diff --git a/csv2tle.c b/csv2tle.c
index f949544..4414a72 100644
--- a/csv2tle.c
+++ b/csv2tle.c
@@ -169,7 +169,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/deproject.c b/deproject.c
index c808ebc..1b652e5 100644
--- a/deproject.c
+++ b/deproject.c
@@ -83,16 +83,17 @@ int main(int argc,char *argv[])
   jpg=read_jpg(jpgfile);
 
   out.nx=3000;
-  out.ny=2000;
+  out.ny=6000;
   out.nz=3;
 
+  /*
   img.x0*=4.0;
   img.y0*=4.0;
   img.a[1]/=4.0;
   img.a[2]/=4.0;
   img.b[1]/=4.0;
   img.b[2]/=4.0;
-
+  */
   out.z=(float *) malloc(sizeof(float)*out.nx*out.ny*out.nz);
 
   for (i=0;i<out.nx;i++) {
diff --git a/detect.c b/detect.c
index e0b9354..314dc2a 100644
--- a/detect.c
+++ b/detect.c
@@ -171,7 +171,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/fakeiod.c b/fakeiod.c
index 5298997..099c2f4 100644
--- a/fakeiod.c
+++ b/fakeiod.c
@@ -395,7 +395,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/faketle.c b/faketle.c
index a5cf18c..31c440d 100644
--- a/faketle.c
+++ b/faketle.c
@@ -216,7 +216,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/ferror.c b/ferror.c
index 05090bf..7d878e5 100644
--- a/ferror.c
+++ b/ferror.c
@@ -19,11 +19,11 @@ va_list arg_ptr;
     vfprintf(stderr, format, arg_ptr);
     va_end(arg_ptr);
 
-    fprintf(stderr, "\nNow terminating the program...\n");
-    fflush(stderr);
-
-    exit(5);
+    //fprintf(stderr, "\nNow terminating the program...\n");
+    //    fflush(stderr);
 
+    //    exit(5);
+    return;
 }
 
 /* ===================================================================== */
diff --git a/jpg2fits.c b/jpg2fits.c
index d843261..ef58436 100644
--- a/jpg2fits.c
+++ b/jpg2fits.c
@@ -428,7 +428,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/launchtle.c b/launchtle.c
index ea0761e..b230031 100644
--- a/launchtle.c
+++ b/launchtle.c
@@ -30,7 +30,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/measure.c b/measure.c
index 3369d26..a90be35 100644
--- a/measure.c
+++ b/measure.c
@@ -236,7 +236,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/normal.c b/normal.c
index 3481ce8..31814ed 100644
--- a/normal.c
+++ b/normal.c
@@ -227,6 +227,10 @@ orbit_t rv2el(int satno,double mjd,xyz_t r0,xyz_t v0)
       orb[i+1].eqinc=0.0;
   }
 
+  orb[i].mnan=modulo(orb[i].mnan,2.0*M_PI);
+  orb[i].ascn=modulo(orb[i].ascn,2.0*M_PI);
+  orb[i].argp=modulo(orb[i].argp,2.0*M_PI);
+
   return orb[i];
 }
 
@@ -296,7 +300,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/pass.c b/pass.c
index ebff9bb..839926b 100644
--- a/pass.c
+++ b/pass.c
@@ -468,7 +468,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/pgm2fits.c b/pgm2fits.c
index 6f8c135..dd9f7f0 100644
--- a/pgm2fits.c
+++ b/pgm2fits.c
@@ -507,7 +507,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/posmatch.c b/posmatch.c
index af49667..ef8e1b5 100644
--- a/posmatch.c
+++ b/posmatch.c
@@ -51,7 +51,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/posvel.c b/posvel.c
index a94c5d4..995c5cc 100644
--- a/posvel.c
+++ b/posvel.c
@@ -90,7 +90,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/propagate.c b/propagate.c
index e080834..20eb2eb 100644
--- a/propagate.c
+++ b/propagate.c
@@ -217,6 +217,10 @@ orbit_t rv2el(int satno,double mjd,xyz_t r0,xyz_t v0)
       orb[i+1].eqinc=0.0;
   }
 
