Misc updates

makefile

@@ -14,14 +14,17 @@ exec_prefix = $(prefix)
 bindir = $(exec_prefix)/bin
 
 all:
-	make rfedit rfplot rffft rfpng rffit rffind
+	make rfedit rfplot rffft rfpng rffit rffind rfdop
 
 rffit: rffit.o sgdp4.o satutl.o deep.o ferror.o dsmin.o simplex.o versafit.o
 	gfortran -o rffit rffit.o sgdp4.o satutl.o deep.o ferror.o dsmin.o simplex.o versafit.o $(LFLAGS)
 
 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)
 
+rfdop: rfdop.o rftrace.o rfio.o rftime.o sgdp4.o satutl.o deep.o ferror.o
+	$(CC) -o rfdop rfdop.o rftrace.o rfio.o rftime.o sgdp4.o satutl.o deep.o ferror.o -lm
+
 rfedit: rfedit.o rfio.o rftime.o
 	$(CC) -o rfedit rfedit.o rfio.o rftime.o -lm
 

rfdop.c

@@ -0,0 +1,115 @@
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <stdlib.h>
+#include <getopt.h>
+#include "rfio.h"
+#include "rftime.h"
+#include "rftrace.h"
+
+void usage(void)
+{
+  printf("rfdop m:t:c:l:d:g:s:o:hi:\n\n");
+  printf("m    date/time (MJD)\n");
+  printf("t    date/time (yyyy-mm-ddThh:mm:ss.sss)\n");
+  printf("i    NORAD number\n");
+  printf("c    TLE catalog file\n");
+  printf("s    site (COSPAR)\n");
+  printf("d    time step [default: 10s]\n");
+  printf("l    trail length [default: 900s]\n");
+  printf("g    enable GRAVES\n");
+  printf("H    skip high orbits (<10 revs/day)\n");
+  printf("o    output file name\n");
+  printf("h    this help\n");
+
+  return;
+}
+
+int main(int argc, char *argv[])
+{
+  int i, n, satno=0;
+  double mjd0, *mjd, t;
+  float length=900.0, dt=10.0;
+  int site_id;
+  int graves=0,skiphigh=0;
+  char *env,tlefile[256],nfd[64],outfname[256]="out.dat";
+  int arg=0;
+
+  // Get defaults
+  env=getenv("ST_COSPAR");
+  if (env!=NULL) {
+    site_id=atoi(env);
+  }
+
+  if (argc>1) {
+    while ((arg=getopt(argc,argv,"m:t:c:i:l:d:gs:o:hH"))!=-1) {
+      switch (arg) {
+      case 't':
+	strcpy(nfd,optarg);
+	mjd0=nfd2mjd(nfd);
+	break;
+
+      case 'm':
+	mjd0=(double) atof(optarg);
+	break;
+
+      case 'c':
+	strcpy(tlefile,optarg);
+	break;
+
+      case 'i':
+	satno=atoi(optarg);
+	break;
+
+      case 'l':
+	length=atof(optarg);
+	break;
+
+      case 'd':
+	dt=atof(optarg);
+	break;
+
+      case 's':
+	site_id=atoi(optarg);
+	break;
+
+      case 'g':
+	graves=1;
+	break;
+
+      case 'H':
+	skiphigh=1;
+	break;
+
+      case 'o':
+	strcpy(outfname,optarg);
+	break;
+
+      case 'h':
+	usage();
+	return 0;
+
+      default:
+	usage();
+	return 0;
+      }
+    }
+  } else {
+    usage();
+    return 0;
+  }
+
+  // Generate MJDs
+  n = (int) floor(length/dt);
+  mjd = (double *) malloc(sizeof(double)*n);
+  for (i=0;i<n;i++) 
+    mjd[i] = mjd0+(double) i*dt/86400.0;
+
+  // Compute Doppler curves
+  compute_doppler(tlefile,mjd,n,site_id,satno,graves,skiphigh,outfname);
+
+  // Free
+  free(mjd);
+
+  return 0;
+}

rfplot.c

@@ -56,7 +56,7 @@ int main(int argc,char *argv[])
   float width=1500;
   float x,y,x0,y0;
   char c;
