More clean up

pull/10/head
Cees Bassa 2018-02-13 17:35:51 +01:00
parent 1db20f9cb8
commit cd68383bcd
7 changed files with 78 additions and 986 deletions

View File

64
data/sites.txt 100644
View File

@ -0,0 +1,64 @@
# No ID Latitude Longitude Elev Observer
1111 RL 38.9478 -104.5614 2073 Ron Lee
4171 CB 52.8344 6.3785 10 Cees Bassa
4172 LB 52.3713 5.2580 -3 Leo Barhorst
4553 CB 53.3210 -2.2330 86 Cees Bassa
0001 MM 30.3340 -97.7610 160 Mike McCants
0002 MM 30.3138 -97.8661 280 Mike McCants
0030 IA 47.4857 18.9727 400 Ivan Artner
0031 IA 47.5612 18.7366 149 Ivan Artner
0070 BC 53.2233 -0.6003 30 Bob Christy
0100 SG 59.4627 17.9138 0 Sven Grahn
0110 LK 32.5408 -96.8906 200 Lyn Kennedy
0433 GR -33.9406 18.5129 10 Greg Roberts
0434 IR -26.1030 27.9288 1646 Ian Roberts
0435 GR -33.9369 18.5101 25 Greg Roberts
0710 LS 52.3261 10.6756 85 Lutz Schindler
0899 FM -34.8961 -56.1227 30 Fernando Mederos
1056 MK 57.0122 23.9833 4 Martins Keruss
1086 NK 46.4778 30.7556 56 Nikolay Koshkin
1244 AM 44.3932 33.9701 69 Andriy Makeyev
1747 DD 45.7275 -72.3526 191 Daniel Deak
1775 KF 44.6062 -75.6910 200 Kevin Fetter
2018 PW 51.0945 -1.1188 124 Peter Wakelin
2115 MW 51.3286 -0.7950 75 Mike Waterman
2414 DH 50.7468 -1.8800 34 David Hopkins
2420 RE 55.9486 -3.1386 40 Russell Eberst
2563 PN 51.0524 2.4043 10 Pierre Nierinck
2675 DB 52.1358 -2.3264 70 David Brierley
2701 TM 43.6876 -79.3924 230 Ted Molczan
2751 BM 51.3440 -1.9849 125 Bruce MacDonald
2756 AK 56.0907 -3.1623 25 Andy Kirkham
4353 ML 52.1541 4.4908 0 Marco Langbroek
4354 ML 52.1168 4.5602 -2 Marco Langbroek
4355 ML 52.1388 4.4994 -2 Marco Langbroek
4541 AR 41.9639 12.4531 80 Alberto Rango
4542 AR 41.9683 12.4545 80 Alberto Rango
4641 AR 41.1060 16.9010 70 Alberto Rango
5555 SG 56.1019 94.5533 154 Sergey Guryanov
5918 BG 59.2985 18.1045 40 Bjorn Gimle
5919 BG 59.2615 18.6206 30 Bjorn Gimle
6226 SC 28.4861 -97.8194 110 Scott Campbell
7777 BY 38.1656 -2.3267 1608 Brad Young remote
7778 BY -31.2733 149.0644 1122 Brad Young SSO
7779 BY 32.9204 -105.5283 2225 Brad Young NM
8048 ST 49.4175 -123.6420 1. Scott Tilley
8049 ST 49.4348 -123.6685 40. Scott Tilley
8305 PG 26.2431 -98.2163 30 Paul Gabriel
8335 BY 36.1397 -95.9838 205 Brad Young
8336 BY 36.1397 -95.9838 205 Brad Young
8438 SS 38.8108 -119.7750 1545 Sierra Stars
8536 TL 36.8479 -76.5010 4 Tim Luton
8539 SN 39.4707 -79.3388 839 Steve Newcomb
8597 TB -34.9638 138.6333 100 Tony Beresford
8600 PC -32.9770 151.6477 18 Paul Camilleri
0699 PC -14.4733 132.2369 108 Paul Camilleri
8730 EC 30.3086 -97.7279 150 Ed Cannon
8739 DB 37.1133 -121.7028 282 Derek Breit
9461 KO 47.9175 19.89365 938 Konkoly Obs
9633 JN 33.0206 -96.9982 153 Jim Nix
9730 MM 30.3150 -97.8660 280 Mike McCants
6242 JM 42.9453 -2.8284 623 Jon Mikel
6241 JM 42.9565 -2.8203 619 Jon Mikel
4160 BD 51.2793 5.4768 35 Bram Dorreman
9999 GR 47.348 5.5151 100 Graves

