From 74392761006218438429acc82716a376a62c4558 Mon Sep 17 00:00:00 2001
From: Cees Bassa <cbassa@gmail.com>
Date: Mon, 24 Mar 2014 16:41:38 +0100
Subject: [PATCH] Identification added

---
 rfplot.c  | 268 +++++++++++++++++++++++++++++++++++++++++++-----------
 rftrace.c | 160 ++++++++++++++++++++++++++++----
 rftrace.h |   3 +-
 3 files changed, 361 insertions(+), 70 deletions(-)

diff --git a/rfplot.c b/rfplot.c
index 3d9d024..be12735 100644
--- a/rfplot.c
+++ b/rfplot.c
@@ -9,48 +9,18 @@
 #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)
-{
-  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;
-}
+void plot_traces(struct trace *t,int nsat);
+struct trace locate_trace(struct spectrogram s,struct select sel);
 
 int main(int argc,char *argv[])
 {
@@ -61,12 +31,12 @@ int main(int argc,char *argv[])
   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,flag=0;
+  int i,j,k,flag=0,isel=0;
   int redraw=1,mode=0,posn=0,click=0;
   float dt,zzmax,s1,s2;
   int ix=0,iy=0,isub=0;
   int i0,j0,i1,j1,jmax;
-  float width=500;
+  float width=1500;
   float x,y,x0,y0;
   char c;
   char path[128],xlabel[64],ylabel[64];
@@ -77,8 +47,9 @@ int main(int argc,char *argv[])
   int arg=0,nsub=3600,nbin=1;
   double f0=0.0,df0=0.0;
   int foverlay=1;
-  struct trace *t;
-  int nsat;
+  struct trace *t,tf;
+  int nsat,satno;
+  struct select sel;
 
   // Read arguments
   if (argc>1) {
@@ -130,11 +101,9 @@ int main(int argc,char *argv[])
   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);
-  fmin=(s.freq-0.5*s.samp_rate)*1e-6;
-  fmax=(s.freq+0.5*s.samp_rate)*1e-6;
 
   // Compute traces
-  t=compute_trace(s.mjd,s.nsub,4171,fmin,fmax,&nsat);
+  t=compute_trace(s.mjd,s.nsub,4171,s.freq*1e-6,s.samp_rate*1e-6,&nsat);
   printf("Traces for %d objects\n",nsat);
 
   cpgopen("/xs");
@@ -147,7 +116,15 @@ int main(int argc,char *argv[])
   ymin=0.0;
   ymax=(float) s.nchan;
   zmin=0.0;
+
+  // Set trace
+  tf.n=0;
   
+  // Set selection
+  isel=0;
+  sel.n=0;
+  sel.w=50.0;
+
   // Forever loop
   for (;;) {
     if (redraw==1) {
@@ -193,6 +170,30 @@ int main(int argc,char *argv[])
 
       cpgswin(xmin,xmax,ymin,ymax);
 
+      // Plot selection
+      if (sel.n>0) {
+	cpgsci(7);
+	// Plot points
+	for (i=0;i<sel.n;i++) 
+	  cpgpt1(sel.x[i],sel.y[i],17);
+	// Plot upper bound
+	for (i=0;i<sel.n;i++) {
+	  if (i==0)
+	    cpgmove(sel.x[i],sel.y[i]+sel.w);
+	  else
+	    cpgdraw(sel.x[i],sel.y[i]+sel.w);
+	}
+	// Plot lower bound
+	for (i=0;i<sel.n;i++) {
+	  if (i==0)
+	    cpgmove(sel.x[i],sel.y[i]-sel.w);
+	  else
+	    cpgdraw(sel.x[i],sel.y[i]-sel.w);
+	}
+
+	cpgsci(1);
+      }
+
       redraw=0;
     }
 
@@ -203,6 +204,46 @@ int main(int argc,char *argv[])
     if (c=='q')
       break;
 
+    // Select start
+    if (c=='s') {
+      sel.x[isel]=x;
+      sel.y[isel]=y;
+      isel++;
+      sel.n=isel;
+      redraw=1;
+      continue;
+    }
+
+    // Fit
+    if (c=='f') {
+      tf=locate_trace(s,sel);
+      tf.site=4171;
+      continue;
+    }
+
+    // Identify
+    if (c=='i') {
+      identify_trace(tf,0);
+      redraw=1;
+      continue;
+    }
+
+    // Identify
+    if (c=='I') {
+      printf("Provide satno: ");
+      scanf("%d",&satno);
+      identify_trace(tf,satno);
+      redraw=1;
+      continue;
+    }
+
+    // Undo
+    if (c=='u') {
+      isel--;
+      sel.n=isel;
+      redraw=1;
+    }
+
     // Increase
     if (c=='v') {
       zmax*=1.01;
@@ -331,7 +372,7 @@ int main(int argc,char *argv[])
     }
 
     // Toggle overlay
-    if (c=='p') {
+    if (c=='p' || c=='X') {
       if (foverlay==0)
 	foverlay=1;
       else if (foverlay==1)
@@ -372,12 +413,21 @@ int main(int argc,char *argv[])
       continue;
     }
 
+    // Recompute traces
+    if (c=='R') {
+      t=compute_trace(s.mjd,s.nsub,4171,s.freq*1e-6,s.samp_rate*1e-6,&nsat);
+      redraw=1;
+      continue;
+    }
+
     // Reset
     if (c=='r') {
       xmin=0.0;
       xmax=(float) s.nsub;
       ymin=0.0;
       ymax=(float) s.nchan;
+      isel=0;
+      sel.n=0;
       redraw=1;
       continue;
     }
@@ -394,10 +444,10 @@ int main(int argc,char *argv[])
       // Set area
       x=width*(ix+0.5);
       y=width*(iy+0.5);
-      xmin=x-width;
-      xmax=x+width;
-      ymin=y-width;
-      ymax=y+width;
+      xmin=x-0.75*width;
+      xmax=x+0.75*width;
+      ymin=y-0.75*width;
+      ymax=y+0.75*width;
 
       // Increment
       iy++;
@@ -443,8 +493,16 @@ int main(int argc,char *argv[])
   // Free
   free(s.z);
   free(s.mjd);
-
-
+  if (tf.n>0) {
+    free(tf.mjd);
+    free(tf.freq);
+    free(tf.za);
+  }
+  for (i=0;i<nsat;i++) {
+    free(t[i].mjd);
+    free(t[i].freq);
+    free(t[i].za);
+  }
 
   return 0;
 }
@@ -565,3 +623,109 @@ void usage(void)
 
   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 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);
+
+      // 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 4171\n",s.mjd[i],f,sigma);
+	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 d9aaf59..4c8f849 100644
--- a/rftrace.c
+++ b/rftrace.c
@@ -1,8 +1,10 @@
 #include <stdio.h>
 #include <string.h>
 #include <math.h>
+#include <stdlib.h>
 #include "sgdp4h.h"
 #include "satutl.h"
+#include "rftrace.h"
 
 #define LIM 80
 #define D2R M_PI/180.0
@@ -21,12 +23,6 @@ struct site {
   float alt;
   char observer[64];
 };
-struct trace {
-  int satno,n;
-  double *mjd;
-  double *freq;
-  float *za;
-};
 
 // Return x modulo y [0,y)
 double modulo(double x,double y)
@@ -149,19 +145,141 @@ struct site get_site(int site_id)
   return s;
 }
 
