From 51c4fac2d9eb47d63d198db6b0cf8661a641e27e Mon Sep 17 00:00:00 2001
From: Cees Bassa <bassa@topaz.(none)>
Date: Wed, 19 Mar 2014 23:06:26 +0100
Subject: [PATCH] Misc changes

---
 data/camera.txt |   2 +
 detect.c        |  12 +++-
 imgstat.c       | 154 +++++++++++++++++++++++++++++++++++++++++++++++-
 launchtle.c     |  13 ++--
 measure.c       |   2 +-
 reduce.c        |   1 -
 tleinfo.c       |   2 +-
 7 files changed, 174 insertions(+), 12 deletions(-)

diff --git a/data/camera.txt b/data/camera.txt
index fc7491b..9874af4 100644
--- a/data/camera.txt
+++ b/data/camera.txt
@@ -2,4 +2,6 @@ W50H   7.21 5.31
 D85V   9.97 14.96
 W12H  29.9  22.8
 D85H  14.96  9.97
+D50H  25.01 16.68
+D50V  16.68 25.01
 FULL  360  240
diff --git a/detect.c b/detect.c
index c2e67b6..3a1560a 100644
--- a/detect.c
+++ b/detect.c
@@ -600,9 +600,15 @@ int main(int argc,char *argv[])
   img.nx=ff.naxis1;
   img.ny=ff.naxis2;
   img.z=(float *) malloc(img.nx*img.ny*sizeof(float));
