diff --git a/.gitignore b/.gitignore index 7afe4d2..3ca5119 100644 --- a/.gitignore +++ b/.gitignore @@ -16,10 +16,12 @@ *.exe *.out *.app + /rfdop /rfedit /rffft /rffind /rffit +/rfinfo /rfplot /rfpng diff --git a/Makefile.osx b/Makefile.osx index bfbb04d..9c38776 100644 --- a/Makefile.osx +++ b/Makefile.osx @@ -19,7 +19,8 @@ LFLAGS = -L$(prefix)/lib -lcpgplot -lpgplot -lX11 -lpng -lm -lgsl -lgslcblas # NOTE: STRF will not compile or link correctly with the system gcc (which is actually clang) # It's best to build with gcc provided by Macports or Homebrew # Under Macports, this is provided as gcc-mp-7, as below, on Homebrew this may be different. -CC = $(prefix)/bin/gcc-mp-7 +# Newer Macports releases use gcc-mp-10, therefor changed the default to gcc-mp-10 +CC = $(prefix)/bin/gcc-mp-10 # Installation INSTALL_PROGRAM = install -m 557 @@ -44,8 +45,8 @@ rffind: rffind.o rfio.o rftime.o rftrack: rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o $(CC) -o rftrack rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o -lm -rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o - $(CC) -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o $(LFLAGS) +rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o + $(CC) -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o $(LFLAGS) rffft: rffft.o rftime.o $(CC) -o rffft rffft.o rftime.o -lfftw3f -lm $(LFLAGS) diff --git a/makefile b/makefile index ba7ba82..62bb0af 100644 --- a/makefile +++ b/makefile @@ -34,8 +34,8 @@ rffind: rffind.o rfio.o rftime.o rftrack: rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o $(CC) -o rftrack rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o -lm -rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o - $(CC) -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o $(LFLAGS) +rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o + $(CC) -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o $(LFLAGS) rffft: rffft.o rftime.o $(CC) -o rffft rffft.o rftime.o -lfftw3f -lm diff --git a/rffft.c b/rffft.c index 37b4fe8..d2127f4 100644 --- a/rffft.c +++ b/rffft.c @@ -13,6 +13,7 @@ void usage(void) printf("rffft: FFT RF observations\n\n"); printf("-i Input file (can be fifo) [stdin]\n"); printf("-p Output prefix\n"); + printf("-o Output filename [default: YYYY-MM-DDTHH:MM:SS.sss_XXXXXX.bin]\n"); printf("-f Center frequency (Hz)\n"); printf("-s Sample rate (Hz)\n"); printf("-c Channel size [100Hz]\n"); @@ -22,6 +23,7 @@ void usage(void) printf("-F Input format char, int, float [int]\n"); printf("-T YYYY-MM-DDTHH:MM:SSS.sss\n"); printf("-R Frequency range to store (Hz)\n"); + printf("-I Invert frequencies\n"); printf("-b Digitize output to bytes [off]\n"); printf("-q Quiet mode, no output [off]\n"); printf("-h This help\n"); @@ -31,11 +33,11 @@ void usage(void) int main(int argc,char *argv[]) { - int i,j,k,l,m,nchan,nint=1,arg=0,nbytes,nsub=60,flag,nuse=1,realtime=1,quiet=0,imin,imax,partial=0; + int i,j,k,l,m,nchan,nint=1,arg=0,nbytes,nsub=60,flag,nuse=1,realtime=1,quiet=0,imin,imax,partial=0,useoutput=0; fftwf_complex *c,*d; fftwf_plan fft; FILE *infile,*outfile; - char infname[128]="",outfname[128]="",path[64]=".",prefix[32]=""; + char infname[128]="",outfname[128]="",path[64]=".",prefix[32]="",output[128]=""; char informat='i',outformat='f'; int16_t *ibuf; char *cbuf; @@ -45,10 +47,11 @@ int main(int argc,char *argv[]) double freq,samp_rate,mjd,freqmin=-1,freqmax=-1; struct timeval start,end; char tbuf[30],nfd[32],header[256]=""; + int sign=1; // Read arguments if (argc>1) { - while ((arg=getopt(argc,argv,"i:f:s:c:t:p:n:hm:F:T:bqR:"))!=-1) { + while ((arg=getopt(argc,argv,"i:f:s:c:t:p:n:hm:F:T:bqR:o:I"))!=-1) { switch(arg) { case 'i': @@ -59,6 +62,11 @@ int main(int argc,char *argv[]) strcpy(path,optarg); break; + case 'o': + strcpy(output,optarg); + useoutput=1; + break; + case 'f': freq=(double) atof(optarg); break; @@ -109,6 +117,10 @@ int main(int argc,char *argv[]) realtime=0; break; + case 'I': + sign=-1; + break; + case 'h': usage(); return 0; @@ -191,7 +203,11 @@ int main(int argc,char *argv[]) // Forever loop for (m=0;;m++) { // File name - sprintf(outfname,"%s/%s_%06d.bin",path,prefix,m); + if (useoutput==0) { + sprintf(outfname,"%s/%s_%06d.bin",path,prefix,m); + } else { + sprintf(outfname,"%s/%s_%06d.bin",path,output,m); + } outfile=fopen(outfname,"w"); // Loop over subints to dump @@ -225,17 +241,17 @@ int main(int argc,char *argv[]) if (informat=='i') { for (i=0;i0.95) printf("%5d %d/%d %.4f\n",orb.satno,nalt,nsel,frac); } rewind(fp); @@ -736,7 +736,7 @@ int main(int argc,char *argv[]) } else { rms=fit_curve(orb,ia); //printf("rms: %.3f kHz, cf: %.6f MHz, TCA: %s\n",rms,d.ffit/1000.0,nfd); - printf("%05d %.3f %.3f %s %04d %02d%010.6f\n",orb.satno,d.ffit/1000.0,rms,nfd,d.p[0].site_id,orb.ep_year-2000,orb.ep_day); + printf("%05d %.6f %.3f %s %04d %02d%010.6f\n",orb.satno,d.ffit/1000.0,rms,nfd,d.p[0].site_id,orb.ep_year-2000,orb.ep_day); redraw=1; } } diff --git a/rfplot.c b/rfplot.c index e370dca..6e2fdff 100644 --- a/rfplot.c +++ b/rfplot.c @@ -17,6 +17,10 @@ struct select { float x[NMAX],y[NMAX],w; float sigma; }; +struct data { + int n; + double *x, *y, *ym; +} gd; 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); @@ -30,6 +34,9 @@ void quadfit(float x[],float y[],int n,float a[]); struct trace locate_trace(struct spectrogram s,struct trace t,int site_id,float sigmamin,float wmin,float wmax,int graves); void filter(struct spectrogram s,int site_id,float sigma,int graves); void peakfind(struct spectrogram s,int site_id,int i0,int i1,int j0,int j1); +void versafit(int m,int n,double *a,double *da,double (*func)(double *),double dchisq,double tol,char *opt); +double chisq_gaussian(double a[]); +float fit_gaussian_point(struct spectrogram s,float x,float y,struct select sel,int site_id,int graves); int main(int argc,char *argv[]) { @@ -54,7 +61,7 @@ int main(int argc,char *argv[]) int ix=0,iy=0,isub=0; int i0,j0,i1,j1,jmax; float width=1500; - float x=0.0,y=0.0,x0=0.0,y0=0.0; + float x=0.0,y=0.0,x0=0.0,y0=0.0,yfit; char c; char path[128],xlabel[128],ylabel[64],filename[32],tlefile[128]; int sec,lsec,ssec; @@ -83,9 +90,15 @@ int main(int argc,char *argv[]) env=getenv("ST_TLEDIR"); sprintf(tlefile,"%s/bulk.tle",env); + // Set selection + isel=0; + sel.n=0; + sel.w=100.0; + sel.sigma=5.0; + // Read arguments if (argc>1) { - while ((arg=getopt(argc,argv,"p:f:w:s:l:b:z:hc:C:gm:o:"))!=-1) { + while ((arg=getopt(argc,argv,"p:f:w:s:l:b:z:hc:C:gm:o:S:W:"))!=-1) { switch (arg) { case 'p': @@ -111,6 +124,14 @@ int main(int argc,char *argv[]) case 'w': df0=(double) atof(optarg); break; + + case 'W': + sel.w=atof(optarg); + break; + + case 'S': + sel.sigma=atof(optarg); + break; case 'z': zmax=atof(optarg); @@ -180,12 +201,6 @@ int main(int argc,char *argv[]) // Set trace tf.n=0; - // Set selection - isel=0; - sel.n=0; - sel.w=100.0; - sel.sigma=5.0; - // Forever loop for (;;) { if (redraw==1) { @@ -315,7 +330,9 @@ int main(int argc,char *argv[]) printf("t Locate traces\n"); printf("s Select track points\n"); printf("u Undo track point selection\n"); + printf("S Set track selection width and sigma\n"); printf("f Fits selected track points (generates out.dat)\n"); + printf("F Fits selected wideband track points (generates out.dat)\n"); printf("g Filter points above a certain sigma limit (generates filter.dat)\n"); printf("G Find peaks in current window (generates peakfind.dat)\n"); printf("i Identify selected track points\n"); @@ -324,6 +341,7 @@ int main(int argc,char *argv[]) printf("b Increase dynamic range\n"); printf("Z Provide manual dynamic range limits\n"); printf("D/mid Mark point (appends to mark.dat)\n"); + printf("w Fit wideband point (appends to mark.dat)\n"); printf("c Center view\n"); printf("C Toggle through color maps\n"); printf("p/right Toggle overlays\n"); @@ -400,6 +418,13 @@ int main(int argc,char *argv[]) continue; } + // Fit + if (c=='F') { + tf=fit_gaussian_trace(s,sel,site_id,graves); + tf.site=site_id; + continue; + } + // Identify if (c=='i') { if (graves==0) @@ -530,6 +555,26 @@ int main(int argc,char *argv[]) fclose(file); } + // Fit single wideband point + if (c=='w') { + file=fopen("mark.dat","a"); + + // Fit point + yfit=fit_gaussian_point(s,x,y,sel,site_id,graves); + cpgpt1(x,yfit,17); + i=(int) floor(x); + f=s.freq-0.5*s.samp_rate+(double) yfit*s.samp_rate/(double) s.nchan; + if (s.mjd[i]>1.0) { + 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); + } + + // Mark if (c=='m') { i0=(int) floor(xmin); @@ -733,6 +778,9 @@ int main(int argc,char *argv[]) free(t[i].freq); free(t[i].za); } + free(gd.x); + free(gd.y); + free(gd.ym); return 0; } @@ -850,6 +898,8 @@ void usage(void) printf("-f Frequency to zoom into (Hz)\n"); printf("-w Bandwidth to zoom into (Hz)\n"); printf("-o Frequency offset to apply (Hz) [0]\n"); + printf("-W Track selection width (in pixels) [100]\n"); + printf("-S Track selection significance [5]\n"); printf("-C Site ID\n"); printf("-c TLE catalog\n"); printf("-g GRAVES data\n"); @@ -986,6 +1036,99 @@ struct trace fit_trace(struct spectrogram s,struct select sel,int site_id,int gr return t; } +// Fit gaussian trace +struct trace fit_gaussian_trace(struct spectrogram s,struct select sel,int site_id,int graves) +{ + int i,j,k,l,sn; + int i0,i1,j0,j1,jmax; + double f,chi2; + double a[4],da[4]; + float x,y,s1,s2,z,za,zs,zm,sigma,rms; + struct trace t; + FILE *file; + + // Set up trace + t.satno=99999; + t.n=(int) ceil(sel.x[sel.n-1]-sel.x[0]); + t.mjd=(double *) malloc(sizeof(double)*t.n); + t.freq=(double *) malloc(sizeof(double)*t.n); + t.za=(float *) malloc(sizeof(float)*t.n); + + // Set up data + gd.n=(int) 2 * sel.w; + gd.x=(double *) malloc(sizeof(double) * gd.n); + gd.y=(double *) malloc(sizeof(double) * gd.n); + gd.ym=(double *) malloc(sizeof(double) * gd.n); + + // Open file + file=fopen("out.dat","w"); + + // Loop over selected regions + for (k=0,l=0;k=s.nchan) + j1=s.nchan; + + // Loop over points + for (j=j0,l=0;jsel.sigma && s.mjd[i]>1.0) { + f=s.freq-0.5*s.samp_rate+(double) a[0]*s.samp_rate/(double) s.nchan; + 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) a[0], 17); + cpgmove((float) i,(float) a[0]-a[1]); + cpgdraw((float) i,(float) a[0]+a[1]); + t.mjd[l]=s.mjd[i]; + t.freq[l]=f; + t.za[l]=0.0; + l++; + } + } + } + t.n=l; + + // Close file + fclose(file); + + return t; +} + // Discrete convolution void convolve(float *y,int n,float *w,int m,float *z) { @@ -1286,3 +1429,71 @@ void peakfind(struct spectrogram s,int site_id,int i0,int i1,int j0,int j1) return; } + +// Gaussian point +double chisq_gaussian(double a[]) +{ + int i; + double arg,amp,chisq; + + // Compute chisq + for (i=0,chisq=0.0;i=s.nchan) + j1=s.nchan; + + // Loop over points + for (j=j0,l=0;j