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

[Noel Zinn]

Interesting, Peter.

Don't have a solution for you, but I have duplicated your MuPAD computations 
with similar numerical integrations in Matlab.  So, that's cool.  But I've 
also dusted off a meridional arc function from an old TM code and duplicated 
Proj.4's computation as well.  So, that's cool, too.

So why the difference?  I don't doubt Proj.4 or my TM code.  Too much usage 
behind them.  Rather, my intuition suggests that either MuPAD and Matlab are 
integrating over a prolate (egg-shaped) ellipsoid (not oblate) or that the 
zero latitude is interpreted differently (vertically?).  Don't know why that 
would be so, but worth investigating.  Too late tonight, however.

Good luck,
Noel


----- Original Message ----- 
From: "Peter Albert" <peter.albert at gmx.de>
To: <proj at xserve.flids.com>
Sent: Tuesday, July 12, 2005 10:24 AM
Subject: [Proj] wrong computation of meridinal distance / pj mlfn.c


> 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 +++
> _______________________________________________
> Proj mailing list
> Proj at xserve.flids.com
> http://xserve.flids.com/mailman/listinfo/proj
>