View File

@ -1,391 +0,0 @@
38860 2214.487
24753 2214.663
33446 2216.526
31490 2216.527
36110 2216.529
37165 2216.530
39363 2216.892
29107 2217.506
33456 2220.485
35951 2222.493
35951 2222.495
29522 2222.510
31598 2229.238
37216 2229.243
31598 2229.494
33412 2229.498
37216 2229.499
37179 2229.603
31598 2229.750
33412 2229.754
37216 2229.755
31598 2230.007
33412 2230.010
37216 2230.010
31598 2230.263
33412 2230.266
37216 2230.267
31598 2230.518
33412 2230.522
37216 2230.523
31598 2230.774
33412 2230.778
37216 2230.778
33433 2230.796
33412 2231.290
27640 2232.259
28537 2232.501
28537 2232.502
31701 2232.503
37386 2232.505
34505 2233.425
36423 2235.470
28414 2235.847
39209 2240.704
33320 2241.600
33320 2241.601
32378 2242.503
28384 2242.505
28888 2242.510
24680 2242.510
37348 2242.510
90004 2242.518
39088 2243.334
38782 2243.765
35937 2244.095
29522 2244.301
39363 2244.729
28649 2245.686
29710 2245.686
36795 2245.686
32783 2245.688
37849 2248.028
32783 2249.688
28051 2250.010
35951 2275.505
35938 2275.577
35938 2275.580
35938 2275.581
35938 2275.905
35938 2275.908
35938 2275.909
37874 2276.396
35938 2276.601
35938 2276.604
35938 2276.605
35938 2277.625
35938 2277.628
35938 2277.629
26619 2278.407
35938 2278.649
35938 2278.652
35938 2278.653
39209 2243.814
35937 2245.797
32378 2242.759
32378 2242.247
38466 2242.506
39209 2242.894
35937 2243.403
32378 2243.528
28384 2243.528
35951 2245.338
39209 2245.534
35937 2245.451
32378 2245.903
35951 2278.168
32958 2246.396
28051 2246.010
37849 2246.979
37214 2246.406
37849 2246.125
37849 2247.503
28051 2250.010
29710 2249.686
37387 2250.006
28649 2249.686
37849 2249.930
37849 2249.076
37849 2240.688
35937 2240.100
35937 2240.332
23533 2240.576
35951 2240.586
29522 2240.590
32378 2241.480
28384 2241.480
29522 2241.615
37849 2241.736
35937 2241.356
35937 2241.356
26092 2241.607
39209 2241.794
35937 2242.380
37849 2243.833
37930 2243.543
39209 2243.604
29522 2237.517
10033 2237.517
36834 2234.001
36834 2238.001
37180 2243.844
38046 2236.570
35951 2237.517
39150 2235.131
37180 2235.844
37849 2237.542
37849 2234.396
37849 2242.784
37849 2244.882
37849 2245.930
33314 2243.340
24753 2212.616
32289 2212.812
28054 2222.496
33408 2224.566
38038 2225.521
28051 2234.009
39260 2236.892
38038 2241.520
39262 2246.806
35937 2247.492
37849 2233.076
37849 2232.027
37387 2234.005
35937 2234.892
38861 2239.040
35931 2234.004
33433 2246.795
38997 2233.565
33434 2233.995
38773 2232.505
33409 2235.447
37216 2246.010
37216 2246.266
31598 2246.007
31598 2246.263
32376 2246.276
33412 2246.266
38248 2234.005
37849 2232.027
35951 2232.409
37820 2232.220
25991 2237.511
28054 2237.512
23533 2237.528
37180 2251.844
27392 2245.863
33320 2257.600
39260 2252.893
38861 2255.041
39630 2237.516
23533 2246.360
39630 2241.612
29522 2248.398
35951 2241.612
29522 2246.350
29522 2247.374
23533 2247.384
36119 2248.503
35951 2248.409
20963 2242.498
39751 2242.508
36121 2266.322
38012 2269.200
39019 2269.200
29108 2268.457
24753 2213.640
24753 2210.568
24753 2211.592
27391 2211.011
25635 2215.007
38354 2211.850
29505 2213.523
28479 2217.371
27640 2212.511
34388 2218.008
39260 2218.568
29499 2230.013
39150 2225.132
38771 2230.010
39262 2230.806
28537 2234.064
28537 2230.939
35951 2244.312
35951 2246.360
35938 2269.105
35938 2271.481
39678 2270.517
26464 2249.006
35937 2246.473
35937 2248.521
35937 2249.545
31113 2273.522
26934 2242.479
35937 2244.421
39209 2242.493
39232 2242.511
39232 2244.559
26619 2242.408
35931 2250.004
37849 2250.125
24753 2252.507
33434 2249.995
38248 2250.006
35951 2252.505
37849 2255.367
37849 2252.222
29479 2252.224
29479 2255.700
29522 2252.495
36834 2253.841
39634 2254.116
35937 2255.688
35937 2252.616
10033 2252.536
23533 2252.507
39630 2252.505
24753 2252.506
78407 2252.515
78408 2252.515
78409 2252.515
24753 2203.400
24753 2202.376
36508 2201.010
24753 2201.352
39380 2202.513
90020 2242.504
27391 2212.059
27640 2210.937
40118 2212.568
33321 2213.800
27640 2212.253
31140 2210.906
36599 2213.015
27391 2213.111
29107 2213.006
26958 2245.014
39630 2245.338
39630 2262.181
39630 2264.843
25789 2263.151
25789 2264.200
32958 2264.650
29522 2245.327
37849 2243.956
39630 2272.831
39630 2270.168
23533 2243.290
23533 2244.314
23533 2245.338
23533 2243.653
39630 2244.314
40137 2242.404
39086 2280.005
36121 2282.160
39765 2280.005
26619 2282.407
19822 2280.509
37387 2280.008
39012 2279.672
39497 2280.036
39497 2279.972
37941 2280.727
39497 2279.404
39497 2280.604
40301 2280.011
31304 2280.504
31698 2280.000
19460 2282.452
27640 2222.307
33412 2225.147
40013 2224.005
37820 2232.603
35683 2233.338
38709 2233.333
37216 2232.315
40344 2232.5
40344 2277.515
39765 2279.406
36037 2235.009
25544 143.050
35938 2279.671
39497 2276.036
26958 2223.014
36122 2274.479
35938 2272.505
32750 2272.945
29658 2272.945
33244 2272.945
32283 2272.945
31797 2272.945
32750 2273.072
29658 2273.072
33244 2273.072
32283 2273.072
31797 2273.072
35951 2272.842
19822 2279.984
39462 2242.505
35951 2253.531
35951 2254.555
05560 2242.534
39088 2243.334
26958 2243.013
40555 2280.006
35938 2281.725
35938 2280.701
39210 2280.006
27858 2232.711
37849 2235.445
39768 2235.006
37391 2234.069
24753 2235.450
33434 2245.995
38248 2246.006
33434 2245.995
33434 2245.995
29522 2249.423
39630 2247.384
35951 2249.436
35951 2247.388
39630 2246.361
39630 2248.409
39630 2249.433
40143 2248.555
40298 2271.604
25789 2264.989
25476 2264.857
40588 2215.993
23802 2264.978
27640 2215.876
27640 2216.900
32289 2216.817
27640 2215.364
40336 2213.126
29522 2262.196
29522 2264.859
27640 2221.508
40013 2220.005
39766 2220.239
38337 2224.003
35937 2237.258
28054 2235.463
24753 2237.498
29522 2239.570
35951 2239.567
24753 2239.567
39630 2239.567
29522 2238.545
35951 2238.545
24753 2238.545
39630 2238.545
29522 2236.497
35951 2236.497
24753 2236.497
39630 2236.497
28414 2239.846
35937 2238.283
35937 2239.307
37180 2239.844
39159 2238.851

