;   PROGRAM                  PSF.PRO

; Last Updated            17 Mar 2000          Mattia Vaccari

; This program produces different GAIA PSFs and represents their properties
;  in four diagrams, starting from the GAIA PSFs generated by syntfits.f
; The PSFs can be read from an IDL array or from a gzipped fits file.
; The diagrams can be displayed on screen or saved in a ps or eps file whose
;  size can be chosen to be 17*6 cm^2 (1 row layout) or 21*29.7 cm^2
;  (an a4 page, 2*2 layout)

; PARAMETERS
step=9.3
; step of SAG_LL_025's PSFs in the x direction
rf=3.
; rebinning factor
rstep=step/rf
; Step in both directions after rebinning

l1 :
print,'Where is the PSF contained, in an array (a) or in a gzipped '+$
 'fits file (f)?'
af='...' & read,af
if af ne 'a' and af ne 'f' then goto,l1

if af eq 'a' then begin
print,'Enter the name of the array:'
array='...' & read,array
result=execute('psf='+array)
psfdim=size(psf)/[1,3,3*rf,1,1]
endif

if af eq 'f' then begin
l2 :
print,'Enter the name of the file without the .fits.gz extension:'
ffn='...' & read,ffn
psf=readfits('psf/'+ffn+'.fits.gz',/silent)
psfdim=size(psf)
if psfdim(0) eq 0 then goto,l2
endif

l3 :
print,'Is it an optical PSF (o) or a global Astro-2 BBP PSF (g)?
og='...' & read,og
if og ne 'o' and og ne 'g' then goto,l3
if og eq 'o' then begin side=744. & intside=511.5 & endif
if og eq 'g' then begin side=1488. & intside=1023. & endif
; side is the side of the 2-D PSF and 1-D PSFs images
; intside is the side of the integrated PSF image
ran=[side/2,side/2]
intran=[0,intside]
; ran and intran are the ranges of the corresponding images
int=side/step
; int is the number of intervals of step width plotted in the 2-D PSF
;  and 1-D PSFs images
rint=intside/rstep
; rint is the number of intervals of rstep width plotted in the
;  integrated PSF image

l4:
print,'Do you want a single-scan PSF (s) or a multiple-scan PSF (m)?'
sc='s' & read,sc
if sc ne 's' and sc ne 'm' then goto,l4

if sc eq 'm' then begin
print,'How many scans do you want to stack?'
npsf=1 & read,npsf

l5 :
print,'How do you want the position angles to be distributed?'
print,'    1=0-360    2=0-45 or 90-135    3=0-90    4=0-45'
dir=1 & read,dir
if dir eq 1 then posa=randomu(seed,npsf)*360. $
else if dir eq 2 then begin
posa=randomu(seed,npsf)*45.
posa=posa+((-1)^round(randomu(seed,npsf)*1000000)*45.+45.)
endif $
else if dir eq 3 then posa=randomu(seed,npsf)*90. $
else if dir eq 4 then posa=randomu(seed,npsf)*45. $
else goto,l5
endif

l6 :
print,'Do you want to display it on screen (s) or to save it in a file (f)?
scps='s' & read,scps
if scps ne 's' and scps ne 'f' then goto,l6

l7 :
if scps eq 'f' then begin
print,'Do you want a small (17*6 cm^2) or big (a4) file (s/b)?'
fisi='b' & read,fisi
if fisi ne 's' and fisi ne 'b' then goto,l7
endif

l8 :
if scps eq 'f' then begin
print,'Do you want a ps or an eps file (ps/eps)?'
ps='ps' & read,ps
if ps ne 'ps' and ps ne 'eps' then goto,l8
case ps of
'ps' : begin
print,'Enter filename without the .ps extension:'
fn='...' & read,fn
ext='.ps'
end
'eps' : begin
print,'Enter filename without the .eps extension:'
fn='...' & read,fn
ext='.eps'
end
endcase
endif

print,'Generating GAIA PSF. Please wait...'

;if af eq 'f' then begin
;psf=rebin(psf,psfdim(1)*rf,psfdim(2)*3*rf)
;endif else if af eq 'a' then psfdim=size(psf)/[1,3,3*rf,1,1]
;if sc eq 's' then stpsf=psf else if sc eq 'm' then begin
;auxpsf=fltarr(psfdim(1)*rf,psfdim(2)*rf*3)
;stpsf=fltarr(psfdim(1)*rf,psfdim(2)*rf*3)
;for n=0,npsf-1 do begin
;auxpsf=rot(psf,posa(n),1,774,2322,/pivot)
;stpsf=stpsf+auxpsf
;endfor
;endif

