[Proj] wrong computation of meridinal distance / pj mlfn.c
ndzinn at houston.rr.com
Tue Jul 12 22:43:01 EDT 2005
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.
----- 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:
> 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
> 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
> However, if you integrate from 90° - the geocentric latitude up to the
> like this:
> 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,
> 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
More information about the Proj