Skip to content

Commit cc05d3d

Browse files
committed
pov mode corrections + part 1 of The Comment Update
1 parent 7887016 commit cc05d3d

File tree

3 files changed

+90
-49
lines changed

3 files changed

+90
-49
lines changed

src/astro.c

Lines changed: 52 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -24,13 +24,15 @@ SatPass passes[MAX_PASSES];
2424
int num_passes = 0;
2525
Satellite *last_pass_calc_sat = NULL;
2626

27+
/* string extraction (sscanf is a bit too beefy for tight TLE loops) */
2728
static double parse_tle_double(const char *str, int start, int len)
2829
{
2930
char buf[32] = {0};
3031
strncpy(buf, str + start, len);
3132
return atof(buf);
3233
}
3334

35+
/* pulls the system clock and mashes it into our custom YYYYDDD.FFFF format */
3436
double get_current_real_time_epoch(void)
3537
{
3638
time_t now = time(NULL);
@@ -40,11 +42,12 @@ double get_current_real_time_epoch(void)
4042
double day_of_year = gmt->tm_yday + 1.0;
4143
double fraction_of_day = (gmt->tm_hour + gmt->tm_min / 60.0 + gmt->tm_sec / 3600.0) / 24.0;
4244

43-
// modified to return full YYYY format to make global time sim less fuckywucky,
44-
// afterwards we can just use the YY format for SGP4 data.
45+
/* modified to return full YYYY format to make global time sim less fuckywucky,
46+
afterwards we just use the YY format for SGP4 data. */
4547
return (year * 1000.0) + day_of_year + fraction_of_day;
4648
}
4749

50+
/* handles year rollover/underflow so the math doesnt explode when looking at past/future passes */
4851
double normalize_epoch(double epoch)
4952
{
5053
int year = (int)(epoch / 1000.0);
@@ -73,14 +76,14 @@ double normalize_epoch(double epoch)
7376
return (year * 1000.0) + day_of_year;
7477
}
7578

76-
// utility to convert normalized epoch format to unix time for sgp4 math
79+
/* utility to convert normalized epoch format to unix time for sgp4 math */
7780
double get_unix_from_epoch(double epoch)
7881
{
7982
epoch = normalize_epoch(epoch);
8083
int year = (int)(epoch / 1000.0);
8184
double day = fmod(epoch, 1000.0);
8285

83-
// pure mathematical unix conversion to avoid OS-level timegm() quantization and stutter
86+
/* pure mathematical unix conversion to avoid OS-level timegm() quantization and stutter */
8487
int y = year - 1;
8588
int leaps_to_year = (y / 4) - (y / 100) + (y / 400);
8689
int leaps_to_1970 = (1969 / 4) - (1969 / 100) + (1969 / 400);
@@ -90,6 +93,7 @@ double get_unix_from_epoch(double epoch)
9093
return unix_days * 86400.0;
9194
}
9295

96+
/* sidereal time keeps the earth spinning under the sats; without this, everything is static */
9397
double epoch_to_gmst(double epoch)
9498
{
9599
double unix_time = get_unix_from_epoch(epoch);
@@ -101,6 +105,7 @@ double epoch_to_gmst(double epoch)
101105
return gmst;
102106
}
103107

108+
/* pretty-print for the ui so humans can actually read the time */
104109
void epoch_to_datetime_str(double epoch, char *buffer)
105110
{
106111
epoch = normalize_epoch(epoch);
@@ -133,6 +138,7 @@ void epoch_to_datetime_str(double epoch, char *buffer)
133138
sprintf(buffer, "%04d-%02d-%02d %02d:%02d:%02.0f UTC", year, month, day, h, m, seconds);
134139
}
135140