-  for (i=0;i<ff.naxis1*ff.naxis2;i++) 
-    img.z[i]=(ff.zmax[i]-ff.zavg[i])/ff.zstd[i]>sigma;
-
+  for (i=0;i<ff.naxis1;i++) {
+    for (j=0;j<ff.naxis2;j++) {
+      k=i+ff.naxis1*j;
+      img.z[k]=(ff.zmax[k]-ff.zavg[k])/ff.zstd[k];
+      if (img.z[k]>sigma)
+	printf("%d %d %f\n",i,j,ff.dt[(int) ff.znum[k]]);
+    }
+  }
+  return 0;
   // Loop over lines
   for (k=0;;k++) {
     // Generate mask
diff --git a/imgstat.c b/imgstat.c
index 5045c96..b9f19c1 100644
--- a/imgstat.c
+++ b/imgstat.c
@@ -4,6 +4,10 @@
 #include <math.h>
 #include "qfits.h"
 
+#define LIM 128
+#define D2R M_PI/180.0
+#define R2D 180.0/M_PI
+
 struct image {
   char filename[64];
   int naxis1,naxis2;
@@ -16,17 +20,158 @@ struct image {
   char nfd[32];
   int cospar;
 };
+struct map {
+  char datadir[LIM],observer[32];
+  double lng,lat;
+  float alt;
+  int site_id;
+} m;
 struct image read_fits(char *filename);
 
+// Return x modulo y [0,y)
+double modulo(double x,double y)
+{
+  x=fmod(x,y);
+  if (x<0.0) x+=y;
+
+  return x;
+}
+
+// Greenwich Mean Sidereal Time
+double gmst(double mjd)
+{
+  double t,gmst;
+
+  t=(mjd-51544.5)/36525.0;
+
+  gmst=modulo(280.46061837+360.98564736629*(mjd-51544.5)+t*t*(0.000387933-t/38710000),360.0);
+
+  return gmst;
+}
+
+// Precess a celestial position
+void precess(double mjd0,double ra0,double de0,double mjd,double *ra,double *de)
+{
+  double t0,t;
+  double zeta,z,theta;
+  double a,b,c;
+
+  // Angles in radians
+  ra0*=D2R;
+  de0*=D2R;
+
+  // Time in centuries
+  t0=(mjd0-51544.5)/36525.0;
+  t=(mjd-mjd0)/36525.0;
+
+  // Precession angles
+  zeta=(2306.2181+1.39656*t0-0.000139*t0*t0)*t;
+  zeta+=(0.30188-0.000344*t0)*t*t+0.017998*t*t*t;
+  zeta*=D2R/3600.0;
+  z=(2306.2181+1.39656*t0-0.000139*t0*t0)*t;
+  z+=(1.09468+0.000066*t0)*t*t+0.018203*t*t*t;
+  z*=D2R/3600.0;
+  theta=(2004.3109-0.85330*t0-0.000217*t0*t0)*t;
+  theta+=-(0.42665+0.000217*t0)*t*t-0.041833*t*t*t;
+  theta*=D2R/3600.0;
+  
+  a=cos(de0)*sin(ra0+zeta);
+  b=cos(theta)*cos(de0)*cos(ra0+zeta)-sin(theta)*sin(de0);
+  c=sin(theta)*cos(de0)*cos(ra0+zeta)+cos(theta)*sin(de0);
+
+  *ra=(atan2(a,b)+z)*R2D;
+  *de=asin(c)*R2D;
+
+  if (*ra<360.0)
+    *ra+=360.0;
+  if (*ra>360.0)
+    *ra-=360.0;
+
+  return;
+}
+
+// Get observing site
+void get_site(int site_id)
+{
+  int i=0;
+  char line[LIM];
+  FILE *file;
+  int id;
+  double lat,lng;
+  float alt;
+  char abbrev[3],observer[64],filename[LIM];
+
+  sprintf(filename,"%s/data/sites.txt",m.datadir);
+  file=fopen(filename,"r");
+  if (file==NULL) {
+    printf("File with site information not found!\n");
+    return;
+  }
+  while (fgets(line,LIM,file)!=NULL) {
+    // Skip
+    if (strstr(line,"#")!=NULL)
+      continue;
+
+    // Strip newline
+    line[strlen(line)-1]='\0';
+
+    // Read data
+    sscanf(line,"%4d %2s %lf %lf %f",
+	   &id,abbrev,&lat,&lng,&alt);
+    strcpy(observer,line+38);
+
+    // Change to km
+    alt/=1000.0;
+    
+    if (id==site_id) {
+      m.lat=lat;
+      m.lng=lng;
+      m.alt=alt;
+      m.site_id=id;
+      strcpy(m.observer,observer);
+    }
+
+  }
+  fclose(file);
+
+  return;
+}
+
+// Convert equatorial into horizontal coordinates
+void equatorial2horizontal(double mjd,double ra,double de,double *azi,double *alt)
+{
+  double h;
+
+  h=gmst(mjd)+m.lng-ra;
+  
+  *azi=modulo(atan2(sin(h*D2R),cos(h*D2R)*sin(m.lat*D2R)-tan(de*D2R)*cos(m.lat*D2R))*R2D,360.0);
+  *alt=asin(sin(m.lat*D2R)*sin(de*D2R)+cos(m.lat*D2R)*cos(de*D2R)*cos(h*D2R))*R2D;
+
+  return;
+}
+
 int main(int argc,char *argv[])
 {
   int i;
   struct image img;
   float zavg,zstd,sx,sy,wx,wy;
+  double ra,de,alt,azi,mjd0=51544.5;
+  char *env;
+  
+  // Get environment variables
+  env=getenv("ST_DATADIR");
+  if (env!=NULL) {
+    strcpy(m.datadir,env);
+  } else {
+    printf("ST_DATADIR environment variable not found.\n");
+  }
 
   // Read image
   img=read_fits(argv[1]);
 
+  // Get site
+  get_site(img.cospar);
+  
   // Compute statistics
   for (i=0,zavg=0.0;i<img.naxis1*img.naxis2;i++)
     zavg+=img.zavg[i];
@@ -41,7 +186,14 @@ int main(int argc,char *argv[])
   wx=img.naxis1*sx/3600.0;
   wy=img.naxis2*sy/3600.0;
 
-  printf("%s %14.8lf %10.6f %10.6f %.2f %.2f %.2f %.2f %.2f %.2f %.1f %.1f\n",argv[1],img.mjd,img.ra0,img.de0,wx,wy,sx,sy,img.xrms,img.yrms,zavg,zstd);
+  // Precess to epoch of date
+  precess(mjd0,img.ra0,img.de0,img.mjd,&ra,&de);
+
+  // Get horizontal coordinates
+  equatorial2horizontal(img.mjd,ra,de,&azi,&alt);
+  azi=modulo(azi+180,360);
+  //  printf("%s %14.8lf %10.6f %10.6f %.2f %.2f %.2f %.2f %.2f %.2f %.1f %.1f\n",argv[1],img.mjd,img.ra0,img.de0,wx,wy,sx,sy,img.xrms,img.yrms,zavg,zstd);
+  printf("%s %14.8lf %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n",argv[1],img.mjd,img.ra0,img.de0,ra,de,azi,alt);
 
   return 0;
 }
diff --git a/launchtle.c b/launchtle.c
index 2007501..ea0761e 100644
--- a/launchtle.c
+++ b/launchtle.c
@@ -14,11 +14,6 @@
 #define CK2     5.413080e-4   /* (0.5 * XJ2 * AE * AE) */
 extern double SGDP4_jd0;
 
-void usage(void)
-{
-  return;
-}
-
 // Compute Julian Day from Date
 double date2mjd(int year,int month,double day)
 {
@@ -186,6 +181,14 @@ void format_tle(orbit_t orb,char *line1,char *line2)
   return;
 }
 
+void usage(void)
+{
+  printf("launch tle c:i:t:T:I:d:\n\n");
+  printf("-c  reference tle\n-i  reference satno\n-t  reference launch time\n-T  launch time\n-I  output satno\n-d  output desig\n");
+
+  return;
+}
+
 int main(int argc,char *argv[])
 {
   int arg=0,satno=0,satno1=84001;
diff --git a/measure.c b/measure.c
index b804a46..6b3e884 100644
--- a/measure.c
+++ b/measure.c
@@ -518,7 +518,7 @@ int main(int argc,char *argv[])
   for (;;) {
     if (redraw!=0) {
       if (redraw==1) 
-	cpgeras();
+	cpgpage();
       cpgsci(1);
       
       cpgsvp(0.1,0.9,0.1,0.85);
diff --git a/reduce.c b/reduce.c
index 47e39dd..91dd67a 100644
--- a/reduce.c
+++ b/reduce.c
@@ -1005,7 +1005,6 @@ struct image read_fits(char *filename)
   for (i=0;i<img.naxis1*img.naxis2;i++) 
     img.mask[i]=1;
 
-
   return img;
 }
 
diff --git a/tleinfo.c b/tleinfo.c
index 8d79155..69cbd67 100644
--- a/tleinfo.c
+++ b/tleinfo.c
@@ -164,7 +164,7 @@ int main(int argc,char *argv[])
   char *env;
 
   env=getenv("ST_TLEDIR");
-  sprintf(tlefile,"%s/classfd.tle",env);
+  sprintf(tlefile,"%s/bulk.tle",env);
 
   // Decode options
   while ((arg=getopt(argc,argv,"c:i:aH1ft"))!=-1) {