Major update

pull/10/head
Cees Bassa 2015-05-25 23:15:20 +02:00
parent 05703c2f03
commit d10d578f5b
9 changed files with 118 additions and 31 deletions

View File

@ -80,7 +80,7 @@ int main(int argc,char *argv[])
}
// Read data
s=read_spectrogram(path,isub,nsub,f0,df0,nbin);
s=read_spectrogram(path,isub,nsub,f0,df0,nbin,0.0);
// Write data
write_spectrogram(s,outfile);

View File

@ -172,7 +172,7 @@ int main(int argc,char *argv[])
}
// Read data
s=read_spectrogram(path,isub,nsub,f0,df0,1);
s=read_spectrogram(path,isub,nsub,f0,df0,1,0.0);
// Exit on emtpy file
if (s.nsub==0)

15
rffit.c
View File

@ -201,7 +201,7 @@ int identify_satellite(char *catalog,double rmsmax)
flag=1;
}
if (rms<rmsmax) {
printf("%05d %.3f kHz\n",orb.satno,rms);
printf("%05d %.3f kHz %.6f MHz\n",orb.satno,rms,d.ffit/1000.0);
i++;
}
}
@ -319,7 +319,8 @@ int main(int argc,char *argv[])
fclose(fp);
}
freopen("/tmp/stderr.txt","w",stderr);
if (freopen("/tmp/stderr.txt","w",stderr)==NULL)
fprintf(stderr,"Failed to redirect stderr\n");
cpgopen("/xs");
cpgask(0);
@ -566,7 +567,8 @@ int main(int argc,char *argv[])
status=scanf("%i",&i);
if (i>=0 && i<=9) {
printf("\nNew value: ");
fgets(string,64,stdin);
if (fgets(string,64,stdin)==NULL)
fprintf(stderr,"Failed to read string\n");
status=scanf("%s",string);
// if (i==0) strcpy(d.satname,string);
if (i==1) orb.eqinc=RAD(atof(string));
@ -634,7 +636,7 @@ int main(int argc,char *argv[])
printf("No points selected!\n");
} else {
rms=fit_curve(orb,ia);
printf("rms: %.3f kHz, cf: %.3f MHz, TCA: %s\n",rms,d.ffit/1000.0,nfd);
printf("rms: %.3f kHz, cf: %.6f MHz, TCA: %s\n",rms,d.ffit/1000.0,nfd);
redraw=1;
}
}
@ -1282,7 +1284,10 @@ void save_data(float xmin,float ymin,float xmax,float ymax,char *filename)
for (i=0,j=0;i<d.n;i++) {
if (d.p[i].t>xmin && d.p[i].t<xmax && d.p[i].f>ymin && d.p[i].f<ymax && d.p[i].flag==2) {
// fprintf(file,"%s\t%14.3lf\t%8.3f\t%04d\n",d.p[i].timestamp,1000.0*d.p[i].freq,d.p[i].flux,d.p[i].site_id);
fprintf(file,"%lf\t%14.3lf\t%8.3f\t%04d\n",d.p[i].mjd,1000.0*d.p[i].freq,d.p[i].flux,d.p[i].site_id);
if (d.p[i].rsite_id==0)
fprintf(file,"%lf\t%14.3lf\t%8.3f\t%04d\n",d.p[i].mjd,1000.0*d.p[i].freq,d.p[i].flux,d.p[i].site_id);
else
fprintf(file,"%lf\t%14.3lf\t%8.3f\t%04d\t%04d\n",d.p[i].mjd,1000.0*d.p[i].freq,d.p[i].flux,d.p[i].site_id,d.p[i].rsite_id);
j++;
}
}

4
rfio.c
View File

