;   PROGRAM                 STACK_BIS.PRO

; Last Updated               07 Jun 2000           Mattia Vaccari

; This program stacks a set of GAIA BBP simulated observations of an HST WFPC2
;  field generated by SIM_STACK.PRO into a GAIA flux map.

lsc:
print,'Do you want to display GAIA flux map on screen (y/n)?'
sc='...' & read,sc
if sc ne 'y' and sc ne 'n' then begin
print,'Invalid choice!'
goto,lsc
endif

lps:
print,'Do you want to save the flux map in magnitudes to an eps file (y/n)?'
ps='...' & read,ps
if ps ne 'y' and ps ne 'n' then begin
print,'Invalid choice!'
goto,lps
endif
if ps eq 'y' then begin
print,'Enter filename without the .eps extension:'
fn='fluxmap.eps' & read,fn
print,'Enter the side of the image in centimetres:'
fms=16. & read,fms
endif

lst:
print,'Do you want to save the important data into an IDL structure (y/n)?'
st='...' & read,st
if st ne 'y' and st ne 'n' then begin
print,'Invalid choice!'
goto,lst
endif
if st eq 'y' then begin
lstfile:
print,'Do you want to save this structure into a dat file (y/n)?'
stfile='...' & read,stfile
if stfile ne 'y' and stfile ne 'n' then begin
print,'Invalid choice!'
goto,lstfile
endif
if stfile eq 'y' then begin
print,'Enter filename without the .dat extension:'
simfn='sim_data.dat' & read,simfn
endif
endif

; RETRIEVE DATA FROM SIMSTR STRUCTURE
dadsfile=simstr.dadsfile & date=simstr.date & instrume=simstr.instrume
camerastr=simstr.camerastr &targname=simstr.targname & ra_targ=simstr.ra_targ
dec_targ=simstr.dec_targ & pa_v3=simstr.pa_v3 & filtnam1=simstr.filtnam1
filtnam2=simstr.filtnam2 & filtnu=simstr.filtnu & atodgain=simstr.atodgain
gain=simstr.gain & photflam=simstr.photflam & photzpt=simstr.photzpt
calpar=simstr.calpar & vi=simstr.vi & dm=simstr.dm & hpix=simstr.hpix
exptime=simstr.exptime & hstdnc=simstr.hstdnc
; hstelc=simstr.hstelc & hstmagins=simstr.hstmagins
muback=simstr.muback & parhstdn=simstr.parhstdn & hstmagsta=simstr.hstmagsta
hstmag=simstr.hstmag & rdim=simstr.rdim & hstmagcen=simstr.hstmagcen
ss=simstr.ss & ssstr=simstr.ssstr & nobs=simstr.nobs & sd=simstr.sd
posa=simstr.posa & cx=simstr.cx & cy=simstr.cy & gaiaet=simstr.gaiaet
rn=simstr.rn & spfx=simstr.spfx & spfy=simstr.spfy & obs=simstr.obs
ssfx=simstr.ssfx & ssfy=simstr.ssfy & gstep=simstr.gstep & gsize=simstr.gsize
parn=simstr.parn & px=simstr.px & py=simstr.py & psx=simstr.psx
psy=simstr.psy & sx=simstr.sx & sy=simstr.sy & xi=simstr.xi & yi=simstr.yi

print,'Stacking GAIA ',ssstr,' pixels/sample BBP simulated observations.'
print,'Please wait...'

; FLUX MAP GRID IN THE HST FOR
fmgrid=fltarr(2,gsize^2)
k=0l
for i=0,gsize-1 do for j=0,gsize-1 do begin
fmgrid(*,k)=[i-(gsize-1)/2,j-(gsize-1)/2]*gstep
k=k+1
endfor
; fmgrid is a (2,gsize^2) array containing along its 2 columns the x and y
;  coordinates of the gsize^2 flux map grid points in the HST FOR.
;  The first row contains the coordinate of the point at the lower left corner
;  of the grid and the two columns are ordered first bottom/up and then
;  left/right (which is the order in which IDL dipslays 2-D arrays).
;  The essentially 2-D arrays folded into 1-D arrays that will be used in the
;  following for the sake of efficiency will maintain the same ordering.

