198 lines
4.1 KiB
C
198 lines
4.1 KiB
C
#include <stdio.h>
|
|
#include <string.h>
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
#include "sgdp4h.h"
|
|
#include "satutl.h"
|
|
#include <getopt.h>
|
|
|
|
#define LIM 80
|
|
|
|
// Cross product
|
|
xyz_t
|
|
cross (xyz_t a, xyz_t b)
|
|
{
|
|
xyz_t c;
|
|
|
|
c.x = a.y * b.z - a.z * b.y;
|
|
c.y = a.z * b.x - a.x * b.z;
|
|
c.z = a.x * b.y - a.y * b.x;
|
|
|
|
return c;
|
|
}
|
|
|
|
// Read a line of maximum length int lim from file FILE into string s
|
|
int
|
|
fgetline (FILE * file, char *s, int lim)
|
|
{
|
|
int c, i = 0;
|
|
|
|
while (--lim > 0 && (c = fgetc (file)) != EOF && c != '\n')
|
|
s[i++] = c;
|
|
if (c == '\t')
|
|
c = ' ';
|
|
if (c == '\n')
|
|
s[i++] = c;
|
|
s[i] = '\0';
|
|
return i;
|
|
}
|
|
|
|
// Compute Julian Day from Date
|
|
double
|
|
date2mjd (int year, int month, double day)
|
|
{
|
|
int a, b;
|
|
double jd;
|
|
|
|
if (month < 3)
|
|
{
|
|
year--;
|
|
month += 12;
|
|
}
|
|
|
|
a = floor (year / 100.);
|
|
b = 2. - a + floor (a / 4.);
|
|
|
|
if (year < 1582)
|
|
b = 0;
|
|
if (year == 1582 && month < 10)
|
|
b = 0;
|
|
if (year == 1582 && month == 10 && day <= 4)
|
|
b = 0;
|
|
|
|
jd =
|
|
floor (365.25 * (year + 4716)) + floor (30.6001 * (month + 1)) + day + b -
|
|
1524.5;
|
|
|
|
return jd - 2400000.5;
|
|
}
|
|
|
|
// DOY to MJD
|
|
double
|
|
doy2mjd (int year, double doy)
|
|
{
|
|
int month, k = 2;
|
|
double day;
|
|
|
|
if (year % 4 == 0 && year % 400 != 0)
|
|
k = 1;
|
|
|
|
month = floor (9.0 * (k + doy) / 275.0 + 0.98);
|
|
|
|
if (doy < 32)
|
|
month = 1;
|
|
|
|
day =
|
|
doy - floor (275.0 * month / 9.0) + k * floor ((month + 9.0) / 12.0) +
|
|
30.0;
|
|
|
|
return date2mjd (year, month, day);
|
|
}
|
|
|
|
int
|
|
main (int argc, char *argv[])
|
|
{
|
|
int i, arg = 0, imode, n = 100, satno, satname, satno2, imin, flag, status;
|
|
char catalog1[] = "20180121_173045_catalog.txt", catalog2[] =
|
|
"jsc_20180122.txt";
|
|
FILE *file1, *file2;
|
|
orbit_t orb1, orb2;
|
|
double mjd0 = 58140.0, mjd;
|
|
xyz_t *r, *v, r0, v0, d, n1, n2;
|
|
double dr, dv, dn, drmin, dvmin, mjdmin;
|
|
char line0[LIM], line1[LIM], line2[LIM];
|
|
|
|
// Reroute stderr
|
|
freopen ("/tmp/stderr.txt", "w", stderr);
|
|
|
|
// Allocate
|
|
r = (xyz_t *) malloc (sizeof (xyz_t) * n);
|
|
v = (xyz_t *) malloc (sizeof (xyz_t) * n);
|
|
|
|
// Open file
|
|
file1 = fopen (catalog1, "r");
|
|
file2 = fopen (catalog2, "r");
|
|
|
|
// Loop over classfd catalog
|
|
while (read_twoline (file1, 0, &orb1) == 0)
|
|
{
|
|
// Compute positions
|
|
imode = init_sgdp4 (&orb1);
|
|
for (i = 0; i < n; i++)
|
|
satpos_xyz (mjd0 + 2400000.5 + (double) i / 86400.0, &r[i], &v[i]);
|
|
|
|
// Loop over ISON catalog
|
|
flag = 0;
|
|
rewind (file2);
|
|
while (read_twoline (file2, 0, &orb2) == 0)
|
|
{
|
|
mjd = doy2mjd (orb2.ep_year, orb2.ep_day);
|
|
|
|
imode = init_sgdp4 (&orb2);
|
|
satpos_xyz (mjd + 2400000.5, &r0, &v0);
|
|
|
|
// Find nearest offset
|
|
for (i = 0; i < n; i++)
|
|
{
|
|
d.x = r[i].x - r0.x;
|
|
d.y = r[i].y - r0.y;
|
|
d.z = r[i].z - r0.z;
|
|
dr = sqrt (d.x * d.x + d.y * d.y + d.z * d.z);
|
|
d.x = v[i].x - v0.x;
|
|
d.y = v[i].y - v0.y;
|
|
d.z = v[i].z - v0.z;
|
|
dv = sqrt (d.x * d.x + d.y * d.y + d.z * d.z);
|
|
|
|
if (flag == 0 || dr < drmin)
|
|
{
|
|
drmin = dr;
|
|
dvmin = dv;
|
|
satno = orb2.satno;
|
|
imin = i;
|
|
mjdmin = mjd;
|
|
flag = 1;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Compute normals
|
|
rewind (file2);
|
|
read_twoline (file2, satno, &orb2);
|
|
mjd = doy2mjd (orb2.ep_year, orb2.ep_day);
|
|
imode = init_sgdp4 (&orb2);
|
|
satpos_xyz (mjd + 2400000.5, &r0, &v0);
|
|
n1 = cross (r[0], v[0]);
|
|
n2 = cross (r0, v0);
|
|
d.x = n1.x - n2.x;
|
|
d.y = n1.y - n2.y;
|
|
d.z = n1.z - n2.z;
|
|
dn = sqrt (d.x * d.x + d.y * d.y + d.z * d.z);
|
|
|
|
// Find object name
|
|
rewind (file2);
|
|
while (fgetline (file2, line0, LIM) > 0)
|
|
{
|
|
fgetline (file2, line1, LIM);
|
|
fgetline (file2, line2, LIM);
|
|
sscanf (line1, "1 %d", &satno2);
|
|
if (satno == satno2)
|
|
{
|
|
sscanf (line0, "SO %d", &satname);
|
|
}
|
|
}
|
|
|
|
// printf("%05d %05d %8.2f km %9.5f km/s %6.0lf s %8.0f km\n",orb1.satno,satno,drmin,dvmin,(mjd0+(double) imin/86400.0-mjd)*86400,dn);
|
|
printf ("%6d %05d %s | %5d %8.2f km %9.5f km/s %6.0lf s %8.0f km\n",
|
|
satname, orb1.satno, orb1.desig, satno, drmin, dvmin,
|
|
(mjd0 + (double) imin / 86400.0 - mjd) * 86400, dn);
|
|
}
|
|
fclose (file1);
|
|
fclose (file2);
|
|
|
|
// Free
|
|
free (r);
|
|
free (v);
|
|
|
|
return 0;
|
|
}
|