View File

@ -1,133 +0,0 @@
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlinear.h>
struct data {
size_t n;
double *x,*y,*sy,*w;
};
int gauss_f(const gsl_vector *q,void *d,gsl_vector *f)
{
int i,n=((struct data *) d)->n;
double *x=((struct data *) d)->x;
double *y=((struct data *) d)->y;
double a[5],ym,arg,ex,fac;
// Extract parameters
for (i=0;i<5;i++)
a[i]=gsl_vector_get(q,i);
// Compute values
for (i=0;i<n;i++) {
arg=(x[i]-a[0])/a[1];
ex=exp(-arg*arg);
fac=2.0*a[2]*ex*arg;
ym=a[2]*ex+a[3]+a[4]*x[i];
gsl_vector_set(f,i,ym-y[i]);
}
return GSL_SUCCESS;
}
int gauss_df(const gsl_vector *q,void *d,gsl_matrix *J)
{
int i,n=((struct data *) d)->n;
double *x=((struct data *) d)->x;
double *y=((struct data *) d)->y;
double a[5],ym,arg,ex,fac;
// Extract parameters
for (i=0;i<5;i++)
a[i]=gsl_vector_get(q,i);
// Compute values
for (i=0;i<n;i++) {
arg=(x[i]-a[0])/a[1];
ex=exp(-arg*arg);
fac=2.0*a[2]*ex*arg;
gsl_matrix_set(J,i,0,fac/a[1]);
gsl_matrix_set(J,i,1,fac*arg/a[1]);
gsl_matrix_set(J,i,2,ex);
gsl_matrix_set(J,i,3,1.0);
gsl_matrix_set(J,i,4,x[i]);
}
return GSL_SUCCESS;
}
int fit_gaussian(double *x,double *y,double *sy,int n,double *q,double *sq,int m,double *chisq)
{
int i;
const gsl_multifit_nlinear_type *T=gsl_multifit_nlinear_trust;
gsl_multifit_nlinear_workspace *w;
gsl_multifit_nlinear_fdf fdf;
gsl_multifit_nlinear_parameters fdf_params=gsl_multifit_nlinear_default_parameters();
gsl_vector *f;
gsl_matrix *J;
struct data d;
gsl_rng *r;
int status,info;
const double xtol=1e-8;
const double gtol=1e-8;
const double ftol=0.0;
// Fill struct
d.n=n;
d.x=(double *) malloc(sizeof(double)*d.n);
d.y=(double *) malloc(sizeof(double)*d.n);
d.sy=(double *) malloc(sizeof(double)*d.n);
d.w=(double *) malloc(sizeof(double)*d.n);
for (i=0;i<d.n;i++) {
d.x[i]=x[i];
d.y[i]=y[i];
d.sy[i]=sy[i];
d.w[i]=1.0/(d.sy[i]*d.sy[i]);
}
gsl_vector_view a=gsl_vector_view_array(q,m);
gsl_vector_view wts=gsl_vector_view_array(d.w,d.n);
gsl_matrix *covar=gsl_matrix_alloc(m,m);
// Function definition
fdf.f=gauss_f;
fdf.df=gauss_df;
fdf.fvv=NULL;
fdf.n=d.n;
fdf.p=m;
fdf.params=&d;
// Allocate and initialize
w=gsl_multifit_nlinear_alloc(T,&fdf_params,d.n,m);
gsl_multifit_nlinear_winit(&a.vector,&wts.vector,&fdf,w);
f=gsl_multifit_nlinear_residual(w);
// Solve and compute covariance and chi-squared
status=gsl_multifit_nlinear_driver(20,xtol,gtol,ftol,NULL,NULL,&info,w);
J=gsl_multifit_nlinear_jac(w);
gsl_multifit_nlinear_covar(J,0.0,covar);
gsl_blas_ddot(f,f,chisq);
// Extract parameters
for (i=0;i<m;i++) {
q[i]=gsl_vector_get(w->x,i);
sq[i]=sqrt(gsl_matrix_get(covar,i,i));
}
// Free
gsl_multifit_nlinear_free (w);
gsl_matrix_free (covar);
gsl_rng_free (r);
free(d.x);
free(d.y);
free(d.sy);
free(d.w);
return status;
}