; FLUX MAP GRID IN THE OBSERVATION FOR
ones=replicate(1.,gsize^2)
;obsfmx=[[+cos(posa)],[+sin(posa)]]#fmgrid+$
; diag([[+cos(posa)],[+sin(posa)]]#transpose([[cx],[cy]]))#ones
;obsfmy=[[-sin(posa)],[+cos(posa)]]#fmgrid+$
; diag([[-sin(posa)],[+cos(posa)]]#transpose([[cx],[cy]]))#ones
;obsfmx=[[+cos(posa)],[-sin(posa)]]#fmgrid+cx#transpose(ones)
;obsfmy=[[+sin(posa)],[+cos(posa)]]#fmgrid+cy#transpose(ones)
obsfmx=[[+cos(posa)],[+sin(posa)]]#fmgrid-$
 diag([[+cos(posa)],[+sin(posa)]]#transpose([[cx],[cy]]))#ones
obsfmy=[[-sin(posa)],[+cos(posa)]]#fmgrid-$
 diag([[-sin(posa)],[+cos(posa)]]#transpose([[cx],[cy]]))#ones
; obsfmx and obsfmy are (nobs,gsize^2) arrays containing the x and y
;  coordinates of the flux map elements in the observation for

; OBSERVATIONS' REBINNING
obsreb=rebin(obs,nobs,xi*ssfx,yi*ssfy,/sample)/(ssfx*ssfy)

; ASSIGNMENT OF A SUBSAMPLE VALUE TO EACH FLUX MAP ELEMENT
xin=round(obsfmx/(sx/ssfx)+0.5*(xi*ssfx-1))
yin=round(obsfmy/(sy/ssfy)+0.5*(yi*ssfy-1))
; xin and yin express obsfmx and obsfmy in samples
nfm=fltarr(nobs,gsize^2)
for n=0,nobs-1 do nfm(n,*)=obsreb([intarr(gsize^2)+n],[xin(n,*)],[yin(n,*)])
;nfm=obsreb(indgen(nobs)#replicate(1,gsize^2),$
; round(obsfmx/(sx/ssfx)+0.5*(xi*ssfx-1)),$
; round(obsfmy/(sy/ssfy)+0.5*(yi*ssfy-1)))
; nfm now contains along its nobs columns the sample value of the sample
;  whose centre is nearest to each flux map element for each observation
; In other words, nfm contains along its nobs columns the nobs
;  one-dimensional partial flux maps

nfm=transpose(reform(nfm,nobs,gsize,gsize),[0,2,1])
; nfm now contains the nobs two-dimensional partial flux maps
gfm=total(nfm,1)
; gfm is the two-dimensional GAIA glaobal flux map

; FLUX MAP IN MAGNITUDES
gfmmag=-2.5*alog10(gfm>parn)
gfmmax=max(gfmmag,min=gfmmin)
; GAIA flux map in uncalibrated magnitudes
pargfmdn=parn/gain
gfmdnc=gfm/gain
; GAIA flux map in data number counts
gfmmagsta=-2.5*alog10((gfmdnc>pargfmdn)/(gaiaet*nobs*(sx*sy/1000.)^2))+$
 calpar(0,filtnu)+calpar(1,filtnu)*vi+calpar(2,filtnu)*vi^2+dm+0.1
; GAIA flux map in stamag standard magnitudes

; FLUX MAP ON SCREEN
;if sc eq 'y' then begin
;mag=2.
;dispsize=mag*gsize
;set_plot,'x'
;loadct,0
;window,!d.window+1,xsize=dispsize,ysize=dispsize
;pfx=[0,1,1,0] & pfy=[0,0,1,1]
;for i=0,gsize-1 do for j=0,gsize-1 do begin
;polyfill,mag*(pfx+i),mag*(pfy+j),color=round((gfmmag(i,j)-gfmmin)/$
; (gfmmax-gfmmin)*255.),/device
;endfor
;endif

if sc eq 'y' then begin
mag=2.
dispsize=mag*gsize
gfmmagsc=rebin(gfmmag,dispsize,dispsize,/sample)
set_plot,'x'
loadct,0
window,!d.window+1,xsize=dispsize,ysize=dispsize
tvscl,gfmmagsc
endif

; SAVE FLUX MAP INTO AN EPS FILE
if ps eq 'y' then begin
set_plot,'ps'
device,filename=fn+'.eps',xsize=fms,ysize=fms,xoffset=(21.-fms)/2.,$
 yoffset=(29.7-fms)/2.,/encapsulated,bits_per_pixel=8
tvscl,gfmmag,xsize=fms,ysize=fms,/centimeters
device,/close
set_plot,'x'
endif
; Note that the image is centered on an a4 page

; SAVE IMPORTANT DATA INTO AN IDL STRUCTURE
if st eq 'y' then begin
; Simulation Structure
result=execute('simstr_bis={dadsfile:dadsfile,date:date,instrume:instrume'+$
 ',camerastr:camerastr,targname:targname,ra_targ:ra_targ,dec_targ:dec_targ'+$
 ',pa_v3:pa_v3,filtnam1:filtnam1,filtnam2:filtnam2,filtnu:filtnu'+$
 ',atodgain:atodgain,gain:gain,photflam:photflam,photzpt:photzpt'+$
 ',calpar:calpar,vi:vi,dm:dm,hpix:hpix,exptime:exptime,hstdnc:hstdnc'+$
; ',hstelc:hstelc,hstmagins:hstmagins'+$
 ',muback:muback,parhstdn:parhstdn,hstmagsta:hstmagsta,hstmag:hstmag'+$
 ',rdim:rdim,hstmagcen:hstmagcen'+$
 ',ss:ss,ssstr:ssstr,nobs:nobs,sd:sd,posa:posa,cx:cx,cy:cy,gaiaet:gaiaet'+$
 ',rn:rn,spfx:spfx,spfy:spfy,obs:obs'+$
 ',ssfx:ssfx,ssfy:ssfy,gstep:gstep,gsize:gsize,nfm:nfm,gfm:gfm,parn:parn'+$
 ',gfmmag:gfmmag,gfmmagsta:gfmmagsta'+$
 ',px:px,py:py,psx:psx,psy:psy,sx:sx,sy:sy,xi:xi,yi:yi}')
endif

; SAVE IMPORTANT DATA INTO A DAT FILE
if stfile eq 'y' then begin
save,filename=simfn+'.dat',/verbose,simstr_bis
endif

end