+  orb[i].mnan=modulo(orb[i].mnan,2.0*M_PI);
+  orb[i].ascn=modulo(orb[i].ascn,2.0*M_PI);
+  orb[i].argp=modulo(orb[i].argp,2.0*M_PI);
+
   return orb[i];
 }
 
@@ -286,7 +290,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/reduce.c b/reduce.c
index 91dd67a..f594068 100644
--- a/reduce.c
+++ b/reduce.c
@@ -1096,7 +1096,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/residuals.c b/residuals.c
index 374b855..fcbdc3e 100644
--- a/residuals.c
+++ b/residuals.c
@@ -510,7 +510,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/runsched.c b/runsched.c
index 006d34f..089a95b 100644
--- a/runsched.c
+++ b/runsched.c
@@ -129,14 +129,14 @@ void send_position(char *sra,char *sde)
   close(skt); 
  
   // Set restart
-  file=fopen("/data/satobs/control/state.txt","w");
+  file=fopen("/media/storage/satobs/control/state.txt","w");
   if (file!=NULL) {
     fprintf(file,"restart");
     fclose(file);
   }
 
   // Set position
-  file=fopen("/data/satobs/control/position.txt","w");
+  file=fopen("/media/storage/satobs/control/position.txt","w");
   if (file!=NULL) {
     fprintf(file,"%s %s\n",sra,sde);
     fclose(file);
diff --git a/rv2tle.c b/rv2tle.c
index 3e60199..5be1d41 100644
--- a/rv2tle.c
+++ b/rv2tle.c
@@ -224,11 +224,12 @@ orbit_t rv2el(int satno,double mjd,xyz_t r0,xyz_t v0)
       orb[i+1].ecc=0.0;
     if (orb[i+1].eqinc<0.0)
       orb[i+1].eqinc=0.0;
-    orb[i+1].mnan=modulo(orb[i+1].mnan,2.0*M_PI);
-    orb[i+1].ascn=modulo(orb[i+1].ascn,2.0*M_PI);
-    orb[i+1].argp=modulo(orb[i+1].argp,2.0*M_PI);
   }
 
+  orb[i].mnan=modulo(orb[i].mnan,2.0*M_PI);
+  orb[i].ascn=modulo(orb[i].ascn,2.0*M_PI);
+  orb[i].argp=modulo(orb[i].argp,2.0*M_PI);
+
   return orb[i];
 }
 
@@ -298,7 +299,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/satfit.c b/satfit.c
index 1965f41..bc80e3f 100644
--- a/satfit.c
+++ b/satfit.c
@@ -1168,7 +1168,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
@@ -1453,15 +1453,7 @@ double chisq(double a[])
   }
   if (nsel>0) 
     rms=sqrt(rms/(float) nsel);
-  /*
-  // Perigee match
-  satpos_xyz(57022.84238426+2400000.5,&satpos,&satvel);
-  dx=-1557.4323-satpos.x;
-  dy=+3594.7374-satpos.y;
-  dz=-7566.5506-satpos.z;
-  chisq+=1000*(dx*dx+dy*dy+dz*dz);
-  printf("%f %f %f\n",dx,dy,dz);
-  */
+  
   d.chisq=chisq;
   d.rms=rms;
   d.nsel=nsel;
@@ -1888,5 +1880,9 @@ orbit_t rv2el(int satno,double mjd,xyz_t r0,xyz_t v0)
       orb[i+1].eqinc=0.0;
   }
 
+  orb[i].mnan=modulo(orb[i].mnan,2.0*M_PI);
+  orb[i].ascn=modulo(orb[i].ascn,2.0*M_PI);
+  orb[i].argp=modulo(orb[i].argp,2.0*M_PI);
+
   return orb[i];
 }
diff --git a/satmap.c b/satmap.c
index 9a19a38..d66f520 100644
--- a/satmap.c
+++ b/satmap.c
@@ -766,7 +766,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/satorbit.c b/satorbit.c
index a126bb9..140e0fb 100644
--- a/satorbit.c
+++ b/satorbit.c
@@ -301,7 +301,7 @@ void plot_track(void)
     
     Isat=orb.satno;
     imode=init_sgdp4(&orb);
