Horizons: Distance to Horizon
04-11-2019, 09:10 PM (This post was last modified: 04-15-2019 09:25 PM by StephenG1CMZ.)
Post: #3
 StephenG1CMZ Senior Member Posts: 918 Joined: May 2015
RE: Horizons: Distance to Horizon
Version 0.2:
Improves cas solve implementation (using h instead of HT)
Improves handling of extreme values
Implements inverse surface area calculations (although these may be approximate and need more checking).
Note: The I at the start of the code is a cut-and-paste error - please delete it.
Code:
 // If you have previously downloaded this, please delete this I: // I   LOCAL CRID:="Horizons V0.2 © 2019 StephenG1CMZ";  //This version has simple geometric horizons  //No refraction    LOCAL HT;//FOR SOLVE //Global h used    LOCAL AU2KM:=149597871;//GOOGLE  LOCAL RADIUS:=6378140;//m //MUST BE >=0  //MOST OF MY TEST VALUES ASSUME THIS VALUE. CHOOSE YOUR OWN.  //EARTH:WORKING IN m ACCURACY IS NORMALLY WITHIN 10cm (OF LATEST DATA)  LOCAL NA:=0;  // NA OR 0 => ONE OF:   //INPUT<0: NOTHING VISIBLE   //INPUT=0: NOTHING VISIBLE   //INPUT>0.5:OVER 1/2 OF SPHERE CANNOT BE VISIBLE   //INPUT=0.5 (INFINITE HEIGHT LIMITING CASE, use UpperLimits for values)   //AVOID /0 (HOPEFULLY ONE OF THE ABOVE)  //Definitions  //DD   Direct Distance       (eye-horizon)  //HT   Height above Sphere   (<0=NA)  //SLen Curved Surface Length (paw-horizon)  //Distances and heights in same units    GetSphereSA(RADIUS)  //RETURN SPHERE SURFACE AREA  BEGIN   RETURN 4*π*RADIUS^2;  END;  EXPORT GetDLen(HT,RADIUS) //Direct distance  //RETURN D is the straight geometrical eye-hZ dist  //TIMING 58us  BEGIN    IF HT*RADIUS<0 THEN     RETURN NA;    END;    RETURN √(HT*(2*RADIUS+HT));  END;  EXPORT GetHeightDD(DD,RADIUS)  //INVERSE  //RETURN HT TO SEE DD AWAY  BEGIN   IF DD<0 THEN     //NO UPPER LIMIT TO EYE HEIGHT OR DIST     //(Imp Limit DD>1ᴇ12 WITH EARTH RADIUS NOT CHECKED)    RETURN NA;   END;   HT:=CAS("solve((h*(2*RADIUS+h))-DD^2=0,h,h≥0)");   RETURN HT;  END;  EXPORT GetHeightSAF(SAF,RADIUS)  //INVERSE  //Get Height for Surface Area (SAF Fraction of sphere visible, 0..0.5)  //OR NA IF IMPOSSIBLE  BEGIN    IF SAF≤0 OR SAF≥0.5 THEN //INCLUDE SAF=0 TRAP     RETURN NA;    END;    HT:=CAS("solve(h/(2*(h+RADIUS))-SAF=0,h,h≥0)");    RETURN HT; //MAX(HT,0);     //MAX CLIPS SMALL NEGATIVES: NEEDED WHEN SAF=0 UNLESS TRAPPED  END;  //NB INVERSE SA RESULTS ARE APPROXIMATE  //BUT NEED CHECKING TBD  EXPORT GetHeightSA(SA,RADIUS)  //INVERSE  //Get Height for Surface Area (SA in units^2)  //OR NA IF IMPOSSIBLE  BEGIN   LOCAL SAF;   IF RADIUS==0 THEN    RETURN NA;   END;   SAF:=SA/GetSphereSA(RADIUS);   PRINT(100*SAF);//diagnostic   HT:=GetHeightSAF(SAF,RADIUS);   //HT:=RADIUS*((SA*SAF)-1);   RETURN HT;  END;  EXPORT GetHeightSL(SurfLen,RADIUS)  //INVERSE  //RETURN HT TO SEE SURFLEN AWAY (OR NA IF SURFLEN IMPOSSIBLE)  BEGIN   IF SurfLen<0 OR SurfLen≥π*RADIUS THEN    RETURN NA;   END;   HT:=(RADIUS/COS(SurfLen/RADIUS))-RADIUS;   RETURN HT;   //RETURN {(RADIUS/COS(SD/RADIUS))-RADIUS,CAS("solve((RADIUS/COS(SD/RADIUS))-RADIUS-HT=0,HT)")};  END;  EXPORT GetSurfLenDD(DD,RADIUS) //Surface Length from DD  //D is the straight geometrical eye-hor dist  //RETURN S IS THE SURFACE DISTANCE   //TIMING 55us  BEGIN    IF RADIUS==0 THEN     RETURN NA;    END;    RETURN RADIUS*ATAN(DD/RADIUS);  END;  //BOTH THESE S DISTANCES SHOULD BE THE SAME (EXCEPT ROUNDING)  //THERE IS NO SIGNIFICANT DIFFERENCE  //HINT:   //S(HT) IS QUICKER IF YOU ONLY NEED S  //S(DD) IS QUICKER IF YOU ALSO NEED D (OR KNOW D)  EXPORT GetSurfLenHT(HT,RADIUS)  //RETURN S IS THE SURFACE DISTANCE  //TIMING 109us QUICKER THAN D+SD  BEGIN    IF HT*RADIUS<0 OR (RADIUS+HT==0) THEN      RETURN NA;    END;    RETURN RADIUS*ACOS(RADIUS/(RADIUS+HT));  END;  //Surfacevisible at Height  //FP Fraction  //PC Percent  //SA Surface Area (same units as radius)  EXPORT GetSurfAreaHTFP(HT,RADIUS)  //Derived From: http://www.hpmuseum.org/forum/thread-12717.html?highlight=Spher  //Height and Radius in same units. Ht is above surface   //Return.Proportion Surface Visible at Ht (0-0.5,HT≥0)  //For example:  //Radius of Earth = 6378 km  //Distance of Apollo 17 from Earth's surface at 1972-12-07T10:39Z = 29000 km  //Visible percentage of Earth's surface in 1972 "Blue Marble" photograph = 40.99% (41%)  BEGIN   //Guard HT<0:WHAT IS REALLY VISIBLE   IF HT<0 OR (HT+RADIUS==0)THEN    RETURN NA;//0 OR 1   END;   RETURN (HT/(2*(HT+RADIUS)));  END;  GetSurfAreaHTPC(HT,RADIUS)  BEGIN   RETURN 100*GetSurfAreaHTFP(HT,RADIUS);  END;  EXPORT PrintReport(HT,RADIUS)  //User-Readable-Report  //Hint: See EX   BEGIN   LOCAL PLACES:=3;//ROUNDING   LOCAL DX:=GetDLen(HT,RADIUS);   LOCAL FP:=GetSurfAreaHTFP(HT,RADIUS);   PRINT("At Height "+HT+" Radius "+RADIUS);   PRINT("Geometric visible horizon: "+ROUND(DX,PLACES));   PRINT("Surface horizon: "+ROUND(GetSurfLenDD(DX,RADIUS),PLACES));   PRINT("Surface Area: "+ROUND(FP*GetSphereSA(RADIUS/1000),PLACES));  END;  EXPORT EX(HT,RADIUS)  //Work an example  //Exhibits program functionality  //Hint: See Report  BEGIN   LOCAL DX:=GetDLen(HT,RADIUS);   LOCAL FP:=GetSurfAreaHTFP(HT,RADIUS);   PRINT("HT "+HT+" R "+RADIUS);   PRINT("Direct Length "+DX);   PRINT("SurfaceLenDD "+GetSurfLenDD(DX,RADIUS));   PRINT("SurfaceLenHT "+GetSurfLenHT(HT,RADIUS));   PRINT("SurfArea "+STRING({100*FP+"%",FP*GetSphereSA(RADIUS)}));   //PRINT("Verify: HT SD "+GetHeightSL(GetSurfLenHT(HT,RADIUS),RADIUS));   PRINT("Verify: HT DD "+GetHeightDD(GetDLen     (HT,RADIUS),RADIUS));   PRINT("Verify: HT SAF"+GetHeightSAF(FP,RADIUS));   //PRINT("Verify: HT SA "+GetHeightSA(SA,RADIUS));   PRINT(" ");  END;  EXPORT UpperLimits(RADIUS)  //Upper Limits = Half Sphere = Infinite Height  //Appply to SurfLen,SurfArea. Not DD and HT  BEGIN   LOCAL SpC,SpA;    SpC:=π*RADIUS;    SpA:=2*π*RADIUS^2;    RETURN {SpC,SpA};//HALF SPHERE  END;  EXAMPLES()   BEGIN   // a feW examples : most test data elsWhere   //HT,D,S (D=STRAIGHT S=CURVED)   PRINT(); PRINT(CRID); PRINT(" ");WAIT(1);   PRINT("Sphere: R "+RADIUS/1000+" km. Circ & Area:");   PRINT(2*UpperLimits(RADIUS/1000));//km   PRINT("Upper limits: SurfLen & SurfArea");   PRINT(STRING(UpperLimits(RADIUS))+" m,m^2");//m   PRINT(" ");   //HT≤0 RETURNS 0    //EX(−1,RADIUS);   EX(0,RADIUS);   PRINT("//MAN at 2m (D 5.1km)");   EX(2,RADIUS);   PRINT("//ISS at 408 km (around 2373 km)");   EX(427200,RADIUS);   PRINT("//Eg APOLLO at 29000 km");   EX(29000ᴇ3,6378000);   PRINT("Verify+:"+GetHeightSA(2.09514275019ᴇ14,RADIUS));   PRINT("Note the differences in these calculations...");   PRINT("My inverse surface area may be inaccurate");   PRINT(" ");   PRINT("//40.5 AU Voyager PALE BLUE DOT");   EX(40.5*AU2KM,RADIUS);   //PRINT("//145 AU Voyager 2019");   //EX(145*AU2KM,RADIUS);   PRINT("//Viewed as Exoplanet from Proxima Centauri");   EX(268395.131*AU2KM,RADIUS);   //This last verifies, except solve doesnt      //Test Traps    PRINT(GetHeightDD(2ᴇ12,RADIUS));   PRINT(GetHeightSL(π*RADIUS,RADIUS));   PRINT(GetHeightSAF(0.50,RADIUS));//undef unless trapped   //PRINT("READY");  END;  EXPORT HORIZONS()  BEGIN   LOCAL DX;   LOCAL HT:=2;   PRINT();   PRINT(CRID);   EXAMPLES();   //PRINT("TEVAL DStraight "+TEVAL(DX:=HZDistanceD(HT,RADIUS)));// 58 CALC D    //PRINT("TEVAL SurfaceD "+TEVAL(HZDistanceSD(DX,RADIUS)));    // 54 S ADD TO ABOVE    //PRINT("TEVAL SurfaceH "+TEVAL(HZDistanceSH(HT,RADIUS)));    //109 VS S   END;
Note: Some of my examples are incorrect - I mangled km and m. In particular, the AU2KM AU2M conversion.

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

 Messages In This Thread Horizons: Distance to Horizon - StephenG1CMZ - 04-09-2019, 05:55 PM RE: Horizons: Distance to Horizon - StephenG1CMZ - 04-09-2019, 06:00 PM RE: Horizons: Distance to Horizon - StephenG1CMZ - 04-11-2019 09:10 PM RE: Horizons: Distance to Horizon - StephenG1CMZ - 04-15-2019, 09:28 PM

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