offset=0.0
CIC=1

;Create two arrays 9*25
e=fltarr(9,25)
s=fltarr(9,25)

;set plotting for postscript file
set_plot,'PS'

;generate multiple plots on same page 2*2
!p.multi=[0,2,2]
;Reset plotting to 1 plot per page
!p.multi=0

;load a particular color table
loadct,39

;set particular qualities of plot
device,/color
device,ysize=8.5,/inches
device,xsize=7.0,/inches
device,yoffset=1.0,/INCHES


;for q=0,8 do begin

;declare tt a variable =0
tt=0

;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



;if q lt 10 then filename='output_'+string(q,format='(I1)')
;if q ge 10 then filename='output_'+string(q,format='(I2)')

;read in data from output_0 and output_1 so they can be plotted
filename='./lsst_sim/output_0'
rdfloat,filename,x,y,z,vx,vy,vz,xo,yo,zo,vxo,vyo,vzo,n

;if q ne 0 then goto,jump2
jump1:

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

if tt eq 0 then begin
plot,[-11000*7.0/8.5/2.,11000*7.0/8.5/2.],[-2000,9000],/NODATA,/xstyle,/ystyle
endif else begin
plot,[-2000*7.0/8.5/2.,2000*7.0/8.5/2.],[3400,5400],/NODATA,/xstyle,/ystyle
endelse

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

    uu=randomu(seed,1)

    if (uu(0) lt (1000./float(N_elements(x))) ) then begin

    oplot,[x(i),xo(i)],[z(i),zo(i)],color=cc(n(i))*CIC
   
endif


endfor


aaa=[4180,1625,2533,800,800,550,530,385,385,365,365,320]
bbb=[2613,800, 0 ,   0,   0, 0,   0,  0,  0,  0,  0,  0]
ccc=[19200,6032,8577,2739,3803,5198,2058,5630,5630,3625,17192,0]
ddd=[0,6.034,6.034-1.822-4.436,6.035-1.822-4.436+4.038,6.035-1.822-4.436+4.038+0.068,6.035-1.822-4.436+4.038+0.068+0.508,6.035-1.822-4.436+4.038+0.068+0.508+0.03,6.035-1.822-4.436+4.038+0.068+0.508+0.03+0.368,6.035-1.822-4.436+4.038+0.068+0.508+0.03+0.368+0.016,6.035-1.822-4.436+4.038+0.068+0.508+0.03+0.368+0.016+0.045,6.035-1.822-4.436+4.038+0.068+0.508+0.03+0.368+0.016+0.045+0.06,6.035-1.822-4.436+4.038+0.068+0.508+0.03+0.368+0.016+0.045+0.06+0.025]*1e3

for iii=0,N_elements(aaa)-1 do begin

if ccc(iii) ne 0 then begin
xq=findgen(1000)/1000*(aaa(iii)-bbb(iii))+bbb(iii)
yq=-xq
zq=ddd(iii)+xq*xq/2/(ccc(iii))
oplot,xq,zq,thick=3
oplot,yq,zq,thick=3
xq=findgen(1000)/1000*aaa(iii)*2-aaa(iii)
oplot,xq,fltarr(N_elements(xq))+zq(N_elements(zq)-1),linestyle=2
xq=findgen(1000)/1000*bbb(iii)*2-bbb(iii)
oplot,xq,fltarr(N_elements(xq))+zq(0),linestyle=2
endif else begin

xq=findgen(1000)/1000*(aaa(iii)-bbb(iii))+bbb(iii)
yq=-xq
zq=fltarr(N_elements(xq))+ddd(iii)
oplot,xq,zq,thick=3
oplot,yq,zq,thick=3
endelse

endfor










if tt eq 1 then goto,jump2

tt=1
goto,jump1


jump2:
;device,/close
;END
;if q eq 0 then tt='In focus'
;if q eq 1 then tt='Primary +10 microns'
;if q eq 2 then tt='Secondary +10 microns'
;if q eq 3 then tt='Tertiary +10 microns'
;if q eq 4 then tt='L1 +10 microns'
;if q eq 5 then tt='L2 +10 microns'
;if q eq 6 then tt='Filter +10 microns'
;if q eq 7 then tt='L3 +10 microns'
;if q eq 8 then tt='Focal Plane +10 microns'


!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=xoff*0.6
yoff=yoff*0.6
;xoff=0.0
;yoff=0.0

ff=1825.75
;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



for qqq=0L,N_elements(spotx)-1 do begin
oplot,[(spotx(qqq)-xf)*1000],[(spoty(qqq)-yf)*1000.],color=250.*(spotc(qqq)-2613)/(4180-2613)*CIC,psym=8
endfor

;nx=fix(2*size/10)
;ny=fix(2*size*8.5/7.0/10)
;image=fltarr(nx,ny)
;for qqq=0L,N_elements(spotx)-1 do begin
;    xv=(spotx(qqq)-xf)*1000/10.
;    yv=(spoty(qqq)-yf)*1000/10.
;    if (xv ge 0 and xv lt nx and $
;        yv ge 0 and yv lt ny) then image(xv,yv)=image(xv,yv)+1
;endfor

;!p.multi=0
;plot,[0,1],[0,1],/nodata
;tvscl,alog10(image>1)

xyouts,-0.6*size,0.9*size,'!4r!3 = '+string(rms,format='(F6.2)')
xyouts,-0.6*size,0.7*size, 'e = '+string(ellip,format='(F6.3)')
xyouts,-0.6*size,0.5*size,'!4u!3 = '+string(pa,format='(F6.2)')



endfor




;endfor

device,/close
set_plot,'X'
!x.margin=[10,3]
!y.margin=[4,2]
END

!p.multi=[0,1,3]
!x.margin=[10,3]
!y.margin=[4,2]
plot,[0,10],[0.,30.0],title='PSF Size',xtitle='Optic Number',ytitle='RMS (microns)',/nodata
xx=findgen(9)+1.
for kk=0,24 do oplot,xx,s(*,kk),psym=4,color=(float(kk)/24*240.+10)
for kk=0,8 do oplot,[kk+1],[mean(s(kk,*))],psym=5,symsize=2.0

plot,[0,10],[1e-2,1.0],/ylog,title='PSF Ellipticity',ytitle='Root Quadropole Moment',xtitle='Optic Number',/nodata
for kk=0,24 do oplot,xx,e(*,kk),psym=4,color=(float(kk)/24*240.+10)
for kk=0,8 do oplot,[kk+1],[mean(e(kk,*))],psym=5,symsize=2.0


plot,[0,10],[0.,6.0],title='PSF Ellipticity*Size',xtitle='Optic Number',ytitle='Sigma*Ellipticity (microns)',/nodata
xx=findgen(9)+1.
for kk=0,24 do oplot,xx,e(*,kk)*s(*,kk),psym=4,color=(float(kk)/24*240.+10)
for kk=0,8 do oplot,[kk+1],[mean(e(kk,*)*s(kk,*))],psym=5,symsize=2.0

!p.multi=0

jumpe:
device,/close
set_plot,'X'
!x.margin=[10,3]
!y.margin=[4,2]
END