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.
Search
Member List
Calendar
Help
