# # # IRAF cl script to demonstrate the capability for specholucy to fit # point source spectra and background for a complex scene # # Author; J. R. Walsh, ST-ECF. jwalsh@eso.org # # REQUIREMENTS: # An image display tool should be active. # The package stecf.specres must be loaded. # The package artdata needs to be loaded for producing PSF images # The packages stsdas.ttools and stsdas.graphic need to be loaded (for # STSDAS tables access and plotting of tables) # set imtype="fits" # # This script produces a simulated image consisting of a background # formed from a constant value and an elliptical brightness distribution # together with three superposed point source spectra with a range of # 10 in flux level. Two of the spectra are tilted. Noise is added to # the image and statistical error and data quality files are also # formed. # The simulated image is then fitted by a background and three point # sources using the routine stecf.specres.specholucy. # print " " print " *** Start of specres.specholucy demo *** " print " [artdata, stsdas.graphics, stsdas.ttools and stecf.specres need to be loaded] " print " " # # First delete the output files from previous efforts # print " " print " Tidying up files from any previous trials " print " " # imdel *.fits ! rm demopsfs.lis ! rm demowav.tab ! rm demopos.tab ! rm test1t*.tab # # Get 128x128, 256x256 and 256x64 images and set to zero # imcopy dev$pix[1:128,1:128] testbas.fits imarith operand1="testbas.fits" op="*" operand2=0.0 result="test128.fits" imdel testbas.fits imcopy dev$pix[1:256,1:256] testbas.fits imarith operand1="testbas.fits" op="*" operand2=0.0 result="test256.fits" # # Form a noiseless Gaussian star image of FWHM 4.0 pixels # print " " print " Creating three identical stellar Gaussian images (FWHM=4pix) " print " " # mkobjects input="test128.fits" output="psf1.fits" backgro=0. objects="demopsf.dat" \ xoffset=0. yoffset=0. star="gaussian" radius=2.0 ar=1.0 pa=0.0 magzero=7.0 gain=1. \ rdnoise=0. poisson=no seed=13 comment=yes imcopy psf1.fits psf2.fits imcopy psf1.fits psf3.fits # # Display the stellar point source image # print " " print " Displaying the stellar point source " print " " # display image="psf1.fits" frame=1 z1=0.001 z2=0.02 # # Run specpsf to form the longslit spectrum for the point sources # ls -1 psf*.fits > demopsfs.lis # # # Create a table of wavelengths of PSF's # tcreate table="demowav.tab" cdfile="col1.col" datafile="demowav.dat" nskip=0 # print " " print " Forming a PSF longslit spectrum from 3 stellar images " print " Slitwidth 0.5arcsec, pixel size 0.05arcsec, centred at Y=128.0 " print " Same pixel sampling as point source " print " " # specpsf refima="test256.fits" psfwavs="demowav.tab" \ psfims="@demopsfs.lis" sliwid=0.5 pixwid=0.050 subsam=1. cenoff=0.0 \ psfpos=128.0 outpsf="rub1_01.fits" imcopy rub1_01.fits testpsf1.fits # # Display the longslit point source spectrum # print " " print " Displaying the longslit stellar point source spectrum " print " " # display image="rub1_01.fits" frame=1 z1=0.0 z2=0.10 # # Position and rotate the longslit point source spectra # print " " print " 1st PSF spectrum: Rotation 0deg; flux 5000, Centre Y=120.5 " print " " # rotate input="rub1_01.fits" output="rub1_02.fits" rotation=0.0 ncols=256 nlines=256 imarith operand1="rub1_02.fits" op="*" operand2=5000.0 result="rub1_03.fits" imshift input="rub1_03.fits" output="rub1_04.fits" xshift=0.0 yshift=-7.50 \ interp="spline3" # print " " print " 2nd PSF spectrum: Rotation 2.5deg; flux 500, Centre Y=199.33 " print " " # rotate input="rub1_01.fits" output="rub1_05.fits" rotation=2.5 ncols=256 \ nlines=256 interpo="spline3" imarith operand1="rub1_05.fits" op="*" operand2=500.0 result="rub1_06.fits" imshift input="rub1_06.fits" output="rub1_07.fits" xshift=0.0 yshift=71.33 \ interp="spline3" # print " " print " 3rd PSF spectrum: Rotation -2.5deg; flux 1000, Centre Y=88.67 " print " " # rotate input="rub1_01.fits" output="rub1_08.fits" rotation=-2.5 ncols=256 \ nlines=256 interpo="spline3" imarith operand1="rub1_08.fits" op="*" operand2=1000.0 result="rub1_09.fits" imshift input="rub1_09.fits" output="rub1_10.fits" xshift=0.0 yshift=-39.11 \ interp="spline3" # # Add the three longslit PSF spectra, a simulated sky spectrum with # a Gaussian profile along the slit (i.e. in Y) and a background # of 50 counts # imarith operand1="rub1_04.fits" op="+" operand2="rub1_07.fits" result="rub1_11.fits" imarith operand1="rub1_11.fits" op="+" operand2="rub1_10.fits" result="rub1_12.fits" # # Display the three summed longslit point source spectra # print " " print " Displaying the 3 longslit stellar point source spectra " print " " # display image="rub1_12.fits" frame=1 z1=0.0 z2=500.0 # # # Form a 2D long slit spectrum from the 1-D template skyback.fits, # scale to a maximum of 100 counts and rotate to appear like a sky spectrum # print " " print " Sum the 3 PSF longslit spectra with a spatially varying longslit " print " background spectrum (sky) " print " " mk2dspec input="test256.fits" output="rub1_13.fits" models="template2d.lis" imarith operand1="rub1_13.fits" op="*" operand2=400.0 result="rub1_14.fits" rotate input="rub1_14.fits" output="rub1_15.fits" rotation=-90.0 ncols=256 \ nlines=256 interpo="nearest" imarith operand1="rub1_12.fits" op="+" operand2="rub1_15.fits" result="rub1_16.fits" imarith operand1="rub1_16.fits" op="+" operand2=50.0 result="rub1_17.fits" # # Display the summed background # print " " print " Displaying the background image " print " " # display image="rub1_17.fits" frame=1 z1=0.0 z2=200.0 # # Add Gaussian noise to the simulated image # print " " print " Adding Gaussian noise to the model image (RON=6 e; Gain=1) " print " " mknoise input="rub1_17.fits" output="testspec1.fits" backgro=0.0 gain=1.0 \ rdnoise=6.0 poisson=yes seed=13 # # Add a bad pixel # imreplace images="testspec1.fits[105:105,178:178]" value=1000 lower=0 upper=500 radius=0 # # Display the noisy simulated image # print " " print " Displaying the noisy simulated longslit point source + background image " print " " # display image="testspec1.fits" frame=1 z1=0.0 z2=200.0 print " " print " Adding wavelength information to header for 2-D spectrum " print " " hedit image="testspec1.fits" fields="CRPIX1" value=128.0 add=yes delete=no \ verify=no show=yes update=yes hedit image="testspec1.fits" fields="CRVAL1" value=8000.0 add=yes delete=no \ verify=no show=yes update=yes hedit image="testspec1.fits" fields="CTYPE1" value="LAMBDA" add=yes delete=no \ verify=no show=yes update=yes hedit image="testspec1.fits" fields="CD1_1" value=5.0 add=yes delete=no \ verify=no show=yes update=yes # # Form the statistical error image of the input image # print " " print " Forming the statistical error image " print " " # imexpr expr="sqrt(a)" output="testspec1er.fits" a="testspec1.fits" # # Form the data quality image of the input image # print " " print " Forming the data quality image and setting some BAD pixels " print " " imcopy test256.fits testspec1dq.fits imreplace images="testspec1dq.fits[128:128,128:128]" value=100 lower=0 upper=0 radius=0 imreplace images="testspec1dq.fits[35:39,58:58]" value=100 lower=0 upper=0 radius=0 imreplace images="testspec1dq.fits[86:86,100:109]" value=100 lower=0 upper=0 radius=0 imreplace images="testspec1dq.fits[157:157,45:45]" value=100 lower=0 upper=0 radius=0 imreplace images="testspec1dq.fits[31:31,157:158]" value=100 lower=0 upper=0 radius=0 imreplace images="testspec1dq.fits[77:77,87:95]" value=100 lower=0 upper=0 radius=0 imreplace images="testspec1dq.fits[167:172,78:79]" value=100 lower=0 upper=0 radius=0 imreplace images="testspec1dq.fits[200:200,49:49]" value=100 lower=0 upper=0 radius=0 imreplace images="testspec1dq.fits[34:34,216:216]" value=100 lower=0 upper=0 radius=0 imreplace images="testspec1dq.fits[209:209,234:234]" value=100 lower=0 upper=0 radius=0 imreplace images="testspec1dq.fits[105:105,178:178]" value=100 lower=0 upper=0 radius=0 # # Form the statistical noise image for the background from the two simulated # background images for later comparison of the noise on the fitted backgrpund # image to that expected # print " " print " Forming the noise image on the background " print " " # imarith operand1="test256" op="+" operand2=50.0 result="rub1_18.fits" imarith operand1="rub1_18.fits" op="+" operand2="rub1_15.fits" result="rub1_19.fits" imcopy rub1_19.fits modelback1.fits mknoise input="modelback1.fits" output="testback1.fits" backgro=0.0 gain=1.0 \ rdnoise=6.0 poisson=yes seed=13 imarith operand1="testback1.fits" op="-" operand2="modelback1.fits" \ result="testnois1.fits" # # Create a table of positions of point sources # tcreate table="demopos.tab" cdfile="col3.col" datafile="demopos.dat" nskip=0 # # Print the contents of the table specifying the positions of the point # source spectra # print " " print " specholucy input table for 3 point source spectra: " print " " print " X (pix) Y (pix) Slope (/pix) " tprint table="demopos.tab" # # Run specholucy on the input data for the three point sources: # print " " print " PACKAGE = stecf " print " TASK = specholucy " print " " print " inpima = testspec1 Name of input flux file " print " errima = testspec1er Name of input statistical error file " print " dqima = testspec1dq Name of input data quality file " print " poitab = demopos.tab Name of table of X;Y pixel centres and slope of sp " print " " print " OUTPUT tables and images" print " Rootab = test1t Root name for output point source extracted spectr " print " poi_back= no Output fitted point sources+background (Y) or back " print " bakima = testout1 Name of output fitted background file " print " bakerr = testout1er Name of output fitted background statistical error " print " bakdq = testout1dq Name of output fitted background data quality file " print " " print " POINT SPREAD FUNCTION Spectrum and Position " print " psfima = testpsf1 Name of input Point Spread Function file " print " psfmeth = g Method for fitting PSF peak (C=Centroid; G=Gaussia " print " psfb = 1 Subsampling factor (Y-direction) for PSF image " print " interpol = spline3 Interpolation method used for shifting PSFs " print " " print " POSITIONS of Point Source Spectra" print " posmeth = c Method for fitting input point source peak (E=Exac " print " icsum = 10 Number of columns of input data to sum for cross-c " print " " print " BACKGROUND CHANNEL parameter" print " skernel = 5.0 Sigma of Gaussian for spatial resolution kernel " print " " print " ITERATIONS and Stopping Criteria" print " niter = 50 Number of iterations for point source and 1D background restoration " print " epspoi = 0.0005 Stopping criterion for fractional change in point " print " epsbac1 = 0.0005 Stopping criterion for fractional change in backgr " print " " print " ntrial = 1 No. of Monte Carlo trials for error estimation " print " (seed = 13) Seed for random number generator" print " (dqlim = 10 ) Limiting data quality value above which to ignore " print " ( verbose = yes) Report parameters of restoration during execution " print " " # specholucy inpima="testspec1.fits" errima="testspec1er.fits" dqima="testspec1dq.fits" \ poitab="demopos.tab" Rootab="test1t" poi_back=no bakima="testout1.fits" \ bakerr="testout1er.fits" bakdq="testout1dq.fits" psfima="testpsf1.fits" \ psfmeth="g" psfb=1 interpol="spline3" posmeth="c" icsum=10 skernel=5.0 \ niter=50 epspoi=5.0E-5 epsbac=5.0E-4 ntrial=1 seed=13 dqlim=10 verbose=yes # # Plot and give statistics on the three extracted point source spectra # print " " print " Point source 1 at 120.5 with 5000 counts " print " " axispar.wb=0.0 axispar.wt=5400. dvpar.append = no pltpar.pointmode = yes sgraph "test1t_1.tab WAVELENGTH FLUX" errcolu="STATISTICAL_ERROR" tstat intable="test1t_1.tab" column="FLUX" rows=1-256 # print " " print " Point source 2 at 199.33 with 500 counts " print " " axispar.wb=0.0 axispar.wt=600. dvpar.append = no pltpar.pointmode = yes sgraph "test1t_2.tab WAVELENGTH FLUX" errcolu="STATISTICAL_ERROR" tstat intable="test1t_2.tab" column="FLUX" rows=1-256 # print " " print " Point source 3 at 88.67 with 1000 counts " print " " axispar.wb=0.0 axispar.wt=1200. dvpar.append = no pltpar.pointmode = yes sgraph "test1t_3.tab WAVELENGTH FLUX" errcolu="STATISTICAL_ERROR" tstat intable="test1t_3.tab" column="FLUX" rows=1-256 # # Display the fitted background # print " " print " Displaying the fitted background image (range 0-160 counts) " print " " display image="testout1.fits" frame=1 z1=0 z2=180 # # Finally perform statistics on the noise of the restored background # image compared with that on the input background image # print " " print " Noise on the fitted background " print " " imarith operand1="testout1.fits" op="-" operand2="modelback1.fits" \ result="testout1nois.fits" imstat image="testout1nois.fits" print " " print " Compared to that on the input background image" print " " imstat image="testnois1.fits" # # End of demo # print " " print " *** End of stecf.specres.specholucy demo *** " print " " #