Horizons: Distance to Horizon
04-09-2019, 06:00 PM (This post was last modified: 04-10-2019 08:48 AM by StephenG1CMZ.)
Post: #2
 StephenG1CMZ Senior Member Posts: 918 Joined: May 2015
RE: Horizons: Distance to Horizon
Version 0.1 implements simple geometric solutions for spheres, taking no account of refraction.

It includes functions for distance to horizon at a given height, and inverse functions.

It also demonstrates some unexpected cas solve behaviour:
The annotated example Print, if included, often seems OK on the 1st run (returning undef), but on later runs all solves return undef, suggestive of an initialisation issue. Earlier versions of the program sometimes generated a large real instead of undef on the 1st run (and on one occasion a large negative, despite the ">=0").
Update: Changing HT to h in the cas solve seems to fix the 2nd run undef problem.
Code:
  LOCAL CRID:="Horizons V0.1 © 2019 StephenG1CMZ";  //This version has simple geometric horizons  //No refraction    LOCAL HT;//FOR SOLVE    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", VALUE RETURNED WHEN EG HT<0;  //Definitions  //DD   Direct Distance       (eye-horizon)  //HT   Height above Sphere   (<0=NA)  //SLen Curved Surface Length (paw-horizon)  //Distances and heights in same units    SpSphere(RADIUS)  //Return Sphere CAV  BEGIN   LOCAL SpC,SpA,SpV;   SpC:=2*π*RADIUS;   SpA:=4*π*RADIUS^2;   SpV:=((4/3)*π*RADIUS^3);   RETURN {SpC,SpA,SpV};   END;  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      RETURN CAS("solve((HT*(2*RADIUS+HT))-DD^2=0,HT,HT≥0)");  END;  EXPORT GetHeightSD(SurfLen,RADIUS)  //INV   //RETURN HT TO SEE SLEN (SD) AWAY  // RETURN HT (OR 0 IF SD IMPOSSIBLE TBD)  BEGIN      RETURN (RADIUS/COS(SurfLen/RADIUS))-RADIUS;   //RETURN {(RADIUS/COS(SD/RADIUS))-RADIUS,CAS("solve((RADIUS/COS(SD/RADIUS))-RADIUS-HT=0,HT)")};  END;  EXPORT GetHeightSAF(SAF,RADIUS)  //INV  //Get Height for Surface Area (SA in units EG m^2,NOT %)//correction saf as a fraction  BEGIN    HT:=CAS("solve(HT/(2*(HT+RADIUS))-SAF=0,HT,HT≥0)");    RETURN 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;    GetSurfaceVisibleHtPC(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/1000)}));//km   //PRINT("Verify: HT SD "+GetHeightSD(GetSurfLenHT(HT,RADIUS),RADIUS));   PRINT("Verify: HT DD "+GetHeightDD(GetDLen     (HT,RADIUS),RADIUS));   PRINT("Verify: HT SAF"+GetHeightSAF(FP,RADIUS));   PRINT(" ");  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. C:A:V:");   PRINT(SpSphere(RADIUS/1000));//km   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("//40.5 AU Voyager PALE BLUE DOT");   EX(40.5*AU2KM,RADIUS);   //PRINT("//145 AU Voyager 2019");   //EX(145*AU2KM,RADIUS);   PRINT("//as Exoplanet of Proxima Centauri");   EX(268395.131*AU2KM,RADIUS);   //This last verifies, except solve doesnt      //CAUTION:   //INCLUDING THIS PRINT HAS THIS EFFECT:   //1ST RUN AS EXPECTED (THIS VALUE undef)   //2ND RUN ALL solve RETURN undef   //PRINT(GetHeightSAF(0.50,RADIUS));//sometimes undef   //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;

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)