1 #include <stdio.h>
2 #include <vector>
3 #include <math.h>
4 #include <unistd.h>
5 #include <stdlib.h>
6 #include <predict/predict.h>
7 #include "testcase_reader.h"
8 #include <iostream>
9 using namespace std;
10 
11 int runtest(const char *filename);
12 
main(int argc,char ** argv)13 int main(int argc, char **argv)
14 {
15 	// Check arguments
16 	if (argc < 2) {
17 		cout << "Usage: " << argv[0] << " <testfiles>" << endl;
18 		return -1;
19 	}
20 
21 	// Test all provided test files
22 	int retval = 0;
23 	for (int i = 1; i < argc; ++i) {
24 		if (runtest(argv[i])) {
25 			cout << argv[i] << ": failed" << endl;
26 			retval = -1;
27 		} else {
28 			cout << argv[i] << ": OK" << endl;
29 		}
30 	}
31 
32 	return retval;
33 }
34 
35 /**
36  * Predict AOS/LOS for several passes after the start time, and check that the times are consistent with each other, that elevation is non-zero in the center
37  * and that the elevation rates are correct at each point.
38  *
39  * \param orbital_elements Orbital elements
40  * \param observer Observer
41  * \param start_time Start time
42  * \return 0 on success, -1 on failure
43  **/
44 int aoslos_timepoint_consistency_test(predict_orbital_elements_t *orbital_elements, predict_observer_t *observer, double start_time);
45 
runtest(const char * filename)46 int runtest(const char *filename)
47 {
48 	// Load testcase
49 	TestCaseReader testcase;
50 	testcase.loadFromFile(filename);
51 	if (!(testcase.containsValidData() && (testcase.containsValidQth()) && (testcase.containsValidTLE()))) {
52 		fprintf(stderr, "Failed to load testfile: %s\n", filename);
53 		return -1;
54 	}
55 
56 	// Get TLE
57 	char *tle[2];
58 	testcase.getTLE(tle);
59 
60 	// Create orbit object
61 	predict_orbital_elements_t *orbital_elements = predict_parse_tle(tle[0], tle[1]);
62 
63 	// Create observer object
64 	predict_observer_t *obs = predict_create_observer("test", testcase.latitude()*M_PI/180.0, testcase.longitude()*M_PI/180.0, testcase.altitude());
65 	if (!obs) {
66 		fprintf(stderr, "Failed to initialize observer!");
67 		return -1;
68 	}
69 
70 	// Use first available time as start time for AOS/LOS finding
71 	double start_time = predict_to_julian(testcase.data()[0][0]);
72 
73 	//check whether the pass makes sense wrt elevation and/or elevation rate at start, end and middle of pass
74 	if (aoslos_timepoint_consistency_test(orbital_elements, obs, start_time) != 0) {
75 		return -1;
76 	}
77 
78 	return 0;
79 }
80 
81 //The tolerance used for fine-tuning the elevations in the aoslos-functions
82 #define ELEVATION_ZERO_TOLERANCE 0.3
83 
84 //Expected time between two passes should at least be 20 minutes
85 #define PASS_TIME_THRESH (20.0/(24*60))
86 
aoslos_timepoint_consistency_test(predict_orbital_elements_t * orbital_elements,predict_observer_t * observer,double start_time)87 int aoslos_timepoint_consistency_test(predict_orbital_elements_t *orbital_elements, predict_observer_t *observer, double start_time)
88 {
89 	if (predict_is_geosynchronous(orbital_elements)) {
90 		return 0;
91 	}
92 
93 	struct predict_observation aos, los;
94 	aos = predict_next_aos(observer, orbital_elements, start_time);
95 	double aos_time = aos.time;
96 	los = predict_next_los(observer, orbital_elements, aos_time);
97 	double los_time = los.time;
98 
99 	struct predict_position orbit;
100 	struct predict_observation observation;
101 
102 	const int NUM_PASSES = 10;
103 	for (int i=0; i < NUM_PASSES; i++) {
104 		if (los_time <= aos_time) {
105 			fprintf(stderr, "los time not strictly larger than aos time: %f %f\n", aos_time, los_time);
106 			return -1;
107 		}
108 
109 		predict_orbit(orbital_elements, &orbit, aos_time);
110 		predict_observe_orbit(observer, &orbit, &observation);
111 		double elevation_rate_aos = observation.elevation_rate;
112 
113 		predict_orbit(orbital_elements, &orbit, los_time);
114 		predict_observe_orbit(observer, &orbit, &observation);
115 		double elevation_rate_los = observation.elevation_rate;
116 
117 		if ((elevation_rate_aos <= 0) || (elevation_rate_los >= 0)) {
118 			fprintf(stderr, "Elevation rates do not have correct signs for aos and los times: %f (%f) %f (%f)\n", elevation_rate_aos, aos_time, elevation_rate_los, los_time);
119 			return -1;
120 		}
121 
122 		double midpass = (los_time - aos_time)/2.0 + aos_time;
123 		predict_orbit(orbital_elements, &orbit, midpass);
124 		predict_observe_orbit(observer, &orbit, &observation);
125 
126 		if (observation.elevation <= 0.0) {
127 			fprintf(stderr, "Elevation is negative during the middle of a pass: %f %f %f\n", aos_time, midpass, los_time);
128 			return -1;
129 		}
130 
131 		double prev_los_time = los_time;
132 		double prev_aos_time = aos_time;
133 		struct predict_observation aos, los;
134 		aos = predict_next_aos(observer, orbital_elements, los_time);
135 		aos_time = aos.time;
136 		los = predict_next_los(observer, orbital_elements, aos_time);
137 		los_time = los.time;
138 
139 		if (!((los_time > aos_time) && (los_time > prev_los_time) && (los_time > prev_aos_time) && (aos_time > prev_los_time))) {
140 			fprintf(stderr, "New AOS/LOS not strictly larger than previous AOS/LOS\n");
141 			return -1;
142 		}
143 
144 		if ((aos_time - prev_los_time) < PASS_TIME_THRESH) {
145 			fprintf(stderr, "Time between passes not significantly different: %f minutes\n", (aos_time - prev_los_time)*24*60);
146 			return -1;
147 		}
148 
149 		//check that time between LOS time and next AOS time produces negative elevations
150 		double timestep = 1.0/(24.0*50);
151 		double curr_time = prev_los_time + timestep;
152 		while (curr_time < aos_time - timestep) {
153 			predict_orbit(orbital_elements, &orbit, curr_time);
154 			predict_observe_orbit(observer, &orbit, &observation);
155 			if (observation.elevation*180.0/M_PI >= ELEVATION_ZERO_TOLERANCE) {
156 				fprintf(stderr, "AOS/LOS failed consistency check between two passes: %f %f %f %f\n", curr_time, observation.elevation, prev_los_time, aos_time);
157 				return -1;
158 			}
159 			curr_time += timestep;
160 		}
161 	}
162 
163 	return 0;
164 }
165