Post Reply 
Horizons: Distance to Horizon
04-15-2019, 09:28 PM (This post was last modified: 04-15-2019 09:50 PM by StephenG1CMZ.)
Post: #4
RE: Horizons: Distance to Horizon
New in Version 0.3:
GetAngular - Gets angular size of a distant sphere (suitable for celestial use)
GetDip2Horizon (from eye height)

Improves: Examples, PrintReport

Code:


 LOCAL CRID:="Horizons V0.3 © 2019 StephenG1CMZ";
 //This version has simple geometric horizons
 //No refraction
 
 LOCAL HT;//FOR SOLVE //Global h used
 
 LOCAL AU2KM:=149597871;//GOOGLE
 LOCAL AU2M:=1000*AU2KM;
 LOCAL KR2D:=180/π;
 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 RSUN:=109*RADIUS;

 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 OR √(<0) (HOPEFULLY ONE OF THE ABOVE)
  //AVOID ASIN>1 (ie RADIUS>DistC)

 //Definitions //Distances and heights in same units
 //DD    Direct Distance       (eye-horizon)
 //DistC DistanceToC=HT+RADIUS (eye to centre of sphere)
 //HT    Height above Sphere   (eye to surface. <0=NA) 
 //RADIUS                      (of Sphere)
 //SLen  Curved Surface Length (paw-horizon)

 //AU    Astronomical Unit
 //mas   milliarcseconds
 
 GetSphereSA(RADIUS)
 //RETURN SPHERE SURFACE AREA
 BEGIN
  RETURN 4*π*RADIUS^2;
 END;

 GetAngularRadians(DistC,RADIUS)
 //The observed angular diameter at the eye (ignoring physiology)
 ///This DistC ASIN implement is suitable for celestial use (observing 3D spheres)
 ///(An alternate using HT ATAN is known for observing 2D disks)
 BEGIN
  LOCAL DIAMETER:=2*RADIUS;
  //LOCAL AS;
  IF DistC==0 THEN
   RETURN NA;//IF INSTEAD HT=0 WERE TRAPPED return π
  END;
 
  IF RADIUS>DistC THEN //TOO WIDE TO SEE??
   //IF NOT TRAPPED OUTPUT JUST STOPS!
   MSGBOX("HORIZONS:GetAngular(RADIUS>DX) => ASIN>1 ERR");
   MSGBOX({RADIUS,DistC}); 
   RETURN NA; //π/2; //BEST VALUE TO RETURN TBC
  END; 
  RETURN 2*ASIN(DIAMETER/(2*DistC)); 
 END;
 EXPORT GetAngularDegrees(DistC,RADIUS)
 BEGIN
  RETURN GetAngularRadians(DistC,RADIUS)*KR2D;
 END;

 GetDipRadiansDD(DD,RADIUS)
 //Dip to horizon
 BEGIN
  IF RADIUS==0 THEN
   RETURN NA;
  END;
  RETURN ATAN(DD/RADIUS);
 END;
 EXPORT GetDipDegreesDD(DD,RADIUS)
 BEGIN
  RETURN GetDipRadiansDD(DD,RADIUS)*KR2D;
 END;

 GetDipRadiansHT(HT,RADIUS)
 //Dip to horizon
 BEGIN
  IF RADIUS+HT==0 THEN
   RETURN NA;
  END;
  RETURN ACOS(RADIUS/(RADIUS+HT));
  //OR RETURN ATAN(DD/RADIUS);
 END;
 EXPORT GetDipDegreesHT(HT,RADIUS)
 BEGIN
  RETURN GetDipRadiansHT(HT,RADIUS)*KR2D;
 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 % INACCURATE
  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
 //UNITS IS ONLY FOR READABILITY
 //Hint: See EX 
 BEGIN
  LOCAL PLACES:=3;//ROUNDING
  LOCAL UNITS:=" ";
  LOCAL DX:=GetDLen(HT,RADIUS);
  LOCAL FP:=GetSurfAreaHTFP(HT,RADIUS);
 
  PRINT("At Height "+HT+UNITS+" Radius "+RADIUS+UNITS);
  PRINT("Eye-horizon : "+ROUND(DX,PLACES)+UNITS);
  PRINT("Paw-horizon : "+ROUND(GetSurfLenDD(DX,RADIUS),PLACES)+UNITS);
  PRINT("Dip-horizon : "+ROUND(GetDipDegreesHT(HT,RADIUS),PLACES)+" °");
  PRINT("Surface Area: "+ROUND(FP*GetSphereSA(RADIUS/1000),PLACES));//m TO km
 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("Dip: "+GetDipDegreesHT(HT,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 elsewhere
  //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 about 408 km ( ~2373 km)"); 
  
  EX(427200,RADIUS);
  PRINT("//51ᴇ3 mas,140 °");
  PRINT(60*60*1000*GetAngularDegrees(427200,108.5/2)+" mas"); //ISS fom Earth
  //The next yields an ASIN error unless trapped.
  //PRINT(           GetAngularDegrees(427200,       RADIUS )+" °");   //Earth from ISS - Wrong Parameters?
  // OK
  PRINT(           GetAngularDegrees(427200+RADIUS,RADIUS )+" °");   //Earth from ISS

  PRINT("//GPS SATELLITES (38%)");
  EX(20000ᴇ3,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("//Sun 31-32");
  PRINT(60*GetAngularDegrees(AU2M,RSUN));

  PRINT("//40.5 AU Voyager PALE BLUE DOT");
  EX(40.5*AU2M,RADIUS);
  //PRINT("//145 AU Voyager 2019");
  //EX(145*AU2M,RADIUS);

  PRINT("//Viewed as Exoplanet from Proxima Centauri");
  EX(268395.131*AU2M,RADIUS); 
  //This last verifies, except solve doesnt
 
  PRINT("Proxima Centauri from here:");
  PRINT("//angular: (1.02 mas)"); 
  PRINT(1000*60*60*GetAngularDegrees(268395.131*AU2M,0.1542*RSUN)+" mas");//Proxima DX and radius
  //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 DD;
  LOCAL HT:=2;
  PRINT();
  PRINT(CRID);
  EXAMPLES();

  //PRINT("TEVAL DStraight "+TEVAL(DD:=HZDistanceD(HT,RADIUS)));// 58 CALC D 
  //PRINT("TEVAL SurfaceD "+TEVAL(HZDistanceSD(DD,RADIUS)));    // 54 S ADD TO ABOVE 
  //PRINT("TEVAL SurfaceH "+TEVAL(HZDistanceSH(HT,RADIUS)));    //109 VS S 
 END;

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-15-2019 09:28 PM



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