Formatting IOD output

.gitignore

@@ -1,6 +1,9 @@
 # Object files
 *.o
 
+# Python compiled objects
+*.pyc
+
 # Emacs backup files
 *~
 

python/extract_tracks.py

@@ -21,13 +21,30 @@ def pos2rel(ff,x,y):
 
     return world[0,0],world[0,1]
 
-def dec2sex(angle):
-    if angle<0.0:
+# IOD position format 2: RA/DEC = HHMMmmm+DDMMmm MX   (MX in minutes of arc)
+def format_position(ra,de):
+    ram=60.0*ra/15.0
+    rah=int(np.floor(ram/60.0))
+    ram-=60.0*rah
+    
+    des=np.sign(de)
+    dem=60.0*np.abs(de)
+    ded=int(np.floor(dem/60.0))
+    dem-=60.0*ded
+
+    if des==-1:
         sign="-"
     else:
         sign="+"
-    angle=60.0*np.abs(angle)
-    
+
+    return ("%02d%06.3f%c%02d%05.2f"%(rah,ram,sign,ded,dem)).replace(".","")
+
+# Format IOD line
+def format_iod_line(norad,cospar,site_id,t,ra,de):
+    pstr=format_position(ra,de)
+    tstr=t.replace("-","").replace("T","").replace(":","").replace(".","")
+
+    return "%05d %-9s %04d G %s 17 25 %s 37 S"%(norad,cospar,site_id,tstr,pstr)
 
 # Settings
 
@@ -78,10 +95,12 @@ for line in lines:
     x,y,t,sig=ff.significant(trksig,id.x0,id.y0,id.dxdt,id.dydt,trkrmin)
 
     # Fit tracks
-    if len(x)>ntrkmin:
+    if len(t)>ntrkmin:
         obs=observation(ff,x,y,t,sig)
 
-        obs.ra,obs.de=pos2rel(ff,obs.x0,obs.y0)
+        print("%s"%format_iod_line(id.norad,"18 555A",ff.site_id,obs.nfd,obs.ra,obs.de))
+
+        
         
         # Plot
         ppgplot.pgopen("/xs")
@@ -112,12 +131,12 @@ for line in lines:
         
         ppgplot.pgend()
 
-        plt.figure()
-        plt.plot(t,x,'.')
-        plt.plot(t,y,'.')
-        plt.plot([obs.tmin,obs.tmax],[obs.xmin,obs.xmax])
-        plt.plot([obs.tmin,obs.tmax],[obs.ymin,obs.ymax])
-        plt.show()
+#        plt.figure()
+#        plt.plot(t,x,'.')
+#        plt.plot(t,y,'.')
+#        plt.plot([obs.tmin,obs.tmax],[obs.xmin,obs.xmax])
+#        plt.plot([obs.tmin,obs.tmax],[obs.ymin,obs.ymax])
+#        plt.show()
 
         # Track and stack
     else:

python/stio.py

@@ -3,6 +3,7 @@
 import numpy as np
 from astropy.io import fits
 from astropy.time import Time
+from astropy import wcs
 
 class observation:
     """Satellite observation"""
@@ -30,6 +31,11 @@ class observation:
         self.ymin=self.y0+self.dydt*(self.tmin-self.tmid)
         self.xmax=self.x0+self.dxdt*(self.tmax-self.tmid)
         self.ymax=self.y0+self.dydt*(self.tmax-self.tmid)
+
+        # Compute ra/dec
+        world=ff.w.wcs_pix2world(np.array([[self.x0,self.y0]]),1)
+        self.ra=world[0,0]
+        self.de=world[0,1]
         
 class satid:
     """Satellite identifications"""
@@ -109,6 +115,7 @@ class fourframe:
             self.cunit=[hdu[0].header['CUNIT1'],hdu[0].header['CUNIT2']]
             self.crres=np.array([hdu[0].header['CRRES1'],hdu[0].header['CRRES2']])
 
+        # Compute image properties
         self.sx=np.sqrt(self.cd[0,0]**2+self.cd[1,0]**2)
         self.sy=np.sqrt(self.cd[0,1]**2+self.cd[1,1]**2)
         self.wx=self.nx*self.sx
@@ -116,6 +123,14 @@ class fourframe:
         self.vmin=np.mean(self.zmax)-2.0*np.std(self.zmax)
         self.vmax=np.mean(self.zmax)+6.0*np.std(self.zmax)
 
+        # Setup WCS
+        self.w=wcs.WCS(naxis=2)
+        self.w.wcs.crpix=self.crpix
+        self.w.wcs.crval=self.crval
+        self.w.wcs.cd=self.cd
+        self.w.wcs.ctype=self.ctype
+        self.w.wcs.set_pv([(2,1,45.0)])
+        
 
     def significant(self,sigma,x0,y0,dxdt,dydt,rmin=10.0):
         """Extract significant points"""