-    
+
     if(imode == SGDP4_ERROR) continue;
 
     jd=m.mjd+2400000.5;
@@ -1371,7 +1371,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/skymap.c b/skymap.c
index 73f758c..bb4d600 100644
--- a/skymap.c
+++ b/skymap.c
@@ -68,6 +68,14 @@ struct observation {
   float dt,st,dr,sr,dx,dy,t;
   int flag;
 };
+struct coeff_lr {
+  int nd,nm,nm1,nf;
+  double sa,ca;
+} clr[60];
+struct coeff_b {
+  int nd,nm,nm1,nf;
+  double sa;
+} cb[60];
 int fgetline(FILE *,char *,int);
 double modulo(double,double);
 void reverse(double,double,double *,double *);
@@ -103,6 +111,10 @@ double doy2mjd(int year,double doy);
 struct observation decode_iod_observation(char *iod_line);
 void plot_iod(char *filename);
 void get_site(int site_id);
+void lunpos_xyz(double mjd,xyz_t *pos,double *ra,double *de);
+void ecliptical2equatorial(double l,double b,double *ra,double *de);
+void skymap_plotsun(void);
+void skymap_plotmoon(void);
 
 void usage()
 {
@@ -131,7 +143,8 @@ void usage()
 void init_skymap(void)
 {
   int i;
-  char *env;
+  char *env,filename[128];
+  FILE *file;
 
   // Default Map parameters
   m.azi0=0;
@@ -191,6 +204,20 @@ void init_skymap(void)
   }
   sprintf(m.tlefile,"%s/classfd.tle",m.tledir);
 
+  // Read LR coefficients
+  sprintf(filename,"%s/data/moonLR.dat",m.datadir);
+  file=fopen(filename,"r");
+  for (i=0;i<60;i++)
+    fscanf(file,"%d %d %d %d %lf %lf",&clr[i].nd,&clr[i].nm,&clr[i].nm1,&clr[i].nf,&clr[i].sa,&clr[i].ca);
+  fclose(file);
+
+  // Read B coefficients
+  sprintf(filename,"%s/data/moonB.dat",m.datadir);
+  file=fopen(filename,"r");
+  for (i=0;i<60;i++)
+    fscanf(file,"%d %d %d %d %lf",&cb[i].nd,&cb[i].nm,&cb[i].nm1,&cb[i].nf,&cb[i].sa);
+  fclose(file);
+
   return;
 }
 
@@ -1606,7 +1633,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
@@ -2018,6 +2045,30 @@ double parallactic_angle(double mjd,double ra,double de)
   return q;
 }
 
+void skymap_plotmoon(void)
+{
+  double rx,ry,ra,de,azi,alt;
+  xyz_t lunpos;
+
+  // Get moon position
+  lunpos_xyz(m.mjd,&lunpos,&ra,&de);
+  equatorial2horizontal(m.mjd,ra,de,&azi,&alt);
+
+  if (strcmp(m.orientation,"horizontal")==0) 
+    forward(azi,alt,&rx,&ry);
+  else if (strcmp(m.orientation,"equatorial")==0) 
+    forward(ra,de,&rx,&ry);
+
+  cpgsci(14);
+  if (m.w>60.0)
+    cpgcirc((float) rx,(float) ry,2.0);
+  else
+    cpgcirc((float) rx,(float) ry,0.25);
+  cpgsci(1);
+
+  return;
+}
+
 void skymap_plotsun(void)
 {
   double rx,ry;
@@ -2205,6 +2256,7 @@ int plot_skymap(void)
 	equatorial2horizontal(m.mjd,m.ra0,m.de0,&m.azi0,&m.alt0);
       }
       skymap_plotsun();
+      skymap_plotmoon();
 
       plot_apex(35786.0,0.0);
       plot_apex(39035,63.4);
@@ -2885,5 +2937,94 @@ void plot_iod(char *filename)
   cpgsci(1);
   
 
