procedure display4k (image) string image ="" {prompt="name of MEF image"} string outimg ="tmp_display4k.fits" {prompt="name of output mosaic image"} string biassec ="[1022:1036,2:2031]" {prompt="overscan (overrides header)"} real gain10 =1.731 {prompt="gain zeropoint value, amp 3 [e-/ADU]"} real gain11 =-1.525e-6 {prompt="1st-order non-linearity term, amp 3"} real gain12 =9.346e-12 {prompt="2nd-order non-linearity term, amp 3"} real gain20 =1.782 {prompt="gain zeropoint value, amp 4 [e-/ADU]"} real gain21 =-1.388e-6 {prompt="1st-order non-linearity term, amp 4"} real gain22 =1.904e-11 {prompt="2nd-order non-linearity term, amp 4"} bool nonlcor =yes {prompt="apply non-linearity correction?"} bool telcomm =yes {prompt="is telemetry info present in header?"} bool dodispl =yes {prompt="display the output mosaic image?"} real zlow =INDEF {prompt="minimum display level [ADU]"} real zhigh =INDEF {prompt="maximum display level [ADU]"} bool zlog =no {prompt="display using logarithmic stretch?"} bool verbose =yes {prompt="verbose?"} # @(#) display4k Author: R.A. Jansen -- Oct 10 2007 # @(#) # @(#) task to combine the dual-amplifier images of the VATT 4K CCD, that # @(#) are written to disk in Multi-Extension FITS (MEF) format, into a # @(#) single flat FITS mosaic image, to display that image in an already # @(#) open 'ds9' and start 'imexam' for interactive examination of the # @(#) image. If "dodispl=no", the single flat FITS image is created but # @(#) not displayed. # @(#) # @(#) Note: this task requires tasks 'chkimg', 'cparse', 'rpparse ' and # @(#) 'rpbuild' from the "rjtools" external package, available at: # @(#) http://www.public.asu.edu/~rjansen/iraf/rjtools.html # @(#) If "rjtools" is not installed as a package, copy the above # @(#) tasks to your current working directory and declare them as # @(#) tasks individually, e.g.: task chkimg = chkimg.cl # @(#) # @(#) Last updated: # @(#) May 3 2008 -- updated the gathering of header info after version # @(#) change in 'AzCamtool' for May 3-7 2008 run. [RAJ] # @(#) May 21 2008 -- added 'biassec' task parameter to specify overscan # @(#) section that differs from the nominal one recorded # @(#) in the image header; rewrote code to generate the # @(#) output header; task can now serve as a first stage # @(#) in the reduction of VATT data. [RAJ] # @(#) Apr 2 2010 -- re-implemented a non-linearity correction assuming # @(#) a polynomial dependence of the gain on the average # @(#) signal level of order <=3. The gain zeropoints are # @(#) assumed to be appropriate for a level of 32750 ADU # @(#) and the output image will always have units of e-. # @(#) Package 'stsdas' is now required to be loaded; all # @(#) work-arounds using 'imarith' have been removed. Do # @(#) not use this version while observing at the VATT. begin string img, mos, osec, ima, imb, tmp, seca, secb string typ, hdr, bsec1, bsec2, cstr, ctmp int nb, nc, nr, mnc, mnr, tz, mhc, mhr, imot, il real eadu1[3], eadu2[3], eadu, zlo, zhi, rtmp real rlev1, rlev2, blev1, blev2 string ctyp, cdate, ctime, cobj, uts, ut1, ut2, cra, cdec, cst, cha string cf1, cf2, cf3, cf4, cfilt, co1, co2, co3, co4, observer real texp, tdrk, ct, dt, reqnx, repch, rjuln, ram, razim, relev real rutc, dgain, sky, sig # Make sure that the required package and tasks are loaded ... cstr = "Please load package '%s' first.\n" if ( !defpac("stsdas") ) { printf(cstr, "stsdas") ; return } if ( !defpac("artdata") ) { printf(cstr, "artdata") ; return } if ( !defpac("rjtools") ) { printf(cstr, "rjtools") ; return } # Initialize... img = image mos = outimg osec = biassec eadu1[1]= gain10 eadu1[2]= gain11 eadu1[3]= gain12 eadu2[1]= gain20 eadu2[2]= gain21 eadu2[3]= gain22 eadu = 0.5*(eadu1[1]+eadu2[1]) zlo = zlow zhi = zhigh ima = "tmp_display4k_ima.fits" imb = "tmp_display4k_imb.fits" tmp = "tmp_display4k_tmp.fits" hdr = "tmp_display4k.hdr" # Check whether the input MEF file exists... unlearn chkimg ; chkimg (img, "access", imtype="", verbose+) if (!chkimg.ok) { return } img = chkimg.impath//chkimg.imroot typ = chkimg.imtype # Check that output mosaic FITS file does not exist already; any file # with the name 'tmp_display4k', however, is silently overwritten... if ( mos == "tmp_display4k" || mos == "tmp_display4k.fits" ) { chkimg (mos, "noaccess", imtype="", verbose-) mos = chkimg.impath//chkimg.imroot//chkimg.imtype if (!chkimg.ok) { delete (mos, yes, verify-) } } else { chkimg (mos, "noaccess", imtype="", verbose+) if (!chkimg.ok) { return } } # Check that any user-supplied overscan section is a valid one... rpparse (str(osec), verbose=yes) if ( !rpparse.ok || rpparse.imrpstr=="[*,*]" ) { osec = "" } # If the task previously crashed, remove the debris... if ( access(ima) ) { delete (ima, yes, verify-) } if ( access(imb) ) { delete (imb, yes, verify-) } if ( access(hdr) ) { delete (hdr, yes, verify-) } if (verbose) { print ("DISPLAY4K[v2010.2]: image=\""//img//typ//"\" ...") } # Determine the dimensions of the individual MEF extensions... cstr = "i_naxis[12]" hselect (img//typ//"[1]", cstr, yes) | scan (nc,nr) # Determine the location of the overscan sections and number of # overscan columns... hselect (img//typ//"[1]", "biassec", yes) | scan (bsec1) hselect (img//typ//"[2]", "biassec", yes) | scan (bsec2) rpparse (bsec1, verbose=no) nb = rpparse.xhi - rpparse.xlo + 1 # Use user-specified overscan section instead header-declared one? if ( osec != "" ) { bsec1 = osec rpparse (osec, verbose=no) rpbuild ((nc - rpparse.xhi + 1), (rpparse.xhi - rpparse.xlo + 1), rpparse.ylo, rpparse.yhi, pln=0) bsec2 = rpbuild.imrpstr } # Construct the section specifiers for each of the two amplifier # images, given their relative orientations... rpbuild (1, (nc-nb), 1, nr, pln=0) seca = rpbuild.imrpstr rpbuild ((nc-nb+1), (2*(nc-nb)), 1, nr, pln=0) secb = rpbuild.imrpstr if (verbose) { printf(" Amp1 section =%20s Amp2 section=%20s\n", "'"//seca//"'", "'"//secb//"'") } # Extract the individual amplifier frames from the MEF file... if (verbose) { printf(" Extracting individual amplifier frames... ") } if (verbose) { printf(" [1]... ") } imcopy (img//typ//"[1][*,*]", ima, verbose-) if (verbose) { printf(" [2]... ") } imcopy (img//typ//"[2][-*,*]", imb, verbose-) if (verbose) { print ("done.") } sleep (1) # Measure the mean signal levels... if (verbose) { printf(" Measuring mean signal levels... ") } unlearn imstat if (verbose) { printf(" [1]... ") } imstat (ima, fields="midpt", nclip=5, lsigma=3., usigma=3., format-) | scan (rlev1) if (verbose) { printf(" [2]... ") } imstat (imb, fields="midpt", nclip=5, lsigma=3., usigma=3., format-) | scan (rlev2) if (verbose) { printf("\n") } printf(" Amp1 mean level = %9.2f ADU ", rlev1) printf(" Amp2 mean level = %9.2f ADU\n", rlev2) sleep (1) # Measure the bias levels in the overscan region... if (verbose) { printf(" Measuring bias levels... ") } if (verbose) { printf(" [1]... ") } imstat (ima//bsec1, fields="midpt", nclip=5, lsigma=3., usigma=3., format-) | scan (blev1) if (verbose) { printf(" [2]... ") } imstat (imb//bsec2, fields="midpt", nclip=5, lsigma=3., usigma=3., format-) | scan (blev2) if (verbose) { printf("\n") } printf(" Amp1 bias level = %9.2f ADU ", blev1) printf(" Amp2 bias level = %9.2f ADU\n", blev2) sleep (1) # Subtract the measured bias levels... if (verbose) { printf(" Subtracting bias level... ") } if (verbose) { printf(" [1]... ") } cstr = "im1 - "//blev1 imcalc (ima, tmp, cstr, pixtype="real", nullval=-33000., verbose-) delete (ima, yes, verify-) rename (tmp, ima, field="all") if (verbose) { printf(" [2]... ") } cstr = "im1 - "//blev2 imcalc (imb, tmp, cstr, pixtype="real", nullval=-33000., verbose-) delete (imb, yes, verify-) rename (tmp, imb, field="all") if (verbose) { print ("done.") } sleep (1) # Apply the specified polynomial non-linearity correction or simply # convert the image from ADU's to e- ... if (nonlcor) { if (verbose) { printf(" Applying non-linearity correction... ") } } else { if (verbose) { printf(" Converting ADU to e-... ") } } if (verbose) { printf(" [1]... ") } if (nonlcor) { rlev1 = rlev1 - 32750. cstr = "im1*("//str(eadu1[1])//" + " cstr = cstr//str(eadu1[2])//"*"//str(rlev1)//" + " cstr = cstr//str(eadu1[3])//"*("//str(rlev1)//")**2)" } else { cstr = "im1*"//str(eadu1[1]) } imcalc (ima, tmp, cstr, pixtype="real", nullval=-33000., verbose-) sleep (1) delete (ima, yes, verify-) rename (tmp, ima, field="all") if (verbose) { printf(" [2]... ") } if (nonlcor) { rlev2 = rlev2 - 32750. cstr = "im1*("//str(eadu2[1])//" + " cstr = cstr//str(eadu2[2])//"*"//str(rlev2)//" + " cstr = cstr//str(eadu2[3])//"*("//str(rlev2)//")**2)" } else { cstr = "im1*"//str(eadu2[1]) } imcalc (imb, tmp, cstr, pixtype="real", nullval=-33000., verbose-) sleep (1) delete (imb, yes, verify-) rename (tmp, imb, field="all") if (verbose) { print ("done.") } sleep (1) # Create an empty output mosaic image ... mnr = nr mnc = 2*(nc-nb) if (verbose) { printf(" Creating empty frame to hold mosaic... ") } mknoise (mos, output="", title=mos, ncols=mnc, nlines=mnr, header="artdata$stdheader.dat", background=0., rdnoise=0., ncosrays=0, comments=no) if (verbose) { print ("done.") } sleep (1) # Copy the individual amplifier sections into this mosaic... if (verbose) { printf(" Copying individual amplifier sections into mosaic... ") } rpbuild (1, (nc-nb), 1, nr, pln=0) if (verbose) { printf("[1]... ") } imcopy (ima//rpbuild.imrpstr, mos//seca, verbose-) rpbuild ((nb+1), nc, 1, nr, pln=0) if (verbose) { printf("[2]... ") } imcopy (imb//rpbuild.imrpstr, mos//secb, verbose-) if (verbose) { print ("done.") } sleep (1) # Attempt to get the orientation right (for Feb 2009/Mar 2010 runs)... imcopy (mos//"[*,-*]", mos//"[*,*]", verbose-) sleep (1) # Write a correctly populated and formatted FITS header (in case this # task is used with 'dodispl=no' to produce a mosaic for later use)... if (verbose) { printf(" Retrieving information from FITS header... ") } hselect (img//typ//"[0]", "imagetyp", yes) | scan (ctyp) hselect (img//typ//"[0]", "object", yes) | scan (cobj) hselect (img//typ//"[0]", "timesys", yes) | scan (uts) hselect (img//typ//"[0]", "timezone", yes) | scan (tz) hselect (img//typ//"[0]", "date-obs", yes) | scan (cdate) hselect (img//typ//"[0]", "time-obs", yes) | scan (ctime) hselect (img//typ//"[0]", "ut", yes) | scan (ut1) hselect (img//typ//"[0]", "exptime", yes) | scan (texp) hselect (img//typ//"[0]", "darktime", yes) | scan (tdrk) if (telcomm) { hselect (img//typ//"[0]", "ra", yes) | scan (cra) hselect (img//typ//"[0]", "dec", yes) | scan (cdec) hselect (img//typ//"[0]", "equinox", yes) | scan (reqnx) hselect (img//typ//"[0]", "epoch", yes) | scan (repch) hselect (img//typ//"[0]", "st", yes) | scan (cst) hselect (img//typ//"[0]", "ha", yes) | scan (cha) hselect (img//typ//"[0]", "julian", yes) | scan (rjuln) hselect (img//typ//"[0]", "airmass", yes) | scan (ram) hselect (img//typ//"[0]", "azimuth", yes) | scan (razim) hselect (img//typ//"[0]", "elevat", yes) | scan (relev) hselect (img//typ//"[0]", "motion", yes) | scan (imot) cf1 = "" ; cf2 = "" ; cf3 = "" ; cf4 = "" hselect (img//typ//"[0]", "filter", yes) | \ translit ("-", "\"", "", del+) | scan (cf1,cf2,cf3,cf4) } else { cra = '00:00:00.00' cdec = '+00:00:00.0' reqnx = 2000.0 repch = 2000.0 cst = '00:00:00.00' cha = '+00:00:00' ram = 1.0000 razim = 180.00 relev = 90.000 imot = 0 cf1 = 'TOP' cf2 = '1' cf3 = 'BOT' cf4 = '1' # The following will populate the header with the current JD; not # that when the frame was written at the telescope. Beware! gdate(utc=yes) rjuln = gdate.jd } co1 = "" ; co2 = "" ; co3 = "" ; co4 = "" hselect (img//typ//"[0]", "observer", yes) | \ translit ("-", "\"", "", del+) | scan (co1,co2,co3,co4) hselect (img//typ//"[0]", "detgain", yes) | scan (dgain) hselect (img//typ//"[0]", "camtemp", yes) | scan (ct) hselect (img//typ//"[0]", "dewtemp", yes) | scan (dt) if (verbose) { print ("done.") } # Prepare the header card image... if (verbose) { printf(" Writing correctly populated FITS header... ") } cstr = "Image type" printf("IMAGETYP= %-20s / %s\n", "'"//substr(ctyp//" ",1,8)//"'", cstr, > hdr) cstr = "Object name or title" il = strlen(cobj) if ( il < 8 ) { il = 8 } printf("OBJECT = %-20s / %s\n", "'"//substr(cobj//" ",1,il)//"'", cstr, >> hdr) cstr = "Time system" printf("TIMESYS = %-20s / %s\n", "'"//substr(uts//" ",1,8)//"'", cstr, >> hdr) cstr = "Local time zone (UT-" printf("TIMEZONE= %20d / %s%d)\n", tz, cstr, tz, >> hdr) cstr = "UT date and time at start of exposure" printf("OBSDATE = %-20s / %s\n", "'"//cdate//"T"//ctime//"'", cstr, >> hdr) cstr = "UT date at start of exposure" printf("DATE-OBS= %-20s / %s\n", "'"//cdate//"'", cstr, >> hdr) cstr = "UT time at start of exposure" printf("TIME-OBS= %-20s / %s\n", "'"//ctime//"'", cstr, >> hdr) cstr = "UT time at start of exposure" printf("UTOPEN = %-20s / %s\n", "'"//ut1//"'", cstr, >> hdr) cstr = "UT time at end of exposure" rutc = (real(ut1) + texp/3600.) if (rutc < 10.) { printf("0%-13.2h\n", rutc) | scan (ut2) } else { printf("%-013.2h\n", rutc) | scan (ut2) } printf("UTCLOSE = %-20s / %s\n", "'"//ut2//"'", cstr, >> hdr) cstr = "[s] Exposure time" printf("EXPTIME = %20.2f / %s\n", texp, cstr, >> hdr) cstr = "[s] Dark time" printf("DARKTIME= %20.2f / %s\n", tdrk, cstr, >> hdr) cstr = "[hh:mm:ss.ss] Right Ascension" printf("RA = %-20s / %s\n", "'"//cra//"'", cstr, >> hdr) cstr = "[Sdd:mm:ss.s] Declination" printf("DEC = %-20s / %s\n", "'"//cdec//"'", cstr, >> hdr) cstr = "Epoch of RA and DEC" printf("EPOCH = %20.1f / %s\n", repch, cstr, >> hdr) cstr = "Equinox of RA and DEC" printf("EQUINOX = %20.1f / %s\n", reqnx, cstr, >> hdr) cstr = "Local sidereal time at start of exposure" printf("LST = %-20s / %s\n", "'"//cst//"'", cstr, >> hdr) cstr = "Hour Angle at start of exposure" printf("HA = %-20s / %s\n", "'"//cha//"'", cstr, >> hdr) cstr = "Airmass at start of exposure" printf("AIRMASS = %20.3f / %s\n", ram, cstr, >> hdr) cstr = "Julian day at start of exposure" printf("JULIAN = %20.2f / %s\n", rjuln, cstr, >> hdr) cstr = "Modified Julian day at start of exposure" printf("MJD = %20.2f / %s\n", (rjuln-2400000.), cstr, >> hdr) cstr = "[deg] Azimuth at start of exposure" printf("AZIMUTH = %20.2f / %s\n", razim, cstr, >> hdr) cstr = "[deg] Elevation at start of exposure" printf("ELEVAT = %20.2f / %s\n", relev, cstr, >> hdr) cstr = "Telescope motion flag" printf("MOTION = %20d / %s\n", imot, cstr, >> hdr) cstr = "Observatory" printf("OBSERVAT= %-20s / %s\n", "'Vatican Observatory, Mt.Graham'", cstr, >> hdr) cstr = "IAU Observatory Code number" printf("IAU_ID = %20d / %s\n", 290, cstr, >> hdr) cstr = "Telescope name" printf("TELESCOP= %-20s / %s\n", "'VATT '", cstr, >> hdr) cstr = "[Sddd:mm:ss.ss] Site longitude, west of zero" printf("SITELONG= %-20s / %s\n", "'+109:53:31.25'", cstr, >> hdr) cstr = "[Sdd:mm:ss.ss] Site latitude" printf("SITELAT = %-20s / %s\n", "'+32:42:04.69'", cstr, >> hdr) cstr = "[m] Site elevation" printf("SITEELEV= %20.1f / %s\n", 3191.0, cstr, >> hdr) cstr = "Telescope focus" printf("TELFOCUS= %20d / %s\n", 0, cstr, >> hdr) cstr = "Telescope filter" printf("TELFILTE= %-20s / %s\n", "' '", cstr, >> hdr) cstr = "Instrument name" printf("INSTRUME= %-20s / %s\n", "'VATT4K '", cstr, >> hdr) cstr = "Instrument filter" cfilt = cf1//" "//cf2//" "//cf3//" "//cf4 printf("FILTER = %-20s / %s\n", "'"//cfilt//"'", cstr, >> hdr) cstr = "Number of detectors" printf("NCCDS = %20d / %s\n", 1, cstr, >> hdr) cstr = "Detector ID" printf("DETECTOR= %-20s / %s\n", "'STA0500A'", cstr, >> hdr) cstr = "Number of amplifiers" printf("NAMPS = %20d / %s\n", 2, cstr, >> hdr) cstr = "[um] Detector pixel size" printf("PIXSIZE = %20d / %s\n", 15, cstr, >> hdr) cstr = "[arcsec] Detector pixel size" printf("SECPIX = %20.4f / %s\n", 0.3754, cstr, >> hdr) cstr = "Commanded detector gain selection" printf("DETGAIN = %20d / %s\n", dgain, cstr, >> hdr) cstr = "Commanded CCD pixel binning" printf("CCDSUM = %-20s / %s\n", "'2 2 '", cstr, >> hdr) cstr = "[e-/ADU] Detector gain" printf("GAIN = %20.2f / %s\n", 1.00, cstr, >> hdr) cstr = "[e-] Read-noise" printf("RDNOISE = %20.2f / %s\n", 4.93, cstr, >> hdr) cstr = "[us] Sample integration time" printf("DWELL = %20d / %s\n", 0, cstr, >> hdr) cstr = "Camera temperature [degrees C]" printf("CAMTEMP = %20.1f / %s\n", ct, cstr, >> hdr) cstr = "Dewar temperature [degrees C]" printf("DEWTEMP = %20.1f / %s\n", dt, cstr, >> hdr) cstr = "Detector size" printf("DETSIZE = %-20s / %s\n", "'[1:2032,1:2032]'", cstr, >> hdr) cstr = "CCD section" printf("CCDSEC = %-20s / %s\n", "'[1:2032,1:2032]'", cstr, >> hdr) cstr = "Data section" printf("DATASEC = %-20s / %s\n", "'[1:2032,1:2032]'", cstr, >> hdr) cstr = "Trim section" printf("TRIMSEC = %-20s / %s\n", "'[1:2032,1:2032]'", cstr, >> hdr) cstr = "Bias section" printf("BIASSEC = %-20s / %s\n", "'[2033:2040,1:2032]'", cstr, >> hdr) cstr = "Program title" printf("PROGRAM = %-20s / %s\n", "' '", cstr, >> hdr) cstr = "Principle Investigator" printf("PROG_PI = %-20s / %s\n", "'R.A. Jansen'", cstr, >> hdr) cstr = "Observers" observer = co1 if ( co2 != "" ) { observer = observer//" "//co2 } if ( co3 != "" ) { observer = observer//" "//co3 } if ( co4 != "" ) { observer = observer//" "//co4 } printf("OBSERVER= %-20s / %s\n", "'"//observer//"'", cstr, >> hdr) cstr = "Original FITS file name" printf("ORIGFILE= %-20s / %s\n", "'"//img//typ//"'", cstr, >> hdr) # Write the header unit of the mosaic FITS image... mkheader (mos, hdr, append-, verbose-) if (verbose) { print ("done.") } # Display the output mosaic... mhc = 0.5*mnc ; mhr = 0.5*mnr rpbuild ((mhc-mhc/2.), (mhc+mhc/2.), (mhr-mhr/2.), (mhr+mhr/2.)) imstat (mos//rpbuild.imrpstr, fields="midpt,stddev", nclip=5, lsigma=4., usigma=3., format-) | scan (sky, sig) printf(" Mean signal level in mosaic: %-9s e-\n",str(sky)) if ( sky > 57500.*eadu ) { print ("WARNING: Frame may be OVEREXPOSED!") } if ( sky > 64500.*eadu ) { print ("WARNING: Frame is SATURATED! Bad astronomer!!") } if ( dodispl ) { print (" Displaying output mosaic image \""//mos//"\"...") if ( zlo == INDEF || zhi == INDEF ) { if ( zlo == INDEF ) { zlo = sky - 5*sig if ( sig > 5000. ) { zlo = sky - 5*5000. } } if ( zhi == INDEF ) { zhi = sky + 30*sig if ( sig > 5000. ) { zlo = sky + 30*5000. } } } if ( zlog ) { display (mos, 1, erase+, border_erase+, fill+, contrast=0., zrange-, zscale-, z1=zlo, z2=zhi, ztrans="log") } else { display (mos, 1, erase+, border_erase+, fill+, contrast=0., zrange-, zscale-, z1=zlo, z2=zhi, ztrans="linear") } #imexam print (" Good bye.") } else { print (" Created output mosaic image \""//mos//"\". Good bye.") } # Clean-up... delete (ima//","//imb//","//hdr, yes, verify-) sleep (1) flprcache sleep (1) end