Geodesy API
10-17-2016, 09:49 PM
Post: #1
 StephenG1CMZ Senior Member Posts: 918 Joined: May 2015
Geodesy API
A collection of useful geodesy routines.
Disclaimer: Not for navigational use. No liability accepted.

Stephen Lewkowicz (G1CMZ)
10-17-2016, 10:11 PM (This post was last modified: 10-23-2016 06:05 PM by StephenG1CMZ.)
Post: #2
 StephenG1CMZ Senior Member Posts: 918 Joined: May 2015
RE: Geodesy API
Anyone using Version 0.1 or Version 0.2 should now replace it.

Stephen Lewkowicz (G1CMZ)
10-23-2016, 06:12 PM (This post was last modified: 11-22-2016 06:37 PM by StephenG1CMZ.)
Post: #3
 StephenG1CMZ Senior Member Posts: 918 Joined: May 2015
RE: Geodesy API
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.

Stephen Lewkowicz (G1CMZ)
11-22-2016, 06:34 PM (This post was last modified: 11-22-2016 09:42 PM by StephenG1CMZ.)
Post: #4
 StephenG1CMZ Senior Member Posts: 918 Joined: May 2015
RE: Geodesy API
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.

Stephen Lewkowicz (G1CMZ)
11-23-2016, 02:09 PM (This post was last modified: 12-14-2016 09:30 AM by StephenG1CMZ.)
Post: #5
 StephenG1CMZ Senior Member Posts: 918 Joined: May 2015
RE: Geodesy API
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;
Version 0.31 worked on my calculator but used code in backups]...After deleting backups once the code was here, it didn't, and was replaced by V 0.32.

A compiler option to ignore backups would be useful.

Stephen Lewkowicz (G1CMZ)
11-23-2016, 03:03 PM
Post: #6
 StephenG1CMZ Senior Member Posts: 918 Joined: May 2015
RE: Geodesy API
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;

Stephen Lewkowicz (G1CMZ)
12-13-2016, 08:57 PM
Post: #7
 StephenG1CMZ Senior Member Posts: 918 Joined: May 2015
RE: Geodesy API
[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.3​7″),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]

Stephen Lewkowicz (G1CMZ)
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)