mirror of
https://github.com/SatDump/SatDump
synced 2026-08-13 17:47:30 -04:00
288 lines
8.1 KiB
C
288 lines
8.1 KiB
C
#include "predict.h"
|
|
#include "unsorted.h"
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include "defs.h"
|
|
#include "sun.h"
|
|
|
|
/**
|
|
* This function reduces angles greater than two pi by subtracting two pi
|
|
* from the angle
|
|
*
|
|
* \copyright GPLv2+
|
|
**/
|
|
double FixAngle(double x)
|
|
{
|
|
|
|
while ( x > 2*M_PI )
|
|
x-=2*M_PI;
|
|
|
|
return x;
|
|
}
|
|
|
|
/**
|
|
* This function is used in the FindMoon() function.
|
|
*
|
|
* \copyright GPLv2+
|
|
**/
|
|
double PrimeAngle(double x)
|
|
{
|
|
x=x-360.0*floor(x/360.0);
|
|
return x;
|
|
}
|
|
|
|
/**
|
|
* Output struct for predict_moon().
|
|
**/
|
|
struct moon {
|
|
///Julian day
|
|
double jd;
|
|
///Related to ecliptic longitude
|
|
double lm;
|
|
///Related to ecliptic latitude
|
|
double b;
|
|
///Parallax-related
|
|
double p;
|
|
///Related to siderial time
|
|
double teg;
|
|
///Range approximation?
|
|
double dx;
|
|
};
|
|
|
|
/**
|
|
* Predict absolute, observer-independent properties of the moon.
|
|
*
|
|
* \param time Time
|
|
* \param moon Output struct
|
|
* \copyright GPLv2+
|
|
**/
|
|
void predict_moon(double time, struct moon *moon)
|
|
{
|
|
double jd, t, t2, t3, l1, m, teg, l, b, w1, w2, bt, p, lm, m1, d, ff, om, ss, ex;
|
|
|
|
|
|
jd = time + JULIAN_TIME_DIFF;
|
|
moon->jd = jd;
|
|
|
|
t=(jd-2415020.0)/36525.0;
|
|
t2=t*t;
|
|
t3=t2*t;
|
|
l1=270.434164+481267.8831*t-0.001133*t2+0.0000019*t3;
|
|
m=358.475833+35999.0498*t-0.00015*t2-0.0000033*t3;
|
|
m1=296.104608+477198.8491*t+0.009192*t2+0.0000144*t3;
|
|
d=350.737486+445267.1142*t-0.001436*t2+0.0000019*t3;
|
|
ff=11.250889+483202.0251*t-0.003211*t2-0.0000003*t3;
|
|
om=259.183275-1934.142*t+0.002078*t2+0.0000022*t3;
|
|
om = om * M_PI/180.0;
|
|
|
|
/* Additive terms */
|
|
|
|
l1=l1+0.000233*sin((51.2+20.2*t)*M_PI/180.0);
|
|
ss=0.003964*sin((346.56+132.87*t-0.0091731*t2)*M_PI/180.0);
|
|
l1=l1+ss+0.001964*sin(om);
|
|
m=m-0.001778*sin((51.2+20.2*t)*M_PI/180.0);
|
|
m1=m1+0.000817*sin((51.2+20.2*t)*M_PI/180.0);
|
|
m1=m1+ss+0.002541*sin(om);
|
|
d=d+0.002011*sin((51.2+20.2*t)*M_PI/180.0);
|
|
d=d+ss+0.001964*sin(om);
|
|
ff=ff+ss-0.024691*sin(om);
|
|
ff=ff-0.004328*sin(om+(275.05-2.3*t)*M_PI/180.0);
|
|
ex=1.0-0.002495*t-0.00000752*t2;
|
|
om=om*M_PI/180.0;
|
|
|
|
l1=PrimeAngle(l1);
|
|
m=PrimeAngle(m);
|
|
m1=PrimeAngle(m1);
|
|
d=PrimeAngle(d);
|
|
ff=PrimeAngle(ff);
|
|
om=PrimeAngle(om);
|
|
|
|
m=m*M_PI/180.0;
|
|
m1=m1*M_PI/180.0;
|
|
d=d*M_PI/180.0;
|
|
ff=ff*M_PI/180.0;
|
|
|
|
/* Ecliptic Longitude */
|
|
|
|
l=l1+6.28875*sin(m1)+1.274018*sin(2.0*d-m1)+0.658309*sin(2.0*d);
|
|
l=l+0.213616*sin(2.0*m1)-ex*0.185596*sin(m)-0.114336*sin(2.0*ff);
|
|
l=l+0.058793*sin(2.0*d-2.0*m1)+ex*0.057212*sin(2.0*d-m-m1)+0.05332*sin(2.0*d+m1);
|
|
l=l+ex*0.045874*sin(2.0*d-m)+ex*0.041024*sin(m1-m)-0.034718*sin(d);
|
|
l=l-ex*0.030465*sin(m+m1)+0.015326*sin(2.0*d-2.0*ff)-0.012528*sin(2.0*ff+m1);
|
|
|
|
l=l-0.01098*sin(2.0*ff-m1)+0.010674*sin(4.0*d-m1)+0.010034*sin(3.0*m1);
|
|
l=l+0.008548*sin(4.0*d-2.0*m1)-ex*0.00791*sin(m-m1+2.0*d)-ex*0.006783*sin(2.0*d+m);
|
|
|
|
l=l+0.005162*sin(m1-d)+ex*0.005*sin(m+d)+ex*0.004049*sin(m1-m+2.0*d);
|
|
l=l+0.003996*sin(2.0*m1+2.0*d)+0.003862*sin(4.0*d)+0.003665*sin(2.0*d-3.0*m1);
|
|
|
|
l=l+ex*0.002695*sin(2.0*m1-m)+0.002602*sin(m1-2.0*ff-2.0*d)+ex*0.002396*sin(2.0*d-m-2.0*m1);
|
|
|
|
l=l-0.002349*sin(m1+d)+ex*ex*0.002249*sin(2.0*d-2.0*m)-ex*0.002125*sin(2.0*m1+m);
|
|
|
|
l=l-ex*ex*0.002079*sin(2.0*m)+ex*ex*0.002059*sin(2.0*d-m1-2.0*m)-0.001773*sin(m1+2.0*d-2.0*ff);
|
|
|
|
l=l+ex*0.00122*sin(4.0*d-m-m1)-0.00111*sin(2.0*m1+2.0*ff)+0.000892*sin(m1-3.0*d);
|
|
|
|
l=l-ex*0.000811*sin(m+m1+2.0*d)+ex*0.000761*sin(4.0*d-m-2.0*m1)+ex*ex*.000717*sin(m1-2.0*m);
|
|
|
|
l=l+ex*ex*0.000704*sin(m1-2.0*m-2.0*d)+ex*0.000693*sin(m-2.0*m1+2.0*d)+ex*0.000598*sin(2.0*d-m-2.0*ff)+0.00055*sin(m1+4.0*d);
|
|
|
|
l=l+0.000538*sin(4.0*m1)+ex*0.000521*sin(4.0*d-m)+0.000486*sin(2.0*m1-d);
|
|
|
|
l=l-0.001595*sin(2.0*ff+2.0*d);
|
|
|
|
/* Ecliptic latitude */
|
|
|
|
b=5.128189*sin(ff)+0.280606*sin(m1+ff)+0.277693*sin(m1-ff)+0.173238*sin(2.0*d-ff);
|
|
b=b+0.055413*sin(2.0*d+ff-m1)+0.046272*sin(2.0*d-ff-m1)+0.032573*sin(2.0*d+ff);
|
|
|
|
b=b+0.017198*sin(2.0*m1+ff)+9.266999e-03*sin(2.0*d+m1-ff)+0.008823*sin(2.0*m1-ff);
|
|
b=b+ex*0.008247*sin(2.0*d-m-ff)+0.004323*sin(2.0*d-ff-2.0*m1)+0.0042*sin(2.0*d+ff+m1);
|
|
|
|
b=b+ex*0.003372*sin(ff-m-2.0*d)+ex*0.002472*sin(2.0*d+ff-m-m1)+ex*0.002222*sin(2.0*d+ff-m);
|
|
|
|
b=b+0.002072*sin(2.0*d-ff-m-m1)+ex*0.001877*sin(ff-m+m1)+0.001828*sin(4.0*d-ff-m1);
|
|
|
|
b=b-ex*0.001803*sin(ff+m)-0.00175*sin(3.0*ff)+ex*0.00157*sin(m1-m-ff)-0.001487*sin(ff+d)-ex*0.001481*sin(ff+m+m1)+ex*0.001417*sin(ff-m-m1)+ex*0.00135*sin(ff-m)+0.00133*sin(ff-d);
|
|
|
|
b=b+0.001106*sin(ff+3.0*m1)+0.00102*sin(4.0*d-ff)+0.000833*sin(ff+4.0*d-m1);
|
|
|
|
b=b+0.000781*sin(m1-3.0*ff)+0.00067*sin(ff+4.0*d-2.0*m1)+0.000606*sin(2.0*d-3.0*ff);
|
|
|
|
b=b+0.000597*sin(2.0*d+2.0*m1-ff)+ex*0.000492*sin(2.0*d+m1-m-ff)+0.00045*sin(2.0*m1-ff-2.0*d);
|
|
|
|
b=b+0.000439*sin(3.0*m1-ff)+0.000423*sin(ff+2.0*d+2.0*m1)+0.000422*sin(2.0*d-ff-3.0*m1);
|
|
|
|
b=b-ex*0.000367*sin(m+ff+2.0*d-m1)-ex*0.000353*sin(m+ff+2.0*d)+0.000331*sin(ff+4.0*d);
|
|
|
|
b=b+ex*0.000317*sin(2.0*d+ff-m+m1)+ex*ex*0.000306*sin(2.0*d-2.0*m-ff)-0.000283*sin(m1+3.0*ff);
|
|
|
|
w1=0.0004664*cos(om*M_PI/180.0);
|
|
w2=0.0000754*cos((om+275.05-2.3*t)*M_PI/180.0);
|
|
bt=b*(1.0-w1-w2);
|
|
|
|
/* Parallax calculations */
|
|
|
|
p=0.950724+0.051818*cos(m1)+0.009531*cos(2.0*d-m1)+0.007843*cos(2.0*d)+0.002824*cos(2.0*m1)+0.000857*cos(2.0*d+m1)+ex*0.000533*cos(2.0*d-m)+ex*0.000401*cos(2.0*d-m-m1);
|
|
|
|
p=p+0.000173*cos(3.0*m1)+0.000167*cos(4.0*d-m1)-ex*0.000111*cos(m)+0.000103*cos(4.0*d-2.0*m1)-0.000084*cos(2.0*m1-2.0*d)-ex*0.000083*cos(2.0*d+m)+0.000079*cos(2.0*d+2.0*m1);
|
|
|
|
p=p+0.000072*cos(4.0*d)+ex*0.000064*cos(2.0*d-m+m1)-ex*0.000063*cos(2.0*d+m-m1);
|
|
|
|
p=p+ex*0.000041*cos(m+d)+ex*0.000035*cos(2.0*m1-m)-0.000033*cos(3.0*m1-2.0*d);
|
|
|
|
p=p-0.00003*cos(m1+d)-0.000029*cos(2.0*ff-2.0*d)-ex*0.000029*cos(2.0*m1+m);
|
|
|
|
p=p+ex*ex*0.000026*cos(2.0*d-2.0*m)-0.000023*cos(2.0*ff-2.0*d+m1)+ex*0.000019*cos(4.0*d-m-m1);
|
|
|
|
b=bt*M_PI/180.0;
|
|
lm=l*M_PI/180.0;
|
|
|
|
/* Find siderial time in radians */
|
|
|
|
t=(jd-2451545.0)/36525.0;
|
|
teg=280.46061837+360.98564736629*(jd-2451545.0)+(0.000387933*t-t*t/38710000.0)*t;
|
|
|
|
while (teg>360.0)
|
|
teg-=360.0;
|
|
|
|
//output
|
|
moon->b = b;
|
|
moon->lm = lm;
|
|
moon->p = p;
|
|
moon->dx = 3.0/(M_PI*p);
|
|
moon->teg = teg;
|
|
}
|
|
|
|
/**
|
|
* Calculate RA and dec for the moon.
|
|
*
|
|
* \param time Time
|
|
* \param ra Right ascension
|
|
* \param dec Declination
|
|
* \copyright GPLv2+
|
|
**/
|
|
void predict_moon_ra_dec(predict_julian_date_t time, double *ra, double *dec)
|
|
{
|
|
struct moon moon;
|
|
predict_moon(time, &moon);
|
|
|
|
/* Semi-diameter calculation */
|
|
/* sem=10800.0*asin(0.272488*p*M_PI/180.0)/pi; */
|
|
/* Convert ecliptic coordinates to equatorial coordinates */
|
|
|
|
double z=(moon.jd-2415020.5)/365.2422;
|
|
double ob=23.452294-(0.46845*z+5.9e-07*z*z)/3600.0;
|
|
ob=ob*M_PI/180.0;
|
|
*dec=asin(sin(moon.b)*cos(ob)+cos(moon.b)*sin(ob)*sin(moon.lm));
|
|
*ra=acos(cos(moon.b)*cos(moon.lm)/cos(*dec));
|
|
|
|
if (moon.lm > M_PI)
|
|
*ra = 2*M_PI - *ra;
|
|
}
|
|
|
|
void predict_observe_moon(const predict_observer_t *observer, double time, struct predict_observation *obs)
|
|
{
|
|
struct moon moon;
|
|
predict_moon(time, &moon);
|
|
|
|
double ra, dec;
|
|
predict_moon_ra_dec(time, &ra, &dec);
|
|
|
|
double n = observer->latitude; /* North latitude of tracking station */
|
|
double e = observer->longitude; /* East longitude of tracking station */
|
|
|
|
|
|
double th = FixAngle(moon.teg*M_PI/180.0 + e);
|
|
double h=th-ra;
|
|
|
|
double az=atan2(sin(h),cos(h)*sin(n)-tan(dec)*cos(n))+M_PI;
|
|
double el=asin(sin(n)*sin(dec)+cos(n)*cos(dec)*cos(h));
|
|
|
|
/* Radial velocity approximation. This code was derived
|
|
from "Amateur Radio Software", by John Morris, GM4ANB,
|
|
published by the RSGB in 1985. */
|
|
|
|
double mm=FixAngle(1.319238+time*0.228027135); /* mean moon position */
|
|
double t2=0.10976;
|
|
double t1=mm+t2*sin(mm);
|
|
double dv=0.01255*moon.dx*moon.dx*sin(t1)*(1.0+t2*cos(mm));
|
|
dv=dv*4449.0;
|
|
t1=6378.0;
|
|
t2=384401.0;
|
|
double t3=t1*t2*(cos(dec)*cos(n)*sin(h));
|
|
t3=t3/sqrt(t2*t2-t2*t1*sin(el));
|
|
|
|
double moon_dv=dv+t3*0.0753125;
|
|
|
|
obs->time = time;
|
|
obs->azimuth = az;
|
|
obs->elevation = el;
|
|
obs->range = moon.dx;
|
|
obs->range_rate = moon_dv;
|
|
}
|
|
|
|
double predict_moon_ra(predict_julian_date_t time)
|
|
{
|
|
double ra, dec;
|
|
predict_moon_ra_dec(time, &ra, &dec);
|
|
return ra;
|
|
}
|
|
|
|
double predict_moon_declination(predict_julian_date_t time)
|
|
{
|
|
double ra, dec;
|
|
predict_moon_ra_dec(time, &ra, &dec);
|
|
return dec;
|
|
}
|
|
|
|
double predict_moon_gha(predict_julian_date_t time)
|
|
{
|
|
struct moon moon;
|
|
predict_moon(time, &moon);
|
|
double moon_gha=moon.teg-predict_moon_ra(time)*180.0/M_PI;
|
|
|
|
if (moon_gha<0.0) moon_gha+=360;
|
|
return moon_gha*M_PI/180.0;
|
|
}
|