+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP Prime Software Library (/forum-15.html)
+--- Thread: Horizons: Distance to Horizon (/thread-12769.html)
Horizons: Distance to Horizon - StephenG1CMZ - 04-09-2019 05:55 PM
After recently seeing this thread http://www.hpmuseum.org/forum/thread-12717.html?highlight=Spher
(implementing surface area visible at height)
I have implemented distance to horizon functionality on the HP Prime.
Note: The HP Prime may be capable of better accuracy than some of my examples might suggest...
Some examples probably started out with a different radius.
RE: Horizons: Distance to Horizon - StephenG1CMZ - 04-09-2019 06:00 PM
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;RE: Horizons: Distance to Horizon - StephenG1CMZ - 04-11-2019 09:10 PM
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;RE: Horizons: Distance to Horizon - StephenG1CMZ - 04-15-2019 09:28 PM
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;