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