1#ifndef _PREDICT_H_
2#define _PREDICT_H_
3
4#ifdef __cplusplus
5extern "C" {
6#endif
7
8#include <time.h>
9#include <stdbool.h>
10
11#define PREDICT_VERSION_MAJOR @PROJECT_VERSION_MAJOR@
12#define PREDICT_VERSION_MINOR @PROJECT_VERSION_MINOR@
13#define PREDICT_VERSION_PATCH @PROJECT_VERSION_PATCH@
14#define PREDICT_VERSION (PREDICT_VERSION_MAJOR * 100 * 100 + PREDICT_VERSION_MINOR * 100 + PREDICT_VERSION_PATCH)
15#define PREDICT_VERSION_STRING "@PROJECT_VERSION@"
16
17/**
18 * Get the major version number of the library
19 *
20 * \return Major version number
21 */
22int predict_version_major();
23
24/**
25 * Get the minor version number of the library
26 *
27 * \return Minor version number
28 */
29int predict_version_minor();
30
31/**
32 * Get the patch version number of the library
33 *
34 * \return Patch version number
35 */
36int predict_version_patch();
37
38/**
39 * Get the version number of the library
40 *
41 * \return 2 digit major, 2 digit minor, 2 digit patch decimal version number
42 */
43int predict_version();
44
45/**
46 * Get the version number string of the library
47 *
48 * \return Version number string ("major.minor.patch")
49 */
50char *predict_version_string();
51
52/**
53 * The representation of time used by libpredict: The number of days since 31Dec79 00:00:00 UTC.
54 **/
55typedef double predict_julian_date_t;
56
57/**
58 * Convert time_t in UTC to Julian date in UTC.
59 *
60 * \param time Time in UTC
61 * \return Julian day in UTC
62 **/
63predict_julian_date_t predict_to_julian(time_t time);
64
65/**
66 * Convert Julian date in UTC back to a time_t in UTC.
67 *
68 * \param date Julian date in UTC
69 * \return Time in UTC
70 **/
71time_t predict_from_julian(predict_julian_date_t date);
72
73/**
74 * Simplified perturbation models used in modeling the satellite orbits.
75 **/
76enum predict_ephemeris {
77  EPHEMERIS_SGP4 = 0,
78  EPHEMERIS_SDP4 = 1,
79  EPHEMERIS_SGP8 = 2,
80  EPHEMERIS_SDP8 = 3
81};
82
83/**
84 * Container for processed TLE data from TLE strings.
85 **/
86typedef struct {
87	///Satellite number (line 1, field 2)
88	int satellite_number;
89	///Element number (line 1, field 13)
90	long element_number;
91	///International designator (line 1, fields 4, 5, 6)
92	char designator[10];
93	///Epoch year (last two digits) (line 1, field 7)
94	int epoch_year;
95	///Epoch day (day of year and fractional portion of day, line 1, field 8)
96	double epoch_day;
97	///Inclination (line 2, field 3)
98	double inclination;
99	///Right Ascension of the Ascending Node [Degrees] (line 2, field 4)
100	double right_ascension;
101	///Eccentricity (line 2, field 5)
102	double eccentricity;
103	///Argument of Perigee [Degrees] (line 2, field 6)
104	double argument_of_perigee;
105	///Mean Anomaly [Degrees] (line 2, field 7)
106	double mean_anomaly;
107	///Mean Motion [Revs per day] (line 2, field 8)
108	double mean_motion;
109	///First Time Derivative of the Mean Motion divided by two (line 1, field 9)
110	double derivative_mean_motion;
111	///Second Time Derivative of Mean Motion divided by six (line 1, field 10)
112	double second_derivative_mean_motion;
113	///BSTAR drag term (decimal point assumed, line 1, field 11)
114	double bstar_drag_term;
115	///Number of revolutions around Earth at epoch (line 2, field 9)
116	int revolutions_at_epoch;
117
118	///Which perturbation model to use
119	enum predict_ephemeris ephemeris;
120	///Ephemeris data structure pointer
121	void *ephemeris_data;
122} predict_orbital_elements_t;
123
124/**
125 * Create predict_orbital_elements_t from TLE strings.
126 *
127 * \param tle_line_1 First line of NORAD two-line element set string
128 * \param tle_line_2 Second line of NORAD two-line element set string
129 * \return Processed TLE parameters
130 * \copyright GPLv2+
131 **/
132predict_orbital_elements_t* predict_parse_tle(const char *tle_line_1, const char *tle_line_2);
133
134/**
135 * Free memory allocated in orbital elements structure.
136 * \param orbital_elements Orbit to free
137 **/
138void predict_destroy_orbital_elements(predict_orbital_elements_t *orbital_elements);
139
140/**
141 * Predicted orbital values for satellite at a given time.
142 **/
143struct predict_position {
144	///Timestamp for last call to orbit_predict
145	predict_julian_date_t time;
146
147	///Whether the orbit has decayed
148	bool decayed;
149
150	///ECI position in km
151	double position[3];
152	///ECI velocity in km/s
153	double velocity[3];
154
155	///Latitude in radians, northing/easting
156	double latitude;
157	///Longitude in radians, northing/easting
158	double longitude;
159	///Altitude in km
160	double altitude;
161	///Footprint diameter in km
162	double footprint;
163	///Whether satellite is eclipsed by the earth
164	int eclipsed;
165	///Eclipse depth
166	double eclipse_depth;
167	///Orbital phase (mean anomaly)
168	double phase;
169	///The current number of revolutions around Earth
170	long revolutions;
171
172	///Current inclination (from xinck within sgp4/sdp4)
173	double inclination;
174	///Current right ascension of the ascending node (from xnodek within sgp4/sdp4)
175	double right_ascension;
176	///Current argument of perigee (from omgadf within sgp4/sdp4)
177	double argument_of_perigee;
178};
179
180/**
181 * Main prediction function. Predict satellite orbit at given time.
182 * \param orbital_elements Orbital elements
183 * \param x Predicted orbit
184 * \param time Julian day in UTC
185 * \return 0 if everything went fine
186 * \copyright GPLv2+
187 **/
188int predict_orbit(const predict_orbital_elements_t *orbital_elements, struct predict_position *x, predict_julian_date_t time);
189
190/**
191 * Find whether an orbit is geosynchronous.
192 *
193 * This function uses the definition of geosynchronous orbits found in
194 * "Classification of geosynchronous objects", Issue 17, 28 March 2015, from the
195 * European Space Agency:
196 *
197 * - Eccentricity smaller than 0.2
198 * - Mean motion between 0.9 and 1.1
199 * - Inclination lower than 70 degrees
200 *
201 * The function is mainly used internally for avoiding long iteration loops in
202 * functions like predict_at_max_elevation() and predict_next_aos(). The wider
203 * definition of a geosynchronous orbits is appropriate here. The definition of
204 * a geostationary satellite would be stricter, but is not considered here.
205 *
206 * \param orbital_elements Orbital elements
207 * \return true if orbit is geosynchronous, false otherwise
208 **/
209bool predict_is_geosynchronous(const predict_orbital_elements_t *orbital_elements);
210
211/**
212 * Get apogee of satellite orbit.
213 *
214 * \param x Orbital elements
215 * \return Apogee of orbit
216 * \copyright GPLv2+
217 **/
218double predict_apogee(const predict_orbital_elements_t *x);
219
220/**
221 * Get perigee of satellite orbit.
222 *
223 * \param x Orbital elements
224 * \return Perigee of orbit
225 * \copyright GPLv2+
226 **/
227double predict_perigee(const predict_orbital_elements_t *x);
228
229/**
230 * Find whether an AOS can ever happen on the given latitude.
231 *
232 * \param x Orbital elements
233 * \param latitude Latitude of ground station in radians
234 * \return true if AOS can happen, otherwise false
235 * \copyright GPLv2+
236 **/
237bool predict_aos_happens(const predict_orbital_elements_t *x, double latitude);
238
239/**
240 * Observation point/ground station (QTH).
241 **/
242typedef struct {
243	///Observatory name
244	char name[128];
245	///Latitude (WGS84, radians)
246	double latitude;
247	///Longitude (WGS84, radians)
248	double longitude;
249	///Altitude (WGS84, meters)
250	double altitude;
251} predict_observer_t;
252
253/**
254 * Data relevant for a relative observation of an orbit or similar with respect to an observation point.
255 **/
256struct predict_observation {
257	///UTC time
258	predict_julian_date_t time;
259	///Azimuth angle (rad)
260	double azimuth;
261	///Azimuth angle rate (rad/s)
262	double azimuth_rate;
263	///Elevation angle (rad)
264	double elevation;
265	///Elevation angle rate (rad/s)
266	double elevation_rate;
267	///Range (km)
268	double range;
269	///Range vector
270	double range_x, range_y, range_z;
271	///Range velocity (km/s)
272	double range_rate;
273	///Visibility status, whether satellite can be seen by optical means.
274	///The satellite is defined to be visible if:
275	//	- The satellite is in sunlight
276	//	- The satellite is above the horizon
277	//	- The sky is dark enough (sun elevation is below a fixed threshold)
278	bool visible;
279};
280
281/**
282 * Create observation point (QTH).
283 *
284 * \param name Name of observation point
285 * \param lat Latitude in radians (easting/northing)
286 * \param lon Longitude in radians (easting/northing)
287 * \param alt Altitude in meters
288 * \return Allocated observation point
289 **/
290predict_observer_t *predict_create_observer(const char *name, double lat, double lon, double alt);
291
292/**
293 * Free observer.
294 *
295 * \param obs Observer to be freed.
296 **/
297void predict_destroy_observer(predict_observer_t *obs);
298
299/**
300 * Find relative position of satellite with respect to an observer. Calculates range, azimuth, elevation and relative velocity.
301 *
302 * \param observer Point of observation
303 * \param orbit Satellite orbit
304 * \param obs Return of object for position of the satellite relative to the observer.
305 * \copyright GPLv2+
306 **/
307void predict_observe_orbit(const predict_observer_t *observer, const struct predict_position *orbit, struct predict_observation *obs);
308
309/**
310 * Estimate relative position of the moon.
311 *
312 * \param observer Point of observation
313 * \param time Time of observation
314 * \param obs Return object for position of the moon relative to the observer
315 * \copyright GPLv2+
316 **/
317void predict_observe_moon(const predict_observer_t *observer, predict_julian_date_t time, struct predict_observation *obs);
318
319/**
320 * Calculate right ascension of the moon.
321 *
322 * \param time Time
323 * \return RA in radians
324 **/
325double predict_moon_ra(predict_julian_date_t time);
326
327/**
328 * Calculate declination of the moon.
329 *
330 * \param time Time
331 * \return Declination in radians
332 **/
333double predict_moon_declination(predict_julian_date_t time);
334
335/**
336 * Calculate the greenwich hour angle (longitude) of the moon.
337 *
338 * \param time Time
339 * \return GHA in radians
340 * \copyright GPLv2+
341 **/
342double predict_moon_gha(predict_julian_date_t time);
343
344/**
345 * Estimate relative position of the sun.
346 *
347 * \param observer Point of observation
348 * \param time Time of observation
349 * \param obs Return object for position of the sun relative to the observer
350 * \copyright GPLv2+
351 **/
352void predict_observe_sun(const predict_observer_t *observer, predict_julian_date_t time, struct predict_observation *obs);
353
354/**
355 * Calculate right ascension of the sun.
356 *
357 * \param observer Point of observation
358 * \param time Time of observation
359 * \return RA in radians
360 **/
361double predict_sun_ra(predict_julian_date_t time);
362
363/**
364 * Calculate declination of the sun.
365 *
366 * \param observer Point of observation
367 * \param time Time of observation
368 * \return Declination in radians
369 **/
370double predict_sun_declination(predict_julian_date_t time);
371
372/**
373 * Calculate the greenwich hour angle (longitude) of the sun.
374 *
375 * \param time Time
376 * \return GHA in radians
377 * \copyright GPLv2+
378 **/
379double predict_sun_gha(predict_julian_date_t time);
380
381/**
382 * Find next acquisition of signal (AOS) of satellite (when the satellite rises above the horizon). Ignores previous AOS of current pass if the satellite is in range at the start time.
383 *
384 * \param observer Point of observation
385 * \param orbital_elements Orbital elements
386 * \param start_time Start time for AOS search
387 * \return Observation of the AOS
388 * \copyright GPLv2+
389 **/
390struct predict_observation predict_next_aos(const predict_observer_t *observer, const predict_orbital_elements_t *orbital_elements, predict_julian_date_t start_time);
391
392/**
393 * Find next loss of signal (LOS) of satellite (when the satellite goes below the horizon). Finds LOS of the current pass if the satellite currently is in range, finds LOS of next pass if not.
394 *
395 * \param observer Point of observation
396 * \param orbital_elements Orbital elements
397 * \param start_time Start time for LOS search
398 * \return Observation of the LOS
399 * \copyright GPLv2+
400 **/
401struct predict_observation predict_next_los(const predict_observer_t *observer, const predict_orbital_elements_t *orbital_elements, predict_julian_date_t start_time);
402
403/**
404 * Find maximum elevation of next or current pass.
405 *
406 * \param observer Ground station
407 * \param orbital_elements Orbital elements of satellite
408 * \param start_time Search time. If elevation is negative, max elevation is sought from the start_time and on. If elevation is positive, max elevation is searched for within the current pass
409 * \return Observed properties at maximum elevation
410 **/
411struct predict_observation predict_at_max_elevation(const predict_observer_t *observer, const predict_orbital_elements_t *orbital_elements, predict_julian_date_t start_time);
412
413/**
414 * Calculate doppler shift of a given downlink frequency with respect to an observer.
415 *
416 * \param observation Observation of a satellite orbit
417 * \param downlink_frequency Downlink frequency of the satellite
418 * \return The frequency difference from the original frequency
419 * \copyright GPLv2+
420 **/
421double predict_doppler_shift(const struct predict_observation *observation, double downlink_frequency);
422
423/**
424 * Calculate squint angle for satellite, i.e. angle between the satellite antenna and the QTH antenna.
425 *
426 * \param observer Point of observation
427 * \param orbit Current state of satellite orbit
428 * \param alon Attitude longitude in radians (describes orientation of the satellite at apogee)
429 * \param alat Attidue latitude in radians (see above)
430 * \return Squint angle in radians. Will output nan if the satellite is not an SDP4 satellite
431 * \copyright GPLv2+
432 **/
433double predict_squint_angle(const predict_observer_t *observer, const struct predict_position *orbit, double alon, double alat);
434
435/*!
436 * \brief Calculate refraction angle.
437 *
438 * This function assumes atmospheric pressure of 101.0kPa and temperature 10deg celsius.
439 *
440 * \param el True elevation angle (rad).
441 *
442 * \return Refraction angle (rad).
443 */
444double predict_refraction(double el);
445
446/*!
447 * \brief Calculate refraction angle.
448 *
449 * Corrects for different atmospheric pressure and temperature.
450 *
451 * \param el True elevation angle in rads.
452 * \param pressure Atmospheric pressure in kPa.
453 * \param temp Temperature in deg celsius.
454 *
455 * \return Refraction angle (rad).
456 */
457double predict_refraction_ext(double el, double pressure, double temp);
458
459/*!
460 * \brief Calculate refraction angle from apparent elevation.
461 *
462 * This function assumes atmospheric pressure of 101.0kPa and temperature 10deg celsius.
463 *
464 * \param apparent_el Apparent elevation angle (rad).
465 *
466 * \return Refraction angle (rad).
467 */
468double predict_refraction_from_apparent(double apparent_el);
469
470/*!
471 * \brief Calculate refraction angle from apparent elevation.
472 *
473 * Corrects for different atmospheric pressure and temperature.
474 *
475 * \param apparent_el Apparent elevation angle (rad).
476 * \param pressure Atmospheric pressure in kPa.
477 * \param temp Temperature in deg celsius.
478 *
479 * \return Refraction angle (rad).
480 */
481double predict_refraction_from_apparent_ext(double apparent_el, double pressure, double temp);
482
483/*!
484 * \brief Calculate refraction rate of change.
485 *
486 * \param el True elevation angle (rad).
487 * \param el_rate Rate of change of true elevation angle (rad/s).
488 *
489 * \return Refraction rate of change (rad/s).
490 */
491double predict_refraction_rate(double el, double el_rate);
492
493/*!
494 * \brief Calculate refraction rate of change.
495 *
496 * Corrects for different atmospheric pressure and temerature.
497 *
498 * \param el True elevation angle (rad).
499 * \param el_rate Rate of change of true elevation angle (rad/s).
500 * \param pressure Atmospheric pressure in kPa.
501 * \param temp Temperature in deg celsius.
502 *
503 * \return Apparent elevation (rad).
504 */
505double predict_refraction_rate_ext(double el, double el_rate, double pressure, double temp);
506
507/*!
508 * \brief Calculate apparent elevation from true elevation.
509 *
510 * \param el True elevation angle (rad).
511 *
512 * \return Apparent elevation (rad).
513 */
514double predict_apparent_elevation(double el);
515
516/*!
517 * \brief Calculate apparent elevation from true elevation.
518 *
519 * Corrects for different atmospheric pressures and temperatures.
520 *
521 * \param el True elevation angle (rad).
522 * \param pressure Atmospheric pressure (kPa).
523 * \param temp Temperature (deg C).
524 *
525 * \return Apparent elevation (rad).
526 */
527double predict_apparent_elevation_ext(double el, double pressure, double temp);
528
529/*!
530 * \brief Calculate apparent elevation rate.
531 *
532 * \param el True elevation angle (rad).
533 * \param el_rate Rate of change of true elevation angle (rad/s).
534 *
535 * \return Rate of change of apparent elevation (rad/s).
536 */
537double predict_apparent_elevation_rate(double el, double el_rate);
538
539/*!
540 * \brief Calculate apparent elevation rate.
541 *
542 * Corrects for different atmospheric pressures and temperatures.
543 *
544 * \param el True elevation angle (rad).
545 * \param el_rate Rate of change of true elevation angle (rad/s).
546 * \param pressure Atmospheric pressure (kPa).
547 * \param temp Temperature (deg C).
548 *
549 * \return Rate of change of apparent elevation (rad/s).
550 */
551double predict_apparent_elevation_rate_ext(double el, double el_rate, double pressure, double temp);
552
553#ifdef __cplusplus
554}
555#endif
556
557#endif //_PREDICT_H_
558