if sc eq 's' then stpsf=rebin(psf,psfdim(1)*rf,psfdim(2)*3*rf) else $
if sc eq 'm' then begin
stpsf=fltarr(psfdim(1)*rf,psfdim(2)*rf*3)
if af eq 'a' then begin
for n=0,npsf-1 do begin
stpsf=stpsf+rot(psf,posa(n),1,774,2322,/pivot)
endfor
endif else if af eq 'f' then begin
for n=0,npsf-1 do begin
stpsf=stpsf+rot(rebin(psf,psfdim(1)*rf,psfdim(2)*3*rf),posa(n),1,774,2322,/pivot)
endfor
endif
endif
; Note that rebinning is performed only if the PSF is entered as a file,
;  meaning that if the PSF is entered as an IDL array, it is assumed that
;  this is a (psfdim(1)*rf,psfdim(2)*3*rf) array, obtained i.e. from a
;  (517,517) PSF of the kind specified in SAG_LL_025 via rebinning
;  (of a factor rf along the x axis and 3*rf along the y axis)

; The PSF is now sampled with a step/rf step along both directions

x1=(psfdim(1)-1)/2*rf-int*rf/2
x2=(psfdim(1)-1)/2*rf+int*rf/2
y1=(psfdim(2)-1)/2*3*rf-int*rf/2
y2=(psfdim(2)-1)/2*3*rf+int*rf/2
stpsf=stpsf(x1:x2,y1:y2)
stdim=size(stpsf)
stpsfmag=2.5*alog10(stpsf)
stpsfmag=max(stpsfmag)-stpsfmag

