voffset=0.0
CIC=1
e=fltarr(9,25)
s=fltarr(9,25)

;elf=fltarr(10)
;e1ff=fltarr(10)
;e2ff=fltarr(10)


set_plot,'PS'
!p.multi=[0,2,2]
!p.multi=0
loadct,39
device,/color
device,ysize=8.5,/inches
device,xsize=7.0,/inches
device,yoffset=1.0,/INCHES
device,filename='idl_x.ps'

ntrials=81
e1f=fltarr(ntrials)
e2f=fltarr(ntrials)
pxf=fltarr(ntrials)
pyf=fltarr(ntrials)
paf=fltarr(ntrials)
rmsf=fltarr(ntrials)

;for trial=0,ntrials-1 do begin
for trial=0,80 do begin

ngridpoints=1
;name='_0_30000points'
if trial lt 10 then name='_'+string(trial,format='(I1)')
if trial ge 10 and trial lt 100 then name='_'+string(trial,format='(I2)')
if trial ge 100 and trial lt 1000 then name='_'+string(trial,format='(I3)')

;centerx=0.2*(-4+(trial mod 9))
;centery=1.5*0.0+0.2*(4-(trial / 9))
;print,(-4+(trial mod 9)),(4-(trial / 9))

centerx=0.0+0.30*(-4+(trial mod 9))
centery=0.0+0.30*(4-(trial / 9))
print,(-4+(trial mod 9)),(4-(trial / 9))
print,centerx,centery

;gridstep=0.6
;centerx=1.5*0.0+gridstep*(-2+(trial mod 5))
;centery=gridstep*(2-(trial / 5))
;print,(-2+(trial mod 5)),(2-(trial / 5))



tt=0
CIC=1
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




filename='output'+name+'.fits'
;rdfloat,filename,x,y,z,vx,vy,vz,xo,yo,zo,vxo,vyo,vzo,n,origx,origy
data=mrdfits(filename,1)
x=data.xi
y=data.yi
z=data.zi
vx=data.vxi
vy=data.vyi
vz=data.vzi
xo=data.xf
yo=data.yf
zo=data.zf
vxo=data.vxf
vyo=data.vyf
vzo=data.vzf
origx=data.vxo
origy=data.vyo
n=data.surface

;rdfloat,filename,x,y,z,vx,vy,vz,xo,yo,zo,vxo,vyo,vzo,n,origx,origy

jump1:

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


if trial eq 0 then begin

if tt eq 0 then begin
plot,[-11000*7.0/8.5/2.,11000*7.0/8.5/2.],[-3000,8000],/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+07.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

endif










if tt eq 1 then goto,jump2

tt=1
goto,jump1





jump2:

if (trial eq 0) then begin

npp=sqrt(ntrials)
!p.multi=[0,npp,npp]
!x.margin=[0,0]
!y.margin=[0,0]
endif

spotx=fltarr(N_elements(x))
spoty=fltarr(N_elements(x))
spotox=fltarr(N_elements(x))
spotoy=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)
      spotox(j)=origx(i)
      spotoy(j)=origy(i)
      spotc(j)=currz
     j=j+1
    endif
    
endfor
spotx=spotx(0:j-1)
spoty=spoty(0:j-1)
spotox=spotox(0:j-1)
spotoy=spotoy(0:j-1)
spotc=spotc(0:j-1)



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



for kk=0,ngridpoints*ngridpoints-1 do begin


xoff=-(ngridpoints/2.0-0.5)+fix(kk mod ngridpoints)
yoff=(ngridpoints/2.0-0.5)-fix(kk/ngridpoints)



xoff=centerx+xoff*1./360.
yoff=centery+yoff*1./360.



;ff=182.575
ff=!Pi/180.

size=13.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



xtol=min(abs(spotox-xoff*ff))
ytol=min(abs(spotoy-yoff*ff))

g1=where(abs(spotox-xoff*ff) le xtol and abs(spotoy-yoff*ff) le ytol)

if (N_elements(g1) gt 1 and xtol lt !Pi/180/3600. and ytol lt !Pi/180/3600.) then begin

print,N_elements(g1)
xf=median(spotx(g1))
yf=median(spoty(g1))
xv=(spotx(g1)-xf)*1000.
yv=(spoty(g1)-yf)*1000.



medx=median(xv)
medy=median(yv)
seeing=28.0

for ty=0,99 do begin
;weighted sum
;seeing in microns


