1 #include "sun.h"
2 #include "unsorted.h"
3 #include "defs.h"
4
5 /**
6 * The function Delta_ET has been added to allow calculations on the position of the sun. It provides the difference between UT (approximately the same as UTC) and ET (now referred to as TDT). This function is based on a least squares fit of data from 1950 to 1991 and will need to be updated periodically. Values determined using data from 1950-1991 in the 1990 Astronomical Almanac. See DELTA_ET.WQ1 for details.
7 *
8 * \copyright GPLv2+
9 **/
Delta_ET(double year)10 double Delta_ET(double year)
11 {
12 double delta_et;
13
14 delta_et=26.465+0.747622*(year-1950)+1.886913*sin(2*M_PI*(year-1975)/33);
15
16 return delta_et;
17 }
18
19 /**
20 * Returns angle in radians from argument in degrees.
21 *
22 * \copyright GPLv2+
23 **/
Radians(double arg)24 double Radians(double arg)
25 {
26 /* Returns angle in radians from argument in degrees */
27 return (arg*M_PI/180.0);
28 }
29
30 /**
31 * Returns angle in degrees from argument in radians.
32 *
33 * \copyright GPLv2+
34 **/
Degrees(double arg)35 double Degrees(double arg)
36 {
37 /* Returns angle in degrees from argument in radians */
38 return (arg*180.0/M_PI);
39 }
40
sun_predict(double time,double position[3])41 void sun_predict(double time, double position[3])
42 {
43 double jul_utc = time + JULIAN_TIME_DIFF;
44 double mjd = jul_utc - 2415020.0;
45 double year = 1900 + mjd / 365.25;
46 double T = (mjd + Delta_ET(year) / SECONDS_PER_DAY) / 36525.0;
47 double M = Radians(fmod(358.47583+fmod(35999.04975*T,360.0)-(0.000150+0.0000033*T)*Sqr(T),360.0));
48 double L = Radians(fmod(279.69668+fmod(36000.76892*T,360.0)+0.0003025*Sqr(T),360.0));
49 double e = 0.01675104-(0.0000418+0.000000126*T)*T;
50 double C = Radians((1.919460-(0.004789+0.000014*T)*T)*sin(M)+(0.020094-0.000100*T)*sin(2*M)+0.000293*sin(3*M));
51 double O = Radians(fmod(259.18-1934.142*T,360.0));
52 double Lsa = fmod(L+C-Radians(0.00569-0.00479*sin(O)), 2*M_PI);
53 double nu = fmod(M+C, 2*M_PI);
54 double R = 1.0000002*(1.0-Sqr(e))/(1.0+e*cos(nu));
55 double eps = Radians(23.452294-(0.0130125+(0.00000164-0.000000503*T)*T)*T+0.00256*cos(O));
56 R = ASTRONOMICAL_UNIT_KM*R;
57
58 position[0] = R*cos(Lsa);
59 position[1] = R*sin(Lsa)*cos(eps);
60 position[2] = R*sin(Lsa)*sin(eps);
61 }
62
predict_observe_sun(const predict_observer_t * observer,double time,struct predict_observation * obs)63 void predict_observe_sun(const predict_observer_t *observer, double time, struct predict_observation *obs)
64 {
65
66 // Find sun position
67 double solar_vector[3];
68 sun_predict(time, solar_vector);
69
70 /* Zero vector for initializations */
71 double zero_vector[3] = {0,0,0};
72
73 /* Solar observed azimuth and elevation vector */
74 vector_t solar_set;
75
76 geodetic_t geodetic;
77 geodetic.lat = observer->latitude;
78 geodetic.lon = observer->longitude;
79 geodetic.alt = observer->altitude / 1000.0;
80 geodetic.theta = 0.0;
81
82 double jul_utc = time + JULIAN_TIME_DIFF;
83 Calculate_Obs(jul_utc, solar_vector, zero_vector, &geodetic, &solar_set);
84
85 double sun_azi = solar_set.x;
86 double sun_ele = solar_set.y;
87
88 double sun_range = 1.0+((solar_set.z-ASTRONOMICAL_UNIT_KM)/ASTRONOMICAL_UNIT_KM);
89 double sun_range_rate = 1000.0*solar_set.w;
90
91 obs->time = time;
92 obs->azimuth = sun_azi;
93 obs->elevation = sun_ele;
94 obs->range = sun_range;
95 obs->range_rate = sun_range_rate;
96 }
97
98 /**
99 * Calculate RA and dec for the sun.
100 *
101 * \param time Time
102 * \param ra Right ascension
103 * \param dec Declination
104 * \copyright GPLv2+
105 **/
predict_sun_ra_dec(predict_julian_date_t time,double * ra,double * dec)106 void predict_sun_ra_dec(predict_julian_date_t time, double *ra, double *dec)
107 {
108 //predict absolute position of the sun
109 double solar_vector[3];
110 sun_predict(time, solar_vector);
111
112 //prepare for radec calculation
113 double jul_utc = time + JULIAN_TIME_DIFF;
114 double zero_vector[3] = {0,0,0};
115 vector_t solar_rad;
116
117 //for some reason, RADec requires QTH coordinates, though
118 //the properties to be calculated are observer-independent.
119 //Pick some coordinates, will be correct anyway.
120 geodetic_t geodetic;
121 geodetic.lat = 10;
122 geodetic.lon = 10;
123 geodetic.alt = 10 / 1000.0;
124 geodetic.theta = 0.0;
125
126 //calculate right ascension/declination
127 Calculate_RADec(jul_utc, solar_vector, zero_vector, &geodetic, &solar_rad);
128 *ra = solar_rad.x;
129 *dec = solar_rad.y;
130 }
131
predict_sun_ra(predict_julian_date_t time)132 double predict_sun_ra(predict_julian_date_t time)
133 {
134 double ra, dec;
135 predict_sun_ra_dec(time, &ra, &dec);
136 return ra;
137 }
138
predict_sun_declination(predict_julian_date_t time)139 double predict_sun_declination(predict_julian_date_t time)
140 {
141 double ra, dec;
142 predict_sun_ra_dec(time, &ra, &dec);
143 return dec;
144 }
145
predict_sun_gha(predict_julian_date_t time)146 double predict_sun_gha(predict_julian_date_t time)
147 {
148 //predict absolute position of sun
149 double solar_vector[3];
150 sun_predict(time, solar_vector);
151
152 //convert to lat/lon/alt
153 geodetic_t solar_latlonalt;
154 Calculate_LatLonAlt(time, solar_vector, &solar_latlonalt);
155
156 //return longitude as the GHA
157 double sun_lon = 360.0-Degrees(solar_latlonalt.lon);
158 return sun_lon*M_PI/180.0;
159 }
160