+  return;
+}
+
+// Moon position
+void lunpos_xyz(double mjd,xyz_t *pos,double *ra,double *de)
+{
+  int i;
+  double t,t2,t3,t4;
+  double l1,d,m,m1,f,a1,a2,a3,e,ef;
+  double suml,sumb,sumr,arglr,argb;
+  double l,b,r;
+
+  // Julian Centuries
+  t=(mjd-51544.5)/36525.0;
+
+  // Powers of t
+  t2=t*t;
+  t3=t2*t;
+  t4=t3*t;
+
+  // angles
+  l1=modulo(218.3164477+481267.88123421*t-0.0015786*t2+t3/538841.0-t4/65194000.0,360.0)*D2R;
+  d=modulo(297.8501921+445267.1114034*t-0.0018819*t2+t3/545868.0-t4/113065000.0,360.0)*D2R;
+  m=modulo(357.5291092+35999.0502909*t-0.0001536*t2+t3/24490000.0,360.0)*D2R;
+  m1=modulo(134.9633964+477198.8675055*t+0.0087417*t2+t3/69699.0-t4/14712000.0,360.0)*D2R;
+  f=modulo(93.2720950+483202.0175233*t-0.0036539*t2-t3/3526000.0+t4/86331000.0,360.0)*D2R;
+  a1=modulo(119.75+131.849*t,360.0)*D2R;
+  a2=modulo(53.09+479264.290*t,360.0)*D2R;
+  a3=modulo(313.45+481266.484*t,360.0)*D2R;
+  e=1.0-0.002516*t-0.0000074*t2;
+
+  // Compute sums
+  for (i=0,suml=sumb=sumr=0.0;i<60;i++) {
+    // Arguments
+    arglr=clr[i].nd*d+clr[i].nm*m+clr[i].nm1*m1+clr[i].nf*f;
+    argb=cb[i].nd*d+cb[i].nm*m+cb[i].nm1*m1+cb[i].nf*f;
+    
+    // E multiplication factor
+    if (abs(clr[i].nm)==1)
+      ef=e;
+    else if (abs(clr[i].nm)==2)
+      ef=e*e;
+    else
+      ef=1.0;
+
+    // Sums
+    suml+=clr[i].sa*sin(arglr)*ef;
+    sumr+=clr[i].ca*cos(arglr)*ef;
+
+    // E multiplication factor
+    if (abs(cb[i].nm)==1)
+      ef=e;
+    else if (abs(cb[i].nm)==2)
+      ef=e*e;
+    else
+      ef=1.0;
+
+    // Sums
+    sumb+=cb[i].sa*sin(argb)*ef;
+  }
+
+  // Additives
+  suml+=3958*sin(a1)+1962*sin(l1-f)+318*sin(a2);
+  sumb+=-2235*sin(l1)+382*sin(a3)+175*sin(a1-f)+175*sin(a1+f)+127*sin(l1-m1)-115*sin(l1+m1);
+
+  // Ecliptic longitude, latitude and distance
+  l=modulo(l1*R2D+suml/1000000.0,360.0);
+  b=sumb/1000000.0;
+  r=385000.56+sumr/1000.0;
+
+  // Equatorial
+  ecliptical2equatorial(l,b,ra,de);
+
+  // Position
+  pos->x=r*cos(*de*D2R)*cos(*ra*D2R);
+  pos->y=r*cos(*de*D2R)*sin(*ra*D2R);
+  pos->z=r*sin(*de*D2R);
+
+  return;
+}
+
+// Convert ecliptical into equatorial coordinates
+void ecliptical2equatorial(double l,double b,double *ra,double *de)
+{
+  double eps=23.4392911;
+
+  *ra=modulo(atan2(sin(l*D2R)*cos(eps*D2R)-tan(b*D2R)*sin(eps*D2R),cos(l*D2R))*R2D,360.0);
+  *de=asin(sin(b*D2R)*cos(eps*D2R)+cos(b*D2R)*sin(eps*D2R)*sin(l*D2R))*R2D;
+
   return;
 }
