Post Reply 
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
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)
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
RE: Horizons: Distance to Horizon - StephenG1CMZ - 04-11-2019 09:10 PM



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