weight=exp(-((xv-medx)^2+(yv-medy)^2)/2/seeing/seeing)/((sqrt(2*!Pi*seeing*seeing))^2)
covxy=total(weight*xv*yv)/total(weight)-total(weight*xv)*total(weight*yv)/total(weight)/total(weight)
resultx=total(weight*xv*xv)/total(weight)-total(weight*xv)*total(weight*xv)/total(weight)/total(weight)
resulty=total(weight*yv*yv)/total(weight)-total(weight*yv)*total(weight*yv)/total(weight)/total(weight)


medx=total(weight*xv)/total(weight)
medy=total(weight*yv)/total(weight)
seeing=sqrt(resultx+resulty)

endfor


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

print,rms,ellip,pa,e1,e2,resultx,resulty

;if (trial eq 0 ) then begin
for qqq=0L,N_elements(spotx)-1 do begin
    uu=randomu(seed,1)

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

oplot,[(spotx(qqq)-xf)*1000],[(spoty(qqq)-yf)*1000.],color=250.*(spotc(qqq)-2613)/(4180-2613)*CIC,psym=8

endif

endfor
;endif


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


e1f(trial)=e1
e2f(trial)=e2
paf(trial)=pa
rmsf(trial)=rms
pxf(trial)=xoff
pyf(trial)=yoff

endif


endfor


endfor






ii=0L
rr=fltarr(1000000)
;sc=fltarr(100000)
;scv=fltarr(100000)
e1diff=fltarr(1000000)
e2diff=fltarr(1000000)
ediff=fltarr(1000000)

;for kk=0,ngridpoints*ngridpoints-1 do begin
;for jj=0,ngridpoints*ngridpoints-1 do begin
for kk=0,ntrials-1 do begin
for jj=0,ntrials-1 do begin

;if kk ne jj then begin


e1diff(ii)=abs(e1f(jj)-e1f(kk))
e2diff(ii)=abs(e2f(jj)-e2f(kk))
ediff(ii)=sqrt((e1f(jj)-e1f(kk))^2+(e2f(jj)-e2f(kk))^2)
;e1a=e1f(kk)/4.0/(1-e1f(kk)*e1f(kk)-e2f(kk)*e2f(kk))
;e1b=e1f(jj)/4.0/(1-e1f(jj)*e1f(jj)-e2f(jj)*e2f(jj))
;e2a=e2f(kk)/4.0/(1-e1f(kk)*e1f(kk)-e2f(kk)*e2f(kk))
;e2b=e2f(jj)/4.0/(1-e1f(jj)*e1f(jj)-e2f(jj)*e2f(jj))

rr(ii)=sqrt((pxf(kk)-pxf(jj))^2+(pyf(kk)-pyf(jj))^2)

;sc(ii)=(e1a*e1b+e2a*e2b)/sqrt(e1a^2+e2a^2)/sqrt(e1b^2+e2b^2)
;scv(ii)=(e1a*e1b+e2a*e2b)
       
ii=ii+1

;endif

endfor
endfor
rr=rr(0:ii-1)
;sc=sc(0:ii-1)
;scv=scv(0:ii-1)
e1diff=e1diff(0:ii-1)
e2diff=e2diff(0:ii-1)
ediff=ediff(0:ii-1)


nbins=20

bins=max(rr(where(ediff ne 0)))/nbins
zz=fltarr(nbins)
zze=fltarr(nbins)
zzv=fltarr(nbins)
zzev=fltarr(nbins)
zzu=fltarr(nbins)
zzeu=fltarr(nbins)
xz=fltarr(nbins)
zzn=fltarr(nbins)

for q=0,nbins-1 do begin
y=0L
for kk=0L,N_elements(rr)-1 do begin
    if (rr(kk) ge float(q)*bins-bins/2. and rr(kk) lt float(q)*bins+bins/2.) then begin

        if (y eq 0) then tt=e1diff(kk)
        if (y gt 0) then tt=[tt,e1diff(kk)]
        if (y eq 0) then ttv=e2diff(kk)
        if (y gt 0) then ttv=[ttv,e2diff(kk)]
        if (y eq 0) then ttu=ediff(kk)
        if (y gt 0) then ttu=[ttu,ediff(kk)]
        y=y+1
    endif
endfor

if (y gt 1) then begin
result=moment(tt)     
zz(q)=result(0)
zze(q)=sqrt(result(1))

xz(q)=bins*(float(q))

result=moment(ttv)     
zzv(q)=result(0)
zzev(q)=sqrt(result(1))

result=moment(ttu)     
zzu(q)=result(0)
zzeu(q)=sqrt(result(1))


zzn(q)=1
endif else begin
zzn(q)=0
xz(q)=bins*(float(q))
endelse
   

endfor