diff --git a/slewto.c b/slewto.c
index 0904ff8..ce60b46 100644
--- a/slewto.c
+++ b/slewto.c
@@ -271,7 +271,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/tle2rv.c b/tle2rv.c
index 65217b2..21a81c1 100644
--- a/tle2rv.c
+++ b/tle2rv.c
@@ -171,7 +171,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/tleinfo.c b/tleinfo.c
index e938166..f2594f8 100644
--- a/tleinfo.c
+++ b/tleinfo.c
@@ -35,7 +35,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
diff --git a/vadd.c b/vadd.c
new file mode 100644
index 0000000..64bd3d9
--- /dev/null
+++ b/vadd.c
@@ -0,0 +1,453 @@
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include "cpgplot.h"
+#include "sgdp4h.h"
+#include "satutl.h"
+#include <getopt.h>
+
+#define LIM 128
+#define XKE 0.07436680 // Guassian Gravitational Constant
+#define XKMPER 6378.135
+#define AE 1.0
+#define XMNPDA 1440.0
+#define R2D 180.0/M_PI
+#define D2R M_PI/180.0
+
+extern double SGDP4_jd0;
+
+// Dot product
+float dot(xyz_t a,xyz_t b)
+{
+  return a.x*b.x+a.y*b.y+a.z*b.z;
+}
+
+// Return x modulo y [0,y)
+double modulo(double x,double y)
+{
+  x=fmod(x,y);
+  if (x<0.0) x+=y;
+
+  return x;
+}
+
+// Magnitude
+double magnitude(xyz_t a)
+{
+  return sqrt(dot(a,a));
+}
+
+// Cross product
+xyz_t cross(xyz_t a,xyz_t b)
+{
+  xyz_t c;
+
+  c.x=a.y*b.z-a.z*b.y;
+  c.y=a.z*b.x-a.x*b.z;
+  c.z=a.x*b.y-a.y*b.x;
+
+  return c;
+}
+
+// Compute Date from Julian Day
+void mjd2date(double mjd,int *year,int *month,double *day)
+{
+  double f,jd;
+  int z,alpha,a,b,c,d,e;
+
+  jd=mjd+2400000.5;
+  jd+=0.5;
+
+  z=floor(jd);
+  f=fmod(jd,1.);
+
+  if (z<2299161)
+    a=z;
+  else {
+    alpha=floor((z-1867216.25)/36524.25);
+    a=z+1+alpha-floor(alpha/4.);
+  }
+  b=a+1524;
+  c=floor((b-122.1)/365.25);
+  d=floor(365.25*c);
+  e=floor((b-d)/30.6001);
+
+  *day=b-d-floor(30.6001*e)+f;
+  if (e<14)
+    *month=e-1;
+  else
+    *month=e-13;
+
+  if (*month>2)
+    *year=c-4716;
+  else
+    *year=c-4715;
+
+  return;
+}
+
+// MJD to DOY
+double mjd2doy(double mjd,int *yr)
+{
+  int year,month,k=2;
+  double day,doy;
+  
+  mjd2date(mjd,&year,&month,&day);
+
+  if (year%4==0 && year%400!=0)
+    k=1;
+
+  doy=floor(275.0*month/9.0)-k*floor((month+9.0)/12.0)+day-30;
+
+  *yr=year;
+
+  return doy;
+}
+
+// Clasical elements
+orbit_t classel(int ep_year,double ep_day,xyz_t r,xyz_t v)
+{
+  int i;
+  double rm,vm,vm2,rvm,mu=1.0;;
+  double chi,xp,yp,sx,cx,b,ee;
+  double a,ecc,incl,node,peri,mm,n;
+  xyz_t h,e,kk,nn;
+  orbit_t orb;
+
+  r.x/=XKMPER;
+  r.y/=XKMPER;
+  r.z/=XKMPER;
+  v.x/=(XKE*XKMPER/AE*XMNPDA/86400.0);
+  v.y/=(XKE*XKMPER/AE*XMNPDA/86400.0);
+  v.z/=(XKE*XKMPER/AE*XMNPDA/86400.0);
+
+  rm=magnitude(r);
+  vm2=dot(v,v);
+  rvm=dot(r,v);
+  h=cross(r,v);
+  chi=dot(h,h)/mu;
+
+  e.x=(vm2/mu-1.0/rm)*r.x-rvm/mu*v.x;
+  e.y=(vm2/mu-1.0/rm)*r.y-rvm/mu*v.y;
+  e.z=(vm2/mu-1.0/rm)*r.z-rvm/mu*v.z;
+
+  a=pow(2.0/rm-vm2/mu,-1);
+  ecc=magnitude(e);
+  incl=acos(h.z/magnitude(h))*R2D;
+  
+  kk.x=0.0;
+  kk.y=0.0;
+  kk.z=1.0;
+  nn=cross(kk,h);
+  if (nn.x==0.0 && nn.y==0.0)
+    nn.x=1.0;
+  node=atan2(nn.y,nn.x)*R2D;
+  if (node<0.0)
+    node+=360.0;
+
+  peri=acos(dot(nn,e)/(magnitude(nn)*ecc))*R2D;
+  if (e.z<0.0)
+    peri=360.0-peri;
+  if (peri<0.0)
+    peri+=360.0;
+
+  // Elliptic motion
+  if (ecc<1.0) {
+    xp=(chi-rm)/ecc;
+    yp=rvm/ecc*sqrt(chi/mu);
+    b=a*sqrt(1.0-ecc*ecc);
+    cx=xp/a+ecc;
+    sx=yp/b;
+    ee=atan2(sx,cx);
+    n=XKE*sqrt(mu/(a*a*a));
+    mm=(ee-ecc*sx)*R2D;
+  } 
+  if (mm<0.0)
+    mm+=360.0;
+
+  // Fill
+  orb.satno=0;
+  orb.eqinc=incl*D2R;
+  orb.ascn=node*D2R;
+  orb.argp=peri*D2R;
+  orb.mnan=mm*D2R;
+  orb.ecc=ecc;
+  orb.rev=XKE*pow(a,-3.0/2.0)*XMNPDA/(2.0*M_PI);
+  orb.bstar=0.0;
+  orb.ep_year=ep_year;
+  orb.ep_day=ep_day;
+  orb.norb=0;
+
+  return orb;
+}
+
+orbit_t rv2el(int satno,double mjd,xyz_t r0,xyz_t v0)
+{
+  int i,imode;
+  orbit_t orb[5],orb1[5];
+  xyz_t r,v;
+  kep_t kep;
+  char line1[70],line2[70];
+  int ep_year;
+  double ep_day;
+
+  // Epoch
+  ep_day=mjd2doy(mjd,&ep_year);
+
+  // Initial guess
+  orb[0]=classel(ep_year,ep_day,r0,v0);
+  orb[0].satno=satno;
+  
+  for (i=0;i<4;i++) {
+    // Propagate
+    imode=init_sgdp4(&orb[i]);
+    imode=satpos_xyz(mjd+2400000.5,&r,&v);
+
+    // Compute initial orbital elements
+    orb1[i]=classel(ep_year,ep_day,r,v);
+
+    // Adjust
+    orb[i+1].rev=orb[i].rev+orb[0].rev-orb1[i].rev;
+    orb[i+1].ascn=orb[i].ascn+orb[0].ascn-orb1[i].ascn;
+    orb[i+1].argp=orb[i].argp+orb[0].argp-orb1[i].argp;
+    orb[i+1].mnan=orb[i].mnan+orb[0].mnan-orb1[i].mnan;
+    orb[i+1].eqinc=orb[i].eqinc+orb[0].eqinc-orb1[i].eqinc;
+    orb[i+1].ecc=orb[i].ecc+orb[0].ecc-orb1[i].ecc;
+    orb[i+1].ep_year=orb[i].ep_year;
+    orb[i+1].ep_day=orb[i].ep_day;
+    orb[i+1].satno=orb[i].satno;
+    orb[i+1].norb=orb[i].norb;
+    orb[i+1].bstar=orb[i].bstar;
+
+    // Keep in range
+    if (orb[i+1].ecc<0.0)
+      orb[i+1].ecc=0.0;
+    if (orb[i+1].eqinc<0.0)
+      orb[i+1].eqinc=0.0;
+  }
+
+  orb[i].mnan=modulo(orb[i].mnan,2.0*M_PI);
+  orb[i].ascn=modulo(orb[i].ascn,2.0*M_PI);
+  orb[i].argp=modulo(orb[i].argp,2.0*M_PI);
+
+  return orb[i];
+}
+
+// Format TLE
+void format_tle(orbit_t orb,char *line1,char *line2)
+{
+  int i,csum;
+  char sbstar[]=" 00000-0",bstar[13];
+
+  // Format Bstar term
+  if (fabs(orb.bstar)>1e-9) {
+    sprintf(bstar,"%11.4e",10*orb.bstar);
+    sbstar[0] = bstar[0];  sbstar[1] = bstar[1];  sbstar[2] = bstar[3];  sbstar[3] = bstar[4];
+    sbstar[4] = bstar[5];  sbstar[5] = bstar[6];  sbstar[6] = bstar[8];  sbstar[7] = bstar[10];  sbstar[8] = '\0';
+  }
+  // Print lines
+  sprintf(line1,"1 %05dU %-8s %2d%012.8f  .00000000  00000-0 %8s 0    0",orb.satno,orb.desig,orb.ep_year-2000,orb.ep_day,sbstar);
+  sprintf(line2,"2 %05d %8.4f %8.4f %07.0f %8.4f %8.4f %11.8f    0",orb.satno,DEG(orb.eqinc),DEG(orb.ascn),1E7*orb.ecc,DEG(orb.argp),DEG(orb.mnan),orb.rev);
+
+  // Compute checksums
+  for (i=0,csum=0;i<strlen(line1);i++) {
+    if (isdigit(line1[i]))
+      csum+=line1[i]-'0';
+    else if (line1[i]=='-')
+      csum++;
+  }
+  sprintf(line1,"%s%d",line1,csum%10);
+  for (i=0,csum=0;i<strlen(line2);i++) {
+    if (isdigit(line2[i]))
+      csum+=line2[i]-'0';
+    else if (line2[i]=='-')
+      csum++;
+  }
+  sprintf(line2,"%s%d",line2,csum%10);
+
+  return;
+}
+
+// Present nfd
+void nfd_now(char *s)
+{
+  time_t rawtime;
+  struct tm *ptm;
+
+  // Get UTC time
+  time(&rawtime);
+  ptm=gmtime(&rawtime);
+    
+  sprintf(s,"%04d-%02d-%02dT%02d:%02d:%02d",ptm->tm_year+1900,ptm->tm_mon+1,ptm->tm_mday,ptm->tm_hour,ptm->tm_min,ptm->tm_sec);
+  
+  return;
+}
+
+// Compute Julian Day from Date
+double date2mjd(int year,int month,double day)
+{
+  int a,b;
+  double jd;
+
+  if (month<3) {
+    year--;
+    month+=12;
+  }
+
+  a=floor(year/100.);
+  b=2.-a+floor(a/4.);
+
+  if (year<1582) b=0;
+  if (year==1582 && month<10) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
+
+  jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
+
+  return jd-2400000.5;
+}
+
+// nfd2mjd
+double nfd2mjd(char *date)
+{
+  int year,month,day,hour,min,sec;
+  double mjd,dday;
+
+  sscanf(date,"%04d-%02d-%02dT%02d:%02d:%02d",&year,&month,&day,&hour,&min,&sec);
+  dday=day+hour/24.0+min/1440.0+sec/86400.0;
+
+  mjd=date2mjd(year,month,dday);
+
+  return mjd;
+}
+
+void usage(void)
+{
+  printf("propagate c:i:t:m:\n\nPropagates orbital elements to a new epoch using the SGP4/SDP4 model.\nDefault operation propagates classfd.tle to now,\n\n-c  input catalog\n-i  Satellite number\n-t  New epoch (YYYY-MM-DDTHH:MM:SS)\n-m  New epoch (MJD)\n");
+
+  return;
+}
+
+
+int main(int argc,char *argv[])
+{
+  int imode,satno=0,arg,satnomin,flag=0,satnonew=-1;
+  FILE *file;
+  orbit_t orb;
+  xyz_t r,v,n,dv;
+  char tlefile[LIM],nfd[32];
+  char line1[80],line2[80],desig[20];
+  double mjd,ra,de,dr,drmin;
+  float vadd=0.0;
+  char direction[16]="radial";
+  char *env;
+
+  // Get environment variable
+  env=getenv("ST_TLEDIR");
+  sprintf(tlefile,"%s/classfd.tle",env);
+
+  // Set date
+  nfd_now(nfd);
+  mjd=nfd2mjd(nfd);
+
+  // Decode options
+  while ((arg=getopt(argc,argv,"c:i:t:m:hv:d:I:"))!=-1) {
+    switch (arg) {
+
+    case 't':
+      strcpy(nfd,optarg);
+      mjd=nfd2mjd(nfd);
+      break;
+
+    case 'm':
+      mjd=(double) atof(optarg);
+      break;
+      
+    case 'c':
+      strcpy(tlefile,optarg);
+      break;
+
+    case 'i':
+      satno=atoi(optarg);
+      break;
+
+    case 'h':
+      usage();
+      return 0;
+      break;
+
+    case 'v':
+      vadd=atof(optarg);
+      break;
+
+    case 'd':
+      strcpy(direction,optarg);
+      break;
+
+    case 'I':
+      satnonew=atoi(optarg);
+      break;
+
+    default:
+      usage();
+      return 0;
+    }
+  }
+
+  // Reloop stderr
+  freopen("/tmp/stderr.txt","w",stderr);
+
+  // Open file
+  file=fopen(tlefile,"r");
+  while (read_twoline(file,satno,&orb)==0) {
+    format_tle(orb,line1,line2);
+    strcpy(desig,orb.desig);
+
+    // Propagate
+    imode=init_sgdp4(&orb);
+    imode=satpos_xyz(mjd+2400000.5,&r,&v);
+    
+    // Compute normal
+    n=cross(r,v);
+
+    // Add velocity
+    if (strcmp(direction,"prograde")==0) {
+      dv.x=vadd*v.x/magnitude(v);
+      dv.y=vadd*v.y/magnitude(v);
+      dv.z=vadd*v.z/magnitude(v);
+    } else if (strcmp(direction,"radial")==0) {
+      dv.x=vadd*r.x/magnitude(r);
+      dv.y=vadd*r.y/magnitude(r);
+      dv.z=vadd*r.z/magnitude(r);
+    } else if (strcmp(direction,"normal")==0) {
+      dv.x=vadd*n.x/magnitude(n);
+      dv.y=vadd*n.y/magnitude(n);
+      dv.z=vadd*n.z/magnitude(n);
+    } else {
+      dv.x=0.0;
+      dv.y=0.0;
+      dv.z=0.0;
+    }
+    v.x+=dv.x/1000.0;
+    v.y+=dv.y/1000.0;
+    v.z+=dv.z/1000.0;
+
+    // Convert
+    orb=rv2el(orb.satno,mjd,r,v);
+
+    if (satnonew==-1) {
+      strcpy(orb.desig,desig);
+    } else {
+      strcpy(orb.desig,"15999A");
+      orb.satno=satnonew;
+    }
+
+    // Print tle
+    format_tle(orb,line1,line2);
+    printf("%s\n%s\n# %05d + %g m/s %s\n",line1,line2,satno,vadd,direction);
+  }
+  fclose(file);
+
+
+  return 0;
+}
diff --git a/xyz2tle.c b/xyz2tle.c
index 66a2e50..45d3504 100644
--- a/xyz2tle.c
+++ b/xyz2tle.c
@@ -77,7 +77,7 @@ double date2mjd(int year,int month,double day)
 
   if (year<1582) b=0;
   if (year==1582 && month<10) b=0;
-  if (year==1852 && month==10 && day<=4) b=0;
+  if (year==1582 && month==10 && day<=4) b=0;
 
   jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
 
@@ -509,6 +509,10 @@ orbit_t rv2el(int satno,double mjd,xyz_t r0,xyz_t v0)
       orb[i+1].eqinc=0.0;
   }
 
+  orb[i].mnan=modulo(orb[i].mnan,2.0*M_PI);
+  orb[i].ascn=modulo(orb[i].ascn,2.0*M_PI);
+  orb[i].argp=modulo(orb[i].argp,2.0*M_PI);
+
   return orb[i];
 }