-
Notifications
You must be signed in to change notification settings - Fork 39
/
sat_id.cpp
1773 lines (1602 loc) · 65.2 KB
/
sat_id.cpp
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
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* Copyright (C) 2018, Project Pluto. See LICENSE. */
/*
sat_id.cpp 8 March 2003, with updates as listed below
An example 'main' function illustrating how to find which satellite(s)
are within a given radius of a given RA/dec, as seen from a given
point. The code reads in a file of observations in MPC format (name
provided as the first command-line argument). For example:
sat_id mpc_obs.txt
would hunt through the file 'mpc_obs.txt' for MPC-formatted
observations. It would then read the file 'alldat.tle', looking
for corresponding satellites within .2 degrees of said observations.
It then spits out the original file, with satellite IDs added
(when found) after each observation line. For each IDed satellite,
the international and NORAD designations are given, along with
its angular distance from the search point, position angle of
motion, and apparent angular rate of motion in arcminutes/second
(or, equivalently, degrees/minute). */
/* 2 July 2003: fixed the day/month/year to JD part of 'get_mpc_data()'
so it will work for all years >= 0 (previously, it worked for years
2000 to 2099... plenty for the practical purpose of ID'ing recently-found
satellites, but this is also an 'example' program.) */
/* 3 July 2005: revised the check on the return value for parse_elements().
Now elements with bad checksums won't be rejected. */
/* 23 June 2006: after comment from Eric Christensen, revised to use
names 'ObsCodes.html' or 'ObsCodes.htm', with 'stations.txt' being a
third choice. Also added the '-a' command line switch to cause the program
to show all lines from input (default is now that only MPC astrometric
input gets echoed.) */
/* 30 June 2006: further comment from Eric Christensen: when computing
object motion from two consecutive observations, if the second one has
a date/time preceding the first, you get a negative rate of motion that's
off by 180 degrees. Fixed this. */
/* 17 Nov 2006: artificial satellite data is now being provided in a
file named 'ALL_TLE.TXT'. I've modified the default TLE to match.
(Note : since modified to be read from 'tle_list.txt', as found in
the https://www.github.com/Bill-Gray/tles repository. This allows
for multiple TLEs to be read.) */
/* 22 Oct 2012: minor cosmetic changes, such as making constant variables
of type 'const', updating URL for the MPC station code file, adding a
comment or two. */
/* 7 Jan 2013: revised output to show satellite name if available, plus
the eccentricity, orbital period, and inclination. */
/* 2013 Dec 8: revised to pay attention to "# MJD" and "#Ephem start"
lines, for files that contain many TLEs covering different time spans
for the same object. I sometimes create such files; when that happens,
for each observation, only the TLE(s) covering that observation's time
should be used, and the others are suppressed. */
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <ctype.h>
#include <stdlib.h>
#include <assert.h>
#if defined( _WIN32) || defined( __WATCOMC__)
#include <malloc.h> /* for alloca() prototype */
#else
#include <unistd.h>
#endif
#if defined(_MSC_VER) && _MSC_VER < 1900
/* For older MSVCs, we have to supply our own */
/* snprintf(). See snprintf.cpp for details. */
int snprintf( char *string, const size_t max_len, const char *format, ...);
#endif
#if defined( __has_include)
#if !__has_include(<watdefs.h>)
#error \
"You need the 'lunar' library (https://www.github.com/Bill-Gray/lunar).\
See https://github.com/Bill-Gray/sat_code/issues/2 for a fix to this."
/* This would normally be followed by dozens of errors. GCC, */
/* at least, stops completely when it can't find a system */
/* include file. */
#include <stop_cascading_errors>
#endif
#endif
#include "norad.h"
#include "observe.h"
#include "watdefs.h"
#include "mpc_func.h"
#include "afuncs.h"
#include "date.h"
#include "sat_util.h"
#include "stringex.h"
#define OBSERVATION struct observation
OBSERVATION
{
char text[81];
double jd, ra, dec;
double observer_loc[3];
};
typedef struct
{
double dist, ra, dec, motion_rate, motion_pa;
int norad_number;
char intl_desig[9];
char text[54];
bool in_shadow;
} match_t;
typedef struct
{
OBSERVATION *obs;
size_t idx1, idx2, n_obs, n_matches;
double speed;
match_t *matches;
} object_t;
/* When we encounter a line for a spacecraft-based observation's offset
from the geocenter (a "second line" in the MPC's punched-card format system),
we store that offset in an offset_t structure. Once the file has been
read in, we go back and match the offsets to their corresponding observations.
If the lines aren't in proper order (which happens), we'll fix it, and
can detect errors where we get one line but not the other. */
typedef struct
{
double jd, posn[3];
char mpc_code[4];
} offset_t;
#define PI 3.1415926535897932384626433832795028841971693993751058209749445923
#define TIME_EPSILON (1./86400.)
#define EARTH_MAJOR_AXIS 6378137.
#define EARTH_MINOR_AXIS 6356752.
#define is_power_of_two( X) (!((X) & ((X) - 1)))
/* For testing purposes, I sometimes want _only_ "my" TLEs (not the
'classified' or Space-Watch TLEs) to be used. This is invoked with -s. */
static bool my_tles_only = false;
#if !defined( ON_LINE_VERSION) && !defined( _WIN32)
#define CSI "\x1b["
#define OSC "\x1b]"
#define REVERSE_VIDEO CSI "0m" CSI "7m"
#define NORMAL_VIDEO CSI "0m"
#else
#define REVERSE_VIDEO
#define NORMAL_VIDEO
#endif
/* Sputnik 1 was launched on 1957 Oct 4. There is no point in
checking for satellites before that date. TLEs, and therefore
this program, won't work after 2057 Jan 1. */
const double oct_4_1957 = 2436115.5;
const double jan_1_2057 = 2472364.5;
static int get_mpc_data( OBSERVATION *obs, const char *buff)
{
obs->jd = extract_date_from_mpc_report( buff, NULL);
if( obs->jd < oct_4_1957 || obs->jd > jan_1_2057)
return( -1); /* not an 80-column MPC record */
if( get_ra_dec_from_mpc_report( buff, NULL, &obs->ra, NULL,
NULL, &obs->dec, NULL))
if( 's' != buff[14] && 'v' != buff[14]) /* satellite offsets and */
return( -2); /* roving obs lat/lons won't have RA/decs */
assert( strlen( buff) < sizeof( obs->text));
strncpy( obs->text, buff, sizeof( obs->text));
obs->text[80] = '\0';
return( 0);
}
static int extract_csv_value( char *obuff, const char *ibuff, int idx, size_t obuff_len)
{
while( idx && *ibuff)
if( *ibuff++ == ',')
idx--;
if( !*ibuff)
return( -1);
if( *ibuff == '"')
ibuff++;
obuff_len--;
while( obuff_len-- && *ibuff && *ibuff != ',' && *ibuff != '"')
*obuff++ = *ibuff++;
*obuff = '\0';
return( 0);
}
/* If we encounter 'field' data -- we're determining which artsats
are in a field, not which ones match moving objects -- the output
is very different : */
static bool field_mode = false;
/* An imaging field is specified as an image name, date/time, RA and
dec in decimal degrees, and MPC obscode, all comma-separated, such as
MyImage,2023 nov 17 03:14:15.9,271.818,-14.142,T05
"Another image",2023-11-17.314159,21 16 58.21,"+31 41 59.26","V00"
In the first, the RA and dec are in decimal degrees; in the second,
sexagesimal notation is used and the RA is assumed to be in hours,
minutes, and seconds. You need the commas, but can skip the
quotation marks. */
static int get_field_data( OBSERVATION *obs, const char *buff)
{
int rval = -1;
size_t n_commas, i;
for( i = n_commas = 0; buff[i] && n_commas < 5; i++)
if( buff[i] == ',')
n_commas++;
if( n_commas == 4 && strlen( buff) < 80)
{
char tbuff[80];
int bytes_read;
extract_csv_value( tbuff, buff, 1, sizeof( tbuff));
obs->jd = get_time_from_string( 0., tbuff, 0, NULL);
if( obs->jd > oct_4_1957 && obs->jd < jan_1_2057)
{
extract_csv_value( tbuff, buff, 4, sizeof( tbuff));
if( strlen( tbuff) != 3) /* MPC code must be three characters */
return( -1);
strcpy( obs->text + 77, tbuff);
extract_csv_value( tbuff, buff, 2, sizeof( tbuff));
obs->ra = get_ra_from_string( tbuff, &bytes_read);
if( bytes_read <= 0)
return( -1);
extract_csv_value( tbuff, buff, 3, sizeof( tbuff));
obs->dec = get_dec_from_string( tbuff, &bytes_read);
if( bytes_read <= 0)
return( -1);
extract_csv_value( obs->text, buff, 0, sizeof( obs->text));
field_mode = true;
rval = 0;
}
}
return( rval);
}
/* This loads up the file 'ObsCodes.html' into memory on its first call.
Then, given an MPC code, it finds the corresponding line and copies
it into 'station_code_data'. It looks in several places for the file;
if you've installed Find_Orb, it should be able to get it from the
~/.find_orb directory. It also checks for the truncated 'ObsCodes.htm'
version of the file. */
int verbose = 0;
static bool check_all_tles = false;
static int get_station_code_data( char *station_code_data,
const char *mpc_code)
{
static char *cached_data, *cached_ptr;
if( !mpc_code) /* freeing memory */
{
if( cached_data)
free( cached_data);
cached_data = cached_ptr = NULL;
return( 0);
}
*station_code_data = '\0';
if( !cached_data)
{
const char *filenames[2] = { "ObsCodes.html", "ObsCodes.htm" };
FILE *ifile = NULL;
size_t size;
int i;
for( i = 0; !ifile && i < 2; i++)
ifile = local_then_config_fopen( filenames[i], "rb");
if( !ifile)
{
printf( "Failed to find MPC station list 'ObsCodes.html'\n");
printf( "This can be downloaded at:\n\n");
printf( "http://www.minorplanetcenter.org/iau/lists/ObsCodes.html\n");
exit( -3);
}
fseek( ifile, 0L, SEEK_END);
size = (size_t)ftell( ifile);
fseek( ifile, 0L, SEEK_SET);
cached_data = (char *)malloc( size + 1);
if( fread( cached_data, 1, size, ifile) != size)
{
printf( "Failed to read station file\n");
exit( -4);
}
fclose( ifile);
if( verbose)
printf( "Station codes: %u bytes read\n", (unsigned)size);
ifile = local_then_config_fopen( "rovers.txt", "rb");
if( ifile)
{ /* if 'rovers.txt' is available, */
size_t size2; /* append its data to ObsCodes.htm */
fseek( ifile, 0L, SEEK_END);
size2 = (size_t)ftell( ifile);
fseek( ifile, 0L, SEEK_SET);
cached_data = (char *)realloc( cached_data, size + size2 + 1);
if( fread( cached_data + size, 1, size2, ifile) != size2)
{
printf( "Failed to read 'rovers.txt'\n");
exit( -4);
}
fclose( ifile);
size += size2;
}
cached_data[size] = '\0';
}
if( !cached_ptr || memcmp( cached_ptr, mpc_code, 3))
{
char search_buff[5];
snprintf_err( search_buff, sizeof( search_buff), "\n%.3s", mpc_code);
cached_ptr = strstr( cached_data, search_buff);
if( cached_ptr)
cached_ptr++;
}
if( cached_ptr)
{
size_t i;
for( i = 0; cached_ptr[i] >= ' '; i++)
station_code_data[i] = cached_ptr[i];
station_code_data[i] = '\0';
}
else
{
static char *codes_already_reported = NULL;
static size_t n_reported = 0;
if( codes_already_reported && strstr( codes_already_reported, mpc_code))
return( -1); /* we've already reported this as 'unknown' */
n_reported++;
codes_already_reported = (char *)realloc( codes_already_reported,
n_reported * 4 + 1);
if( n_reported == 1) /* realloc( ) doesn't initialize to zero */
*codes_already_reported = '\0';
strcat( codes_already_reported, mpc_code);
strcat( codes_already_reported, " ");
printf( "Station code '%s' not found.\n", mpc_code);
if( n_reported == 1) /* only really need the following text once */
{
#ifdef ON_LINE_VERSION
printf( "If this is a new MPC code, it could be that this service needs to be\n");
printf( "updated to know about it. Please contact pluto at projectpluto.com so\n");
printf( "I can fix this.\n");
#else
printf( "If this is a new MPC code, you may need to get this file:\n");
printf( "http://www.minorplanetcenter.org/iau/lists/ObsCodes.html\n");
printf( "and replace the existing ObsCodes.html.\n");
#endif
}
}
return( cached_ptr ? 0 : -1);
}
/* Spacecraft-based observations have the spacecraft position given on a
second line. See https://www.minorplanetcenter.net/iau/info/SatelliteObs.html
for details. That documentation is a little outdated; the decimal point
can actually be anywhere within the twelve-byte field. */
static double get_spacecraft_coord( const char *iptr)
{
char tbuff[13];
size_t i = 0;
double rval;
strncpy( tbuff, iptr, 13);
tbuff[12] = '\0';
while( tbuff[i] != '+' && tbuff[i] != '-' && tbuff[i])
i++;
if( !tbuff[i])
rval = 0.;
else
{
rval = atof( tbuff + i + 1);
if( tbuff[i] == '-')
rval = -rval;
}
return( rval);
}
/* Code to check if a 'second line' (v for roving observer or s for
spacecraft observation) matches a 'first line' (the one that has the actual
astrometry). If the date and obscode match, and the observation doesn't
already have a non-geocentric position set, then we have a match. */
static bool offset_matches_obs( const offset_t *offset, const OBSERVATION *obs)
{
if( offset->jd == obs->jd && !obs->observer_loc[0]
&& !strcmp( offset->mpc_code, obs->text + 77)
&& !obs->observer_loc[1] && !obs->observer_loc[2])
return( true);
else
return( false);
}
/* (XXX) locations are specified with text such as
COM Long. 239 18 45 E, Lat. 33 54 11 N, Alt. 100m, Google Earth */
static char xxx_location[80];
static int set_observer_location( OBSERVATION *obs)
{
mpc_code_t code_data;
int rval;
if( memcmp( obs->text + 77, "XXX", 3))
{
char station_data[100];
rval = get_station_code_data( station_data, obs->text + 77);
if( !rval)
get_mpc_code_info( &code_data, station_data);
}
else
{
static bool first_time = true;
rval = get_xxx_location_info( &code_data, xxx_location);
if( rval && first_time)
printf( "WARNING: (XXX) observations won't be handled, because the\n"
"observatory location was either missing or incorrectly formatted.\n"
"See https://www.projectpluto.com/xxx.htm for information on how\n"
"this line should be handled.\n");
first_time = false;
}
if( !rval)
{
observer_cartesian_coords( obs->jd, code_data.lon,
code_data.rho_cos_phi, code_data.rho_sin_phi, obs->observer_loc);
j2000_to_epoch_of_date( obs->jd, &obs->ra, &obs->dec);
}
return( rval);
}
/* Loads up MPC-formatted 80-column observations from a file. Makes
a pass to find out how many observations there are, allocates space
for them, then reads them again to actually load the observations. */
static const char *_target_desig;
static OBSERVATION *get_observations_from_file( FILE *ifile, size_t *n_found,
const double t_low, const double t_high)
{
OBSERVATION *rval = NULL, obs;
void *ades_context = init_ades2mpc( );
char buff[400];
size_t count = 0, n_allocated = 0, n_offsets = 0, i;
offset_t *offsets = NULL;
int n_errors_found = 0;
assert( ades_context);
memset( &obs, 0, sizeof( OBSERVATION));
while( fgets_with_ades_xlation( buff, sizeof( buff), ades_context, ifile))
if( (!get_mpc_data( &obs, buff) || !get_field_data( &obs, buff))
&& obs.jd > t_low && obs.jd < t_high
&& (!_target_desig || strstr( buff, _target_desig)))
{
if( buff[14] == 's' || buff[14] == 'v')
{ /* satellite obs or roving observer */
offset_t toff;
toff.jd = obs.jd;
strlcpy_error( toff.mpc_code, buff + 77);
if( buff[14] == 's')
{
for( i = 0; i < 3; i++)
{
toff.posn[i] = get_spacecraft_coord( buff + 34 + i * 12);
if( buff[32] == '2') /* posn is actually in AU */
toff.posn[i] /= AU_IN_KM;
if( !toff.posn[i])
{
if( n_errors_found++ < 10)
fprintf( stderr, REVERSE_VIDEO "Malformed satellite offset\n%s\n"
NORMAL_VIDEO, buff);
i = 3;
}
}
}
else
{
double lat, lon, alt_in_meters, rho_sin_phi, rho_cos_phi;
if( sscanf( buff + 34, "%lf %lf %lf", &lon, &lat, &alt_in_meters) != 3)
if( n_errors_found++ < 10)
fprintf( stderr, "Couldn't parse roving observer:\n%s\n", buff);
lon *= PI / 180.;
lat *= PI / 180.;
lat_alt_to_parallax( lat, alt_in_meters,
&rho_cos_phi, &rho_sin_phi,
EARTH_MAJOR_AXIS, EARTH_MINOR_AXIS);
observer_cartesian_coords( obs.jd, lon, rho_cos_phi,
rho_sin_phi, toff.posn);
}
n_offsets++;
offsets = (offset_t *)realloc( offsets, n_offsets * sizeof( offset_t));
offsets[n_offsets - 1] = toff;
}
else if( !set_observer_location( &obs))
{
if( count == n_allocated)
{
n_allocated += 10 + n_allocated / 2;
rval = (OBSERVATION *)realloc( rval,
(n_allocated + 1) * sizeof( OBSERVATION));
}
rval[count] = obs;
count++;
}
}
else if( !strncmp( buff, "COM verbose ", 12))
verbose = atoi( buff + 12) + 1;
else if( !strcmp( buff, "COM check tles"))
check_all_tles = true;
else if( !memcmp( buff, "COM Long.", 9))
strlcpy_error( xxx_location, buff);
else if( !strcmp( buff, "COM ignore obs"))
while( fgets_trimmed( buff, sizeof( buff), ifile))
if( !strcmp( buff, "COM end ignore obs"))
break;
*n_found = count;
free_ades2mpc_context( ades_context);
/* for each spacecraft offset, look for the corresponding
observation. If we don't find one, emit a warning. */
for( i = 0; i < n_offsets; i++)
{
size_t j = 0;
while( j < count && !offset_matches_obs( offsets + i, rval + j))
j++;
if( j == count)
{
if( n_errors_found++ < 10)
fprintf( stderr, REVERSE_VIDEO "Unmatched artsat obs for JD %f\n"
NORMAL_VIDEO, offsets[i].jd);
}
else
memcpy( rval[j].observer_loc, offsets[i].posn, 3 * sizeof( double));
}
/* Warn about obs with no offset data, too. */
for( i = 0; i < count; i++)
if( !rval[i].observer_loc[0] && !rval[i].observer_loc[1] && !rval[i].observer_loc[2])
if( n_errors_found++ < 10)
fprintf( stderr, REVERSE_VIDEO "No position for this observation :\n%s\n"
NORMAL_VIDEO, rval[i].text);
free( offsets);
if( n_errors_found >= 10)
fprintf( stderr, "Showing first ten of %d errors\n", n_errors_found);
return( rval);
}
static int id_compare( const OBSERVATION *a, const OBSERVATION *b)
{
return( memcmp( a->text, b->text, 12));
}
static int compare_obs( const void *a, const void *b, void *context)
{
const OBSERVATION *aptr = (const OBSERVATION *)a;
const OBSERVATION *bptr = (const OBSERVATION *)b;
int rval = id_compare( aptr, bptr);
INTENTIONALLY_UNUSED_PARAMETER( context);
if( !rval) /* same IDs? Then sort by JD of observation */
rval = (aptr->jd > bptr->jd ? 1 : -1);
return( rval);
}
/* Copied straight from 'shellsor.cpp' in Find_Orb. See comments there. */
void shellsort_r( void *base, const size_t n_elements, const size_t elem_size,
int (*compare)(const void *, const void *, void *), void *context)
{
#if (defined _GNU_SOURCE || defined __GNU__ || defined __linux)
qsort_r( base, n_elements, elem_size, compare, context);
#else
size_t gap = 250104703;
char *data = (char *)base;
char *pivot = (char *)alloca( elem_size);
while( gap < n_elements)
gap = gap * 8 / 3 + 1;
while( (gap = gap * 3 / 8) != 0)
{
size_t j;
const size_t spacing = elem_size * gap;
for( j = gap; j < n_elements; j++)
{
char *tptr = data + j * elem_size;
char *tptr2 = tptr - spacing;
if( (compare)( tptr2, tptr, context) > 0)
{
memcpy( pivot, tptr, elem_size);
memcpy( tptr, tptr2, elem_size);
tptr = tptr2;
tptr2 -= spacing;
while( tptr2 >= base && (compare)( tptr2, pivot, context) > 0)
{
memcpy( tptr, tptr2, elem_size);
tptr = tptr2;
tptr2 -= spacing;
}
memcpy( tptr, pivot, elem_size);
}
}
}
#endif
}
/* Given a unit vector, this creates a perpendicular xi_vect
in the xy plane and an eta_vect perpendicular to them both. */
static void create_orthogonal_vects( const double *v, double *xi_vect, double *eta_vect)
{
const double tval = sqrt( v[0] * v[0] + v[1] * v[1]);
xi_vect[2] = 0.;
if( !tval) /* 'mid' is directly at a celestial pole */
{
xi_vect[0] = 1.;
xi_vect[1] = 0.;
}
else
{
xi_vect[0] = v[1] / tval;
xi_vect[1] = -v[0] / tval;
}
vector_cross_product( eta_vect, v, xi_vect);
}
/* relative_motion() is intended for situations where you've
computed that an object moved from p1 to p2, and want to
compare that motion vector to an observed motion from p3 to
p4. Initial use cases are in Sat_ID (we have a computed motion
from TLEs and two observed RA/decs) and astcheck (similar, but
the computed positions are from orbital elements). In each case,
you're trying to determine if the observed and computed motions
match to within some threshhold.
I used to do this by computing (delta_RA * cos_dec, delta_dec)
for both observed and computed motions. That works well as long as
the two points are very close together. As they're separated, you
get distortions. This should be more accurate.
Given four points
(ra_dec[0], ra_dec[1]) = starting point object 1
(ra_dec[2], ra_dec[3]) = ending point object 1
(ra_dec[4], ra_dec[5]) = starting point object 2
(ra_dec[6], ra_dec[7]) = ending point object 2
we compute their (xi, eta) sky plane coordinates, using a plane
tangent to the midpoint of the starting locations. That way, any
distortion will affect both ends equally. Then we compute how far
each object moved in the sky plane coordinates. Then we compute
the differences in speed. */
double relative_motion( const double *ra_dec)
{
double mid[3], xi_vect[3], eta_vect[3];
double xi[4], eta[4], v[4][3];
double delta_xi, delta_eta;
int i;
for( i = 0; i < 4; i++)
polar3_to_cartesian( v[i], ra_dec[i + i], ra_dec[i + i + 1]);
for( i = 0; i < 3; i++)
mid[i] = v[0][i] + v[1][i] + v[2][i] + v[3][i];
normalize_vect3( mid);
create_orthogonal_vects( mid, xi_vect, eta_vect);
for( i = 0; i < 4; i++)
{
const double dist = dot_product( mid, v[i]);
xi[i] = dot_product( xi_vect, v[i]) / dist;
eta[i] = dot_product( eta_vect, v[i]) / dist;
}
delta_xi = (xi[0] - xi[1]) - (xi[2] - xi[3]);
delta_eta = (eta[0] - eta[1]) - (eta[2] - eta[3]);
return( sqrt( delta_xi * delta_xi + delta_eta * delta_eta));
}
static double angular_sep( const double delta_ra, const double dec1,
const double dec2, double *posn_ang)
{
double p1[2], p2[2], dist;
p1[0] = 0.;
p1[1] = dec1;
p2[0] = delta_ra;
p2[1] = dec2;
calc_dist_and_posn_ang( p1, p2, &dist, posn_ang);
if( posn_ang)
*posn_ang *= 180. / PI;
return( dist);
}
/* Out of all observations for a given object, this function will pick
two that "best" describe the object's motion. For that purpose, we look
for a pair closest to 'optimal_dist' apart. We also limit the separation
in time to 'max_time_sep'; that's to avoid a situation where the observations
are really close to the optimal distance apart, but are actually from
different orbits.
This code thinks in terms of pairs of observations. If somebody insists
on providing a single observation, we duplicate it. (Unless the '-1'
switch is specified, in which case single observations are just dropped.)
*/
static bool include_singletons = true;
static double find_good_pair( OBSERVATION *obs, const size_t n_obs,
size_t *idx1, size_t *idx2)
{
size_t a, b;
double speed = 0., dt;
const double max_time_sep = 0.1; /* .1 day = 2.4 hr */
double best_score = 1e+30;
*idx1 = *idx2 = 0;
for( b = 0; b < n_obs; b++)
for( a = b + 1; a < n_obs
&& (dt = obs[a].jd - obs[b].jd) < max_time_sep; a++)
if( !memcmp( &obs[a].text[77], &obs[b].text[77], 3))
{
const double optimal_dist = PI / 180.; /* one degree */
const double dist = angular_sep( obs[b].ra - obs[a].ra,
obs[b].dec, obs[a].dec, NULL);
const double score = fabs( dist - optimal_dist);
assert( dt >= .0);
if( best_score > score)
{
best_score = score;
*idx2 = a;
*idx1 = b;
speed = dist / dt;
}
}
speed *= 180. / PI; /* cvt speed from radians/day to deg/day */
speed /= hours_per_day; /* ...then to deg/hour = arcmin/minute */
return( speed);
}
/* Aberration from the Ron-Vondrak method, from Meeus'
_Astronomical Algorithms_, p 153, just the leading terms */
static void compute_aberration( const double t_cen, double *ra, double *dec)
{
const double l3 = 1.7534703 + 628.3075849 * t_cen;
const double sin_l3 = sin( l3), cos_l3 = cos( l3);
const double sin_2l3 = 2. * sin_l3 * cos_l3;
const double cos_2l3 = 2. * cos_l3 * cos_l3 - 1.;
const double x = -1719914. * sin_l3 - 25. * cos_l3
+6434. * sin_2l3 + 28007 * cos_2l3;
const double y = 25. * sin_l3 + 1578089 * cos_l3
+25697. * sin_2l3 - 5904. * cos_2l3;
const double z = 10. * sin_l3 + 684185. * cos_l3
+11141. * sin_2l3 - 2559. * cos_2l3;
const double c = 17314463350.; /* speed of light is 173.1446335 AU/day */
const double sin_ra = sin( *ra), cos_ra = cos( *ra);
*ra -= (y * cos_ra - x * sin_ra) / (c * cos( *dec));
*dec += ((x * cos_ra + y * sin_ra) * sin( *dec) - z * cos( *dec)) / c;
}
static void error_exit( const int exit_code)
{
printf(
"sat_id takes the name of an input file of MPC-formatted (80-column)\n\
astrometry as a command-line argument. It searches for matches between\n\
the observation data and satellites in TLEs specified in 'tle_list.txt'.\n\
By default, matches within .2 degrees are shown.\n\n\
Additional command-line arguments are:\n\
-a YYYYMMDD Only use observations after this time\n\
-b YYYYMMDD Only use observations before this time\n\
-c Check all TLEs for existence\n\
-m (nrevs) Only consider objects with fewer # revs/day (default=6)\n\
-n (NORAD) Only consider objects with this NORAD identifier\n\
-r (radius) Only show matches within this radius in degrees (default=4)\n\
-t (fname) Get TLEs from this filename\n\
-v Verbose output. '-v2' gets still more verboseness.\n\
-y Set tolerance for apparent motion mismatch\n\
-z (rate) Only consider observations above 'rate' deg/hr (default=.001)\n\
\n\
See 'sat_id.txt' for additional details.\n");
exit( exit_code);
}
static int compute_artsat_ra_dec( double *ra, double *dec, double *dist,
const OBSERVATION *optr, tle_t *tle,
const double *sat_params, bool *in_shadow)
{
double pos[3]; /* Satellite position vector */
double sun_xyzr[4], tval;
double t_since = (optr->jd - tle->epoch) * minutes_per_day;
const double j2000 = 2451545.; /* JD 2451545 = 2000 Jan 1.5 */
int sxpx_rval;
if( select_ephemeris( tle))
sxpx_rval = SDP4( t_since, tle, sat_params, pos, NULL);
else
sxpx_rval = SGP4( t_since, tle, sat_params, pos, NULL);
if( sxpx_rval == SXPX_WARN_PERIGEE_WITHIN_EARTH)
sxpx_rval = 0;
if( verbose > 2 && sxpx_rval)
printf( "TLE failed for JD %f: %d\n", optr->jd, sxpx_rval);
get_satellite_ra_dec_delta( optr->observer_loc, pos, ra, dec, dist);
compute_aberration( (optr->jd - j2000) / 36525., ra, dec);
lunar_solar_position( optr->jd, NULL, sun_xyzr);
ecliptic_to_equatorial( sun_xyzr);
tval = dot_product( sun_xyzr, pos);
if( tval < 0. && in_shadow) /* elongation greater than 90 degrees; */
{ /* may be in earth's shadow */
double tvect[3], r2 = 0.;
const double earth_r = EARTH_MAJOR_AXIS / 1000.; /* in km */
size_t i;
tval /= sun_xyzr[3] * sun_xyzr[3];
for( i = 0; i < 3; i++)
{
tvect[i] = pos[i] - sun_xyzr[i] * tval;
r2 += tvect[i] * tvect[i];
}
*in_shadow = (r2 < earth_r * earth_r);
}
else if( in_shadow)
*in_shadow = false;
return( sxpx_rval);
}
static bool is_in_range( const double jd, const double tle_start,
const double tle_range)
{
return( !tle_range || !tle_start ||
(jd >= tle_start && jd <= tle_start + tle_range));
}
/* Determines if we have _any_ observations between the given JDs. If we
don't, we can skip an individual TLE or an entire file. */
static bool got_obs_in_range( const object_t *objs, size_t n_objects,
const double jd_start, const double jd_end)
{
while( n_objects--)
{
OBSERVATION *obs = objs->obs;
size_t i;
if( obs[0].jd < jd_end && obs[objs->n_obs - 1].jd > jd_start)
for( i = 0; i < objs->n_obs; i++)
if( obs[i].jd > jd_start && obs[i].jd < jd_end)
return( true);
objs++;
}
return( false);
}
static int _pack_intl_desig( char *desig_out, const char *desig)
{
size_t i = 0;
int rval = 0;
while( desig[i] && isdigit( desig[i]))
i++;
memset( desig_out, ' ', 8);
desig_out[8] = '\0';
if( i == 5 && isupper( desig[5])) /* already in YYYNNAaa form */
{
memcpy( desig_out, desig, 5);
desig += 5;
}
else if( i == 4 && desig[4] == '-' && isdigit( desig[5])
&& isdigit( desig[6]) && isdigit( desig[7]) && isupper( desig[8]))
{
desig_out[0] = desig[2]; /* desig is in YYYY-NNNAaa form */
desig_out[1] = desig[3];
desig_out[2] = desig[5];
desig_out[3] = desig[6];
desig_out[4] = desig[7];
desig += 8;
}
else /* not a recognized intl desig form */
rval = -1;
if( !rval)
{
desig_out[5] = *desig++;
if( isupper( *desig))
desig_out[6] = *desig++;
if( isupper( *desig))
desig_out[7] = *desig++;
}
return( rval);
}
static int _compare_intl_desigs( const char *desig1, const char *desig2)
{
char odesig1[9], odesig2[9];
_pack_intl_desig( odesig1, desig1);
_pack_intl_desig( odesig2, desig2);
return( strcmp( odesig1, odesig2));
}
/* Code to look through 'sat_xref.txt', if available, and assign NORAD
and international designations to TLEs with only default designations.
See 'sat_xref.txt' and 'eph2tle.cpp' in the Find_Orb repository. */
static int look_up_extended_identifiers( const char *line0, tle_t *tle)
{
static char buff[100];
int match_found = !strcmp( line0, buff + 21);
static bool got_sat_xref_txt = true; /* until proven otherwise */
if( !match_found && got_sat_xref_txt)
{
FILE *ifile = local_then_config_fopen( "sat_xref.txt", "rb");
if( !ifile) /* don't look for it again */
got_sat_xref_txt = false;
else
{
while( !match_found && fgets_trimmed( buff, sizeof( buff), ifile))
match_found = !strcmp( line0, buff + 21);
fclose( ifile);
}
}
if( match_found)
{
tle->norad_number = atoi( buff);
memcpy( tle->intl_desig, buff + 12, 8);
}
return( match_found);
}
/* The international (YYYY-NNNletter(s)) designation and the NORAD
five-digit designation are always shown for a match. If they're in
the name as well, we should remove them so that more of the actual
name gets displayed. */
static void remove_redundant_desig( char *name, const char *desig)
{
size_t len = strlen( desig), i;
while( len && desig[len - 1] == ' ')
len--;
for( i = 0; name[i]; i++)
if( !memcmp( name + i, desig, len))
{
size_t n = len;
if( i >= 3 && name[i - 2] == '=' && name[i - 1] == ' '
&& name[i - 3] == ' ')
{
i -= 3; /* remove preceding '=' */
n += 3;
}
else if( name[i + n] == ' ' && name[i + n + 1] == '='
&& name[i + n + 2] == ' ')
n += 3;
memmove( name + i, name + i + n, strlen( name + i + n) + 1);
i--;
}
}
/* We check astrometry first against TLEs from github.com/Bill-Gray/tles,
then some other sources such as the amateur community's TLEs, and
only then against Space-Track TLEs. If we've already checked an
object against the previous sources, then we really ought not to
check the Space-Track TLEs. For one thing, the mere fact that we've
gone to the effort of computing "our own" TLEs means the Space-Track
TLEs are unreliable (some have poor accuracy; others are not consistently
available).
Therefore, we maintain a list of '# ID:' numbers from tle_list.txt, and when we
get to '# ID off', we figure we're in Space-Track TLE territory. If we
encounter a NORAD number that's in our list, we essentially say :
we already found this object and have handled it. */
static bool already_found_desig( const int curr_norad, size_t n_norad_ids, const int *norad_ids)