offset=0.0
CIC=1

e=fltarr(9,25)
s=fltarr(9,25)
;set plotting for postscript
set_plot,'PS'

loadct,39
device,/color
device,ysize=8.5,/inches
device,xsize=7.0,/inches
device,yoffset=1.0,/inches
device,filename='/afs/slac/u/ki/lances/idl/psf.ps'

;for q=0,8 do begin

;declare tt a variable =0
tt=0


filename='./lsst_sim/output_0'
rdfloat,filename,x,y,z,vx,vy,vz,xo,yo,zo,vxo,yo,vzo,n

r=sqrt(x*x+y*y)*abs(y)/y
ro=sqrt(xo*xo+yo*yo)*abs(yo)/yo

;declare variable called CIC
CIC=1

;make array of 15 elements 
cc=fltarr(15)
cc(1)=60
cc(2)=240
cc(3)=160
cc(4)=240
cc(5)=60
cc(6)=240
cc(7)=60
cc(8)=240
cc(9)=60
cc(10)=240
cc(11)=60
cc(12)=240
cc(13)=200
cc(14)=200

!p.multi=[0,5,5]
!x.margin=[0,0]
!y.margin=[0,0]

aa=findgen(11)/10*2*!Pi
usersym,0.1*cos(aa),0.1*sin(aa),/fill

for kk=0,24 do begin

if kk mod 5 eq 0 then xoff=-2.0
if kk mod 5 eq 1 then xoff=-1.0
if kk mod 5 eq 2 then xoff=0.0
if kk mod 5 eq 3 then xoff=1.0
if kk mod 5 eq 4 then xoff=2.0
if kk/5 eq 0 then yoff=-2.0
if kk/5 eq 1 then yoff=-1.0
if kk/5 eq 2 then yoff=0.0
if kk/5 eq 3 then yoff=1.0
if kk/5 eq 4 then yoff=2.0


;xoff=0.0+xoff*0.6
;yoff=yoff*0.6
xoff=0.0
yoff=0.0

ff=182.575

size=8.0
plot,[-size,size],[-size*8.5/7.0,size*8.5/7.0],/NODATA,/xstyle,/ystyle

xq=fix((size*8.5/7.0))
if (xq mod 2 eq 0) then xq=xq+1
for qq=-xq,xq,2 do begin
    xq=findgen(1000)/1000.*size*8.5/7.0*2-size*8.5/7.0
    oplot,xq,fltarr(N_elements(xq))+qq*5.0,linestyle=1
    oplot,fltarr(N_elements(xq))+qq*5.0,xq,linestyle=1

endfor

spotx=fltarr(N_elements(x))
spoty=fltarr(N_elements(x))
spotc=fltarr(N_elements(x))

j=0L
for i=0L,N_elements(x)-1 do begin

    if (n(i) eq 1) then begin
       currz=sqrt((x(i)*x(i)+y(i)*y(i)))
    endif

    if n(i) eq 12 then begin
      spotx(j)=xo(i)
      spoty(j)=yo(i)
      spotc(j)=currz
     j=j+1
    endif
    



endfor

spotx=spotx(0:j-1)
spoty=spoty(0:j-1)
spotc=spotc(0:j-1)

xf=median(spotx(where(abs(spotx-xoff*ff) lt 2 and abs(spoty-yoff*ff) lt 2)))
yf=median(spoty(where(abs(spotx-xoff*ff) lt 2 and abs(spoty-yoff*ff) lt 2)))

print,xf,yf

;xf=(max(spotx)-min(spotx))/2.
;yf=(max(spoty)-min(spoty))/2.

;xf=0
;yf=0


gg=where(abs(spotx-xf) lt size/1000. and abs(spoty-yf) lt size/1000.)
xv=(spotx(gg)-xf)*1000.
yv=(spoty(gg)-yf)*1000.

;if (xoff eq 0.0 and yoff eq 0.0 and q eq 1) then goto, jumpe

resultx=moment(xv)
resulty=moment(yv)
covxy=correlate(xv,yv,/covariance)

ellip=((resultx(1)-resulty(1))^2+(2.0*covxy)^2)/(resultx(1)+resulty(1))^2
ellip=sqrt(ellip)
rms=sqrt(resultx(1)+resulty(1))
pa=0.5*atan((2.0*covxy),(resultx(1)-resulty(1)))

print,rms,ellip,pa

;e(q,kk)=ellip
;s(q,kk)=rms

endfor

END