diff --git a/python/calibrate.py b/python/calibrate.py
new file mode 100644
index 0000000..3125adf
--- /dev/null
+++ b/python/calibrate.py
@@ -0,0 +1,142 @@
+#!/usr/bin/env python
+
+import os
+from astropy.io import fits
+import matplotlib.pyplot as plt
+import numpy as np
+from astropy import wcs
+from scipy import optimize
+
+def residual(a,x,y,z):
+    return z-(a[0]+a[1]*x+a[2]*y)
+
+def offsets(ra0,de0,ra,de):
+    w=wcs.WCS(naxis=2)
+    w.wcs.crpix=np.array([0.0,0.0])
+    w.wcs.crval=np.array([ra0,de0])
+    w.wcs.cd=np.array([[1.0/3600.0,0.0],[0.0,1.0/3600.0]])
+    w.wcs.ctype=["RA---TAN","DEC--TAN"]
+    w.wcs.set_pv([(2,1,45.0)])
+
+    world=np.stack((ra,de),axis=-1)
+    pix=w.wcs_world2pix(world,1)
+    return pix[:,0],pix[:,1]
+    
+
+def read_pixel_catalog(fname):
+    d=np.loadtxt(fname)
+    x,y,mag=d[:,0],d[:,1],d[:,2]
+    
+    return x,y,mag
+
+def read_tycho2_catalog(maxmag=9.0):
+    hdu=fits.open("/home/bassa/code/python/satellite/astrometry-tycho2/build/tyc2.fits")
+
+    ra=hdu[1].data.field('RA')
+    de=hdu[1].data.field('DEC')
+    mag=hdu[1].data.field('MAG_VT')
+
+    c=mag<maxmag
+    print("Selected %d stars from Tycho-2 catalog with m<%.1f"%(np.sum(c),maxmag))
+    return ra[c],de[c],mag[c]
+
+def match_catalogs(ra,de,x,y,xs,ys,rmin):
+    res=[]
+    for i in range(len(x)):
+        dx,dy=xs-x[i],ys-y[i]
+        r=np.sqrt(dx*dx+dy*dy)
+        if np.min(r)<rmin:
+            j=np.argmin(r)
+            res.append([ra[i],de[i],xs[j],ys[j]])
+            
+    return np.array(res)
+
+def fit_wcs(res,x0,y0):
+    ra,de,x,y=res[:,0],res[:,1],res[:,2],res[:,3]
+    ra0,de0=np.mean(ra),np.mean(de)
+    dx,dy=x-x0,y-y0
+
+    for i in range(5):
+        w=wcs.WCS(naxis=2)
+        w.wcs.crpix=np.array([0.0,0.0])
+        w.wcs.crval=np.array([ra0,de0])
+        w.wcs.cd=np.array([[1.0,0.0],[0.0,1.0]])
+        w.wcs.ctype=["RA---TAN","DEC--TAN"]
+        w.wcs.set_pv([(2,1,45.0)])
+
+        world=np.stack((ra,de),axis=-1)
+        pix=w.wcs_world2pix(world,1)
+        rx,ry=pix[:,0],pix[:,1]
+
+        ax,cov_q,infodict,mesg,ierr=optimize.leastsq(residual,[0.0,0.0,0.0],args=(dx,dy,rx),full_output=1)
+        ay,cov_q,infodict,mesg,ierr=optimize.leastsq(residual,[0.0,0.0,0.0],args=(dx,dy,ry),full_output=1)
+
+        ra0,de0=w.wcs_pix2world([[ax[0],ay[0]]],1)[0]
+
+    w=wcs.WCS(naxis=2)
+    w.wcs.crpix=np.array([x0,y0])
+    w.wcs.crval=np.array([ra0,de0])
+    w.wcs.cd=np.array([[ax[1],ax[2]],[ay[1],ay[2]]])
+    w.wcs.ctype=["RA---TAN","DEC--TAN"]
+    w.wcs.set_pv([(2,1,45.0)])
+
+    whdr={"CRPIX1":x0,"CRPIX2":y0,"CRVAL1":ra0,"CRVAL2":de0,"CD1_1":ax[1],"CD1_2":ax[2],"CD2_1":ay[1],"CD2_2":ay[2],"CTYPE1":"RA---TAN","CTYPE2":"DEC--TAN","CUNIT1":"DEG","CUNIT2":"DEG"}
+
+    return whdr
+
+# FITS file
+fname="2018-02-12T17:39:55.270.fits"
+
+# Read fits
+hdu=fits.open(fname)
+
+# Fix to force WCS with two axes only
+hdu[0].header["NAXIS"]=2
+w=wcs.WCS(hdu[0].header)
+
+# Read data
+data=hdu[0].data
+zavg,zstd,zmax,znum=data
+ny,nx=zavg.shape
+
+# Get extrema
+vmin=np.mean(zavg)-2.0*np.std(zavg)
+vmax=np.mean(zavg)+3.0*np.std(zavg)
+
+# Read pixel catalog
+xs,ys,ms=read_pixel_catalog(fname+".cat")
+
+# Read Tycho2 catalog
+ra,de,mag=read_tycho2_catalog(9.0)
+
+# Convert to pixel coordinates
+xc,yc=w.all_world2pix(ra,de,1)
+c=(xc>0) & (xc<nx) & (yc>0) & (yc<ny)
+xt,yt=xc[c],yc[c]
+rat,det=ra[c],de[c]
+            
+# Match catalogs
+res=match_catalogs(rat,det,xt,yt,xs,ys,10.0)
+
+whdr=fit_wcs(res,nx//2,ny//2)
+
+    
+
+
+#hdu=fits.PrimaryHDU(header=w.to_header(),data=data)
+#hdr=fits.Header()
+hdr=hdu[0].header
+for k,v in whdr.items():
+    hdr[k]=v
+hdu=fits.PrimaryHDU(header=hdr,data=data)
+hdu.writeto("test.fits",overwrite=True)
+
+# Convert catalog stars
+plt.figure(figsize=(10,8))
+plt.imshow(zavg,origin="lower",aspect=1.0,interpolation="None",vmin=vmin,vmax=vmax)
+plt.scatter(xs,ys,s=80,edgecolors="r",facecolors="none")
+plt.scatter(xt,yt,s=150,edgecolors="y",facecolors="none")
+#plt.scatter(xt[c],yt[c],s=50,edgecolors="k",facecolors="none")
+#plt.xlim(0,720)
+#plt.ylim(0,576)
+plt.show()
diff --git a/python/extract_tracks.py b/python/extract_tracks.py
new file mode 100644
index 0000000..239987e
--- /dev/null
+++ b/python/extract_tracks.py
@@ -0,0 +1,130 @@
+#!/usr/bin/env python
+
+import sys
+from stio import fourframe,satid,observation
+import numpy as np
+import ppgplot
+import matplotlib.pyplot as plt
+from astropy import wcs
+from astropy.coordinates import SkyCoord
+
+def pos2rel(ff,x,y):
+    # Setup wcs
+    w=wcs.WCS(naxis=2)
+    w.wcs.crpix=ff.crpix
+    w.wcs.crval=ff.crval
+    w.wcs.cd=ff.cd
+    w.wcs.ctype=ff.ctype
+    w.wcs.set_pv([(2,1,45.0)])
+
+    world=w.wcs_pix2world(np.array([[x,y]]),1)
+
+    return world[0,0],world[0,1]
+
+def dec2sex(angle):
+    if angle<0.0:
+        sign="-"
+    else:
+        sign="+"
+    angle=60.0*np.abs(angle)
+    
+
+# Settings
+
+# Minimum predicted velocity (pixels/s)
+drdtmin=10.0
+
+# Track selection region around prediction (pixels)
+trkrmin=10.0
+
+# Track selection sigma
+trksig=5.0
+
+# Minimum track points
+ntrkmin=10
+
+fname=sys.argv[1]
+
+# Read four frame
+ff=fourframe(fname)
+
+# Read satelite IDs
+try:
+    f=open(fname+".id")
+except OSError:
+    print("Cannot open",fname+".id")
+else:
+    lines=f.readlines()
+    f.close()
+
+# ppgplot arrays
+tr=np.array([-0.5,1.0,0.0,-0.5,0.0,1.0])
+heat_l=np.array([0.0,0.2,0.4,0.6,1.0])
+heat_r=np.array([0.0,0.5,1.0,1.0,1.0])
+heat_g=np.array([0.0,0.0,0.5,1.0,1.0])
+heat_b=np.array([0.0,0.0,0.0,0.3,1.0])
+
+# Loop over identifications
+for line in lines:
+    # Decode
+    id=satid(line)
+
+    # Skip slow moving objects
+    drdt=np.sqrt(id.dxdt**2+id.dydt**2)
+    if drdt<drdtmin:
+        continue
+
+    # Extract significant pixels
+    x,y,t,sig=ff.significant(trksig,id.x0,id.y0,id.dxdt,id.dydt,trkrmin)
+
+    # Fit tracks
+    if len(x)>ntrkmin:
+        obs=observation(ff,x,y,t,sig)
+
+        obs.ra,obs.de=pos2rel(ff,obs.x0,obs.y0)
+        
+        # Plot
+        ppgplot.pgopen("/xs")
+        ppgplot.pgpap(0.0,1.0)
+        ppgplot.pgsvp(0.1,0.95,0.1,0.8)
+
+        ppgplot.pgsch(0.8)
+        ppgplot.pgmtxt("T",6.0,0.0,0.0,"UT Date: %.23s  COSPAR ID: %04d"%(ff.nfd,ff.site_id))
+        ppgplot.pgmtxt("T",4.8,0.0,0.0,"R.A.: %10.5f (%4.1f'') Decl.: %10.5f (%4.1f'')"%(ff.crval[0],3600.0*ff.crres[0],ff.crval[1],3600.0*ff.crres[1]))
+        ppgplot.pgmtxt("T",3.6,0.0,0.0,"FoV: %.2f\\(2218)x%.2f\\(2218) Scale: %.2f''x%.2f'' pix\\u-1\\d"%(ff.wx,ff.wy,3600.0*ff.sx,3600.0*ff.sy))
+        ppgplot.pgmtxt("T",2.4,0.0,0.0,"Stat: %5.1f+-%.1f (%.1f-%.1f)"%(np.mean(ff.zmax),np.std(ff.zmax),ff.vmin,ff.vmax))
+        
+        ppgplot.pgsch(1.0)
+        ppgplot.pgwnad(0.0,ff.nx,0.0,ff.ny)
+        ppgplot.pglab("x (pix)","y (pix)"," ")
+        ppgplot.pgctab(heat_l,heat_r,heat_g,heat_b,5,1.0,0.5)
+        
+        ppgplot.pgimag(ff.zmax,ff.nx,ff.ny,0,ff.nx-1,0,ff.ny-1,ff.vmax,ff.vmin,tr)
+        ppgplot.pgbox("BCTSNI",0.,0,"BCTSNI",0.,0)
+
+        ppgplot.pgsci(3)
+#        ppgplot.pgpt(x,y,4)
+
+        ppgplot.pgsci(4)
+        ppgplot.pgpt(np.array([obs.x0]),np.array([obs.y0]),4)
+        ppgplot.pgmove(obs.xmin,obs.ymin)
+        ppgplot.pgdraw(obs.xmax,obs.ymax)
+        
+        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()
+
+        # Track and stack
+    else:
+        ztrk=ff.track(id.dxdt,id.dydt,0.0)
+    
+#        plt.figure(figsize=(15,10))
+#        plt.imshow(ztrk,origin='lower')
+#        plt.scatter(id.x0,id.y0,s=150,edgecolors="y",facecolors="none")
+#        plt.show()
+
diff --git a/python/hough.py b/python/hough.py
new file mode 100644
index 0000000..fcdc0f3
--- /dev/null
+++ b/python/hough.py
@@ -0,0 +1,49 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from astropy.io import fits
+import cv2
+import sys
+
+hdu=fits.open(sys.argv[1])
+
+data=hdu[0].data
+zavg,zstd,zmax,znum=data
+
+zsig=(zmax-zavg)/(zstd+1e-9)
+
+ny,nx=zsig.shape
+tmp=np.zeros(nx*ny*3,dtype='uint8').reshape(ny,nx,3)
+tmp[:,:,0]=zmax.astype('uint8')
+tmp[:,:,1]=0.1*zmax.astype('uint8')
+tmp[:,:,2]=0.1*zmax.astype('uint8')
+img=cv2.cvtColor(tmp,cv2.COLOR_BGR2RGB)
+
+mask=np.array(zsig>5.0,dtype='uint8')
+tmp=np.zeros(nx*ny*3,dtype='uint8').reshape(ny,nx,3)
+tmp[:,:,0]=255*mask
+tmp[:,:,1]=255*mask
+tmp[:,:,2]=255*mask
+
+mat=cv2.cvtColor(tmp,cv2.COLOR_BGR2GRAY)
+lines=cv2.HoughLines(mat,1,np.pi/180.0,200)
+if lines is not None:
+    for rho,theta in lines[0]:
+        a = np.cos(theta)
+        b = np.sin(theta)
+        x0 = a*rho
+        y0 = b*rho
+        x1 = int(x0 + 1000*(-b))
+        y1 = int(y0 + 1000*(a))
+        x2 = int(x0 - 1000*(-b))
+        y2 = int(y0 - 1000*(a))
+
+        cv2.line(img,(x1,y1+10),(x2,y2+10),(0,0,255),2)
+
+minLineLength = 100
+maxLineGap = 10
+lines = cv2.HoughLinesP(mat,1,np.pi/180,100,minLineLength,maxLineGap)
+if lines is not None:
+    for x1,y1,x2,y2 in lines[0]:
+        cv2.line(img,(x1,y1),(x2,y2),(0,255,0),2)
+
+    cv2.imwrite(sys.argv[1]+'.jpg',img)
diff --git a/python/iod.py b/python/iod.py
new file mode 100644
index 0000000..609c8bc
--- /dev/null
+++ b/python/iod.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+
+class iod:
+    """IOD Observation"""
+
+    def __init__(self,line):
+        s=line.split()
+        self.iodline=line
+        self.norad=int(s[0])
+        self.cospar="%s %s"%(s[1],s[2])
+        self.site=int(s[3])
+        self.conditions=s[4]
+        self.datestring=s[5]
+
+
+        
+    def __repr__(self):
+        return "%05d %s %04d %s %s"%(self.norad,self.cospar,self.site,self.conditions,self.datestring)
+        
+
+iodlines=["10529 77 112D   4553 G 20120511214446298 17 25 2034286+472232 37 R",
+          "10529 77 112D   4553 G 20120511214452938 17 25 2038653+454274 37 S",
+          "23862 96 029D   4553 G 20120511221706217 17 25 2121225+480700 37 S",
+          "23862 96 029D   4553 G 20120511221726241 17 25 2122590+453398 37 R"]
+
+obslist=[iod(line) for line in iodlines]
+    
+print(type(obslist),type(obslist[0]))
+for obs in obslist:
+    print(obs)
diff --git a/python/plot.py b/python/plot.py
new file mode 100644
index 0000000..d1118d7
--- /dev/null
+++ b/python/plot.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+import numpy as np
+import sys
+from astropy.io import fits
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+
+fname="2018-02-12T17:39:55.270.fits"
+
+# Read fits
+hdu=fits.open(fname)
+
+# Read data
+data=hdu[0].data
+zavg,zstd,zmax,znum=data
+
+# Read header
+ny,nx=zavg.shape
+nz=hdu[0].header['NFRAMES']
+dt=[hdu[0].header['DT%04d'%i] for i in range(nz)]
+
+# Sigma frame
+zsig=(zmax-zavg)/(zstd+1e-9)
+
+# Position
+xm,ym=np.meshgrid(np.arange(nx),np.arange(ny))
+
+# Selection width
+w=10
+sigma=5.0
+
+# Selection
+c=(zsig>sigma) & (xm>w) & (xm<nx-w) & (ym>w) & (ym<ny-w)
+
+x=np.ravel(xm[c])
+y=np.ravel(ym[c])
+inum=np.ravel(znum[c]).astype('int')
+sig=np.ravel(zsig[c])
+t=np.array([dt[i] for i in inum])
+
+# Load predictions
+d=np.loadtxt(fname+".id",usecols=(1,2,3,4,5,6))
+
+x0=d[:,0]
+y0=d[:,1]
+t0=np.zeros_like(x0)
+x1=d[:,2]
+y1=d[:,3]
+t1=d[:,4]
+
+for i in range(len(x0)):
+    xr=x0[i]+(x1[i]-x0[i])*(t-t0[i])/(t1[i]-t0[i])
+    yr=y0[i]+(y1[i]-y0[i])*(t-t0[i])/(t1[i]-t0[i])
+    r=np.sqrt((x-xr)**2+(y-yr)**2)
+    c=r<10.0
+    print(i,np.sum(c))
+    plt.figure()
+    plt.plot([x0[i],x1[i]],[y0[i],y1[i]],'r')
+    plt.scatter(x[c],y[c],c=t[c],s=10.0*(sig[c]-sigma))
+    plt.scatter(x,y,c=t,s=1)
+    plt.show()
+
+#fig=plt.figure()
+#ax=fig.add_subplot(111,projection='3d')
+#ax.scatter(x,y,t,c=t,s=10*(sig-5))
+
+#for i in range(len(x0)):
+#    ax.plot([x0[i],x1[i]],[y0[i],y1[i]],[t0[i],t1[i]],'k')
+#plt.show()
diff --git a/python/stio.py b/python/stio.py
new file mode 100644
index 0000000..e0a9056
--- /dev/null
+++ b/python/stio.py
@@ -0,0 +1,173 @@
+#!/usr/bin/env python
+
+import numpy as np
+from astropy.io import fits
+from astropy.time import Time
+
+class observation:
+    """Satellite observation"""
+
+    def __init__(self,ff,x,y,t,sig):
+        """Fit a satellite track"""
+
+        # Get times
+        self.tmin=np.min(t)
+        self.tmid=np.mean(t)
+        self.tmax=np.max(t)
+        self.mjd=ff.mjd+self.tmid/86400.0
+        self.nfd=Time(self.mjd,format='mjd',scale='utc').isot
+        
+        # Very simple polynomial fit; no weighting, no cleaning
+        px=np.polyfit(t-self.tmid,x,1)
+        py=np.polyfit(t-self.tmid,y,1)
+
+        # Store
+        self.x0=px[1]
+        self.y0=py[1]
+        self.dxdt=px[0]
+        self.dydt=py[0]
+        self.xmin=self.x0+self.dxdt*(self.tmin-self.tmid)
+        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)
+        
+class satid:
+    """Satellite identifications"""
+
+    def __init__(self,line):
+        s=line.split()
+        self.nfd=s[0]
+        self.x0=float(s[1])
+        self.y0=float(s[2])
+        self.t0=0.0
+        self.x1=float(s[3])
+        self.y1=float(s[4])
+        self.t1=float(s[5])
+        self.norad=int(s[6])
+        self.catalog=s[7]
+        self.state=s[8]
+        self.dxdt=(self.x1-self.x0)/(self.t1-self.t0)
+        self.dydt=(self.y1-self.y0)/(self.t1-self.t0)
+
+    def __repr__(self):
+        return "%s %f %f %f -> %f %f %f %d %s %s"%(self.nfd,self.x0,self.y0,self.t0,self.x1,self.y1,self.t1,self.norad,self.catalog,self.state)
+        
+class fourframe:
+    """Four frame class"""
+
+    def __init__(self,fname=None):
+        if fname==None:
+            # Initialize empty fourframe
+            self.nx=0
+            self.ny=0
+            self.nz=0
+            self.mjd=-1
+            self.nfd=None
+            self.zavg=None
+            self.zstd=None
+            self.zmax=None
+            self.znum=None
+            self.dt=None
+            self.site_id=0
+            self.observer=None
+            self.texp=0.0
+            self.fname=None
+            self.crpix=np.array([0.0,0.0])
+            self.crval=np.array([0.0,0.0])
+            self.cd=np.array([[1.0,0.0],[0.0,1.0]])
+            self.ctype=["RA---TAN","DEC--TAN"]
+            self.cunit=np.array(["deg","deg"])
+            self.crres=np.array([0.0,0.0])
+        else:
+            # Read FITS file
+            hdu=fits.open(fname)
+
+            # Read image planes
+            self.zavg,self.zstd,self.zmax,self.znum=hdu[0].data
+
+            # Frame properties
+            self.ny,self.nx=self.zavg.shape
+            self.nz=hdu[0].header['NFRAMES']
+    
+            # Read frame time oselfsets
+            self.dt=np.array([hdu[0].header['DT%04d'%i] for i in range(self.nz)])
+    
+            # Read header
+            self.mjd=hdu[0].header['MJD-OBS']
+            self.nfd=hdu[0].header['DATE-OBS']
+            self.site_id=hdu[0].header['COSPAR']
+            self.observer=hdu[0].header['OBSERVER']
+            self.texp=hdu[0].header['EXPTIME']
+            self.fname=fname
+
+            # Astrometry keywords
+            self.crpix=np.array([hdu[0].header['CRPIX1'],hdu[0].header['CRPIX2']])
+            self.crval=np.array([hdu[0].header['CRVAL1'],hdu[0].header['CRVAL2']])
+            self.cd=np.array([[hdu[0].header['CD1_1'],hdu[0].header['CD1_2']],
+                              [hdu[0].header['CD2_1'],hdu[0].header['CD2_2']]])
+            self.ctype=[hdu[0].header['CTYPE1'],hdu[0].header['CTYPE2']]
+            self.cunit=[hdu[0].header['CUNIT1'],hdu[0].header['CUNIT2']]
+            self.crres=np.array([hdu[0].header['CRRES1'],hdu[0].header['CRRES2']])
+
+        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
+        self.wy=self.ny*self.sy
+        self.vmin=np.mean(self.zmax)-2.0*np.std(self.zmax)
+        self.vmax=np.mean(self.zmax)+6.0*np.std(self.zmax)
+
+
+    def significant(self,sigma,x0,y0,dxdt,dydt,rmin=10.0):
+        """Extract significant points"""
+
+        # Generate sigma frame
+        zsig=(self.zmax-self.zavg)/(self.zstd+1e-9)
+
+        # Select
+        c=(zsig>sigma)
+
+        # Positions
+        xm,ym=np.meshgrid(np.arange(self.nx),np.arange(self.ny))
+        x,y=np.ravel(xm[c]),np.ravel(ym[c])
+        inum=np.ravel(self.znum[c]).astype('int')
+        sig=np.ravel(zsig[c])
+        t=np.array([self.dt[i] for i in inum])
+
+        # Predicted positions
+        xr=x0+dxdt*t
+        yr=y0+dydt*t
+        r=np.sqrt((x-xr)**2+(y-yr)**2)
+        c=r<rmin
+
+        return x[c],y[c],t[c],sig[c]
+
+    def track(self,dxdt,dydt,tref):
+        """Track and stack"""
+        # Empty frame
+        ztrk=np.zeros_like(self.zavg)
+        
+        # Loop over frames
+        for i in range(self.nz):
+            dx=int(np.round(dxdt*(self.dt[i]-tref)))
+            dy=int(np.round(dydt*(self.dt[i]-tref)))
+            if dx>=0:
+                i1min,i1max=dx,self.nx-1
+                i2min,i2max=0,self.nx-dx-1
+            else:
+                i1min,i1max=0,self.nx+dx-1
+                i2min,i2max=-dx,self.nx-1
+            if dy>=0:
+                j1min,j1max=dy,self.ny-1
+                j2min,j2max=0,self.ny-dy-1
+            else:
+                j1min,j1max=0,self.ny+dy-1
+                j2min,j2max=-dy,self.ny-1
+            zsel=np.where(self.znum==i,self.zmax,0.0)
+            ztrk[j2min:j2max,i2min:i2max]+=zsel[j1min:j1max,i1min:i1max]
+
+        return ztrk
+                
+    def __repr__(self):
+       return "%s %dx%dx%d %s %.3f %d %s"%(self.fname,self.nx,self.ny,self.nz,self.nfd,self.texp,self.site_id,self.observer)
+
+   
diff --git a/python/sun.py b/python/sun.py
new file mode 100644
index 0000000..4468550
--- /dev/null
+++ b/python/sun.py
@@ -0,0 +1,34 @@
+#!/usr/bin/env python
+from astropy.coordinates import SkyCoord,EarthLocation,AltAz,get_sun
+from astropy.time import Time
+import astropy.units as u
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import optimize,interpolate
+
+
+
+loc=EarthLocation(lat=52.8344*u.deg,lon=6.3785*u.deg,height=10*u.m)
+t0=Time.now()
+refalt=-6.0
+
+# Compute solar position every 10 minutes
+t=Time(t0.mjd+np.arange(144)/144.0,format='mjd',scale='utc')
+psun=get_sun(t)
+hor=psun.transform_to(AltAz(obstime=t,location=loc))
+
+# Interpolating function
+alt_diff=hor.alt.deg-refalt
+f=interpolate.interp1d(t.mjd,alt_diff)
+
+# Find sunrise/sunset
+sign=alt_diff[1:]*alt_diff[:-1]
+idx=np.where(sign<0.0)[0]
+for i in idx:
+    # Set
+    if alt_diff[i]>alt_diff[i+1]:
+        tset=Time(optimize.bisect(f,t[i].mjd,t[i+1].mjd),format='mjd',scale='utc')
+    # Rise
+    else:
+        trise=Time(optimize.bisect(f,t[i].mjd,t[i+1].mjd),format='mjd',scale='utc')
+
diff --git a/python/vis.py b/python/vis.py
new file mode 100644
index 0000000..5a05209
--- /dev/null
+++ b/python/vis.py
@@ -0,0 +1,11 @@
+#!/usr/bin/env python
+from skyfield.positionlib import ICRF
+from skyfield.api import load,Topos
+
+ts=load.timescale()
+t=ts.now()
+observer=Topos(latitude_degrees=52.8344,longitude_degrees=6.3785,elevation_m=10.0)
+p=ICRF([0.0,0.0,0.0],observer_data=observer,t=t)
+q=p.from_altaz(alt_degrees=30.0,az_degrees=30.0)
+
+print(p,q)