diff --git a/makefile b/makefile
index 452df99..fed7c96 100644
--- a/makefile
+++ b/makefile
@@ -8,7 +8,7 @@ LFLAGS = -lcpgplot -lpgplot -lX11 -lpng -lm -lgsl -lgslcblas
 CC = gcc
 
 all:
-	make rfedit rfplot rffft rfpng
+	make rfedit rfplot rffft rfpng rffind
 
 rfpng: rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o
 	gfortran -o rfpng rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o $(LFLAGS)
diff --git a/rffind.c b/rffind.c
index 4f35883..40b424b 100644
--- a/rffind.c
+++ b/rffind.c
@@ -2,128 +2,101 @@
 #include <string.h>
 #include <stdlib.h>
 #include <math.h>
-#include <cpgplot.h>
 #include <getopt.h>
 #include "rftime.h"
 #include "rfio.h"
-#include <gsl/gsl_multifit.h>
 
 #define LIM 128
 #define NMAX 64
 
-
 void usage(void)
 {
 }
 
-float fit_polynomial(float *x,float *y,int n,int m,float *a)
+void filter(struct spectrogram s,int site_id,float sigma,char *filename)
 {
+  int i,j,k,l;
+  float s1,s2,avg,std,dz;
+  FILE *file;
+  double f;
+  int *mask;
 
-  int i,j;
-  double chisq;
-  gsl_matrix *X,*cov;
-  gsl_vector *yy,*w,*c;
-
-
-  X=gsl_matrix_alloc(n,m);
-  yy=gsl_vector_alloc(n);
-  w=gsl_vector_alloc(n);
-
-  c=gsl_vector_alloc(m);
-  cov=gsl_matrix_alloc(m,m);
-
-  // Fill matrices
-  for(i=0;i<n;i++) {
-    for (j=0;j<m;j++)
-      gsl_matrix_set(X,i,j,pow(x[i],j));
-    
-    gsl_vector_set(yy,i,y[i]);
-    gsl_vector_set(w,i,1.0);
-  }
-
-  // Do fit
-  gsl_multifit_linear_workspace *work=gsl_multifit_linear_alloc(n,m);
-  gsl_multifit_wlinear(X,w,yy,c,cov,&chisq,work);
-  gsl_multifit_linear_free(work);
-
-  // Save parameters
-  for (i=0;i<m;i++)
-    a[i]=gsl_vector_get(c,(i));
-
-  gsl_matrix_free(X);
-  gsl_vector_free(yy);
-  gsl_vector_free(w);
-  gsl_vector_free(c);
-  gsl_matrix_free(cov);
-
-  return chisq;
-}
+  mask=(int *) malloc(sizeof(int)*s.nchan);
 