!p.multi=0
!x.margin=[10,3]
!y.margin=[4,2]
!p.multi=[0,1,3]

g=where(zzn ne 0)

xz=xz(g)
zz=zz(g)
zze=zze(g)
zzev=zzev(g)
zzv=zzv(g)


g=where(e1f ne 0)
meany=mean(sqrt(e1f(g)^2+e2f(g)^2))


plot,xz,zz,psym=4,xtitle='Separation (degrees)',ytitle='Abs(E!D1!N Residual)',yr=[0,meany*2.0]
errplot,xz,zz-zze,zz+zze,width=0.0

xx=findgen(1000)/1000.*max(xz)*2.0
oplot,xx,fltarr(1000)+meany*0.81,linestyle=1

plot,xz,zzv,psym=4,xtitle='Separation (degrees)',ytitle='Abs(E!D2!N Residual)',yr=[0,meany*2.0]
errplot,xz,zzv-zzev,zzv+zzev,width=0.0

oplot,xx,fltarr(1000)+meany*0.81,linestyle=1

plot,xz,zzu,psym=4,xtitle='Separation (degrees)',ytitle='Ellipticity Residual',yr=[0,meany*2.0]
errplot,xz,zzu-zzeu,zzu+zzeu,width=0.0

oplot,xx,fltarr(1000)+meany*1.28,linestyle=1

;elf(trial)=ellip
;e1ff(trial)=e1
;e2ff(trial)=e2



;result=moment(elf)
;print,result(0),sqrt(result(1))
;result=moment(e1ff)
;print,result(0),sqrt(result(1))
;result=moment(e2ff)
;print,result(0),sqrt(result(1))

g=where(e1f ne 0)

!p.multi=[0,1,3]
q=histogram(paf(g),min=-!Pi/2.,bin=0.2)
xx=findgen(N_elements(q))*0.2-!Pi/2.
plot,xx,q,psym=10,title='Position Angle'

q=histogram(sqrt(e1f(g)^2+e2f(g)^2),min=0,bin=max(sqrt(e1f(g)^2+e2f(g)^2))/10.)
xx=findgen(N_elements(q))*max(sqrt(e1f(g)^2+e2f(g)^2))/10.
plot,xx,q,psym=10,title='Ellipticity'

q=histogram(rmsf(g),min=min(rmsf(g)),bin=(max(rmsf(g))-min(rmsf(g)))/10.)
xx=findgen(N_elements(q))*(max(rmsf(g))-min(rmsf(g)))/10.+min(rmsf(g))
plot,xx,q,psym=10,title='RMS Size'

!p.multi=0
sqnt=fix(sqrt(ntrials))
e1v=fltarr(sqnt,sqnt)
e2v=fltarr(sqnt,sqnt)
k=0
for j=sqnt-1,0,-1 do begin
for i=0,sqnt-1 do begin
;    e1v(i,j)=e1f(k)
;   e2v(i,j)=e2f(k)
    e1v(i,j)=sqrt((e1f(k))^2+(e2f(k))^2)*cos(paf(k)+!Pi/2.)
    e2v(i,j)=sqrt((e1f(k))^2+(e2f(k))^2)*sin(paf(k)+!Pi/2.)
    k=k+1
endfor
endfor
minx=min(pxf)
maxx=max(pxf)
miny=min(pyf)
maxy=max(pyf)
rx=maxx-minx
ry=maxy-miny
dx=rx/(sqnt-1)
dy=ry/(sqnt-1)
plot,[minx-dx,maxx+dx],[miny-dy,maxy+dy],/xstyle,/ystyle,/nodata,$
     xr=[minx-dx,maxx+dx],yr=[miny-dy,maxy+dy]

;plot,[-0.5,sqnt-0.5],[-0.5,sqnt-0.5],/nodata,xr=[-0.5,sqnt-0.5],$
;     yr=[-0.5,sqnt-0.5],/xstyle,/ystyle
;velovect,e1v,e2v
maxe=max(sqrt(e1v^2+e2v^2))
e1v=e1v/0.3
e2v=e2v/0.3
for j=sqnt-1,0,-1 do begin
for i=0,sqnt-1 do begin
    oplot,[i*dx+minx-e1v(i,j)/2.*dx,i*dx+minx+e1v(i,j)/2.*dx],$
          [j*dy+miny-e2v(i,j)/2.*dy,j*dy+miny+e2v(i,j)/2.*dy]
endfor
endfor

device,/close
set_plot,'X'

print,mean(rmsf(g)),mean(sqrt(e1f(g)^2+e2f(g)^2))

END