[Proj] wrong computation of meridinal distance / pj mlfn.c

[Clifford J Mugnier]

Peter:

Geocentric latitude is not the same as geodetic latitude.  Is that possibly
your conundrum?  True meridianal distance on the surface of the ellipsoid
is measured with respect to geodetic latitude.

Geodetic latidude is defined as the angle from the equatorial plane to the
ellipsoid normal terminated by the semi-major axis.  That ellipsoid normal
does NOT intersect the center of the ellipsoid of revolution.

Cliff Mugnier
LOUISIANA STATE UNIVERSITY

--------------------------------

Dear all,

I am new to proj and also relatively new to projection in general;however,
I think I found a problem in how the meridinal distance is approximated in
proj (using the subroutine in pj_mlfn.c).

Actually, it seems to me as if the meridinal distance on a given ellipsoid
is not calculated from the equator to the given (geodetic) latitude but
between the Pole and (90° - the appropriate geocentric latitude)!

Let me give you an example: For the sake of simplicity, let the major axis
of our ellipsoid be a = 1.0 and the squared eccentricity es = 0.007.

Let the geodetic latitude be 10N, then the appropriate geocentric latitude
calculates to 9.931398465N.

(Following this transformation:
float(arctan((a-es*a)/a*tan(10*PI/180))*180/PI)
)

Now, if I understand the sinusoidal projection correctly, then the
y-coordinate of a given point in this projection should just equal its
meridinal distance from the equator, in our case:

echo "0 10" | proj +proj=sinu +a=1.0 +es=0.007 -f %.8f | cut -f 2

The result is 0.17332956, which generally makes sense given the major axis
being set to 1. (If I choose a realistic ellipsoid, the resulting value is
just 0.17332956 times a.)

Just to be sure, let's cross-check my calculation of geocentric latitude by
adding "+geoc" to the previous call and using the geocentric latitude
instead of 10N:

echo "0 9.931398465" | proj +geoc +proj=sinu +a=1.0 +es=0.007 -f %.8f |
cut -f 2

gives 0.17332956, no problems so far.

However, now I tried to calculate the elliptical integral myself, i.e.
int(sqrt(a-es*(cos(t)^2)) dt) from the equator up to the given geocentric
latitude. I did this using mupad, a computer algebra system:

numeric::int(sqrt(1.0-0.007*(cos(x))^2), x=0..9.931398465*PI/180)

The result now is 0.1727339231, which significantly differs from the proj
output.

However, if you integrate from 90° - the geocentric latitude up to the pole
like this:

numeric::int(sqrt(1.0-0.007*(cos(x))^2),
x=(90-9.931398465)*PI/180..90*PI/180)

then you get 0.1733295629, which perfectly matches proj's results.



I would greatly appreciate if someone could tell me whether I really got
something wrong here, or if there actually is a problem in proj. I looked
into the code, and the meridinal distance is unfortunately approximated
using a weird 8th order polynomial which I can't say anything about ...

Thanks in advance and kind regards,

Peter


P.S.: BTW, I am using proj-4.4.9