; Radially Integrated PSF
xc=(findgen(stdim(1))-(stdim(1)-1)/2)*rstep
r=findgen(rint+1)*rstep
aux=fltarr(stdim(1))+1.
;xco=xc#aux
;yco=aux#xc
;cco=complex(xco,yco)
cco=complex(xc#aux,aux#xc)
intpsf=fltarr(rint+1)
for n=0,rint do begin
index=where(abs(cco) lt r(n),count)
if count eq 0 then intpsf(n)=0 else intpsf(n)=total(stpsf(index))
endfor
intpsf=intpsf/total(stpsf)

as=1
; axis style

if scps eq 's' then begin

set_plot,'x'

window,!d.window+1,xsize=500,ysize=500
contour,stpsfmag,xc,xc,levels=[0.5,1.5,2.5,3.5,4.5,5.5,6.5],$
 c_labels=intarr(7)+1,xrange=ran,yrange=ran,xstyle=as,ystyle=as,$
 xticks=1,yticks=1,xtickname=[' ',' '],ytickname=[' ',' '],$
 title='2-D PSF (flux/arcsec^2)',xtitle='mas',ytitle='mas',$
 xmargin=[5,5],ymargin=[5,5]
axis,0,0,xaxis=0,/data,xrange=ran,yrange=ran,ticklen=0.01,$
 xstyle=as,ystyle=as,xmargin=[5,5],ymargin=[5,5]
axis,0,0,yaxis=0,/data,xrange=ran,yrange=ran,ticklen=0.01,$
 xstyle=as,ystyle=as,xmargin=[5,5],ymargin=[5,5]

window,!d.window+1,xsize=500,ysize=500
plot,r,intpsf,title='Encircled Energy',xtitle='r (mas)',$
 ytitle='Normalized Integrated Flux inside r',$
 xrange=intran,yrange=[0,1],xstyle=as,ystyle=as

window,!d.window+1,xsize=500,ysize=500
plot,xc,stpsf(*,int*rf/2)/max(stpsf),$
 title='1-D PSF in the x direction at y=0',$
 xtitle='mas',ytitle='Normalized Flux/arcsec^2',$
 xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 yticks=1,ytickname=[' ',' ']
axis,0,0,yaxis=0,/data,xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 ticklen=0.01

window,!d.window+1,xsize=500,ysize=500
plot,xc,stpsf(int*rf/2,*)/max(stpsf),$
 title='1-D PSF in the y direction at x=0',$
 xtitle='mas',ytitle='Normalized Flux/arcsec^2',$
 xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 yticks=1,ytickname=[' ',' ']
axis,0,0,yaxis=0,/data,xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 ticklen=0.01

endif

if scps eq 'f' then begin

set_plot,'ps'

if ps eq 'eps' then device,/encapsulated

if fisi eq 's' then begin

!p.multi=[0,4,1]

device,filename=fn+ext,xsize=17.,ysize=6.,$
 xoffset=(21.-17.)/2.,yoffset=(29.7-6.)/2.

contour,stpsfmag,xc,xc,levels=[0.5,1.5,2.5,3.5,4.5,5.5,6.5],$
 c_labels=intarr(7),xrange=ran,yrange=ran,xstyle=as,ystyle=as,$
 title='2-D PSF',xtitle='mas',ytitle='mas',$
 xticks=4,yticks=4,charsize=1.5,$
 position=[1.5,1.5,4.5,4.5]/[17,6,17,6]
axis,0,0,xaxis=0,/data,xrange=ran,yrange=ran,ticklen=0.01,$
 xstyle=as,ystyle=as,xticks=1,xtickname=[' ',' ']
axis,0,0,yaxis=0,/data,xrange=ran,yrange=ran,ticklen=0.01,$
 xstyle=as,ystyle=as,yticks=1,ytickname=[' ',' ']

plot,r,intpsf,title='Encircled Energy',xtitle='r (mas)',$
 xrange=intran,yrange=[0,1],xstyle=as,ystyle=as,$
 xticks=4,yticks=4,charsize=1.5,$
 position=[5.5,1.5,8.5,4.5]/[17,6,17,6]

plot,xc,stpsf(*,int*rf/2)/max(stpsf),title='1-D X Axis PSF',$
 xtitle='mas',xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 xticks=4,yticks=4,charsize=1.5,$
 position=[9.5,1.5,12.5,4.5]/[17,6,17,6]
axis,0,0,yaxis=0,/data,xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 yticks=1,ytickname=[' ',' ']

plot,xc,stpsf(int*rf/2,*)/max(stpsf),title='1-D Y Axis PSF',$
 xtitle='mas',xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 xticks=4,yticks=4,charsize=1.5,$
 position=[13.5,1.5,16.5,4.5]/[17,6,17,6]
axis,0,0,yaxis=0,/data,xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 yticks=1,ytickname=[' ',' ']

endif

if fisi eq 'b' then begin

!p.multi=[0,2,2]

device,filename=fn+ext,xsize=19.,ysize=27.7,xoffset=1.,yoffset=1.

yxr=float(!d.y_size)/float(!d.x_size)

contour,stpsfmag,xc,xc,levels=[0.5,1.5,2.5,3.5,4.5,5.5,6.5],$
 c_labels=intarr(7)+1,xrange=ran,yrange=ran,$
 xstyle=as,ystyle=as,xticks=1,yticks=1,$
 xtickname=[' ',' '],ytickname=[' ',' '],$
 title='2-D PSF (flux/arcsec^2)',xtitle='mas',ytitle='mas',$
 position=[0.06,0.9-0.43/yxr,0.49,0.9]
axis,0,0,xaxis=0,/data,xrange=ran,yrange=ran,ticklen=0.01,$
 xstyle=as,ystyle=as
axis,0,0,yaxis=0,/data,xrange=ran,yrange=ran,ticklen=0.01,$
 xstyle=as,ystyle=as

plot,r,intpsf,title='Encircled Energy',xtitle='r (mas)',$
 ytitle='Normalized Integrated Flux inside r',$
 xrange=intran,yrange=[0,1],xstyle=as,ystyle=as,$
 position=[0.56,0.9-0.43/yxr,0.99,0.9]

plot,xc,stpsf(*,int*rf/2)/max(stpsf),$
 title='1-D PSF in the x direction at y=0',$
 xtitle='mas',ytitle='Normalized Flux/arcsec^2',$
 xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 yticks=1,ytickname=[' ',' '],position=[0.06,0.1,0.49,0.1+0.43/yxr]
axis,0,0,yaxis=0,/data,xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 ticklen=0.01

plot,xc,stpsf(int*rf/2,*)/max(stpsf),$
 title='1-D PSF in the y direction at x=0',$
 xtitle='mas',ytitle='Normalized Flux/arcsec^2',$
 xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 yticks=1,ytickname=[' ',' '],position=[0.56,0.1,0.99,0.1+0.43/yxr]
axis,0,0,yaxis=0,/data,xrange=ran,yrange=[0,1],xstyle=as,ystyle=as,$
 ticklen=0.01

endif

device,/close

!p.multi=0

set_plot,'x'

endif

end