G.3 STACK_TRIS.PRO

This program stacks a set of GAIA BBP simulated observations generated by SIM_STACK.PRO into a flux map using the drizzling technique mentioned at the end of Section 6.4 and described in Appendix F.

  1 ;   PROGRAM               STACK_TRIS.PRO
  2 
  3 ; Last Updated             09 Jun 2000         Mattia Vaccari
  4 
  5 ; This program stacks a set of GAIA BBP simulated observations of an HST WFPC2
  6 ;  field generated by SIM_STACK.PRO into a GAIA flux map.
  7 ; The technique used in the stacking is inspired by the drizzling technique
  8 ;  devised for HST Deep Field North (Williams et al. 1996)
  9 
 10 lsc:
 11 print,'Do you want to display GAIA flux map on screen (y/n)?'
 12 sc='...' & read,sc
 13 if sc ne 'y' and sc ne 'n' then begin
 14     print,'Invalid choice!'
 15     goto,lsc
 16     endif
 17 
 18 lps:
 19 print,'Do you want to save the flux map in magnitudes to an eps file (y/n)?'
 20 ps='...' & read,ps
 21 if ps ne 'y' and ps ne 'n' then begin
 22     print,'Invalid choice!'
 23     goto,lps
 24     endif
 25 if ps eq 'y' then begin
 26     print,'Enter filename without the .eps extension:'
 27     fn='fluxmap.eps' & read,fn
 28     print,'Enter the side of the image in centimetres:'
 29     fms=16. & read,fms
 30     endif
 31 
 32 lst:
 33 print,'Do you want to save the important data into an IDL structure (y/n)?'
 34 st='...' & read,st
 35 if st ne 'y' and st ne 'n' then begin
 36     print,'Invalid choice!'
 37     goto,lst
 38     endif
 39 if st eq 'y' then begin
 40     lstfile:
 41     print,'Do you want to save this structure into a dat file (y/n)?'
 42     stfile='...' & read,stfile
 43     if stfile ne 'y' and stfile ne 'n' then begin
 44         print,'Invalid choice!'
 45         goto,lstfile
 46         endif
 47     if stfile eq 'y' then begin
 48         print,'Enter filename without the .dat extension:'
 49         simfn='sim_data.dat' & read,simfn
 50         endif
 51     endif
 52 
 53 ; RETRIEVE DATA FROM SIMSTR STRUCTURE
 54 dadsfile=simstr.dadsfile & date=simstr.date & instrume=simstr.instrume
 55 camerastr=simstr.camerastr &targname=simstr.targname & ra_targ=simstr.ra_targ
 56 dec_targ=simstr.dec_targ & pa_v3=simstr.pa_v3 & filtnam1=simstr.filtnam1
 57 filtnam2=simstr.filtnam2 & filtnu=simstr.filtnu & atodgain=simstr.atodgain
 58 gain=simstr.gain & photflam=simstr.photflam & photzpt=simstr.photzpt
 59 calpar=simstr.calpar & vi=simstr.vi & dm=simstr.dm & hpix=simstr.hpix
 60 exptime=simstr.exptime & hstdnc=simstr.hstdnc
 61 ; hstelc=simstr.hstelc & hstmagins=simstr.hstmagins
 62 muback=simstr.muback & parhstdn=simstr.parhstdn & hstmagsta=simstr.hstmagsta
 63 hstmag=simstr.hstmag & rdim=simstr.rdim & hstmagcen=simstr.hstmagcen
 64 ss=simstr.ss & ssstr=simstr.ssstr & nobs=simstr.nobs & sd=simstr.sd
 65 posa=simstr.posa & cx=simstr.cx & cy=simstr.cy & gaiaet=simstr.gaiaet
 66 rn=simstr.rn & spfx=simstr.spfx & spfy=simstr.spfy & obs=simstr.obs
 67 ssfx=simstr.ssfx & ssfy=simstr.ssfy & gstep=simstr.gstep & gsize=simstr.gsize
 68 parn=simstr.parn & px=simstr.px & py=simstr.py & psx=simstr.psx
 69 psy=simstr.psy & sx=simstr.sx & sy=simstr.sy & xi=simstr.xi & yi=simstr.yi
 70 
 71 print,'Stacking GAIA ',ssstr,' pixels/sample BBP simulated observations.'
 72 print,'Please wait...'
 73 
 74 ; ARRAYS' CREATION
 75 sscx=reform(transpose(((findgen(xi)-0.5*(xi-1))*sx)#replicate(1,ssfx)),xi*ssfx)+$
 76 reform(((findgen(ssfx)-0.5*(ssfx-1))*(0.5*sx/ssfx))#replicate(1,xi),xi*ssfx)
 77 sscy=reform(transpose(((findgen(yi)-0.5*(yi-1))*sy)#replicate(1,ssfy)),yi*ssfy)+$
 78 reform(((findgen(ssfy)-0.5*(ssfy-1))*(0.5*sy/ssfy))#replicate(1,yi),yi*ssfy)
 79 ; sscx and sscy contain the coordinates of subsamples' centers
 80 ;  in the observation FOR
 81 fmgrid=(findgen(gsize)-0.5*(gsize-1))*gstep
 82 ; Coordinates of the centers of the flux map elements in the HST FOR
 83 pfm=fltarr(gsize,gsize)
 84 ; pfm will contain the partial flux map
 85 nfm=fltarr(nobs,gsize,gsize)
 86 ; nfm will contain the nobs partial flux maps
 87 gfm=fltarr(gsize,gsize)
 88 ; gfm will contain the global flux map
 89 
 90 for n=0,nobs-1 do begin
 91     
 92     ; RETRIEVE NTH OBSERVATION
 93     ob=reform(obs(n,*,*))
 94     
 95     ; SUBSAMPLING OF GAIA OBSERVATION
 96     ssob=rebin(ob,xi*ssfx,yi*ssfy,/sample)/(ssfx*ssfy)
 97     ; subsampling of each sample into ssfx*ssfy subsamples
 98     
 99     ; COORDINATES OF SUBSAMPLES' CENTERS IN THE HST FOR
100     auxx   =+(sscx#(fltarr(yi*ssfy)+1.))*cos(posa(n))$
101     -((fltarr(xi*ssfx)+1.)#sscy)*sin(posa(n))$
102     +cx(n)
103     auxy   =+(sscx#(fltarr(yi*ssfy)+1.))*sin(posa(n))$
104     +((fltarr(xi*ssfx)+1.)#sscy)*cos(posa(n))$
105     +cy(n)
106     
107     ; REBINNING OF GAIA PARTIAL FLUX MAP
108     for i=0,gsize-1 do begin
109         xin=where(abs(auxx-fmgrid(i)) lt gstep/2.,wrx)
110         for j=0,gsize-1 do begin
111             if wrx ne 0 then yin=where(abs(auxy(xin)-fmgrid(j)) lt gstep/2.,wry)$
112             else wry=0
113             if wry ne 0 then pfm(i,j)=total(ssob(xin(yin))) else pfm(i,j)=0.
114             endfor
115         endfor
116     ; Rebinning of GAIA subsamples into GAIA partial flux map
117     
118     ; NTH PARTIAL FLUX MAP
119     nfm(n,*,*)=pfm
120     
121     endfor
122 
123 ; GAIA GLOBAL FLUX MAP
124 gfm=total(nfm,1)
125 
126 ; FLUX MAP IN MAGNITUDES
127 gfmmag=-2.5*alog10(gfm>parn)
128 gfmmax=max(gfmmag,min=gfmmin)
129 ; GAIA flux map in uncalibrated magnitudes
130 pargfmdn=parn/gain
131 gfmdnc=gfm/gain
132 ; GAIA flux map in data number counts
133 gfmmagsta=-2.5*alog10((gfmdnc>pargfmdn)/(gaiaet*nobs*(sx*sy/1000.)^2))+$
134 calpar(0,filtnu)+calpar(1,filtnu)*vi+calpar(2,filtnu)*vi^2+dm+0.1
135 ; GAIA flux map in stamag standard magnitudes
136 
137 ; FLUX MAP ON SCREEN
138 ;if sc eq 'y' then begin
139 ;mag=2.
140 ;dispsize=mag*gsize
141 ;set_plot,'x'
142 ;loadct,0
143 ;window,!d.window+1,xsize=dispsize,ysize=dispsize
144 ;pfx=[0,1,1,0] & pfy=[0,0,1,1]
145 ;for i=0,gsize-1 do for j=0,gsize-1 do begin
146 ;polyfill,mag*(pfx+i),mag*(pfy+j),color=round((gfmmag(i,j)-gfmmin)/$
147 ; (gfmmax-gfmmin)*255.),/device
148 ;endfor
149 ;endif
150 
151 if sc eq 'y' then begin
152     mag=2.
153     dispsize=mag*gsize
154     gfmmagsc=rebin(gfmmag,dispsize,dispsize,/sample)
155     set_plot,'x'
156     loadct,0
157     window,!d.window+1,xsize=dispsize,ysize=dispsize
158     tvscl,gfmmagsc
159     endif
160 
161 ; SAVE FLUX MAP INTO AN EPS FILE
162 if ps eq 'y' then begin
163     set_plot,'ps'
164     device,filename=fn+'.eps',xsize=fms,ysize=fms,xoffset=(21.-fms)/2.,$
165     yoffset=(29.7-fms)/2.,/encapsulated,bits_per_pixel=8
166     tvscl,gfmmag,xsize=fms,ysize=fms,/centimeters
167     device,/close
168     set_plot,'x'
169     endif
170 ; Note that the image is centered on an a4 page
171 
172 ; SAVE IMPORTANT DATA INTO AN IDL STRUCTURE
173 if st eq 'y' then begin
174     ; Simulation Structure
175     result=execute('simstr_tris={dadsfile:dadsfile,date:date,instrume:instrume'+$
176     ',camerastr:camerastr,targname:targname,ra_targ:ra_targ,dec_targ:dec_targ'+$
177     ',pa_v3:pa_v3,filtnam1:filtnam1,filtnam2:filtnam2,filtnu:filtnu'+$
178     ',atodgain:atodgain,gain:gain,photflam:photflam,photzpt:photzpt'+$
179     ',calpar:calpar,vi:vi,dm:dm,hpix:hpix,exptime:exptime,hstdnc:hstdnc'+$
180     ; ',hstelc:hstelc,hstmagins:hstmagins'+$
181     ',muback:muback,parhstdn:parhstdn,hstmagsta:hstmagsta,hstmag:hstmag'+$
182     ',rdim:rdim,hstmagcen:hstmagcen'+$
183     ',ss:ss,ssstr:ssstr,nobs:nobs,sd:sd,posa:posa,cx:cx,cy:cy,gaiaet:gaiaet'+$
184     ',rn:rn,spfx:spfx,spfy:spfy,obs:obs'+$
185     ',ssfx:ssfx,ssfy:ssfy,gstep:gstep,gsize:gsize,nfm:nfm,gfm:gfm,parn:parn'+$
186     ',gfmmag:gfmmag,gfmmagsta:gfmmagsta'+$
187     ',px:px,py:py,psx:psx,psy:psy,sx:sx,sy:sy,xi:xi,yi:yi}')
188     endif
189 
190 ; SAVE IMPORTANT DATA INTO A DAT FILE
191 if stfile eq 'y' then begin
192     save,filename=simfn+'.dat',/verbose,simstr_tris
193     endif
194 
195 end