GetAngular - Gets angular size of a distant sphere (suitable for celestial use)
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;
Search
Member List
Calendar
Help