-// Compute trace
-struct trace *compute_trace(double *mjd,int n,int site_id,float fmin,float fmax,int *nsat)
+// Identify trace
+void identify_trace(struct trace t,int satno)
 {
-  int i,j,imode,flag,satno,tflag;
+  int i,imode,flag=0;
+  struct point *p;
+  struct site s;
+  double *v;
+  orbit_t orb;
+  xyz_t satpos,satvel;
+  FILE *file;
+  double dx,dy,dz,dvx,dvy,dvz,r,za;
+  double sum1,sum2,beta,freq0,rms,mjd0;
+  char nfd[32],nfdmin[32],text[16];
+  int satnomin;
+  double rmsmin,freqmin;
+  
+  // Reloop stderr
+  freopen("/tmp/stderr.txt","w",stderr);
+
+  // Get site
+  s=get_site(t.site);
+
+  // Allocate
+  p=(struct point *) malloc(sizeof(struct point)*t.n);
+  v=(double *) malloc(sizeof(double)*t.n);
+
+  // Get observer position
+  for (i=0;i<t.n;i++) 
+    obspos_xyz(t.mjd[i],s.lng,s.lat,s.alt,&p[i].obspos,&p[i].obsvel);
+
+  printf("Fitting trace:\n");
+
+  // 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,sum1=0.0,sum2=0.0;i<t.n;i++) {
+      // Get satellite position
+      satpos_xyz(t.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[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;
+    }
+    freq0=sum1/sum2;
+
+    // Compute residuals
+    for (i=0,rms=0.0;i<t.n;i++) 
+      rms+=pow(t.freq[i]-(1.0-v[i]/C)*freq0,2);
+    rms=sqrt(rms/(double) t.n);
+
+    // Find TCA
+    for (i=0,mjd0=0.0;i<t.n-1;i++) 
+      if (v[i]*v[i-1]<0.0)
+	mjd0=t.mjd[i];
+    
+    if (mjd0>0.0)
+      mjd2nfd(mjd0,nfd);
+    else
+      strcpy(nfd,"0000-00-00T00:00:00");
+    
+    if (rms<1000) {
+      printf("%05d: %s  %8.3f MHz %8.3f kHz\n",orb.satno,nfd,1e-6*freq0,1e-3*rms);
+      if (flag==0 || rms<rmsmin) {
+	satnomin=orb.satno;
+	strcpy(nfdmin,nfd);
+	freqmin=freq0;
+	rmsmin=rms;
+	flag=1;
+      }
+    }
+  }
+  fclose(file);
+  fclose(stderr);
+
+  if (flag==1) {
+    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);
+    if (text[0]=='y') {
+      file=fopen("frequencies.txt","a");
+      fprintf(file,"%05d %8.3f\n",satnomin,1e-6*freqmin);
+      fclose(file);
+      file=fopen("log.txt","a");
+      fprintf(file,"%05d %8.3f %.3f %.19s\n",satnomin,1e-6*freqmin,1e-3*rmsmin,nfdmin);
+      fclose(file);
+      printf("Frequency stored\n\n");
+    }
+  } else {
+    printf("\nTrace not identified..\n");
+  }
+
+  // Free
+  free(p);
+  free(v);
+
+  return;
+}
+
+// Compute trace
+struct trace *compute_trace(double *mjd,int n,int site_id,float freq,float bw,int *nsat)
+{
+  int i,j,imode,flag,satno,tflag,m;
   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;
+  double freq0;
   char line[LIM],text[8];
   struct trace *t;
+  float fmin,fmax;
+
+  // Frequency limits
+  fmin=freq-0.5*bw;
+  fmax=freq+0.5*bw;
+
+  // Reloop stderr
+  freopen("/tmp/stderr.txt","w",stderr);
 
   // Find number of satellites in frequency range
   infile=fopen("frequencies.txt","r");
@@ -176,6 +294,12 @@ struct trace *compute_trace(double *mjd,int n,int site_id,float fmin,float fmax,
   fclose(infile);
   *nsat=i;
 
+  // Valid MJDs
+  for (i=0;i<n;i++)
+    if (mjd[i]==0.0)
+      break;
+  m=i;
+
   // Allocate traces
   t=(struct trace *) malloc(sizeof(struct trace)* *nsat);
 
@@ -183,10 +307,10 @@ struct trace *compute_trace(double *mjd,int n,int site_id,float fmin,float fmax,
   s=get_site(site_id);
 
   // Allocate
-  p=(struct point *) malloc(sizeof(struct point)*n);
+  p=(struct point *) malloc(sizeof(struct point)*m);
 
   // Get observer position
-  for (i=0;i<n;i++) 
+  for (i=0;i<m;i++) 
     obspos_xyz(mjd[i],s.lng,s.lat,s.alt,&p[i].obspos,&p[i].obsvel);
 
   infile=fopen("frequencies.txt","r");
@@ -199,10 +323,11 @@ struct trace *compute_trace(double *mjd,int n,int site_id,float fmin,float fmax,
 
     // Allocate
     t[j].satno=satno;
-    t[j].n=n;
-    t[j].mjd=(double *) malloc(sizeof(double)*n);
-    t[j].freq=(double *) malloc(sizeof(double)*n);
-    t[j].za=(float *) malloc(sizeof(float)*n);
+    t[j].site=site_id;
+    t[j].n=m;
+    t[j].mjd=(double *) malloc(sizeof(double)*m);
+    t[j].freq=(double *) malloc(sizeof(double)*m);
+    t[j].za=(float *) malloc(sizeof(float)*m);
 
     sprintf(text," %d",satno);
     // Loop over TLEs
@@ -214,7 +339,7 @@ struct trace *compute_trace(double *mjd,int n,int site_id,float fmin,float fmax,
 	printf("Error\n");
       
       // Loop over points
-      for (i=0,flag=0,tflag=0;i<n;i++) {
+      for (i=0,flag=0,tflag=0;i<m;i++) {
 	// Get satellite position
 	satpos_xyz(mjd[i]+2400000.5,&satpos,&satvel);
 	
@@ -241,6 +366,7 @@ struct trace *compute_trace(double *mjd,int n,int site_id,float fmin,float fmax,
     j++;
   }
   fclose(infile);
+  fclose(stderr);
 
   // Free
   free(p);
diff --git a/rftrace.h b/rftrace.h
index 6eebb5e..7d399b1 100644
--- a/rftrace.h
+++ b/rftrace.h
@@ -1,7 +1,8 @@
 struct trace {
-  int satno,n;
+  int satno,n,site;
   double *mjd;
   double *freq;
   float *za;
 };
 struct trace *compute_trace(double *mjd,int n,int site_id,float fmin,float fmax,int *nsat);
+void identify_trace(struct trace t,int satno);