HP Forums
Geodesy API - Printable Version

+- 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;
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.


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.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]