G.1 SIM_STACK.PRO

This program simulates a given number of GAIA BBP simulated observations of a given HST WFPC2 image and stacks them into a flux map following the procedures described in Sections 6.3 and 6.4, respectively.

  1 ;   PROGRAM             SIM_STACK.PRO
  2 
  3 ; Last Updated           05 Jun 2000          Mattia Vaccari
  4 
  5 ; This program generates a set of GAIA BBP simulated observations
  6 ;  of an HST WFPC2 field and stacks them into a GAIA flux map.
  7 
  8 ; PARAMETERS
  9 hpix=[45.5,96.6]                     ; HST pixel in mas (PC and WFCs)
 10 spfx=2l                              ; Subpixeling factors
 11 spfy=2l
 12 px=37.2                              ; GAIA pixel in mas
 13 py=111.6
 14 psx=[6,6,6,1,4,4]                    ; GAIA sample in pixels
 15 psy=[8,4,2,8,8,4]
 16 sx=psx*px                            ; GAIA sample in mas
 17 sy=psy*py
 18 obsside=[24998.4,55353.6]
 19 ; Side of GAIA observation obtained from a WPPC2 field in mas
 20 ;  (PC and WFCs respectively)
 21 xi=round(obsside#(1./sx))
 22 yi=round(obsside#(1./sy))
 23 ; GAIA observation in samples for different combinations of
 24 ;  WFPC2 cameras and BBP sample sizes
 25 gaiaet=0.8618
 26 ; GAIA observation exposure time in s
 27 rn=[5.44289,6.74395,8.78593,10.4357,6.12805,7.83178]
 28 ; GAIA BBP readnoise in e-/sample rms, calculated following SAG_MV_04
 29 ;  assuming full CCD readout
 30 ssfx=psx                             ; Subsampling factors
 31 ssfy=psy*3
 32 ; The subsamples are squares of side 37.2 mas, which is also the sample size
 33 ;  in the along scan direction
 34 gstep=hpix/[1.,2.]                   ; GAIA flux map step in mas
 35 ; gstep equals 45.5 and 48.3 mas for PC and WFCs respectively
 36 gsize=[350l,320l]
 37 ; Side of GAIA flux map in steps for PC and WFCs respectively
 38 ; The side of GAIA flux map is thus 15925 and 15456 mas for PC and WFCs
 39 ;  respectively
 40 
 41 print,'The following .fits.gz files are available in the ./hstdata directory:'
 42 spawn,'cd hstdata ; ls -C1 *.fits.gz | awk ''{print $1}''',fitsfiles,$
 43 count=n_fits
 44 
 45 for i=0,n_fits-1 do begin
 46     print,i+1,strmid(fitsfiles(i),0,strlen(fitsfiles(i))-8),format='(i3,6x,a)'
 47     endfor
 48 
 49 print,'Enter the number of the file you want to open:'
 50 fnu=0 & read,fnu
 51 print,'Opening file number ',fnu
 52 hstdnc=readfits('hstdata/'+fitsfiles(fnu-1),header,/silent)
 53 himdim=size(hstdnc)
 54 
 55 f=n_elements(header)
 56 
 57 for j=0,f-2 do begin
 58     r=header(j)
 59     i=strpos(r,'DADSFILE=')
 60     if i ne -1 then d=execute(strcompress(strmid(r,0,30),/remove_all)) else begin
 61         i=strpos(r,'ROOTNAME=')
 62         if i ne -1 then begin
 63             d=execute(strcompress(strmid(r,0,30),/remove_all))
 64             DADSFILE=ROOTNAME
 65             endif
 66         endelse
 67     i=strpos(r,'INSTRUME=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 68     /remove_all))
 69     i=strpos(r,'DATE    =')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 70     /remove_all))
 71     i=strpos(r,'TARGNAME=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 72     /remove_all))
 73     i=strpos(r,'RA_TARG =')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 74     /remove_all))
 75     i=strpos(r,'DEC_TARG=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 76     /remove_all))
 77     i=strpos(r,'PA_V3   =')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 78     /remove_all))
 79     i=strpos(r,'FILTNAM1=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 80     /remove_all))
 81     i=strpos(r,'FILTNAM2=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 82     /remove_all))
 83     i=strpos(r,'EXPTIME =')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 84     /remove_all))
 85     i=strpos(r,'ATODGAIN=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 86     /remove_all))
 87     i=strpos(r,'PHOTFLAM=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 88     /remove_all))
 89     i=strpos(r,'PHOTZPT =')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
 90     /remove_all))
 91     endfor
 92 
 93 print,'DADSFILE',dadsfile,format='(a,15x,a)'
 94 print,'INSTRUMENT',instrume,format='(a,13x,a)'
 95 print,'OBSERVATION DATE',date,format='(a,7x,a)'
 96 print,'TARGET NAME',targname,format='(a,12x,a)'
 97 print,'RIGHT ASCENSCION',ra_targ,format='(a,1x,a)'
 98 print,'DECLINATION',dec_targ,format='(a,6x,a)'
 99 print,'V3 POSITION ANGLE',pa_v3
