predictlib.c 5.15 KB
Newer Older
Jeroen Vreeken's avatar
Jeroen Vreeken committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
#include "sgp4sdp4.h"
#include "predictlib.h"

struct predict {
	sat_t sat;
	double lat;
	double lon;
	double alt;
};

/** \brief SGP4SDP4 driver for doing AOS/LOS calculations.
 *  \param sat Pointer to the satellite data.
 *  \param qth Pointer to the QTH data.
 *  \param t The time for calculation (Julian Date)
 *
 */
void
predict_calc (struct predict *predict, double t)
{
	sat_t *sat = &predict->sat;

    obs_set_t     obs_set;
    geodetic_t    sat_geodetic;
    geodetic_t    obs_geodetic;
    double        age;


    obs_geodetic.lon = predict->lon * de2ra;
    obs_geodetic.lat = predict->lat * de2ra;
    obs_geodetic.alt = predict->alt / 1000.0;
    obs_geodetic.theta = 0;


    sat->jul_utc = t;
    sat->tsince = (sat->jul_utc - sat->jul_epoch) * xmnpda;

    /* call the norad routines according to the deep-space flag */
    if (sat->flags & DEEP_SPACE_EPHEM_FLAG)
        SDP4 (sat, sat->tsince);
    else
        SGP4 (sat, sat->tsince);

    Convert_Sat_State (&sat->pos, &sat->vel);

    /* get the velocity of the satellite */
    Magnitude (&sat->vel);
    sat->velo = sat->vel.w;
    Calculate_Obs (sat->jul_utc, &sat->pos, &sat->vel, &obs_geodetic, &obs_set);
    Calculate_LatLonAlt (sat->jul_utc, &sat->pos, &sat_geodetic);

    while (sat_geodetic.lon < -pi)
        sat_geodetic.lon += twopi;

    while (sat_geodetic.lon > (pi))
        sat_geodetic.lon -= twopi;

    sat->az = Degrees (obs_set.az);
    sat->el = Degrees (obs_set.el);
    sat->range = obs_set.range;
    sat->range_rate = obs_set.range_rate;
    sat->ssplat = Degrees (sat_geodetic.lat);
    sat->ssplon = Degrees (sat_geodetic.lon);
    sat->alt = sat_geodetic.alt;
    sat->ma = Degrees (sat->phase);
    sat->ma *= 256.0/360.0;
    sat->phase = Degrees (sat->phase);

    /* same formulas, but the one from predict is nicer */
    //sat->footprint = 2.0 * xkmper * acos (xkmper/sat->pos.w);
    sat->footprint = 12756.33 * acos (xkmper / (xkmper+sat->alt));
    age = sat->jul_utc - sat->jul_epoch;
    sat->orbit = (long) floor((sat->tle.xno * xmnpda/twopi +
                    age * sat->tle.bstar * ae) * age +
                    sat->tle.xmo/twopi) + sat->tle.revnum - 1;
}

void gtk_sat_data_init_sat (struct predict *predict)
{
    geodetic_t obs_geodetic;
    obs_set_t obs_set;
    geodetic_t sat_geodetic;
    double jul_utc, age;
	sat_t *sat = &predict->sat;

        sat->flags = 0;

        select_ephemeris (sat);

        /* initialise variable fields */
        sat->jul_utc = 0.0;
        sat->tsince = 0.0;
        sat->az = 0.0;
        sat->el = 0.0;
        sat->range = 0.0;
        sat->range_rate = 0.0;
        sat->ra = 0.0;
        sat->dec = 0.0;
        sat->ssplat = 0.0;
        sat->ssplon = 0.0;
        sat->alt = 0.0;
        sat->velo = 0.0;
        sat->ma = 0.0;
        sat->footprint = 0.0;
        sat->phase = 0.0;
        sat->aos = 0.0;
        sat->los = 0.0;

    jul_utc = Julian_Date_of_Epoch (sat->tle.epoch); // => tsince = 0.0
    sat->jul_epoch = jul_utc;

    /* initialise observer location */
    obs_geodetic.lon = 0.0;
    obs_geodetic.lat = 0.0;
    obs_geodetic.alt = 0.0;
    obs_geodetic.theta = 0;


    /* execute computations */
    if (sat->flags & DEEP_SPACE_EPHEM_FLAG)
        SDP4 (sat, 0.0);
    else
        SGP4 (sat, 0.0);

    /* scale position and velocity to km and km/sec */
    Convert_Sat_State (&sat->pos, &sat->vel);

    /* get the velocity of the satellite */
    Magnitude (&sat->vel);
    sat->velo = sat->vel.w;
    Calculate_Obs (jul_utc, &sat->pos, &sat->vel, &obs_geodetic, &obs_set);
    Calculate_LatLonAlt (jul_utc, &sat->pos, &sat_geodetic);

    while (sat_geodetic.lon < -pi)
        sat_geodetic.lon += twopi;
    
    while (sat_geodetic.lon > (pi))
        sat_geodetic.lon -= twopi;

    sat->az = Degrees (obs_set.az);
    sat->el = Degrees (obs_set.el);
    sat->range = obs_set.range;
    sat->range_rate = obs_set.range_rate;
    sat->ssplat = Degrees (sat_geodetic.lat);
    sat->ssplon = Degrees (sat_geodetic.lon);
    sat->alt = sat_geodetic.alt;
    sat->ma = Degrees (sat->phase);
    sat->ma *= 256.0/360.0;
    sat->footprint = 2.0 * xkmper * acos (xkmper/sat->pos.w);
    age = 0.0;
    sat->orbit = (long) floor((sat->tle.xno * xmnpda/twopi +
                               age * sat->tle.bstar * ae) * age +
                              sat->tle.xmo/twopi) + sat->tle.revnum - 1;

    /* orbit type */
//    sat->otype = get_orbit_type (sat);
}

struct predict *predict_create(char *tle, double lat, double lon, double alt)
{
	struct predict *predict;
	char *tle_dup = strdup(tle);
	int i;
	
	predict = calloc(sizeof(struct predict), 1);
	
	for (i = 0; tle[i] != '1' && i < strlen(tle); i++);
	tle += i;
	
        Convert_Satellite_Data (tle, &predict->sat.tle);

	predict->lat = lat;
	predict->lon = lon;
	predict->alt = alt;
	
	gtk_sat_data_init_sat (predict);

	free(tle_dup);
	
	return predict;
}

bool predict_tle_check(char *tle)
{
	return !Good_Elements(tle);
}

void predict_free(struct predict *predict)
{
	free(predict);
}

/* calculate for a given julian date */
int predict_calc_azel(struct predict *predict, double t, double *az, double *el)
{
	predict_calc(predict, t);

	*az = predict->sat.az;
	*el = predict->sat.el;
	
	return 0;
}