Reduction notes:
Here are some basic steps in the reduction of the imaging data. In what follows I am using the format of IDL commands. I admit that it’s a bit of a stream of consciousness, but is at least somewhat complete.
1) Start IDL
sandusky> idl
(starts command line version of IDL)
sandusky> idlde
(starts the gui version)
This gives you the IDL prompt
IDL>
This means you are no longer at the command level, though (if you use a UNIX based OS) you can execute commands by typing $ and the command. For example.
IDL>$ls
Generates a directory listing just as if you were at the command level of the OS.
2) Build the Flat field and Bias.
You will find that there are several bias and flat field frames taken in any given night. The idea here is to use them collectively to make a generalized flat field and bias. So, your first step is to open the bias (or flat field) files.
IDL> b1=readfits(‘firstbias.fits’,h1)
(where b1 is the bias image and h1 is the header)
Repeat as needed
IDL> bn=readfits(‘nthbias.fits’,hn)
For the torus imager, all of the images will be floating point arrays that are 512 by 512. You can find their dimensions with this command.
IDL> help, bias
(or any other vector)
This will return
bias=FLTARR(512,512)
Once you have done this you can try to merge them into a single bias. There are two ways to do this. Here they are.
A) Addition and filtering: This approach involves averaging them all together and then filtering the image to remove defects. Use THIS one for the flat fields. It’s optional for the bias.
IDL> btot=(b1+b2+b3+…+bn)/n
The filtering is achieved by using a program called sigma_filter.
IDL> bfilt=sigma_filter(btot,4,n_sigma=3,/iterate,/all_pixels)
(where the 4 is the width of a region in the image over which you get an median value, n_sigma is the statistical deviation you will allow inside that box, such that anything further off than that value is replaced by the median, /iterate means the program should run recursively until it finds no further outliers, and /all is used to make sure that the median filter doesn’t skip the pixels on the edges of the image. Check out the web site idlastro.gsfc.nasa.gov/contents.html for more information.)
IMPORTANT NOTE! You may find that when you run sigma_filter that it tells you the procedure isn’t there. Don’t worry. It just means that the path isn’t set right for automatic compiling (readfits is the same way, so you may already have run into this problem). In this case you just need to compile it yourself, which is just SOOOOO hard.
IDL>.run sigma_filter
Got that? If this doesn’t work, then check your spelling.
B) Voting: This is a means of building a bias that goes pixel by pixel to come up with a median image that gets rid of any spurious counts. It actually involves programming, so I’ll include an inefficient way of doing this that showcases some common idl programming commands.
IDL> bfilt=fltarr(512,512)
(here you must create the final image)
IDL> .run
-for n=0, 511 do begin
-for nn=0,511 do begin
-tmp=[b1(n,nn),b2(n,nn),b3(n,nn),….,bn(n,nn)]
(note the square brackets. They indicate that you are creating array values. The curved brackets indicate array elements in b1 and b2 that you are putting into the tmp array)
-bfilt(n,nn)=median(tmp)
(median is an idl command that returns the median value of the array)
-endfor
-endfor
-end
IDL>
(and there you go….)
Finally you should save the image.
IDL>writefits,bfilt,’filteredbias.fits’,h
(The header is optional here)
Now, I say to use the averaging method for the flat field, because the exposure lengths and sky brightness change from frame to frame in the flat fields. This is unlike the bias, which is always at very close to the same level. If the flats were all taken to the same level of brightness you could use voting there too, but for this project I don’t recommend it.
Just for the sake of clarity I recommend at this point that you close out the IDL session and start over for the actual data analysis.
IDL>exit
3) Open the data files
sandusky>IDL
(start the session-if needed)
IDL> bias=readfits(‘filteredbias.fits’,hb)
(bias is the averaged flat field and h is the header-if you saved one)
IDL> flat=readfits(‘filteredflat.fits’,hf)
(flat is the averaged flat field, hf the header-if you saved one)
IDL> onb=readfits(‘onbandfile.fits’,hon)
(onb is the onband, SII, image, hon the header)
IDL> offb=readfits(‘offbandfile.fits’,hoff)
(offb is the offband, scattered light, image, hoff the header)
So, how do you find which is the one is on band vs. off band? Well, you should look at the header. Some of the files have VERY minimal headers. Their names are typically written in UPPER CASE, for example TORUS00.FIT. The ones in lower case have been minimally processed by Ron (torus00.fits) and will have a FAR more expansive header. To see it you just ask to print out the vector.
IDL>print,h
This gives you something that looks like...
SIMPLE = T / Written by IDL: 24-Jan-2001 16:59:46.00
BITPIX = 16
NAXIS = 2
NAXIS1 = 512
NAXIS2 = 512
BUNIT = 'ADU ' /Pixel value units
ORIGIN = 'Wisconsin Space Physics, NASA Goddard Planetary Science' /
OBSERVER= 'Woodward' /
TELESCOP= 'McMath-Pierce West Auxiliary' /
INSTRUME= 'Torus CCD Imager' /
COMMENT Sky clear but hazy
FBASE = 'torus00 ' /Base filename image in .fit, header in .hdr
IMAGTYPE= 'OBJECT ' /Image type (BIAS, COMP, FLAT, OBJECT)
OBJECT = 'west torus' /
FILTER = '6735/9 ' /
HOURANGL= '3:30 ' /
AIRMASS = 1.51000 /
COMMENT Moons finally out of way
TIMEBASE= 1000 /Units of EXPTIME in ms
EXPTIME = 600 /Exposure time in in units of TIMEBASE
DATE = '11/01/2001' /Observation date (dd/mm/yyyy) (UTC)
TIME = '07:31:19' /Observation time (UTC)
CCDTSET = -100.000 /CCD temperature setpoint (deg C)
CCDTEMP = -90.7200 /CCD temperature (deg C)
STARTIME= '07:31:19' /Start time of exposure (UT)
FIRSTCOL= 0 /First column of CCD subimage
FIRSTROW= 0 /First row of CCD subimage
PARBIN = 1 /Parallel binning (columns)
SERBIN = 1 /Serial binning (rows)
ENDTIME = '07:41:33' /End time of exposure (UT)
SAVETIME= '07:41:33' /UT when file was saved
PEXPTIME= 200 /pvc dialog box exposure time (prefer EXPTIME)
CAMNAME = 'Camera A (16)' /Camera name
CAMSN = 258 /Camera head serial number
CCDNAME = 'S200/AIA SI502B' /CCD name
PVCAMVER= 544 /CCD Controller software version
HISTORY 31 keywords added from .HDR file
In the various fields you will get everything, including time of observation, filter (on vs. off) used, airmass, exposure time, date, etc.. In the above case the filter listed is ‘6735/9’, which is the SII on band filter.
4) ‘Flatten’ the flat field.
The flat field is supposed to measure out variations in the sensitivity of the optics with location across the CCD. It usually looks brightest in the center and faintest in the corners. Thus, you need to use the flat field to correct the data image. There are 3 steps to creating a true ‘flat’ correction image.
IDL>flat1=flat-bias
(This removes the bias level from the flat field, leaving only the data.)
IDL>avg=total(flat1)/n_elements(flat1)
(ok, now we’ve used two new commands here. Total(x) returns the added value of every element in array ‘x’, while n_elements(x) counts how many elements there are in ‘x’. Note that x can be any array or any PORTION of an array. For example IDL>total(flat(250:450,125:375)) returns the total flux in pixels covering the xrange from pixel 250 to 450 and the yrange from pixel 125 to 375. This is a general property of IDL that works for all commands, like the median() command used above.)
IDL>flat2=flat1/avg
(this will create a true flat field that can be used to correct the data image).
5) Process each image:
There’s actually not much to this. We probably won’t need to work with dark frames, which deletes a step. Here’s what you do.
IDL>onb1=onb-bias
(again, we just want to get rid of the bias level)
IDL>onb2=onb1/flat2
(this step takes out the variations in the average throughput.)
The same steps apply to the offb image.
Next you will probably want to get rid of the cosmic rays. For now the only program I recommend is sigma filter. Be careful with this one though, because it can have a BIG effect on the scattered light distribution if you push it too hard. I recommend going with n_sigma=5 or higher to catch the cosmic rays.
IDL>onimg=sigma_filter(onb2,3,n_sigma=5,/iterate,/all)
You also will want to keep the box width small to about 3 or 4. There will be some cosmic rays you just can’t catch this way, but you’ll get most and the filtered image will not be as affected by the slope of the scattered light.
6) Image display:
Now, you may want to take a look at the images at this point (or earlier). Here are some useful commands to help you with this.
IDL>window,xsize=x,ysize=y
(this creates a window of dimensions x,y….you will need to define x and y as variables before this OR just type in the numbers)
another option might be this
IDL>win,x,y,n
(this is a macro I wrote that you may have access to that’s just a shortcut of the above. The n value assigns a value to the window that you will see in the header strip of the display and looks like ‘IDL2’)
To switch between windows is simple.
IDL>wset,n
This command is very useful if you want to have several images up at once, for example the on band, off band, and the subtraction.
The most basic command for displaying an image is as follows.
IDL>tvscl,image
A note of caution, this will display the entire image, whether the image window matches it or not! So, if you have a 2048 x 2048 image and a 512 x 512 window, then you will only have the lower left corner (0:511,0:511) displayed.
Stretching the display. You will find that it helps you see what’s going on if you use a nifty tool in IDL called xloadct.
IDL>xloadct
xloadct allows you to change the color table being used as well as the gamma, minimum value, and maximum value. The options and functions buttons let you do things like reverse the table and non-linearly affect the color-intensity relationship. DO NOT use the ‘replace” table button under options.
Your display may still not be so good even with xloadct. That may be unavoidable if there are still cosmic rays that you haven’t gotten rid of, however one way to address this is to use a log stretch of the image.
IDL>tvscl,alog10(image)
Keep in mind that you need to be sure that there are no pixel values that are ≤ 0 if you do this.
7) Difference the images:
What you are trying to do here is to subtract the scattered light from Jupiter. To do that you should subtract the off band image from the on band image. Why? Because Jupiter’s contribution is pure continuum, not line emission, so the off band image distribution will be the same as the on band image. The caveat is that there won’t be the same amount of scattered light in the on and off band images. The first difference is in the widths of the filters. The on band is 9Å wide, while the off band is 30Å. A good starting point is to just subtract the off band filter after dividing by the ratio of their widths (3.3).
IDL>diff=onimg-offimg/3.3
Now, this will most likely not work because you haven’t accounted for the difference in filter transmission between the two, the shape of the transmission curve (it’s not a knife-edge that drops to zero on either side of the width!), and the solar spectrum, which has Fraunhofer lines. THAT tweaking is what YOU need to do.
7) Save the result:
Once you are done you should save your work! There are two ways to do this.
A)FITS: This is what you’d expect. Just type.
IDL>writefits,data,’mydata.fits’,header
B)IDLsave: IDL has a function that allows you to save a collection of arrays for recall in a future IDL session. I like using this, but it’s a question of taste. In this case you can save more than just your data.
IDL>save,diff,onimg,offimg,hdiff,honb,hoffb,factor,filename=’mydata.sav’
And, there you go! Only one command left!
IDL>exit