+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP Prime Software Library (/forum-15.html)
+--- Thread: Geodesy API (/thread-7062.html)
Geodesy API - StephenG1CMZ - 10-17-2016 09:49 PM
A collection of useful geodesy routines.
Disclaimer: Not for navigational use. No liability accepted.
RE: Geodesy API - StephenG1CMZ - 10-17-2016 10:11 PM
Anyone using Version 0.1 or Version 0.2 should now replace it.
RE: Geodesy API - StephenG1CMZ - 10-23-2016 06:12 PM
Version 0.3 implements LatLong-to-distance calculations.
Spherical calculations use Haversine.
Ellipsoidal calculations use Vincenty.
These routines have not been tested apart from a handful of values.
Code:
// GEODESY © 2016 StephenG1CMZ
LOCAL VER:="0.3";
LOCAL NL:=CHAR(10);
LOCAL SP:=" ";
G_Use(); //FORWARD
G_RadiusB(); //FORWARD
G_MeanRadius(); //FORWARD
//We assume RA≥RB also Trig in Radians
//Geoid Data
EXPORT WGS84a_m:=6378137.0;
EXPORT WGS84a_km:=WGS84a_m/1000;
LOCAL WGS84fINV:=298.257223563;
EXPORT WGS84f:=1/WGS84fINV;
LOCAL WGS84b_m:=G_RadiusB(WGS84a_m,WGS84f);
LOCAL WGS84b_km:=WGS84b_m/1000;
LOCAL RA:=WGS84a_m; //RADIUS EQ
LOCAL FLATTEN:=WGS84f;
LOCAL RB:=G_RadiusB(RA,FLATTEN);
LOCAL RR:=G_MeanRadius(RA,RB);
LOCAL EPSILON_V:=3ᴇ−11;
LOCAL LoopSafe:=140;
//FORMAT:
//RA RB F
EXPORT G_XGNAME:=1;
EXPORT G_XRA:=2;
EXPORT G_XF:=4;
EXPORT G_EarthABF:={
{"WGS84_m", WGS84a_m,0,WGS84f},
{"WGS84_km", WGS84a_m/1000,0,WGS84f},
{"Bessel_m", 6377397.155,0,1/299.1528128},
{"International_m",6378388.000,0,1/297},
{"Sphere", 1,0,0}};
EXPORT G_ThisEll:=G_EarthABF(1);
EXPORT G_XLAT:=1;
EXPORT G_XLON:=2;
LOCAL ESS_UNITS:="m";
EXPORT G_About()
BEGIN
PRINT();
PRINT("Geodesy API V 0.3 © 2016 StephenG1CMZ");
PRINT("NOT FOR NAVIGATIONAL USE");
PRINT("No liability accepted");
PRINT("Expect the ordering of Vincenty results to change in later versions");
END;
// STD MATHS
D2R(DEGS)
BEGIN
RETURN DEGS*(π/180);
END;
R2D(RADS)
BEGIN
RETURN RADS*(180/π);
END;
ATAN2YX(YY,XX)
BEGIN //MY ATAN2
IF XX==0 AND YY==0 THEN
RETURN 0; //CONVENIENT;
END;
//-π..π THE SPECIAL CHAR BELOW IS i=sqrt(−1), NOT CHINESE
RETURN ARG(XX+*YY);//BUILT-IN HANDLES MOST. END;
END;
HAV(THETA)
//TRIG: HAVERSINE=SIN^2(THETA/2)
BEGIN
LOCAL SN:=SIN(THETA/2);
RETURN SN^2;
END;
// END STD MATHS
EXPORT G_RadiusB(RadiusA,Flatten)
//RETURN POLAR RADIUS
BEGIN
RETURN (1-Flatten)*RadiusA;//RADIUS POLAR
END;
EXPORT G_Use(ST)
//A RADIUS EQU
//B RADIUS POL
//AT PRESENT RB IS ALWAYS DERIVED FROM F
//LATER ALLOW F TO BE DERIVED
BEGIN
LOCAL II,POS;
FOR II FROM 1 TO SIZE(G_EarthABF) DO
POS:=INSTRING(G_EarthABF(G_XGNAME),ST);
IF POS THEN
G_ThisEll:=G_EarthABF;
RA:=G_ThisEll(G_XRA);
FLATTEN:=G_ThisEll(G_XF);
RB:=G_RadiusB(RA,FLATTEN);
RETURN POS;
END;
END;
RETURN POS;//0=>UNCHANGED
END;
EXPORT G_DxVincentyDo(RA,RB,FLATTEN,Lat1,Lon1,Lat2,Lon2)
// GEODESY INVERSE PROBLEM ELLIPSOIDAL
// USE: RA AND FLATTEN FROM ThisEll
// RETURNS DIST IN SAME UNITS AS R
// LAT AND LON IN DEGREES
BEGIN
LOCAL NUM,DEN;
LOCAL LAT1:=HMS→(Lat1);
LOCAL LAT2:=HMS→(Lat2);
LOCAL LON1:=HMS→(Lon1);
LOCAL LON2:=HMS→(Lon2);
LOCAL UU1,SINUU1,COSUU1; //U1 AUX LATITUDE
LOCAL UU2,SINUU2,COSUU2; //U2 AUX LATITUDE
LOCAL USQ;
LOCAL AAA,BBB;
LOCAL LONDIFFradians; //L=DIFF IN LON
LOCAL PHI1,PHI2;
LOCAL LAMBDA,SINLAMBDA,COSLAMBDA,LAMBDAOLD;
LOCAL SIGMA,SIN_SIGMA,COS_SIGMA,DELTA_SIGMA,SINSQ_SIGMA; //AUX ARC LENGTH
LOCAL SIN_ALPHA,COSSQ_ALPHA;
LOCAL COS_2SIGMAM,COSSQ_2SIGMAM;
LOCAL CC;
LOCAL BRACKET_TERM;
LOCAL LOOPCOUNT;
LOCAL DISTS;
LOCAL AZ_FWD,AZ_FNL,AZ_REV;
// CALC L
LONDIFFradians:=D2R(LON2-LON1); //DIFF IN LON
// CALC U1 U2
PHI1:=D2R(LAT1);
PHI2:=D2R(LAT2);
UU1:=ATAN((1-FLATTEN)*TAN(PHI1));//CORRECT
UU2:=ATAN((1-FLATTEN)*TAN(PHI2));//CORRECT
//TO CALC LAMBDA:
// PRE-CALC U1 TRIG ONCEONLY
SINUU1:=SIN(UU1); COSUU1:=COS(UU1);
// PRE-CALC U2 TRIG EACH POSITION
SINUU2:=SIN(UU2); COSUU2:=COS(UU2);
//ITERATE CALCULATING LAMBDA
LOOPCOUNT:=0;
LAMBDA:=LONDIFFradians; //V13: DIFF IN LON ON AUX. FIRST APPROX
REPEAT
LOOPCOUNT:=LOOPCOUNT+1;
SINLAMBDA:=SIN(LAMBDA); COSLAMBDA:=COS(LAMBDA); //PRECALC
//V14
BRACKET_TERM:=(COSUU1*SINUU2 - SINUU1*COSUU2*COSLAMBDA);
SINSQ_SIGMA:=(COSUU2*SINLAMBDA)^2 + (BRACKET_TERM)^2;
IF SINSQ_SIGMA==0 THEN
//PRINT("SINσ 0: NO DISTANCE");
//0 DISTANCE MIGHT SUGGEST UNDERFLOW
RETURN {1,0,0,0,0,0};//VALID ZERO
END;
SIN_SIGMA:=√(SINSQ_SIGMA);
//V15
COS_SIGMA:=SINUU1*SINUU2 + COSUU1*COSUU2*COSLAMBDA;
//V16
SIGMA:=ATAN2YX(SIN_SIGMA,COS_SIGMA);
//V17
SIN_ALPHA:=(COSUU1*COSUU2*SINLAMBDA)/SIN_SIGMA;
COSSQ_ALPHA:=1-(SIN_ALPHA*SIN_ALPHA);
//V18
IF COSSQ_ALPHA==0 THEN
//PRINT("AVOID COSSQα 0");
COS_2SIGMAM:=0;
ELSE
COS_2SIGMAM:=COS_SIGMA - (2*SINUU1*SINUU2/COSSQ_ALPHA);
END;
CC:=(FLATTEN/16)*COSSQ_ALPHA*(4 + FLATTEN*(4 - 3*COSSQ_ALPHA));
LAMBDAOLD:=LAMBDA;
LAMBDA:=LONDIFFradians + (1-CC)*FLATTEN*SIN_ALPHA*(SIGMA + CC*SIN_SIGMA*(COS_2SIGMAM + CC*COS_SIGMA*(−1+2*(COS_2SIGMAM*COS_2SIGMAM))));
//PRINT(LAMBDA);
UNTIL (ABS(LAMBDA-LAMBDAOLD)≤EPSILON_V) OR (LOOPCOUNT≥LoopSafe);
//V19
USQ:=COSSQ_ALPHA*((RA^2 - RB^2)/(RB^2));
//(V3)
AAA:=1 + (USQ/16384)*(4096 + USQ*(−768 + USQ*(320 - 175*USQ)));
//(V4)
BBB:= (USQ/ 1024)*( 256 + USQ*(−128 + USQ*( 74 - 47*USQ)));
//(V6)
COSSQ_2SIGMAM:=(COS_2SIGMAM)^2;
DELTA_SIGMA:=(1/6)*BBB*COS_2SIGMAM*(−3+4*(SIN_ALPHA)^2)*(−3+4*COSSQ_2SIGMAM) ;
DELTA_SIGMA:=(1/4)*BBB*(COS_SIGMA*(−1+2*COSSQ_2SIGMAM)) - DELTA_SIGMA ;
DELTA_SIGMA:=BBB*SIN_SIGMA*(COS_2SIGMAM + DELTA_SIGMA);
//PRINT("Δσ "+DELTA_SIGMA);
DISTS:=AAA*(SIGMA - DELTA_SIGMA);//ANGULAR ≤π Radians
DISTS:=RB*DISTS; //same units as radius
//V20:FORWARD AZIMUTH
//CHECKED SINLAMBDA:ORIGINAL (LON2-LON1) OR LOOPED VALUE? USE LOOPED
AZ_FWD:=ATAN2YX((COSUU2*SINLAMBDA),BRACKET_TERM);
//V21:FINAL AND REVERSE AZIMUTH
NUM:=(COSUU1*SINLAMBDA);
DEN:=(−SINUU1*COSUU2 + COSUU1*SINUU2*COSLAMBDA);
AZ_FNL:=ATAN2YX(NUM,DEN);
AZ_REV:=ATAN2YX(−NUM,−DEN);
//RETURN RESULTS
DISTS:=ROUND(DISTS,3);//VINCENTY CLAIMS mm,
AZ_FWD:=R2D(AZ_FWD);
AZ_FWD:=IFTE(AZ_FWD<0,((AZ_FWD+360) MOD 360),AZ_FWD);//NO NEGATIVE BEARINGS
AZ_FNL:=R2D(AZ_FNL);
AZ_FNL:=IFTE(AZ_FNL<0,((AZ_FNL+360) MOD 360),AZ_FNL);//NO NEGATIVE BEARINGS
AZ_REV:=R2D(AZ_REV);
AZ_REV:=IFTE(AZ_REV<0,((AZ_REV+360) MOD 360),AZ_REV);
// VALIDITY: CONVERGENCE FAILURE
// DISTS: SAME UNITS AS RADIUS
// AZIMUTHS: POSITIVE 0..360
// LOOPCOUNT: CONVERGENCE INFO (LOWER BETTER):DIAGNOSTIC ONLY
// EXPECT SEQUENCE ORDER TO CHANGE
// ESPECIALLY AZIMUTHS
RETURN {(LOOPCOUNT≠LoopSafe),DISTS,(AZ_FWD),(AZ_FNL),(AZ_REV),LOOPCOUNT};
END;
EXPORT G_DxVincenty(Elevation,Point1,Point2)
// GEODESY INVERSE PROBLEM ELLIPSOIDAL
// USE: RA AND FLATTEN FROM ThisEll TO DERIVE RB
// RETURNS DIST IN SAME UNITS AS R
// LAT AND LON IN DEGREES
BEGIN
LOCAL LAT1:=Point1(G_XLAT);
LOCAL LON1:=Point1(G_XLON);
LOCAL LAT2:=Point2(G_XLAT);
LOCAL LON2:=Point2(G_XLON);
LOCAL RA:=G_ThisEll(G_XRA)+Elevation;
LOCAL FLATTEN:=G_ThisEll(G_XF);
LOCAL RB:=G_RadiusB(RA,FLATTEN);
RETURN G_DxVincentyDo(RA,RB,FLATTEN,LAT1,LON1,LAT2,LON2);
END;
// Main Spherical Geodesy routines-
EXPORT G_MeanRadius (RA,RB)
//In Geodesy, this is usually used as a radius for spherical
BEGIN
RETURN (2*RA+RB)/3;
END;
EXPORT G_DistHav(Point1,Point2)
BEGIN
LOCAL LAT1:=HMS→(Point1(G_XLAT));
LOCAL LAT2:=HMS→(Point2(G_XLAT));
LOCAL LON1:=HMS→(Point1(G_XLON));
LOCAL LON2:=HMS→(Point2(G_XLON));
LOCAL DLON:=D2R(LON2-LON1);
LOCAL DLAT:=D2R(LAT2-LAT1);
LOCAL AA:=HAV(DLAT) + (COS(D2R(LAT1))*COS(D2R(LAT2))*HAV(DLON));
LOCAL CC:=2*ASIN(MIN(1,√(AA)));//MIN 1 AVOIDS ROUNDING OUTRANGE
LOCAL DD:=RR*CC;
RETURN ROUND(DD,−6);
END;
// End Main Spherical Routines
G_EllStats(PR)
//PR:PRINT ELSE TEXTOUT
BEGIN
LOCAL II;
LOCAL ST:={};
ST(1):=G_ThisEll(G_XGNAME);
ST(2):="RA: "+RA;
ST(3):="RB: "+RB;
ST(4):="f: "+(FLATTEN);
IF FLATTEN≠0 THEN
ST(5):="1/f: "+(1/FLATTEN);
END;
IF PR THEN
FOR II FROM 1 TO SIZE(ST) DO
PRINT(ST(II));
END;
ELSE
FOR II FROM 1 TO SIZE(ST) DO
TEXTOUT_P(ST(II),0,20*II+100);
END;
END;
END;
G_EX(LAT1,LON1,LAT2,LON2)
//An example displaying results
BEGIN
LOCAL TM,RESULT1,RESULT2,DISTS;
PRINT(" ");
PRINT("From LAT "+HMS→(LAT1)+" "+→HMS(LAT1));
PRINT("From LON "+HMS→(LON1)+" "+→HMS(LON1));
PRINT("To LAT "+HMS→(LAT2)+" "+→HMS(LAT2));
PRINT("To LON "+HMS→(LON2)+" "+→HMS(LON2));
TM:=TICKS;
RESULT1:=G_DistHav({LAT1,LON1},{LAT2,LON2});
RESULT2:=G_DxVincenty(0,{LAT1,LON1},{LAT2,LON2});
TM:=TICKS-TM;
PRINT("Elapsed "+(TM/1000)+" s");
PRINT("Spherical Haversine: "+(RESULT1)+" m");
PRINT("Ellipsoidal Vincenty: "+(RESULT2(2))+" m");
PRINT(" AZ "+(→HMS(RESULT2(3))));
PRINT(" AZ "+(→HMS(RESULT2(4))));
PRINT(" AZ "+(→HMS(RESULT2(5))));
PRINT("Validity: " +(RESULT2(1))+" "+(RESULT2(6)));
//PRINT(RESULT2);
//RESULT:=Z_GEODESY(LAT1,LON1,LAT2,LON2);
END;
EXPORT G_NewEll(LST)
BEGIN
G_EarthABF(0):=LST;
RETURN 1;
END;
EXPORT ExampleUse()
//A short example of key procedures
BEGIN
LOCAL OKC,DIST,RESULT,EXAMPLEPOINT,XLAT,XLON;
OKC:=G_NewEll({"Add your own",{RA,0,FLATTEN}});//Insert into List. Use use to use it.
OKC:=G_Use("WGS84_m");//DEFAULT
IF OKC THEN //Found in list
DIST:=G_DistHav({1,2},{3,4});
//RESULT:=G_DxVincentyDo(RA,RB,FLATTEN,LAT1,LON1,LAT2,LON2);
RESULT:=G_DxVincenty(0,{1,2},{3,4});
//Variables prefixed X index into data structures
EXAMPLEPOINT:={"From Here",1,2};
XLAT:=2; XLON:=3;
XLAT:=1; XLON:=2; //REINSTATE CURRENT DEFAULT
EXAMPLEPOINT:={1,2};
//and similarly for geoid structure
END;
END;
EXPORT G_DrawEll(Flatten)
//The stats are omitted because flattening drawn may differ
//makes illustrate changing f easier
//for now use myRadiusA of 100 pixels
BEGIN
LOCAL XX,YY;
LOCAL DX,DY;
LOCAL MyRadiusA:=100;
LOCAL HEIGHT:=100;
LOCAL WIDTH:=100;
//LOCAL SCALE:=RadiusA/WIDTH;
LOCAL EWIDTH:=MyRadiusA;
LOCAL EHEIGHT:=((1-Flatten)*MyRadiusA);//RADIUS POLAR
RECT();
//G_EllStats(0);
FOR YY FROM −HEIGHT TO HEIGHT DO
//PRINT(YY);
FOR XX FROM −WIDTH TO WIDTH DO
DX:=XX/(EWIDTH);
DY:=YY/(EHEIGHT);
IF (DX^2+DY^2)≤1 THEN
//PRINT({WIDTH+10+XX,HEIGHT+20+YY});
PIXON_P(WIDTH+100+XX,HEIGHT+20+YY,#009900h);
END;//IF
END;//FOR
END;//FOR
//G_EllStats(0);
WAIT;
END;
G_TEST()
//A few tests. More needed.
BEGIN
G_EX(50°03′58.76″,−005°42′53.10″,58°38′38.48″,−003°04′12.34″);//50PERCENT
PRINT("EXPECT 969 954_m");
G_EX(−37°57′3.72030″,144°25′29.52440″,−37°39′10.15610″,143°55′35.38390″);
PRINT("EXPECTED: 54972.271m");//VINCENTY AUSTRALIA
PRINT("306d52 05 * 127d10 25 * 307d10 25.07 ");
G_EX(31.8300167,35.0662833,31.83000000,35.0708167);
PRINT("PYTHON: CIRCLE 428.4 VIN 429.1677 m");
//G_EX(0,0,0.5,179.5);
G_EX(0,0,0.5,179.7);//WILL NOT CONVERGE
//G_EX(0,0,0,90);
//G_EX(0,0,0.5,0.5);
G_EX(2,2,2,2);
G_EX(50°03′58.76″,−005°42′53.10″,58°38′38.48″,−003°04′12.34″);//50PERCENT
PRINT("EXPECT 969 954_m");
//IF 0 THEN
G_EX(50°3′59″,−5°42′53″,58°38′38″,−3°4′12″);
PRINT("RECHECK: EXPECTED 968.9km 009DEG0711");
IF 0 THEN
G_EX(−37°57′3.72030″,144°25′29.52440″,−37°39′10.15610″,143°55′35.38390″);//EXPECTS 54972m
PRINT("EXPECTED: 54972.271m ");
PRINT("306d52:5.37 / 127d10:25.07");//VINCENTY AUSTRALIA
//IF 1 THEN
G_EX(31.8300167,35.0662833,31.83000000,35.0708167);
PRINT("PYTHON: CIRCLE 428.4 VIN 429.1677 m");
//IF 0 THEN
G_EX(0,0,0.5,179.5);
G_EX(0,0,0.5,179.7);//WILL NOT CONVERGE
//IF 1 THEN
G_EX(0,0,0,90);
G_EX(0,0,0.5,0.5);
G_EX(2,2,2,2);
END;
//PRINT("SEMI: "+π*ESS(1));//MAX
END;
EXPORT GEODESY()
BEGIN
G_About();
G_EllStats(1);
PRINT("Esc");WAIT;
G_TEST();
END;Update: Version 0.3 Vincenty is deprecated as it has bugs detailed below.
RE: Geodesy API - StephenG1CMZ - 11-22-2016 06:34 PM
Version 0.3 Vincenty is deprecated as it has the following typo:
In the calculation of DELTA_SIGMA there is a reference to SIN_ALPHA which should read SIN_SIGMA.
RE: Geodesy API - StephenG1CMZ - 11-23-2016 02:09 PM
Version 0.31 fixes the typo in the inverse Vincenty (LatLon to distance).
It also adds more predefined ellipsoids (including the recently announced roundest object in nature) and a few routines.
Note that the Direct Vincenty routine (LatLon and distance to LatLon2) is included for debugging only It's results as implemented here are known to be incorrect by about 100 km, or worse.
Code:
// GEODESY © 2016 StephenG1CMZ
LOCAL VER:="0.31";
LOCAL NL:=CHAR(10);
LOCAL SP:=" ";
//FORWARD
G_Flattening();//Derive Flattening from radius RA RB
G_MeanRadius();//Derive mean radius from radius RA RB
G_RadiusB(); //Derive RB from RA and flattening
G_Use(); //Use this ellipsoid
//END FORWARD
//We assume RA≥RB also Trig in Radians
// Ell=ElIpsoid
// Sp=Spherical
//Geoid Data
EXPORT WGS84a_m:=6378137.0;
EXPORT WGS84a_km:=WGS84a_m/1000;
LOCAL WGS84fINV:=298.257223563;
EXPORT WGS84f:=1/WGS84fINV;
LOCAL WGS84b_m:=G_RadiusB(WGS84a_m,WGS84f);
LOCAL WGS84b_km:=WGS84b_m/1000;
LOCAL RA:=WGS84a_m; //RADIUS EQ
LOCAL FLATTEN:=WGS84f;
LOCAL RB:=G_RadiusB(RA,FLATTEN);
LOCAL RR:=G_MeanRadius(RA,RB);
LOCAL EPSILON_V:=1ᴇ−12;
LOCAL LoopSafe:=140;
//FORMAT:
//RA RB F
EXPORT G_XGNAME:=1;
EXPORT G_XRA:=2;
EXPORT G_XF:=4;
EXPORT G_EllABF:={
{"WGS84_m", WGS84a_m, 0, WGS84f},
{"Bessel_m", 6377397.155, 0, 1/299.1528128},
{"International_m", 6378388.000, 0, 1/297},
{"Sphere", 1,0,0},
{"Mercury_km", 2439.7, 2439.7, 0.0000},
{"Venus_km", 6051.8, 6051.8, 0.0000},
{"Earth WGS84_km", WGS84a_km, 0, WGS84f},
{"Moon_km", 1738.1, 1736.0, 0.0012},
{"Mars_km", 3396.2, 3376.2, 0.00589},
{"Jupiter_km", 71492, 66854, 0.06487},
{"Saturn_km", 60268, 54364, 0.09796},
{"Uranus_km", 22559, 24973, 0.02293},
{"Neptune_km", 24764, 24341, 0.01708},
{"Pluto_km", 1187, 1187, 0.0000},
{"Sun", 695508, 0, 0.00005},
{"Kepler 11145123 _km", 1.5ᴇ6, 1.5ᴇ6-3, G_Flattening(1.5ᴇ6,1.5ᴇ6-3)}};
// Kepler 11145123 is roundest object in nature as at 2016-11
EXPORT G_ThisEll:=G_EllABF(1);
EXPORT G_XLAT:=1;
EXPORT G_XLON:=2;
LOCAL ESS_UNITS:="m";
EXPORT G_About()
BEGIN
PRINT();
PRINT("Geodesy API V "+VER+" © 2016 StephenG1CMZ");
PRINT("NOT FOR NAVIGATIONAL USE");
PRINT("No liability accepted");
PRINT("Expect the ordering of Vincenty results to change in later versions");
END;
// STD MATHS
D2R(DEGS)
BEGIN
RETURN DEGS*(π/180);
END;
R2D(RADS)
BEGIN
RETURN RADS*(180/π);
END;
ATAN2YX(YY,XX)
BEGIN //MY ATAN2
IF XX==0 AND YY==0 THEN
RETURN 0; //CONVENIENT;
END;
//-π..π THE SPECIAL CHAR BELOW IS i=sqrt(−1), NOT CHINESE
RETURN ARG(XX+*YY);//BUILT-IN HANDLES MOST. END;
END;
HAV(THETA)
//TRIG: HAVERSINE=SIN^2(THETA/2)
BEGIN
LOCAL SN:=SIN(THETA/2);
RETURN SN^2;
END;
// END STD MATHS
EXPORT G_Flattening(RadiusA,RadiusB)
BEGIN
RETURN (RadiusA-RadiusB)/RadiusA;
END;
EXPORT G_RadiusB(RadiusA,Flatten)
//RETURN POLAR RADIUS
BEGIN
RETURN (1-Flatten)*RadiusA;//RADIUS POLAR
END;
EXPORT G_DxVincentyDo(RA,RB,FLATTEN,Lat1,Lon1,Lat2,Lon2)
// GEODESY INVERSE PROBLEM ELLIPSOIDAL
// USE: RA AND FLATTEN FROM ThisEll
// RETURNS DIST IN SAME UNITS AS R
// LAT AND LON IN DEGREES
BEGIN
LOCAL NUM,DEN;
LOCAL LAT1:=HMS→(Lat1);
LOCAL LAT2:=HMS→(Lat2);
LOCAL LON1:=HMS→(Lon1);
LOCAL LON2:=HMS→(Lon2);
LOCAL UU1,SINUU1,COSUU1; //U1 AUX LATITUDE
LOCAL UU2,SINUU2,COSUU2; //U2 AUX LATITUDE
LOCAL USQ;
LOCAL AAA,BBB;
LOCAL LONDIFFradians; //L=DIFF IN LON
LOCAL PHI1,PHI2;
LOCAL LAMBDA,SINLAMBDA,COSLAMBDA,LAMBDAOLD;
LOCAL SIGMA,SIN_SIGMA,COS_SIGMA,DELTA_SIGMA,SINSQ_SIGMA; //AUX ARC LENGTH
LOCAL SIN_ALPHA,COSSQ_ALPHA;
LOCAL COS_2SIGMAM,COSSQ_2SIGMAM;
LOCAL CC;
LOCAL BRACKET_TERM;
LOCAL LOOPCOUNT;
LOCAL DISTS;
LOCAL AZ_FWD,AZ_FNL,AZ_REV;
// CALC L
LONDIFFradians:=D2R(LON2-LON1); //DIFF IN LON
// CALC U1 U2
PHI1:=D2R(LAT1);
PHI2:=D2R(LAT2);
UU1:=ATAN((1-FLATTEN)*TAN(PHI1));//CORRECT
UU2:=ATAN((1-FLATTEN)*TAN(PHI2));//CORRECT
//TO CALC LAMBDA:
// PRE-CALC U1 TRIG ONCEONLY
SINUU1:=SIN(UU1); COSUU1:=COS(UU1);
// PRE-CALC U2 TRIG EACH POSITION
SINUU2:=SIN(UU2); COSUU2:=COS(UU2);
//ITERATE CALCULATING LAMBDA
LOOPCOUNT:=0;
LAMBDA:=LONDIFFradians; //V13: DIFF IN LON ON AUX. FIRST APPROX
REPEAT
LOOPCOUNT:=LOOPCOUNT+1;
SINLAMBDA:=SIN(LAMBDA); COSLAMBDA:=COS(LAMBDA); //PRECALC
//V14
BRACKET_TERM:=(COSUU1*SINUU2 - SINUU1*COSUU2*COSLAMBDA);
SINSQ_SIGMA:=(COSUU2*SINLAMBDA)^2 + (BRACKET_TERM)^2;
IF SINSQ_SIGMA==0 THEN
//PRINT("SINσ 0: NO DISTANCE");
//0 DISTANCE MIGHT SUGGEST UNDERFLOW
RETURN {1,0,0,0,0,0};//VALID ZERO
END;
SIN_SIGMA:=√(SINSQ_SIGMA);
//V15
COS_SIGMA:=SINUU1*SINUU2 + COSUU1*COSUU2*COSLAMBDA;
//V16
SIGMA:=ATAN2YX(SIN_SIGMA,COS_SIGMA);
//V17
SIN_ALPHA:=(COSUU1*COSUU2*SINLAMBDA)/SIN_SIGMA;
COSSQ_ALPHA:=1-(SIN_ALPHA*SIN_ALPHA);
//V18
IF COSSQ_ALPHA==0 THEN
//PRINT("AVOID COSSQα 0");
COS_2SIGMAM:=0;
ELSE
COS_2SIGMAM:=COS_SIGMA - (2*SINUU1*SINUU2/COSSQ_ALPHA);
END;
CC:=(FLATTEN/16)*COSSQ_ALPHA*(4 + FLATTEN*(4 - 3*COSSQ_ALPHA));
LAMBDAOLD:=LAMBDA;
LAMBDA:=LONDIFFradians + (1-CC)*FLATTEN*SIN_ALPHA*(SIGMA + CC*SIN_SIGMA*(COS_2SIGMAM + CC*COS_SIGMA*(−1+2*(COS_2SIGMAM*COS_2SIGMAM))));
//PRINT(LAMBDA);
UNTIL (ABS(LAMBDA-LAMBDAOLD)≤EPSILON_V) OR (LOOPCOUNT≥LoopSafe);
//V19
USQ:=COSSQ_ALPHA*((RA^2 - RB^2)/(RB^2));
//(V3)
AAA:=1 + (USQ/16384)*(4096 + USQ*(−768 + USQ*(320 - 175*USQ)));
//(V4)
BBB:= (USQ/ 1024)*( 256 + USQ*(−128 + USQ*( 74 - 47*USQ)));
//(V6)
COSSQ_2SIGMAM:=(COS_2SIGMAM)^2;
//DELTA_SIGMA:=(1/6)*BBB*COS_2SIGMAM*(−3+4*(SIN_ALPHA)^2)*(−3+4*COSSQ_2SIGMAM) ;//V0.3
DELTA_SIGMA:=(1/6)*BBB*COS_2SIGMAM*(−3+4*(SIN_SIGMA)^2)*(−3+4*COSSQ_2SIGMAM) ; //V0.4
DELTA_SIGMA:=(1/4)*BBB*(COS_SIGMA*(−1+2*COSSQ_2SIGMAM)) - DELTA_SIGMA ;
DELTA_SIGMA:=BBB*SIN_SIGMA*(COS_2SIGMAM + DELTA_SIGMA);
//PRINT("Δσ "+DELTA_SIGMA);
DISTS:=AAA*(SIGMA - DELTA_SIGMA);//ANGULAR ≤π Radians
DISTS:=RB*DISTS; //same units as radius
//V20:FORWARD AZIMUTH
//CHECKED SINLAMBDA:ORIGINAL (LON2-LON1) OR LOOPED VALUE? USE LOOPED
AZ_FWD:=ATAN2YX((COSUU2*SINLAMBDA),BRACKET_TERM);
//V21:FINAL AND REVERSE AZIMUTH
NUM:=(COSUU1*SINLAMBDA);
DEN:=(−SINUU1*COSUU2 + COSUU1*SINUU2*COSLAMBDA);
AZ_FNL:=ATAN2YX(NUM,DEN);
AZ_REV:=ATAN2YX(−NUM,−DEN);
//RETURN RESULTS
DISTS:=ROUND(DISTS,3);//VINCENTY CLAIMS mm,
AZ_FWD:=R2D(AZ_FWD);
AZ_FWD:=IFTE(AZ_FWD<0,((AZ_FWD+360) MOD 360),AZ_FWD);//NO NEGATIVE BEARINGS
AZ_FNL:=R2D(AZ_FNL);
AZ_FNL:=IFTE(AZ_FNL<0,((AZ_FNL+360) MOD 360),AZ_FNL);//NO NEGATIVE BEARINGS
AZ_REV:=R2D(AZ_REV);
AZ_REV:=IFTE(AZ_REV<0,((AZ_REV+360) MOD 360),AZ_REV);
// VALIDITY: CONVERGENCE FAILURE
// DISTS: SAME UNITS AS RADIUS
// AZIMUTHS: POSITIVE 0..360
// LOOPCOUNT: CONVERGENCE INFO (LOWER BETTER):DIAGNOSTIC ONLY
// EXPECT SEQUENCE ORDER TO CHANGE
// ESPECIALLY AZIMUTHS
RETURN {(LOOPCOUNT≠LoopSafe),DISTS,(AZ_FWD),(AZ_FNL),(AZ_REV),LOOPCOUNT};
END;
EXPORT G_LL_VincentyDo(RA,RB,FLATTEN,Lat1,Lon1,GdDist,Bearing)
//DIRECT: GIVEN LL DIST BEAR, RETURN LL2
BEGIN
LOCAL USQ,NUM,DEN,TMP;
LOCAL LAT1:=HMS→(Lat1);
LOCAL LON1:=HMS→(Lon1);
LOCAL FNL,REV;
LOCAL SIGMA1;
LOCAL SIGMA,SIN_SIGMA,COS_SIGMA;
LOCAL TWOSIGMAM,TWOSIGMA1,COS_2SIGMAM,COSSQ_2SIGMAM;
LOCAL DELTA_SIGMA,SIGMAOLD;
LOCAL SIN_ALPHA;
LOCAL AAA,BBB,CCC,LLL,ALPHA2;
LOCAL LAMBDA,PHI2,LON2;
LOCAL S_BA;//CALC BEFORE ITERATING
LOCAL PHI1:= D2R(LAT1);//GEODETIC LAT
LOCAL LAMBDA1:=D2R(LON1);
LOCAL ALPHA1:=D2R(Bearing);//INITIAL BEARING
LOCAL SIN_ALPHA1:=SIN(ALPHA1);
LOCAL COS_ALPHA1:=COS(ALPHA1);
LOCAL SINSQ_ALPHA,COSSQ_ALPHA;
MSGBOX("DIRECT NOT WORKING"+NL+"ERRORS>100km"+NL+"Can you spot the bug?");
//CALCULATE UU1
LOCAL TANUU1:=(1-FLATTEN)*TAN(PHI1);//CORRECT
LOCAL UU1:=ATAN(TANUU1);//CORRECT
//LOCAL COSUU1:=COS(UU1);
LOCAL COSUU1:=1/√(1+TANUU1^2);//TRY
//LOCAL SINUU1:=SIN(UU1);// = COS*TAN SAVES A TRIG
LOCAL SINUU1:=COSUU1*TANUU1;//TRY
IF GdDist==0 AND 1 THEN
RETURN {LAT1,LON1,0};//AVOID OR REVEAL ROUNDING ERROR
END;
// V1
SIGMA1:=ATAN2YX(TANUU1,COS_ALPHA1);
// V2
SIN_ALPHA:=COSUU1*SIN_ALPHA1;
SINSQ_ALPHA:=(SIN_ALPHA)^2;
//V3 V4
COSSQ_ALPHA:=1-SINSQ_ALPHA;
USQ:=COSSQ_ALPHA*((RA^2 - RB^2)/(RB^2));
//(V3)
AAA:=1 + (USQ/16384)*(4096 + USQ*(−768 + USQ*(320 - 175*USQ)));
//(V4)
BBB:= (USQ/ 1024)*( 256 + USQ*(−128 + USQ*( 74 - 47*USQ)));
//ITERATE
S_BA:=(GdDist/(RB*AAA));
SIGMA:=S_BA;//1ST EST
REPEAT
// V5
TWOSIGMAM:=2*SIGMA1 + SIGMA;
//V6 PREP
COS_2SIGMAM:=COS(TWOSIGMAM);
SIN_SIGMA:=SIN(SIGMA);
COS_SIGMA:=COS(SIGMA);
//(V6)
COSSQ_2SIGMAM:=(COS_2SIGMAM)^2;
DELTA_SIGMA:=(1/6)*BBB*COS_2SIGMAM*(−3+4*(SIN_SIGMA)^2)*(−3+4*COSSQ_2SIGMAM) ;// V0.4
DELTA_SIGMA:=(1/4)*BBB*(COS_SIGMA*(−1+2*COSSQ_2SIGMAM)) - DELTA_SIGMA ;
DELTA_SIGMA:=BBB*SIN_SIGMA*(COS_2SIGMAM + DELTA_SIGMA);
//PRINT("Δσ "+DELTA_SIGMA);
// V7
SIGMAOLD:=SIGMA;
SIGMA:=S_BA+DELTA_SIGMA;
UNTIL ABS(SIGMA-SIGMAOLD)≤EPSILON_V;
//UPDATE?
//SOME SOURCES OMIT AND USE OLDER VALUES
//SO MAY INHIBIT FOR COMATIBILITY TESTS
IF 1 THEN
SIN_SIGMA:=SIN(SIGMA);
COS_SIGMA:=COS(SIGMA);
END;
//V8 LATITUDE
TMP:=(SINUU1*SIN_SIGMA - COSUU1*COS_SIGMA*COS_ALPHA1);
NUM:=(SINUU1*COS_SIGMA) + (COSUU1*SIN_SIGMA*COS_ALPHA1);
DEN:=((SIN_ALPHA)^2+(TMP)^2);
DEN:=(1-FLATTEN)*√(DEN);
PHI2:=ATAN2YX(NUM,DEN);
//V9 LAMBDA FOR LON
NUM:=SIN_SIGMA*SIN_ALPHA1;
DEN:= (COSUU1*COS_SIGMA) - (SINUU1*SIN_SIGMA*COS_ALPHA1);
LAMBDA:=ATAN2YX(NUM,DEN);
//V10 CCC FOR LON
CCC:=(FLATTEN/16)*COSSQ_ALPHA*(4+FLATTEN*(4-3*COSSQ_ALPHA));
//V11 LLL:=DIFF IN LON
LLL:=(SIGMA + CCC*SIN_SIGMA*(COS_2SIGMAM + CCC*COS_SIGMA*(−1 + 2*COSSQ_2SIGMAM)));
LLL:=LAMBDA - (1 - CCC)*FLATTEN*SIN_ALPHA*LLL;
//V12 ALPHA2:=GEODETIC LATITUDE
//DEN:=(−SINUU1*SIN_SIGMA + COSUU1*COS_SIGMA*COS_ALPHA1);
FNL:=ATAN2YX(SIN_ALPHA,−TMP);
REV:=ATAN2YX(−SIN_ALPHA,TMP);
// RETURN
// PHI2: GEODETIC LATITUDE −90..90
// LONGITUDE: NORMALISE TO −180..180
// BEARING2: FINAL BEARING: NORMALISE TO 0..360
// NORMALISATIONS TBD
PHI2:=R2D(PHI2);
LON2:=LON1+R2D(LLL);
// ALPHA2 == FNL BEARING
FNL:=R2D(FNL);
//FNL:=IFTE(FNL<0,((FNL+360) MOD 360),FNL);//NO NEGATIVE BEARINGS
REV:=R2D(REV);
//REV:=IFTE(REV<0,((REV+360) MOD 360),REV);//NO NEGATIVE BEARINGS
RETURN {PHI2,LON2,REV,FNL};
END;
EXPORT G_DxVincenty(Elevation,Point1,Point2)
// GEODESY INVERSE PROBLEM ELLIPSOIDAL
// USE: RA AND FLATTEN FROM ThisEll TO DERIVE RB
// RETURNS DIST IN SAME UNITS AS R
// LAT AND LON IN DEGREES
BEGIN
LOCAL LAT1:=Point1(G_XLAT);
LOCAL LON1:=Point1(G_XLON);
LOCAL LAT2:=Point2(G_XLAT);
LOCAL LON2:=Point2(G_XLON);
LOCAL RA:=G_ThisEll(G_XRA)+Elevation;
LOCAL FLATTEN:=G_ThisEll(G_XF);
LOCAL RB:=G_RadiusB(RA,FLATTEN);
RETURN G_DxVincentyDo(RA,RB,FLATTEN,LAT1,LON1,LAT2,LON2);
END;
EXPORT G_LL_Vincenty(Elevation,Point1,GdDist,Bearing)
// GEODESY DIRECT PROBLEM ELLIPSOIDAL
// USE: RA AND FLATTEN FROM ThisEll TO DERIVE RB
// Point1: Lat Lon
// GdDistance: Geodesic Distance
// Bearing:
// LAT AND LON IN DEGREES
BEGIN
LOCAL LAT1:=Point1(G_XLAT);
LOCAL LON1:=Point1(G_XLON);
LOCAL RA:=G_ThisEll(G_XRA)+Elevation;
LOCAL FLATTEN:=G_ThisEll(G_XF);
LOCAL RB:=G_RadiusB(RA,FLATTEN);
RETURN G_LL_VincentyDo(RA,RB,FLATTEN,LAT1,LON1,GdDist,Bearing);
END;
G_UseThis(POSN)
//VINCENTY RECALCULATES THESE==DUPLICATION
BEGIN
IF POSN>0 THEN
G_ThisEll:=G_EarthABF(POSN);
RA:=G_ThisEll(G_XRA);
FLATTEN:=G_ThisEll(G_XF);
RB:=G_RadiusB(RA,FLATTEN);
//ELSE.UNCHANGED
END;
RETURN POSN;
END;
EXPORT G_UseN(POSN)
//N==0: USER CHOOSE
//N>0: CHOOSE ENTRY N
BEGIN
LOCAL OKC;
LOCAL NN:=POSN;
IF NN==0 OR NN>SIZE(G_EllABF) THEN
OKC:=CHOOSE(NN,"Ellipsoid",G_EarthABF);
IF OKC THEN
MSGBOX("DETAILS TBD");
RETURN G_UseN(NN);
END;
END;
RETURN G_UseThis(NN);
END;
EXPORT G_UseNew(LST)
//CREATE NEW ELLIPSOID (NO CHECK FOR DUPLICATES)
//LST FORMAT AS IN EXAMPLES
//IDEALLY RETURN 0 IF LIST TOO BIG
BEGIN
G_EllABF(0):=LST;
RETURN G_UseThis(SIZE(G_EllABF));
END;
EXPORT G_UseST(ST)
//SEARCH FOR STRING
//IF FOUND USE THAT ELLIPSOID
//A RADIUS EQU
//B RADIUS POL
//AT PRESENT RB IS ALWAYS DERIVED FROM F
//LATER ALLOW F TO BE DERIVED
BEGIN
LOCAL II,POS;
FOR II FROM 1 TO SIZE(G_EllABF) DO
POS:=INSTRING(G_EllABF(G_XGNAME),ST);
IF POS THEN
G_UseThis(POS);
RETURN POS;
END;
END;
RETURN POS;//0=>UNCHANGED
END;
G_Use(PP)
//NOT YET WORKING-(USING MENU RUN ROUTINE INTERACTIVELY)
//PROGRAM LOGIC LOOKS OK BUT TYPE HANDLING RUNTIME ERR
//INTENDED TO PROVIDE A SINGLE EXPORT TO REPLACE 3 EXPORT ROUTINES
//WHICH COULD THEN BE MADE LOCAL
//(NUM): SELECT ELLIPSOID N FROM LIST (0=SHOW CHOOSE LIST O USER)
// EG 1
//(STR): SELECT 1st ELLIPSOID CONTAINING STRING
// EG "WGS84" or "WGS84_m"
//(LST): CREATE (AND USE) NEW ELLIPSOID
// EG {"My Planet",1,0,0.1}
BEGIN
//MSGBOX(TYPE(PP));
CASE
IF TYPE(PP)==0 THEN //NUM
RETURN G_UseN(PP);
END;
IF TYPE(PP)==2 THEN //STRING:SEARCH FOR THAT ELLIPSOID NAME
RETURN G_UseST(PP);
END;
IF TYPE(PP)==6 THEN //LIST: CREATE AND USE NEW ELLIPSOID
RETURN G_UseNew(PP);
END;
DEFAULT
RETURN MSGBOX("G_Use ERR: NUM OR STR OR LST");
END;
END;
// Main Spherical Geodesy routines-
EXPORT G_SpMeanRadius (RA,RB)
//In Geodesy, this is usually used as a radius for spherical
BEGIN
RETURN (2*RA+RB)/3;
END;
EXPORT G_SpDistHav(Point1,Point2)
BEGIN
LOCAL LAT1:=HMS→(Point1(G_XLAT));
LOCAL LAT2:=HMS→(Point2(G_XLAT));
LOCAL LON1:=HMS→(Point1(G_XLON));
LOCAL LON2:=HMS→(Point2(G_XLON));
LOCAL DLON:=D2R(LON2-LON1);
LOCAL DLAT:=D2R(LAT2-LAT1);
LOCAL AA:=HAV(DLAT) + (COS(D2R(LAT1))*COS(D2R(LAT2))*HAV(DLON));
LOCAL CC:=2*ASIN(MIN(1,√(AA)));//MIN 1 AVOIDS ROUNDING OUTRANGE
LOCAL DD:=RR*CC;
RETURN ROUND(DD,−6);
END;
// End Main Spherical Routines
// Useful routines
EXPORT G_SpSurfaceArea(Radius)
// RETURN SPHERICAL SURFACE AREA IN SAME UNITS AS RR
BEGIN
RETURN 4*π*(Radius^2);
END;
//
EXPORT G_FaceArea(RadiusA,RadiusB)
BEGIN
RETURN π*RadiusA*RadiusB;
END;
G_EllStats(PR)
//PR:PRINT ELSE TEXTOUT
BEGIN
LOCAL II;
LOCAL ST:={};
ST(1):=G_ThisEll(G_XGNAME);
ST(2):="RA: "+RA;
ST(3):="RB: "+RB;
ST(4):="f: "+(FLATTEN);
IF FLATTEN≠0 THEN
ST(5):="1/f: "+(1/FLATTEN);
END;
IF PR THEN
FOR II FROM 1 TO SIZE(ST) DO
PRINT(ST(II));
END;
ELSE
FOR II FROM 1 TO SIZE(ST) DO
TEXTOUT_P(ST(II),0,20*II+100);
END;
END;
END;
G_EX(LAT1,LON1,LAT2,LON2)
//An example displaying results
BEGIN
LOCAL TM,RESULT1,RESULT2,DISTS;
PRINT(" ");
PRINT("From LAT "+HMS→(LAT1)+" "+→HMS(LAT1));
PRINT("From LON "+HMS→(LON1)+" "+→HMS(LON1));
PRINT("To LAT "+HMS→(LAT2)+" "+→HMS(LAT2));
PRINT("To LON "+HMS→(LON2)+" "+→HMS(LON2));
TM:=TICKS;
RESULT1:=G_DistHav({LAT1,LON1},{LAT2,LON2});
RESULT2:=G_DxVincenty(0,{LAT1,LON1},{LAT2,LON2});
TM:=TICKS-TM;
PRINT("Elapsed "+(TM/1000)+" s");
PRINT("Spherical Haversine: "+(RESULT1)+" m");
PRINT("Ellipsoidal Vincenty: "+(RESULT2(2))+" m");
PRINT(" AZ "+(→HMS(RESULT2(3))));
PRINT(" AZ "+(→HMS(RESULT2(4))));
PRINT(" AZ "+(→HMS(RESULT2(5))));
PRINT("Validity: " +(RESULT2(1))+" "+(RESULT2(6)));
//PRINT(RESULT2);
//RESULT:=Z_GEODESY(LAT1,LON1,LAT2,LON2);
//WAIT;
END;
EXPORT G_ExampleUse()
//A short example of key procedures
BEGIN
LOCAL OKC,DIST,RESULT,EXAMPLEPOINT,XLAT,XLON;
OKC:=G_UseNew({"Add your own",{RA,0,FLATTEN}});//Insert into List. Use use to use it.
OKC:=G_UseST("WGS84_m");//DEFAULT
IF OKC THEN //Found in list
DIST:=G_DistHav({1,2},{3,4});
//RESULT:=G_DxVincentyDo(RA,RB,FLATTEN,LAT1,LON1,LAT2,LON2);
RESULT:=G_DxVincenty(0,{1,2},{3,4});
//Variables prefixed X index into data structures
EXAMPLEPOINT:={"From Here",1,2};
XLAT:=2; XLON:=3;
XLAT:=1; XLON:=2; //REINSTATE CURRENT DEFAULT
EXAMPLEPOINT:={1,2};
//and similarly for geoid structure
END;
END;
EXPORT G_DrawEll(Flatten)
//The stats are omitted because flattening drawn may differ
//makes illustrate changing f easier
//for now use myRadiusA of 100 pixels
BEGIN
LOCAL XX,YY;
LOCAL DX,DY;
LOCAL MyRadiusA:=100;
LOCAL HEIGHT:=100;
LOCAL WIDTH:=100;
//LOCAL SCALE:=RadiusA/WIDTH;
LOCAL EWIDTH:=MyRadiusA;
LOCAL EHEIGHT:=((1-Flatten)*MyRadiusA);//RADIUS POLAR
RECT();
//G_EllStats(0);
FOR YY FROM −HEIGHT TO HEIGHT DO
//PRINT(YY);
FOR XX FROM −WIDTH TO WIDTH DO
DX:=XX/(EWIDTH);
DY:=YY/(EHEIGHT);
IF (DX^2+DY^2)≤1 THEN
//PRINT({WIDTH+10+XX,HEIGHT+20+YY});
PIXON_P(WIDTH+100+XX,HEIGHT+20+YY,#009900h);
END;//IF
END;//FOR
END;//FOR
//G_EllStats(0);
WAIT;
END;
G_TEST()
//A few tests. More needed.
BEGIN
G_EX(50°03′58.76″,−005°42′53.10″,58°38′38.48″,−003°04′12.34″);//50PERCENT
PRINT("EXPECT 969 954_m");
G_EX(−37°57′3.72030″,144°25′29.52440″,−37°39′10.15610″,143°55′35.38390″);
PRINT("EXPECTED: 54972.271m");//VINCENTY AUSTRALIA
PRINT("306d52 05 * 127d10 25 * 307d10 25.07 ");
G_EX(31.8300167,35.0662833,31.83000000,35.0708167);
PRINT("PYTHON: CIRCLE 428.4 VIN 429.1677 m");
//G_EX(0,0,0.5,179.5);
G_EX(0,0,0.5,179.7);//WILL NOT CONVERGE
//G_EX(0,0,0,90);
//G_EX(0,0,0.5,0.5);
G_EX(2,2,2,2);
G_EX(50°03′58.76″,−005°42′53.10″,58°38′38.48″,−003°04′12.34″);//50PERCENT
PRINT("EXPECT 969 954_m");
//IF 0 THEN
G_EX(50°3′59″,−5°42′53″,58°38′38″,−3°4′12″);
PRINT("RECHECK: EXPECTED 968.9km 009DEG0711");
IF 0 THEN
PRINT(G_LL_Vincenty(0,{45,30},0,0));
PRINT(→HMS(G_LL_Vincenty(0,{−37°57′3.72030″,144°25′29.52440″},54972.271,306°52′5.37″)));//AU
PRINT("EXPECTED:"+→HMS({−37°39′10.5610″,143°55′35.38390″,127°10′25.07″,307°10′2507″}));
END;
IF 0 THEN
G_EX(−37°57′3.72030″,144°25′29.52440″,−37°39′10.15610″,143°55′35.38390″);//EXPECTS 54972m
PRINT("EXPECTED: 54972.271m ");
PRINT("306d52:5.37 / 127d10:25.07");//VINCENTY AUSTRALIA
//IF 1 THEN
G_EX(31.8300167,35.0662833,31.83000000,35.0708167);
PRINT("PYTHON: CIRCLE 428.4 VIN 429.1677 m");
//IF 0 THEN
G_EX(0,0,0.5,179.5);
G_EX(0,0,0.5,179.7);//WILL NOT CONVERGE
//IF 1 THEN
G_EX(0,0,0,90);
G_EX(0,0,0.5,0.5);
G_EX(2,2,2,2);
END;
END;
EXPORT GEODESY()
BEGIN
G_About();
G_EllStats(1);
PRINT("Esc");WAIT;
IF MSGBOX("Run Test?") THEN
G_TEST();
END;
END;A compiler option to ignore backups would be useful.
RE: Geodesy API - StephenG1CMZ - 11-23-2016 03:03 PM
Version 0.32
Code:
// GEODESY © 2016 StephenG1CMZ
LOCAL VER:="0.32";
LOCAL NL:=CHAR(10);
LOCAL SP:=" ";
//FORWARD
G_Flattening(); //Derive Flattening from radius RA RB
G_SpMeanRadius();//Derive mean radius from radius RA RB
G_RadiusB(); //Derive RB from RA and flattening
G_Use(); //Use this ellipsoid
//END FORWARD
//We assume RA≥RB also Trig in Radians
// Ell=ElIpsoid
// Sp=Spherical
//Geoid Data
EXPORT WGS84a_m:=6378137.0;
EXPORT WGS84a_km:=WGS84a_m/1000;
LOCAL WGS84fINV:=298.257223563;
EXPORT WGS84f:=1/WGS84fINV;
LOCAL WGS84b_m:=G_RadiusB(WGS84a_m,WGS84f);
LOCAL WGS84b_km:=WGS84b_m/1000;
LOCAL RA:=WGS84a_m; //RADIUS EQ
LOCAL FLATTEN:=WGS84f;
LOCAL RB:=G_RadiusB(RA,FLATTEN);
LOCAL RR:=G_SpMeanRadius(RA,RB);
LOCAL EPSILON_V:=1ᴇ−12;
LOCAL LoopSafe:=140;
//FORMAT:
//RA RB F
EXPORT G_XGNAME:=1;
EXPORT G_XRA:=2;
EXPORT G_XF:=4;
EXPORT G_EllABF:={
{"WGS84_m", WGS84a_m, 0, WGS84f},
{"Bessel_m", 6377397.155, 0, 1/299.1528128},
{"International_m", 6378388.000, 0, 1/297},
{"Sphere", 1,0,0},
{"Mercury_km", 2439.7, 2439.7, 0.0000},
{"Venus_km", 6051.8, 6051.8, 0.0000},
{"Earth WGS84_km", WGS84a_km, 0, WGS84f},
{"Moon_km", 1738.1, 1736.0, 0.0012},
{"Mars_km", 3396.2, 3376.2, 0.00589},
{"Jupiter_km", 71492, 66854, 0.06487},
{"Saturn_km", 60268, 54364, 0.09796},
{"Uranus_km", 22559, 24973, 0.02293},
{"Neptune_km", 24764, 24341, 0.01708},
{"Pluto_km", 1187, 1187, 0.0000},
{"Sun", 695508, 0, 0.00005},
{"Kepler 11145123 _km", 1.5ᴇ6, 1.5ᴇ6-3, G_Flattening(1.5ᴇ6,1.5ᴇ6-3)}};
// Kepler 11145123 is roundest object in nature as at 2016-11
EXPORT G_ThisEll:=G_EllABF(1);
EXPORT G_XLAT:=1;
EXPORT G_XLON:=2;
LOCAL ESS_UNITS:="m";
EXPORT G_About()
BEGIN
PRINT();
PRINT("Geodesy API V "+VER+" © 2016 StephenG1CMZ");
PRINT("NOT FOR NAVIGATIONAL USE");
PRINT("No liability accepted");
PRINT("Expect the ordering of Vincenty results to change in later versions");
END;
// STD MATHS
D2R(DEGS)
BEGIN
RETURN DEGS*(π/180);
END;
R2D(RADS)
BEGIN
RETURN RADS*(180/π);
END;
ATAN2YX(YY,XX)
BEGIN //MY ATAN2
IF XX==0 AND YY==0 THEN
RETURN 0; //CONVENIENT;
END;
//-π..π THE SPECIAL CHAR BELOW IS i=sqrt(−1), NOT CHINESE
RETURN ARG(XX+*YY);//BUILT-IN HANDLES MOST. END;
END;
HAV(THETA)
//TRIG: HAVERSINE=SIN^2(THETA/2)
BEGIN
LOCAL SN:=SIN(THETA/2);
RETURN SN^2;
END;
// END STD MATHS
EXPORT G_Flattening(RadiusA,RadiusB)
BEGIN
RETURN (RadiusA-RadiusB)/RadiusA;
END;
EXPORT G_RadiusB(RadiusA,Flatten)
//RETURN POLAR RADIUS
BEGIN
RETURN (1-Flatten)*RadiusA;//RADIUS POLAR
END;
EXPORT G_DxVincentyDo(RA,RB,FLATTEN,Lat1,Lon1,Lat2,Lon2)
// GEODESY INVERSE PROBLEM ELLIPSOIDAL
// USE: RA AND FLATTEN FROM ThisEll
// RETURNS DIST IN SAME UNITS AS R
// LAT AND LON IN DEGREES
BEGIN
LOCAL NUM,DEN;
LOCAL LAT1:=HMS→(Lat1);
LOCAL LAT2:=HMS→(Lat2);
LOCAL LON1:=HMS→(Lon1);
LOCAL LON2:=HMS→(Lon2);
LOCAL UU1,SINUU1,COSUU1; //U1 AUX LATITUDE
LOCAL UU2,SINUU2,COSUU2; //U2 AUX LATITUDE
LOCAL USQ;
LOCAL AAA,BBB;
LOCAL LONDIFFradians; //L=DIFF IN LON
LOCAL PHI1,PHI2;
LOCAL LAMBDA,SINLAMBDA,COSLAMBDA,LAMBDAOLD;
LOCAL SIGMA,SIN_SIGMA,COS_SIGMA,DELTA_SIGMA,SINSQ_SIGMA; //AUX ARC LENGTH
LOCAL SIN_ALPHA,COSSQ_ALPHA;
LOCAL COS_2SIGMAM,COSSQ_2SIGMAM;
LOCAL CC;
LOCAL BRACKET_TERM;
LOCAL LOOPCOUNT;
LOCAL DISTS;
LOCAL AZ_FWD,AZ_FNL,AZ_REV;
// CALC L
LONDIFFradians:=D2R(LON2-LON1); //DIFF IN LON
// CALC U1 U2
PHI1:=D2R(LAT1);
PHI2:=D2R(LAT2);
UU1:=ATAN((1-FLATTEN)*TAN(PHI1));//CORRECT
UU2:=ATAN((1-FLATTEN)*TAN(PHI2));//CORRECT
//TO CALC LAMBDA:
// PRE-CALC U1 TRIG ONCEONLY
SINUU1:=SIN(UU1); COSUU1:=COS(UU1);
// PRE-CALC U2 TRIG EACH POSITION
SINUU2:=SIN(UU2); COSUU2:=COS(UU2);
//ITERATE CALCULATING LAMBDA
LOOPCOUNT:=0;
LAMBDA:=LONDIFFradians; //V13: DIFF IN LON ON AUX. FIRST APPROX
REPEAT
LOOPCOUNT:=LOOPCOUNT+1;
SINLAMBDA:=SIN(LAMBDA); COSLAMBDA:=COS(LAMBDA); //PRECALC
//V14
BRACKET_TERM:=(COSUU1*SINUU2 - SINUU1*COSUU2*COSLAMBDA);
SINSQ_SIGMA:=(COSUU2*SINLAMBDA)^2 + (BRACKET_TERM)^2;
IF SINSQ_SIGMA==0 THEN
//PRINT("SINσ 0: NO DISTANCE");
//0 DISTANCE MIGHT SUGGEST UNDERFLOW
RETURN {1,0,0,0,0,0};//VALID ZERO
END;
SIN_SIGMA:=√(SINSQ_SIGMA);
//V15
COS_SIGMA:=SINUU1*SINUU2 + COSUU1*COSUU2*COSLAMBDA;
//V16
SIGMA:=ATAN2YX(SIN_SIGMA,COS_SIGMA);
//V17
SIN_ALPHA:=(COSUU1*COSUU2*SINLAMBDA)/SIN_SIGMA;
COSSQ_ALPHA:=1-(SIN_ALPHA*SIN_ALPHA);
//V18
IF COSSQ_ALPHA==0 THEN
//PRINT("AVOID COSSQα 0");
COS_2SIGMAM:=0;
ELSE
COS_2SIGMAM:=COS_SIGMA - (2*SINUU1*SINUU2/COSSQ_ALPHA);
END;
CC:=(FLATTEN/16)*COSSQ_ALPHA*(4 + FLATTEN*(4 - 3*COSSQ_ALPHA));
LAMBDAOLD:=LAMBDA;
LAMBDA:=LONDIFFradians + (1-CC)*FLATTEN*SIN_ALPHA*(SIGMA + CC*SIN_SIGMA*(COS_2SIGMAM + CC*COS_SIGMA*(−1+2*(COS_2SIGMAM*COS_2SIGMAM))));
//PRINT(LAMBDA);
UNTIL (ABS(LAMBDA-LAMBDAOLD)≤EPSILON_V) OR (LOOPCOUNT≥LoopSafe);
//V19
USQ:=COSSQ_ALPHA*((RA^2 - RB^2)/(RB^2));
//(V3)
AAA:=1 + (USQ/16384)*(4096 + USQ*(−768 + USQ*(320 - 175*USQ)));
//(V4)
BBB:= (USQ/ 1024)*( 256 + USQ*(−128 + USQ*( 74 - 47*USQ)));
//(V6)
COSSQ_2SIGMAM:=(COS_2SIGMAM)^2;
//DELTA_SIGMA:=(1/6)*BBB*COS_2SIGMAM*(−3+4*(SIN_ALPHA)^2)*(−3+4*COSSQ_2SIGMAM) ;//V0.3
DELTA_SIGMA:=(1/6)*BBB*COS_2SIGMAM*(−3+4*(SIN_SIGMA)^2)*(−3+4*COSSQ_2SIGMAM) ; //V0.4
DELTA_SIGMA:=(1/4)*BBB*(COS_SIGMA*(−1+2*COSSQ_2SIGMAM)) - DELTA_SIGMA ;
DELTA_SIGMA:=BBB*SIN_SIGMA*(COS_2SIGMAM + DELTA_SIGMA);
//PRINT("Δσ "+DELTA_SIGMA);
DISTS:=AAA*(SIGMA - DELTA_SIGMA);//ANGULAR ≤π Radians
DISTS:=RB*DISTS; //same units as radius
//V20:FORWARD AZIMUTH
//CHECKED SINLAMBDA:ORIGINAL (LON2-LON1) OR LOOPED VALUE? USE LOOPED
AZ_FWD:=ATAN2YX((COSUU2*SINLAMBDA),BRACKET_TERM);
//V21:FINAL AND REVERSE AZIMUTH
NUM:=(COSUU1*SINLAMBDA);
DEN:=(−SINUU1*COSUU2 + COSUU1*SINUU2*COSLAMBDA);
AZ_FNL:=ATAN2YX(NUM,DEN);
AZ_REV:=ATAN2YX(−NUM,−DEN);
//RETURN RESULTS
DISTS:=ROUND(DISTS,3);//VINCENTY CLAIMS mm,
AZ_FWD:=R2D(AZ_FWD);
AZ_FWD:=IFTE(AZ_FWD<0,((AZ_FWD+360) MOD 360),AZ_FWD);//NO NEGATIVE BEARINGS
AZ_FNL:=R2D(AZ_FNL);
AZ_FNL:=IFTE(AZ_FNL<0,((AZ_FNL+360) MOD 360),AZ_FNL);//NO NEGATIVE BEARINGS
AZ_REV:=R2D(AZ_REV);
AZ_REV:=IFTE(AZ_REV<0,((AZ_REV+360) MOD 360),AZ_REV);
// VALIDITY: CONVERGENCE FAILURE
// DISTS: SAME UNITS AS RADIUS
// AZIMUTHS: POSITIVE 0..360
// LOOPCOUNT: CONVERGENCE INFO (LOWER BETTER):DIAGNOSTIC ONLY
// EXPECT SEQUENCE ORDER TO CHANGE
// ESPECIALLY AZIMUTHS
RETURN {(LOOPCOUNT≠LoopSafe),DISTS,(AZ_FWD),(AZ_FNL),(AZ_REV),LOOPCOUNT};
END;
EXPORT G_LL_VincentyDo(RA,RB,FLATTEN,Lat1,Lon1,GdDist,Bearing)
//DIRECT: GIVEN LL DIST BEAR, RETURN LL2
BEGIN
LOCAL USQ,NUM,DEN,TMP;
LOCAL LAT1:=HMS→(Lat1);
LOCAL LON1:=HMS→(Lon1);
LOCAL FNL,REV;
LOCAL SIGMA1;
LOCAL SIGMA,SIN_SIGMA,COS_SIGMA;
LOCAL TWOSIGMAM,TWOSIGMA1,COS_2SIGMAM,COSSQ_2SIGMAM;
LOCAL DELTA_SIGMA,SIGMAOLD;
LOCAL SIN_ALPHA;
LOCAL AAA,BBB,CCC,LLL,ALPHA2;
LOCAL LAMBDA,PHI2,LON2;
LOCAL S_BA;//CALC BEFORE ITERATING
LOCAL PHI1:= D2R(LAT1);//GEODETIC LAT
LOCAL LAMBDA1:=D2R(LON1);
LOCAL ALPHA1:=D2R(Bearing);//INITIAL BEARING
LOCAL SIN_ALPHA1:=SIN(ALPHA1);
LOCAL COS_ALPHA1:=COS(ALPHA1);
LOCAL SINSQ_ALPHA,COSSQ_ALPHA;
MSGBOX("DIRECT NOT WORKING"+NL+"ERRORS>100km"+NL+"Can you spot the bug?");
//CALCULATE UU1
LOCAL TANUU1:=(1-FLATTEN)*TAN(PHI1);//CORRECT
LOCAL UU1:=ATAN(TANUU1);//CORRECT
//LOCAL COSUU1:=COS(UU1);
LOCAL COSUU1:=1/√(1+TANUU1^2);//TRY
//LOCAL SINUU1:=SIN(UU1);// = COS*TAN SAVES A TRIG
LOCAL SINUU1:=COSUU1*TANUU1;//TRY
IF GdDist==0 AND 1 THEN
RETURN {LAT1,LON1,0};//AVOID OR REVEAL ROUNDING ERROR
END;
// V1
SIGMA1:=ATAN2YX(TANUU1,COS_ALPHA1);
// V2
SIN_ALPHA:=COSUU1*SIN_ALPHA1;
SINSQ_ALPHA:=(SIN_ALPHA)^2;
//V3 V4
COSSQ_ALPHA:=1-SINSQ_ALPHA;
USQ:=COSSQ_ALPHA*((RA^2 - RB^2)/(RB^2));
//(V3)
AAA:=1 + (USQ/16384)*(4096 + USQ*(−768 + USQ*(320 - 175*USQ)));
//(V4)
BBB:= (USQ/ 1024)*( 256 + USQ*(−128 + USQ*( 74 - 47*USQ)));
//ITERATE
S_BA:=(GdDist/(RB*AAA));
SIGMA:=S_BA;//1ST EST
REPEAT
// V5
TWOSIGMAM:=2*SIGMA1 + SIGMA;
//V6 PREP
COS_2SIGMAM:=COS(TWOSIGMAM);
SIN_SIGMA:=SIN(SIGMA);
COS_SIGMA:=COS(SIGMA);
//(V6)
COSSQ_2SIGMAM:=(COS_2SIGMAM)^2;
DELTA_SIGMA:=(1/6)*BBB*COS_2SIGMAM*(−3+4*(SIN_SIGMA)^2)*(−3+4*COSSQ_2SIGMAM) ;// V0.4
DELTA_SIGMA:=(1/4)*BBB*(COS_SIGMA*(−1+2*COSSQ_2SIGMAM)) - DELTA_SIGMA ;
DELTA_SIGMA:=BBB*SIN_SIGMA*(COS_2SIGMAM + DELTA_SIGMA);
//PRINT("Δσ "+DELTA_SIGMA);
// V7
SIGMAOLD:=SIGMA;
SIGMA:=S_BA+DELTA_SIGMA;
UNTIL ABS(SIGMA-SIGMAOLD)≤EPSILON_V;
//UPDATE?
//SOME SOURCES OMIT AND USE OLDER VALUES
//SO MAY INHIBIT FOR COMATIBILITY TESTS
IF 1 THEN
SIN_SIGMA:=SIN(SIGMA);
COS_SIGMA:=COS(SIGMA);
END;
//V8 LATITUDE
TMP:=(SINUU1*SIN_SIGMA - COSUU1*COS_SIGMA*COS_ALPHA1);
NUM:=(SINUU1*COS_SIGMA) + (COSUU1*SIN_SIGMA*COS_ALPHA1);
DEN:=((SIN_ALPHA)^2+(TMP)^2);
DEN:=(1-FLATTEN)*√(DEN);
PHI2:=ATAN2YX(NUM,DEN);
//V9 LAMBDA FOR LON
NUM:=SIN_SIGMA*SIN_ALPHA1;
DEN:= (COSUU1*COS_SIGMA) - (SINUU1*SIN_SIGMA*COS_ALPHA1);
LAMBDA:=ATAN2YX(NUM,DEN);
//V10 CCC FOR LON
CCC:=(FLATTEN/16)*COSSQ_ALPHA*(4+FLATTEN*(4-3*COSSQ_ALPHA));
//V11 LLL:=DIFF IN LON
LLL:=(SIGMA + CCC*SIN_SIGMA*(COS_2SIGMAM + CCC*COS_SIGMA*(−1 + 2*COSSQ_2SIGMAM)));
LLL:=LAMBDA - (1 - CCC)*FLATTEN*SIN_ALPHA*LLL;
//V12 ALPHA2:=GEODETIC LATITUDE
//DEN:=(−SINUU1*SIN_SIGMA + COSUU1*COS_SIGMA*COS_ALPHA1);
FNL:=ATAN2YX(SIN_ALPHA,−TMP);
REV:=ATAN2YX(−SIN_ALPHA,TMP);
// RETURN
// PHI2: GEODETIC LATITUDE −90..90
// LONGITUDE: NORMALISE TO −180..180
// BEARING2: FINAL BEARING: NORMALISE TO 0..360
// NORMALISATIONS TBD
PHI2:=R2D(PHI2);
LON2:=LON1+R2D(LLL);
// ALPHA2 == FNL BEARING
FNL:=R2D(FNL);
FNL:=IFTE(FNL<0,((FNL+360) MOD 360),FNL);//NO NEGATIVE BEARINGS
REV:=R2D(REV);
REV:=IFTE(REV<0,((REV+360) MOD 360),REV);//NO NEGATIVE BEARINGS
RETURN {PHI2,LON2,REV,FNL};
END;
EXPORT G_DxVincenty(Elevation,Point1,Point2)
// GEODESY INVERSE PROBLEM ELLIPSOIDAL
// USE: RA AND FLATTEN FROM ThisEll TO DERIVE RB
// RETURNS DIST IN SAME UNITS AS R
// LAT AND LON IN DEGREES
BEGIN
LOCAL LAT1:=Point1(G_XLAT);
LOCAL LON1:=Point1(G_XLON);
LOCAL LAT2:=Point2(G_XLAT);
LOCAL LON2:=Point2(G_XLON);
LOCAL RA:=G_ThisEll(G_XRA)+Elevation;
LOCAL FLATTEN:=G_ThisEll(G_XF);
LOCAL RB:=G_RadiusB(RA,FLATTEN);
RETURN G_DxVincentyDo(RA,RB,FLATTEN,LAT1,LON1,LAT2,LON2);
END;
EXPORT G_LL_Vincenty(Elevation,Point1,GdDist,Bearing)
// GEODESY DIRECT PROBLEM ELLIPSOIDAL
// USE: RA AND FLATTEN FROM ThisEll TO DERIVE RB
// Point1: Lat Lon
// GdDistance: Geodesic Distance
// Bearing:
// LAT AND LON IN DEGREES
BEGIN
LOCAL LAT1:=Point1(G_XLAT);
LOCAL LON1:=Point1(G_XLON);
LOCAL RA:=G_ThisEll(G_XRA)+Elevation;
LOCAL FLATTEN:=G_ThisEll(G_XF);
LOCAL RB:=G_RadiusB(RA,FLATTEN);
RETURN G_LL_VincentyDo(RA,RB,FLATTEN,LAT1,LON1,GdDist,Bearing);
END;
G_UseThis(POSN)
//VINCENTY RECALCULATES THESE==DUPLICATION
BEGIN
IF POSN>0 THEN
G_ThisEll:=G_EllABF(POSN);
RA:=G_ThisEll(G_XRA);
FLATTEN:=G_ThisEll(G_XF);
RB:=G_RadiusB(RA,FLATTEN);
//ELSE.UNCHANGED
END;
RETURN POSN;
END;
EXPORT G_UseN(POSN)
//N==0: USER CHOOSE
//N>0: CHOOSE ENTRY N
BEGIN
LOCAL OKC;
LOCAL NN:=POSN;
IF NN==0 OR NN>SIZE(G_EllABF) THEN
OKC:=CHOOSE(NN,"Ellipsoid",G_EllABF);
IF OKC THEN
MSGBOX("DETAILS TBD");
RETURN G_UseN(NN);
END;
END;
RETURN G_UseThis(NN);
END;
EXPORT G_UseNew(LST)
//CREATE NEW ELLIPSOID (NO CHECK FOR DUPLICATES)
//LST FORMAT AS IN EXAMPLES
//IDEALLY RETURN 0 IF LIST TOO BIG
BEGIN
G_EllABF(0):=LST;
RETURN G_UseThis(SIZE(G_EllABF));
END;
EXPORT G_UseST(ST)
//SEARCH FOR STRING
//IF FOUND USE THAT ELLIPSOID
//A RADIUS EQU
//B RADIUS POL
//AT PRESENT RB IS ALWAYS DERIVED FROM F
//LATER ALLOW F TO BE DERIVED
BEGIN
LOCAL II,POS;
FOR II FROM 1 TO SIZE(G_EllABF) DO
POS:=INSTRING(G_EllABF(G_XGNAME),ST);
IF POS THEN
G_UseThis(POS);
RETURN POS;
END;
END;
RETURN POS;//0=>UNCHANGED
END;
G_Use(PP)
//NOT YET WORKING-(USING MENU RUN ROUTINE INTERACTIVELY)
//PROGRAM LOGIC LOOKS OK BUT TYPE HANDLING RUNTIME ERR
//INTENDED TO PROVIDE A SINGLE EXPORT TO REPLACE 3 EXPORT ROUTINES
//WHICH COULD THEN BE MADE LOCAL
//(NUM): SELECT ELLIPSOID N FROM LIST (0=SHOW CHOOSE LIST O USER)
// EG 1
//(STR): SELECT 1st ELLIPSOID CONTAINING STRING
// EG "WGS84" or "WGS84_m"
//(LST): CREATE (AND USE) NEW ELLIPSOID
// EG {"My Planet",1,0,0.1}
BEGIN
//MSGBOX(TYPE(PP));
CASE
IF TYPE(PP)==0 THEN //NUM
RETURN G_UseN(PP);
END;
IF TYPE(PP)==2 THEN //STRING:SEARCH FOR THAT ELLIPSOID NAME
RETURN G_UseST(PP);
END;
IF TYPE(PP)==6 THEN //LIST: CREATE AND USE NEW ELLIPSOID
RETURN G_UseNew(PP);
END;
DEFAULT
RETURN MSGBOX("G_Use ERR: NUM OR STR OR LST");
END;
END;
// Main Spherical Geodesy routines-
EXPORT G_SpMeanRadius(RA,RB)
//In Geodesy, this is usually used as a radius for spherical
BEGIN
RETURN (2*RA+RB)/3;
END;
EXPORT G_SpDistHav(Point1,Point2)
BEGIN
LOCAL LAT1:=HMS→(Point1(G_XLAT));
LOCAL LAT2:=HMS→(Point2(G_XLAT));
LOCAL LON1:=HMS→(Point1(G_XLON));
LOCAL LON2:=HMS→(Point2(G_XLON));
LOCAL DLON:=D2R(LON2-LON1);
LOCAL DLAT:=D2R(LAT2-LAT1);
LOCAL AA:=HAV(DLAT) + (COS(D2R(LAT1))*COS(D2R(LAT2))*HAV(DLON));
LOCAL CC:=2*ASIN(MIN(1,√(AA)));//MIN 1 AVOIDS ROUNDING OUTRANGE
LOCAL DD:=RR*CC;
RETURN ROUND(DD,−6);
END;
// End Main Spherical Routines
// Useful routines
EXPORT G_SpSurfaceArea(Radius)
// RETURN SPHERICAL SURFACE AREA IN SAME UNITS AS RR
BEGIN
RETURN 4*π*(Radius^2);
END;
//
EXPORT G_FaceArea(RadiusA,RadiusB)
BEGIN
RETURN π*RadiusA*RadiusB;
END;
G_EllStats(PR)
//PR:PRINT ELSE TEXTOUT
BEGIN
LOCAL II;
LOCAL ST:={};
ST(1):=G_ThisEll(G_XGNAME);
ST(2):="RA: "+RA;
ST(3):="RB: "+RB;
ST(4):="f: "+(FLATTEN);
IF FLATTEN≠0 THEN
ST(5):="1/f: "+(1/FLATTEN);
END;
IF PR THEN
FOR II FROM 1 TO SIZE(ST) DO
PRINT(ST(II));
END;
ELSE
FOR II FROM 1 TO SIZE(ST) DO
TEXTOUT_P(ST(II),0,20*II+100);
END;
END;
END;
G_EX(LAT1,LON1,LAT2,LON2)
//An example displaying results
BEGIN
LOCAL TM,RESULT1,RESULT2,DISTS;
PRINT(" ");
PRINT("From LAT "+HMS→(LAT1)+" "+→HMS(LAT1));
PRINT("From LON "+HMS→(LON1)+" "+→HMS(LON1));
PRINT("To LAT "+HMS→(LAT2)+" "+→HMS(LAT2));
PRINT("To LON "+HMS→(LON2)+" "+→HMS(LON2));
TM:=TICKS;
RESULT1:=G_SpDistHav({LAT1,LON1},{LAT2,LON2});
RESULT2:=G_DxVincenty(0,{LAT1,LON1},{LAT2,LON2});
TM:=TICKS-TM;
PRINT("Elapsed "+(TM/1000)+" s");
PRINT("Spherical Haversine: "+(RESULT1)+" m");
PRINT("Ellipsoidal Vincenty: "+(RESULT2(2))+" m");
PRINT(" AZ "+(→HMS(RESULT2(3))));
PRINT(" AZ "+(→HMS(RESULT2(4))));
PRINT(" AZ "+(→HMS(RESULT2(5))));
PRINT("Validity: " +(RESULT2(1))+" "+(RESULT2(6)));
//PRINT(RESULT2);
//RESULT:=Z_GEODESY(LAT1,LON1,LAT2,LON2);
//WAIT;
END;
EXPORT G_ExampleUse()
//A short example of key procedures
BEGIN
LOCAL OKC,DIST,RESULT,EXAMPLEPOINT,XLAT,XLON;
OKC:=G_UseNew({"Add your own",{RA,0,FLATTEN}});//Insert into List. Use use to use it.
OKC:=G_UseST("WGS84_m");//DEFAULT
IF OKC THEN //Found in list
DIST:=G_SpDistHav({1,2},{3,4});
//RESULT:=G_DxVincentyDo(RA,RB,FLATTEN,LAT1,LON1,LAT2,LON2);
RESULT:=G_DxVincenty(0,{1,2},{3,4});
//Variables prefixed X index into data structures
EXAMPLEPOINT:={"From Here",1,2};
XLAT:=2; XLON:=3;
XLAT:=1; XLON:=2; //REINSTATE CURRENT DEFAULT
EXAMPLEPOINT:={1,2};
//and similarly for geoid structure
END;
END;
EXPORT G_DrawEll(Flatten)
//The stats are omitted because flattening drawn may differ
//makes illustrate changing f easier
//for now use myRadiusA of 100 pixels
BEGIN
LOCAL XX,YY;
LOCAL DX,DY;
LOCAL MyRadiusA:=100;
LOCAL HEIGHT:=100;
LOCAL WIDTH:=100;
//LOCAL SCALE:=RadiusA/WIDTH;
LOCAL EWIDTH:=MyRadiusA;
LOCAL EHEIGHT:=((1-Flatten)*MyRadiusA);//RADIUS POLAR
RECT();
//G_EllStats(0);
FOR YY FROM −HEIGHT TO HEIGHT DO
//PRINT(YY);
FOR XX FROM −WIDTH TO WIDTH DO
DX:=XX/(EWIDTH);
DY:=YY/(EHEIGHT);
IF (DX^2+DY^2)≤1 THEN
//PRINT({WIDTH+10+XX,HEIGHT+20+YY});
PIXON_P(WIDTH+100+XX,HEIGHT+20+YY,#009900h);
END;//IF
END;//FOR
END;//FOR
//G_EllStats(0);
WAIT;
END;
G_TEST()
//A few tests. More needed.
BEGIN
G_EX(50°03′58.76″,−005°42′53.10″,58°38′38.48″,−003°04′12.34″);//50PERCENT
PRINT("EXPECT 969 954_m");
G_EX(−37°57′3.72030″,144°25′29.52440″,−37°39′10.15610″,143°55′35.38390″);
PRINT("EXPECTED: 54972.271m");//VINCENTY AUSTRALIA
PRINT("306d52 05 * 127d10 25 * 307d10 25.07 ");
G_EX(31.8300167,35.0662833,31.83000000,35.0708167);
PRINT("PYTHON: CIRCLE 428.4 VIN 429.1677 m");
//G_EX(0,0,0.5,179.5);
G_EX(0,0,0.5,179.7);//WILL NOT CONVERGE
//G_EX(0,0,0,90);
//G_EX(0,0,0.5,0.5);
G_EX(2,2,2,2);
G_EX(50°03′58.76″,−005°42′53.10″,58°38′38.48″,−003°04′12.34″);//50PERCENT
PRINT("EXPECT 969 954_m");
//IF 0 THEN
G_EX(50°3′59″,−5°42′53″,58°38′38″,−3°4′12″);
PRINT("RECHECK: EXPECTED 968.9km 009DEG0711");
IF 0 THEN //DIRECT TESTS
PRINT(G_LL_Vincenty(0,{45,30},0,0));
PRINT(→HMS(G_LL_Vincenty(0,{−37°57′3.72030″,144°25′29.52440″},54972.271,306°52′5.37″)));//AU
PRINT("EXPECTED:"+→HMS({−37°39′10.5610″,143°55′35.38390″,127°10′25.07″,307°10′2507″}));
END;
IF 0 THEN
G_EX(−37°57′3.72030″,144°25′29.52440″,−37°39′10.15610″,143°55′35.38390″);//EXPECTS 54972m
PRINT("EXPECTED: 54972.271m ");
PRINT("306d52:5.37 / 127d10:25.07");//VINCENTY AUSTRALIA
//IF 1 THEN
G_EX(31.8300167,35.0662833,31.83000000,35.0708167);
PRINT("PYTHON: CIRCLE 428.4 VIN 429.1677 m");
//IF 0 THEN
G_EX(0,0,0.5,179.5);
G_EX(0,0,0.5,179.7);//WILL NOT CONVERGE
//IF 1 THEN
G_EX(0,0,0,90);
G_EX(0,0,0.5,0.5);
G_EX(2,2,2,2);
END;
END;
EXPORT GEODESY()
BEGIN
G_About();
G_EllStats(1);
PRINT("Esc");WAIT;
IF MSGBOX("Run Test?") THEN
G_TEST();
END;
END;RE: Geodesy API - StephenG1CMZ - 12-13-2016 08:57 PM
[version 0.33 is a backup containing investigative test data.
If I am understanding the results correctly, two identical calls of the Direct Vincenty routine return different results . First 41 degrees latitude is returned, then 47 degrees. Direct spherical also returns 47. Not sure when I will get the time to look into this, so thought I ought to mention the mistake in case anyone is using V0.32.
Code:
// GEODESY © 2016 StephenG1CMZ
LOCAL VER:="0.33";
LOCAL NL:=CHAR(10);
LOCAL SP:=" ";
//FORWARD
G_Flattening(); //Derive Flattening from radius RA RB
G_SpMeanRadius();//Derive mean radius from radius RA RB
G_RadiusB(); //Derive RB from RA and flattening
G_Use(); //Use this ellipsoid
//END FORWARD
//We assume RA≥RB also Trig in Radians
// Ell=Ellipsoid
// Dx=Distance
// Gd=Geodesic
// LL=Lat_Lon
// Sph=Spherical
//Geoid Data
EXPORT WGS84a_m:=6378137.0;
EXPORT WGS84a_km:=WGS84a_m/1000;
LOCAL WGS84fINV:=298.257223563;
EXPORT WGS84f:=1/WGS84fINV;
LOCAL WGS84b_m:=G_RadiusB(WGS84a_m,WGS84f);
LOCAL WGS84b_km:=WGS84b_m/1000;
LOCAL RA:=WGS84a_m; //RADIUS EQ
LOCAL FLATTEN:=WGS84f;
LOCAL RB:=G_RadiusB(RA,FLATTEN);
LOCAL RR:=G_SpMeanRadius(RA,RB);
LOCAL EPSILON_V:=1ᴇ−12;
LOCAL LoopSafe:=140;
//FORMAT:
//RA RB F
EXPORT G_XGNAME:=1;
EXPORT G_XRA:=2;
EXPORT G_XF:=4;
EXPORT G_EllABF:={
{"WGS84_m", WGS84a_m, 0, WGS84f},
{"Bessel_m", 6377397.155, 0, 1/299.1528128},
{"International_m", 6378388.000, 0, 1/297},
{"Sphere", 1,0,0},
{"Mercury_km", 2439.7, 2439.7, 0.0000},
{"Venus_km", 6051.8, 6051.8, 0.0000},
{"Earth WGS84_km", WGS84a_km, 0, WGS84f},
{"Moon_km", 1738.1, 1736.0, 0.0012},
{"Mars_km", 3396.2, 3376.2, 0.00589},
{"Jupiter_km", 71492, 66854, 0.06487},
{"Saturn_km", 60268, 54364, 0.09796},
{"Uranus_km", 22559, 24973, 0.02293},
{"Neptune_km", 24764, 24341, 0.01708},
{"Pluto_km", 1187, 1187, 0.0000},
{"Sun", 695508, 0, 0.00005},
{"Kepler 11145123 _km", 1.5ᴇ6, 1.5ᴇ6-3, G_Flattening(1.5ᴇ6,1.5ᴇ6-3)}};
// Kepler 11145123 is roundest object in nature as at 2016-11
EXPORT G_ThisEll:=G_EllABF(1);
EXPORT G_XLAT:=1;
EXPORT G_XLON:=2;
LOCAL ESS_UNITS:="m";
EXPORT G_About()
BEGIN
PRINT();
PRINT("Geodesy API V "+VER+" © 2016 StephenG1CMZ");
PRINT("NOT FOR NAVIGATIONAL USE");
PRINT("No liability accepted");
PRINT("Initial results suggest 14% error (Direct)");
END;
// STD MATHS
D2R(DEGS)
BEGIN
RETURN DEGS*(π/180);
END;
R2D(RADS)
BEGIN
RETURN RADS*(180/π);
END;
ATAN2YX(YY,XX)
BEGIN //MY ATAN2
IF XX==0 AND YY==0 THEN
RETURN 0; //CONVENIENT;
END;
//-π..π THE SPECIAL CHAR BELOW IS i=sqrt(−1), NOT CHINESE
RETURN ARG(XX+*YY);//BUILT-IN HANDLES MOST. END;
END;
HAV(THETA)
//TRIG: HAVERSINE=SIN^2(THETA/2)
BEGIN
LOCAL SN:=SIN(THETA/2);
RETURN SN^2;
END;
// END STD MATHS
EXPORT G_Flattening(RadiusA,RadiusB)
BEGIN
RETURN (RadiusA-RadiusB)/RadiusA;
END;
EXPORT G_RadiusB(RadiusA,Flatten)
//RETURN POLAR RADIUS
BEGIN
RETURN (1-Flatten)*RadiusA;//RADIUS POLAR
END;
EXPORT EPC(FNL,NN,EPCVAL)
//PRINT EXPECTED VAL AND ERR
BEGIN
LOCAL LST:=FNL;
PRINT("EPC Ret"+LST(NN));
PRINT("EPC Exp: "+HMS→(EPCVAL));
PRINT((100*(LST(NN)/EPCVAL))+"%");
END;
EXPORT G_DxVincentyDo(RA,RB,FLATTEN,Lat1,Lon1,Lat2,Lon2)
// GEODESY INVERSE PROBLEM ELLIPSOIDAL
// USE: RA AND FLATTEN FROM ThisEll
// RETURNS DIST IN SAME UNITS AS R
// LAT AND LON IN DEGREES
BEGIN
LOCAL NUM,DEN;
LOCAL LAT1:=HMS→(Lat1);
LOCAL LAT2:=HMS→(Lat2);
LOCAL LON1:=HMS→(Lon1);
LOCAL LON2:=HMS→(Lon2);
LOCAL UU1,SINUU1,COSUU1; //U1 AUX LATITUDE
LOCAL UU2,SINUU2,COSUU2; //U2 AUX LATITUDE
LOCAL USQ;
LOCAL AAA,BBB;
LOCAL LONDIFFradians; //L=DIFF IN LON
LOCAL PHI1,PHI2;
LOCAL LAMBDA,SINLAMBDA,COSLAMBDA,LAMBDAOLD;
LOCAL SIGMA,SIN_SIGMA,COS_SIGMA,DELTA_SIGMA,SINSQ_SIGMA; //AUX ARC LENGTH
LOCAL SIN_ALPHA,COSSQ_ALPHA;
LOCAL COS_2SIGMAM,COSSQ_2SIGMAM;
LOCAL CC;
LOCAL BRACKET_TERM;
LOCAL LOOPCOUNT;
LOCAL DISTS;
LOCAL AZ_FWD,AZ_FNL,AZ_REV;
// CALC L
LONDIFFradians:=D2R(LON2-LON1); //DIFF IN LON
// CALC U1 U2
PHI1:=D2R(LAT1);
PHI2:=D2R(LAT2);
UU1:=ATAN((1-FLATTEN)*TAN(PHI1));//CORRECT
UU2:=ATAN((1-FLATTEN)*TAN(PHI2));//CORRECT
//TO CALC LAMBDA:
// PRE-CALC U1 TRIG ONCEONLY
SINUU1:=SIN(UU1); COSUU1:=COS(UU1);
// PRE-CALC U2 TRIG EACH POSITION
SINUU2:=SIN(UU2); COSUU2:=COS(UU2);
//ITERATE CALCULATING LAMBDA
LOOPCOUNT:=0;
LAMBDA:=LONDIFFradians; //V13: DIFF IN LON ON AUX. FIRST APPROX
REPEAT
LOOPCOUNT:=LOOPCOUNT+1;
SINLAMBDA:=SIN(LAMBDA); COSLAMBDA:=COS(LAMBDA); //PRECALC
//V14
BRACKET_TERM:=(COSUU1*SINUU2 - SINUU1*COSUU2*COSLAMBDA);
SINSQ_SIGMA:=(COSUU2*SINLAMBDA)^2 + (BRACKET_TERM)^2;
IF SINSQ_SIGMA==0 THEN
//PRINT("SINσ 0: NO DISTANCE");
//0 DISTANCE MIGHT SUGGEST UNDERFLOW
RETURN {1,0,0,0,0,0};//VALID ZERO
END;
SIN_SIGMA:=√(SINSQ_SIGMA);
//V15
COS_SIGMA:=SINUU1*SINUU2 + COSUU1*COSUU2*COSLAMBDA;
//V16
SIGMA:=ATAN2YX(SIN_SIGMA,COS_SIGMA);
//V17
SIN_ALPHA:=(COSUU1*COSUU2*SINLAMBDA)/SIN_SIGMA;
COSSQ_ALPHA:=1-(SIN_ALPHA*SIN_ALPHA);
//V18
IF COSSQ_ALPHA==0 THEN
//PRINT("AVOID COSSQα 0");
COS_2SIGMAM:=0;
ELSE
COS_2SIGMAM:=COS_SIGMA - (2*SINUU1*SINUU2/COSSQ_ALPHA);
END;
CC:=(FLATTEN/16)*COSSQ_ALPHA*(4 + FLATTEN*(4 - 3*COSSQ_ALPHA));
LAMBDAOLD:=LAMBDA;
LAMBDA:=LONDIFFradians + (1-CC)*FLATTEN*SIN_ALPHA*(SIGMA + CC*SIN_SIGMA*(COS_2SIGMAM + CC*COS_SIGMA*(−1+2*(COS_2SIGMAM*COS_2SIGMAM))));
//PRINT(LAMBDA);
UNTIL (ABS(LAMBDA-LAMBDAOLD)≤EPSILON_V) OR (LOOPCOUNT≥LoopSafe);
//V19
USQ:=COSSQ_ALPHA*((RA^2 - RB^2)/(RB^2));
//(V3)
AAA:=1 + (USQ/16384)*(4096 + USQ*(−768 + USQ*(320 - 175*USQ)));
//(V4)
BBB:= (USQ/ 1024)*( 256 + USQ*(−128 + USQ*( 74 - 47*USQ)));
//(V6)
COSSQ_2SIGMAM:=(COS_2SIGMAM)^2;
//DELTA_SIGMA:=(1/6)*BBB*COS_2SIGMAM*(−3+4*(SIN_ALPHA)^2)*(−3+4*COSSQ_2SIGMAM) ;//V0.3
DELTA_SIGMA:=(1/6)*BBB*COS_2SIGMAM*(−3+4*(SIN_SIGMA)^2)*(−3+4*COSSQ_2SIGMAM) ; //V0.4
DELTA_SIGMA:=(1/4)*BBB*(COS_SIGMA*(−1+2*COSSQ_2SIGMAM)) - DELTA_SIGMA ;
DELTA_SIGMA:=BBB*SIN_SIGMA*(COS_2SIGMAM + DELTA_SIGMA);
//PRINT("Δσ "+DELTA_SIGMA);
DISTS:=AAA*(SIGMA - DELTA_SIGMA);//ANGULAR ≤π Radians
DISTS:=RB*DISTS; //same units as radius
//V20:FORWARD AZIMUTH
//CHECKED SINLAMBDA:ORIGINAL (LON2-LON1) OR LOOPED VALUE? USE LOOPED
AZ_FWD:=ATAN2YX((COSUU2*SINLAMBDA),BRACKET_TERM);
//V21:FINAL AND REVERSE AZIMUTH
NUM:=(COSUU1*SINLAMBDA);
DEN:=(−SINUU1*COSUU2 + COSUU1*SINUU2*COSLAMBDA);
AZ_FNL:=ATAN2YX(NUM,DEN);
AZ_REV:=ATAN2YX(−NUM,−DEN);
//RETURN RESULTS
DISTS:=ROUND(DISTS,3);//VINCENTY CLAIMS mm,
AZ_FWD:=R2D(AZ_FWD);
AZ_FWD:=IFTE(AZ_FWD<0,((AZ_FWD+360) MOD 360),AZ_FWD);//NO NEGATIVE BEARINGS
AZ_FNL:=R2D(AZ_FNL);
AZ_FNL:=IFTE(AZ_FNL<0,((AZ_FNL+360) MOD 360),AZ_FNL);//NO NEGATIVE BEARINGS
AZ_REV:=R2D(AZ_REV);
AZ_REV:=IFTE(AZ_REV<0,((AZ_REV+360) MOD 360),AZ_REV);
// VALIDITY: CONVERGENCE FAILURE
// DISTS: SAME UNITS AS RADIUS
// AZIMUTHS: POSITIVE 0..360
// LOOPCOUNT: CONVERGENCE INFO (LOWER BETTER):DIAGNOSTIC ONLY
// EXPECT SEQUENCE ORDER TO CHANGE
// ESPECIALLY AZIMUTHS
RETURN {(LOOPCOUNT≠LoopSafe),DISTS,(AZ_FWD),(AZ_FNL),(AZ_REV),LOOPCOUNT};
END;
EXPORT G_LL_VincentyDo(RA,RB,FLATTEN,Lat1,Lon1,GdDist,Bearing)
//DIRECT: GIVEN LL DIST BEAR, RETURN LL2
BEGIN
LOCAL USQ,NUM,DEN,TMP;
LOCAL LAT1:=HMS→(Lat1);
LOCAL LON1:=HMS→(Lon1);
LOCAL FNL,REV;
LOCAL SIGMA1;
LOCAL SIGMA,SIN_SIGMA,COS_SIGMA;
LOCAL TWOSIGMAM,TWOSIGMA1,COS_2SIGMAM,COSSQ_2SIGMAM;
LOCAL DELTA_SIGMA,SIGMAOLD;
LOCAL SIN_ALPHA;
LOCAL AAA,BBB,CCC,LLL,ALPHA2;
LOCAL LAMBDA,PHI2,LON2;
LOCAL S_BA;//CALC BEFORE ITERATING
LOCAL PHI1:= D2R(LAT1);//GEODETIC LAT
LOCAL LAMBDA1:=D2R(LON1);
LOCAL ALPHA1:=D2R(Bearing);//INITIAL BEARING
LOCAL SIN_ALPHA1:=SIN(ALPHA1);
LOCAL COS_ALPHA1:=COS(ALPHA1);
LOCAL SINSQ_ALPHA,COSSQ_ALPHA;
PRINT("In:");
PRINT({RA,RB,FLATTEN});
PRINT({Lat1,Lon1,GdDist,Bearing});
//PRINT("DEBUG:DO NOT USE!");
//PRINT("Some VincentyDirect results wrong!");WAIT(0.1);
//Karney OK. Australia>100km err
//CALCULATE UU1
LOCAL TANUU1:=(1-FLATTEN)*TAN(PHI1);//CORRECT
LOCAL UU1:=ATAN(TANUU1);//CORRECT
//LOCAL COSUU1:=COS(UU1);
LOCAL COSUU1:=1/√(1+TANUU1^2);//TRY
//LOCAL SINUU1:=SIN(UU1);// = COS*TAN SAVES A TRIG
LOCAL SINUU1:=COSUU1*TANUU1;//TRY
IF GdDist==0 AND 1 THEN
RETURN {LAT1,LON1,0};//AVOID OR REVEAL ROUNDING ERROR
END;
// V1
//IF EQUATORIAL (PHI1=0 AND BEARING=π): USE SIGMA1=0
SIGMA1:=ATAN2YX(TANUU1,COS_ALPHA1);
// V2
SIN_ALPHA:=COSUU1*SIN_ALPHA1;
SINSQ_ALPHA:=(SIN_ALPHA)^2;
//V3 V4
COSSQ_ALPHA:=1-SINSQ_ALPHA;
USQ:=COSSQ_ALPHA*((RA^2 - RB^2)/(RB^2));
//(V3)
AAA:=1 + (USQ/16384)*(4096 + USQ*(−768 + USQ*(320 - 175*USQ)));
//(V4)
BBB:= (USQ/ 1024)*( 256 + USQ*(−128 + USQ*( 74 - 47*USQ)));
//ITERATE
S_BA:=(GdDist/(RB*AAA));
SIGMA:=S_BA;//1ST EST
REPEAT
// V5
TWOSIGMAM:=2*SIGMA1 + SIGMA;
//V6 PREP
COS_2SIGMAM:=COS(TWOSIGMAM);
SIN_SIGMA:=SIN(SIGMA);
COS_SIGMA:=COS(SIGMA);
//(V6)
COSSQ_2SIGMAM:=(COS_2SIGMAM)^2;
DELTA_SIGMA:=(1/6)*BBB*COS_2SIGMAM*(−3+4*(SIN_SIGMA)^2)*(−3+4*COSSQ_2SIGMAM) ;// V0.4
DELTA_SIGMA:=(1/4)*BBB*(COS_SIGMA*(−1+2*COSSQ_2SIGMAM)) - DELTA_SIGMA ;
DELTA_SIGMA:=BBB*SIN_SIGMA*(COS_2SIGMAM + DELTA_SIGMA);
//PRINT("Δσ "+DELTA_SIGMA);
// V7
SIGMAOLD:=SIGMA;
SIGMA:=S_BA+DELTA_SIGMA;
UNTIL ABS(SIGMA-SIGMAOLD)≤EPSILON_V;
//UPDATE?
//SOME SOURCES OMIT AND USE OLDER VALUES
//SO MAY INHIBIT FOR COMATIBILITY TESTS
IF 1 THEN
SIN_SIGMA:=SIN(SIGMA);
COS_SIGMA:=COS(SIGMA);
END;
//V8 LATITUDE
TMP:=(SINUU1*SIN_SIGMA - COSUU1*COS_SIGMA*COS_ALPHA1);
NUM:=(SINUU1*COS_SIGMA) + (COSUU1*SIN_SIGMA*COS_ALPHA1);
DEN:=((SIN_ALPHA)^2+(TMP)^2);
DEN:=(1-FLATTEN)*√(DEN);
PHI2:=ATAN2YX(NUM,DEN);
//V9 LAMBDA FOR LON
NUM:=SIN_SIGMA*SIN_ALPHA1;
DEN:= (COSUU1*COS_SIGMA) - (SINUU1*SIN_SIGMA*COS_ALPHA1);
LAMBDA:=ATAN2YX(NUM,DEN);
//V10 CCC FOR LON
CCC:=(FLATTEN/16)*COSSQ_ALPHA*(4+FLATTEN*(4-3*COSSQ_ALPHA));
//V11 LLL:=DIFF IN LON
LLL:=(SIGMA + CCC*SIN_SIGMA*(COS_2SIGMAM + CCC*COS_SIGMA*(−1 + 2*COSSQ_2SIGMAM)));
LLL:=LAMBDA - (1 - CCC)*FLATTEN*SIN_ALPHA*LLL;
//V12 ALPHA2:=GEODETIC LATITUDE
//DEN:=(−SINUU1*SIN_SIGMA + COSUU1*COS_SIGMA*COS_ALPHA1);
FNL:=ATAN2YX(SIN_ALPHA,−TMP);
REV:=ATAN2YX(−SIN_ALPHA,TMP);
// RETURN
// PHI2: GEODETIC LATITUDE −90..90
// LONGITUDE: NORMALISE TO −180..180
// BEARING2: FINAL BEARING: NORMALISE TO 0..360
// NORMALISATIONS TBD
PHI2:=R2D(PHI2);
LON2:=LON1+R2D(LLL);
// ALPHA2 == FNL BEARING
FNL:=R2D(FNL);
FNL:=IFTE(FNL<0,((FNL+360) MOD 360),FNL);//NO NEGATIVE BEARINGS
REV:=R2D(REV);
REV:=IFTE(REV<0,((REV+360) MOD 360),REV);//NO NEGATIVE BEARINGS
PRINT("Out:");
PRINT({PHI2,LON2});
RETURN {PHI2,LON2,REV,FNL};
END;
EXPORT G_DxVincenty(Elevation,Point1,Point2)
// GEODESY INVERSE PROBLEM ELLIPSOIDAL
// USE: RA AND FLATTEN FROM ThisEll TO DERIVE RB
// RETURNS DIST IN SAME UNITS AS R
// LAT AND LON IN DEGREES
BEGIN
LOCAL LAT1:=Point1(G_XLAT);
LOCAL LON1:=Point1(G_XLON);
LOCAL LAT2:=Point2(G_XLAT);
LOCAL LON2:=Point2(G_XLON);
LOCAL RA:=G_ThisEll(G_XRA)+Elevation;
LOCAL FLATTEN:=G_ThisEll(G_XF);
LOCAL RB:=G_RadiusB(RA,FLATTEN);
RETURN G_DxVincentyDo(RA,RB,FLATTEN,LAT1,LON1,LAT2,LON2);
END;
EXPORT G_LL_Vincenty(Elevation,Point1,GdDist,Bearing)
// GEODESY DIRECT PROBLEM ELLIPSOIDAL
// USE: RA AND FLATTEN FROM ThisEll TO DERIVE RB
// Point1: Lat Lon
// GdDistance: Geodesic Distance
// Bearing:
// LAT AND LON IN DEGREES
BEGIN
LOCAL LAT1:=Point1(G_XLAT);
LOCAL LON1:=Point1(G_XLON);
LOCAL RA:=G_ThisEll(G_XRA)+Elevation;
LOCAL FLATTEN:=G_ThisEll(G_XF);
LOCAL RB:=G_RadiusB(RA,FLATTEN);
RETURN G_LL_VincentyDo(RA,RB,FLATTEN,LAT1,LON1,GdDist,Bearing);
END;
G_UseThis(POSN)
//VINCENTY RECALCULATES THESE==DUPLICATION
BEGIN
IF POSN>0 THEN
G_ThisEll:=G_EllABF(POSN);
RA:=G_ThisEll(G_XRA);
FLATTEN:=G_ThisEll(G_XF);
RB:=G_RadiusB(RA,FLATTEN);
RR:=G_SpMeanRadius(RA,RB);
//ELSE.UNCHANGED
END;
RETURN POSN;
END;
EXPORT G_UseN(POSN)
//N==0: USER CHOOSE
//N>0: CHOOSE ENTRY N
BEGIN
LOCAL OKC;
LOCAL NN:=POSN;
IF NN==0 OR NN>SIZE(G_EllABF) THEN
OKC:=CHOOSE(NN,"Ellipsoid",G_EllABF);
IF OKC THEN
MSGBOX("DETAILS TBD");
RETURN G_UseN(NN);
END;
END;
RETURN G_UseThis(NN);
END;
EXPORT G_UseNew(LST)
//CREATE NEW ELLIPSOID (NO CHECK FOR DUPLICATES)
//LST FORMAT AS IN EXAMPLES
//IDEALLY RETURN 0 IF LIST TOO BIG
BEGIN
G_EllABF(0):=LST;
RETURN G_UseThis(SIZE(G_EllABF));
END;
EXPORT G_UseST(ST)
//SEARCH FOR STRING
//IF FOUND USE THAT ELLIPSOID
//A RADIUS EQU
//B RADIUS POL
//AT PRESENT RB IS ALWAYS DERIVED FROM F
//LATER ALLOW F TO BE DERIVED
BEGIN
LOCAL II,POS;
FOR II FROM 1 TO SIZE(G_EllABF) DO
POS:=INSTRING(G_EllABF(G_XGNAME),ST);
IF POS THEN
G_UseThis(POS);
RETURN POS;
END;
END;
RETURN POS;//0=>UNCHANGED
END;
G_Use(PP)
//NOT YET WORKING-(USING MENU RUN ROUTINE INTERACTIVELY)
//PROGRAM LOGIC LOOKS OK BUT TYPE HANDLING RUNTIME ERR
//INTENDED TO PROVIDE A SINGLE EXPORT TO REPLACE 3 EXPORT ROUTINES
//WHICH COULD THEN BE MADE LOCAL
//(NUM): SELECT ELLIPSOID N FROM LIST (0=SHOW CHOOSE LIST O USER)
// EG 1
//(STR): SELECT 1st ELLIPSOID CONTAINING STRING
// EG "WGS84" or "WGS84_m"
//(LST): CREATE (AND USE) NEW ELLIPSOID
// EG {"My Planet",1,0,0.1}
BEGIN
//MSGBOX(TYPE(PP));
CASE
IF TYPE(PP)==0 THEN //NUM
RETURN G_UseN(PP);
END;
IF TYPE(PP)==2 THEN //STRING:SEARCH FOR THAT ELLIPSOID NAME
RETURN G_UseST(PP);
END;
IF TYPE(PP)==6 THEN //LIST: CREATE AND USE NEW ELLIPSOID
RETURN G_UseNew(PP);
END;
DEFAULT
RETURN MSGBOX("G_Use ERR: NUM OR STR OR LST");
END;
END;
// Main Spherical Geodesy routines-
EXPORT G_SpMeanRadius(RA,RB)
//In Geodesy, this is usually used as a radius for spherical
BEGIN
RETURN (2*RA+RB)/3;
END;
EXPORT G_SpDistHav(Point1,Point2)
BEGIN
LOCAL LAT1:=HMS→(Point1(G_XLAT));
LOCAL LAT2:=HMS→(Point2(G_XLAT));
LOCAL LON1:=HMS→(Point1(G_XLON));
LOCAL LON2:=HMS→(Point2(G_XLON));
LOCAL DLON:=D2R(LON2-LON1);
LOCAL DLAT:=D2R(LAT2-LAT1);
LOCAL AA:=HAV(DLAT) + (COS(D2R(LAT1))*COS(D2R(LAT2))*HAV(DLON));
LOCAL CC:=2*ASIN(MIN(1,√(AA)));//MIN 1 AVOIDS ROUNDING OUTRANGE
LOCAL DD:=RR*CC;
RETURN ROUND(DD,−6);
END;
EXPORT G_LL_Sp(Elevation,Point1,GdDist,Bearing)
//Sph Direct LL1 Dist to LL2
//https://uk.answers.yahoo.com/question/index?qid=20161205044438AAMKcRG
//BASED ON MOVABLE TYPE
BEGIN
LOCAL LAT1:=HMS→(Point1(1));
LOCAL LON1:=HMS→(Point1(2)); LOCAL ALPHA1:=D2R(Bearing);
LOCAL ALPHA:=D2R(Bearing);
LOCAL COSALPHA:=COS(ALPHA);
LOCAL DELTA:=GdDist/RR;
LOCAL PHI1:=D2R(LAT1);
LOCAL SINPHI1:=SIN(PHI1);
LOCAL COSPHI1:=COS(PHI1);
LOCAL SINDELTA:=SIN(DELTA);
LOCAL COSDELTA:=COS(DELTA);
LOCAL SINPHI2:=(SINPHI1*COSDELTA + COSPHI1*SINDELTA*COSALPHA);
LOCAL PHI2:=ASIN(SINPHI2);
LOCAL LAT2=R2D(PHI2);
LOCAL NUM:=SIN(ALPHA1)*SIN(GdDist/RR)*COSPHI1;
LOCAL DEN:=COS(GdDist/RR) - (SINPHI1*SINPHI2);
LOCAL LAMBDA:=ATAN2YX(NUM,DEN);
LOCAL LON2:=LON1+R2D(LAMBDA);
IF Elevation≠0 THEN
MSGBOX("TBD");
END;
PRINT("Sph In:");
PRINT(RR);
PRINT({Point1,GdDist,Bearing});
//NORMALISE LAT2 LON2 TBD
PRINT("Sph Out: "+{LAT2,LON2});
RETURN {LAT2,LON2};
END;
// End Main Spherical Routines
// Useful routines
EXPORT G_SpSurfaceArea(Radius)
// RETURN SPHERICAL SURFACE AREA IN SAME UNITS AS RR
BEGIN
RETURN 4*π*(Radius^2);
END;
//
EXPORT G_FaceArea(RadiusA,RadiusB)
BEGIN
RETURN π*RadiusA*RadiusB;
END;
G_EllStats(PR)
//PR:PRINT ELSE TEXTOUT
BEGIN
LOCAL II;
LOCAL ST:={};
ST(1):=G_ThisEll(G_XGNAME);
ST(2):="RA: "+RA;
ST(3):="RB: "+RB;
ST(4):="f: "+(FLATTEN);
IF FLATTEN≠0 THEN
ST(5):="1/f: "+(1/FLATTEN);
END;
IF PR THEN
FOR II FROM 1 TO SIZE(ST) DO
PRINT(ST(II));
END;
ELSE
FOR II FROM 1 TO SIZE(ST) DO
TEXTOUT_P(ST(II),0,20*II+100);
END;
END;
END;
G_EX(LAT1,LON1,LAT2,LON2)
//An example displaying results
BEGIN
LOCAL TM,RESULT1,RESULT2,DISTS;
PRINT(" ");
PRINT("From LAT "+HMS→(LAT1)+" "+→HMS(LAT1));
PRINT("From LON "+HMS→(LON1)+" "+→HMS(LON1));
PRINT("To LAT "+HMS→(LAT2)+" "+→HMS(LAT2));
PRINT("To LON "+HMS→(LON2)+" "+→HMS(LON2));
TM:=TICKS;
RESULT1:=G_SpDistHav({LAT1,LON1},{LAT2,LON2});
RESULT2:=G_DxVincenty(0,{LAT1,LON1},{LAT2,LON2});
TM:=TICKS-TM;
PRINT("Elapsed "+(TM/1000)+" s");
PRINT("Spherical Haversine: "+(RESULT1)+" m");
PRINT("Ellipsoidal Vincenty: "+(RESULT2(2))+" m");
PRINT(" AZ "+(→HMS(RESULT2(3))));
PRINT(" AZ "+(→HMS(RESULT2(4))));
PRINT(" AZ "+(→HMS(RESULT2(5))));
PRINT("Validity: " +(RESULT2(1))+" "+(RESULT2(6)));
//PRINT(RESULT2);
//RESULT:=Z_GEODESY(LAT1,LON1,LAT2,LON2);
//WAIT;
END;
EXPORT G_ExampleUse()
//A short example of key procedures
BEGIN
LOCAL OKC,DIST,RESULT,EXAMPLEPOINT,XLAT,XLON;
OKC:=G_UseNew({"Add your own",{RA,0,FLATTEN}});//Insert into List. Use use to use it.
OKC:=G_UseST("WGS84_m");//DEFAULT
IF OKC THEN //Found in list
DIST:=G_SpDistHav({1,2},{3,4});
//RESULT:=G_DxVincentyDo(RA,RB,FLATTEN,LAT1,LON1,LAT2,LON2);
RESULT:=G_DxVincenty(0,{1,2},{3,4});
//Variables prefixed X index into data structures
EXAMPLEPOINT:={"From Here",1,2};
XLAT:=2; XLON:=3;
XLAT:=1; XLON:=2; //REINSTATE CURRENT DEFAULT
EXAMPLEPOINT:={1,2};
//and similarly for geoid structure
END;
END;
EXPORT G_DrawEll(Flatten)
//The stats are omitted because flattening drawn may differ
//makes illustrate changing f easier
//for now use myRadiusA of 100 pixels
BEGIN
LOCAL XX,YY;
LOCAL DX,DY;
LOCAL MyRadiusA:=100;
LOCAL HEIGHT:=100;
LOCAL WIDTH:=100;
//LOCAL SCALE:=RadiusA/WIDTH;
LOCAL EWIDTH:=MyRadiusA;
LOCAL EHEIGHT:=((1-Flatten)*MyRadiusA);//RADIUS POLAR
RECT();
//G_EllStats(0);
FOR YY FROM −HEIGHT TO HEIGHT DO
//PRINT(YY);
FOR XX FROM −WIDTH TO WIDTH DO
DX:=XX/(EWIDTH);
DY:=YY/(EHEIGHT);
IF (DX^2+DY^2)≤1 THEN
//PRINT({WIDTH+10+XX,HEIGHT+20+YY});
PIXON_P(WIDTH+100+XX,HEIGHT+20+YY,#009900h);
END;//IF
END;//FOR
END;//FOR
//G_EllStats(0);
WAIT;
END;
G_TEST()
//A few tests. More needed.
BEGIN
G_EX(50°03′58.76″,−005°42′53.10″,58°38′38.48″,−003°04′12.34″);//50PERCENT
PRINT("EXPECT 969 954_m");
G_EX(−37°57′3.72030″,144°25′29.52440″,−37°39′10.15610″,143°55′35.38390″);
PRINT("EXPECTED: 54972.271m");//VINCENTY AUSTRALIA
PRINT("306d52 05 * 127d10 25 * 307d10 25.07 ");
G_EX(31.8300167,35.0662833,31.83000000,35.0708167);
PRINT("PYTHON: CIRCLE 428.4 VIN 429.1677 m");
//G_EX(0,0,0.5,179.5);
G_EX(0,0,0.5,179.7);//WILL NOT CONVERGE
//G_EX(0,0,0,90);
//G_EX(0,0,0.5,0.5);
G_EX(2,2,2,2);
G_EX(50°03′58.76″,−005°42′53.10″,58°38′38.48″,−003°04′12.34″);//50PERCENT
PRINT("EXPECT 969 954_m");
//IF 0 THEN
G_EX(50°3′59″,−5°42′53″,58°38′38″,−3°4′12″);
PRINT("RECHECK: EXPECTED 968.9km 009DEG0711");
IF 1 THEN //DIRECT TESTS
PRINT("DIRECT TESTS");
//KARNEY 2012 EG IS WITHIN SECONDS
//Karney, C.F.F. J Geod (2013) 87: 43. doi:10.1007/s00190-012-0578-z
PRINT("Karney Ell");
PRINT(G_LL_Vincenty(0,{40,0},10000*1000,30));//KARNEY 2012 EXAMPLE (LON ASSUMED)
PRINT((G_LL_Vincenty(0,{40,0},10000*1000,30)));//KARNEY 2012 EXAMPLE (LON ASSUMED)
PRINT("EXPECT:");
PRINT(→HMS({41.79331020506,137.84490004377,149.09016931807}));//EXPECTED
PRINT({41.79331020506,137.84490004377,149.09016931807});//EXPECTED
PRINT("Znd");
EPC(G_LL_Vincenty(0,{40,0},1000*1000,30),1,41.79331020506);
PRINT("K Sph");
PRINT(G_LL_Sp(0,{40,0},1000*1000,30));
EPC(G_LL_Sp(0,{40,0},1000*1000,30),1,41.79331020506);
//A 0 DISTANCE GIVES A ROUNDING ERR UNLESS 0 TRAPPED
PRINT("ZERO");
PRINT(G_LL_Vincenty(0,{45,30},0,0));
//AUSTRALIA INV EG IS WRONG
PRINT("Au:");
PRINT(→HMS(G_LL_Sp(0,{−37°57′3.72030″,144°25′29.52440″},54972.271,306°52′5.37″)));//AU
EPC(G_LL_Sp(0,{−37°57′3.72030″,144°25′29.52440″},54972.271,306°52′5.37″),1,−37°39′10.15610″);//AU
PRINT(→HMS(G_LL_Vincenty(0,{−37°57′3.72030″,144°25′29.52440″},54972.271,306°52′5.37″)));//AU
EPC(G_LL_Vincenty(0,{−37°57′3.72030″,144°25′29.52440″},54972.271,306°52′5.37″),1,−37°39′10.15610″);//AU
PRINT("EXPECTED:"+→HMS({−37°39′10.15610″,143°55′35.38564″,127°10′25.07″,307°10′2507″}));//OR306TBC
END;
IF 0 THEN
G_EX(−37°57′3.72030″,144°25′29.52440″,−37°39′10.15610″,143°55′35.38390″);//EXPECTS 54972m
PRINT("EXPECTED: 54972.271m ");
PRINT("306d52:5.37 / 127d10:25.07");//VINCENTY AUSTRALIA
//IF 1 THEN
G_EX(31.8300167,35.0662833,31.83000000,35.0708167);
PRINT("PYTHON: CIRCLE 428.4 VIN 429.1677 m");
//IF 0 THEN
G_EX(0,0,0.5,179.5);
G_EX(0,0,0.5,179.7);//WILL NOT CONVERGE
//IF 1 THEN
G_EX(0,0,0,90);
G_EX(0,0,0.5,0.5);
G_EX(2,2,2,2);
END;
END;
EXPORT GEODESY()
BEGIN
G_About();
G_EllStats(1);
PRINT("Esc");WAIT;
IF MSGBOX("Run Test?") THEN
G_TEST();
END;
END;
[/u]