@ -4,7 +4,7 @@
#include "rftime.h"
#include "rfio.h"
struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,double df0,int nbin)
struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,double df0,int nbin,double foff)
{
int i,j,k,l,flag=0,status,msub,ibin,nadd;
char filename[128],header[256],nfd[32];
@ -32,6 +32,8 @@ struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,dou
status=fread(header,sizeof(char),256,file);
status=sscanf(header,"HEADER\nUTC_START %s\nFREQ %lf Hz\nBW %lf Hz\nLENGTH %f s\nNCHAN %d\n",s.nfd0,&s.freq,&s.samp_rate,&length,&nch);
s.freq+=foff;
// Close file
fclose(file);

2
rfio.h
View File

@ -6,5 +6,5 @@ struct spectrogram {
float *z;
char nfd0[32];
};
struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,double df0,int nbin);
struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,double df0,int nbin,double foff);
void write_spectrogram(struct spectrogram s,char *prefix);

View File

@ -20,7 +20,7 @@ void dec2sex(double x,char *s,int f,int len);
void time_axis(double *mjd,int n,float xmin,float xmax,float ymin,float ymax);
void usage(void);
void plot_traces(struct trace *t,int nsat);
struct trace fit_trace(struct spectrogram s,struct select sel,int site_id);
struct trace fit_trace(struct spectrogram s,struct select sel,int site_id,int graves);
void convolve(float *y,int n,float *w,int m,float *z);
float gauss(float x,float w);
void quadfit(float x[],float y[],int n,float a[]);
@ -189,7 +189,7 @@ void filter(struct spectrogram s,int site_id)
void peakfind(struct spectrogram s,int site_id,int i0,int i1,int j0,int j1)
{
int i,j,k,l,m=21,n;
float *w,*y,*sy,*a,*b,*c,d[3],dx[3],dw=2.0,x0,c0=-0.005;
float *w,*y,*sy,*a,*b,*c,d[3],dx[3],dw=1.0,x0,c0=-0.000008;
double f;
FILE *file;
@ -273,7 +273,7 @@ int main(int argc,char *argv[])
float heat_b[] = {0.0, 0.0, 0.0, 0.3, 1.0};
float xmin,xmax,ymin,ymax,zmin,zmax=8.0;
int i,j,k,flag=0,isel=0,sn;
int redraw=1,mode=0,posn=0,click=0,graves=0;
int redraw=1,mode=0,posn=0,click=0,graves=0,grid=0;
float dt,zzmax,s1,s2,z,za,sigma,zs,zm;
int ix=0,iy=0,isub=0;
int i0,j0,i1,j1,jmax;
@ -294,6 +294,8 @@ int main(int argc,char *argv[])
char *env;
int site_id=0;
int cmap=0;
double foff=0.0,mjdgrid=0.0;
int jj0,jj1;
// Get site
env=getenv("ST_COSPAR");
@ -307,7 +309,7 @@ int main(int argc,char *argv[])
// Read arguments
if (argc>1) {
while ((arg=getopt(argc,argv,"p:f:w:s:l:b:z:hc:C:g"))!=-1) {
while ((arg=getopt(argc,argv,"p:f:w:s:l:b:z:hc:C:gm:o:"))!=-1) {
switch (arg) {
case 'p':
@ -346,6 +348,16 @@ int main(int argc,char *argv[])
usage();
return 0;
case 'm':
cmap=atoi(optarg);
if (cmap>2)
cmap=0;
break;
case 'o':
foff=(double) atof(optarg);
break;
case 'c':
strcpy(tlefile,optarg);
break;
@ -365,8 +377,8 @@ int main(int argc,char *argv[])
}
// Read data
s=read_spectrogram(path,isub,nsub,f0,df0,nbin);
s=read_spectrogram(path,isub,nsub,f0,df0,nbin,foff);
printf("Read spectrogram\n%d channels, %d subints\nFrequency: %g MHz\nBandwidth: %g MHz\n",s.nchan,s.nsub,s.freq*1e-6,s.samp_rate*1e-6);
// Exit on empty data
@ -474,6 +486,31 @@ int main(int argc,char *argv[])
cpgsci(1);
}
// Plot grid
if (grid==1) {
cpgsci(2);
for (i=0,flag=0;i<s.nsub-1;i++) {
dt=86400.0*(s.mjd[i]-mjdgrid);
jj1=(int) (floor) (dt/2.3);
if (i==0)
jj0=jj1;
if (jj1-jj0>0.0) {
flag=0;
jj0=jj1;
}
if (jj0%2==0)
cpgsls(1);
else
cpgsls(2);
if (flag==0) {
cpgmove((float) i,ymin);
cpgdraw((float) i,ymax);
flag=1;
}
}
cpgsci(1);
cpgsls(1);
}
redraw=0;
}
@ -484,6 +521,16 @@ int main(int argc,char *argv[])
if (c=='q')
break;
// Toggle grid
if (c=='k') {
if (grid==0)
grid=1;
else
grid=0;
mjdgrid=s.mjd[(int) floor(x)];
redraw=1;
}
// Track
if (c=='t') {
for (i=0;i<nsat;i++) {
@ -523,7 +570,7 @@ int main(int argc,char *argv[])
// Fit
if (c=='f') {
tf=fit_trace(s,sel,site_id);
tf=fit_trace(s,sel,site_id,graves);
tf.site=site_id;
continue;
}
@ -627,7 +674,10 @@ int main(int argc,char *argv[])
j=(int) floor(y);
f=s.freq-0.5*s.samp_rate+(double) j*s.samp_rate/(double) s.nchan;
if (s.mjd[i]>1.0) {
fprintf(file,"%lf %lf %f %d\n",s.mjd[i],f,s.z[i+s.nsub*j],site_id);
if (graves==0)
fprintf(file,"%lf %lf %f %d\n",s.mjd[i],f,s.z[i+s.nsub*j],site_id);
else
fprintf(file,"%lf %lf %f %d 9999\n",s.mjd[i],f,s.z[i+s.nsub*j],site_id);
printf("%lf %lf %f %d\n",s.mjd[i],f,s.z[i+s.nsub*j],site_id);
}
fclose(file);
@ -672,7 +722,10 @@ int main(int argc,char *argv[])
f=s.freq-0.5*s.samp_rate+(double) jmax*s.samp_rate/(double) s.nchan;
if (sigma>5.0 && s.mjd[i]>1.0) {
fprintf(file,"%lf %lf %f %d\n",s.mjd[i],f,zzmax,site_id);
if (graves==0)
fprintf(file,"%lf %lf %f %d\n",s.mjd[i],f,zzmax,site_id);
else
fprintf(file,"%lf %lf %f %d 9999\n",s.mjd[i],f,zzmax,site_id);
cpgpt1((float) i,(float) jmax,17);
}
}
@ -1029,7 +1082,7 @@ void plot_traces(struct trace *t,int nsat)
}
// Fit trace
struct trace fit_trace(struct spectrogram s,struct select sel,int site_id)
struct trace fit_trace(struct spectrogram s,struct select sel,int site_id,int graves)
{
int i,j,k,l,sn;
int i0,i1,j0,j1,jmax;
@ -1085,7 +1138,10 @@ struct trace fit_trace(struct spectrogram s,struct select sel,int site_id)
// Store
if (sigma>5.0 && s.mjd[i]>1.0) {
f=s.freq-0.5*s.samp_rate+(double) jmax*s.samp_rate/(double) s.nchan;
fprintf(file,"%lf %lf %f %d\n",s.mjd[i],f,sigma,site_id);
if (graves==0)
fprintf(file,"%lf %lf %f %d\n",s.mjd[i],f,sigma,site_id);
else
fprintf(file,"%lf %lf %f %d 9999\n",s.mjd[i],f,sigma,site_id);
cpgpt1((float) i,(float) jmax,17);
t.mjd[l]=s.mjd[i];
t.freq[l]=f;

View File

@ -103,7 +103,7 @@ int main(int argc,char *argv[])
for (isub=0;;isub+=15) {
// Read data
s=read_spectrogram(path,isub,nsub,f0,df0,nbin);
s=read_spectrogram(path,isub,nsub,f0,df0,nbin,0.0);
if (s.mjd[0]<54000)
break;
@ -113,7 +113,7 @@ int main(int argc,char *argv[])
printf("Read spectrogram\n%d channels, %d subints\nFrequency: %g MHz\nBandwidth: %g MHz\n",s.nchan,s.nsub,s.freq*1e-6,s.samp_rate*1e-6);
// Compute traces
t=compute_trace(tlefile,s.mjd,s.nsub,site_id,s.freq*1e-6,s.samp_rate*1e-6,&nsat);
t=compute_trace(tlefile,s.mjd,s.nsub,site_id,s.freq*1e-6,s.samp_rate*1e-6,&nsat,0);
printf("Traces for %d objects for location %d\n",nsat,site_id);
printf("%s\n",pngfile);

10
rfsim.c
View File

@ -16,6 +16,7 @@
struct point {
xyz_t obspos,obsvel;
xyz_t grvpos,grvvel;
};
struct site {
int id;
@ -36,7 +37,7 @@ int main(int argc,char *argv[])
int imode;
double *mjd,mjd0=57028.0;
struct point *p;
struct site s;
struct site s,g;
FILE *file;
orbit_t orb;
double dx,dy,dz,dvx,dvy,dvz,za;
@ -52,18 +53,21 @@ int main(int argc,char *argv[])
// 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++)
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,19822,&orb)==0) {
while (read_twoline(file,25544,&orb)==0) {
// Initialize
imode=init_sgdp4(&orb);
if (imode==SGDP4_ERROR)

View File

@ -94,6 +94,19 @@ void obspos_xyz(double mjd,double lng,double lat,float alt,xyz_t *pos,xyz_t *vel
return;
}
// Convert equatorial into horizontal coordinates
void equatorial2horizontal(double mjd,double ra,double de,double lng,double lat,double *azi,double *alt)
{
double h;
h=gmst(mjd)+lng-ra;
*azi=modulo(atan2(sin(h*D2R),cos(h*D2R)*sin(lat*D2R)-tan(de*D2R)*cos(lat*D2R))*R2D,360.0);
*alt=asin(sin(lat*D2R)*sin(de*D2R)+cos(lat*D2R)*cos(de*D2R)*cos(h*D2R))*R2D;
return;
}
// Get observing site
struct site get_site(int site_id)
{
@ -167,7 +180,8 @@ void identify_trace(char *tlefile,struct trace t,int satno)
sprintf(freqlist,"%s/data/frequencies.txt",env);
// Reloop stderr
freopen("/tmp/stderr.txt","w",stderr);
if (freopen("/tmp/stderr.txt","w",stderr)==NULL)
fprintf(stderr,"Failed to redirect stderr\n");
// Get site
s=get_site(t.site);
@ -208,6 +222,7 @@ void identify_trace(char *tlefile,struct trace t,int satno)
r=sqrt(dx*dx+dy*dy+dz*dz);
v[i]=(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;
beta=(1.0-v[i]/C);
sum1+=beta*t.freq[i];
sum2+=beta*beta;
@ -283,18 +298,18 @@ struct trace *compute_trace(char *tlefile,double *mjd,int n,int site_id,float fr
struct trace *t;
float fmin,fmax;
char *env,freqlist[LIM];
double ra,de,azi,alt;
env=getenv("ST_DATADIR");
sprintf(freqlist,"%s/data/frequencies.txt",env);
printf("%d\n",graves);
// Frequency limits
fmin=freq-0.5*bw;
fmax=freq+0.5*bw;
// Reloop stderr
freopen("/tmp/stderr.txt","w",stderr);
if (freopen("/tmp/stderr.txt","w",stderr)==NULL)
fprintf(stderr,"Failed to redirect stderr\n");
// Find number of satellites in frequency range
infile=fopen(freqlist,"r");
@ -395,9 +410,14 @@ struct trace *compute_trace(char *tlefile,double *mjd,int n,int site_id,float fr
dvz=satvel.z-p[i].grvel.z;
r=sqrt(dx*dx+dy*dy+dz*dz);
vg=(dvx*dx+dvy*dy+dvz*dz)/r;
// Graves frequency
ra=modulo(atan2(dy,dx)*R2D,360.0);
de=asin(dz/r)*R2D;
equatorial2horizontal(mjd[i],ra,de,sg.lng,sg.lat,&azi,&alt);
t[j].freq[i]=(1.0-v/C)*(1.0-vg/C)*freq0;
if (!((azi<90.0 || azi>270.0) && alt>15.0 && alt<40.0))
t[j].za[i]=100.0;
}
}
}