Low Distortion Projection help needed - solving for "k" aka Grid Scale Factor at a Point

I always use TM for my LDPs, Lamberts don't play nice with many available programs. To get K I simply compare ground to ellipsoid distances. I keep everything NAD83.

FWIW,

Mr Dennis derived his equation from Stemƒ??s manual ( https://www.ngs.noaa.gov/PUBS_LIB/ManualNOSNGS5.pdf ). Stemƒ??s manual solves for ƒ??kƒ? in parts; I think it is easier to code that way.

BTW, ƒ??lnƒ? refers to natural (not base 10) logarithms.

Also you indicate a value for ƒ??eƒ? which is the value for ƒ??e^2ƒ?. In Mr Stemƒ??s manual he uses an e^2 of 0.0066943800229034. This is a derived value and not a defining one for GRS80.

I cannot edit my post.

I just wanted to add, on the issue of the value of ƒ??e^2ƒ? that Dr Milbert, former NGS Chief Geodesist, iteratively computed the value in a post to this site.

He reports: ƒ??e2=6.694380022903415749574948586457761E-0003ƒ? He says he used quadruple precision. See #post-280228

Dr Milbertƒ??s post included the other derived quantities for GRS80.?ÿ

If, in addition to having the satisfaction of completing it yourself, you wanted an answer by other means, Leica's current Infinity software will show k for any point in its Points tab. LGO didn't do this (at least not for custom projections), but Infinity will. Maybe even a demo copy of Infinity could help you out.

I took a different approach for LDPs. If you look at Stem's page 28, the formulas for k0?ÿand R0 can be rewritten to bypass a lot of the computation overhead. In designing an LDP, you are going to assign a value to k0 so you don't need to solve for it.

Stem's publication is here:?( https://www.ngs.noaa.gov/PUBS_LIB/ManualNOSNGS5.pdf /p>

What you do need is R0, the mapping radius of the central parallel and K, the mapping radius of the equator. I get R0 by rewriting Stem's formula for k0 this way:

R0 = k0 * a/(w0 * tan(phi0)), where phi0 is the latitude of the central parallel.

Then K = R0 * exp(Q0 * sin(phi0)).

Q0 and w0 are defined on page 27. You'll see a lot of the elements of Michael's formula in the formula for Q.

Then k can be calculated by Stem's formula for k on page 28.

My eyes are old and rheumy and my fingers are stiff and numb, so getting parentheses in the right places to make a complicated formula work in a spreadsheet is just too difficult. Breaking it down into intermediate pieces helps me get it right.

I hope it helps you, too.

?ÿ

?ÿ

?ÿ

Math teacher has it right, you chose the number, keep it simple, think of it in terms of ppm

I donƒ??t whether there are any Matlab users but I just found out today that Clive Moler is working to develop a quadruple precision capability!

?ÿ

GeeOddMike, MightyMoe, MathTeacher, FrozenNorth : thank you all for your input, clarifications, hints, reminders, and re-written formulas. I truly appreciate it all and everything is helpful.?ÿ

?ÿ

And, I'm having a fun time figuring this out 🙂

Just testing Greek letters and subscripts. I can see the Greek letter phi on my screen but not the result of the subscript code. I'm trying to make a better representation of phi0,

??₀

?ÿ

?ÿ

?ÿ

So how do I get the subscript 0 in the reply?

?ÿ@MathTeacher what is w0, or even w for that matter? I don't see it on page 27 of the manual, but I do see the formula for Ws.. Thanks :)?ÿ

???ÿ

Alt +227

@kjypls?ÿw?ÿis the denominator in the formula for the length of a line from the semi-minor axis perpendicular (normal) to the surface of the ellipsoid. Every point has one; ws is at the latitude of the southern standard parallel, so the formula uses the sine of that latitude. At the central parallel, w0 uses the sine of the central parallel.

If you look at the formulas for?ÿk0 and?ÿk on page 28, you can gain some insight. Note that?ÿk0?ÿuses the tangent function while?ÿk?ÿuses the sine of one latitude and the cosine of another. The formula for?ÿk0 uses w0 while the formula for?ÿk uses the definition for w. They are the same formula, but, in the case of?ÿk0, the sine and cosine are functions of the same angle and their quotient is the tangent of that angle. Not sure why the w in one and the definition of w in the other.

Hope this helps and good luck in the development.

GeeOddMike, Mathematica gives you any precision you want. Here is pi to 50 places (in less than 1ms)

pi = 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 3751

and its cheaper than Matlab.

Program step in Mathematica is?ÿ N[pi,50]?ÿ?ÿ?ÿ If ones needs more places, say 10,000 then N[pi,10000]

A small correction for my post above. Inside the square brackets in my post " pi " the first letter "p" needs to be

a Capital letter P.?ÿ It should read N[Pi,50]?ÿ or?ÿ N[Pi,10000].

?ÿ

JOHN NOLTON

?ÿ

As a longtime user of?ÿ Matlab I can whip out a script in short order for most tasks of interest to me. Like spoken languages, my ability to deal with more than one language has proven increasingly problematic.?ÿ

?ÿ

As far as cost goes, the Matlab home edition is only $149. The Mathematica Home Editions ( http://www.wolfram.com/mathematica-home-edition/ ) is $320. Both tools try to get you to sign up for service plans. BTW, there are even a few free tools like Scilab that are readily available..?ÿ?ÿ

Both tools are a great improvement over the days when coding from scratch using FORTRAN, C, or C++.?ÿ I never much cared for coding problems using spreadsheets.

Thanks for the info.?ÿ

DMM

?ÿ