-void filter(float *x,float *y,int n,int m,float sigma,int *mask)
-{
-  int i,j,k,nn;
-  float *xx,*yy,*a,chi2,ym;
-  float rms;
-
-  // Allocate
-  a=(float *) malloc(sizeof(float)*m);
-  xx=(float *) malloc(sizeof(float)*n);
-  yy=(float *) malloc(sizeof(float)*n);
-
-  // Set intial mask
-  for (i=0;i<n;i++) {
-    if (y[i]<2)
-      mask[i]=1;
-    else
-      mask[i]=0;
-  }
+  // Open file
+  file=fopen(filename,"a");
 
-  // Iterations
-  for (k=0;k<10;k++) {
-    // Apply mask
-    for (i=0,j=0;i<n;i+=100) {
-      if (mask[i]!=1)
-	continue;
-      xx[j]=x[i];
-      yy[j]=y[i];
-      j++;
+  // Loop over subints
+  for (i=0;i<s.nsub;i++) {
+    // Set mask
+    for (j=0;j<s.nchan;j++)
+      mask[j]=1;
+
+    // Iterate to remove outliers
+    for (k=0;k<10;k++) {
+
+      // Find average
+      for (j=0,s1=s2=0.0;j<s.nchan;j++) {
+	if (mask[j]==1) {
+	  s1+=s.z[i+s.nsub*j];
+	  s2+=1.0;
+	}
+      }
+      avg=s1/s2;
+      
+      // Find standard deviation
+      for (j=0,s1=s2=0.0;j<s.nchan;j++) {
+	if (mask[j]==1) {
+	  dz=s.z[i+s.nsub*j]-avg;
+	  s1+=dz*dz;
+	  s2+=1.0;
+	}
+      }
+      std=sqrt(s1/s2);
+
+      // Update mask
+      for (j=0,l=0;j<s.nchan;j++) {
+	if (fabs(s.z[i+s.nsub*j]-avg)>sigma*std) {
+	  mask[j]=0;
+	  l++;
+	}
+      }
     }
-    nn=j;
-    
-    // Fit polynomial
-    chi2=fit_polynomial(xx,yy,nn,m,a);
-    
-    // Scale
-    for (i=0,rms=0.0,nn=0;i<n;i++) {
-      for (j=0,ym=0.0;j<m;j++)
-	ym+=a[j]*pow(x[i],j);
-      yy[i]=y[i]/ym-1.0;
-      if (mask[i]==1) {
-	rms+=yy[i]*yy[i];
-	nn++;
+       // Reset mask
+    for (j=0;j<s.nchan;j++) {
+      if (s.z[i+s.nsub*j]-avg>sigma*std) 
+	mask[j]=1;
+      else
+	mask[j]=0;
+    }    
+
+    // Find maximum when points are adjacent
+    for (j=0;j<s.nchan-1;j++) {
+      if (mask[j]==1 && mask[j+1]==1) {
+	if (s.z[i+s.nsub*j]<s.z[i+s.nsub*(j+1)])
+	  mask[j]=0;
+      }
+    }
+    for (j=s.nchan-2;j>=0;j--) {
+      if (mask[j]==1 && mask[j-1]==1) {
+	if (s.z[i+s.nsub*j]<s.z[i+s.nsub*(j-1)])
+	  mask[j]=0;
       }
     }
-    rms=sqrt(rms/(float) (nn-1));
 
-    // Update mask
-    for (i=0;i<n;i++) {
-      if (fabs(yy[i])>sigma*rms)
-	mask[i]=0;
+    // Mark points
+    for (j=0;j<s.nchan;j++) {
+      if (mask[j]==1) {
+	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);
+      }
     }
-  }
+  } 
 
-  // Recompute final mask
-  for (i=0;i<n;i++) {
-    if (yy[i]>sigma*rms)
-      mask[i]=0;
-    else
-      mask[i]=1;
-  }
+  fclose(file);
 
-  // Free
-  free(a);
-  free(xx);
-  free(yy);
+  free(mask);
 
   return;
 }
@@ -139,8 +112,6 @@ int main(int argc,char *argv[])
   float avg,std;
   int arg=0;
   float sigma=5.0;
-  float *x,*y;
-  int *mask;
   FILE *file;
   double f,f0,df0;
   char filename[128]="find.dat";
@@ -155,7 +126,7 @@ int main(int argc,char *argv[])
 
   // Read arguments
   if (argc>1) {
-    while ((arg=getopt(argc,argv,"p:f:w:s:l:hc:o:"))!=-1) {
+    while ((arg=getopt(argc,argv,"p:f:w:s:l:hc:o:S:"))!=-1) {
       switch (arg) {
 	
       case 'p':
@@ -178,6 +149,10 @@ int main(int argc,char *argv[])
 	f0=(double) atof(optarg);
 	break;
 	
+      case 'S':
+	sigma=atof(optarg);
+	break;
+
       case 'w':
 	df0=(double) atof(optarg);
 	break;
@@ -199,42 +174,18 @@ int main(int argc,char *argv[])
   // Read data
   s=read_spectrogram(path,isub,nsub,f0,df0,1);
 
-  x=(float *) malloc(sizeof(float)*s.nchan);
-  y=(float *) malloc(sizeof(float)*s.nchan);
-  mask=(int *) malloc(sizeof(int)*s.nchan);
+  // Exit on emtpy file
+  if (s.nsub==0)
+    return 0;
 
   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);
 
-  file=fopen(filename,"w");
-
-  // Loop over subints
-  for (i=0;i<s.nsub;i++) {
-    // Fill array
-    for (j=0;j<s.nchan;j++) {
-      x[j]=-0.5*s.samp_rate*1e-6+s.samp_rate*1e-6*(float) j/(float) (s.nchan-1);
-      y[j]=s.z[i+s.nsub*j];
-    }
-
-    // Filter data
-    filter(x,y,s.nchan,m,sigma,mask);
-
-    // Output
-    for (j=0;j<s.nchan;j++) {
-      if (mask[j]==0) {
-	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);
-      }
-    }
-  }
-  fclose(file);
+  // Filter
+  filter(s,site_id,sigma,filename);
 
   // Free
   free(s.z);
   free(s.mjd);
-  free(x);
-  free(y);
-  free(mask);
 
   return 0;
 }
diff --git a/rfio.c b/rfio.c
index 12586f9..11f1eb8 100644
--- a/rfio.c
+++ b/rfio.c
@@ -6,7 +6,7 @@
 
 struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,double df0,int nbin)
 {
-  int i,j,k,l,flag=0,status,msub,ibin;
+  int i,j,k,l,flag=0,status,msub,ibin,nadd;
   char filename[128],header[256],nfd[32];
   FILE *file;
   struct spectrogram s;
@@ -23,6 +23,7 @@ struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,dou
   file=fopen(filename,"r");
   if (file==NULL) {
     printf("%s does not exist\n",filename);
+    s.nsub=0;
     return s;
   }
 
@@ -40,8 +41,12 @@ struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,dou
     j0=(int) ((f0-0.5*df0-s.freq+0.5*s.samp_rate)*(float) nch/s.samp_rate);
     j1=(int) ((f0+0.5*df0-s.freq+0.5*s.samp_rate)*(float) nch/s.samp_rate);
     
-    if (j0<0 || j1>nch) 
+    if (j0<0 || j1>nch) {
       fprintf(stderr,"Requested frequency range out of limits\n");
+      s.nsub=0;
+      s.nchan=0;
+      return s;
+    }
   } else {
     s.nchan=nch;
     j0=0;
@@ -64,7 +69,7 @@ struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,dou
     s.mjd[j]=0.0;
 
   // Loop over files
-  for (k=0,i=0,l=0,ibin=0;l<nsub;k++) {
+  for (k=0,i=0,l=0,ibin=0,nadd=0;l<nsub;k++) {
     // Generate filename
     sprintf(filename,"%s_%06d.bin",prefix,k+isub);
 
@@ -80,11 +85,13 @@ struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,dou
     for (;l<nsub;l++,ibin++) {
       // Read header
       status=fread(header,sizeof(char),256,file);
+
       if (status==0)
 	break;
       status=sscanf(header,"HEADER\nUTC_START    %s\nFREQ         %lf Hz\nBW           %lf Hz\nLENGTH       %f s\nNCHAN        %d\n",nfd,&freq,&samp_rate,&length,&nchan);
       s.mjd[i]+=nfd2mjd(nfd)+0.5*length/86400.0;
       s.length[i]+=length;
+      nadd++;
 
       // Read buffer
       status=fread(z,sizeof(float),nch,file);
@@ -98,12 +105,13 @@ struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,dou
       // Increment
       if (l%nbin==nbin-1) {
 	// Scale
-	s.mjd[i]/=(float) nbin;
-	
+	s.mjd[i]/=(float) nadd;
+
 	for (j=0;j<s.nchan;j++) 
-	  s.z[i+s.nsub*j]/=(float) nbin;
+	  s.z[i+s.nsub*j]/=(float) nadd;
 
 	ibin=0;
+	nadd=0;
 	i++;
       }
     }
@@ -111,6 +119,11 @@ struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,dou
     // Close file
     fclose(file);
   }
+  // Scale last subint
+  s.mjd[i]/=(float) nadd;
+
+  for (j=0;j<s.nchan;j++) 
+    s.z[i+s.nsub*j]/=(float) nadd;
 
   // Swap frequency range
   if (f0>0.0 && df0>0.0) {
diff --git a/rfplot.c b/rfplot.c
index 11c1306..35cb0b4 100644
--- a/rfplot.c
+++ b/rfplot.c
@@ -22,6 +22,95 @@ void usage(void);
 void plot_traces(struct trace *t,int nsat);
 struct trace locate_trace(struct spectrogram s,struct select sel,int site_id);
 
+void filter(struct spectrogram s,int site_id)
+{
+  int i,j,k,l,jmax,zmax;
+  float s1,s2,avg,std,dz;
+  FILE *file;
+  double f;
+  int *mask;
+  float sigma=5;
+
+  mask=(int *) malloc(sizeof(int)*s.nchan);
+
+  // Open file
+  file=fopen("filter.dat","w");
+
+  // Loop over subints
+  for (i=0;i<s.nsub;i++) {
+    // Set mask
+    for (j=0;j<s.nchan;j++)
+      mask[j]=1;
+
+    // Iterate to remove outliers
+    for (k=0;k<10;k++) {
+
+      // Find average
+      for (j=0,s1=s2=0.0;j<s.nchan;j++) {
+	if (mask[j]==1) {
+	  s1+=s.z[i+s.nsub*j];
+	  s2+=1.0;
+	}
+      }
+      avg=s1/s2;
+      
+      // Find standard deviation
+      for (j=0,s1=s2=0.0;j<s.nchan;j++) {
+	if (mask[j]==1) {
+	  dz=s.z[i+s.nsub*j]-avg;
+	  s1+=dz*dz;
+	  s2+=1.0;
+	}
+      }
+      std=sqrt(s1/s2);
+
+      // Update mask
+      for (j=0,l=0;j<s.nchan;j++) {
+	if (fabs(s.z[i+s.nsub*j]-avg)>sigma*std) {
+	  mask[j]=0;
+	  l++;
+	}
+      }
+    }
+    // Reset mask
+    for (j=0;j<s.nchan;j++) {
+      if (s.z[i+s.nsub*j]-avg>sigma*std) 
+	mask[j]=1;
+      else
+	mask[j]=0;
+    }    
+
+    // Find maximum when points are adjacent
+    for (j=0;j<s.nchan-1;j++) {
+      if (mask[j]==1 && mask[j+1]==1) {
+	if (s.z[i+s.nsub*j]<s.z[i+s.nsub*(j+1)])
+	  mask[j]=0;
+      }
+    }
+    for (j=s.nchan-2;j>=0;j--) {
+      if (mask[j]==1 && mask[j-1]==1) {
+	if (s.z[i+s.nsub*j]<s.z[i+s.nsub*(j-1)])
+	  mask[j]=0;
+      }
+    }
+
+    // Mark points
+    for (j=0;j<s.nchan;j++) {
+      if (mask[j]==1) {
+	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);
+	cpgpt1((float) i+0.5,(float) j+0.5,17);
+      }
+    }
+  }
+  fclose(file);
+
+  free(mask);
+
+  return;
+}
+
 int main(int argc,char *argv[])
 {
   struct spectrogram s;
@@ -48,7 +137,7 @@ int main(int argc,char *argv[])
   double f0=0.0,df0=0.0;
   int foverlay=1;
   struct trace *t,tf;
-  int nsat,satno;
+  int nsat,satno,status;
   struct select sel;
   char *env;
   int site_id=0;
@@ -120,9 +209,13 @@ int main(int argc,char *argv[])
 
   // Read data
   s=read_spectrogram(path,isub,nsub,f0,df0,nbin);
-
+ 
   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
+  if (s.nsub==0)
+    return 0;
+
   // Compute traces
   t=compute_trace(tlefile,s.mjd,s.nsub,site_id,s.freq*1e-6,s.samp_rate*1e-6,&nsat);
   printf("Traces for %d objects for location %d\n",nsat,site_id);
@@ -237,6 +330,9 @@ int main(int argc,char *argv[])
       continue;
     }
 
+    if (c=='g')
+      filter(s,site_id);
+
     // Fit
     if (c=='f') {
       tf=locate_trace(s,sel,site_id);
@@ -254,7 +350,7 @@ int main(int argc,char *argv[])
     // Identify
     if (c=='I') {
       printf("Provide satno: ");
-      scanf("%d",&satno);
+      status=scanf("%d",&satno);
       identify_trace(tlefile,tf,satno);
       redraw=1;
       continue;
@@ -399,7 +495,7 @@ int main(int argc,char *argv[])
 	j1=s.nchan-1;
 
       printf("Provide filename: ");
-      scanf("%s",filename);
+      status=scanf("%s",filename);
       
       file=fopen(filename,"a");
       // Loop over image
diff --git a/rfpng.c b/rfpng.c
new file mode 100644
index 0000000..04953c3
--- /dev/null
+++ b/rfpng.c
@@ -0,0 +1,422 @@
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include <cpgplot.h>
+#include <getopt.h>
+#include "rftime.h"
+#include "rfio.h"
+#include "rftrace.h"
+
+#define LIM 128
+#define NMAX 64
+
+struct select {
+  int flag,n;
+  float x[NMAX],y[NMAX],w;
+};
+
+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 locate_trace(struct spectrogram s,struct select sel,int site_id);
+
+int main(int argc,char *argv[])
+{
+  struct spectrogram s;
+  float tr[]={-0.5,1.0,0.0,-0.5,0.0,1.0};
+  float cool_l[]={-0.5,0.0,0.17,0.33,0.50,0.67,0.83,1.0,1.7};
+  float cool_r[]={0.0,0.0,0.0,0.0,0.6,1.0,1.0,1.0,1.0};
+  float cool_g[]={0.0,0.0,0.0,1.0,1.0,1.0,0.6,0.0,1.0};
+  float cool_b[]={0.0,0.3,0.8,1.0,0.3,0.0,0.0,0.0,1.0};
+  float xmin,xmax,ymin,ymax,zmin,zmax=8.0;
+  int i,j,k;
+  float dt,zzmax,s1,s2;
+  int ix=0,iy=0,isub=0;
+  int i0,j0,i1,j1,jmax;
+  float width=1500;
+  float x,y,x0,y0;
+  char c;
+  char path[128],xlabel[64],ylabel[64],filename[32],tlefile[128],pngfile[128];
+  int sec,lsec,ssec;
+  char stime[16];
+  double fmin,fmax,fcen,f;
+  FILE *file;
+  int arg=0,nsub=900,nbin=1;
+  double f0=0.0,df0=0.0,dy=2500;
+  int foverlay=1;
+  struct trace *t,tf;
+  int nsat,satno;
+  struct select sel;
+  char *env;
+  int site_id=0;
+
+  // Get site
+  env=getenv("ST_COSPAR");
+  if (env!=NULL) {
+    site_id=atoi(env);
+  } else {
+    printf("ST_COSPAR environment variable not found.\n");
+  }
+  env=getenv("ST_TLEDIR");
+  sprintf(tlefile,"%s/bulk.tle",env); 
+
+  // Read arguments
+  if (argc>1) {
+    while ((arg=getopt(argc,argv,"p:f:w:s:l:b:z:hc:C:"))!=-1) {
+      switch (arg) {
+	
+      case 'p':
+	strcpy(path,optarg);
+	break;
+	
+      case 'f':
+	f0=(double) atof(optarg);
+	break;
+	
+      case 'w':
+	df0=(double) atof(optarg);
+	break;
+	
+      case 'z':
+	zmax=atof(optarg);
+	break;
+
+      case 'c':
+	strcpy(tlefile,optarg);
+	break;
+
+      case 'h':
+	usage();
+	return 0;
+
+      default:
+	usage();
+	return 0;
+      }
+    }
+  } else {
+    usage();
+    return 0;
+  }
+
+  for (isub=0;;isub+=15) {
+    // Read data
+    s=read_spectrogram(path,isub,nsub,f0,df0,nbin);
+    if (s.mjd[0]<54000)
+      break;
+
+    // Output filename
+    sprintf(pngfile,"%.19s_%8.3f.png/png",s.nfd0,s.freq*1e-6);
+
+    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);
+    printf("Traces for %d objects for location %d\n",nsat,site_id);
+
+    printf("%s\n",pngfile);
+    cpgopen(pngfile);
+    cpgctab(cool_l,cool_r,cool_g,cool_b,9,1.0,0.5);
+    cpgsch(0.8);
+    cpgask(1);
+    
+    // Default limits
+    xmin=0.0;
+    xmax=(float) s.nsub;
+    ymin=0;
+    ymax=(float) s.nchan;
+    zmin=0.0;
+    
+    cpgsci(1);
+    cpgpage();
+    cpgsvp(0.1,0.95,0.1,0.95);
+    cpgswin(xmin,xmax,ymin,ymax);
+  
+    cpgimag(s.z,s.nsub,s.nchan,1,s.nsub,1,s.nchan,zmin,zmax,tr);
+    
+    // Pixel axis
+    cpgbox("CTSM1",0.,0,"CTSM1",0.,0);
+    
+    // Time axis
+    cpgbox("B",0.,0,"",0.,0);
+    time_axis(s.mjd,s.nsub,xmin,xmax,ymin,ymax);
+    
+    // Freq axis
+    fmin=s.freq-0.5*s.samp_rate+ymin*s.samp_rate/(float) s.nchan;
+    fmax=s.freq-0.5*s.samp_rate+ymax*s.samp_rate/(float) s.nchan;
+    fmin*=1e-6;
+    fmax*=1e-6;
+    
+    // Plot traces
+    cpgswin(xmin,xmax,fmin,fmax);
+    cpgsci(3);
+    plot_traces(t,nsat);
+    cpgsci(1);
+    
+    // Human readable frequency axis
+    fcen=0.5*(fmax+fmin);
+    fcen=floor(1000*fcen)/1000.0;
+    sprintf(ylabel,"Frequency - %.3f MHz",fcen);
+    fmin-=fcen;
+    fmax-=fcen;
+    cpgswin(xmin,xmax,fmin,fmax);
+    cpgbox("",0.,0,"BTSN",0.,0);
+    
+    sprintf(xlabel,"UT Date: %.10s",s.nfd0);
+    cpglab(xlabel,ylabel," ");
+    
+    cpgswin(xmin,xmax,ymin,ymax);
+   
+    cpgend();
+
+    // Free
+    free(s.z);
+    free(s.mjd);
+
+    for (i=0;i<nsat;i++) {
+      free(t[i].mjd);
+      free(t[i].freq);
+      free(t[i].za);
+    }
+    //    free(tf.mjd);
+    //    free(tf.freq);
+    //    free(tf.za);
+
+  }
+
+  return 0;
+}
+
+// Convert Decimal into Sexagesimal
+void dec2sex(double x,char *s,int f,int len)
+{
+  int i;
+  double sec,deg,min;
+  char sign;
+  char *form[]={"::",",,","hms","  "};
+
+  sign=(x<0 ? '-' : ' ');
+  x=3600.*fabs(x);
+
+  sec=fmod(x,60.);
+  x=(x-sec)/60.;
+  min=fmod(x,60.);
+  x=(x-min)/60.;
+  //  deg=fmod(x,60.);
+  deg=x;
+
+  if (len==7) sprintf(s,"%c%02i%c%02i%c%07.4f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
+  if (len==6) sprintf(s,"%c%02i%c%02i%c%06.3f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
+  if (len==5) sprintf(s,"%c%02i%c%02i%c%05.2f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
+  if (len==4) sprintf(s,"%c%02i%c%02i%c%04.1f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
+  if (len==2) sprintf(s,"%c%02i%c%02i%c%02i%c",sign,(int) deg,form[f][0],(int) min,form[f][1],(int) floor(sec),form[f][2]);
+
+  return;
+}
+
+void time_axis(double *mjd,int n,float xmin,float xmax,float ymin,float ymax)
+{
+  int i,imin,imax;
+  double mjdt,mjdmin,mjdmax;
+  float dt,t,tmin,tmax;
+  int lsec,ssec,sec;
+  char stime[16];
+
+  // Find extrema
+  for (i=0;i<n;i++) {
+    if (i==0) {
+      mjdmin=mjd[i];
+      mjdmax=mjd[i];
+    } else {
+      if (mjd[i]>mjdmax) mjdmax=mjd[i];
+    }
+  }
+  dt=(float) 86400*(mjdmax-mjdmin);
+
+  // Choose tickmarks
+  if (dt>43000) {
+    lsec=10800;
+    ssec=3600;
+  } else if (dt>21600) {
+    lsec=10800;
+    ssec=3600;
+  } else if (dt>7200) {
+    lsec=1800;
+    ssec=300;
+  } else if (dt>3600) {
+    lsec=600;
+    ssec=120;
+  } else if (dt>900) {
+    lsec=300;
+    ssec=60;
+  } else {
+    lsec=60;
+    ssec=10;
+  }
+  // Override
+  lsec=60;
+  ssec=10;
+
+  // Extrema
+  tmin=86400.0*(mjdmin-floor(mjdmin));
+  tmax=tmin+dt;
+  tmin=lsec*floor(tmin/lsec);
+  tmax=lsec*ceil(tmax/lsec);
+
+  // Large tickmarks
+  for (t=tmin;t<=tmax;t+=lsec) {
+    mjdt=floor(mjdmin)+t/86400.0;
+    if (mjdt>=mjdmin && mjdt<mjdmax) {
+      for (i=0;i<n-1;i++)
+	if (mjdt>=mjd[i] && mjdt<mjd[i+1])
+	  break;
+      sec=(int) floor(fmod(t,86400.0));
+      dec2sex(((float) sec+0.1)/3600.0,stime,0,2);
+      stime[6]='\0';
+      cpgtick(xmin,ymin,xmax,ymin,((float) i-xmin)/(xmax-xmin),0.5,0.5,0.3,0.0,stime);
+    }
+  }
+
+  // Small tickmarks
+  for (t=tmin;t<=tmax;t+=ssec) {
+    mjdt=floor(mjdmin)+t/86400.0;
+    if (mjdt>=mjdmin && mjdt<mjdmax) {
+      for (i=0;i<n-1;i++)
+	if (mjdt>=mjd[i] && mjdt<mjd[i+1])
+	  break;
+      sec=(int) floor(t);
+      cpgtick(xmin,ymin,xmax,ymin,((float) i-xmin)/(xmax-xmin),0.25,0.25,1.0,1.0,"");
+    }
+  }
+
+  return;
+}
+
+void usage(void)
+{
+  printf("rfplot: plot RF observations\n\n");
+  printf("-p <path>    Input path to file /a/b/c_??????.bin\n");
+  printf("-s <start>   Number of starting subintegration [0]\n");
+  printf("-l <length>  Number of subintegrations to plot [3600]\n");
+  printf("-b <nbin>    Number of subintegrations to bin [1]\n");
+  printf("-z <zmax>    Image scaling upper limit [8.0]\n");
+  printf("-f <freq>    Frequency to zoom into (Hz)\n");
+  printf("-w <bw>      Bandwidth to zoom into (Hz)\n");
+  printf("-h           This help\n");
+
+  return;
+}
+
+void plot_traces(struct trace *t,int nsat)
+{
+  int i,j,flag,textflag;
+  char text[8];
+
+  // Loop over objects
+  for (i=0;i<nsat;i++) {
+    sprintf(text," %d",t[i].satno);
+
+    // Plot label at start of trace
+    if (t[i].za[0]<=90.0)
+	cpgtext(0.0,(float) t[i].freq[0],text);
+
+    // Loop over trace
+    for (j=0,flag=0,textflag=0;j<t[i].n;j++) {
+      // Plot label for rising sources
+      if (j>0 && t[i].za[j-1]>90.0 && t[i].za[j]<=90.0)
+	cpgtext((float) j,(float) t[i].freq[j],text);
+
+      // Plot line
+      if (flag==0) {
+	cpgmove((float) j,(float) t[i].freq[j]);
+	flag=1;
+      } else {
+	cpgdraw((float) j,(float) t[i].freq[j]);
+      }
+
+      // Below horizon
+      if (t[i].za[j]>90.0)
+	flag=0;
+      else
+	flag=1;
+    }	  
+  }
+
+  return;
+}
+
+// Locate trace
+struct trace locate_trace(struct spectrogram s,struct select sel,int site_id)
+{
+  int i,j,k,l,sn;
+  int i0,i1,j0,j1,jmax;
+  double f;
+  float x,y,s1,s2,z,za,zs,zm,sigma;
+  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);
+
+  // Open file
+  file=fopen("out.dat","w");
+
+  // Loop over selected regions
+  for (k=0,l=0;k<sel.n-1;k++) {
+    for (x=sel.x[k];x<=sel.x[k+1];x+=1.0) {
+      y=(x-sel.x[k])/(sel.x[k+1]-sel.x[k])*(sel.y[k+1]-sel.y[k])+sel.y[k];
+      i=(int) floor(x);
+      j0=(int) floor(y-sel.w);
+      j1=(int) floor(y+sel.w);
+
+      // Keep in range
+      if (j0<0)
+	j0=0;
+      if (j1>=s.nchan)
+	j1=s.nchan;
+
+      // Find maximum and significance
+      zm=0.0;
+      jmax=0;
+      s1=0.0;
+      s2=0.0;
+      sn=0;
+      for (j=j0;j<j1;j++) {
+	z=s.z[i+s.nsub*j];
+	s1+=z;
+	s2+=z*z;
+	sn++;
+	if (z>zm) {
+	  zm=z;
+	  jmax=j;
+	}
+      }
+      za=s1/(float) sn;
+      zs=sqrt(s2/(float) sn-za*za);
+      sigma=(zm-za)/zs;
+
+      // 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);
+	cpgpt1((float) i,(float) jmax,17);
+	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;
+}
diff --git a/rftrace.c b/rftrace.c
index cf74975..66daf46 100644
--- a/rftrace.c
+++ b/rftrace.c
@@ -96,7 +96,7 @@ void obspos_xyz(double mjd,double lng,double lat,float alt,xyz_t *pos,xyz_t *vel
 // Get observing site
 struct site get_site(int site_id)
 {
-  int i=0;
+  int i=0,status;
   char line[LIM];
   FILE *file;
   int id;
@@ -123,7 +123,7 @@ struct site get_site(int site_id)
     line[strlen(line)-1]='\0';
 
     // Read data
-    sscanf(line,"%4d %2s %lf %lf %f",
+    status=sscanf(line,"%4d %2s %lf %lf %f",
 	   &id,abbrev,&lat,&lng,&alt);
     strcpy(observer,line+38);
 
@@ -148,7 +148,7 @@ struct site get_site(int site_id)
 // Identify trace
 void identify_trace(char *tlefile,struct trace t,int satno)
 {
-  int i,imode,flag=0;
+  int i,imode,flag=0,status;
   struct point *p;
   struct site s;
   double *v;
@@ -246,7 +246,7 @@ void identify_trace(char *tlefile,struct trace t,int satno)
     printf("\nBest fitting object:\n");
     printf("%05d: %s  %8.3f MHz %8.3f kHz\n",satnomin,nfdmin,1e-6*freqmin,1e-3*rmsmin);
     printf("Store frequency? [y/n]\n");
-    scanf("%s",text);
+    status=scanf("%s",text);
     if (text[0]=='y') {
       file=fopen(freqlist,"a");
       fprintf(file,"%05d %8.3f\n",satnomin,1e-6*freqmin);
@@ -270,7 +270,7 @@ void identify_trace(char *tlefile,struct trace t,int satno)
 // Compute trace
 struct trace *compute_trace(char *tlefile,double *mjd,int n,int site_id,float freq,float bw,int *nsat)
 {
-  int i,j,imode,flag,satno,tflag,m;
+  int i,j,imode,flag,satno,tflag,m,status;
   struct point *p;
   struct site s;
   FILE *file,*infile;
@@ -299,7 +299,7 @@ struct trace *compute_trace(char *tlefile,double *mjd,int n,int site_id,float fr
   for (i=0;;) {
     if (fgetline(infile,line,LIM)<=0)
       break;
-    sscanf(line,"%d %lf",&satno,&freq0);
+    status=sscanf(line,"%d %lf",&satno,&freq0);
 
     if (freq0>=fmin && freq0<=fmax)
       i++;
@@ -334,7 +334,7 @@ struct trace *compute_trace(char *tlefile,double *mjd,int n,int site_id,float fr
   for (j=0;;) {
     if (fgetline(infile,line,LIM)<=0)
       break;
-    sscanf(line,"%d %lf",&satno,&freq0);
+    status=sscanf(line,"%d %lf",&satno,&freq0);
     if (freq0<fmin || freq0>fmax)
       continue;