Geodetic Integral Formulation [Long: FWF Req'd]
Kaimbridge M. GoldChild
In article <3A6B08...@bc.sympatico.ca>,
Te...@bc.sympatico.ca wrote:
> Kaimbridge wrote:
> > Sure there are, you just have to calculate the two integrals
>
> What are the two integrals? Where can I find them written out?
> I have been looking for fundamental formulae for this purpose.
<snip>
> I would appreciate knowing the integral formulations.
Okay, let's start off with a few definitions, redefinitions, and
clarifications.
At the very onset, it should be made clear that, except for the Greek
variables noted below--or designated proxies (i.e., "ğ", "ã",
"Ğ")--and the "internal variables" (Vincenty's A, B, C; Sodano's A_o,
B_o, C_o, D_o and A', B', C'; Rainsford's A_0, A_2, A_4, A_6), all of
the integratory assignments/designations (T{IP}, L{IP}, Fm{MN,UN},
etc.) are purely arbitrary, but mnemonically designed creations of
this writer (with some designations still in contention--e.g.,
previously defining integral kernels T{IP}, L{IP}, etc., as T'{IP},
L'{IP}, etc.)--however, the mathematical integrity of the equations
themselves is of the highest caliber.
In an attempt to keep the vagaries and complexities of calculus to
a minimum, the definite integrals involved will be identified and
worked in mean value kernel format: Definite integral = F{IP} * BD,
integrand = F{P_tn}. While kernel integration is defined here using
"natural" quadrature (P_tn = LB to UB in an increment of Iota), the
ideal means of general, applicational integration is via Gaussian
Quadrature, which is discussed at the end of this presentation.
As for analyzing and experimenting with the various equations presented
in this document, an excellent desktop "workbench" is UBasic, a free,
math friendly, traditional BASIC, which can be downloaded at its
author's site:
ftp://rkmath.rikkyo.ac.jp/pub/ubibm/
A major reconceptualization rarely acknowledged or used in geodesy,
but believed rudimentary by this presenter, is the angularization of
the elliptic parameters:
a,b = equatorial, polar radii;
cos{Oz} = b/a; sin{Oz}^2 = e^2 = [a^2 - b^2]/a^2; -:[1]:-
tan{Oz}^2 = e'^2 = [a^2 - b^2]/b^2;
ver{Oz} = 2 * sin{.5*Oz}^2 = f = [a - b]/a;
The online formula being referenced here is Thaddeus Vincenty's "Direct
And Inverse Solutions Of Geodesics On The Ellipsoid With Applications
Of Nested Equations", originally published in Survey Review (April
1975, Vol. XXIII, No. 176, pp. 88-93), as presented by the (Australian)
Intergovernmental Committee on Surveying and Mapping:
http://www.anzlic.org.au/icsm/gdatm/chapter4.htm
Let ğ = "sigma"; ã = "alpha"; Ğ = "Delta" (difference triangle);
LB,UB = ğ_1,ğ_2 ("angular distance on a sphere, -:[2]:-
from the equator to point_1,point_2");
BD = UB - LB = ğ = ğ_2 - ğ_1; BS = UB + LB = 2ğ_m;
AP = ã ("azimuth of the geodesic at the equator"),
= asin{cos{reduced Lat_1} * |sin{ã_(1-2)}|},
= asin{cos{reduced Lat_2} * |sin{ã_(2-1)}|};
(ã_(1-2),ã_(1-2) = forward,reverse local azimuths)
IP = integrated value of function's point range from LB to UB;
(i.e., the average, or the mean value theorem's "c")
There are two interrelated benchmark integrands used in the inverse and
direct principle problems of geodesy: L{P_tn} for longitudinal
spheroidification and T{P_tn} for the actual distance. Since L{P_tn}
utilizes T{P_tn}, it is the distance integral that will be looked at
first.
In finding geodetic distance, one must first convert the geographical
coordinates/waypoints into spheroidal values, to create an "auxiliary"
(the "conformal"?) sphere. While the latitudinal conversion is a simple
closed function reduction, the longitudinal spheroidification is
considerably more complex, as it involves both a definite integral and
anti-differentiation, by which it can be found either iteratively or
non-iteratively: Given the complexity of the non-iterative (involving
the Maclaurin series), the iterative method is usually preferred,
particularly if significant ellipticity and/or angular distance is
involved.
To the precision that A and Ğğ are convergent to their respective
series, as well as to their approximative accuracy, Vincenty's
equation for distance works out to the following:
s = b * A * [ğ - Ğğ], -:[3]:-
= b * [(sin{IPg} * cos{AP} * tan{Oz})^2 + 1]^.5 * BD,
= a * cos{Oz} * [1 + (sin{IPg} * cos{AP} * tan{Oz})^2]^.5 * BD,
= a * [cos{Oz}^2 + (sin{IPg} * cos{AP} * sin{Oz})^2]^.5 * BD,
= a * T{IPg} * BD;
Therefore, let
-:[4]:-
T{P_tn} = [cos{Oz}^2 + (sin{P_tn} * cos{AP} * sin{Oz})^2]^.5,
= cos{Oz} * [(sin{P_tn} * cos{AP} * tan{Oz})^2 + 1]^.5;
UT=oo T{P_tn}
T{IPg} = sum (=======): P_1 = LB, P_ut = UB;
TN=1 UT (the closer UT is to oo [infinity], the
greater the precision of integration)
For longitudinal spheroidification:
M = w = L = Long_2 - Long_1; -:[5]:-
Y = "Lamda" = M + Å{Y} (spheroidal M);
Å{Y} = sin{AP} * L{IPy} * BD = [1 - C] * f * sin{ã} * [....];
sin{Oz}^2
L{P_tn} = ===========;
1 + T{P_tn}
sin{Oz}^2 UT=oo L{P_tn}
L{IPy} = ========== = sum (=======): P_1 = LB, P_ut = UB;
1 + T{IPy} TN=1 UT
(of course, though, L{IPy} = sin{Oz}^2 / [1 + T{IPy}],
NOT sin{Oz}^2 / [1 + T{IPg}]), as L{IPy} is an
"inverse harmonic" integral, meaning dp/F{p} is
found from sum(1/F{P_tn}), NOT 1/sum(F{P_tn}), thus
you can't just calculate L{IPy} as a simple closed
function, using T{IPg}--hence, IP is split into
distinct IPg, IPy identities!)
Given the relationship between L{P_tn} and T{P_tn}, these can be
considered the benchmark integrand forms of the needed integrals.
A working online execution of Vincenty's A_k, B_k, C formulation (A_k,
B_k being Helmert's enhanced, k_1 based A, B values discussed later in
this paper) can be found at:
http://www.best.com/~williams/gccalc.htm
But this is just the "bare bones" explanation, suitable for a
platonic understanding, using numerical integration: The structure
and theory of T{P_tn} is actually more complex, including a more
direct approach involving the AP reduction of P_tn [RP_tn] to create
the alternatively defined, but fully equivalent, T{IPr}, as well
as there being an alternate, more fundamental, elliptically based
kernel--R{P_tn} -> R{IPg}--where the reduction is more properly
assigned to Oz [Rz], thus T{IPg/IPr} can be considered the "transferred
reduction" integral kernel (though, further complicating its anatomical
layout, using the here defined standard form of the integrand, the
reduction is retransferred back to Oz in modified form [Tz]!).
Furthermore, most of the popular fornmularies--including Vincenty's and
Emanuel Sodano's--use binomial series expansions (or approximations of)
for finding the integrals.
So now let's take a second, closer look.
-:[6]:-
tan{Rz}^2 = cos{AP}^2 * tan{Oz}^2,
= u^2 = cos{ã}^2 * [a^2 - b^2]/b^2 = [Ra^2 - b^2]/b^2;
(i.e., Rz = atan{cos{AP}*tan{Oz}})
Ra = a * [1 - (sin{AP} * sin{Oz})^2]^.5 = b * sec{Rz};
R{P_tn} = cos{Rz} * [(sin{P_tn} * tan{Rz})^2 + 1]^.5,
= [1 - (cos{P_tn} * sin{Rz})^2]^.5,
= [cos{Rz}^2 + (sin{P_tn} * sin{Rz})^2]^.5;
sin{Tz} = cos{AP} * sin{Oz};
sin{RP_tn} = cos{AP} * sin{P_tn};
cos{RP_tn} = [1 - sin{RP_tn}^2]^.5;
T{RP_tn} = cos{Oz} * [(sin{RP_tn} * tan{Oz})^2 + 1]^.5,
= [1 - (cos{RP_tn} * sin{Oz})^2]^.5,
= [cos{Oz}^2 + (sin{RP_tn} * sin{Oz})^2]^.5,
= T{P_tn} = [cos{Oz}^2 + (sin{P_tn} * sin{Tz})^2]^.5;
Thus R{P_tn} = cos{Rz} * sec{Oz} * T{(R)P_tn};
T{(R)P_tn} = cos{Oz} * sec{Rz} * R{P_tn};
(BUT, cos{RP_tn} * sin{Oz} x=x cos{P_tn} * sin{Tz})
T{IPr} = T{IPg} = cos{Oz} * sec{Rz} * R{IPg};
and Cr = a * T{IPg} = a * T{IPr}, -:[7]:-
= Ra * R{IPg} = b * sec{Rz} * R{IPg};
s = Cr * BD;
While the numerical integrand forms may seem a bit complicated, the
binomial series expansions can be even more confusing and convoluted,
particularly as defined in most presentations!
NB = TN_bs (binomial series term number) = 0 or 1 to UN;
UN = UT_nb (upper term number limit of the series);
MN = TN_cm (cosine multiple term number);
"Binomial Coefficients", "Binomial Series": -:[8]:-
P{17,5} 13 * 14 * 15 * 16 * 17
(e.g.) BC{17,5} = ======= = ======================;
5! 1 * 2 * 3 * 4 * 5
P{.5,3} -1.5 * -.5 * .5
BC{.5,3} = ======= = ===============;
3! 1 * 2 * 3
1 -1 1
H1 = BC{.5,1} = =; H2 = BC{.5,2} = =; H3 = BC{.5,3} = ==;
2 8 16
[1 + (-Q)]^.5 = 1 + [H1*(-Q)^1 ] + [H2*(-Q)^2] + [H3*(-Q)^3] + ...,
= [1 - Q]^.5 = 1 - [H1*Q^1] + [H2*Q^2] - [H3*Q^3] + ...,
1 1 1
= 1 - =Q - =Q^2 - ==Q^3 + ...;
2 8 16
(unless noted otherwise, all closed functions -> binomial series
throughout this paper are expanded out to a UN = 3 equivalency)
--------------------------------------------------------------------
cos{Q}^2 = [1 + cos{2*Q}]/2^1; -:[9]:-
sin{Q}^2 = [1 - cos{2*Q}]/2^1;
cos{Q}^4 = [3 + (4 * cos{2*Q}) + cos{4*Q}]/2^3;
sin{Q}^4 = [3 - (4 * cos{2*Q}) + cos{4*Q}]/2^3;
cos{Q}^6 = [10 + (15 * cos{2*Q}) + (6 * cos{4*Q}) + cos{6*Q}]/2^5;
sin{Q}^6 = [10 - (15 * cos{2*Q}) + (6 * cos{4*Q}) - cos{6*Q}]/2^5;
cos{Q}^(2*NB) = [ (.5 * BC{2*NB,NB-0} * cos{2*0*Q})
NB
+ SUM BC{2*NB,NB-MN} * cos{2*MN*Q}]/2^([2*NB]-1);
MN=1
sin{Q}^(2*NB) = [ (.5 * BC{2*NB,NB-0} * cos{2*0*Q})
NB
+ SUM BC{-1,NB} * BC{2*NB,NB-MN}
MN=1 * cos{2*MN*Q}]/2^([2*NB]-1);
and, rearranging by cosine multiples and adding sin{E} factors in an
inverted half-pyramid, where:
Fm{MN,UN}= cos{2*MN*Q} term's inverted pyramid function;
H1*1 H2*3 H3*10
CSm{0,3} = 1 - ====sin{E}^2 + ====sin{E}^4 - =====sin{E}^6; -:[10]:-
2^1 2^3 2^5
H1*1 H2*4 H3*15
CSm{1,3} = - ====sin{E}^2 + ====sin{E}^4 - =====sin{E}^6;
2^1 2^3 2^5
H2*1 H3*6
CSm{2,3} = + ====sin{E}^4 - ====sin{E}^6;
2^3 2^5
H3*1
CSm{3,3} = - ====sin{E}^6;
2^5
Or, doing out, combining and simplifying all of the coefficients,
1 3 5
CSm{0,3} = 1 - =sin{E}^2 - ==sin{E}^4 - ===sin{E}^6; -:[11]:-
4 64 256
1 1 15
CSm{1,3} = - =sin{E}^2 - ==sin{E}^4 - ===sin{E}^6;
4 16 512
1 3
CSm{2,3} = - ==sin{E}^4 - ===sin{E}^6;
64 256
1
CSm{3,3} = - ===sin{E}^6;
512
(likewise, the sin{Q}^XP equivalent inverted pyramid uses
the same exact terms as the cosine pyramid, except that
the polarities are all positive for even MN rows and
negative for odd ones and sin{E}^XP changes to tan{E}^XP:
STm{0,3} = 1 + CA_1*tan{E}^2 + CA_2*tan{E}^4 + CA_3*tan{E}^6;
STm{1,3} = - CB_1*tan{E}^2 - CB_2*tan{E}^4 - CB_3*tan{E}^6;
STm{2,3} = + CC_2*tan{E}^4 + CC_3*tan{E}^6;
STm{3,3} = - CD_3*tan{E}^6;)
from which [1 - (cos{Q} * sin{E})^2]^.5
= 1 - [H1 * cos{Q}^2 * sin{E}^2]
+ [H2 * cos{Q}^4 * sin{E}^4]
- [H3 * cos{Q}^6 * sin{E}^6],
= CSm{0,3} + [CSm{1,3} * cos{2*1*Q}]
+ [CSm{2,3} * cos{2*2*Q}]
+ [CSm{3,3} * cos{2*3*Q}];
and [(sin{Q} * tan{E})^2 + 1]^.5 -:[12]:-
= 1 + [H1 * sin{Q}^2 * tan{E}^2]
+ [H2 * sin{Q}^4 * tan{E}^4]
+ [H3 * sin{Q}^6 * tan{E}^6],
= STm{0,3} + [STm{1,3} * cos{2*1*Q}]
+ [STm{2,3} * cos{2*2*Q}]
+ [STm{3,3} * cos{2*3*Q}];
(and [(sin{Q} * tan{E})^2 + 1]^.5
= sec{E} * [1 - (cos{Q} * sin{E})^2]^.5)
--------------------------------------------------------------------
The calculus modification of the binomial series involves a redefining
of cos{2*MN*Q}:
2 * sin{MN*BD} * cos{MN*BS}
CM'{MN} = cos'{2*MN*IP_mn} = ===========================; -:[13]:-
2 * MN * BD
While the references seen by this author express the integrand as
"cos{2*MN*P}" (with maybe one or two sources denoting it as
"sin{2*MN*P}", which, given the "sin{MN*BD} * cos{MN*BS}" element in
CM'{MN}, supports the notion of the actual identity of the integrand
as being sin{2*MN*P}--i.e., cos'{2*MN*P}--rather than cos{2*MN*P}),
if one tries to integrate cos{2*MN*P} numerically (i.e., mean value),
cos{2*MN*IP} x=x cos{2*MN*IP_mn}, as it should--this writer is not sure
if this is because either it is not a "real" definite integral or it is
a special subtype: In either case, anyone with a clarification or
explanation regarding this discrepancy is more than welcome to respond.
It is the latter binomial series that Sodano uses for sec{Rz} * R{IPg}
in his first of two major papers ("A Rigorous Non-iterative Procedure
For Rapid Inverse Solution Of Very Long Geodesics", Bulletin Geodesique,
1958, No.s 47-48, pp. 13-25), replacing tan{E} with tan{Rz} and shifting
CM'{MN}'s MN divisor from CM'{MN} to A_o, B_o, C_o, D_o (as well as some
polarity inversions--i.e., + -> - and vice versa--which were necessary
due to both H2 being negative as well as some quirky arrangements he
gave for some of the spherical components used, but which I've corrected
below):
(for purposes of neatness and spatial compression, let
OS = sin{Oz}^2, OT = tan{Oz}^2, OV = ver{Oz},
RS = sin{Rz}^2, RT = tan{Rz}^2; CP = cos{AP}^2);
SI_1 = sin{IP_1}^2 = [1 - CM'{1}]/2^1; -:[14]:-
SI_2 = sin{IP_2}^4 = [3 - (4 * CM'{1}) + CM'{2}]/2^3;
SI_3 = sin{IP_3}^6 = [10 - (15 * CM'{1}) + (6 * CM'{2}) - CM'{3}]/2^5;
(CI_nb = cos{IP_nb}^(2*NB) = [...]/2^([2*NB]-1))
*** Since the CM'{MN} based series expansions used to find the
exponentially unique sine-cosine functions of IPg and IPr are also
used to find IPy's and the yet to be introduced IPf's, the P based
integrated kernels for a specific binomial term will be denoted as
IP_nb and the RP based as IR_nb (e.g., cos{IPg}^6 = cos{IP_3}^6 and
cos{IPy}^6 = cos{IPr}^6 = cos{IR_3}^6). In turn, the reverse applies
to analytical but noncalculable conversions and cancellations (e.g.,
CI_3 x=x [1-SI_1]^3 and SI_3/SI_1 x=x SI_2, thus they will be
expressed as CI^3 = [1-SI]^3 and SI^3/SI = SI^2). ***
H1*1 H2*3 H3*10
A_o = RTm{0,3} = 1 + ====RT^1 + ====RT^2 + =====RT^3; -:[15]:-
2^1 2^3 2^5
H1*1 H2*4 H3*15
1 * B_o = RTm{1,3} = - ====RT^1 - ====RT^2 - =====RT^3;
2^1 2^3 2^5
H2*1 H3*6
2 * C_o = RTm{2,3} = + ====RT^2 + ====RT^3;
2^3 2^5
H3*1
3 * D_o = RTm{3,3} = - ====RT^3;
2^5
In Vincenty's formulation, however, more than adequate (for earthly
applications) approximations for C_o and D_o are used, instead:
-:[16]:-
A = A_o = RTm{0,3},
= [(sin{IPg} * tan{Rz})^2 + 1]^.5: P_1 = 0, P_2 = 90°;
(i.e., AP adjusted elliptic integral of the second
kind, to an expansion out to UN = 3)
B = -B_o / A = -RTm{1,3} / RTm{0,3};
By introducing an approximation of RTm{MN>0,UN}, one gets:
RTa{MN>0,UN} = 2 * BC{.5,MN} * RTm{0,UN}^(1-MN) * RTm{1,UN}^MN;
-:[17]:-
RTm{1,3} = 1 * B_o = RTa{1,3} = 2 * H1 * A * (-B)^1;
RTm{2,3} = 2 * C_o ~=~ RTa{2,3} = 2 * H2 * A * (-B)^2;
RTm{3,3} = 3 * D_o ~=~ RTa{3,3} = 2 * H3 * A * (-B)^3;
and [(sin{IPg} * tan{Rz})^2 + 1]^.5 (P_1 = LB, P_2 = UB) -:[18]:-
= 1 + [H1 * SI_1 * RT^1] + [H2 * SI_2 * RT^2]
+ [H3 * SI_3 * RT^3],
= RTm{0,3} + [RTm{1,3} * CM'{1}] + [RTm{2,3} * CM'{2}]
+ [RTm{3,3} * CM'{3}],
= A_o + [(1 * B_o) * CM'{1}] + [(2 * C_o) * CM'{2}]
+ [(3 * D_o) * CM'{3}];
~=~ RTm{0,3} + [RTa{1,3} * CM'{1}] + [RTa{2,3} * CM'{2}]
+ [RTa{3,3} * CM'{3}],
= A + [2 * H1 * A * (-B)^1 * CM'{1}]
+ [2 * H2 * A * (-B)^2 * CM'{2}]
+ [2 * H3 * A * (-B)^3 * CM"{3}],
= A * [ğ - Ğğ]/ğ
There is a much more efficient series expansion for R{IPg},
introduced in Friedrich Helmert's benchmark classic, "The
Mathematical And Physical Theory Of Higher Geodesy" (originally
published in German, 1880; translated by USAF in 1964), based on
the integrand model cos{.5*Rz}^2 * [1 - (2 * cos{2*P_tn} * tan{.5*Rz}^2)
+ tan{.5*Rz}^4]^.5. However, besides the unidentifed function
involving CM'{MN} being more complex, there is no equivalent
integrand model for L{IPy}.
It *is* perfectly useful, though, for increasing the integrative
efficiency of A (and, thus, is the ideal model for finding the elliptic
integral of the second kind) and B, which Vincenty presented in a
subsequent addendum to his original article (Survey Review, April 1976,
pg.294), using Helmert's k_1 and notated here as A_k and B_k,
respectively:
[1 + u^2]^.5 - 1
tan{.5*Rz}^2 = k_1 = ================; -:[19]:-
[1 + u^2]^.5 + 1
A_k = RHm{0,1} = sec{Rz} * cos{.5*Rz}^2 * [1 + (H1^2 * tan{.5*Rz}^4)];
B_k = -RHm{1,1}/RHm{0,1} = tan{.5*Rz}^2 * [1 - (.375 * tan{.5*Rz}^4)];
While the series expansion of [(sin{IPr} * tan{Oz})^2 + 1]^.5 is
identical to that of [(sin{IPg} * tan{Rz})^2 + 1]^.5, where only
the assignment of cos{AP} changes,
SR_nb = CP^NB * SI_nb = sin{IR_nb}^(2*NB);
-:[20]:-
SR_nb * OT^NB = SI_nb * RT^NB;
the complexity of cos{IR_nb}^(2*NB) makes the expansion of
[1 - (cos{IPr} * sin{Oz})^2]^.5 much more laborous:
CR_nb = cos{IR_nb}^(2*NB) = [1 - sin{IR}^2]^NB;
*** VIA BINOMIAL EXPANSION, using SR_nb ***
-:[21]:-
CR_1 = cos{IR_1}^2 = 1 - SR_1;
CR_2 = cos{IR_2}^4 = 1 - [2 * SR_1] + SR_2;
CR_3 = cos{IR_3}^6 = 1 - [3 * SR_1] + [3 * SR_2] - SR_3,
= 1 - [3 * CP * (1 - CM'{1})/2^1]
+ [3 * CP^2 * (3 - <4 * CM'{1}> + CM'{2})/2^3]
- [CP^3 * (10 - <15 * CM'{1}> + <6 * CM'{2}> - CM'{3})/2^5],
3*1 3*3 1*10
= 1 - [===CP - ===CP^2 + ====CP^3]
2^1 2^3 2^5
3*1 3*4 1*15
+ [(===CP - ===CP^2 + ====CP^3) * CM'{1}]
2^1 2^3 2^5
3*1 1*6
+ [(===CP^2 - ===CP^3) * CM'{2}]
2^3 2^5
1*1
+ [(===CP^3) * CM'{3}];
2^5
Likewise, assigning TSm{MN,3} as the layer function, the inverted
pyramid is just as involved:
1*1
(e.g.) TSm{1,3} = - [H1 * (===CP) * OS^1] -:[22]:-
2^1
2*1 1*4
+ [H2 * (===CP - ===CP^2) * OS^2]
2^1 2^3
3*1 3*4 1*15
- [H3 * (===CP - ===CP^2 + ====CP^3) * OS^3];
2^1 2^3 2^5
Hence,
-:[23]:-
[1 - (cos{IPr} * sin{Oz})^2]^.5
= 1 - [H1 * CR_1 * OS] + [H2 * CR_2 * OS^2]
- [H3 * CR_3 * OS^3],
= TSm{0,3} + [TSm{1,3} * CM'{1}] + [TSm{2,3} * CM'{2}]
+ [TSm{3,3} * CM'{3}];
As the convergency rate of CR_nb * OS^NB is slower than that of
CI_nb * RS^NB, the obvious question is why would one bother going
to the trouble of finding CR_nb: For simply finding T{IPr}, one
probably wouldn't (just let T{IPr} = cos{Oz} * sec{Rz} * R{IPg}).
However, in finding L{IPy}, it becomes quite necessary.
As T{0} = cos{Oz} and sin{Oz}^2 = 4 * sin{.5*Oz}^2 * cos{.5*Oz}^2
= ver{Oz} * [2 * cos{.5*Oz}^2],
sin{Oz}^2 sin{Oz}^2 sin{Oz}^2
ver{Oz} = ================ = =========== = ========= = L{0}:
2 * cos{.5*Oz}^2 1 + cos{Oz} 1 + T{0}
sin{Oz}^2 1 - T{IPy}
L{IPy} = ========== = ==========, -:[24]:-
1 + T{IPy} cos{IR}^2
1 - [1 - (H1 * CR * OS) + (H2 * CR^2 * OS^2) - ...]
= ===================================================,
CR
= [H1 * OS] - [H2 * CR_1 * OS^2] + [H3 * CR_2 * OS^3],
= LSm{0,3} + [LSm{1,3} * CM'{1}] + [LSm{2,3} * CM'{2}];
(since CR_(nb-1) is used in term NB,
the inverted pyramid only extends out
to LSm{UN-1,UN}, instead of LSm{UN,UN})
1
where LSm{0,3} = + [H1 * OS] - [H2 * (1 - ===CP) * OS^2] -:[25]:-
2^1
2*1 3
+ [H3 * (1 - ===CP + ===CP^2) * OS^3];
2^1 2^3
1
LSm{1,3} = - [H2 * ( + ===CP) * OS^2]
2^1
2*1 4
+ [H3 * ( + ===CP - ===CP^2) * OS^3];
2^1 2^3
1
LSm{2,3} = + [H3 * ( + ===CP^2) * OS^3];
2^3
While this is the most direct model, there is another, ver{Oz} based,
binomial series expansion presented by Hume Rainsford ("Long Geodesics
On The Ellipsoid", Bulletin Geodesique, Sept 1955, No.37, pp.12-21),
from which Vincenty extracts his "C" based approximations (though the
actual integrand model and integrating function--"SV_nb"--has yet to
be identified by this writer):
SV_1 = SV{IR_1}^2 = SR_1; -:[26]:-
SV_2 = SV{IR_2}^4 = [4 * SR_1] - [4 * SR_2];
SV_3 = SV{IR_3}^6 = [8 * SR_1] - [18 * SR_2] + [10 * SR_3];
SV_nb = SV{IR_nb}^(2*NB) = ???;
1*1
A_0 = LVm{0,3} = 1 - [H1 * (===CP) * OV^1] -:[27]:-
2^1
4*1 4*3
+ [H2 * (===CP - ===CP^2) * OV^2]
2^1 2^3
8*1 18*3 10*10
- [H3 * (===CP - ====CP^2 + =====CP^3) * OV^3];
2^1 2^3 2^5
~=~ LVm{0,2} = 1 - C = LVm{0,3} + [H3 * (...) * OV^3];
1*1
1 * A_2 = LVm{1,3} = + [H1 * (===CP) * OV^1]
2^1
4*1 4*4
- [H2 * (===CP - ===CP^2) * OV^2]
2^1 2^3
8*1 18*4 10*15
+ [H3 * (===CP - ====CP^2 + =====CP^3) * OV^3];
2^1 2^3 2^5
~=~ [1 - C] * C^1 = LVa{1,2} = LVm{0,2} * [1 - LVm{0,2}]^1;
4*1
2 * A_4 = LVm{2,3} = + [H2 * ( - ===CP^2) * OV^2]
2^3
18*1 10*6
- [H3 * ( - ====CP^2 + ====CP^3) * OV^3];
2^3 2^5
~=~ [1 - C] * C^2 = LVa{2,2} = LVm{0,2} * [1 - LVm{0,2}]^2;
10*1
3 * A_6 = LVm{3,3} = + [H3 * ( + ====CP^3) * OV^3];
2^5
L{IPy} = ver{Oz} * [1 - (SV{IPy}^2 * ver{Oz})]^.5, -:[28]:-
= ver{Oz} * [1 - (H1 * SV_1 * OV^1) + (H2 * SV_2 * OV^2)
- (H3 * SV_3 * OV^3)],
= ver{Oz} * [LVm{0,3} + (LVm{1,3} * CM'{1})
+ (LVm{2,3} * CM'{2})
+ (LVm{3,3} * CM'{3})],
= ver{Oz} * [A_0 + (1 * A_2 *CM'{1}) + (2 * A_4 * CM'{2})
+ (3 * A_6 * CM'{3})];
~=~ ver{Oz} * [LVm{0,2} + (LVa{1,2} * CM'{1})
+ (LVa{2,2} * CM'{2})],
= ver{Oz} * [1 - C] * [1 + (C^1 * CM'{1})
+ (C^2 * CM'{2})];
= [1 - C] * f * [ğ + C * sin{ğ} * (....)]/ğ
The alternative longitudinal integral, presented by Sodano and based
on Helmert's "k_1" series expansion model (also not yet identified
by this writer), is actually geodetically flawed (thus designated
here as "fL{IPf}"), but, like Vincenty's approximations, is adequate
for earthly applications:
sin{Oz}^2 1 - T{IPy}
While L{IPy} = ========== = ==========,
1 + T{IPy} cos{IR}^2
sin{Oz}^2
fL{IPf} = [2 * sin{.5*Oz}^4] + ======================, -:[29]:-
[sec{Oz} * T{IPf}] + 1
sin{Oz}^2 * [(<SR * OT> + 1)^.5 - 1]
= [.5 * ver{Oz}^2] + ====================================,
SR * OT
sin{Oz}^2 * [(<SI * RT> + 1)^.5 - 1]
= [.5 * ver{Oz}^2] + ====================================,
SI * RT
OS * [(1 + <H1*SI*RT> + <H2*SI^2*RT^2>
+ <H3*SI^3*RT^3> + ...) - 1]
= [.5 * OV^2] + ======================================,
SI * RT
= ([.5 * OV^2] + [OS * H1]) + [OS * H2 * SI_1 * RT^1]
+ [OS * H3 * SI_2 * RT^2],
= OV + [OS * H2 * SI_1 * RT^1]
+ [OS * H3 * SI_2 * RT^2],
= OS * [(.5 * sec{.5*Oz}^2) + (H2 * SI_1 * RT^1)
+ (H3 * SI_2 * RT^2)],
= OS * [LTm{0,3} + (LTm{1,3} * CM'{1})
+ (LTm{2,3} * CM'{2})],
= A' + [1 * B' * CM'{1}] + [2 * C' * CM'{2}];
~=~ OS * [LTm{0,3} + (LTa{1,3} * CM'{1})
+ (LTa{2,3} * CM'{2})];
A' H2*1 H3*3
where == = LTm{0,3} = [.5 * sec{.5*Oz}^2] + ====RT^1 + ====RT^2;
OS 2^1 2^3
-:[30]:-
1 * B' H2*1 H3*4
====== = LTm{1,3} = - ====RT^1 - ====RT^2;
OS 2^1 2^3
2 * C' H3*1
====== = LTm{2,3} = + ====RT^2;
OS 2^3
and LTm{1,3} ~=~ LTa{1,3} = LTm{0,3} * -:[31]:-
[1 - (2 * cos{.5*Oz}^2 * LTm{0,3})]^1;
LTm{2,3} ~=~ LTa{2,3} = LTm{0,3} *
[1 - (2 * cos{.5*Oz}^2 * LTm{0,3})]^2;
LTm{MN,UN} ~=~ LTa{MN,UN},
= LTm{0,UN} * [1 - (2 * cos{.5*Oz}^2 * LTm{0,UN})]^MN;
sin{Oz}^2 sin{Oz}^2
Just as L{P_tn} = =========== = =================================,
1 + T{P_tn} 1 + [cos{Oz} * sec{Rz} * R{P_tn}]
one *can* find L{P_tn} from fL{P_tn}:
sin{Oz}^2
L{P_tn} = =================================================;
sin{Oz}^2 -:[32]:-
1 + [cos{Oz} * (=========================== - 1)]
fL{P_tn} - [.5 * ver{Oz}^2]
However, just as L{IPy} can't be found as a closed function using
T{IPr}, T{IPg} or R{IPg}, neither can it be converted from fL{IPf},
for the same reason:
1 1 + 2 + 3 + 4 + 5
Harmonic Mean = ======================= x=x =================;
1 1 1 1 1 1 5
= * [= + = + = + = + =]
5 1 2 3 4 5
When one attains even a minimal understanding of the modified binomial
series expansions of the needed integral kernels, the obvious question
is: Why even consider binomial series expansion, given the more
conceptually explicit, direct and "natural" quadrature--i.e., numerical
integration:
UM = .5 * [UB + MB] = [.75 * UB] + [.25 * LB];
MB = .5 * [UB + LB] = [.50 * UB] + [.50 * LB];
LM = .5 * [MB + LB] = [.25 * UB] + [.75 * LB];
UT - TN
CAEn_tn = =======; -:[33]:-
UT - 1
F{IP}: P_tn = LB + [CAEn_tn * BD] or UB - [CAEn_tn * BD];
Better: F{IP} = .5 * [F{UI} + F{LI}];
F{UI}: UP_tn = MB + [CAEn_tn * (.5 * BD)];
F{LI}: LP_tn = MB - [CAEn_tn * (.5 * BD)];
-:[34]:-
Even Better: F{UI} = .5 * [F{UUI} + F{UMI}];
F{UUI}: UUP_tn = UM + [CAEn_tn * (.25 * BD)];
F{UMI}: UMP_tn = UM - [CAEn_tn * (.25 * BD)];
F{LI} = .5 * [F{LMI} + F{LLI}];
F{LMI}: LMP_tn = LM + [CAEn_tn * (.25 * BD)];
F{LLI}: LLP_tn = LM - [CAEn_tn * (.25 * BD)];
Better Still (QP_tn = P_tn, UP_tn, etc., QI = IP, UI, etc.),
UT F{QP_tn} .5 * [F{QP_1} + F{QP_ut}]
F{QI} = sum (========) - =========================; -:[35]:-
TN=1 UT - 1 UT - 1
(the trapezoidal mean)
While CAEn_tn is the simplest and purest form of quadrature, it is
the least efficient (even using four points per CAEn_tn with the
trapezoidal modification--though it is pretty effective for the
elliptic integral of the second kind).
The ideal form of quadrature, at least for elliptic integrals, is
Gaussian Quadrature, which uses "nodes"--CAEg_ng--and "weights".
Recalling the caveat expressed early on about this writer's creative
license regarding notation, the Gaussian coefficients have been
repackaged/reformatted into more palpable identities:
NG = TNg; UG = UTg;
+/-node_ng = CAEg_ng = cos{AEg_ng}, where AE is the integratorial
amplitude;
weight_ng = OE_ng/UG, where OE is a second, external ("outside")
amplitude;
Gaussian Quadrature requires a minimum of both a positive and negative
use of each node, thus
UP_ng = MB + [cos{AEg_ng} * (.5 * BD)];
UG F{UP_ng}
F{UI} = sum OE_ng * (========);
NG=1 UG
-:[36]:-
LP_ng = MB - [cos{AEg_ng} * (.5 * BD)];
UG F{LP_ng}
F{LI} = sum OE_ng * (========);
NG=1 UG
F{IP} = .5 * [F{UI} + F{LI}];
The whole nature of Gaussian Quadrature is that the nodes are not
evenly incremented, meaning UP_1 x=x UB and LP_1 x=x LB.
The degree of integrational precision for a given node/weight set is
hard to define, using just the two point, UP-LP, per node method.
Besides a prodigious increase in the per point efficiency, using the
four point, UUP-UMP-LMP-LLP, per node method (where UUP_ng, F{UUI}_ug,
etc., are composed and calculated the same way as UP_ng and F{UI}_ug,
etc.--in comparison to UP_tn and F{UI}_ut--with OE_ng staying as is,
unmodified) also provides a definite measure of minumum precision for a
given set (at least for elliptic integrals): Minumum precision is found
at LB = +/-45°, BD = 180°.
-:[37]:-
(Minimum precision of T{IP} for four points per node)
UG = 3 (10^-12 @ cos{Oz} = .9964; 10^-20 @ cos{Oz} = >.9999)
-----------------------------------------------------------------------
AEg_1 = 21.176901 236399 478395 29°; OE_1 = 29.448511 021643 089397 18°;
AEg_2 = 48.607827 177484 556609 05°; OE_2 = 62.010346 638476 673610 18°;
AEg_3 = 76.194942 064285 687087 31°; OE_3 = 80.428480 879127 199623 04°;
UG = 10 (10^-12 @ cos{Oz} = .4279; 10^-20 @ cos{Oz} = .7391)
-----------------------------------------------------------------------
AEg_1 = 6.720618 907402 419355 92°; OE_1 = 10.092082 693867 176810 75°;
AEg_2 = 15.426618 657841 930044 15°; OE_2 = 23.262905 697588 601378 88°;
AEg_3 = 24.184017 406962 473608 77°; OE_3 = 35.908438 629843 510815 43°;
AEg_4 = 32.953008 836254 641760 51°; OE_4 = 47.714058 239468 106744 17°;
AEg_5 = 41.726372 490316 420428 34°; OE_5 = 58.401656 707906 707988 09°;
AEg_6 = 50.501822 122206 501088 41°; OE_6 = 67.720478 429191 205165 23°;
AEg_7 = 59.278402 422230 505723 13°; OE_7 = 75.452031 929620 389953 51°;
AEg_8 = 68.055636 223185 121417 58°; OE_8 = 81.415073 491728 602071 69°;
AEg_9 = 76.833249 489037 493575 14°; OE_9 = 85.469825 422423 159420 14°;
AEg_10= 85.611062 359897 584696 30°; OE_10= 87.521243 889185 748420 09°;
The actual calculations of AEg_ng and OE_ng involve Legendre Polynomials
and their first derivatives which, while the subroutine is direct and
relatively easy, is beyond the scope of this presentation (as is a
subexpansion of AEg_ng, which essentially just increases the number of
points per node): However, UG_3 should be quite sufficient for most
terrestrial navigation and geodesy purposes, with UG_10 providing
additional integrational power for extended precision and more elliptic
bodies (e.g., Jupiter's cos{Oz} = .9351, Saturn's = .9020, and its moon
Prometheus = .6800 [outer facing] and .4595 [inner facing]), though, as
ellipticity increases, antipodal events become more of an issue.
While this response is undoubtably longer and more involved than was
sought, it will provide a rich resource for any future usenet inquiries
and searches for such information, as most of the different integral
models that have been documented here are the ones used in most of the
popular formularies in circulation today--both on and off line--which
are usually given without any identification of the actual integrals
(just identified as "internal variables" or something similarly vague)
or the coefficient structures to extend the various series out further.
~Kaimbridge M. GoldChild~