100 print,'FIRST FILTER NAME',filtnam1,format='(a,6x,a)'
101 print,'SECOND FILTER NAME',filtnam2,format='(a,5x,a)'
102 print,'EXPOSURE TIME',exptime,format='(a,4x,a)'
103 print,'ATODGAIN',atodgain,format='(a,9x,a)'
104 
105 ; HST DATA NUMBER COUNTS
106 if himdim(0) eq 2 then begin
107     print,'This file contains one image.'
108     imagenu=1
109     endif
110 
111 if himdim(0) eq 3 then begin
112     print,'This file contains ',himdim(3),' images. Which do you want to take?'$
113     ,format='(a,i1,a)'
114     imagenu=1 & read,imagenu
115     endif
116 
117 imagenu=imagenu-1
118 hsti=imagenu<1
119 hpix=hpix(hsti)
120 camerastr=['PC','WFC2','WFC3','WFC4'] & camerastr=camerastr(imagenu)
121 print,'Creating HST WFPC2 '+camerastr+' image.'
122 hstdnc=hstdnc(*,*,imagenu)
123 ; Note that if the .fits.gz file containing the original data consists of one
124 ;  image only, the program assumes that this is a PC image. Otherwise, it
125 ;  assumes the file is in the format used by the HDA for the WFPC2 data files,
126 ;  i.e. that it contains, in this order the four images obtained with the
127 ;  PC, WFC2, WFC3 and WFC4.
128 
129 ; A/D GAIN
130 gain=[[7.12,7.12,6.90,7.10],[13.99,14.50,13.95,13.95]]
131 gain=gain(imagenu,round(atodgain/7.)-1)
132 ; HST WFPC2 A/D conversion gain, from "WFPC2 Instrument Handbook",
133 ;  Version 4.0, page 82
134 
135 lss:
136 print,'What BBP sample size do you want to use?'
137 print,'1=6x8  2=6x4  3=6x2  4=1x8  5=4x8  6=4x4'
138 ss=1 & read,ss
139 if total([1.,2.,3.,4.,5.,6.] eq ss) eq 0 then begin
140     print,'Invalid choice!'
141     goto,lss
142     endif
143 ss=ss-1
144 
145 ; PARAMETERS' CHOICE
146 psx=psx(ss)
147 psy=psy(ss)
148 sx=sx(ss)
149 sy=sy(ss)
150 xi=xi(hsti,ss)
151 yi=yi(hsti,ss)
152 rn=rn(ss)
153 ssfx=ssfx(ss)
154 ssfy=ssfy(ss)
155 gstep=gstep(hsti)
156 gsize=gsize(hsti)
157 
158 ; HST IMAGE PHOTOMETRIC CALIBRATION
159 vi=1.2
160 ; Average V-I Galaxy Color Index from Prugniel and H\'eraudeau 1998
161 
162 calpar=[[21.729,-0.051,+0.009],$ ; F555W
163 [22.093,+0.254,+0.012],$ ; F606W
164 [20.920,-0.124,+0.028]]  ; F814W
165 
166 case filtnam1 of
167     'F555W' : begin
168         filtnu=0
169         stamag='V'
170         end
171     'F606W' : begin
172         filtnu=1
173         stamag='V'
174         end
175     'F814W' : begin
176         filtnu=2
177         stamag='I'
178         end
179     else : begin
180         print,'The ',filtnam1,' filter is not supported! ',$
181         'Using F555W-V transformations!'
182         filtnu=0
183         stamag='V'
184         end
185     endcase
186 
187 case stamag of
188     'V' : muback=22.06
189     'I' : muback=21.46
190     endcase
191 ; Sky background in mag/arcsec^2. From SAG_GG_05.
192 ; Note that the choice of this parameter does not affect in any way
193 ;  the processing of the data, but their representation only
194 
195 gr=[1.987,2.003,2.006,1.955]
196 ; Gain ratio from Holtzman et al. 1995b
197 dm=2.5*alog10(gr)
198 if atodgain eq 7. then dm=dm(imagenu) else dm=0.
199 ; Since most photometric calibration observations were done with atodgain=15,
200 ;  there is a correction dm to be applied for exposures carried out with
201 ;  atodgain=7
202 
203 ; HST IMAGE IN UNCALIBRATED MAGNITUDES
204 parhstdn=(exptime*(hpix/1000.)^2)*10^(0.4*(calpar(0,filtnu)+$
205 calpar(1,filtnu)*vi+calpar(2,filtnu)*vi^2+dm+0.1-muback))
206 ; Sky background parameter in data numbers
207 ;parhstdn=4.0
208 ; An alternative sky background parameter can be used to obtain a better
209 ;  image visibility (i.e. an higher contrast). In the case of m100_900,
210 ;  for instance, only 132 out of 122500 pixels (0.1%) in the [350,350] central
211 ;  region of the PC image "contained" less than 4 data numbers.
212 hstmag=-2.5*alog10(hstdnc>parhstdn)
213 rdim=[350l,160l] & rdim=rdim(hsti)
214 hstmagcen=hstmag((himdim(1)-rdim)/2:(himdim(1)+rdim)/2-1,$
215 (himdim(2)-rdim)/2:(himdim(2)+rdim)/2-1)
216 ; hstmagcen covers the central part of hstmag corresponding to a square area
217 ;  of side 15925 and 15456 mas, for PC and WFC respectively, i.e. exactly the
218 ;  portion of the HST image that is covered by GAIA flux map.
219 
220 ; HST IMAGE IN WFPC2 INSTRUMENTAL MAGNITUDES
221 ; See "WFPC2 Instrument Handbook", Version 4.0, Section 8.7
222 ;zp=-2.5*alog10(photflam)+photzpt
223 ; Zero-point of the "FILTNAM1" magnitude scale
224 ;hstmagins=-2.5*alog10((hstdnc>parhstdn)/(exptime*(hpix/1000.)^2))+zp
225 ; HST image in "FILTNAM1" magnitudes
226 
227 ; HST IMAGE IN STANDARD MAGNITUDES
228 hstmagsta=-2.5*alog10((hstdnc>parhstdn)/(exptime*(hpix/1000.)^2))+$
229 calpar(0,filtnu)+calpar(1,filtnu)*vi+calpar(2,filtnu)*vi^2+dm+0.1
230 ; HST image in stamag standard magnitudes
231 ; The transformations between WFPC2 Data Number counts to standard magnitudes
232 ;  are carried out using Equation 9 and Table 10 from Holtzman et al. 1995b,
233 ;  with a constant correction of 0.1 mag/arcsec^2 for infinite aperture
234 
235 ; HST ELECTRON COUNTS
236 ;hstelc=gain*hstdnc
237 ; The atodgain keyword can be equal to 7 or 15, but the true value of the
238 ;  A/D conversion gain (given by the variable gain) is actually very
239 ;  near to 7 and 14, respectively, and is different in WFPC2's different
240 ;  detectors. Please note that we do not use the gain value in the
241 ;  conversion to a magnitude scale. Actually, we only use it to give an
242 ;  estimate of GAIA's electron counts, which we have to combine with
243 ;  readnoise estimates (given in electrons) to properly take into
244 ;  account the noise.
245 ;parhstel=gain*parhstdn     ; sky background parameter in electrons
246 ; parhstel=30.
247 ; An alternative sky background parameter can be used to obtain a better
248 ;  image visibility (i.e. an higher contrast). In the case of m100_900,
249 ;  for instance, only 171 out of 122500 pixels (0.14%) in the central region
250 ;  of the PC image "contained" less than 30 electrons.
251 
252 ; GENERATION OF CENTERS AND POSITION ANGLES OF OBSERVATIONS
253 print,'How many simulated observations do you want to generate?'
254 nobs=1 & read,nobs
255 ; if nobs eq 1 then posa=0. & cx=0. & cy=0. else begin ..... endif
256 lsd:
257 print,'How do you want the position angles of the scan directions '+$
258 'to be distributed?'
259 print,'          1=0-360    2=0-45 or 90-135    3=0-90    4=0-45'
260 sd=1 & read,sd
261 if total([1.,2.,3.,4.] eq sd) eq 0 then begin
262     print,'Invalid choice!'
263     goto,lsd
264     endif
265 ; The value of sd goes in increasing order of non-randomness
266 ;  of the scan directions
267 case sd of
268     1 : posa=randomu(seed,nobs)*(2*!pi)
269     2 : begin
270         posa=randomu(seed,nobs)*(!pi/4)
271         posa=posa+((-1)^round(randomu(seed,nobs)*1000000.)*(!pi/4)+(!pi/4))
272         end
273     3 : posa=randomu(seed,nobs)*(!pi/2)
274     4 : posa=randomu(seed,nobs)*(!pi/4)
275     endcase
276 ; posa contains the position angles (measured counterclockwise)
277 ;  of GAIA FORs with respect to the HST FOR
278 cx=(randomu(seed,nobs)-0.5)*sx
279 cy=(randomu(seed,nobs)-0.5)*sy
280 ; cx and cy contain the coordinates of the centers of GAIA FORs
281 ;  (i.e. of GAIA observations) in the HST FOR
282 
283 ; In other words, the GAIA FOR for the nth observation is obtained
284 ;  starting from the HST FOR by a translation of a vector [cx(n),cy(n)]
285 ;  followed by a rotation around the center of the FOR thus obtained
286 ;  of an angle posa(n) measured counterclockwise
287 
288 lsc:
289 print,'Do you want to display GAIA flux map on screen (y/n)?'
290 sc='...' & read,sc
291 if sc ne 'y' and sc ne 'n' then begin
292     print,'Invalid choice!'
293     goto,lsc
294     endif
295 
296 lps:
297 print,'Do you want to save the flux map in magnitudes to an eps file (y/n)?'
298 ps='...' & read,ps
299 if ps ne 'y' and ps ne 'n' then begin
300     print,'Invalid choice!'
301     goto,lps
302     endif
303 if ps eq 'y' then begin
304     print,'Enter filename without the .eps extension:'
305     fn='fluxmap.eps' & read,fn
306     print,'Enter the side of the image in centimetres:'
307     fms=16. & read,fms
308     endif
309 
310 lst:
311 print,'Do you want to save the important data into an IDL structure (y/n)?'
312 st='...' & read,st
313 if st ne 'y' and st ne 'n' then begin
314     print,'Invalid choice!'
315     goto,lst
316     endif
317 if st eq 'y' then begin
318     lstfile:
319     print,'Do you want to save this structure into a dat file (y/n)?'
320     stfile='...' & read,stfile
321     if stfile ne 'y' and stfile ne 'n' then begin
322         print,'Invalid choice!'
323         goto,lstfile
324         endif
325     if stfile eq 'y' then begin
326         print,'Enter filename without the .dat extension:'
327         simfn='sim_data.dat' & read,simfn
328         endif
329     endif
330 
331 ssstr=['6x8','6x4','6x2','1x8','4x8','4x4'] & ssstr=ssstr(ss)
332 print,'Generating and stacking GAIA ',ssstr,' pixels/sample BBP '+$
333 'simulated observations.'
334 print,'Please wait...'
335 
336 ; EXPOSURE TIME SCALING OF HST ELECTRON COUNTS
337 ;hgelc=(gain*gaiaet/exptime)*hstdnc
338 ; HST electron counts scaled to GAIA exposure time
339 ; We assume that the electron counts per unit time per unit area of the two
340 ;  instruments are equal, which is conservative for most combinations of
341 ;  HST-GAIA filters.
342 
343 ; SUBPIXELING OF HST IMAGE
344 sphim=rebin(hstdnc,himdim(1)*spfx,himdim(2)*spfy,/sample)*$
345 ((gain*gaiaet)/(exptime*spfx*spfy))
346 ; Subpixeling of each HST pixel into 4 subpixels
347 ; The rebin function with the sample keyword performs a nearest neighbour
348 ;  sampling, thus preserving the total number of electrons with great accuracy
349 
350 ; PSF
351 psfname=['68','64','62','18','48','44']
352 psf=readfits('psf/'+'glo_10_'+psfname(ss)+'.fits.gz',/silent)
353 ;psf=readfits('psf/'+psfname+'.fits.gz',/silent)
354 ; (517,517) PSF obtained following SAG_LL_025, using syntfits.f and
355 ; assuming V-I=1.2 and Field Point 10
356 ; The PSF is sampled with a step of 1/4 pixel along both directions
357 ; The total area mapped by the PSF is thus of
358 ;  (517,517)*(1/4)*(px,py)=(4.8081,14.4243) arcsec^2
359 xtot=[21,21,129,33,33,21] & xtot=xtot(ss)
360 xmed=[10,10,64,16,16,10]  & xmed=xmed(ss)
361 dx=[24,24,4,16,16,24]     & dx=dx(ss)
362 xin=(findgen(xtot)-xmed)*dx+258
363 ytot=[33,17,17,33,17,65] & ytot=ytot(ss)
364 ymed=[16,8,8,16,8,32]    & ymed=ymed(ss)
365 dy=[16,32,32,16,32,8]    & dy=dy(ss)
366 yin=(findgen(ytot)-ymed)*dy+258
367 psf=psf(xin,*)
368 psf=psf(*,yin)
369 ; PSF resampling with a step equal to sample size
370 case ss of
371     0 : psf=psf(6:14,10:22)
372     1 : psf=psf(6:14,5:11)
373     2 : psf=psf(40:88,5:11)
374     3 : psf=psf(10:22,10:22)
375     4 : psf=psf(10:22,5:11)
376     5 : psf=psf(6:14,20:44)
377     endcase
378 ; Reduced PSF with a step equal to sample size
379 ; In other words, we truncate the PSF to a central region of about
380 ; (50,50) pixel = (1860,5580) arcsec containing about 99% of the
381 ;  complete PSF (i.e. of the energy)
382 psf=psf/total(psf)
383 ; Normalization
384 
385 ; ARRAYS' CREATION
386 scx=(findgen(xi)-0.5*(xi-1))*sx
387 scy=(findgen(yi)-0.5*(yi-1))*sy
388 ; scx and scy contain the coordinates of samples' centers
389 ;  in the observation FOR
390 spcx=( findgen(himdim(1)*spfx)-0.5*(himdim(1)*spfx-1) )*(hpix/spfx)
391 spcy=( findgen(himdim(2)*spfy)-0.5*(himdim(2)*spfy-1) )*(hpix/spfy)
392 ; spcx and spcy contain the coordinates of subpixels' centers
393 ;  in the HST FOR
394 ob=fltarr(xi,yi)
395 ; ob will contain the nth observation during processing
396 obs=fltarr(nobs,xi,yi)
397 ; obs will contain the nobs observations
398 sscx=(findgen(xi*ssfx)-0.5*(xi*ssfx-1))*(sx/ssfx)
399 sscy=(findgen(yi*ssfy)-0.5*(yi*ssfy-1))*(sy/ssfy)
400 ; sscx and sscy contain the coordinates of subsamples' centers
401 ;  in the observation FOR
402 fmgrid=(findgen(gsize)-0.5*(gsize-1))*gstep
403 ; Coordinates of the centers of the flux map elements in the HST FOR
404 pfm=fltarr(gsize,gsize)
405 ; pfm will contain the partial flux map
406 nfm=fltarr(nobs,gsize,gsize)
407 ; nfm will contain the nobs partial flux maps
408 gfm=fltarr(gsize,gsize)
409 ; gfm will contain the global flux map
410 
411 for n=0,nobs-1 do begin
412     
413     ; COORDINATES OF SUBPIXELS' CENTERS IN THE OBSERVATION FOR
414     auxx   =+((spcx-cx(n))#(fltarr(himdim(2)*spfy)+1.))*cos(posa(n))$
415     +((fltarr(himdim(1)*spfx)+1.)#(spcy-cy(n)))*sin(posa(n))
416     auxy   =-((spcx-cx(n))#(fltarr(himdim(2)*spfy)+1.))*sin(posa(n))$
417     +((fltarr(himdim(1)*spfx)+1.)#(spcy-cy(n)))*cos(posa(n))
418     
419     ; REBINNING OF GAIA OBSERVATION
420     for i=0,xi-1 do begin
421         xin=where(abs(auxx-scx(i)) lt sx/2.,wrx)
422         for j=0,yi-1 do begin
423             if wrx ne 0 then yin=where(abs(auxy(xin)-scy(j)) lt sy/2.,wry)$
424             else wry=0
425             if wry ne 0 then ob(i,j)=total(sphim(xin(yin))) else ob(i,j)=0.
426             endfor
427         endfor
428     ; Rebinning of HST subpixels into GAIA samples
429     
430     ; CONVOLUTION WITH PSF
431     ob=kconvol(ob,psf)
432     ; !!! Note the use of kconvol instead of convol !!!
433     
434     ; NOISE
435     ;obmed=median(ob)
436     ;ob=poidev(ob>obmed,seed=seed)
437     ; The use of the > operator overestimates the photon noise, but avoids
438     ;  possible problems connected with the use of the poidev function with
439     ;  an argument containing mostly very small numbers
440     ob=poidev(ob,seed=seed)
441     ob=ob+randomn(seed,xi,yi)*rn
442     ; Addition of photon noise and readnoise
443     
444     ; NTH OBSERVATION
445     obs(n,*,*)=ob
446     
447     ; SUBSAMPLING OF GAIA OBSERVATION
448     ssob=rebin(ob,xi*ssfx,yi*ssfy,/sample)/(ssfx*ssfy)
449     ; subsampling of each sample into ssfx*ssfy subsamples
450     
451     ; COORDINATES OF SUBSAMPLES' CENTERS IN THE HST FOR
452     auxx   =+(sscx#(fltarr(yi*ssfy)+1.))*cos(posa(n))$
453     -((fltarr(xi*ssfx)+1.)#sscy)*sin(posa(n))$
454     +cx(n)
455     auxy   =+(sscx#(fltarr(yi*ssfy)+1.))*sin(posa(n))$
456     +((fltarr(xi*ssfx)+1.)#sscy)*cos(posa(n))$
457     +cy(n)
458     
459     ; REBINNING OF GAIA PARTIAL FLUX MAP
460     for i=0,gsize-1 do begin
461         xin=where(abs(auxx-fmgrid(i)) lt gstep/2.,wrx)
462         for j=0,gsize-1 do begin
463             if wrx ne 0 then yin=where(abs(auxy(xin)-fmgrid(j)) lt gstep/2.,wry)$
464             else wry=0
465             if wry ne 0 then pfm(i,j)=total(ssob(xin(yin))) else pfm(i,j)=0.
466             endfor
467         endfor
468     ; Rebinning of GAIA subsamples into GAIA partial flux map
469     
470     ; NTH PARTIAL FLUX MAP
471     nfm(n,*,*)=pfm
472     
473     endfor
474 
475 ; GAIA GLOBAL FLUX MAP
476 gfm=total(nfm,1)
477 
478 ; FLUX MAP IN MAGNITUDES
479 parn=(nobs*gaiaet/exptime)*(gstep^2/hpix^2)*(gain*parhstdn)
480 gfmmag=-2.5*alog10(gfm>parn)
481 gfmmax=max(gfmmag,min=gfmmin)
482 ; GAIA flux map in uncalibrated magnitudes
483 pargfmdn=parn/gain
484 gfmdnc=gfm/gain
485 ; GAIA flux map in data number counts
486 gfmmagsta=-2.5*alog10((gfmdnc>pargfmdn)/(gaiaet*nobs*(sx*sy/1000.)^2))+$
487 calpar(0,filtnu)+calpar(1,filtnu)*vi+calpar(2,filtnu)*vi^2+dm+0.1
488 ; GAIA flux map in stamag standard magnitudes
489 
490 ; FLUX MAP ON SCREEN
491 ;if sc eq 'y' then begin
492 ;mag=2.
493 ;dispsize=mag*gsize
494 ;set_plot,'x'
495 ;loadct,0
496 ;window,!d.window+1,xsize=dispsize,ysize=dispsize
497 ;pfx=[0,1,1,0] & pfy=[0,0,1,1]
498 ;for i=0,gsize-1 do for j=0,gsize-1 do begin
499 ;polyfill,mag*(pfx+i),mag*(pfy+j),color=round((gfmmag(i,j)-gfmmin)/$
500 ; (gfmmax-gfmmin)*255.),/device
501 ;endfor
502 ;endif
503 
504 if sc eq 'y' then begin
505     mag=2.
506     dispsize=mag*gsize
507     gfmmagsc=rebin(gfmmag,dispsize,dispsize,/sample)
508     set_plot,'x'
509     loadct,0
510     window,!d.window+1,xsize=dispsize,ysize=dispsize
511     tvscl,gfmmagsc
512     endif
513 
514 ; SAVE FLUX MAP INTO AN EPS FILE
515 if ps eq 'y' then begin
516     set_plot,'ps'
517     device,filename=fn+'.eps',xsize=fms,ysize=fms,xoffset=(21.-fms)/2.,$
518     yoffset=(29.7-fms)/2.,/encapsulated,bits_per_pixel=8
519     tvscl,gfmmag,xsize=fms,ysize=fms,/centimeters
520     device,/close
521     set_plot,'x'
522     endif
523 ; Note that the image is centered on an a4 page
524 
525 ; SAVE IMPORTANT DATA INTO AN IDL STRUCTURE
526 if st eq 'y' then begin
527     ; Simulation Structure
528     result=execute('simstr={dadsfile:dadsfile,date:date,instrume:instrume'+$
529     ',camerastr:camerastr,targname:targname,ra_targ:ra_targ,dec_targ:dec_targ'+$
530     ',pa_v3:pa_v3,filtnam1:filtnam1,filtnam2:filtnam2,filtnu:filtnu'+$
531     ',atodgain:atodgain,gain:gain,photflam:photflam,photzpt:photzpt'+$
532     ',calpar:calpar,vi:vi,dm:dm,hpix:hpix,exptime:exptime,hstdnc:hstdnc'+$
533     ; ',hstelc:hstelc,hstmagins:hstmagins'+$
534     ',muback:muback,parhstdn:parhstdn,hstmagsta:hstmagsta,hstmag:hstmag'+$
535     ',rdim:rdim,hstmagcen:hstmagcen'+$
536     ',ss:ss,ssstr:ssstr,nobs:nobs,sd:sd,posa:posa,cx:cx,cy:cy,gaiaet:gaiaet'+$
537     ',rn:rn,spfx:spfx,spfy:spfy,obs:obs'+$
538     ',ssfx:ssfx,ssfy:ssfy,gstep:gstep,gsize:gsize,nfm:nfm,gfm:gfm,parn:parn'+$
539     ',gfmmag:gfmmag,gfmmagsta:gfmmagsta'+$
540     ',px:px,py:py,psx:psx,psy:psy,sx:sx,sy:sy,xi:xi,yi:yi}')
541     endif
542 
543 ; SAVE IMPORTANT DATA INTO A DAT FILE
544 if stfile eq 'y' then begin
545     save,filename=simfn+'.dat',/verbose,simstr
546     endif
547 
548 end