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

[Strebe at aol.com]

At this very moment I happen to be sitting in the first session of Map 
Projections oral presentations at the ICC biannual conference of the International 
Cartographic Association, so my attention and resources are a little limited. 
However, it seems to me that you are computing the wrong integral. The integral 
you give calculates the arc length from the pole to 10 degrees south of it, 
not from the equator. That seems quite apparent from the way you set it up. You 
need to subtract your result from the complete elliptic integral.

I have checked the proj result with my own completely independent 
implementation. The proj result is correct.

Regards,

daan Strebe
Chair, Commission on Map Projections
International Cartographic Association

----- 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
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20050713/7c1439eb/attachment.html