;   PROGRAM             SIM_STACK.PRO

; Last Updated           05 Jun 2000          Mattia Vaccari

; This program generates a set of GAIA BBP simulated observations
;  of an HST WFPC2 field and stacks them into a GAIA flux map.

; PARAMETERS
hpix=[45.5,96.6]                     ; HST pixel in mas (PC and WFCs)
spfx=2l                              ; Subpixeling factors
spfy=2l
px=37.2                              ; GAIA pixel in mas
py=111.6
psx=[6,6,6,1,4,4]                    ; GAIA sample in pixels
psy=[8,4,2,8,8,4]
sx=psx*px                            ; GAIA sample in mas
sy=psy*py
obsside=[24998.4,55353.6]
; Side of GAIA observation obtained from a WPPC2 field in mas
;  (PC and WFCs respectively)
xi=round(obsside#(1./sx))
yi=round(obsside#(1./sy))
; GAIA observation in samples for different combinations of
;  WFPC2 cameras and BBP sample sizes
gaiaet=0.8618
; GAIA observation exposure time in s
rn=[5.44289,6.74395,8.78593,10.4357,6.12805,7.83178]
; GAIA BBP readnoise in e-/sample rms, calculated following SAG_MV_04
;  assuming full CCD readout
ssfx=psx                             ; Subsampling factors
ssfy=psy*3
; The subsamples are squares of side 37.2 mas, which is also the sample size
;  in the along scan direction
gstep=hpix/[1.,2.]                   ; GAIA flux map step in mas
; gstep equals 45.5 and 48.3 mas for PC and WFCs respectively
gsize=[350l,320l]
; Side of GAIA flux map in steps for PC and WFCs respectively
; The side of GAIA flux map is thus 15925 and 15456 mas for PC and WFCs
;  respectively

print,'The following .fits.gz files are available in the ./hstdata directory:'
spawn,'cd hstdata ; ls -C1 *.fits.gz | awk ''{print $1}''',fitsfiles,$
 count=n_fits

for i=0,n_fits-1 do begin
print,i+1,strmid(fitsfiles(i),0,strlen(fitsfiles(i))-8),format='(i3,6x,a)'
endfor

print,'Enter the number of the file you want to open:'
fnu=0 & read,fnu
print,'Opening file number ',fnu
hstdnc=readfits('hstdata/'+fitsfiles(fnu-1),header,/silent)
himdim=size(hstdnc)

f=n_elements(header)

for j=0,f-2 do begin
r=header(j)
i=strpos(r,'DADSFILE=')
if i ne -1 then d=execute(strcompress(strmid(r,0,30),/remove_all)) else begin
i=strpos(r,'ROOTNAME=')
if i ne -1 then begin
d=execute(strcompress(strmid(r,0,30),/remove_all))
DADSFILE=ROOTNAME
endif
endelse
i=strpos(r,'INSTRUME=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
i=strpos(r,'DATE    =')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
i=strpos(r,'TARGNAME=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
i=strpos(r,'RA_TARG =')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
i=strpos(r,'DEC_TARG=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
i=strpos(r,'PA_V3   =')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
i=strpos(r,'FILTNAM1=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
i=strpos(r,'FILTNAM2=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
i=strpos(r,'EXPTIME =')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
i=strpos(r,'ATODGAIN=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
i=strpos(r,'PHOTFLAM=')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
i=strpos(r,'PHOTZPT =')& if i ne -1 then d=execute(strcompress(strmid(r,0,30),$
/remove_all))
endfor

print,'DADSFILE',dadsfile,format='(a,15x,a)'
print,'INSTRUMENT',instrume,format='(a,13x,a)'
print,'OBSERVATION DATE',date,format='(a,7x,a)'
print,'TARGET NAME',targname,format='(a,12x,a)'
print,'RIGHT ASCENSCION',ra_targ,format='(a,1x,a)'
print,'DECLINATION',dec_targ,format='(a,6x,a)'
print,'V3 POSITION ANGLE',pa_v3
print,'FIRST FILTER NAME',filtnam1,format='(a,6x,a)'
print,'SECOND FILTER NAME',filtnam2,format='(a,5x,a)'
print,'EXPOSURE TIME',exptime,format='(a,4x,a)'
print,'ATODGAIN',atodgain,format='(a,9x,a)'

; HST DATA NUMBER COUNTS
if himdim(0) eq 2 then begin
print,'This file contains one image.'
imagenu=1
endif

if himdim(0) eq 3 then begin
print,'This file contains ',himdim(3),' images. Which do you want to take?'$
 ,format='(a,i1,a)'
imagenu=1 & read,imagenu
endif

imagenu=imagenu-1
hsti=imagenu<1
hpix=hpix(hsti)
camerastr=['PC','WFC2','WFC3','WFC4'] & camerastr=camerastr(imagenu)
print,'Creating HST WFPC2 '+camerastr+' image.'
hstdnc=hstdnc(*,*,imagenu)
; Note that if the .fits.gz file containing the original data consists of one
;  image only, the program assumes that this is a PC image. Otherwise, it
;  assumes the file is in the format used by the HDA for the WFPC2 data files,
;  i.e. that it contains, in this order the four images obtained with the
;  PC, WFC2, WFC3 and WFC4.

; A/D GAIN
gain=[[7.12,7.12,6.90,7.10],[13.99,14.50,13.95,13.95]]
gain=gain(imagenu,round(atodgain/7.)-1)
; HST WFPC2 A/D conversion gain, from "WFPC2 Instrument Handbook",
;  Version 4.0, page 82

lss:
print,'What BBP sample size do you want to use?'
print,'1=6x8  2=6x4  3=6x2  4=1x8  5=4x8  6=4x4'
ss=1 & read,ss
if total([1.,2.,3.,4.,5.,6.] eq ss) eq 0 then begin
print,'Invalid choice!'
goto,lss
endif
ss=ss-1

; PARAMETERS' CHOICE
psx=psx(ss)
psy=psy(ss)
sx=sx(ss)
sy=sy(ss)
xi=xi(hsti,ss)
yi=yi(hsti,ss)
rn=rn(ss)
ssfx=ssfx(ss)
ssfy=ssfy(ss)
gstep=gstep(hsti)
gsize=gsize(hsti)

; HST IMAGE PHOTOMETRIC CALIBRATION
vi=1.2
; Average V-I Galaxy Color Index from Prugniel and H\'eraudeau 1998

calpar=[[21.729,-0.051,+0.009],$ ; F555W
        [22.093,+0.254,+0.012],$ ; F606W
        [20.920,-0.124,+0.028]]  ; F814W

case filtnam1 of
'F555W' : begin
filtnu=0
stamag='V'
end
'F606W' : begin
filtnu=1
stamag='V'
end
'F814W' : begin
filtnu=2
stamag='I'
end
else : begin
print,'The ',filtnam1,' filter is not supported! ',$
 'Using F555W-V transformations!'
filtnu=0
stamag='V'
end
endcase

case stamag of
'V' : muback=22.06
'I' : muback=21.46
endcase
; Sky background in mag/arcsec^2. From SAG_GG_05.
; Note that the choice of this parameter does not affect in any way
;  the processing of the data, but their representation only

gr=[1.987,2.003,2.006,1.955]
; Gain ratio from Holtzman et al. 1995b
dm=2.5*alog10(gr)
if atodgain eq 7. then dm=dm(imagenu) else dm=0.
; Since most photometric calibration observations were done with atodgain=15,
;  there is a correction dm to be applied for exposures carried out with
;  atodgain=7

; HST IMAGE IN UNCALIBRATED MAGNITUDES
parhstdn=(exptime*(hpix/1000.)^2)*10^(0.4*(calpar(0,filtnu)+$
 calpar(1,filtnu)*vi+calpar(2,filtnu)*vi^2+dm+0.1-muback))
; Sky background parameter in data numbers
;parhstdn=4.0
; An alternative sky background parameter can be used to obtain a better
;  image visibility (i.e. an higher contrast). In the case of m100_900,
;  for instance, only 132 out of 122500 pixels (0.1%) in the [350,350] central
;  region of the PC image "contained" less than 4 data numbers.
hstmag=-2.5*alog10(hstdnc>parhstdn)
rdim=[350l,160l] & rdim=rdim(hsti)
hstmagcen=hstmag((himdim(1)-rdim)/2:(himdim(1)+rdim)/2-1,$
 (himdim(2)-rdim)/2:(himdim(2)+rdim)/2-1)
; hstmagcen covers the central part of hstmag corresponding to a square area
;  of side 15925 and 15456 mas, for PC and WFC respectively, i.e. exactly the
;  portion of the HST image that is covered by GAIA flux map.

; HST IMAGE IN WFPC2 INSTRUMENTAL MAGNITUDES
; See "WFPC2 Instrument Handbook", Version 4.0, Section 8.7
;zp=-2.5*alog10(photflam)+photzpt
; Zero-point of the "FILTNAM1" magnitude scale
;hstmagins=-2.5*alog10((hstdnc>parhstdn)/(exptime*(hpix/1000.)^2))+zp
; HST image in "FILTNAM1" magnitudes

; HST IMAGE IN STANDARD MAGNITUDES
hstmagsta=-2.5*alog10((hstdnc>parhstdn)/(exptime*(hpix/1000.)^2))+$
 calpar(0,filtnu)+calpar(1,filtnu)*vi+calpar(2,filtnu)*vi^2+dm+0.1
; HST image in stamag standard magnitudes
; The transformations between WFPC2 Data Number counts to standard magnitudes
;  are carried out using Equation 9 and Table 10 from Holtzman et al. 1995b,
;  with a constant correction of 0.1 mag/arcsec^2 for infinite aperture

; HST ELECTRON COUNTS
;hstelc=gain*hstdnc
; The atodgain keyword can be equal to 7 or 15, but the true value of the
;  A/D conversion gain (given by the variable gain) is actually very
;  near to 7 and 14, respectively, and is different in WFPC2's different
;  detectors. Please note that we do not use the gain value in the
;  conversion to a magnitude scale. Actually, we only use it to give an
;  estimate of GAIA's electron counts, which we have to combine with
;  readnoise estimates (given in electrons) to properly take into
;  account the noise.
;parhstel=gain*parhstdn     ; sky background parameter in electrons
; parhstel=30.
; An alternative sky background parameter can be used to obtain a better
;  image visibility (i.e. an higher contrast). In the case of m100_900,
;  for instance, only 171 out of 122500 pixels (0.14%) in the central region
;  of the PC image "contained" less than 30 electrons.

; GENERATION OF CENTERS AND POSITION ANGLES OF OBSERVATIONS
print,'How many simulated observations do you want to generate?'
nobs=1 & read,nobs
; if nobs eq 1 then posa=0. & cx=0. & cy=0. else begin ..... endif
lsd:
print,'How do you want the position angles of the scan directions '+$
 'to be distributed?'
print,'          1=0-360    2=0-45 or 90-135    3=0-90    4=0-45'
sd=1 & read,sd
if total([1.,2.,3.,4.] eq sd) eq 0 then begin
print,'Invalid choice!'
goto,lsd
endif
; The value of sd goes in increasing order of non-randomness
;  of the scan directions
case sd of
1 : posa=randomu(seed,nobs)*(2*!pi)
2 : begin
posa=randomu(seed,nobs)*(!pi/4)
posa=posa+((-1)^round(randomu(seed,nobs)*1000000.)*(!pi/4)+(!pi/4))
end
3 : posa=randomu(seed,nobs)*(!pi/2)
4 : posa=randomu(seed,nobs)*(!pi/4)
endcase
; posa contains the position angles (measured counterclockwise)
;  of GAIA FORs with respect to the HST FOR
cx=(randomu(seed,nobs)-0.5)*sx
cy=(randomu(seed,nobs)-0.5)*sy
; cx and cy contain the coordinates of the centers of GAIA FORs
;  (i.e. of GAIA observations) in the HST FOR

; In other words, the GAIA FOR for the nth observation is obtained
;  starting from the HST FOR by a translation of a vector [cx(n),cy(n)]
;  followed by a rotation around the center of the FOR thus obtained
;  of an angle posa(n) measured counterclockwise

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

ssstr=['6x8','6x4','6x2','1x8','4x8','4x4'] & ssstr=ssstr(ss)
print,'Generating and stacking GAIA ',ssstr,' pixels/sample BBP '+$
 'simulated observations.'
print,'Please wait...'

; EXPOSURE TIME SCALING OF HST ELECTRON COUNTS
;hgelc=(gain*gaiaet/exptime)*hstdnc
; HST electron counts scaled to GAIA exposure time
; We assume that the electron counts per unit time per unit area of the two
;  instruments are equal, which is conservative for most combinations of
;  HST-GAIA filters.

; SUBPIXELING OF HST IMAGE
sphim=rebin(hstdnc,himdim(1)*spfx,himdim(2)*spfy,/sample)*$
 ((gain*gaiaet)/(exptime*spfx*spfy))
; Subpixeling of each HST pixel into 4 subpixels
; The rebin function with the sample keyword performs a nearest neighbour
;  sampling, thus preserving the total number of electrons with great accuracy

; PSF
psfname=['68','64','62','18','48','44']
psf=readfits('psf/'+'glo_10_'+psfname(ss)+'.fits.gz',/silent)
;psf=readfits('psf/'+psfname+'.fits.gz',/silent)
; (517,517) PSF obtained following SAG_LL_025, using syntfits.f and
; assuming V-I=1.2 and Field Point 10
; The PSF is sampled with a step of 1/4 pixel along both directions
; The total area mapped by the PSF is thus of
;  (517,517)*(1/4)*(px,py)=(4.8081,14.4243) arcsec^2
xtot=[21,21,129,33,33,21] & xtot=xtot(ss)
xmed=[10,10,64,16,16,10]  & xmed=xmed(ss)
dx=[24,24,4,16,16,24]     & dx=dx(ss)
xin=(findgen(xtot)-xmed)*dx+258
ytot=[33,17,17,33,17,65] & ytot=ytot(ss)
ymed=[16,8,8,16,8,32]    & ymed=ymed(ss)
dy=[16,32,32,16,32,8]    & dy=dy(ss)
yin=(findgen(ytot)-ymed)*dy+258
psf=psf(xin,*)
psf=psf(*,yin)
; PSF resampling with a step equal to sample size
case ss of
0 : psf=psf(6:14,10:22)
1 : psf=psf(6:14,5:11)
2 : psf=psf(40:88,5:11)
3 : psf=psf(10:22,10:22)
4 : psf=psf(10:22,5:11)
5 : psf=psf(6:14,20:44)
endcase
; Reduced PSF with a step equal to sample size
; In other words, we truncate the PSF to a central region of about
; (50,50) pixel = (1860,5580) arcsec containing about 99% of the
;  complete PSF (i.e. of the energy)
psf=psf/total(psf)
; Normalization

; ARRAYS' CREATION
scx=(findgen(xi)-0.5*(xi-1))*sx
scy=(findgen(yi)-0.5*(yi-1))*sy
; scx and scy contain the coordinates of samples' centers
;  in the observation FOR
spcx=( findgen(himdim(1)*spfx)-0.5*(himdim(1)*spfx-1) )*(hpix/spfx)
spcy=( findgen(himdim(2)*spfy)-0.5*(himdim(2)*spfy-1) )*(hpix/spfy)
; spcx and spcy contain the coordinates of subpixels' centers
;  in the HST FOR
ob=fltarr(xi,yi)
; ob will contain the nth observation during processing
obs=fltarr(nobs,xi,yi)
; obs will contain the nobs observations
sscx=(findgen(xi*ssfx)-0.5*(xi*ssfx-1))*(sx/ssfx)
sscy=(findgen(yi*ssfy)-0.5*(yi*ssfy-1))*(sy/ssfy)
; sscx and sscy contain the coordinates of subsamples' centers
;  in the observation FOR
fmgrid=(findgen(gsize)-0.5*(gsize-1))*gstep
; Coordinates of the centers of the flux map elements in the HST FOR
pfm=fltarr(gsize,gsize)
; pfm will contain the partial flux map
nfm=fltarr(nobs,gsize,gsize)
; nfm will contain the nobs partial flux maps
gfm=fltarr(gsize,gsize)
; gfm will contain the global flux map

for n=0,nobs-1 do begin

; COORDINATES OF SUBPIXELS' CENTERS IN THE OBSERVATION FOR
auxx   =+((spcx-cx(n))#(fltarr(himdim(2)*spfy)+1.))*cos(posa(n))$
        +((fltarr(himdim(1)*spfx)+1.)#(spcy-cy(n)))*sin(posa(n))
auxy   =-((spcx-cx(n))#(fltarr(himdim(2)*spfy)+1.))*sin(posa(n))$
        +((fltarr(himdim(1)*spfx)+1.)#(spcy-cy(n)))*cos(posa(n))

; REBINNING OF GAIA OBSERVATION
for i=0,xi-1 do begin
xin=where(abs(auxx-scx(i)) lt sx/2.,wrx)
for j=0,yi-1 do begin
if wrx ne 0 then yin=where(abs(auxy(xin)-scy(j)) lt sy/2.,wry)$
 else wry=0
if wry ne 0 then ob(i,j)=total(sphim(xin(yin))) else ob(i,j)=0.
endfor
endfor
; Rebinning of HST subpixels into GAIA samples

; CONVOLUTION WITH PSF
ob=kconvol(ob,psf)
; !!! Note the use of kconvol instead of convol !!!

; NOISE
;obmed=median(ob)
;ob=poidev(ob>obmed,seed=seed)
; The use of the > operator overestimates the photon noise, but avoids
;  possible problems connected with the use of the poidev function with
;  an argument containing mostly very small numbers
ob=poidev(ob,seed=seed)
ob=ob+randomn(seed,xi,yi)*rn
; Addition of photon noise and readnoise

; NTH OBSERVATION
obs(n,*,*)=ob

; SUBSAMPLING OF GAIA OBSERVATION
ssob=rebin(ob,xi*ssfx,yi*ssfy,/sample)/(ssfx*ssfy)
; subsampling of each sample into ssfx*ssfy subsamples

; COORDINATES OF SUBSAMPLES' CENTERS IN THE HST FOR
auxx   =+(sscx#(fltarr(yi*ssfy)+1.))*cos(posa(n))$
        -((fltarr(xi*ssfx)+1.)#sscy)*sin(posa(n))$
        +cx(n)
auxy   =+(sscx#(fltarr(yi*ssfy)+1.))*sin(posa(n))$
        +((fltarr(xi*ssfx)+1.)#sscy)*cos(posa(n))$
        +cy(n)

; REBINNING OF GAIA PARTIAL FLUX MAP
for i=0,gsize-1 do begin
xin=where(abs(auxx-fmgrid(i)) lt gstep/2.,wrx)
for j=0,gsize-1 do begin
if wrx ne 0 then yin=where(abs(auxy(xin)-fmgrid(j)) lt gstep/2.,wry)$
 else wry=0
if wry ne 0 then pfm(i,j)=total(ssob(xin(yin))) else pfm(i,j)=0.
endfor
endfor
; Rebinning of GAIA subsamples into GAIA partial flux map

; NTH PARTIAL FLUX MAP
nfm(n,*,*)=pfm

endfor

; GAIA GLOBAL FLUX MAP
gfm=total(nfm,1)

; FLUX MAP IN MAGNITUDES
parn=(nobs*gaiaet/exptime)*(gstep^2/hpix^2)*(gain*parhstdn)
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={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
endif

end