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

[Peter Albert]

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


-- 
5 GB Mailbox, 50 FreeSMS http://www.gmx.net/de/go/promail
+++ GMX - die erste Adresse für Mail, Message, More +++