;+ ; NAME: ; MIPS70_LINCOR (1.0) ; ; PURPOSE: ; This program applies a linearity correction to a MIPS 70um image. ; It reads in a single FITS image, which is allowed to have extensions: ; the first extension is assumed to be the image, the second, if present, ; is assumed to be the uncertainty, and all other extensions are simply ; copied to the output file (which will have the same name as the input ; file, with '.fits' replaced by '.lincor.fits'). It applies a ; correction to each pixel, putting the result in a new file. Calibrated ; data are allowed in imaging mode, as the script will look for the ; 'FLUXCONV' or 'JANSCALE' keywords to determine whether the data are ; calibrated and will modify the fit coefficient, threshold, and ; uncertainty by the calibration factor (set by 'calfac'). ; This program requires the IDL astronomy library. ; ; CATEGORY: ; MIPS (Custom) ; ; CALLING SEQUENCE: ; MIPS70_LINCOR,input_file ; ; INPUTS: ; input_file : The name of the input FITS file ; ; KEYWORD PARAMETERS: ; overwrite : overwrite output file ; ; OUTPUTS: ; The output is a FITS file with pixels corrected for nonlinearity. ; ; MODIFICATION HISTORY: ; Written by: C. W. Engelbracht (2006-06-20) ; 2007-03-23: updated correction function to MIPS70 units (CWE) ; 2007-03-28: updated to handle FITS files with extensions (KDG) ; 2007-03-28: modified correction coefficients (CWE) ; 2007-06-13: modified correction coefficients, apply slope unc. (CWE) ; 2008-01-23: rename lincor_slope_unc, check for calibrated data (CWE) ; 2008-06-09: modified fit function (CWE) ; 2008-06-23: modified fit coefficient and updated uncertainty (CWE) ; 2008-10-28: change fit function again, make itime corr. optional (CWE) ; 2009-01-21: restore error plane, remove default itime corr. (CWE) ; 2009-01-27: update uncertainty coefficients (CWE) ;- pro mips70_lincor, input_file, overwrite = overwrite, $ itime_corr = itime_corr ; set coefficients for non-calibrated data coeff = [0.219, 0.182] threshold = 10^(-coeff[0] / coeff[1]) ; set uncertainty in coefficients unc_coeff = [0.015, 0.014] ; set calibration factor, in MJy/sr/MIPS70 calfac = 702.0 ; check we got at least 1 parameter if (n_params() lt 1) then begin print, 'Syntax - mips70_lincor, input_file, overwrite = overwrite,' print, ' no_itime_corr = no_itime_corr' return endif ; check input file if (not file_test (input_file)) then begin print, 'Input file (' + input_file + ') does not exist.' retall endif if (not file_test (input_file, /read)) then begin print, 'Input file (' + input_file + ') exists but I cannot read it.' retall endif ; compose output file name len = strlen (input_file) if (strmid (input_file, len-5, len) eq '.fits') then begin output_file = strmid (input_file, 0, len-5) + '.lincor.fits' endif else begin output_file = input_file + '.lincor.fits' endelse ; check output file if (file_test (output_file)) then begin if (not keyword_set (overwrite)) then begin print, 'Output file (' + output_file + ') already exists.' retall endif if (not file_test (output_file, /write)) then begin print, 'Cannot overwrite (' + output_file + ').' retall endif endif ; open files fits_open, input_file, input_fcb fits_open, output_file, output_fcb, /write ; read the primary header/image fits_read, input_fcb, image0, header0, exten_no=0 if (input_fcb.nextend eq 0) then begin image = image0 endif ; check for calibrated data fluxconv = sxpar (header0, 'FLUXCONV') if (fluxconv eq 0) then fluxconv = sxpar (header0, 'JANSCALE') ; check for calibration in 1st extension if ((fluxconv eq 0) AND (input_fcb.nextend GT 0)) then begin fits_read, input_fcb, dummy, theader0, exten_no=1, /header_only, /pdu fluxconv = sxpar (theader0, 'FLUXCONV') endif if (fluxconv gt 0) then begin print, 'Data are calibrated, correcting fit parameters' coeff[0] -= coeff[1] * alog10 (calfac) threshold *= calfac endif ; convert to current exposure time, if requested ; N.B. The asteroid data indicate you *don't* want to use this, so it doesn't ; turn on by default, but it's here if you want to try it. if (keyword_set (itime_corr)) then begin nfrms = sxpar (header0, 'DCE_FRMS') if (nfrms eq 0) then begin print, 'Unable to determine exposure time because DCE_FRMS = 0' print, 'Defaulting to 3-second exposure' nfrms = 28 endif ratio = (float (nfrms) - 4.0) / 24.0 print, 'Performing exposure-time correction: ratio is ', ratio coeff[0] += coeff[1] * alog10 (ratio) threshold /= ratio endif ; now get image and (if present) uncertainty extensions if (input_fcb.nextend gt 0) then begin fits_read, input_fcb, image, header, exten_no = 1 if (input_fcb.nextend GT 1) then begin fits_read, input_fcb, image_unc, header_unc, exten_no = 2 endif endif ; compute correction new_image = image idx = where (image gt threshold, count) if (count gt 0) then begin new_image[idx] *= 10^(coeff[0] + coeff[1] * alog10 (image[idx])) endif ; compute uncertainty if (keyword_set (image_unc)) then begin new_image_unc = image_unc if (count gt 0) then begin new_image_unc[idx] = sqrt (unc_coeff[0]^2*(10^coeff[0]*image[idx]^(coeff[1]+1.0)*alog(10.0))^2 + $ unc_coeff[1]^2*(10^coeff[0]*image[idx]^(coeff[1]+1.0)*alog(image[idx]))^2 + $ image_unc[idx]^2*(10^coeff[0]*(coeff[1]+1.0)*image[idx]^coeff[1])^2) endif endif ; update image header sxaddpar,header0,'LINCOR70','T',' 70um linearity correction performed' sxaddpar,header0,'comment','linearity corr. made using mips70_lincor.pro v1.0'$ ,after='lincor70' ; write results to output file if (input_fcb.nextend GT 0) then begin fits_write,output_fcb,image0,header0 fits_write,output_fcb,new_image,header endif else begin fits_write,output_fcb,new_image,header0 endelse if (input_fcb.nextend GT 1) then begin fits_write,output_fcb,new_image_unc,header_unc endif ; copy the rest of the extensions in the image (if they exist) for i = 3,input_fcb.nextend do begin fits_read,input_fcb,image,header,exten_no=i fits_write,output_fcb,image,header endfor ; close files fits_close,input_fcb fits_close,output_fcb end