Post Reply 
Geodesy API
12-13-2016, 08:57 PM
Post: #7
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)
https://my.numworks.com/python/steveg1cmz
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
Geodesy API - StephenG1CMZ - 10-17-2016, 09:49 PM
RE: Geodesy API - StephenG1CMZ - 10-17-2016, 10:11 PM
RE: Geodesy API - StephenG1CMZ - 10-23-2016, 06:12 PM
RE: Geodesy API - StephenG1CMZ - 11-22-2016, 06:34 PM
RE: Geodesy API - StephenG1CMZ - 11-23-2016, 02:09 PM
RE: Geodesy API - StephenG1CMZ - 11-23-2016, 03:03 PM
RE: Geodesy API - StephenG1CMZ - 12-13-2016 08:57 PM



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