View File

@ -1,225 +0,0 @@
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "sgdp4h.h"
#include "satutl.h"
#include "cpgplot.h"
#define LIM 80
#define D2R M_PI/180.0
#define R2D 180.0/M_PI
#define XKMPER 6378.135 // Earth radius in km
#define XKMPAU 149597879.691 // AU in km
#define FLAT (1.0/298.257)
#define C 299792.458 // Speed of light in km/s
struct point {
xyz_t obspos,obsvel;
};
struct site {
int id;
double lng,lat;
float alt;
char observer[64];
};
// Return x modulo y [0,y)
double modulo(double x,double y)
{
x=fmod(x,y);
if (x<0.0) x+=y;
return x;
}
// 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 == '\n')
// s[i++] = c;
s[i] = '\0';
return i;
}
// Greenwich Mean Sidereal Time
double gmst(double mjd)
{
double t,gmst;
t=(mjd-51544.5)/36525.0;
gmst=modulo(280.46061837+360.98564736629*(mjd-51544.5)+t*t*(0.000387933-t/38710000),360.0);
return gmst;
}
// Greenwich Mean Sidereal Time
double dgmst(double mjd)
{
double t,dgmst;
t=(mjd-51544.5)/36525.0;
dgmst=360.98564736629+t*(0.000387933-t/38710000);
return dgmst;
}
// Observer position
void obspos_xyz(double mjd,double lng,double lat,float alt,xyz_t *pos,xyz_t *vel)
{
double ff,gc,gs,theta,s,dtheta;
s=sin(lat*D2R);
ff=sqrt(1.0-FLAT*(2.0-FLAT)*s*s);
gc=1.0/ff+alt/XKMPER;
gs=(1.0-FLAT)*(1.0-FLAT)/ff+alt/XKMPER;
theta=gmst(mjd)+lng;
dtheta=dgmst(mjd)*D2R/86400;
pos->x=gc*cos(lat*D2R)*cos(theta*D2R)*XKMPER;
pos->y=gc*cos(lat*D2R)*sin(theta*D2R)*XKMPER;
pos->z=gs*sin(lat*D2R)*XKMPER;
vel->x=-gc*cos(lat*D2R)*sin(theta*D2R)*XKMPER*dtheta;
vel->y=gc*cos(lat*D2R)*cos(theta*D2R)*XKMPER*dtheta;
vel->z=0.0;
return;
}
// Get observing site
struct site get_site(int site_id)
{
int i=0;
char line[LIM];
FILE *file;
int id;
double lat,lng;
float alt;
char abbrev[3],observer[64];
struct site s;
char *env,filename[LIM];
env=getenv("ST_DATADIR");
sprintf(filename,"%s/data/sites.txt",env);
file=fopen(filename,"r");
if (file==NULL) {
printf("File with site information not found!\n");
return;
}
while (fgets(line,LIM,file)!=NULL) {
// Skip
if (strstr(line,"#")!=NULL)
continue;
// Strip newline
line[strlen(line)-1]='\0';
// Read data
sscanf(line,"%4d %2s %lf %lf %f",
&id,abbrev,&lat,&lng,&alt);
strcpy(observer,line+38);
// Change to km
alt/=1000.0;
// Copy site
if (id==site_id) {
s.lat=lat;
s.lng=lng;
s.alt=alt;
s.id=id;
strcpy(s.observer,observer);
}
}
fclose(file);
return s;
}
// Plot overlay
void overlay(double *mjd,int n,int site_id)
{
int i,imode,flag,satno,tflag;
struct point *p;
struct site s;
FILE *file,*infile;
orbit_t orb;
xyz_t satpos,satvel;
double dx,dy,dz,dvx,dvy,dvz,r,v,za;
double freq,freq0;
char line[LIM],text[8];
// Get site
s=get_site(site_id);
// Allocate
p=(struct point *) malloc(sizeof(struct point)*n);
// Get observer position
for (i=0;i<n;i++)
obspos_xyz(mjd[i],s.lng,s.lat,s.alt,&p[i].obspos,&p[i].obsvel);
infile=fopen("frequencies.txt","r");
while (fgetline(infile,line,LIM)>0) {
sscanf(line,"%d %lf",&satno,&freq0);
sprintf(text," %d",satno);
// Loop over TLEs
file=fopen("/home/bassa/code/c/satellite/sattools/tle/bulk.tle","r");
while (read_twoline(file,satno,&orb)==0) {
// Initialize
imode=init_sgdp4(&orb);
if (imode==SGDP4_ERROR)
printf("Error\n");
// Loop over points
for (i=0,flag=0,tflag=0;i<n;i++) {
// Get satellite position
satpos_xyz(mjd[i]+2400000.5,&satpos,&satvel);
dx=satpos.x-p[i].obspos.x;
dy=satpos.y-p[i].obspos.y;
dz=satpos.z-p[i].obspos.z;
dvx=satvel.x-p[i].obsvel.x;
dvy=satvel.y-p[i].obsvel.y;
dvz=satvel.z-p[i].obsvel.z;
r=sqrt(dx*dx+dy*dy+dz*dz);
v=(dvx*dx+dvy*dy+dvz*dz)/r;
za=acos((p[i].obspos.x*dx+p[i].obspos.y*dy+p[i].obspos.z*dz)/(r*XKMPER))*R2D;
freq=(1.0-v/C)*freq0;
if (flag==0) {
cpgmove((float) i,(float) freq);
flag=1;
} else {
cpgdraw((float) i,(float) freq);
}
if (za<90.0 && flag==0) {
flag=1;
} else if (za>90.0 && flag==1) {
if (tflag==0) {
tflag=1;
cpgtext((float) i,(float) freq,text);
}
flag=0;
}
}
}
fclose(file);
}
fclose(infile);
// Free
free(p);
return;
}