-  char path[128],xlabel[64],ylabel[64],filename[32],tlefile[128];
+  char path[128],xlabel[128],ylabel[64],filename[32],tlefile[128];
   int sec,lsec,ssec;
   char stime[16];
   double fmin,fmax,fcen,f;

rftrace.c

@@ -635,3 +635,110 @@ struct trace *compute_trace(char *tlefile,double *mjd,int n,int site_id,float fr
 
   return t;
 }
+
+// Compute trace
+void compute_doppler(char *tlefile,double *mjd,int n,int site_id,int satno,int graves, int skiphigh, char *outfname)
+{
+  int i,j,imode,flag,tflag,m,status;
+  struct point *p;
+  struct site s,sg;
+  FILE *file,*outfile;
+  orbit_t orb;
+  xyz_t satpos,satvel;
+  double dx,dy,dz,dvx,dvy,dvz,r,v,rg,vg;
+  double freq0;
+  char line[LIM],text[8];
+  struct trace *t;
+  float fmin,fmax;
+  char *env,freqlist[LIM];
+  double ra,de,azi,alt;
+  double rag,deg,azig,altg;
+
+  // Reloop stderr
+  if (freopen("/tmp/stderr.txt","w",stderr)==NULL)
+    fprintf(stderr,"Failed to redirect stderr\n");
+
+  // 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);
+
+  // Compute Graves positions
+  if (graves==1) {
+    sg=get_site(9999);
+    for (i=0;i<n;i++) 
+      obspos_xyz(mjd[i],sg.lng,sg.lat,sg.alt,&p[i].grpos,&p[i].grvel);
+  }
+
+  // Open output file
+  outfile=fopen(outfname, "w");
+
+  // Print header
+  if (graves==1)
+    fprintf(outfile, "# satno mjd r v azi alt rg vg azig altg\n");
+  else
+    fprintf(outfile, "# satno mjd r v azi alt\n");
+
+  // Loop over TLEs
+  file=fopen(tlefile,"r");
+  while (read_twoline(file,satno,&orb)==0) {
+    // Initialize
+    imode=init_sgdp4(&orb);
+    if (imode==SGDP4_ERROR)
+      printf("Error\n");
+
+    // Skip high satellites
+    if (skiphigh==1 && orb.rev<10.0)
+          continue;
+
+    // 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;
+      ra=modulo(atan2(dy,dx)*R2D,360.0);
+      de=asin(dz/r)*R2D;
+      equatorial2horizontal(mjd[i],ra,de,sg.lng,sg.lat,&azi,&alt);
+
+      // Compute Graves velocity/frequency
+      if (graves==1) {
+	dx=satpos.x-p[i].grpos.x;  
+	dy=satpos.y-p[i].grpos.y;
+	dz=satpos.z-p[i].grpos.z;
+	dvx=satvel.x-p[i].grvel.x;
+	dvy=satvel.y-p[i].grvel.y;
+	dvz=satvel.z-p[i].grvel.z;
+	rg=sqrt(dx*dx+dy*dy+dz*dz);
+	vg=(dvx*dx+dvy*dy+dvz*dz)/rg;
+	rag=modulo(atan2(dy,dx)*R2D,360.0);
+	deg=asin(dz/rg)*R2D;
+	equatorial2horizontal(mjd[i],rag,deg,sg.lng,sg.lat,&azig,&altg);
+	fprintf(outfile,"%05d %14.8lf %f %f %f %f %f %f %f %f\n",orb.satno,mjd[i],r,v,azi,alt,rg,vg,azig,altg);
+      } else {
+	fprintf(outfile,"%05d %14.8lf %f %f %f %f\n",orb.satno,mjd[i],r,v,azi,alt);
+      }
+    }
+  }
+  fclose(file);
+  fclose(outfile);
+
+  fclose(stderr);
+
+  // Free
+  free(p);
+
+  return;
+}

rftrace.h

@@ -7,3 +7,5 @@ struct trace {
 struct trace *compute_trace(char *tlefile,double *mjd,int n,int site_id,float fmin,float fmax,int *nsat,int graves);
 void identify_trace(char *tlefile,struct trace t,int satno);
 void identify_trace_graves(char *tlefile,struct trace t,int satno);
+void compute_doppler(char *tlefile,double *mjd,int n,int site_id,int satno,int graves, int skiphigh,char *outfname);
+int fgetline(FILE *file,char *s,int lim);