;; (field-rotated) images generation for SPHERE-IRDIS ;; 1st part & 2nd part identical to pfb_red_ird_snr and pfb_ird_red_contrast ;; ;; data are saved either in IDL format or in FITS format ; ;; - originally written ;; [marcel.carbillet@unice.fr & isabelle.smith@unice.fr - april 2007] ;; * after Soft.Pack.SPHERE version 2.1: ;; - sky background computation improved. ;; [arthur.vigan@oamp.fr - april 2007] ;; - time-exposure and number-of-exposures computations corr'd. ;; [isabelle.smith@unice.fr & marcel.carbillet@unice.fr - may 2007] ;; * for Test Case DRH data generation: ;; - new input: one background image per exposure ;; [isabelle.smith@unice.fr & marcel.carbillet@unice.fr - may/july 2008] ; common bruits, seed_ph ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; 1st part: input data & parameters ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; CPU time init. t0=systime(/s) ;; set directories for diffractive simulation data & photometric data dir_data = !CAOS_ENV.WORK+!CAOS_ENV.DELIM+'Projects'+!CAOS_ENV.DELIM $ + 'TestCaseDRH'+!CAOS_ENV.DELIM+'TestCaseOutputs'+!CAOS_ENV.DELIM filename = 'dataDRH' dir_photometry = !CAOS_ENV.MODULES+!CAOS_ENV.DELIM+(ird_info()).pack_name $ + !CAOS_ENV.DELIM+'pack_data'+!CAOS_ENV.DELIM+'data-reduction' $ + !CAOS_ENV.DELIM ; OUTPUT: dir_output = !CAOS_ENV.WORK+!CAOS_ENV.DELIM+'Projects'+!CAOS_ENV.DELIM $ + 'TestCaseDRH'+!CAOS_ENV.DELIM+'TestCaseFinalData'+!CAOS_ENV.DELIM ;; read diffractive simulation data & set parameters coronotype='ALC_idealcase_4.0lovD' ;;; 'lyotapo4.5' '4qperfect' 'lyotapo4.0' etc. filtre0='H2' filtre1='H3' filter='IRDIS_'+filtre0+filtre1 ;filename=coronotype ;+'_'+filter ;restore, FI=dir_data+filename+'.sav' restore, FI=dir_data+filename+'_0'+'.sav' ; structure "images" restored, where: ; corotag = on-axis PSFs(dim, dim, filters, obj/ref) ; airytag = off-axis PSF(dim, dim, sep, filters) ; resolution= pixel size [arcsec/px] ; lambda = PSFs wavelengths [m] ; width = PSFs spectral widths [m] ; instrument= SPHERE instrument ('IRDIS' here) ; plan_sep = off-axis separation ["] ; filename = name of the file in which this structure was saved coropsf = images.corotag ; star PSFs airypsf = images.airytag ; planets' PSF lambda = images.lambda*1E6 ; wavelenghts (m -> microns) width = images.width*1E6 ; bandwidths (m -> microns) ;; ONLY ONE PSF AND NOT 5: AT 2" ;status = tag_exist(images, 'plan_sep') ;if status then pos_arcsec = images.plan_sep $ ; ; FOR DATA SIM'D FROM VERSION 2.1 OF SOFT.PACK.SPHERE ;else pos_arcsec = [0.1,0.2,0.5,1.0,2.0] ; ; FOR DATA SIM'D WITH VERSION 2.0 OF SOFT.PACK.SPHERE pos_arcsec = [0.1,0.2,0.5,1.0,2.0] nfiltre = n_elements(lambda) ; nb of filters taille = float((size(images.corotag))[1]) ; nb of x- & y-pixels (should be 1024.) ;; ONLY ONE PSF AND NOT 5: AT 2" ;nplanet = (size(images.airytag))[3] ; nb of planets ;; REDEFINED LATER FROM POS_ARCSEC nplanet = (size(pos_arcsec))[1] ; nb of planets DDout = 160. ; pupil size [px] resol = taille/DDout ; resolution element (lambda/D) size [px] temp = [400, 700, 1200, 2500] ; possible planet temperatures [K] ntemp = n_elements(temp) ; nb of planet temperatures temperature= 2 ; planet temperature cons'd (0 will be for the star) ; [1: 400K, 2: 700K, 3: 1200K, 4: 2500K] Sp = 'M0' ; star spectral type: 'M0' or 'G0' distance = 10. ; star distance [pc] dec = -45 ; star declination [deg.] hahr_init = -2 ; hour angle at begining of observation [hr] dn = 0.5 ; max. nb of pixel between 2 inst. images for the 1"0 planet ffoffset = 1E-3 ; flat-field noise ronamp = 10. ; RON amplitude [e-] seed_ph = 1L ; seed for photon noise seed_ron = 1L ; seed for read-out noise seed_ff = 10000L ; seed for flat-field noise k = 1.0 ; subtraction factor to tweak intensity ratio dynamic = 1E5 ; detector dynamic ;; dynamic not used anymore: images should be saturated dyn_factor = 0.5 ; max. allowed dynamic fraction ;; REDEFINED FROM CAOS-SIMULATION OUTPUT FILES ;tpose = 18. ; single-image exposure time [s] ;Tint = 3600. ; global exposure time [s] ;; ADDED nbposes = 72 ; number of output files (number of exposures = 2*nbposes) tpose = 100 Tint = 2*nbposes*tpose Dtel = 8. ; telescope diameter [m] area = !PI*Dtel^2./4.*(1-.14^2.) ; telescope surface [m^2] tr_tel = 0.88 ; telescope transmission tr_CP = 0.60 ; common path transmission tr_IRDIS = 0.38 ; IRDIS transmission output = 'FITS' ; output data format ('FITS' or 'IDL') outind = 1B ; save also whole set of short-exposure data ? ;distancestr = strmid(strcompress(distance,/r),0,2)+'pc' ;; NEXT PART MOVED IN THE LOOP OVER THE DIFFERENT TIME EXPOSURES: ;; reformatting in 4D arrays if needed (compatibility with version 1.0 code) if (size(coropsf))[0] eq 3. then begin coropsftmp=fltarr(taille,taille,2,nfiltre) coropsftmp[*,*,0,0]=coropsf[*,*,0] coropsftmp[*,*,0,1]=coropsf[*,*,1] coropsftmp[*,*,1,0]=coropsf[*,*,2] coropsftmp[*,*,1,1]=coropsf[*,*,3] coropsf=coropsftmp endif if (size(coropsf))[0] eq 4. then begin coropsftmp=fltarr(taille,taille,2,nfiltre) for l=0,nfiltre-1 do begin coropsftmp[*,*,0,l]=coropsf[*,*,l,0] coropsftmp[*,*,1,l]=coropsf[*,*,l,1] endfor coropsf=coropsftmp endif ;; apodizer case (for further apodizer transmission computation) dirapopup = !CAOS_ENV.MODULES+!CAOS_ENV.DELIM+(ird_info()).pack_name+!CAOS_ENV.DELIM $ + 'pack_data'+!CAOS_ENV.DELIM apofile=dirapopup+'apod-4.0lovD.fits' pupfile=dirapopup+'pupilUT_PALC.fits' ;; define transmission of apodizer * Lyot stop tr_stop=1. if (coronotype eq '4qachro') or (coronotype eq '4qmono') $ or (coronotype eq '4qperfect') or (coronotype eq '4qachrow') then tr_stop=.73 if (coronotype eq 'lyot3') then tr_stop=.65 if (coronotype eq 'lyot4') then tr_stop=.74 if (coronotype eq 'lyot5') then tr_stop=.8 if tr_stop eq 1. then begin ; all other cases = apodized Lyot cases apo=readfits(apofile) & pup=readfits(pupfile) tr_stop=total(apo^2*pup)/total(pup) print, "transmission (apodizer * Lyot stop)=", strtrim(tr_stop,2) endif ;; Transmission atmosphere => detector atmos=fltarr(2,4413) openr,1,dir_photometry+'atmosSW.txt' readf,1,atmos close,1 tr_atmos = fltarr(nfiltre) for l=0,nfiltre-1 do tr_atmos[l] = atmos(1,max(where(atmos[0,*] le lambda[l]))) qe=fltarr(2,190) openr,1,dir_photometry+'QE_hawaii2RG_10nm.txt' readf,1,qe close,1 tr_qe = fltarr(nfiltre) for l=0,nfiltre-1 do tr_qe[l] = qe(1,max(where(qe[0,*] le lambda[l]))) tr=tr_tel*tr_CP*tr_IRDIS*tr_stop*tr_atmos*0.95*tr_qe print,'total transmission',tr ;; resampling IRDIS (lambda/2d at 0.95mic) lambda_sh=0.95 sampin=lambda/lambda[0]*resol sampout=lambda/lambda_sh*2. newtaille=round(sampout[0]/sampin[0]*taille) newresol=resol*newtaille/taille ;; MOVED IN THE LOOP : coro_irdis=fltarr(newtaille,newtaille,2,nfiltre) ;; ONLY ONE PSF AND NOT 5: AT 2" airy_irdis=fltarr(newtaille,newtaille,nfiltre) for i=0,nfiltre-1 do begin coro_irdis[*,*,*,i]=pfb_rescale(double(coropsf[*,*,*,i]),newtaille) airy_irdis[*,*,i]=pfb_rescale(double(airypsf[*,*,0,i]),newtaille) endfor pixel = lambda_sh*1E-6/Dtel/2.*!RADEG*3600. ; IRDIS CCD pixel size ["] ; (twive Shannon @ 0.95 microns) print,'WARNING: pixel size is from detection 0.017"/pix but from Shannon: ',pixel print, 'Nulling rejection with corono (total/max): ', $ total(airy_irdis[*,*,0])/total(coro_irdis[*,*,0,0]), $ max(airy_irdis[*,*,0])/max(coro_irdis[*,*,0,0]) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; 2nd part: photometry ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; flux in photons/s/m2 for R=30 , for 10pc filtrelist=['Y1', 'Y2', 'Y3', 'J1', 'J2', 'J3', 'J4', 'H1', 'H2', 'H3', 'H4','K1','K2'] fluxtab=readfits(dir_photometry+'old/photometrie_DBI_'+Sp+'_10pc_4planets.fits')*10.^2./distance^2. flux=fltarr(ntemp+1,2) flux[*,0]=fluxtab[*,(where(filtrelist eq filtre0))] flux[*,1]=fluxtab[*,(where(filtrelist eq filtre1))] print,'Flux in '+filtre0, ': ', flux[*,0] print,'Flux in '+filtre1, ': ', flux[*,1] ;; NUMBER OF EXPOSURES IS NOW DETERMINED BY THE NUMBER OF INPUT FILES (see input parameters) ;powermaxcoro_f=fltarr(nfiltre) ;for f=0,nfiltre-1 do $ ; powermaxcoro_f[f]=max(coro_irdis[*,*,0,f]/total(airy_irdis[*,*,f])*(flux[0,f]*area*tr[f])) ;powermaxcoro=max(powermaxcoro_f) ;tposemax = dyn_factor*dynamic/powermaxcoro ;print, 'tpose max = ', tposemax ;Dmpla=2.5*alog10(flux[0,0]/flux) ; differences of magnitude ;; time exposure stuff ;if (tpose gt tposemax) then begin ; tpose = tposemax ; print, ' ' ; print,'!!! WARNING: because of saturation, new exposure time becomes ', tpose, 's !!!' ; print, ' ' ;endif ;nbposes = round(Tint/tpose) ;dummy = Tint ;Tint = nbposes*tpose ;if Tint ne dummy then begin ; print, ' ' ; print, '!!! WARNING: total integration time is', Tint, 's !!!' ; print, ' ' ;endif ;; sky background arrays computation n_sky = fltarr(2) mskyarcsec = fltarr(2) skytab = fltarr(newtaille,newtaille,nfiltre) for l=0,nfiltre-1 do begin ; sky background values from ; http://www.eso.org/gen-fac/pubs/astclim/paranal/skybackground/ ;; LAMBDA EN MICRONS if (lambda[l] ge 2.) then mskyarcsec[l]=13. if (lambda[l] ge 1.45) and (lambda[l] lt 2.) then mskyarcsec[l]=14.4 if (lambda[l] le 1.45) then mskyarcsec[l]=16.5 ; if (lambda[l] ge 1.9125d-6) then mskyarcsec[l]=13. ; if (lambda[l] ge 1.4375d-6) and (lambda[l] lt 1.9125d-6) then mskyarcsec[l]=14.4 ; if (lambda[l] ge 1.06d-6) and (lambda[l] lt 1.4375d-6) then mskyarcsec[l]=16.5 ; if (lambda[l] le 1.06d-6) then mskyarcsec[l]=19.71 ; sky background nb of photons/px n_sky[l] = (n_phot(mskyarcsec[l], LAMBDA=lambda[l]*1E-6, WIDTH=width[l]*1E-6))[0] $ * area*tr[l]*tpose*pixel^2 skytab[*,*,l]=replicate(n_sky[l],newtaille,newtaille) endfor ;; MOVED AS BEFORE : ;; post-coronagraph on-axis PSFs objectsky = fltarr(newtaille,newtaille,2,nfiltre) for f=0,nfiltre-1 do begin coro_irdis[*,*,0,f] = coro_irdis[*,*,0,f]/ $ ; object total(airy_irdis[*,*,f])*(flux[0,f]*area*tr[f]*tpose) coro_irdis[*,*,1,f] = coro_irdis[*,*,1,f]/ $ ; reference total(airy_irdis[*,*,f])*(flux[0,f]*area*tr[f]*tpose) objectsky[*,*,0,f] = coro_irdis[*,*,0,f] + skytab[*,*,0] ; object + sky objectsky[*,*,1,f] = coro_irdis[*,*,1,f] + skytab[*,*,1] ; reference + sky endfor ;; post-coronagraph 2.0-arcsec off-axis PSFs (equivalent to pre-coronagraph on-axis PSFs) ;psf = fltarr(newtaille,newtaille,2,nfiltre) ;for f=0,nfiltre-1 do begin ; ; object + sky ; psf[*,*,0,f] = airy_irdis[*,*,0]/total(airy_irdis[*,*,0])*(flux[0,f]*area*tr[f]*Tint) $ ; + skytab[*,*,0] ; ; reference + sky ; psf[*,*,1,f] = airy_irdis[*,*,1]/total(airy_irdis[*,*,0])*(flux[0,f]*area*tr[f]*Tint) $ ; + skytab[*,*,1] ;endfor ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; 3rd part: exosystem image computation ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ;; ONE DIMENSION SUPPRESSED: NO MORE REFERENCE IMAGES totimg = fltarr(newtaille,newtaille,nfiltre) hahr = hahr_init ; init. value of hour angle [hr] angle_rot = 0. ; init. value of planet rotation angle [rad.] dis_planet = pos_arcsec/pixel ; planets separations [px] for ntime=0,nbposes-1 do begin ;; 2 lines ADDED numero1 = string(2*ntime+1,format='(I4.4)') numero2 = string(2*ntime+2,format='(I4.4)') print, "CAOS output file #", strtrim(ntime+1,2), " / ", strtrim(nbposes,2), "..." ;; MOVED FROM AN UPPER PART OF THE PROGRAM: ; next lines commented when a stationary background is wanted for comparison restore, FI=dir_data+filename+'_'+strtrim(ntime,2)+'.sav' coropsf = images.corotag ; star PSFs ; ONLY ONE PSF AND NOT 5: AT 2" ;;airypsf = images.airytag ; planets' PSF ;;window, 1, XS=newtaille, YS=newtaille ;;tvscl, rebin(abs(coropsf(*,*,0,0))^.2,512,512) ; MOVED FROM AN UPPER PART OF THE PROGRAM: ; reformatting in 4D arrays if needed (compatibility with version 1.0 code) coropsftmp=fltarr(taille,taille,2,nfiltre) for l=0,nfiltre-1 do begin coropsftmp[*,*,0,l]=coropsf[*,*,l,0] coropsftmp[*,*,1,l]=coropsf[*,*,l,1] endfor coropsf=coropsftmp ;; MOVED FROM AN UPPER PART OF THE PROGRAM: coro_irdis=fltarr(newtaille,newtaille,2,nfiltre) ;;airy_irdis=fltarr(newtaille,newtaille,nfiltre) for i=0,nfiltre-1 do begin coro_irdis[*,*,*,i]=pfb_rescale(double(coropsf[*,*,*,i]),newtaille) ;; airy_irdis[*,*,i]=pfb_rescale(double(airypsf[*,*,0,i]),newtaille) endfor ;; MOVED FROM AN UPPER PART OF THE PROGRAM: ;; post-coronagraph on-axis PSFs objectsky = fltarr(newtaille,newtaille,2,nfiltre) for f=0,nfiltre-1 do begin coro_irdis[*,*,0,f] = coro_irdis[*,*,0,f]/ $ ; object total(airy_irdis[*,*,f])*(flux[0,f]*area*tr[f]*tpose) coro_irdis[*,*,1,f] = coro_irdis[*,*,1,f]/ $ ; reference total(airy_irdis[*,*,f])*(flux[0,f]*area*tr[f]*tpose) objectsky[*,*,0,f] = coro_irdis[*,*,0,f] + skytab[*,*,0] ; object + sky objectsky[*,*,1,f] = coro_irdis[*,*,1,f] + skytab[*,*,1] ; reference + sky endfor ;; 1 line ADDED, the rest of the loop is just duplicated to use the "star of interest" image as the first time step image, and the "reference" image as the second time step image print, "Exposure #", strtrim(2*ntime+1,2), " / ", strtrim(2*nbposes,2), "..." ;; ADDED FOR INFORMATION: print, "Total number of photons reaching the CCD for f=0:", total(coro_irdis[*,*,0,0]) planets = fltarr(newtaille,newtaille,nfiltre) ;; ONE DIMENSION SUPPRESSED: NO MORE REFERENCE IMAGES images = fltarr(newtaille,newtaille,nfiltre) vit_rot_init=dparang(DEC=dec,HAHR=hahr) ; [deg./s] dt=dn/(abs(vit_rot_init)*dis_planet[3]) ; needed time [s] for the rotation ; of "dn" pixels at radius #3 (1"0) if dt le tpose then begin n_petites_rot=floor(tpose/dt) endif else begin n_petites_rot=1 dt = tpose endelse dt_last = tpose - n_petites_rot*dt for k=0,n_petites_rot-1 do begin ; planets rotation hahr += dt/3600. ; current hour angle vit_rot=dparang(DEC=dec,HAHR=hahr) ; rotation velocity [deg./s] angle_rot += vit_rot*dt*!DTOR ; rotation angle [rad.] for i=0,nplanet-1 do begin for f=0,nfiltre-1 do begin ; planets dummy=fftshift(airy_irdis[*,*,f], $ dis_planet[i]*cos(!PI/4+angle_rot), $ dis_planet[i]*sin(!PI/4+angle_rot) ) planets[*,*,f]+=dt*dummy/total(dummy)*total(airy_irdis[*,*,f]) endfor endfor endfor if dt_last ne 0. then begin dt=dt_last hahr += dt/3600. ; current hour angle vit_rot=dparang(DEC=dec,HAHR=hahr) ; rotation velocity [deg./s] angle_rot += vit_rot*dt*!DTOR ; rotation angle [rad.] for i=0,nplanet-1 do begin for f=0,nfiltre-1 do begin ; planets dummy=fftshift(airy_irdis[*,*,f], $ dis_planet[i]*cos(!PI/4+angle_rot), $ dis_planet[i]*sin(!PI/4+angle_rot) ) planets[*,*,f]+=dt*dummy/total(dummy)*total(airy_irdis[*,*,f]) endfor endfor endif planets/=tpose for f=0,nfiltre-1 do begin ;; planets photometry planets[*,*,f] = planets[*,*,f] / $ total(airypsf[*,*,0,f])*(flux[temperature,f]*area*tr[f]*tpose) ;; photon noise ; tmp_seed=seed_ph images[*,*,f] = pfb_photon_seed(objectsky[*,*,0,f]+planets[*,*,f]) ;; NEXT LINE COMMENTED: PB WITH THE SEEDS? ; images[*,*,1,f] = pfb_photon_seed(objectsky[*,*,1,f]) ; seed_ph=tmp_seed ; non-coronagraphic case... ; images[*,*,f] = pfb_photon_seed(psf[*,*,0,f]+planets[*,*,f]) ; images[*,*,1,f] = pfb_photon_seed(psf[*,*,1,f]) ;; flat-field noise ff_s=randomn(seed_ff,newtaille,newtaille)*ffoffset+1. idx = where(ff_s le 0., dummy) & if (dummy gt 0) then ff_s[idx]=1. images[*,*,f]*=ff_s ; ff_r=randomn(seed_ff,newtaille,newtaille)*ffoffset+1. ; idx = where(ff_r le 0., dummy) & if (dummy gt 0) then ff_r[idx]=1. ; images[*,*,1,f]*=ff_r ;; RON images[*,*,f]+=randomn(seed_ron,newtaille,newtaille)*ronamp ; images[*,*,1,f]+=randomn(seed_ron,newtaille,newtaille)*ronamp endfor ;; put possible negative pixels to zero ;; (can theoretically occur because of gaussian RON) idx = where(images lt 0., dummy) if (dummy gt 0) then images[idx]=0. ; window, 2, XS=newtaille, YS=newtaille ; tvscl, images(*,*,1)^.2 ;; save (also) individual short-exposure images if outind then begin ;; SAVE FILES MODIFIED if output eq 'IDL' then begin save, images, FI=dir_output+'im_filt_'+strtrim(2*ntime+1,2)+'.sav' end else if output eq 'FITS' then begin writefits,dir_output+'im_filt0_'+numero1+'.fits',images(*,*,0) writefits,dir_output+'im_filt1_'+numero1+'.fits',images(*,*,1) ; writefits,dir_output+'im_filt'+'10_'+numero1+'.fits',images(*,*,1,0) ; writefits,dir_output+'im_filt'+'11_'+numero1+'.fits',images(*,*,1,1) end ; if output eq 'IDL' then $ ; save, images, FI=dir_data+!CAOS_ENV.DELIM+'rot-'+strtrim(temp[temperature-1],2) $ ; +!CAOS_ENV.DELIM+filename+'.rot'+strtrim(ntime,2)+'.sav' $ ; else if output eq 'FITS' then $ ; writefits, dir_data+!CAOS_ENV.DELIM+'rot-'+strtrim(temp[temperature-1],2) $ ; +!CAOS_ENV.DELIM+filename+'.rot'+strtrim(ntime,2)+'.fits', images ; endif endif ;; long-exposure image totimg+=images ;; NEW TIME STEP print, "Exposure #", strtrim(2*ntime+2,2), " / ", strtrim(2*nbposes,2), "..." ;; ADDED FOR INFORMATION: print, "Total number of photons reaching the CCD for f=1:", total(coro_irdis[*,*,1,0]) planets = fltarr(newtaille,newtaille,nfiltre) ;; ONE DIMENSION SUPPRESSED: NO MORE REFERENCE IMAGES images = fltarr(newtaille,newtaille,nfiltre) vit_rot_init=dparang(DEC=dec,HAHR=hahr) ; [deg./s] dt=dn/(abs(vit_rot_init)*dis_planet[3]) ; needed time [s] for the rotation ; of "dn" pixels at radius #3 (1"0) if dt le tpose then begin n_petites_rot=floor(tpose/dt) endif else begin n_petites_rot=1 dt = tpose endelse dt_last = tpose - n_petites_rot*dt for k=0,n_petites_rot-1 do begin ; planets rotation hahr += dt/3600. ; current hour angle vit_rot=dparang(DEC=dec,HAHR=hahr) ; rotation velocity [deg./s] angle_rot += vit_rot*dt*!DTOR ; rotation angle [rad.] for i=0,nplanet-1 do begin for f=0,nfiltre-1 do begin ; planets dummy=fftshift(airy_irdis[*,*,f], $ dis_planet[i]*cos(!PI/4+angle_rot), $ dis_planet[i]*sin(!PI/4+angle_rot) ) planets[*,*,f]+=dt*dummy/total(dummy)*total(airy_irdis[*,*,f]) endfor endfor endfor if dt_last ne 0. then begin dt=dt_last hahr += dt/3600. ; current hour angle vit_rot=dparang(DEC=dec,HAHR=hahr) ; rotation velocity [deg./s] angle_rot += vit_rot*dt*!DTOR ; rotation angle [rad.] for i=0,nplanet-1 do begin for f=0,nfiltre-1 do begin ; planets dummy=fftshift(airy_irdis[*,*,f], $ dis_planet[i]*cos(!PI/4+angle_rot), $ dis_planet[i]*sin(!PI/4+angle_rot) ) planets[*,*,f]+=dt*dummy/total(dummy)*total(airy_irdis[*,*,f]) endfor endfor endif planets/=tpose for f=0,nfiltre-1 do begin ;; planets photometry planets[*,*,f] = planets[*,*,f] / $ total(airypsf[*,*,0,f])*(flux[temperature,f]*area*tr[f]*tpose) ;; photon noise ; tmp_seed=seed_ph ;; OBJECT_SKY[*,*,1,f] IF NON STATIONARY BACKGROUND ;; OBJECT_SKY[*,*,0,f] IF STATIONARY BACKGROUND images[*,*,f] = pfb_photon_seed(objectsky[*,*,1,f]+planets[*,*,f]) ;; NEXT LINE COMMENTED: PB WITH THE SEEDS? ; images[*,*,1,f] = pfb_photon_seed(objectsky[*,*,1,f]) ; seed_ph=tmp_seed ; non-coronagraphic case... ; images[*,*,f] = pfb_photon_seed(psf[*,*,0,f]+planets[*,*,f]) ; images[*,*,1,f] = pfb_photon_seed(psf[*,*,1,f]) ;; flat-field noise ff_s=randomn(seed_ff,newtaille,newtaille)*ffoffset+1. idx = where(ff_s le 0., dummy) & if (dummy gt 0) then ff_s[idx]=1. images[*,*,f]*=ff_s ; ff_r=randomn(seed_ff,newtaille,newtaille)*ffoffset+1. ; idx = where(ff_r le 0., dummy) & if (dummy gt 0) then ff_r[idx]=1. ; images[*,*,1,f]*=ff_r ;; RON images[*,*,f]+=randomn(seed_ron,newtaille,newtaille)*ronamp ; images[*,*,1,f]+=randomn(seed_ron,newtaille,newtaille)*ronamp endfor ;; put possible negative pixels to zero ;; (can theoretically occur because of gaussian RON) idx = where(images lt 0., dummy) if (dummy gt 0) then images[idx]=0. window, 2, XS=newtaille, YS=newtaille tvscl, images(*,*,1)^.2 ;; save (also) individual short-exposure images if outind then begin ;; SAVE FILES MODIFIED if output eq 'IDL' then begin save, images, FI=dir_output+'im_filt_'+strtrim(2*ntime+2,2)+'.sav' end if output eq 'FITS' then begin writefits,dir_output+'im_filt0_'+numero2+'.fits',images(*,*,0) writefits,dir_output+'im_filt1_'+numero2+'.fits',images(*,*,1) ; writefits,dir_output+'im_filt'+'10_'+numero2+'.fits',images(*,*,1,0) ; writefits,dir_output+'im_filt'+'11_'+numero2+'.fits',images(*,*,1,1) endif ; if output eq 'IDL' then $ ; save, images, FI=dir_data+!CAOS_ENV.DELIM+'rot-'+strtrim(temp[temperature-1],2) $ ; +!CAOS_ENV.DELIM+filename+'.rot'+strtrim(ntime,2)+'.sav' $ ; else if output eq 'FITS' then $ ; writefits, dir_data+!CAOS_ENV.DELIM+'rot-'+strtrim(temp[temperature-1],2) $ ; +!CAOS_ENV.DELIM+filename+'.rot'+strtrim(ntime,2)+'.fits', images ; endif endif ;; long-exposure image totimg+=images endfor window, /FREE, XS=newtaille, YS=newtaille, TIT='(final image [object, filter#0])^0.1' tvscl, totimg[*,*,0]^.1 window, /FREE, XS=newtaille, YS=newtaille, TIT='(final image [object, filter#0])^0.1' tvscl, totimg[*,*,1]^.1 ;; save total long-exposure image ;; SAVE FILES MODIFIED if output eq 'IDL' then begin ; save, totimg, FI=dir_output+filename+'.totrot-'+strtrim(temp[temperature-1],2)+'K.sav' save, totimg, FI=dir_output+'im_fin.sav' end if output eq 'FITS' then begin writefits,dir_output+'im_fin'+'0'+'.fits',totimg(*,*,0) writefits,dir_output+'im_fin'+'1'+'.fits',totimg(*,*,1) ;writefits,dir_output+'im_fin'+'10'+'.fits',totimg(*,*,1,0) ;writefits,dir_output+'im_fin'+'11'+'.fits',totimg(*,*,1,1) end ;if output eq 'IDL' then $ ; save, totimg, FI=dir_data+filename+'.totrot-'+strtrim(temp[temperature-1],2)+'K.sav' $ ;else if output eq 'FITS' then $ ; writefits, dir_data+filename+'.totrot-'+strtrim(temp[temperature-1],2)+'K.fits', totimg end