227
rfsim.c
View File

@ -1,227 +0,0 @@
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "sgdp4h.h"
#include "satutl.h"
#include "rftime.h"
#define LIM 80
#define D2R M_PI/180.0
#define R2D 180.0/M_PI
#define XKMPER 6378.135 // Earth radius in km
#define XKMPAU 149597879.691 // AU in km
#define FLAT (1.0/298.257)
#define C 299792.458 // Speed of light in km/s
struct point {
xyz_t obspos,obsvel;
xyz_t grvpos,grvvel;
};
struct site {
int id;
double lng,lat;
float alt;
char observer[64];
};
struct site get_site(int site_id);
void obspos_xyz(double mjd,double lng,double lat,float alt,xyz_t *pos,xyz_t *vel);
double gmst(double mjd);
double dgmst(double mjd);
double modulo(double x,double y);
void equatorial2horizontal(double mjd,struct site s,double ra,double de,double *azi,double *alt);
int main(int argc,char *argv[])
{
int i,j,n=86400;
int imode;
double *mjd,mjd0=57028.0;
struct point *p;
struct site s,g;
FILE *file;
orbit_t orb;
double dx,dy,dz,dvx,dvy,dvz,za;
double r0,v0,ra0,de0,azi0,alt0;
xyz_t satpos,satvel;
double freq,freq0=2272.5e6;
char nfd[32];
double t=0.0,a,f;
// Arrays
mjd=(double *) malloc(sizeof(double)*n);
p=(struct point *) malloc(sizeof(struct point)*n);
// Get sites
s=get_site(4171);
g=get_site(9999);
// MJD range
for (i=0;i<n;i++)
mjd[i]=mjd0+(double) i/86400.0;
// Get positions
for (i=0;i<n;i++) {
obspos_xyz(mjd[i],s.lng,s.lat,s.alt,&p[i].obspos,&p[i].obsvel);
obspos_xyz(mjd[i],g.lng,g.lat,g.alt,&p[i].grvpos,&p[i].grvvel);
}
// Loop over objects
file=fopen("/home/bassa/code/c/satellite/sattools/tle/catalog.tle","r");
while (read_twoline(file,25544,&orb)==0) {
// Initialize
imode=init_sgdp4(&orb);
if (imode==SGDP4_ERROR)
printf("Error\n");
// Loop over points
for (i=0;i<n;i++) {
// Get satellite position
satpos_xyz(mjd[i]+2400000.5,&satpos,&satvel);
// Observer position
dx=satpos.x-p[i].obspos.x;
dy=satpos.y-p[i].obspos.y;
dz=satpos.z-p[i].obspos.z;
dvx=satvel.x-p[i].obsvel.x;
dvy=satvel.y-p[i].obsvel.y;
dvz=satvel.z-p[i].obsvel.z;
r0=sqrt(dx*dx+dy*dy+dz*dz);
v0=(dvx*dx+dvy*dy+dvz*dz)/r0;
ra0=modulo(atan2(dy,dx)*R2D,360.0);
de0=asin(dz/r0)*R2D;
za=acos((p[i].obspos.x*dx+p[i].obspos.y*dy+p[i].obspos.z*dz)/(r0*XKMPER))*R2D;
if (za<90.0)
printf("%14.8lf %f %f %10.0lf %f\n",mjd[i],r0,v0,(1.0-v0/C)*freq0,za);
}
}
fclose(file);
// Free
free(mjd);
free(p);
return 0;
}
// Get observing site
struct site get_site(int site_id)
{
int i=0,status;
char line[LIM];
FILE *file;
int id;
double lat,lng;
float alt;
char abbrev[3],observer[64];
struct site s;
char *env,filename[LIM];
env=getenv("ST_DATADIR");
sprintf(filename,"%s/data/sites.txt",env);
file=fopen(filename,"r");
if (file==NULL) {
printf("File with site information not found!\n");
return;
}
while (fgets(line,LIM,file)!=NULL) {
// Skip
if (strstr(line,"#")!=NULL)
continue;
// Strip newline
line[strlen(line)-1]='\0';
// Read data
status=sscanf(line,"%4d %2s %lf %lf %f",
&id,abbrev,&lat,&lng,&alt);
strcpy(observer,line+38);
// Change to km
alt/=1000.0;
// Copy site
if (id==site_id) {
s.lat=lat;
s.lng=lng;
s.alt=alt;
s.id=id;
strcpy(s.observer,observer);
}
}
fclose(file);
return s;
}
// Observer position
void obspos_xyz(double mjd,double lng,double lat,float alt,xyz_t *pos,xyz_t *vel)
{
double ff,gc,gs,theta,s,dtheta;
s=sin(lat*D2R);
ff=sqrt(1.0-FLAT*(2.0-FLAT)*s*s);
gc=1.0/ff+alt/XKMPER;
gs=(1.0-FLAT)*(1.0-FLAT)/ff+alt/XKMPER;
theta=gmst(mjd)+lng;
dtheta=dgmst(mjd)*D2R/86400;
pos->x=gc*cos(lat*D2R)*cos(theta*D2R)*XKMPER;
pos->y=gc*cos(lat*D2R)*sin(theta*D2R)*XKMPER;
pos->z=gs*sin(lat*D2R)*XKMPER;
vel->x=-gc*cos(lat*D2R)*sin(theta*D2R)*XKMPER*dtheta;
vel->y=gc*cos(lat*D2R)*cos(theta*D2R)*XKMPER*dtheta;
vel->z=0.0;
return;
}
// Greenwich Mean Sidereal Time
double gmst(double mjd)
{
double t,gmst;
t=(mjd-51544.5)/36525.0;
gmst=modulo(280.46061837+360.98564736629*(mjd-51544.5)+t*t*(0.000387933-t/38710000),360.0);
return gmst;
}
// Greenwich Mean Sidereal Time
double dgmst(double mjd)
{
double t,dgmst;
t=(mjd-51544.5)/36525.0;
dgmst=360.98564736629+t*(0.000387933-t/38710000);
return dgmst;
}
// Return x modulo y [0,y)
double modulo(double x,double y)
{
x=fmod(x,y);
if (x<0.0) x+=y;
return x;
}
// Convert equatorial into horizontal coordinates
void equatorial2horizontal(double mjd,struct site s,double ra,double de,double *azi,double *alt)
{
double h;
h=gmst(mjd)+s.lng-ra;
*azi=modulo(atan2(sin(h*D2R),cos(h*D2R)*sin(s.lat*D2R)-tan(de*D2R)*cos(s.lat*D2R))*R2D,360.0);
*alt=asin(sin(s.lat*D2R)*sin(de*D2R)+cos(s.lat*D2R)*cos(de*D2R)*cos(h*D2R))*R2D;
return;
}

View File

@ -469,17 +469,21 @@ struct trace *compute_trace(char *tlefile,double *mjd,int n,int site_id,float fr
// Find number of satellites in frequency range
infile=fopen(freqlist,"r");
for (i=0;;) {
if (fgetline(infile,line,LIM)<=0)
break;
status=sscanf(line,"%d %lf",&satno,&freq0);
if (freq0>=fmin && freq0<=fmax)
i++;
if (infile==NULL) {
printf("%s not found\n",freqlist);
return t;
} else {
for (i=0;;) {
if (fgetline(infile,line,LIM)<=0)
break;
status=sscanf(line,"%d %lf",&satno,&freq0);
if (freq0>=fmin && freq0<=fmax)
i++;
}
fclose(infile);
*nsat=i;
}
fclose(infile);
*nsat=i;
// Break out
if (i==0)
return t;