141+
/* rips lines from a TLE file and populates the satellite struct */
136142
bool add_satellite_from_tle(const char* line0, const char* line1, const char* line2)
137143
{
138144
if (sat_count >= MAX_SATELLITES) return false;
@@ -163,9 +169,11 @@ bool add_satellite_from_tle(const char* line0, const char* line1, const char* li
163169
double initial_r[3] = {0};
164170
double initial_v[3] = {0};
165171

172+
/* shove the TLE into the sgp4 state machine */
166173
ConvertTLEToSGP4(&sat->satrec, &parsed_objs[0], 0.0, initial_r, initial_v);
167174
free(parsed_objs);
168175

176+
/* manual scraping for the rest of the struct because we like control */
169177
double raw_epoch = parse_tle_double(line1, 18, 14);
170178
int yy = (int)(raw_epoch / 1000.0);
171179
int year = (yy < 57) ? 2000 + yy : 1900 + yy;
@@ -191,6 +199,7 @@ bool add_satellite_from_tle(const char* line0, const char* line1, const char* li
191199
return false;
192200
}
193201

202+
/* bulk loading of celestial junk from flat files */
194203
void load_tle_data(const char *filename)
195204
{
196205
FILE *file = fopen(filename, "r");
@@ -209,7 +218,7 @@ void load_tle_data(const char *filename)
209218
if (strncmp(line0, "# EPOCH:", 8) == 0)
210219
{
211220
unsigned int mask = 0, cust_mask = 0, ret_mask = 0;
212-
// scan past it, variables are unused here but keeps pointer advanced
221+
/* scan past it, variables are unused here but keeps pointer advanced */
213222
if (strstr(line0, "RET_MASK:"))
214223
{
215224
sscanf(line0, "# EPOCH:%*d MASK:%u CUST_MASK:%u RET_MASK:%u", &mask, &cust_mask, &ret_mask);
@@ -221,7 +230,7 @@ void load_tle_data(const char *filename)
221230
}
222231
else
223232
{
224-
rewind(file); // not a header, restart
233+
rewind(file); /* not a header, restart */
225234
}
226235
}
227236

@@ -242,6 +251,7 @@ void load_tle_data(const char *filename)
242251
fclose(file);
243252
}
244253

254+
/* parsing for strings that were likely copy-pasted in a hurry */
245255
void load_manual_tles(AppConfig *config)
246256
{
247257
for (int i = 0; i < config->manual_tle_count; i++)
@@ -264,7 +274,8 @@ void load_manual_tles(AppConfig *config)
264274
}
265275
}
266276

267-
// precalculated unix time passed down to prevent excessyear/day conversions
277+
/* main sgp4 crank; outputs raw ECI coordinates */
278+
/* precalculated unix time passed down to prevent excessyear/day conversions */
268279
Vector3 calculate_position(Satellite *sat, double current_unix)
269280
{
270281
double tsince = (current_unix - sat->epoch_unix) / 60.0;
@@ -282,6 +293,7 @@ Vector3 calculate_position(Satellite *sat, double current_unix)
282293
return pos;
283294
}
284295

296+
/* projects 3D orbital space onto a 2D equirectangular map plane */
285297
void get_map_coordinates(Vector3 pos, double gmst_deg, float earth_offset, float map_w, float map_h, float *out_x, float *out_y)
286298
{
287299
float r = Vector3Length(pos);
@@ -304,6 +316,7 @@ void get_map_coordinates(Vector3 pos, double gmst_deg, float earth_offset, float
304316
*out_y = (v - 0.5f) * map_h;
305317
}
306318

319+
/* finds where the sat hits the high and low points of its orbit in 2D */
307320
void get_apsis_2d(Satellite *sat, double current_time, bool is_apoapsis, double gmst_deg, float earth_offset, float map_w, float map_h, Vector2 *out)
308321
{
309322
(void)gmst_deg;
@@ -328,6 +341,7 @@ void get_apsis_2d(Satellite *sat, double current_time, bool is_apoapsis, double
328341
get_map_coordinates(pos3d, gmst_target, earth_offset, map_w, map_h, &out->x, &out->y);
329342
}
330343

344+
/* predicts the timestamps for the next perigee and apoapsis */
331345
void get_apsis_times(Satellite *sat, double current_time, double *out_peri_unix, double *out_apo_unix)
332346
{
333347
double current_unix = get_unix_from_epoch(current_time);
@@ -351,6 +365,7 @@ void get_apsis_times(Satellite *sat, double current_time, double *out_peri_unix,
351365
*out_apo_unix = get_unix_from_epoch(t_apo);
352366
}
353367

368+
/* bakes the future orbital path into a vertex buffer so sgp4 isnt re-ran every frame */
354369
void update_orbit_cache(Satellite *sat, double current_epoch)
355370
{
356371
double period_days = (2.0 * PI / sat->mean_motion) / 86400.0;
@@ -364,6 +379,7 @@ void update_orbit_cache(Satellite *sat, double current_epoch)
364379
sat->orbit_cached = true;
365380
}
366381

382+
/* converts raw orbital data into azimuth/elevation for a specific ground station */
367383
void get_az_el(Vector3 eci_pos, double gmst_deg, float obs_lat, float obs_lon, float obs_alt, double *az, double *el)
368384
{
369385
double sat_r = Vector3Length(eci_pos);
@@ -376,7 +392,7 @@ void get_az_el(Vector3 eci_pos, double gmst_deg, float obs_lat, float obs_lon, f
376392

377393
double sat_lat = asin(eci_pos.y / sat_r);
378394
double sat_lon_eci = atan2(-eci_pos.z, eci_pos.x);
379-
double theta = (gmst_deg + 0) * DEG2RAD; // assuming earth_rotation_offset handled befor
395+
double theta = (gmst_deg + 0) * DEG2RAD; /* assuming earth_rotation_offset handled befor */
380396
double sat_lon_ecef = sat_lon_eci - theta;
381397

382398
double s_x = sat_r * cos(sat_lat) * cos(sat_lon_ecef);
@@ -411,6 +427,7 @@ void get_az_el(Vector3 eci_pos, double gmst_deg, float obs_lat, float obs_lon, f
411427
*az += 360.0;
412428
}
413429

430+
/* qsort callback to keep passes chronological */
414431
int compare_passes(const void *a, const void *b)
415432
{
416433
const SatPass *p1 = (const SatPass *)a;
@@ -422,6 +439,7 @@ int compare_passes(const void *a, const void *b)
422439
return 0;
423440
}
424441

442+
/* heavy lifting for pass prediction; brute force search with binary search refinement */
425443
void CalculatePasses(Satellite *sat, double start_epoch)
426444
{
427445
num_passes = 0;
@@ -444,7 +462,7 @@ void CalculatePasses(Satellite *sat, double start_epoch)
444462

445463
get_az_el(calculate_position(current_sat, t_unix), gmst, home_location.lat, home_location.lon, home_location.alt, &az, &el);
446464

447-
// back up if happens to already be in a pass to catch the true start
465+
/* back up if happens to already be in a pass to catch the true start */
448466
if (el > 0)
449467
{
450468
for (int i = 0; i < 30 && el > 0; i++)
@@ -472,7 +490,7 @@ void CalculatePasses(Satellite *sat, double start_epoch)
472490
if (!in_pass)
473491
{
474492
in_pass = true;
475-
// binary search to find exact AOS
493+
/* binary search to find exact AOS because stepping by 1min is too crunchy for radio work */
476494
double t_low = t - coarse_step;
477495
double t_high = t;
478496
for (int b = 0; b < 10; b++)
@@ -503,7 +521,7 @@ void CalculatePasses(Satellite *sat, double start_epoch)
503521
if (in_pass)
504522
{
505523
in_pass = false;
506-
// binary search to find exact LOS crossing
524+
/* binary search to find exact LOS crossing */
507525
double t_low = t - coarse_step;
508526
double t_high = t;
509527
for (int b = 0; b < 10; b++)
@@ -525,7 +543,7 @@ void CalculatePasses(Satellite *sat, double start_epoch)
525543
double step = (current_pass.los_epoch - current_pass.aos_epoch) / 399.0;
526544
if (step > 0)
527545
{
528-
current_pass.max_el = -90.0f; // reset to find true max during high-res pass
546+
current_pass.max_el = -90.0f; /* reset to find true max during high-res pass */
529547
for (int k = 0; k < 400; k++)
530548
{
531549
double pt = current_pass.aos_epoch + k * step;
@@ -535,7 +553,7 @@ void CalculatePasses(Satellite *sat, double start_epoch)
535553
get_az_el(calculate_position(current_sat, pt_unix), p_gmst, home_location.lat, home_location.lon, home_location.alt, &p_az, &p_el);
536554
current_pass.path_pts[current_pass.num_pts++] = (Vector2){(float)p_az, (float)p_el};
537555

538-
// Absolutely ensure max elevation is pinpointed precisely
556+
/* ensure max elevation is pinpointed */
539557
if (p_el > current_pass.max_el)
540558
{
541559
current_pass.max_el = (float)p_el;
@@ -579,10 +597,12 @@ void CalculatePasses(Satellite *sat, double start_epoch)
579597
}
580598
}
581599

582-
// Sort overall passes generated by timeframe chronological arrival order
600+
/* make sure the list actually makes sense chronologically */
601+
/* Sort overall passes generated by timeframe chronological arrival order */
583602
qsort(passes, num_passes, sizeof(SatPass), compare_passes);
584603
}
585604

605+
/* formats the internal epoch into a HH:MM:SS string for quick glancing */
586606
void epoch_to_time_str(double epoch, char *str)
587607
{
588608
time_t t = (time_t)get_unix_from_epoch(epoch);
@@ -597,7 +617,7 @@ void epoch_to_time_str(double epoch, char *str)
597617
}
598618
}
599619

600-
// i truly do hope this doesnt drift or something its so eyeballed istg
620+
/* i truly do hope this doesnt drift or something its so eyeballed istg */
601621
Vector3 calculate_sun_position(double current_time_days)
602622
{
603623
double unix_time = get_unix_from_epoch(current_time_days);
@@ -631,6 +651,7 @@ Vector3 calculate_sun_position(double current_time_days)
631651
return Vector3Normalize(pos);
632652
}
633653

654+
/* basic shadow-cone check; tells us if the sat is in the dark (no visual/solar power) */
634655
bool is_sat_eclipsed(Vector3 pos_km, Vector3 sun_dir_norm)
635656
{
636657
float dot = Vector3DotProduct(pos_km, sun_dir_norm);
@@ -645,37 +666,37 @@ Vector3 calculate_moon_position(double current_time_days)
645666
{
646667
double unix_time = get_unix_from_epoch(current_time_days);
647668
double jd = (unix_time / 86400.0) + 2440587.5;
648-
double D_days = jd - 2451545.0; // days since J2000
669+
double D_days = jd - 2451545.0; /* days since J2000 */
649670

650-
// arguments needed later
671+
/* some arguments */
651672
double L_moon = fmod(218.316 + 13.176396 * D_days, 360.0) * DEG2RAD;
652673
double M_moon = fmod(134.963 + 13.064993 * D_days, 360.0) * DEG2RAD;
653674
double F_moon = fmod(93.272 + 13.229350 * D_days, 360.0) * DEG2RAD;
654675

655-
// solar mean anomaly and lunar elongation hehe for perturbation calculations
676+
/* solar mean anomaly and lunar elongation hehe for perturbation calculations */
656677
double M_sun = fmod(357.528 + 0.9856003 * D_days, 360.0) * DEG2RAD;
657678
double D_elong = fmod(297.850 + 12.190749 * D_days, 360.0) * DEG2RAD;
658679

659-
// she perturbate on my variation till I evect
680+
/* she perturbate on my variation till I evect */
660681
double E = 1.0 - 0.002516 * cos(M_sun);
661682
double evection = 1.274 * DEG2RAD * sin(2.0 * D_elong - M_moon);
662683
double variation = 0.658 * DEG2RAD * sin(2.0 * D_elong);
663684
double annual_eq = -0.186 * DEG2RAD * E * sin(M_sun);
664685
double parallactic = -0.035 * DEG2RAD * sin(D_elong);
665686

666-
// apply perturbations to longitude and distance
687+
/* apply perturbations to longitude and distance */
667688
double lambda = L_moon + evection + variation + annual_eq + parallactic + (6.289 * DEG2RAD) * sin(M_moon) + (-0.059 * DEG2RAD) * sin(2.0 * D_elong + M_moon);
668689

669690
double beta = (5.128 * DEG2RAD) * sin(F_moon) + (0.280 * DEG2RAD) * sin(F_moon + M_moon) + (0.277 * DEG2RAD) * sin(F_moon - M_moon) + (0.173 * DEG2RAD) * sin(2.0 * D_elong - F_moon);
670691

671692
double dist_km = 385001.0 - 20905.0 * cos(M_moon) - 3699.0 * cos(2.0 * D_elong - M_moon) - 2956.0 * cos(2.0 * D_elong);
672693

673-
// ecliptic to ECI because that wouldn't worky whatsoever
694+
/* ecliptic to ECI because that wouldn't worky whatsoever */
674695
double x_ecl = dist_km * cos(beta) * cos(lambda);
675696
double y_ecl = dist_km * cos(beta) * sin(lambda);
676697
double z_ecl = dist_km * sin(beta);
677698

678-
// obliquity of the ecliptic (T is centuries since J2000)
699+
/* obliquity of the ecliptic (T is centuries since J2000) */
679700
double T = D_days / 36525.0;
680701
double eps = (23.439291 - 0.0130042 * T) * DEG2RAD;
681702

@@ -687,22 +708,23 @@ Vector3 calculate_moon_position(double current_time_days)
687708
return pos;
688709
}
689710

711+
/* internal helper to figure out straight-line distance to a satellite */
690712
double get_sat_range(Satellite *sat, double epoch, Marker obs)
691713
{
692714
double t_unix = get_unix_from_epoch(epoch);
693715
double theta = epoch_to_gmst(epoch) * DEG2RAD;
694716

695717
Vector3 eci = calculate_position(sat, t_unix);
696718

697-
// direct cartesian rotation (ECI to ECEF)
719+
/* direct cartesian rotation (ECI to ECEF) */
698720
double cos_t = cos(theta);
699721
double sin_t = sin(theta);
700722

701723
double s_x = eci.x * cos_t - eci.z * sin_t;
702724
double s_y = -eci.x * sin_t - eci.z * cos_t;
703725
double s_z = eci.y;
704726

705-
// observer ECEF calc
727+
/* observer ECEF calc */
706728
double lat_rad = obs.lat * DEG2RAD;
707729
double lon_rad = obs.lon * DEG2RAD;
708730
double obs_rad = EARTH_RADIUS_KM + obs.alt / 1000.0;
@@ -712,23 +734,24 @@ double get_sat_range(Satellite *sat, double epoch, Marker obs)
712734
double o_y = obs_rad * cos_lat * sin(lon_rad);
713735
double o_z = obs_rad * sin(lat_rad);
714736

715-
// dist
737+
/* dist */
716738
double dx = s_x - o_x;
717739
double dy = s_y - o_y;
718740
double dz = s_z - o_z;
719741

720742
return sqrt(dx * dx + dy * dy + dz * dz);
721743
}
722744

745+
/* shifts the frequency based on velocity relative to the observer; essential for tuning */
723746
double calculate_doppler_freq(Satellite *sat, double epoch, Marker obs, double base_freq)
724747
{
725-
// line-of-sight range
726-
double dt = 0.1 / 86400.0; // 0.1 seconds step
748+
/* line-of-sight range */
749+
double dt = 0.1 / 86400.0; /* 0.1 seconds step */
727750
double r1 = get_sat_range(sat, epoch - dt, obs);
728751
double r2 = get_sat_range(sat, epoch + dt, obs);
729-
double range_rate = (r2 - r1) / 0.2; // km/s
752+
double range_rate = (r2 - r1) / 0.2; /* km/s */
730753

731-
double c = 299792.458; // in km/s
754+
double c = 299792.458; /* in km/s */
732755
return base_freq * (c / (c + range_rate));
733756
}
734757

0 commit comments

Comments
 (0)