// ПРЯМАЯ И ОБРАТНАЯ ЗАДАЧИ
// "внешние" переменные заполняемые функцией CalcDelta (см.ниже)
// для анализа в обратной геодезической задаче
function _CalcDelta(Ellps: TEllipsoidProperty;P1:T3dpoint;Dist,Azim:extended;IsCalcK:boolean;var P0: T3dPoint) : T3dPoint;
var T2, N2, U2, v2, nu,_v, _u, L,
v_c, vc2, vc3, t4,
V, cosp1, t : extended;
FCTask : extended;
begin
FCTask:=sqrt(Ellps.one_es)/Ellps.a;
t:=tan(P1.X); t2:=sqr(t); t4:=sqr(t2);
nu:=sqrt(Ellps.es/Ellps.one_es)*cos(P1.X);
N2:=nu*nu;
_v:=Dist*sin(Azim); v2:=sqr(_v);
_u:=Dist*cos(Azim); U2:=sqr(_u);
cosp1:=1/cos(P1.X);
V:=sqrt(1+N2);
v_c:=V*FCTask; vc2:=sqr(v_c); vc3:=power(v_c,3);
// "внешние" переменные для анализа в обратной геодезической задаче
// b1
P0.B := v_c*(1+N2);
// l1
P0.L :=v_c*cosp1;
if iscalcK then exit;
//L5
L:=-P0.L*vc3*t*(1+3*t2+N2)/3;
//dB
result.X:=_u-0.5*FCTask*(V*sqr(_v)*t+(1+N2)*U2*N2*FCTask *(3*t/FCTask-t2*_u+_u))+
FCTask*FCTask/6*((1+N2)*v2*(1+3*T2+N2-9*N2*T2)*(0.25*(1+N2)*t*v2*FCTask-_u)+
P0.B*t*U2*3*(N2*U2*FCTask-(4-13*N2+3*T2*(2-3*N2))*_v/6) +
P0.B*V*v2*_u*((1+15*t2*(2+3*t2)) *v2/4 - (2+15*t2*(1+t2))*U2)*FCTask/5);
result.X:= result.X*P0.B+L*N2;
// dL
result.Y:=_v*(1+v_c*t*_u)+vc2*((1+3*t2+N2)*U2*_v-t2*V2*_v+v_c*t*(2+3*t2+N2)*U2*_u*_v+
0.2*vc2*_v*(t2*(1+3*t2)*V2*V2+(2+15*t2+15*t2*t2) *U2*U2-(1+20*t2+30*t2*t2)*V2*U2))/3;
result.Y:=result.Y*P0.L+L*(V2*_v*_u+N2);
// dA
result.Z:=v_c*_v*(t+ 0.5*v_c*(1+2*t2+N2)*_u)+
vc3*t*((5+6*t2+N2-4*N2*N2) *U2*_v -(1+2*t2+N2)*V2*_v +
v_c *((1.25+7*t2+6*t4+1.5*N2+2*N2*T2) *U2*_u*_v - (0.25+5*t2+6*t4+0.5*N2+2*N2*T2)*(V2*_v*_u+N2)) +
vc2*t*((1+20*t2+24*t4)*V2*V2*_v -(58+280*t2+240*t4) *V2*_v*U2 +(61+180*t2+120*t4) *_v*U2*U2)/20)/6;
end;
// ПРЯМАЯ ЗАДАЧА
function TrueGeoTask(Ellps: TEllpsType;const P1: T3dPoint; Azimuth,Dist: double): T3dPoint;
var P0,D : T3dpoint;
E : TEllipsoidProperty;
begin
result:=P1;
// если дистанция долее половини екватора Земли
if IsNan(Dist) or (Dist>2e7) or (Azimuth>1e2) then exit;
E:=Ellipsoids[Ellps];
D:=_CalcDelta(E,P1,Dist,Azimuth,false,P0);
if cos(P1.X)=0 then D.Y:=Azimuth;
result:=Set3dPoint(adjlon(P1.X+D.X,true),p1.Y*Byte(cos(P1.X)<>0)+D.Y,0);
if abs(P1.X+D.X)>PI/2 then result.L:=result.L+pi;
result.L:=adjlon(result.L);
end;
// ОБРАТНАЯ ЗАДАЧА
// X - dist Y-A12 z-A21
function ReversGeoTask(Ellps: TEllpsType;const P1,P2 : T3dPoint): t3dPoint;
var DeltaL,DeltaB,MinB,
MinL,aLen,dA12 : double;
D,D0,D1,P0 : T3dPoint;
Delta : T3dPoint;
E : TEllipsoidProperty;
I : byte;
begin
E:=Ellipsoids[Ellps];
D0.X:=P2.X-P1.X;
D0.Y:=P2.Y-P1.Y;
MinL:=MaxDouble;
MinB:=MinL;
D.X:=D0.X; D.Y:=D0.Y;
result.Y:=0;
i:=0; P0:=Set3dPoint(0,0,0);
dA12:=0; aLen:=0;
while i<254 do
begin
_CalcDelta(E,P1, aLen, dA12, true,P0);
if D.X=0 then D.X:=D.X+1e-10;
dA12:=Arctan(P0.b*D.Y/(P0.l*D.X));
if cos(dA12)=0 then break;
// aLen:=D.X/(P0.b*cos(dA12));
aLen:=Abs(D.X/(P0.b*cos(dA12))); // <<< + Abs исправил 27/06/2017
inc(i);
D1:=_CalcDelta(E,P1, aLen, dA12,false,P0);
DeltaB:=Abs(D1.X-D0.X);
DeltaL:=Abs(D1.Y-D0.Y);
if (DeltaB<MinB) and (DeltaL<MinL) then
begin
MinB:=DeltaB;
MinL:=DeltaL;
result.X:=aLen;
result.Y:=dA12;
end;
D.X:=D.X+0.5*DeltaB*(1-2*byte(D.X<D1.X));
D.Y:=D.Y+0.5*DeltaL*(1-2*byte(D.Y<D1.Y));
if Abs(D.X)>Abs(5*D0.X) then Break;
end;
if result.X<0 then result.Y:=-result.Y;
result.X:=Abs(result.X);
if d0.x<=0 then result.y:=pi-result.y;
if (d0.x>0) and (result.y<0) then result.y:=2*pi+result.y;
end;
где типы
TEllpsType =
(ellNone,ellMerit, ellSGS85, ellGRS80, ellIAU76, ellAiry, ellAPL49,
ellNWL9D, ellModAriy, ellAndrae, ellAustSA, ellGRS67,
ellBessel, ellBesselNamibia, ellClark66, ellClark80,
ellCPM, ellDelambre, ellEngelis, ellEverest30, ellEverest48,
ellEverest56, ellEverest69, ellEverestSS, ellFischer60,
ellFischer60m, ellFischer68, ellHelmert, ellHough,
ellInternational, ellKrass, ellKaula, ellLerch,
ellMaupertius, ellNewInternational, ellPlessis, ellSEAsia,
ellWalbeck, ellWGS60, ellWGS66, ellWGS72, ellWGS84, ellSphere
);
TDatumParam = array of double;
TDatum = packed record
eDest : TEllpsType;
Value : TDatumParam;
end;
TEllipsoidProperty = packed record
eName : string;
a : double; // major axis or radius if es=0
b : double; // minor axis or radius if es=0
e : double; // eccentricity
es : double; // e ^ 2
one_es : double; // 1 - e^2
rone_es : double; // 1/(1 - e^2)
ra : double; // 1/A
Rf : double; // references
Code : string;
end;
КОНСТАНТЫ ЭЛЛИПСОИДОВ
Ellipsoids : array[TEllpsType] of TEllipsoidProperty = (
(ename:'ПУСТОЙ СФЕРОИД';a:6000000;b:6000000;e: 0;es: 0;one_es: 1;rone_es: 1;ra: 1.66666E-7;rf: 0;code:'Empty'),
(eName:'MERIT 1983';a:6378137.0;b:6356752.29821597;e: 0.08181922;es: 0.00669438499958719;one_es: 0.993305615000413;rone_es: 1.00673950181947;
ra: 0.15678559428874E-6;rf: 0.00335281317789691;code:'MERIT'),
(ename:'Soviet Geodetic System 85';a:6378136.0;b:6356751.30156878;e: 0.08181922;
es: 0.0066943849995883;one_es: 0.993305615000412;rone_es: 1.00673950181947;ra: 0.156785618870466E-6;rf: 0.00335281317789691;code:'SGS85'),
(ename:'GRS 1980(IUGG 1980)';a:6378137.0;b:6356752.31414036;e: 0.08181919;
es: 0.00669438002289957;one_es: 0.9933056199771;rone_es: 1.00673949677548;ra: 0.15678559428874E-6;rf: 0.00335281068118232;code:'GRS80'),
(ename:'IAU 1976';a:6378140.0;b:6356755.28815753;e: 0.08181922;
es: 0.00669438499958741;one_es: 0.993305615000413;rone_es: 1.00673950181947;ra: 0.156785520543607E-6;rf: 0.00335281317789691;code:'IAU76'),
(ename:'Airy 1830';a:6377563.396;b:6356256.91;e: 0.08167337;
es: 0.00667053976159737;one_es: 0.993329460238403;rone_es: 1.00671533466852;ra: 0.156799695731319E-6;rf: 0.00334085052190358;code:'airy'),
(ename:'Appl. Physics. 1965';a:6378137.0;b:6356751.79631182;e: 0.08182018;
es: 0.0066945418545874;one_es: 0.993305458145413;rone_es: 1.00673966079587;ra: 0.15678559428874E-6;rf: 0.00335289186923722;code:'APL4.9'),
(ename:'Naval Weapons Lab 1965';a:6378145.0;b:6356759.76948868;e: 0.08182018;
es: 0.00669454185458873;one_es: 0.993305458145411;rone_es: 1.00673966079587;ra: 0.156785397635206E-6;rf: 0.00335289186923722;code:'NWL9D'),
(ename:'Modified Airy';a:6377340.189;b:6356034.446;e: 0.08167338;
es: 0.00667054060589878;one_es: 0.993329459394101;rone_es: 1.0067153355242;ra: 0.156805183722966E-6;rf: 0.00334085094546921;code:'mod_airy'),
(ename:'Andrae 1876 (Den.Iclnd.)';a:6377104.43;b:6355847.41523333;e: 0.08158159;
es: 0.00665555555555652;one_es: 0.993344444444443;rone_es: 1.00670014876791;ra: 0.156810980747888E-6;rf: 0.00333333333333333;code:'andrae'),
(ename:'Australian Natl & S. Amer. 1969';a:6378160.0;b:6356774.71919531;e: 0.08182018;
es: 0.0066945418545864;one_es: 0.993305458145414;rone_es: 1.00673966079587;ra: 0.156785028911159E-6;rf: 0.00335289186923722;code:'aust_SA'),
(ename:'GRS 67(IUGG 1967)';a:6378160.0;b:6356774.51609071;e: 0.08182057;
es: 0.00669460532856936;one_es: 0.993305394671431;rone_es: 1.00673972512833;ra: 0.156785028911159E-6;rf: 0.00335292371299641;code:'GRS67'),
(ename:'Bessel 1841';a:6377397.155;b:6356078.96281819;e: 0.081696831215255833813584066738272;
es: 0.006674372230614;one_es: 0.993325627769386;rone_es: 1.00671921879917;ra: 0.156803783063123E-6;rf: 0.00334277318217481;code:'bessel'),
(ename:'Bessel 1841 (Namibia)';a:6377483.865;b:6356165.38296633;e: 0.08169683;
es: 0.00667437223180078;one_es: 0.993325627768199;rone_es: 1.00671921879917;ra: 0.156801651116369E-6;rf: 0.00334277318217481;code:'bess_nam'),
(ename:'Clarke 1866';a:6378206.4;b:6356583.8;e: 0.08227185;
es: 0.0067686579972912;one_es: 0.993231342002709;rone_es: 1.00681478494592;ra: 0.156783888335755E-6;rf: 0.00339007530392876;code:'clrk66'),
(ename:'Clarke 1880 mod.';a:6378249.145;b:6356514.96582849;e: 0.08248322;
es: 0.00680348119602181;one_es: 0.993196518803978;rone_es: 1.00685008562476;ra: 0.156782837619931E-6;rf: 0.00340754628384929;code:'clrk80'),
(ename:'Comm. des Poids et Mesures 1799';a:6375738.7;b:6356666.22191211;e: 0.07729088;
es: 0.00597388071841887;one_es: 0.994026119281581;rone_es: 1.00600978244187;ra: 0.156844570810281E-6;rf: 0.00299141463998325;code:'CPM'),
(ename:'Delambre 1810 (Belgium)';a:6376428;b:6355957.92616372;e: 0.08006397;
es: 0.00641023989446932;one_es: 0.993589760105531;rone_es: 1.00645159617364;ra: 0.15682761571212E-6;rf: 0.00321027287319422;code:'delmbr'),
(ename:'Engelis 1985';a:6378136.05;b:6356751.32272154;e: 0.08181928;
es: 0.00669439396253357;one_es: 0.993305606037466;rone_es: 1.00673951090364;ra: 0.15678561764138E-6;rf: 0.00335281767444543;code:'engelis'),
(ename:'Everest 1830';a:6377276.345;b:6356075.41314024;e: 0.08147298;
es: 0.00663784663019951;one_es: 0.9933621533698;rone_es: 1.00668220206264;ra: 0.156806753526376E-6;rf: 0.00332444929666288;code:'evrst30'),
(ename:'Everest 1948';a:6377304.063;b:6356103.03899315;e: 0.08147298;
es: 0.00663784663020128;one_es: 0.993362153369799;rone_es: 1.00668220206264;ra: 0.156806071989232E-6;rf: 0.00332444929666288;code:'evrst48'),
(ename:'Everest 1956';a:6377301.243;b:6356100.2283681;e: 0.08147298;
es: 0.00663784663020017;one_es: 0.9933621533698;rone_es: 1.00668220206264;ra: 0.156806141327829E-6;rf: 0.00332444929666288;code:'evrst56'),
(ename:'Everest 1969';a:6377295.664;b:6356094.6679152;e: 0.08147298;
es: 0.00663784663020106;one_es: 0.993362153369799;rone_es: 1.00668220206264;ra: 0.156806278505327E-6;rf: 0.00332444929666288;code:'evrst69'),
(ename:'Everest (Sabah & Sarawak)';a:6377298.556;b:6356097.5503009;e: 0.08147298;
es: 0.00663784663019851;one_es: 0.993362153369801;rone_es: 1.00668220206264;ra: 0.156806207396259E-6;rf: 0.00332444929666288;code:'evrstSS'),
(ename:'Fischer (Mercury Datum) 1960';a:6378166;b:6356784.28360711;e: 0.08181333;
es: 0.00669342162296482;one_es: 0.993306578377035;rone_es: 1.00673852541468;ra: 0.156784881422026E-6;rf: 0.00335232986925913;code:'fschr60'),
(ename:'Modified Fischer 1960';a:6378155;b:6356773.32048274;e: 0.08181333;
es: 0.00669342162296449;one_es: 0.993306578377036;rone_es: 1.00673852541468;ra: 0.156785151818982E-6;rf: 0.00335232986925913;code:'fschr60m'),
(ename:'Fischer 1968';a:6378150;b:6356768.33724438;e: 0.08181333;
es: 0.00669342162296749;one_es: 0.993306578377033;rone_es: 1.00673852541468;ra: 0.156785274726998E-6;rf: 0.00335232986925913;code:'fschr68'),
(ename:'Helmert 1906';a:6378200;b:6356818.16962789;e: 0.08181333;
es: 0.00669342162296627;one_es: 0.993306578377034;rone_es: 1.00673852541468;ra: 0.156784045655514E-6;rf: 0.00335232986925913;code:'helmert'),
(ename:'Hough';a:6378270.0;b:6356794.34343434;e: 0.08199189;
es: 0.00672267002233429;one_es: 0.993277329977666;rone_es: 1.00676817019722;ra: 0.15678232498781E-6;rf: 0.00336700336700337;code:'hough'),
(ename:'International 1909 (Hayford)';a:6378388.0;b:6356911.94612795;e: 0.08199189;
es: 0.00672267002233207;one_es: 0.993277329977668;rone_es: 1.00676817019722;ra: 0.156779424519173E-6;rf: 0.00336700336700337;code:'intl'),
(ename:'(CK-42) Krassovsky 1942';a:6378245;b:6356863.01877305;e: 0.081813334;
es: 0.00669342162296504;one_es: 0.993306578377035;rone_es: 1.00673852541468;ra: 0.156782939507655E-6;rf: 0.00335232986925913;code:'krass'),
(ename:'Kaula 1961';a:6378163;b:6356776.99208691;e: 0.08182155;
es: 0.0066947659459099;one_es: 0.99330523405409;rone_es: 1.00673988791802;ra: 0.156784955166558E-6;rf: 0.00335300429184549;code:'kaula'),
(ename:'Lerch 1979';a:6378139;b:6356754.29151034;e: 0.08181922;
es: 0.00669438499958852;one_es: 0.993305615000411;rone_es: 1.00673950181947;ra: 0.15678554512531E-6;rf: 0.00335281317789691;code:'lerch'),
(ename:'Maupertius 1738';a:6397300;b:6363806.28272251;e: 0.10219488;
es: 0.0104437926591934;one_es: 0.989556207340807;rone_es: 1.0105540166205;ra: 0.15631594578963E-6;rf: 0.00523560209424084;code:'mprts'),
(ename:'New International 1967';a:6378157.5;b:6356772.2;e: 0.08182023;
es: 0.0066945504730862;one_es: 0.993305449526914;rone_es: 1.00673966953093;ra: 0.156785090365047E-6;rf: 0.00335289619298362;code:'new_intl'),
(ename:'Plessis 1817 (France)';a:6376523;b:6355863;e: 0.08043334;
es: 0.00646952287129587;one_es: 0.993530477128704;rone_es: 1.00651165014081;ra: 0.15682527923133E-6;rf: 0.00324001026891929;code:'plessis'),
(ename:'Southeast Asia';a:6378155.0;b:6356773.3205;e: 0.08181333;
es: 0.00669342161757036;one_es: 0.99330657838243;rone_es: 1.00673852540921;ra: 0.156785151818982E-6;rf: 0.00335232986655221;code:''),
(ename:'Walbeck';a:6376896.0;b:6355834.8467;e: 0.08120682;
es: 0.00659454809019966;one_es: 0.9934054519098;rone_es: 1.00663832484261;ra: 0.156816106143177E-6;rf: 0.00330272805139054;code:'Walbeck'),
(ename:'WGS 60';a:6378165.0;b:6356783.28695944;e: 0.08181333;es: 0.00669342162296482;one_es: 0.993306578377035;rone_es: 1.00673852541468;ra: 0.156784906003529E-6;rf: 0.00335232986925913;code:'WGS60'),
(ename:'WGS 66';a:6378145.0;b:6356759.76948868;e: 0.08182018;es: 0.00669454185458873;one_es: 0.993305458145411;rone_es: 1.00673966079587;ra: 0.156785397635206E-6;rf: 0.00335289186923722;code:'WGS66'),
(ename:'WGS 72';a:6378135.0;b:6356750.52001609;e: 0.08181881;es: 0.00669431777826779;one_es: 0.993305682221732;rone_es: 1.00673943368903;ra: 0.1567856434522E-6;rf: 0.0033527794541675;code:'WGS72'),
(ename:'WGS 84';a:6378137.0;b:6356752.31424518;e: 0.08181919;es: 0.00669437999014111;
one_es: 0.99330562000985889;rone_es:1.0067394967422762251591434067861;ra:0.15678559428874E-6;
rf: 0.00335281066474748;code:'WGS84'),
(ename:'Normal Sphere (r=6370997)';a:6370997.0;b:6370997;e: 0;es: 0;one_es: 1;rone_es: 1;ra: 0.156961304486566E-6;rf: 0;code:'sphere'),
|