Doubts in converting spherical coordinates

I'm trying to convert spherical coordinates (namely latitude and longitude from a GPS device) into Cartesian coordinates. I'm following this simple conversion derived from the polar coordinates conversion equations.

Then I'm calculating the distance between the two point applying the euclidean distance but the value I'm finding is not always the same as the distance I can calculate using the haversine formula. In particular I'm noticing that given different longitudes but same latitudes leads to the same distances computed by the two algorithms, whereas having the same longitude and changing the latitude carries different values.

Here is the C code I am using:

double ComputeDistance(double lat1,double lon1, double lat2, double lon2)
{
    
    double dlon, dlat, a, c;
    dlon = lon2- lon1;
    dlat = lat2 - lat1;
    a = pow(sin(dlat/2),2) + cos(lat1) * cos(lat2) * pow(sin(dlon/2),2);
    c = 2 * atan2(sqrt(a), sqrt(1-a));
    return 6378140 * c;  /* 6378140 is the radius of the Earth in meters*/

    
}
int main (int argc, const char * argv[]) {

    double lat1 = 41.788251028649575;
    double lat2 = 41.788251028649575;
    double long1 = -118.1457209154;
    double long2 = -118.1407209154;//just ~10 meters distant
    
    lat1 = DEGREES_TO_RADIANS(lat1);
    lat2 = DEGREES_TO_RADIANS(lat2);
    long1 = DEGREES_TO_RADIANS(long1);
    long2 = DEGREES_TO_RADIANS(long2);
    
    //transform in cartesian coordinates
    double x = 6378140 * cos(lat1) * cos(long1);
    double y = 6378140 * cos(lat1) * sin(long1);
    
    double x2 = 6378140 * cos(lat2) * cos(long2);
    double y2 = 6378140 * cos(lat2) * sin(long2);

    
    double dist = sqrt(pow(x2 - x, 2) + pow(y2 - y, 2));
    printf("DIST %lf\n", dist);
    printf("NDIST %lf\n", ComputeDistance(lat1, long1, lat2, long2));

    return 0;
}

Am I doing something incorrect or there is some math behind it I am not seeing (and maybe ask this on Mathoverflow boards?). UPDATE It is not necessary to cross boards, as someone correctly pointed this conversion is not meaningful for computing the exact distance between two points (the distance between the two poles is zero). So I am reformulating it as: why at small deltas (0.0001 which corresponds to 10 mts more or less) of latitudes the distance appear to be so different from the haversine formula (20-25%)?

UPDATE 2: As Oli Charlesworth pointed out, not considering the z axis makes this conversion a projection that does not mind the north-south difference. This is also the cause of the difference in deltas I was pointing out. In fact, in the correct transform the z is related to the latitude and if you consider it , then compute the euclidean distance between the two points (in a 3d space now), both latitude and longitude will lead to a good approximation for small deltas. For Example for a degree of latitude the error is ~ 1.41 meters.

From this, there is no 2D map projection where distance is preserved. Calculating the distance from the 2D projection of points is useless.