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

[Noel Zinn]

Peter,

Here's the solution to your original query: why MuPAD doesn't match Proj.4 
for the meridional arc.  (It's been an informative thread in the meantime.)

The answer is that integrating over geocentric latitude is incorrect (as 
would be integrating over geodetic latitude).  Given your parametric 
formulation below, the integration parameter 't' is defined as follows:

t = atan(b*tan(phi)/a), where phi is geodetic latitude

The geocentric latitude is:

t = atan(b^2*tan(phi)/a^2), somewhat different

For phi = 10 degrees, t = 9.96564260327, which gives me the answer that 
matches Proj.4 when your function is numerically integrated in Matlab.

I derived this by differentiating your two parametric equation for dx/dt and 
dy/dt and then solved for dy/dx, which is the horizontal to the ellipsoid. 
The vertical is the negative inverse of that.  The arc tangent of the 
vertical is the angle (t) that we're looking for.

Let me know if you agree.

Regards,
Noel


----- Original Message ----- 
From: "Albert Peter" <Peter.Albert at dwd.de>
To: "'PROJ.4 and general Projections Discussions'" <proj at xserve.flids.com>
Sent: Wednesday, July 13, 2005 10:52 AM
Subject: AW: [Proj] Re: wrong computation of meridinal distance / pj mlfn. c


> Gerald,
>
> thanks for your answer.
>
>> For the math involved see the section in the libproj4 manual (p. 24)
>> on the web site
>> http://members.verizon.net/~vze2hc4d/proj4
> <http://members.verizon.net/~vze2hc4d/proj4>
>
> that helped a bit further: when I plug your formula into MuPAD I can at
> least reproduce proj's results, which is good. So it is not the mysterious
> polynomial approximation anymore. I will go and try to understand where
> "your" formula comes from.
>
>> Note that all references to latitude in the manual are to the geodetic
>> latitude.
>> Geocentric is only of interest to satellite problems and I can't
>> imagine that substituting the geocentric form into the integral would
>> simplify the solution.
>
> Well, I actually did not substitute it but rather started from the
> definition of the ellipse as
>
> x=a*cos(t) and y=b*sin(t),
>
> where for t the geocentric latitude has to be taken.
>
> Then the arc length ds id
>
> ds = sqrt( (dx/dt)^2 + (dy/dt)^2 ) dt
>
> which brought me to the formuala I used.
>
> So, currently I don't see where the mistake in this approach lies.
>
> Best regards,
>
>   Peter
>
>
> _______________________________________________
> Proj mailing list
> Proj at xserve.flids.com
> http://xserve.flids.com/mailman/listinfo/proj
>