diff options
Diffstat (limited to 'noao/onedspec/irsiids')
47 files changed, 5434 insertions, 0 deletions
diff --git a/noao/onedspec/irsiids/addsets.par b/noao/onedspec/irsiids/addsets.par new file mode 100644 index 00000000..1ad84e8b --- /dev/null +++ b/noao/onedspec/irsiids/addsets.par @@ -0,0 +1,8 @@ +# ADDSETS parameter file + +input,s,a,,,,Input image root file name +records,s,a,,,,Range of spectral records +output,s,a,,,,Output file root name for new spectra +start_rec,i,a,1,0,9999,Next starting spectral record +subset,i,h,2,1,,Number of spectra to add together within string +weighting,b,h,yes,,,Apply integration time weighting to calibrated data diff --git a/noao/onedspec/irsiids/batchred.cl b/noao/onedspec/irsiids/batchred.cl new file mode 100644 index 00000000..1bd5a3a3 --- /dev/null +++ b/noao/onedspec/irsiids/batchred.cl @@ -0,0 +1,168 @@ +#{ BATCHRED -- Script file to generate another script file +# which runs several ONEDSPEC tasks in an automated fashion. +# +# Currently the following procedures are automated: +# 1. STANDARD +# 2. SENSFUNC +# 3. BSWITCH +# 4. CALIBRATE +# 5. ADDSETS +# + +{ +# Say hello to the guy on the other side of the screen and check batch file. +print ("\n----B A T C H I I D S / I R S F I L E G E N E R A T O R----\n") + +s2 = "process.cl" # Batch file to be created. +if (access (s2)) { + print ("A batch file already exists - ") + if (query) + delete (s2, verify=no) +} + +# Initialize +rt = input # Root name for spectra +ot = output # Output root name +ttl = ">>&'" + ttylog + "')\n" # Log file for tty output + +out = rt +s1 = "" +st = "" +sns = "" +stat = "" +print ("i = ", start_rec, >>s2) + +if (standard) { # STANDARD? + print ("\n#---STANDARD---\n") + print ("\n#---STANDARD---\n", >>s2) + + st = std # STD file + if (access (st)) { + print (st, " - already exists") + if (query) + delete (st, verify=no) + } + + # Loop over all stars + b1 = yes + while (b1) { + # Check that the last entry was different - otherwise end input + records = "" + s3 = records + if (s3 == "") + b1 = no + else { + print ("standard (input='",rt,"',output='",st,"',",>>s2) + print ("\trecords='",s3,"',",>>s2) + print ("\tstar_name='",star_name,"',beam_switch=yes,",>>s2) + print ("\tsamestar=yes,apertures='',bandwidth=INDEF,",>>s2) + print ("\tbandsep=INDEF,interact=no,",ttl,>>s2) + } + } + + print ("") +} + + +if (sensfunc) { # SENSFUNC? + print ("\n#---SENSFUNC---\n") + print ("\n#---SENSFUNC---\n", >>s2) + + if (st == "") + st = std # STD file + sns = sensitivity # Sensitivity image + stat = stats # Statistics file + + print ("\nsensfunc (standards='",st,"',sensitivity='",sns,"',",>>s2) + print ("\tlogfile='",stat,"',apertures='',ignoreaps=no,",>>s2) + print ("\tfunction='",function,"',order=",order,",",>>s2) + print ("\tinteract=no,",ttl,>>s2) + + print ("") +} + + +if (bswitch) { # BSWITCH? + print ("\n#---BSWITCH---\n") + print ("\n#---BSWITCH---\n", >>s2) + + # Save starting output record number + in = out + out = "b" // ot + wt = weight # Weighting? + if (stat == "") + stat = stats # Statistics file + + # Accumulate records + print ("next_rec = i", >>s2) + b1 = yes + while (b1) { + records = "" + s3 = records + if (s3 == "") + b1 = no + else { + print ("j = next_rec", >> s2) + print ("bswitch (input='",in,"',output='",out,"',",>>s2) + print ("\trecords='",s3,"',stats='",stat,"',",>>s2) + print ("\tweighting=",wt,",subset=",subset,",",>>s2) + print ("\tstart_rec=j,",>>s2) + print ("\twave1=",wave1,",wave2=",wave2,",",ttl,>>s2) + } + } + + # Output records + print ("j = next_rec", >>s2) + s1 = "str (i) // '-' // str(j-1)" + print ("s1 = ", s1, >>s2) + + print ("") +} + +if (calibrate) { # CALIBRATE? + print ("\n#---CALIBRATE---\n") + print ("\n#---CALIBRATE---\n", >>s2) + + in = out + out = "c" // ot + if (sns == "") + sns = sensitivity # Sensivity file name + + if (s1 == "") { + records = "" + s1 = records + print ("s1 = '", s1, "'", >>s2) + } + + print ("calibrate (input='",in,"',output='",out,"',records=s1,",>>s2) + print ("\tignoreaps=no,",>>s2) + print ("\textinct=no,flux=yes,",>>s2) + print ("\tsensitivity='",sns,"',fnu=",fnu,",",ttl,>>s2) + + print ("") +} + +if (addsets) { # ADDSETS? + print ("\n#---ADDSETS---\n") + print ("\n#---ADDSETS---\n", >>s2) + + in = out + out = "a" // ot + if (s1 == "") { + records = "" + s1 = records + print ("s1 = '", s1, "'", >>s2) + } + + print ("addsets (input='",in,"',output='",out,"',records=s1,",>>s2) + print ("\tstart_rec=i,subset=2,",ttl,>>s2) +} + +# All done with generator. Ask whether to execute it. +print ("File generation complete - filename=",s2) +if (proceed == no) # Execute batch file? + bye +} + +# Execute generated batch file +process & diff --git a/noao/onedspec/irsiids/batchred.par b/noao/onedspec/irsiids/batchred.par new file mode 100644 index 00000000..dbb61b8c --- /dev/null +++ b/noao/onedspec/irsiids/batchred.par @@ -0,0 +1,38 @@ +# BATCHRED -- Parameter file for batch reduction prep task + +input,s,a,,,,Input root name for spectra +output,s,a,,,,Output root name for spectra +start_rec,i,a,,0,9999,Starting record for output spectra +ttylog,s,a,"ttylog",,,File name to contain a log of terminal output +standard,b,a,yes,,,Generate commands for STANDARD +sensfunc,b,a,yes,,,Generate commands for SENSFUNC +bswitch,b,a,yes,,,Generate commands for BSWITCH +calibrate,b,a,yes,,,Generate commands for CALIBRATE +addsets,b,a,yes,,,Generate commands for ADDSETS + +std,s,a,"std",,,STANDARD and SENSFUNC standard star file +star_name,s,q,,,,STANDARD star name +stats,s,a,"stats",,,SENSFUNC and BSWITCH statistics file +sensitivity,s,a,sens,,,SENSFUNC and CALIBRATE sensitivity spectra +weight,b,a,no,,,BSWITCH weighted averages? +function,s,h,"chebyshev",,,SENSFUNC fitting function +order,i,h,7,1,,SENSFUNC fitting order + +records,s,q,,,,Record string to process +proceed,b,q,yes,,,Begin batch processing +query,b,q,no,,,Delete files(s)? + +fnu,b,h,no +wave1,r,h,0.0 +wave2,r,h,0.0 +subset,i,h,32767 + +rt,s,h +ot,s,h +in,s,h +out,s,h +stat,s,h +sns,s,h +st,s,h +wt,b,h +ttl,s,h diff --git a/noao/onedspec/irsiids/bplot.cl b/noao/onedspec/irsiids/bplot.cl new file mode 100644 index 00000000..53ce4cfc --- /dev/null +++ b/noao/onedspec/irsiids/bplot.cl @@ -0,0 +1,35 @@ +# BPLOT -- Batch plotting of spectra with SPLOT + +procedure bplot (images, records) + +string images {prompt="List of images to plot"} +string records = "" {prompt="List of records to plot"} +string graphics = "stdgraph" {prompt="Graphics output device"} +string cursor = "onedspec$gcurval.dat" {prompt="Cursor file(s)\n"} + +struct *ilist, *clist + +begin + int line, ap + file ifile, cfile, cur, image + + ifile = mktemp ("bplot") + cfile = mktemp ("bplot") + + names (images, records, >& ifile) + files (cursor, > cfile) + cur = "" + + ilist = ifile; clist = cfile + while (fscan (ilist, image) != EOF) { + if ((cursor != "") && (fscan (clist, cur) == EOF)) { + clist = cfile + line = fscan (clist, cur) + } + splot (image, graphics=graphics, cursor=cur) + } + clist = ""; ilist = "" + + delete (ifile, verify=no) + delete (cfile, verify=no) +end diff --git a/noao/onedspec/irsiids/bswitch.par b/noao/onedspec/irsiids/bswitch.par new file mode 100644 index 00000000..3a8b7e8f --- /dev/null +++ b/noao/onedspec/irsiids/bswitch.par @@ -0,0 +1,15 @@ +# BSWITCH parameter file + +input,s,a,,,,Input spectra file root name +records,s,a,,,,Ranges of spectral records +output,s,a,,,,Output file root name for new spectra +start_rec,i,a,1,0,9999,Next starting spectral record +stats,s,a,"stats",,,File to contain statistics summary +ids_mode,b,h,yes,,,Are data in quadruples +extinct,b,h,yes,,,Apply de-extinction correction +weighting,b,h,no,,,Apply statistical weights during addition +subset,i,h,32767,1,,Generate sums at subset intervals +wave1,r,h,0.0,,,Starting wavelength to accumulate stats +wave2,r,h,0.0,,,Ending wavelength to accumulate stats +observatory,s,h,"kpno",,,Observatory of data +extinction,s,h,)_.extinction,,,Extinction file diff --git a/noao/onedspec/irsiids/coefs.par b/noao/onedspec/irsiids/coefs.par new file mode 100644 index 00000000..dbabec65 --- /dev/null +++ b/noao/onedspec/irsiids/coefs.par @@ -0,0 +1,3 @@ +input,s,a,,,,Input image root file name +records,s,a,,,,Range of spectral records +database,s,h,"database",,,IDENTIFY database diff --git a/noao/onedspec/irsiids/coincor.par b/noao/onedspec/irsiids/coincor.par new file mode 100644 index 00000000..a0f2e0bb --- /dev/null +++ b/noao/onedspec/irsiids/coincor.par @@ -0,0 +1,9 @@ +# COINCOR parameter file + +input,s,a,,,,Input image root file name +records,s,a,,,,Range of spectral records +output,s,a,,,,Output file root name for new spectra +start_rec,i,a,1,0,9999,Next starting spectral record +ccmode,s,h,)_.ccmode,,,Correction mode (photo|iids) +deadtime,r,h,)_.deadtime,,,Deadtime in seconds +power,r,h,)_.power,,,IIDS power law coefficient diff --git a/noao/onedspec/irsiids/coincor.x b/noao/onedspec/irsiids/coincor.x new file mode 100644 index 00000000..df572c02 --- /dev/null +++ b/noao/onedspec/irsiids/coincor.x @@ -0,0 +1,123 @@ +# Coincidence correction options +define CC_PHOTO_MODE 1 # Photometer style correction +define CC_IIDS_MODE 2 # IIDS style +define CC_POWER_MODE 3 # Power law correction +define CC_USER_MODE 4 # User supplies a function + + +# COINCOR -- Coincidence correction for detector deadtime + +procedure coincor (input, output, npts, expo, coflag, dt, power, mode) + +real input[npts] +real output[npts] +real expo +int coflag +real dt +real power +int mode, npts + +begin + # Check that exposure time is legit + if (expo <= 0.0) + return + + # Select the method by which the correction is performed + switch (mode) { + case CC_PHOTO_MODE: + # Photoelectric photometer + call ccphoto (input, output, npts, expo, coflag, dt) + + case CC_IIDS_MODE: + # IIDS style correction + if (coflag == -1) { + call cciids (input, output, npts, expo, coflag, dt) + if (power != 1.0) + call ccpower (output, output, npts, expo, coflag, + power) + } else if ((coflag == 0) && (power != 1.0)) + call ccpower (input, output, npts, expo, coflag, power) + else + call amovr (input, output, npts) + + case CC_USER_MODE: + # Provided by the user + call ccuser (input, output, npts, expo, coflag, dt) + } +end + +# CCPHOTO -- Photoelectric photometer coincidence correction + +procedure ccphoto (input, output, npts, expo, coflag, dt) + +real input[npts], output[npts], expo, dt +int coflag +int npts + +int i + +begin + do i = 1, npts + output[i] = input[i] * exp (input[i] * dt / expo) + coflag = 2 +end + +# CCUSER -- User supplied correction scheme + +procedure ccuser (input, output, npts, expo, coflag, dt) + +real input[npts], output[npts], expo, dt +int coflag +int npts + +begin + coflag = 3 +end + +# CCIIDS -- IIDS style correction scheme +# From Instrumentation for Astronomy III (SPIE Vol 172) p.88 by Larry Goad +# +# Note that only the "Detect" mode of observation is supported. + +procedure cciids (input, output, npts, expo, coflag, dt) + +real input[npts], output[npts], expo, dt +int npts, coflag + +int i +real tsweep, value + +begin + # Allow tsweep to be the deadtime so that a different value + # may be entered for other instruments. + # For the IIDS, tsweep = 1.424e-3 sec + tsweep = dt + + do i = 1, npts { + value = 1 - input[i] / expo * tsweep + if ((value < 0.) || (value > 1.)) + output[i] = input[i] + else + output[i] = -expo * log (value)/ tsweep + } + coflag = 0 +end + +# CCPOWER -- Power law correction +# Power law correction from Massey and De Veny, NOAO Newsletter #6. + +procedure ccpower (input, output, npts, expo, coflag, power) + +real input[npts], output[npts], expo, power +int npts, coflag + +int i + +begin + do i = 1, npts + if (input[i] > 0.) + output[i] = expo * (input[i] / expo) ** power + else + output[i] = input[i] + coflag = 1 +end diff --git a/noao/onedspec/irsiids/conversion.x b/noao/onedspec/irsiids/conversion.x new file mode 100644 index 00000000..9d0b8c15 --- /dev/null +++ b/noao/onedspec/irsiids/conversion.x @@ -0,0 +1,213 @@ +define MAX_CHARS 256 + + +# ASCII_TO_EBCDIC -- Vector procedure to convert ASCII characters to EBCDIC +# characters using the lookup table atoe. + +procedure ascii_to_ebcdic (inbuffer, outbuffer, nchars) + +char inbuffer[ARB] +short outbuffer[ARB], atoe[MAX_CHARS] +int l, nchars + +data (atoe[l], l = 1, 8) / 0b, 1b, 2b, 3b, '7' , '-' , '.' , '/' / +data (atoe[l], l = 9, 16) /26b, 5b, '%' , 13b, 14b, 15b, 16b, 17b / +data (atoe[l], l = 17, 24) /20b, 21b, 22b, 23b, '<' , '=' , '2' , '&' / +data (atoe[l], l = 25, 32) /30b, 31b, '?' , '\'', 34b, 35b, 36b, 37b / +data (atoe[l], l = 33, 40) /'@' , 'O' , 177b, '{' , '[' , 'l' , 'P' , '}' / +data (atoe[l], l = 41, 48) /'M' , ']' , '\\' , 'N' , 'k' , '`' , 'K' , 'a'/ +data (atoe[l], l = 49, 56) /360b, 361b, 362b, 363b, 364b, 365b, 366b, 367b/ +data (atoe[l], l = 57, 64) /370b, 371b, 'z' , '^' , 'L' , '~' , 'n' , 'o' / +data (atoe[l], l = 65, 72) /'|' , 301b, 302b, 303b, 304b, 305b, 306b, 307b/ +data (atoe[l], l = 73, 80) /310b, 311b, 321b, 322b, 323b, 324b, 325b, 326b/ +data (atoe[l], l = 81, 88) /327b, 330b, 331b, 342b, 343b, 344b, 345b, 346b/ +data (atoe[l], l = 89, 96) /347b, 350b, 351b, 'J' , 340b, 'Z' , '_' , 'm' / +data (atoe[l], l = 97, 104) /'y' , 201b, 202b, 203b, 204b, 205b, 206b, 207b/ +data (atoe[l], l = 105, 112) /210b, 211b, 221b, 222b, 223b, 224b, 225b, 226b/ +data (atoe[l], l = 113, 120) /227b, 230b, 231b, 242b, 243b, 244b, 245b, 246b/ +data (atoe[l], l = 121, 128) /247b, 250b, 251b, 300b, 'j' , 320b, 241b, 7b/ +data (atoe[l], l = 129, 136) /' ' , '!' , '"' , '#' , '$' , 25b, 6b, 27b/ +data (atoe[l], l = 137, 144) /'(' , ')' , '*' , '+' , ',' , 11b, 12b, 33b/ +data (atoe[l], l = 145, 152) /'0' , '1' , 32b, '3' , '4' , '5' , '6' , 10b/ +data (atoe[l], l = 153, 160) /'8' , '9' , ':' , ';' , 4b, 24b, '>' , 341b/ +data (atoe[l], l = 161, 168) /'A' , 'B' , 'C' , 'D' , 'E' , 'F' , 'G' , 'H' / +data (atoe[l], l = 169, 176) /'I' , 'Q' , 'R' , 'S' , 'T' , 'U' , 'V' , 'W' / +data (atoe[l], l = 177, 184) /'X' , 'Y' , 'b' , 'c' , 'd' , 'e' , 'f' , 'g' / +data (atoe[l], l = 185, 192) /'h' , 'i' , 'p' , 'q' , 'r' , 's' , 't' , 'u' / +data (atoe[l], l = 193, 200) /'v' , 'w' , 'x' , 200b, 212b, 213b, 214b, 215b/ +data (atoe[l], l = 201, 208) /216b, 217b, 220b, 232b, 233b, 234b, 235b, 236b/ +data (atoe[l], l = 209, 216) /237b, 240b, 252b, 253b, 254b, 255b, 256b, 257b/ +data (atoe[l], l = 217, 224) /260b, 261b, 262b, 263b, 264b, 265b, 266b, 267b/ +data (atoe[l], l = 225, 232) /270b, 271b, 272b, 273b, 274b, 275b, 276b, 277b/ +data (atoe[l], l = 233, 240) /312b, 313b, 314b, 315b, 316b, 317b, 332b, 333b/ +data (atoe[l], l = 241, 248) /334b, 335b, 336b, 337b, 352b, 353b, 354b, 355b/ +data (atoe[l], l = 249, 256) /356b, 357b, 372b, 373b, 374b, 375b, 376b, 377b/ + +begin + call alutcs (inbuffer, outbuffer, nchars, atoe) +end + +# EBCDIC_TO_ASCII -- Vector procedure to convert EBCDIC characters to ASCII +# characters. + +procedure ebcdic_to_ascii (inbuffer, outbuffer, nchars) + +char outbuffer[ARB] +short inbuffer[ARB], etoa[MAX_CHARS] +int l, nchars + +data (etoa[l], l = 1, 8) / 0b, 1b, 2b, 3b, 234b, 11b, 206b, 177b / +data (etoa[l], l = 9, 16) /227b, 215b, 216b, 13b, 14b, 15b, 16b, 17b/ +data (etoa[l], l = 17, 24) /20b, 21b, 22b, 23b, 235b, 205b, 10b, 207b / +data (etoa[l], l = 25, 32) /30b, 31b, 222b, 217b, 34b, 35b, 36b, 37b / +data (etoa[l], l = 33, 40) /200b, 201b, 202b, 203b, 204b, 12b, 27b, 33b/ +data (etoa[l], l = 41, 48) /210b, 211b, 212b, 213b, 214b, 5b, 6b, 7b/ +data (etoa[l], l = 49, 56) /220b, 221b, 26b, 223b, 224b, 225b, 226b, 4b/ +data (etoa[l], l = 57, 64) /230b, 231b, 232b, 233b, 24b, 25b, 236b, 32b/ +data (etoa[l], l = 65, 72) /' ' , 240b, 241b, 242b, 243b, 244b, 245b, 246b/ +data (etoa[l], l = 73, 80) /247b, 250b, '[' , '.' , '<' , '(' , '+' , '!' / +data (etoa[l], l = 81, 88) /'&' , 251b, 252b, 253b, 254b, 255b, 256b, 257b/ +data (etoa[l], l = 89, 96) /260b, 261b, ']' , '$' , '*' , ')' , ';' , '^' / +data (etoa[l], l = 97, 104) /'-' , '/' , 262b, 263b, 264b, 265b, 266b, 267b/ +data (etoa[l], l = 105, 112) /270b, 271b, '|' , ',' , '%' , '_' , '>' , '?' / +data (etoa[l], l = 113, 120) /272b, 273b, 274b, 275b, 276b, 277b, 300b, 301b/ +data (etoa[l], l = 121, 128) /302b, '`' , ':' , '#' , '@' , '\'' , '=' , '"'/ +data (etoa[l], l = 129, 136) /303b, 'a' , 'b' , 'c' , 'd' , 'e' , 'f' , 'g' / +data (etoa[l], l = 137, 144) /'h' , 'i' , 304b, 305b, 306b, 307b, 310b, 311b/ +data (etoa[l], l = 145, 152) /312b, 'j' , 'k' , 'l' , 'm' , 'n' , 'o' , 'p' / +data (etoa[l], l = 153, 160) /'q' , 'r' , 313b, 314b, 315b, 316b, 317b, 320b/ +data (etoa[l], l = 161, 168) /321b, '~' , 's' , 't' , 'u' , 'v' , 'w' , 'x' / +data (etoa[l], l = 169, 176) /'y' , 'z' , 322b, 323b, 324b, 325b, 326b, 327b/ +data (etoa[l], l = 177, 184) /330b, 331b, 332b, 333b, 334b, 335b, 336b, 337b/ +data (etoa[l], l = 185, 192) /340b, 341b, 342b, 343b, 344b, 345b, 346b, 347b/ +data (etoa[l], l = 193, 200) /'{' , 'A' , 'B' , 'C' , 'D' , 'E' , 'F' , 'G' / +data (etoa[l], l = 201, 208) /'H' , 'I' , 350b, 351b, 352b, 353b, 354b, 355b/ +data (etoa[l], l = 209, 216) /'}' , 'J' , 'K' , 'L' , 'M' , 'N' , 'O' , 'P' / +data (etoa[l], l = 217, 224) /'Q' , 'R' , 356b, 357b, 360b, 361b, 362b, 363b/ +data (etoa[l], l = 225, 232) /'\\', 237b, 'S' , 'T' , 'U' , 'V' , 'W' , 'X' / +data (etoa[l], l = 233, 240) /'Y' , 'Z' , 364b, 365b, 366b, 367b, 370b, 371b/ +data (etoa[l], l = 241, 248) /'0' , '1' , '2' , '3' , '4' , '5' , '6' , '7' / +data (etoa[l], l = 249, 256) /'8' , '9' , 372b, 373b, 374b, 375b, 376b, 377b/ + +begin + call alutsc (inbuffer, outbuffer, nchars, etoa) +end + +# IBM_TO_ASCII -- Vector procedure for converting IBM characters to ASCII +# characters. + +procedure ibm_to_ascii (inbuffer, outbuffer, nchars) + +char outbuffer[ARB] +short inbuffer[ARB], ibmtoa[MAX_CHARS] +int l, nchars + +data (ibmtoa[l], l = 1, 8) /0b, 1b, 2b, 3b, 234b, 11b, 206b, 177b / +data (ibmtoa[l], l = 9, 16) /1227b, 215b, 216b, 13b, 14b, 15b, 16b, 17b/ +data (ibmtoa[l], l = 17, 24) /20b, 21b, 22b, 23b, 235b, 205b, 10b, 207b / +data (ibmtoa[l], l = 25, 32) /30b, 31b, 222b, 217b, 34b, 35b, 36b, 37b / +data (ibmtoa[l], l = 33, 40) /200b, 201b, 202b, 203b, 204b, 12b, 27b, 33b/ +data (ibmtoa[l], l = 41, 48) /210b, 211b, 212b, 213b, 214b, 5b, 6b, 7b/ +data (ibmtoa[l], l = 49, 56) /220b, 221b, 26b, 223b, 224b, 225b, 226b, 4b/ +data (ibmtoa[l], l = 57, 64) /230b, 231b, 232b, 233b, 24b, 25b, 236b, 32b/ +data (ibmtoa[l], l = 65, 72) /' ' , 240b, 241b, 242b, 243b, 244b, 245b, 246b/ +data (ibmtoa[l], l = 73, 80) /247b, 250b, 0b, '.' , '<' , '(' , '+' , '|' / +data (ibmtoa[l], l = 81, 88) /'&' , 251b, 252b, 253b, 254b, 255b, 256b, 257b/ +data (ibmtoa[l], l = 89, 96) /260b, 261b, '!' , '$' , '*' , ')' , ';' , '^' / +data (ibmtoa[l], l = 97, 104) /'-' , '/' , 262b, 263b, 264b, 265b, 266b, 267b/ +data (ibmtoa[l], l = 105,112) /270b, 271b, 0b, ',' , '%' , '_' , '>' , '?' / +data (ibmtoa[l], l = 113, 120) /272b, 273b, 274b, 275b, 276b, 277b, 300b, 301b/ +data (ibmtoa[l], l = 121, 128) /302b, '`' , ':' , '#' , '@' , '\'' , '=' , '"'/ +data (ibmtoa[l], l = 129, 136) /303b, 'a' , 'b' , 'c' , 'd' , 'e' , 'f' , 'g' / +data (ibmtoa[l], l = 137, 144) /'h' , 'i' , 304b, 305b, 306b, 307b, 310b, 311b/ +data (ibmtoa[l], l = 145, 152) /312b, 'j' , 'k' , 'l' , 'm' , 'n' , 'o' , 'p' / +data (ibmtoa[l], l = 153, 160) /'q' , 'r' , 313b, 314b, 315b, 316b, 317b, 320b/ +data (ibmtoa[l], l = 161, 168) /321b, '~' , 's' , 't' , 'u' , 'v' , 'w' , 'x' / +data (ibmtoa[l], l = 169, 176) /'y' , 'z' , 322b, 323b, 324b, 325b, 326b, 327b/ +data (ibmtoa[l], l = 177, 184) /330b, 331b, 332b, 333b, 334b, 335b, 336b, 337b/ +data (ibmtoa[l], l = 185, 192) /340b, 341b, 342b, 343b, 344b, 345b, 346b, 347b/ +data (ibmtoa[l], l = 193, 200) /'{' , 'A' , 'B' , 'C' , 'D' , 'E' , 'F' , 'G' / +data (ibmtoa[l], l = 201, 208) /'H' , 'I' , 350b, 351b, 352b, 353b, 354b, 355b/ +data (ibmtoa[l], l = 209, 216) /'}' , 'J' , 'K' , 'L' , 'M' , 'N' , 'O' , 'P' / +data (ibmtoa[l], l = 217, 224) /'Q' , 'R' , 356b, 357b, 360b, 361b, 362b, 363b/ +data (ibmtoa[l], l = 225, 232) /'\\', 237b, 'S' , 'T' , 'U' , 'V' , 'W' , 'X' / +data (ibmtoa[l], l = 233, 240) /'Y' , 'Z' , 364b, 365b, 366b, 367b, 370b, 371b/ +data (ibmtoa[l], l = 241, 248) /'0' , '1' , '2' , '3' , '4' , '5' , '6' , '7' / +data (ibmtoa[l], l = 249, 256) /'8' , '9' , 372b, 373b, 374b, 375b, 376b, 377b / + +begin + call alutsc (inbuffer, outbuffer, nchars, ibmtoa) +end + +# Vector procedure to convert ASCII characters to ibm characters + +procedure ascii_to_ibm (inbuffer, outbuffer, nchars) + +char inbuffer[ARB] +short outbuffer[ARB], atoibm[MAX_CHARS] +int l, nchars + +data (atoibm[l], l = 1, 8) /0b, 1b, 2b, 3b, '7' , '-' , '.' , '/' / +data (atoibm[l], l = 9, 16) /26b, 5b, '%' , 13b, 14b, 15b, 16b, 17b / +data (atoibm[l], l = 17, 24) /20b, 21b, 22b, 23b, '<' , '=' , '2' , '&' / +data (atoibm[l], l = 25, 32) /30b, 31b, '?' , '\'', 34b, 35b, 36b, 37b / +data (atoibm[l], l = 33, 40) /'@' , 'Z' , 177b, '{' , '[' , 'l' , 'P' , '}' / +data (atoibm[l], l = 41, 48) /'M' , ']' , '\\', 'N' , 'k' , '`' , 'K' , 'a' / +data (atoibm[l], l = 49, 56) /360b, 361b, 362b, 363b, 364b, 365b, 366b, 367b/ +data (atoibm[l], l = 57, 64) /370b, 371b, 'z' , '^' , 'L' , '~' , 'n' , 'o' / +data (atoibm[l], l = 65, 72) /'|' , 301b, 302b, 303b, 304b, 305b, 306b, 307b/ +data (atoibm[l], l = 73, 80) /310b, 311b, 321b, 322b, 323b, 324b, 325b, 326b/ +data (atoibm[l], l = 81, 88) /327b, 330b, 331b, 342b, 343b, 344b, 345b, 346b/ +data (atoibm[l], l = 89, 96) /347b, 350b, 351b, 255b, 340b, 275b, '_' , 'm' / +data (atoibm[l], l = 97, 104) /'y' , 201b, 202b, 203b, 204b, 205b, 206b, 207b/ +data (atoibm[l], l = 105, 112) /210b, 211b, 221b, 222b, 223b, 224b, 225b, 226b/ +data (atoibm[l], l = 113, 120) /227b, 230b, 231b, 242b, 243b, 244b, 245b, 246b/ +data (atoibm[l], l = 121, 128) /247b, 250b, 251b, 300b, 'O' , 320b, 241b, 7b/ +data (atoibm[l], l = 129, 136) /' ' , '!' , '"' , '#' , '$' , 25b, 6b, 27b/ +data (atoibm[l], l = 137, 144) /'(' , ')' , '*' , '+' , ',' , 11b, 12b, 33b/ +data (atoibm[l], l = 145, 152) /'0' , '1' , 32b, '3' , '4' , '5' , '6' , 10b/ +data (atoibm[l], l = 153, 160) /'8' , '9' , ':' , ';' , 4b, 24b, '>' , 341b/ +data (atoibm[l], l = 161, 168) /'A' , 'B' , 'C' , 'D' , 'E' , 'F' , 'G' , 'H' / +data (atoibm[l], l = 169, 176) /'I' , 'Q' , 'R' , 'S' , 'T' , 'U' , 'V' , 'W' / +data (atoibm[l], l = 177, 184) /'X' , 'Y' , 'b' , 'c' , 'd' , 'e' , 'f' , 'g' / +data (atoibm[l], l = 185, 192) /'h' , 'i' , 'p' , 'q' , 'r' , 's' , 't' , 'u' / +data (atoibm[l], l = 193, 200) /'v' , 'w' , 'x' , 200b, 212b, 213b, 214b, 215b/ +data (atoibm[l], l = 201, 208) /216b, 217b, 220b, 232b, 233b, 234b, 235b, 236b/ +data (atoibm[l], l = 209, 216) /237b, 240b, 252b, 253b, 254b, 255b, 256b, 257b/ +data (atoibm[l], l = 217, 224) /260b, 261b, 262b, 263b, 264b, 265b, 266b, 267b/ +data (atoibm[l], l = 225, 232) /270b, 271b, 272b, 273b, 274b, 275b, 276b, 277b/ +data (atoibm[l], l = 233, 240) /312b, 313b, 314b, 315b, 316b, 317b, 332b, 333b/ +data (atoibm[l], l = 241, 248) /334b, 335b, 336b, 337b, 352b, 353b, 354b, 355b/ +data (atoibm[l], l = 249, 256) /356b, 357b, 372b, 373b, 374b, 375b, 376b, 377b/ + +begin + call alutcs (inbuffer, outbuffer, nchars, atoibm) +end + +# ALUTSC -- Vector operator to map one set of characters to another using a +# lookup table. + +procedure alutsc (a, b, nchars, lut) + +char b[nchars] +int nchars, i +short a[nchars], lut[ARB] + +begin + do i = 1, nchars, 1 + b[i] = lut[a[i] + 1] +end + +# ALUTCS -- Vector operator to map one set of characters to another using +# a lookup table. + +procedure alutcs (a, b, nchars, lut) + +char a[nchars] +int nchars, i +short b[nchars], lut[ARB] + +begin + do i = nchars, 1, -1 + b[i] = lut[a[i] + 1] +end diff --git a/noao/onedspec/irsiids/doc/addsets.hlp b/noao/onedspec/irsiids/doc/addsets.hlp new file mode 100644 index 00000000..6ce49122 --- /dev/null +++ b/noao/onedspec/irsiids/doc/addsets.hlp @@ -0,0 +1,66 @@ +.help addsets Feb85 noao.imred.iids/noao.imred.irs +.ih +NAME +addsets - Add subsets of a string of spectra +.ih +USAGE +addsets input records +.ih +PARAMETERS +.ls input +The root file name for the input spectra in the string. +.le +.ls records +The range of spectra indicating the elements of the string. +The names of the spectra will be formed by appending the range +elements to the input root name. +.le +.ls output +This is the root file name for the names of the spectra which will +be created by the addset operation. +.le +.ls start_rec = 1 +The starting record number to be appended to the root name of the +created spectra. +.le +.ls subset = 2 +The length of the substring of spectra which will be added together. +For IIDS/IRS data which has been processed through BSWITCH, this +parameter should be 2. This implies that spectra will be taken +2 at a time, added, and the sum written as a new spectrum. +.le +.ls weighting = yes +If set to yes, an average of the substring of spectra is generated +(if flux calibrated) weighted by the integration times of the +individual spectra. If set to no, a simple average is generated. +If not flux calibrated, this parameter has no effect - a simple +sum is generated. +.le +.ih +DESCRIPTION +Every "subset" group of spectra will be accumulated and the sum will be +written as a new spectrum. For example, if the input string contains +100 spectra, and subset=2, then 50 new spectra will be created. Each +new spectrum will be the sum of the consecutive pairs in the original string. + +If there are insufficient spectra to complete a subset accumulation, +the sum is written out anyway and a warning printed. For example, +if the input string contains 23 spectra, and subset=4, there will be +6 new spectra created, but the last one will be based on only 3 spectra. + +Subset may be set to 1 to allow a copy operation although this is not +a very efficient way to do so. +.ih +EXAMPLES +The following three examples are those described above. + +.nf + cl> addsets nite1 2001-2100 + cl> addsets nite1 2001-2023 subset=4 + cl> addsets nite1 2001-2010 subset=1 output=nite2 \ + >>> start_rec=2001 +.fi +.ih +SEE ALSO +bswitch +.endhelp diff --git a/noao/onedspec/irsiids/doc/batchred.hlp b/noao/onedspec/irsiids/doc/batchred.hlp new file mode 100644 index 00000000..9301f8b0 --- /dev/null +++ b/noao/onedspec/irsiids/doc/batchred.hlp @@ -0,0 +1,145 @@ +.help batchred Feb85 noao.imred.iids/noao.imred.irs +.ih +NAME +batchred - Automated processing of IIDS/IRS spectra +.ih +USAGE +batchred +.ih +PARAMETERS +This script task has many parameters, but most are used as +variables internal to the task and are not user parameters. +There are 5 parameters having similar purposes: standard, +sensfunc, bswitch, calibrate, and addsets. Each corresponds +to the ONEDSPEC task of the same name and BATCHRED will generate +the commands necessary to invoke those tasks if the associated +parameter is set to yes (the default in all cases). + +.ls standard = yes +.le +.ls sensfunc = yes +.le +.ls bswitch = yes +.le +.ls calibrate = yes +.le +.ls addsets = yes +.le +.ls fnu = no +This parameter is identical to the fnu parameter for CALIBRATE. +.le +.ls wave1 = 0.0 +This parameter is identical to the wave1 parameter for BSWITCH. +.le +.ls wave2 = 0.0 +This parameter is identical to the wave2 parameter for BSWITCH. +.le +.ls subset = 32767 +This parameter is identical to the subset parameter for BSWITCH. +.le +.ih +DESCRIPTION +Through a question and answer session, a series of commands to +ONEDSPEC is generated which are then processed as a batch job +to reduce "typical" spectra from the IIDS and IRS spectrographs. + +By setting the appropriate hidden parameters, the user may +"turn off" command generation for any of the possible tasks. + +A script task is generated having the name "process.cl" which is +submitted to the CL as the final command of BATCHRED. +All terminal output which would normally appear during the course +of running each of the individual tasks is redirected to a log file +(default=ttylog). + +After the script has been generated, the user may suppress running +the processing task. The script file remains on disk so that subsequent +cases may be appended, such as when +several independent runs of data are to be processed in one +stream (e.g. several nights of data, each to be reduced separately). + +The questions which are asked are described below: + +"Root name for spectra file names:" This is the input root file name +for all spectra which will be run through STANDARD and BSWITCH. + +"Root name for spectra to be created:" This is the output root file +name which all newly created spectra will use. It is also the +input file name for tasks CALIBRATE and ADDSETS since these tasks +operate on spectra created by BSWITCH. + +"Starting record number for spectra to be created:" All created spectra +will have a suffix number starting with this value and incremented +by one for each new spectrum created. + +"File name to contain statistics information:" This file will contain +informative output from SENSFUNC and BSWITCH. (default=stats) + +"File name to contain a log of terminal output:" All tasks talk back +to let you know how things are proceding. The backtalk is saved +in this file. (default=ttylog) + +"File name for output from STANDARD and input to SENSFUNC:" Just +what it says. (default=std) + +"Record string to process:" The spectra are assumed to be representable +by strings (try "help ranges" for details on the formats allowed). +Both STANDARD and BSWITCH expect ranges of spectral record numbers +which are appended to the root given in answer to the first question +above. This question is asked repeatedly so that you can enter as +many strings of spectra as you like and is ended by hitting return +without entering a value. There is a short delay after entering +each string of records while a check is made to verify that all +your spectra actually exist. + +"Standard star name:" For each record string STANDARD expects +the name of the standard star observed, but it must be given in +a manner acceptable to STANDARD. (see STANDARD and LCALIB for +more details). + +"Use weighted averages:" If answered yes, then SENSFUNC and BSWITCH +will use their weighted averaging schemes. + +"Apply magnitude fudging:" If answered yes, then SENSFUNC will +use its "fudge" option. (see SENSFUNC) + +"Solve for grey additive extinction constant:" If answered yes, then +SENSFUNC will solve for this value. + +"File name for sensitivity image file:" This will be the root name +for the output sensitivity spectra from SENSFUNC. + +At anytime during the processing phase, you can inquire about the +progress by listing the latest contents of the file "ttylog" +either by "type ttylog" or by "tail ttylog". The latter command +lists the last 12 lines of the file. + +Be sure to have all your record strings, standard star names, +and options well planned and written down so that you can enter +the answers correctly. The batch reductions are not overly +tolerant of incorrect entries although some preliminary checks +are performed during the entry process. + +.ih +EXAMPLES + +The following invokes the batch reductions using all task options; + + cl> batchred + +The following inhibits the STANDARD and SENSFUNC tasks which must have +been run previously. This is equivalent to the IPPS "autoreduce": + + cl> batchred standard- sensfunc- +.ih +BUGS +If you make an error while entering the requested information, there +is no way to effect repairs other than to (1) start all over, or (2) edit +the generated script file "process.cl" using the system editor. + +If a task encounters an irrecoverable error, the background job +hangs until you kill it using "kill N" where N is the job number. +.ih +SEE ALSO +mkscript, standard, sensfunc, bswitch, calibrate, addsets +.endhelp diff --git a/noao/onedspec/irsiids/doc/bswitch.hlp b/noao/onedspec/irsiids/doc/bswitch.hlp new file mode 100644 index 00000000..a50647b4 --- /dev/null +++ b/noao/onedspec/irsiids/doc/bswitch.hlp @@ -0,0 +1,228 @@ +.help bswitch Sep87 noao.imred.iids/noao.imred.irs +.ih +NAME +bswitch - generate sky-subtracted accumulated spectra +.ih +USAGE +bswitch input records +.ih +PARAMETERS +.ls input +The root name for the input spectra to be beam-switched. +.le +.ls records +The range of spectra to be included in the beam-switch operation. +Each range item will be appended to the root name to form an image +name. For example, if "input" is "nite1" and records is "1011-1018", +then spectra nite1.1011, nite.1012 ... nite1.1018 will be included. +.le +.ls output +New spectra are created by the beam-switch operation. This parameter +specifies the root name to be used for the created spectra. +.le +.ls start_rec = 1 +Each new spectrum created has "output" as its root name and a trailing +number appended. The number begins with start_rec and is incremented +for each new spectrum. For example, if "output" is given as "nite1b" +and start_rec is given as 1001, then new spectra will be created as +nite1b.1001, nite1b.1002 ... +.le +.ls stats = "stats" +A file by this name will have statistical data appended to it, or created +if necessary. If a null file name is given (""), no statistical output +is given. For each aperture, a listing of countrates for each +observation is given relative to the observation with the highest rate. +.le +.ls ids_mode = yes +If the data are taken under the usual IIDS "beam-switch" mode, this +parameter should be set to yes so that accumulations will be performed +in pairs. But if the data are taken where there is no sky observation +or different numbers of sky observations, ids_mode should be set to no. +If weighting is in effect, ids_mode=yes implies weighting of the +object-sky sum; if ids_mode=no, then weighting is applied to the +object and sky independently because then there is no guarantee that +an object and sky observation are related. +.le +.ls extinct = yes +If set to yes, a correction for atmospheric extinction is applied. +The image header must have either a valid entry for AIRMASS or +for hour angle (or right ascension and sidereal time) and declination. +.le +.ls weighting = no +If set to yes, the entire spectrum or a specified region will be used +to obtain a countrate indicative of the statistical weight to be +applied to the spectrum during the accumulations. +.le +.ls subset = 32767 +A subset value larger than the number of independent spectra to be +added indicates that the operation is to produce a single spectrum +for each aperture regardless of how many input spectra are entered. +If subset is a smaller number, say 4, then the accumulations +are written out after every 4 spectra and then re-initialized to zero +for the next 4. +.le +.ls wave1 = 0.0 +If weighting=yes, this parameter indicates the starting point in the +spectrum for the countrate to be assessed. For emission-line objects, +this is particularly useful because the regime of information is then +confined to a narrow spectral region rather than the entire spectrum. +Defaults to the beginning of the spectrum. +.le +.ls wave2 = 0.0 +This provides the ending wavelength for the countrate determination. +Defaults to the endpoint of the spectrum. +.le +.ls observatory = "observatory" +Observatory at which the spectra were obtained if +not specified in the image header by the keyword OBSERVAT. The +observatory may be one of the observatories in the observatory +database, "observatory" to select the observatory defined by the +environment variable "observatory" or the task \fBobservatory\fR, or +"obspars" to select the current parameters set in the \fBobservatory\fR +task. See help for \fBobservatory\fR for additional information. +.le +.ls extinction = ")_.extinction" +The the name of the file containing extinction values. +Required if extinct=yes. +.le +.ih +DESCRIPTION +Data from multiaperture spectrographs are summed according to +aperture number and sky subtracted if sky observations are available. +Data for up to 50 apertures may be simultaneously accumulated. +The accumulated spectra are written to new images. + +The exposure times for each observation may be different. All +internal computations are performed in terms of count rates, +and converted back to counts (for statistical analysis) prior to writing +the new image. Therefore, the time on the sky and object may +be different as well. When these extensions to the normal +mode are required, the flag ids_mode must be set to no. +Then object and sky accumulations are performed totally +independently and a difference is derived at the conclusion +of the operation. + +If ids_mode is set to yes, then the usual IIDS/IRS "beam-switch" +observing mode is assumed. This implies that an equal number of +sky and object spectra are obtained through each aperture +after 2N spectra have been accumulated, where N is the number +of instrument apertures (2 for the IIDS/IRS). It is also assumed +that the object and sky exposure times are equal for each aperture. +Note that the "nebular" mode (where all instrument apertures +point at an extended object simultaneously, and then all apertures +point at sky simultaneously) is an acceptable form for +beam-switched data in ids_mode. + +The accumulations are optionally weighted by the countrate +over a region of the spectrum to improve the statistics during +variable conditions. The user may specify the region of spectrum +by wavelength. In ids_mode, the statistics are obtained from +object-sky differences; otherwise, the statistics are performed +on object+sky and sky spectra separately. + +The spectra may be extinction corrected if this has not already +been performed. +In order to perform either the extinction correction or the +weighting process, the spectra must have been placed on a linear +wavelength scale (or linear in the base 10 logarithm). + +Strings of spectra are accumulated to produce a single +summed spectrum for each observing aperture. But in some cases +it is desirable to produce summed spectra from subsets of the +entire string to evaluate the presence of variations either due +to observing conditions or due to the physical nature of the +object. A subset parameter may be set to the frequency at which +spectra are to be summed. + +In order that the processing occur with minimal user interaction, +elements from the extended image header are used to direct the +flow of operation and to obtain key observing parameters. +The required parameters are: object/sky flag (OFLAG=1/0), exposure +time in seconds (ITM), beam (that is, aperture) number (BEAM-NUM), airmass (AIRMASS) +or alternatively hour angle (HA) and declination (DEC), or +right ascension (RA), sidereal time (ST), declination (DEC), and the +observatory (OBSERVAT), +starting wavelength (W0), and wavelength increment per channel (WPC), +where the names in parenthesis are the expected keywords in the +header. If the observatory is not specified in the image the +observatory parameter is used. See \fBobservatory\fR for further +details on the observatory database. + +The following header flags are used as well: DC_FLAG +for dispersion corrected data (must=0), BS_FLAG for beam-switching +(must not be 1 which indicates the operation was already done), +EX_FLAG for extinction correction (if = 0 extinction is assumed already +done). + +The headers may be listed with the IMHEADER task, setting +the parameter "long" = yes. The values for the parameters follow +the rules used for IIDS and IRS data. + +After the beam-switch operation, the newly created spectra will +have header elements taken from the last object spectrum. +A few parameters will be updated to reflect the operation +(e.g. integration time, processing flags). + +.ih +EXAMPLES +The following example will accumulate a series of 16 spectra obtained +in the normal beam-switched mode and create two new extinction corrected +spectra having names nite1bs.1 and nite1bs.2: + + cl> bswitch nite1 1011-1026 nite1bs 1 + +The following example performs the same functions but accumulates the data +to produce 8 new spectra representing the individual object-sky pairs: + + cl> bswitch nite1 1011-1026 nite1bs 1 subset=4 + +The following example produces an extinction corrected spectrum for every +input spectrum. Note that ids_mode is set to off to generate separate object and +sky sums, and subset is set to 2 so that every pair of spectra (one object and +one sky) are written out as two new spectra: + + cl> bswitch nite1 1011-1026 nite1bs 1 subset=2 ids_mode- + +The next example produces a pair of spectra for each of 3 independent +objects observed, provided that each was observed for the same number +of observations (16 in this case). + +.nf + cl> bswitch nite1 1011-1026,1051-1066,1081-1096 nite1bs 1 \ + >>> subset=16 +.fi + +The next example shows how to use the weighting parameters where +the indicative flux is derived from the region around the emission-line +of 5007A. + +.nf + cl> bswitch nite1 1011-1026 nite1bs 1 weighting- \ + >>> wave1=4990, wave2=5020 +.fi +.ih +TIME REQUIREMENTS +The principle time expenditure goes toward extinction correcting the +data. For IIDS type spectra (length=1024 pixels), approximately 30 cpu +seconds are required to beam-switch a series of 16 spectra. +.ih +BUGS +The number of apertures is restricted to 50 and must be labeled +between 0 and 49 in the image header (the IIDS uses 0 and 1). + +Until an image header editor is available, BSWITCH +can be applied only to data with properly prepared headers +such as IIDS/IRS data read by RIDSMTN, RIDSFILE and some data via RFITS. + +When used to perform the function of extinction correction only (the +third example above), the statistics file fails to note the output +image name for the sky spectrum. + +The data must be on a linear wavelength scale. +The starting wavelength, W0, and a wavelength +per channel, WPC, are required header information, and the DC_FLAG +must be set to 0. +.ih +SEE ALSO +observatory, sensfunc, imheader, lcalib, ridsmtn, ridsfile, rfits +.endhelp diff --git a/noao/onedspec/irsiids/doc/coefs.hlp b/noao/onedspec/irsiids/doc/coefs.hlp new file mode 100644 index 00000000..777933bc --- /dev/null +++ b/noao/onedspec/irsiids/doc/coefs.hlp @@ -0,0 +1,57 @@ +.help coefs May85 noao.imred.iids/noao.imred.irs +.ih +NAME +coefs -- Extract dispersion coefs from mtn HeNeAr headers +.ih +USAGE +coefs input records database +.ih +PARAMETERS +.ls input +The input image root name for the spectral images containing the +dispersion coefficients. +.le +.ls records +The range of records for which the root name applies. +.le +.ls database +The database file name which will contain the coefficients. +.le +.ih +DESCRIPTION +The spectra specified by the combination of the root name +and the records are scanned for the presence of dispersion +coefficients. If present, the coefficients and necessary +information are written to the file indicated by the database +parameter. This file an then be used by the linearization +program DISPCOR to correct any spectra for which the +database is appropriate. + +Each invocation of COEFS appends to the database file, or +creates a new file if necessary. + +The following assumptions are made concerning the coefficients, +which are always correct for IIDS and IRS mountain reduced +data at Kitt Peak. +.ls 5 (1) +The coefficients represent Legendre polynomials. +.le +.ls (2) +The coefficients apply to pixels 1 through 1024 in the original data. +.le +.ih +EXAMPLES +The following example reads the coefficients from the headers +for nite1 arc spectra taken near the beginning and end of the +night and creates a database file called nite1.db: + + cl> coefs nite1 3-4,201-202 nite1.db + +.ih +TIME REQUIREMENTS +Approximately 1 second per spectrum is required. This is primarily +overhead due to file access. +.ih +SEE ALSO +dispcor, identify +.endhelp diff --git a/noao/onedspec/irsiids/doc/coincor.hlp b/noao/onedspec/irsiids/doc/coincor.hlp new file mode 100644 index 00000000..74e002f3 --- /dev/null +++ b/noao/onedspec/irsiids/doc/coincor.hlp @@ -0,0 +1,101 @@ +.help coincor Feb87 noao.imred.iids/noao.imred.irs +.ih +NAME +coincor -- Correct detector count rates +.ih +USAGE +coincor input records +.ih +PARAMETERS +.ls input +The root file name of the input spectra. +.le +.ls records +The range of spectra. +The names of the spectra will be formed by appending the range +elements to the input root name. +.le +.ls output +This is the root file name for the corrected spectra. If no root name +is specified (specified with the null string "") then the operation +is done in place. +.le +.ls start_rec = 1 +The starting record number to be appended to the root name of the +created spectra. +.le +.ls ccmode = )_.ccmode +The mode used to model the detector count rate corrections. +In the following C(obs) is the observed count rate and C(cor) is the +corrected count rate. +.ls "photo" +Photoelectric photometer with discriminator mode. The count rate +correction is + + C(cor) = C(obs) * exp (C(obs) * deadtime) + +where the parameter \fIdeadtime\fR is the representative deadtime in seconds. +.le +.ls "iids" +IIDS correction given by + + C(cor) = (-ln(1-C(obs)*deadtime)/deadtime)**power + +where \fBdeadtime\fR is a parameter related to the sweep time used to +correct for coincidence losses and \fBpower\fR is a power law coefficient. +.le +.le +.ls deadtime = )_.deadtime +For the "photo" mode this parameter is the period, in seconds, during +which no counts can be registered by the detector. Note that this is +based on a per pixel basis. So if the discriminator dead period is of +order 50 nanoseconds and 2000 pixels are observed per readout, the +effective deadtime is about 10E-4 seconds. For the "iids" mode this +parameter defines the sweep time correction and has a value of 1.424E-3 +seconds. +.le +.ls power = )_.power +The IIDS power law coefficient. The standard value is 0.975. +.le +.ih +DESCRIPTION +The input spectra are corrected for detector count rate errors. If no +output root name is given then the operation is done in place. The type +of correction is specified by the parameter \fIccmode\fR. The available +modes are for a general photomultiplier with discriminator coincidence +correction, and the NOAO IIDS. The parameters for these modes are +\fIdeadtime\fR and \fIpower\fR. The exposure time, in seconds, is a +required image header parameter (keyword = EXPOSURE). + +The default mode is for the NOAO IIDS. The IIDS correction includes a +power law correction for a nonlinear effect in the IIDS image tube chain +which is not included by the mountain reduction software at the telescope. +If the spectra have been coincidence corrected at the telescope +then only the nonlinear power law correction is applied. + +The coincidence correction flag may take the values -1 for no correction, +0 for the IIDS correction with \fIpower\fR = 1 (the correction +applied by the mountain reduction software), 1 for the full IIDS +correction, and 2 for the photomuliplier mode correction. +.ih +EXAMPLES +The following example corrects a series of IIDS spectra: + + cl> coincor nite1 1-250 output=nite1cc start_rec=1 + +The following example corrects a series of spectra from the +Lick ITS: + +.nf + cl> coincor its 1-250 output=itscc start=1 ccmode=photo \ + >>> deadtime=2.4E-4 power=1 +.fi +.ih +TIME REQUIREMENTS +\fBCoincor\fR requires approximately 1 second per spectrum of length 1024. +.ih +SEE ALSO +.nf +The \fBimred.iids\fR package is designed for reducing NOAO IIDS spectra. +.fi +.endhelp diff --git a/noao/onedspec/irsiids/doc/extinct.hlp b/noao/onedspec/irsiids/doc/extinct.hlp new file mode 100644 index 00000000..66aca3d6 --- /dev/null +++ b/noao/onedspec/irsiids/doc/extinct.hlp @@ -0,0 +1,49 @@ +.help extinct Apr85 noao.onedspec +.ih +NAME +extinct -- Correct spectra for atmospheric extinction +.ih +USAGE +extinct root records output +.ih +PARAMETERS +.ls root +The root name for the input spectra to be corrected. +.le +.ls records +The range of spectra to be included in the extinction operation. +.le +.ls output +The root name for the output corrected spectra +.le +.ls start_rec +The starting record number for the output corrected spectra. +.le +.ls nr_aps = 2 +The number of instrument apertures for this data set. +.le +.ih +DESCRIPTION +The input spectra are corrected for atmospheric extinction. +EXTINCT redirects the spectra through the task BSWITCH so all +procedures are identical to those described for that task. + +Because BSWITCH attempts to perform a beam-switch operation +unless the subset parameter is equal to the number of +instrument apertures (in which case beam-switching degenerates +to a copy operation), the hidden parameter nr_aps should be set +appropriately to the instrument. For IIDS and IRS data, this +is 2. +.ih +EXAMPLES + + cl> extinct nite1 1001-1032 nite1ex +.ih +BUGS +The input string of spectra must be ordered so that only +one spectrum from each aperture is present among substrings +of length nr_aps. +.ih +SEE ALSO +bswitch +.endhelp diff --git a/noao/onedspec/irsiids/doc/flatdiv.hlp b/noao/onedspec/irsiids/doc/flatdiv.hlp new file mode 100644 index 00000000..e6e8c22e --- /dev/null +++ b/noao/onedspec/irsiids/doc/flatdiv.hlp @@ -0,0 +1,94 @@ +.help flatdiv Dec86 noao.imred.iids/noao.imred.irs +.ih +NAME +flatdiv -- Divide spectra by flat field spectra +.ih +USAGE +flatdiv input records +.ih +PARAMETERS +.ls input +The root file name for the input records to be divided. +.le +.ls records +The range of spectra to be included in the divide operation. +Each range item will be appended to the root name to form an +image file name. +.le +.ls output +New spectra are created by the flatdiv operation. This parameter +specifies the root name to be used for the created spectra. +.le +.ls start_rec +Each new spectrum created as "output" has its root name and a +trailing number appended starting with "start_rec". Subsequent +output images will have an incremented trailing number. +Note that even if an output image is not created because the input +image has already been flattened or the input image is not found the +output record number is still incremented. +.le +.ls flat_file +The root name for the sensitivity spectra as produced by FLATFIT. +Normally with multi-aperture instruments, FLATFIT will produce a +spectrum appropriate to each aperture and the file name will have +"flat_file" as the file name root and the aperture number appended. +.le +.ls coincor = )_.coincor +If set to yes, coincidence correction is applied to the data during +the division, and the following three parameters are required. +For more about this correction see \fBcoincor\fR. +.ls ccmode = )_.ccmode +The mode by which the coincidence correction is to be performed. +This may be "iids" or "photo". +.le +.ls deadtime = )_.deadtime +The detector deadtime in seconds. +.le +.ls power = )_.power +Power law IIDS non-linear correction exponent. +.le +.le +.ih +DESCRIPTION +The input spectra are divided by the flat fields which are +represented by spectra produced by FLATFIT. + +To avoid possible division by zero, any zeroes in the flat field +spectra generated by FLATFIT are replaced by 1.0. + +The input spectra may optionally be corrected for coincidence losses. + +If the input and output spectra (after appending the record numbers) are +the same then the division is performed in-place; i.e. the flattened spectra +replace the original input spectra. +Note that even if an output image is not created because the input +image has already been flattened or the input image is not found the +output record number is still incremented. This is to insure that if +in-place division is desired that the input and output names remain +matched. +.ih +EXAMPLES +The following example divides a series of spectra to produce 20 new +spectra having names nite1.1221 ... nite1.1240. + + cl> flatdiv nite1 1201-1220 nite1 1221 + +The same spectra as above are simultaneously corrected for +coincidence losses. + + cl> flatdiv nite1 1201-1220 nite1 1221 coincor=yes + +The flattened spectra replace the unflattened spectra. + + cl> flatdiv nite1 1201-1220 nite1 1201 + +Note that the input record numbers must be contiguous and the starting +output record number must be the same as the first input record number. +.ih +TIME REQUIREMENTS +Approximately 1 second is required to correct a spectrum of length +1024 points. +.ih +SEE ALSO +coincor, flatfit +.endhelp diff --git a/noao/onedspec/irsiids/doc/flatfit.hlp b/noao/onedspec/irsiids/doc/flatfit.hlp new file mode 100644 index 00000000..af84cb3c --- /dev/null +++ b/noao/onedspec/irsiids/doc/flatfit.hlp @@ -0,0 +1,188 @@ +.help flatfit Dec86 noao.imred.iids/noao.imred.irs +.ih +NAME +flatfit -- Sum and normalize flat field spectra +.ih +USAGE +flatfit root records +.ih +PARAMETERS +.ls root +The root file name for the input names of the flat field +spectra to be accumulated and fit for normalization. +.le +.ls records +The range of spectra indicating the elements of the string. +The names of the spectra will be formed by appending the range +elements to the input root name. +.le +.ls output +This is the root file name for the names of the spectra which will +be created during normalization. The aperture number for the observation +will be appended to the root in form "root.nnnn" where nnnn is the aperture +number with leading 0's. +.le +.ls function = "chebyshev" +The accumulated spectra are fit by this function type - either +chebyshev or legendre polynomials, or spline3 or spline1 interpolators. +.le +.ls order = 4 +The order of the fit using the above function. This should generally be +a low order fit to avoid introduction of high spatial frequency wiggles. +.le +.ls niter = 1 +The number of iterations to reject discrepant pixels upon initial +startup of the solution. +.le +.ls lower = 2.0 +The number of sigmas for which data values less than this cutoff are +rejected. +.le +.ls upper = 2.0 +The number of sigmas for which data values greater than this cutoff are +rejected. +.le +.ls ngrow = 0 +The number of pixels on either side of a rejected pixel to also be rejected. +.le +.ls div_min = 1.0 +During the normalization process, a division by zero will produce +this value as a result. +.le +.ls interact = yes +If set to yes, graphical interaction with the normalization process +is provided for at least the first aperture for which sums are available. +If set to no, no interaction is provided. +.le +.ls all_interact = no +If set to yes, then interaction will be provided for all apertures +for which sums have been accumulated. If set to no then the parameter interact +will determine if the first aperture data is to be interactive. +.le +.ls coincor = )_.coincor +If set to yes, coincidence correction is applied to the data during +the summation process, and the following three parameters are required. +See \fBcoincor\fR for more about this correction. +.ls ccmode = )_.ccmode +The mode by which the coincidence correction is to be performed. +This may be "iids" or "photo". +.le +.ls deadtime = )_.deadtime +The detector deadtime in seconds. +.le +.ls power = )_.power +Power law IIDS non-linear correction exponent. +.le +.le +.ls cursor = "" +Graphics cursor input. When null the standard cursor is used otherwise +the specified file is used. +.le +.ih +DESCRIPTION +The specified spectra are added by aperture number to produce +summations which are then fit by a specified fitting function. +The fitting function is then divided into the sum to produce a +normalized (to 1.0) sum in which the low frequency spatial +response has been removed. + +The resultant normalized images may then be divided into all other +data to remove the pixel-to-pixel variations without introducing +any color terms. The spectra may be used directly if they happen +to be object spectra in which the low frequency response is to be +removed. + +During the accumulation process the spectra may be corrected for +coincidence losses if the detector is subject to the phenomenon. + +After accumulating all input spectra, the pixels in each sum are +fit according to +the specified function. If the interactive switches are set, then +graphical interaction is made available. If only the interact parameter +is set to yes, then only the data from the first aperture will +be available for interaction. Data from subsequent apertures will +be fit using the same parameters and number of iterations as the +first. If the all_interact parameter is also +set, then data from each aperture will be presented for interaction. + +At each step in the fit, pixels which are discrepant by more than +"upper" sigmas above the fit, or "lower" sigmas below the fit, are +rejected. The rejection process may be applied many times (iterations) +to continue rejecting pixels. If the upper and lower sigmas are +not equal, the resulting fit will be biased slightly above the mean +(for lower < upper) or below the mean (upper < lower). This is useful +when the spectrum being fit is that of a star having either absorption +or emission lines. + +A display is presented of the sum and the fit through the data. +A status line is printed containing the fit type, the order of +the fit, the rms residual from the fit, and the number of data +points in the fit after one iteration of rejection. + +The following cursor keystrokes are then active: +.ls ? +Clear the screen and display the active keystrokes +.le +.ls / +Indicate active keystrokes on the status line +.le +.ls e +Change plot mode to an error plot. This display is defined +as the deviation from the fit divided by the data values [ (data - fit)/ data] +at each pixel +.le +.ls f +Change plot mode back to the fit through the data display +.le +.ls o +Change the order of the fit. +.le +.ls l +Change the lower rejection criterion (in units of sigma). +.le +.ls u +Change the upper rejection criterion. +.le +.ls s +Change both rejection criteria to the same value. +.le +.ls r +Reinstate rejected pixels. +.le +.ls i +Iterate one more time. +.le +.ls n +Iterate several more times - the user is prompted for the count. +.le +.ls q +Quit and accept the solution +.le +.ls <CR> +RETURN is the same as 'q' but a confirmation request to exit must be +answered as yes. +.le + +All keystrokes but ?,/,e,f, and q force another iteration which will +reject additional pixels. To fully inhibit pixel rejection, the sigmas +should be set to a large value (e.g. 100). +.ih +EXAMPLES +The following example will accumulate 8 spectra and fit the first +aperture data interactively but not the second, and apply coincidence +corrections to the sums. The upper and lower rejection criteria +have been altered to bias the seventh order fit to a higher level. + + cl> flatfit nite1 1-4,201-204 coin+ low=1.4 up=3 order=7 +.ih +BUGS +For some reason, the error plot is supposed to have a zero level line +drawn, but none appears. + +As in most of the IRAF software, the order of a fit refers to the number +of terms in the fit, so that a fit of order 1 implies a constant and order +2 implies a linear fit. +.ih +SEE ALSO +coincor, flatdiv +.endhelp diff --git a/noao/onedspec/irsiids/doc/powercor.hlp b/noao/onedspec/irsiids/doc/powercor.hlp new file mode 100644 index 00000000..e1f9c70e --- /dev/null +++ b/noao/onedspec/irsiids/doc/powercor.hlp @@ -0,0 +1,62 @@ +.help powercor Oct86 noao.imred.iids/noao.imred.irs +.ih +NAME +powercor -- Apply power law correction to mountain reduced spectra +.ih +USAGE +powercor input records +.ih +PARAMETERS +.ls input +The root file name of the input spectra. +.le +.ls records +The range of spectra. +The names of the spectra will be formed by appending the range +elements to the input root name. +.le +.ls output +This is the root file name for the corrected spectra. +.le +.ls start_rec = 1 +The starting record number to be appended to the root name of the +created spectra. +.le +.ls power = )iids.power +The power law coefficient. +.le +.ih +DESCRIPTION +A power law correction to the IIDS count rates is applied to the input +spectra. The mountain reduction software applies a coincidence correction +to the observed IIDS count rates but does not correct for a nonlinear effect +in the image tube chain. This second correction takes the form of a +power law + + C(out) = C(in) ** power + +where C(in) is the input, coincidence corrected, count rate and C(out) +is the corrected count rate. The power is a parameter of the task +which defaults to the \fBiids\fR package parameter set to the appropriate +value for the IIDS. The exposure time, in seconds, is a required +image header parameter (keyword = EXPOSURE) used to convert the +total counts to count rates. + +Note that if the original raw spectra are being reduced then the either +\fBcoincor\fR or \fBpowercor\fR may be used to apply both the coincidence +correction and the power law correction at the same time. In other words, +the tasks apply the coincidence correction if the coincidence flag (CO-FLAG) is +-1 (uncorrected) and the power law correction alone if the flag is zero +(coincidence corrected only). The flag is 1 when both the coincidence and +nonlinear correction have been applied. + +This task is a script calling \fBcoincor\fR with \fIccmode\fR = "iids". +.ih +EXAMPLES +The following example corrects a series of IIDS spectra: + + cl> powercor nite1 1-250 output=nite1cc start_rec=1 +.ih +SEE ALSO +coincor +.endhelp diff --git a/noao/onedspec/irsiids/doc/process.hlp b/noao/onedspec/irsiids/doc/process.hlp new file mode 100644 index 00000000..5cedcde3 --- /dev/null +++ b/noao/onedspec/irsiids/doc/process.hlp @@ -0,0 +1,20 @@ +.help process Oct85 noao.imred.iids/noao.imred.irs +.ih +NAME +process -- A task generated by BATCHRED +.ih +USAGE +process +.ih +DESCRIPTION +The task \fBbatchred\fR creates a script called process.cl for batch +reductions. \fBBatchred\fR also has an option to automatically run +this script. +.ih +EXAMPLES +The task \fBbatchred\fR is run to setup a set of beam switching operations. +It creates the script \fBprocess.cl\fR which the user runs as a background +process as follows: + + cl> process& +.endhelp diff --git a/noao/onedspec/irsiids/doc/slist1d.hlp b/noao/onedspec/irsiids/doc/slist1d.hlp new file mode 100644 index 00000000..6c7d2702 --- /dev/null +++ b/noao/onedspec/irsiids/doc/slist1d.hlp @@ -0,0 +1,59 @@ +.help slist1d Jan92 noao.imred.irs/iids +.ih +NAME +slist1d -- List spectral header information +.ih +USAGE +slist1d input records +.ih +PARAMETERS +.ls input +The image root name for the spectra to be listed. +.le +.ls records +The record string for the spectra to be listed. The records will be appended +to the root name to form image names of the type "root.xxxx". +.le +.ls long_header = no +If set to yes, then a complete listing of the header elements +is given. If set to no, then a single line per spectrum is given which lists +in the following order: the image name, object or sky spectrum, exposure +time, spectrum length, and image title. +.le +.ih +DESCRIPTION +Each spectrum in the list implied by the root name and the record string +is opened and the header is read. The pixel file is not accessed in order +to save time. The header listing is directed to STDOUT and may be +redirected for printing. + +A warning message is issued if +a requested image is not found, but otherwise proceeds. +.ih +EXAMPLES +The following example lists 8 spectral headers in long form on the printer: + +.nf + cl> slist1d nite1 1001-1008 | lprint +.fi + +The next example lists the same spectral headers but in short form +on the terminal + +.nf + cl> slist1d nite1 1001-1008 long- +.fi +.ih +REVISIONS +.ls SLIST1D V2.10 +This task is the same as V2.9 \fBslist\fR and applies only to the older +IRS/IIDS record extension spectra. In V2.10 \fBslist\fR +has been revised for multiaperture spectra. +.le +.ih +BUGS +SLIST1D does not inform the user if the pixel file can or cannot be read. +.ih +SEE ALSO +slist, imheader +.endhelp diff --git a/noao/onedspec/irsiids/doc/subsets.hlp b/noao/onedspec/irsiids/doc/subsets.hlp new file mode 100644 index 00000000..a9f0ae68 --- /dev/null +++ b/noao/onedspec/irsiids/doc/subsets.hlp @@ -0,0 +1,49 @@ +.help subsets May85 noao.imred.iids/noao.imred.irs +.ih +NAME +subsets - Subtract pairs of spectra in a string +.ih +USAGE +subsets input records +.ih +PARAMETERS +.ls input +The root file name for the input spectra in the string. +.le +.ls records +The range of spectra indicating the elements of the string. +The names of the spectra will be formed by appending the range +elements to the input root name. +.le +.ls output +This is the root file name for the names of the spectra which will +be created by the subtraction operation. +.le +.ls start_rec +The starting record number to be appended to the root name of the +created spectra. +.le +.ih +DESCRIPTION +Pairs of spectra are formed from the input string in the order that +the record numbers would suggest. +The first spectrum in the pair is assumed to be the +principle spectrum and the second spectrum in the pair is subtracted +from the first. The result is written out as a new spectrum. + +No compensation is made for exposure time during the subtraction. +The header from the principle spectrum is assigned to the output +spectrum. + +.ih +EXAMPLES +The following example forms 50 new spectra from nite1.2001-nite1.2002, +nite1.2003-nite1.2004, ... + + cl> subsets nite1 2001-2100 + +The following example creates new spectra from the pairs nite2.2001-nite2.2002, +nite2.2003-nite2.2004 in spite of the order of the record numbers entered. + + cl> subsets nite2 2001,2003,2002,2004 +.endhelp diff --git a/noao/onedspec/irsiids/doc/sums.hlp b/noao/onedspec/irsiids/doc/sums.hlp new file mode 100644 index 00000000..0d8b27e9 --- /dev/null +++ b/noao/onedspec/irsiids/doc/sums.hlp @@ -0,0 +1,44 @@ +.help sums Jul85 noao.imred.iids/noao.imred.irs +.ih +NAME +sums -- Generate sums of the sky and object spectra for each aperture +.ih +USAGE +sums input records +.ih +PARAMETERS +.ls input +The root file name for the input spectra in the string. +.le +.ls records +The range of spectra indicating the elements of the string. +The names of the spectra will be formed by appending the range +elements to the input root name. +.le +.ls output +This is the root file name for the names of the spectra which will +be created by the summation operation. +.le +.ls start_rec +The starting record number to be appended to the root name of the +created spectra. +.le +.ih +DESCRIPTION +All the object spectra for each aperture are summed, and the +sky spectra are also summed to produce two new spectra for +each observing aperture. Exposure times are accumulated. +No tests are made to check whether the object is consistent +among the specified spectra. This could be accomplished by +checking the titles or telescope positions, but it isn't. + +The header parameters OFLAG and BEAM-NUM must be properly +set in the headers. +.ih +EXAMPLES +The following example forms 4 new spectra from nite1.2001-nite1.2002, +nite1.2003-nite1.2004, ... assuming this string is derived from +IIDS spectra. + + cl> sums nite1 2001-2100 +.endhelp diff --git a/noao/onedspec/irsiids/doc/widstape.hlp b/noao/onedspec/irsiids/doc/widstape.hlp new file mode 100644 index 00000000..855f223d --- /dev/null +++ b/noao/onedspec/irsiids/doc/widstape.hlp @@ -0,0 +1,90 @@ +.help widstape Mar85 noao.imred.iids/noao.imred.irs +.ih +NAME +widstape -- Write a Cyber style IDSOUT tape +.ih +USAGE +widstape idsout input records +.ih +PARAMETERS +.ls idsout +The output file name to receive the card-image data. This may be a +magtape specification (e.g. mta, mtb) or disk file name. +.le +.ls input +The input root file name for the spectra to be written +.le +.ls records +The record string to be appended to the root name to create the image +names of the spectra to be written. +.le +.ls new_tape = no +If set to yes, the tape is rewound and output begins at BOT. If no, +output begins at EOT unless an explicit file specification is given +as part of the magtape file name for parameter "idsout" (e.g. mta[2]). +If idsout contains a file specification of [1], then writing begins +at BOT regardless of the value for new_tape. +.le +.ls block_size = 3200 +The tape block size in bytes. This must be an integral factor of 80. +.le +.ls ebcdic = no +The default character code is ASCII, but if this parameter is set to yes, +the output character will be in EBCDIC. +.le +.ih +DESCRIPTION +The specified spectra are copied to the output file in a card-image format +defined in the IPPS-IIDS/IRS Reduction Manual. Values from the extended +image header are used to fill in the observational parameters. + +The basic format consists of 4 - 80 byte header cards, 128 data cards +having 8 data elements per card in 1PE10.3 FORTRAN equivalent format, +and a trailing blank card for a total of 133 cards. +Thus spectra up to 1024 points may be contained in the IDSOUT format. +The format is outlined below: + +.nf + Line Column Type + 1 1-5 Integer Record number within IDSOUT text file + 6-10 Integer Integration time + 11-25 Real Wavelength of first bin + 26-40 Real Dispersion + 41-45 Integer 0 (Index of first pixel) + 46-50 Integer Line length - 1 (Index of last pixel) + 71-80 Integer UT time + 2 1-10 Real Siderial time + 11-25 Real Right Ascension + 26-40 Real Declination + 3 21-35 Real Hour Angle + 36-50 Real Air mass + 51-58 Integer UT date + 60-76 String Image title + 78-80 String END + 4 1-64 String Record label + 78-80 String END +5-132 Real 1024 pixel values, 8 per line + 133 Blank line +.fi + +The data of type real are in exponent format; i.e FORTRAN 'E' format (1.234e3). + +There are no special marks between spectral images, +and when multiple spectra are written with a single command, the first card +of a subsequent spectrum may be within the same physical tape block +as the last card of the previous spectrum. This assures that all tape +blocks (except the very last one in the tape file) are all the same +length. A double end-of-mark is written after the last spectrum. +.ih +EXAMPLES +The following example writes an IDSOUT format tape starting at the +beginning of the tape. + + cl> widstape mta nite1 1001-1200 new_tape+ +.ih +TIME REQUIREMENTS: UNIX/VAX 11/750 +Each spectrum of 1024 points requires about 2 second. +.ih +SEE ALSO +rcardimage, ridsout +.endhelp diff --git a/noao/onedspec/irsiids/extinct.cl b/noao/onedspec/irsiids/extinct.cl new file mode 100644 index 00000000..68c5a2de --- /dev/null +++ b/noao/onedspec/irsiids/extinct.cl @@ -0,0 +1,22 @@ +#{ EXTINCT -- Use the BSWITCH task to perform the correction for +# atmospheric extinction. + +{ +# Root name +rt = root + +# Records +rec = records + +# Output root +out = output + +# Output starting record +strt = start_rec + +# Do operation +# Inhibit weighting and statisitic file generation + +bswitch (input=rt, records=rec, output=out, start_rec=strt, subset=nr_aps, + weighting=no, ids_mode=no, stats="") +} diff --git a/noao/onedspec/irsiids/extinct.par b/noao/onedspec/irsiids/extinct.par new file mode 100644 index 00000000..a605dae0 --- /dev/null +++ b/noao/onedspec/irsiids/extinct.par @@ -0,0 +1,11 @@ +# EXINCT + +root,s,a,,,,Root name for spectra file names +records,s,a,,,,Record string to process +output,s,a,,,,Root name for spectra to be created +start_rec,i,a,1,0,9999,Next starting spectral record +nr_aps,i,h,2 +strt,i,h +rt,s,h +out,s,h +rec,s,h diff --git a/noao/onedspec/irsiids/flatdiv.par b/noao/onedspec/irsiids/flatdiv.par new file mode 100644 index 00000000..84de42d4 --- /dev/null +++ b/noao/onedspec/irsiids/flatdiv.par @@ -0,0 +1,12 @@ + +# FLATDIV parameter file + +input,s,a,,,,Input image root file name +records,s,a,,,,Range of spectral records +output,s,a,,,,Output file root name for new spectra +start_rec,i,a,1,0,9999,Next starting spectral record +flat_file,s,a,,,,Image root name for output flat field spectra +coincor,b,h,)_.coincor,,,Apply coincidence correction to spectra +ccmode,s,h,)_.ccmode,,,Correction mode (photo|iids) +deadtime,r,h,)_.deadtime,,,Deadtime in seconds +power,r,h,)_.power,,,IIDS power law coefficient diff --git a/noao/onedspec/irsiids/flatfit.par b/noao/onedspec/irsiids/flatfit.par new file mode 100644 index 00000000..3167697d --- /dev/null +++ b/noao/onedspec/irsiids/flatfit.par @@ -0,0 +1,24 @@ +# FLATFIT parameter file + +input,s,a,,,,Input image root file name +records,s,a,,,,Range of spectral records +output,s,a,,,,Output file root name for new spectra +function,s,h,"chebyshev",,,Function to fit (chebyshev|legendre|spline3|spline1) +order,i,h,6,1,,Fitting order (number of terms) +niter,i,h,1,1,,Number of rejection iterations +lower,r,h,2,0,,Lower rejection criterion in sigmas +upper,r,h,2,0,,Upper rejection criterion in sigmas +ngrow,i,h,0,,,Growing region +div_min,r,h,1.0,,,Value to use if division by zero occurs +interact,b,h,yes,,,Interact with the first accumulation? +all_interact,b,h,no,,,Interact with all accumulations? +coincor,b,h,)_.coincor,,,Apply coincidence correction to flats +ccmode,s,h,)_.ccmode,,,Correction mode (photo|iids) +deadtime,r,h,)_.deadtime,,,Deadtime in seconds +power,r,h,)_.power,,,IIDS power law coefficient +new_order,i,a,4,1,,enter order +new_lower,r,a,,,,enter nr sigma +new_upper,r,a,,,,enter nr sigma +new_niter,i,a,,,,enter nr of iterations +confirm,b,a,,,,Exit and save solution? +cursor,*gcur,h,"",,,Graphics cursor input diff --git a/noao/onedspec/irsiids/getnimage.x b/noao/onedspec/irsiids/getnimage.x new file mode 100644 index 00000000..c232a85c --- /dev/null +++ b/noao/onedspec/irsiids/getnimage.x @@ -0,0 +1,133 @@ +include <mach.h> + + +# GET_NEXT_IMAGE -- Use root filename and ranges string (if any) to +# generate the next image filename. Return EOF +# when image list is exhausted. + +int procedure get_next_image (infile, records, nrecs, image, sz_name) + +int infile, records[ARB], nrecs, sz_name +char image[sz_name] + +int next_num, stat +int flag1, flag2, flag3 +char image_0[SZ_FNAME] + +int clgfil(), get_next_entry(), strlen() + +common /gnicom/ flag1, flag2 + +data flag3/YES/ + +begin + # Reset initializer, record counter, and get root name + if ((flag1 == YES) || (flag3 == YES)) { + next_num = -1 + call rst_get_entry () + } + + # If no ranges specified, act like template expander + if (nrecs == MAX_INT) { + stat = clgfil (infile, image, sz_name) + + # Otherwise append record numbers to first template expansion + } else { + if (flag1 == YES) { + stat = clgfil (infile, image_0, sz_name) + if (stat == EOF) + return (stat) + } + + stat = get_next_entry (records, next_num) + if (stat != EOF) { + call strcpy (image_0, image, sz_name) + call sprintf (image[strlen(image)+1], sz_name, ".%04d") + call pargi (next_num) + } + } + + flag1 = NO + flag3 = NO + return (stat) +end + + +# Reset the initialization parameter to TRUE + +procedure reset_next_image () + +int flag1, flag2 +common /gnicom/ flag1, flag2 + +begin + flag1 = YES +end + + +# GET_NEXT_ENTRY -- Given a list of ranges and the current file number, +# find and return the next file number in order of entry. +# EOF is returned at the end of the list. + +int procedure get_next_entry (ranges, number) + +int ranges[ARB] # Range array +int number # Both input and output parameter + +int ip, first, last, step, next_number, remainder +int flag1, flag2, flag3 + +common /gnicom/ flag1, flag2 + +data flag3/YES/ + +begin + number = number + 1 + next_number = MAX_INT + if ((flag2 == YES) || (flag3 == YES)) { + ip = 1 + flag2 = NO + flag3 = NO + } + + first = min (ranges[ip], ranges[ip+1]) + last = max (ranges[ip], ranges[ip+1]) + step = ranges[ip+2] + + if (number >= first && number <= last) { + remainder = mod (number - first, step) + if (remainder == 0) + return (number) + if (number - remainder + step <= last) + next_number = number - remainder + step + else + go to 10 + + } else if (first > number) + next_number = min (next_number, first) + + else { +10 ip = ip + 3 + if (ranges[ip] != -1 && ranges[ip+1] !=0 && ranges[ip+2] !=0) + next_number = min (ranges[ip], ranges[ip+1]) + } + + if (next_number == MAX_INT) { + ip = 1 + flag2 = YES + return (EOF) + + } else { + number = next_number + return (number) + } +end + +procedure rst_get_entry () + +int first, flag2 +common /gnicom/ first, flag2 + +begin + flag2 = YES +end diff --git a/noao/onedspec/irsiids/idsmtn.h b/noao/onedspec/irsiids/idsmtn.h new file mode 100644 index 00000000..5fa3c73e --- /dev/null +++ b/noao/onedspec/irsiids/idsmtn.h @@ -0,0 +1,101 @@ +# Definitions for the Mountain format IDS tape reader: + +define MAX_RANGES 100 +define DUMMY 3 # Value returned if DUMMY IDS record is read + +define NBITS_CHAR (SZB_CHAR * NBITS_BYTE) +define SZ_IDS_RECORD (2108 * 16 / NBITS_CHAR) +define NPIX_IDS_REC 1024 +define LEN_USER_AREA 2880 +define SZ_IDS_ID 64 +define NBITS_VN_3WRD_FP 48 +define NBITS_VN_2WRD_FP 32 +define NBITS_VN_LONG_INT 32 +define NBITS_VN_INT 16 +define NBITS_FORTH_CHAR 8 +define DATA_BYTE 9 # First byte of data +define NBYTES_DATA (1024 * 32 / NBITS_BYTE) # Number of data bytes +define NBYTES_INT (NBITS_INT / NBITS_BYTE) +define NBYTES_VN_3WRD_FP 6 +define NBYTES_VN_2WRD_FP 4 +define NBITS_2WRD_HIGH 8 +define WRD2_EXP_OFFSET 25 +define NBITS_2WRD_EXP 6 +define WRD2_MANT_SIGN 24 +define WRD2_EXP_SIGN 31 +define NSIG_VN_BITS 15 +define VN_LONG_SIGN 31 +define WRD3_MANT_SIGN 31 +define WRD3_EXP_SIGN 47 +define MAX_NCOEFF 25 + + +# The control parameter structure is defined below: + +define LEN_CP 10 + SZ_FNAME + 1 + +define IS_REDUCED Memi[$1] +define LONG_HEADER Memi[$1+1] +define PRINT_PIXELS Memi[$1+2] +define MAKE_IMAGE Memi[$1+3] +define OFFSET Memi[$1+4] +define DATA_TYPE Memi[$1+5] +define IRAF_FILE Memc[P2C($1+10)] + + +# The header structure is defined below: + +define LEN_IDS 40 + SZ_IDS_ID + 1 + +define HA Memr[P2R($1)] +define AIRMASS Memr[P2R($1+1)] +define RA Memr[P2R($1+2)] +define DEC Memr[P2R($1+3)] +define W0 Memr[P2R($1+4)] +define WPC Memr[P2R($1+5)] +define LINE Memi[$1+6] +define NP1 Memi[$1+7] +define NP2 Memi[$1+8] +define ITM Memr[P2R($1+9)] +define BEAM Memi[$1+10] +define W Memi[$1+11] +define UT Memr[P2R($1+13)] +define ST Memr[P2R($1+14)] +define DF_FLAG Memi[$1+15] +define SM_FLAG Memi[$1+16] +define QF_FLAG Memi[$1+17] +define DC_FLAG Memi[$1+18] +define QD_FLAG Memi[$1+19] +define EX_FLAG Memi[$1+20] +define BS_FLAG Memi[$1+21] +define CA_FLAG Memi[$1+22] +define CO_FLAG Memi[$1+23] +define OFLAG Memi[$1+24] +define POINT Memi[$1+25] +define DRA Memi[$1+26] +define DDEC Memi[$1+27] +define ALPHA_ID Memc[P2C($1+35)] +define LABEL Memc[P2C($1+40)] + + +# Bit offsets to various IDS header words are defined below: + +define NREC_OFFSET ((0 * 16) + 1) +define RFLAGS_OFFSET ((1 * 16) + 1) +define ITM_OFFSET ((2 * 16) + 1) +define DATA_OFFSET ((4 * 16) + 1) +define W0_OFFSET ((2052 * 16) + 1) +define WPC_OFFSET ((2055 * 16) + 1) +define NP1_OFFSET ((2058 * 16) + 1) +define NP2_OFFSET ((2059 * 16) + 1) +define OFLAG_OFFSET ((2060 * 16) + 1) +define SMODE_OFFSET ((2061 * 16) + 1) +define UT_OFFSET ((2062 * 16) + 1) +define ST_OFFSET ((2064 * 16) + 1) +define BEAM_OFFSET ((2066 * 16) + 1) +define HA_OFFSET ((2067 * 16) + 1) +define RA_OFFSET ((2070 * 16) + 1) +define DEC_OFFSET ((2073 * 16) + 1) +define DRA_OFFSET ((2076 * 16) + 1) +define DDEC_OFFSET ((2077 * 16) + 1) +define LABEL_OFFSET ((2078 * 16) + 1) diff --git a/noao/onedspec/irsiids/irsiids.hd b/noao/onedspec/irsiids/irsiids.hd new file mode 100644 index 00000000..d0a20d98 --- /dev/null +++ b/noao/onedspec/irsiids/irsiids.hd @@ -0,0 +1,18 @@ +# Help directory for the IRS/IIDS tasks. + +$doc = "./doc/" + +addsets hlp=doc$addsets.hlp, src=t_addsets.x +batchred hlp=doc$batchred.hlp, src=batchred.cl +bswitch hlp=doc$bswitch.hlp, src=t_bswitch.x +coefs hlp=doc$coefs.hlp, src=t_coefs.x +coincor hlp=doc$coincor.hlp, src=t_coincor.x +extinct hlp=doc$extinct.hlp, src=extinct.cl +flatdiv hlp=doc$flatdiv.hlp, src=t_flatdiv.x +flatfit hlp=doc$flatfit.hlp, src=t_flatfit.x +powercor hlp=doc$powercor.hlp, src=powercor.cl +process hlp=doc$process.hlp, src=process.cl +slist1d hlp=doc$slist1d.hlp, src=t_slist1d.x +subsets hlp=doc$subsets.hlp, src=t_subsets.x +sums hlp=doc$sums.hlp, src=t_sums.x +widstape hlp=doc$widstape.hlp, src=x_widstape.x diff --git a/noao/onedspec/irsiids/mkpkg b/noao/onedspec/irsiids/mkpkg new file mode 100644 index 00000000..01ee3403 --- /dev/null +++ b/noao/onedspec/irsiids/mkpkg @@ -0,0 +1,22 @@ +# IRS/IIDS Tasks + +$checkout libpkg.a .. +$update libpkg.a +$checkin libpkg.a .. +$exit + +libpkg.a: + coincor.x + conversion.x + getnimage.x <mach.h> + t_addsets.x <error.h> <imhdr.h> + t_bswitch.x <smw.h> <error.h> <imhdr.h> <mach.h> <time.h> + t_coefs.x <error.h> + t_coincor.x <error.h> <imhdr.h> + t_flatdiv.x <error.h> <imhdr.h> + t_flatfit.x <gset.h> <imhdr.h> <math/curfit.h> + t_slist1d.x <error.h> <fset.h> <imhdr.h> <smw.h> + t_subsets.x <error.h> <imhdr.h> + t_sums.x <error.h> <imhdr.h> + t_widstape.x <error.h> <imhdr.h> <mach.h> <smw.h> + ; diff --git a/noao/onedspec/irsiids/powercor.cl b/noao/onedspec/irsiids/powercor.cl new file mode 100644 index 00000000..a89e2478 --- /dev/null +++ b/noao/onedspec/irsiids/powercor.cl @@ -0,0 +1,4 @@ +#{ Apply nonlinear IIDS correction + +coincor (input, records, output, start_rec=start_rec, ccmode="iids", + power=power) diff --git a/noao/onedspec/irsiids/powercor.par b/noao/onedspec/irsiids/powercor.par new file mode 100644 index 00000000..e4e8bb90 --- /dev/null +++ b/noao/onedspec/irsiids/powercor.par @@ -0,0 +1,7 @@ +# POWERCOR parameter file + +input,s,a,,,,Input image root file name +records,s,a,,,,Range of spectral records +output,s,a,,,,Output file root name for new spectra +start_rec,i,a,1,0,9999,Next starting spectral record +power,r,h,)_.power,,,Power law coefficient diff --git a/noao/onedspec/irsiids/slist1d.par b/noao/onedspec/irsiids/slist1d.par new file mode 100644 index 00000000..ae091a20 --- /dev/null +++ b/noao/onedspec/irsiids/slist1d.par @@ -0,0 +1,3 @@ +input,s,a,,,,Input image root file name +records,s,a,,,,Range of spectral records +long_header,b,h,no,,,List header in long format diff --git a/noao/onedspec/irsiids/subsets.par b/noao/onedspec/irsiids/subsets.par new file mode 100644 index 00000000..5368ca06 --- /dev/null +++ b/noao/onedspec/irsiids/subsets.par @@ -0,0 +1,6 @@ +# SUBSETS parameter file + +input,s,a,,,,Input image root file name +records,s,a,,,,Range of spectral records +output,s,a,,,,Output file root name for new spectra +start_rec,i,a,1,0,9999,Next starting spectral record diff --git a/noao/onedspec/irsiids/sums.par b/noao/onedspec/irsiids/sums.par new file mode 100644 index 00000000..bbeab466 --- /dev/null +++ b/noao/onedspec/irsiids/sums.par @@ -0,0 +1,8 @@ + +# SUMS parameter file + +input,s,a,,,,Input image root file name +records,s,a,,,,Range of spectral records +output,s,a,,,,Output file root name for new spectra +start_rec,i,a,1,0,9999,Next starting spectral record +newoutput,s,q,,,,New output file root name diff --git a/noao/onedspec/irsiids/t_addsets.x b/noao/onedspec/irsiids/t_addsets.x new file mode 100644 index 00000000..cd145b4a --- /dev/null +++ b/noao/onedspec/irsiids/t_addsets.x @@ -0,0 +1,195 @@ +include <error.h> +include <imhdr.h> + + +# T_ADDSETS -- Add a series of spectra by subsets. A single spectrum +# is produced for every "subset" number of input spectra. The input +# list is accumulated until "subset" number of spectra have been +# encountered. The result is then written out. +# +# If the input data are calibrated (CA_FLAG = 0) then the result +# is an average over the subset size, but the header exposure +# time is updated. +# +# If the data is uncalibrated then the resulting spectrum is a sum +# of the total counts. + +procedure t_addsets () + +pointer image +pointer recstr, ofile +int root, start_rec, subset +int nrecs +int nrem, ifile, ca_flag +real itm, expo, wt, wtsum +bool weight +pointer sp, recs, im, cur_pix, sp_sum + +real imgetr() +int clpopni(), clgeti(), imgeti() +int get_next_image(), decode_ranges() +bool clgetb() +pointer immap(), imgl1r() + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (ofile, SZ_FNAME, TY_CHAR) + call salloc (recstr, SZ_LINE, TY_CHAR) + call salloc (recs, 300, TY_INT) + + # Open input file name template + root = clpopni ("input") + + # Get range specification if any + call clgstr ("records", Memc[recstr], SZ_LINE) + if (decode_ranges (Memc[recstr], Memi[recs], 100, nrecs) == ERR) + call error (0, "Bad range specification") + + # Get rootname for output files and starting record + call clgstr ("output", Memc[ofile], SZ_FNAME) + start_rec = clgeti ("start_rec") + + # Get subset size + subset = clgeti ("subset") + + # Apply integration time weighting? + weight = clgetb ("weighting") + + # Initialize range decoder + call reset_next_image () + + #Initialize file counter + ifile = 0 + wtsum = 0.0 + + # Loop over all input images by subsets + while (get_next_image (root, Memi[recs], nrecs, Memc[image], + SZ_FNAME) != EOF) { + + # Open image + iferr (im = immap (Memc[image], READ_ONLY, 0)) { + call erract (EA_WARN) + next + } + + # Allocate space for current subset + if (mod (ifile, subset) == 0) { + call calloc (sp_sum, IM_LEN (im,1), TY_REAL) + + # Zero exposure counter + expo = 0.0 + } + + # Add in current spectrum + iferr (itm = imgetr (im, "EXPOSURE")) + iferr (itm = imgetr (im, "ITIME")) + iferr (itm = imgetr (im, "EXPTIME")) + itm = 1 + iferr (ca_flag = imgeti (im, "CA-FLAG")) + ca_flag = -1 + cur_pix = imgl1r (im) + + # Apply integration time weighting + if (weight) + wt = itm + else + wt = 1.0 + + if (ca_flag != 0) + wt = 1.0 + + wtsum = wtsum + wt + call amulkr (Memr[cur_pix], wt, Memr[cur_pix], IM_LEN(im,1)) + call aaddr (Memr[cur_pix], Memr[sp_sum], Memr[sp_sum], IM_LEN(im,1)) + expo = expo + itm + + # Issue status report + call printf ("[%s] added\n") + call pargstr (Memc[image]) + + ifile = ifile + 1 + if (mod (ifile, subset) == 0) { + call wrt_set (Memr[sp_sum], subset, im, Memc[ofile], start_rec, + expo, wtsum, ca_flag) + wtsum = 0.0 + call mfree (sp_sum, TY_REAL) + } else + call imunmap (im) + + } + # Check that there are no remaining spectra in an unfulfilled subset + nrem = mod (ifile, subset) + if (nrem != 0) { + call wrt_set (Memr[sp_sum], nrem, im, Memc[ofile], start_rec, + expo, wtsum, ca_flag) + wtsum = 0.0 + call mfree (sp_sum, TY_REAL) + + call eprintf ("Unfulfilled subset accumulation written - ") + call eprintf ("missing %d spectra\n") + call pargi (subset - nrem) + } + + # Update record number + call clputi ("next_rec", start_rec) + + # Free space + call sfree (sp) + call clpcls (root) +end + +# WRT_SET -- Write spectra ccumulated from the set + +procedure wrt_set (sp_sum, subset, im, ofile, start_rec, expo, wtsum, ca_flag) + +real sp_sum[ARB] +int subset, start_rec, ca_flag +pointer im +char ofile[SZ_FNAME] +real expo, wtsum + +char newfile[SZ_FNAME] +pointer imnew, newpix + +pointer impl1r(), immap() +int strlen() + +begin + # Create new spectrum - first make up a name + call strcpy (ofile, newfile, SZ_FNAME) + call sprintf (newfile[strlen (newfile) + 1], SZ_FNAME, ".%04d") + call pargi (start_rec) + + imnew = immap (newfile, NEW_COPY, im) + + IM_NDIM (imnew) = 1 + IM_LEN (imnew,1) = IM_LEN (im,1) + IM_PIXTYPE (imnew) = TY_REAL + call strcpy (IM_TITLE(im), IM_TITLE(imnew), SZ_LINE) + + call imunmap (im) + + newpix = impl1r (imnew) + + # If this spectrum is calibrated, perform an average + # weighted by integration time and copy new pixels into image + if (ca_flag == 0) + call adivkr (sp_sum, real (wtsum), Memr[newpix], IM_LEN(imnew,1)) + else + call amovr (sp_sum, Memr[newpix], IM_LEN(imnew,1)) + + # Update keyword + call imaddr (imnew, "EXPOSURE", expo) + + # Send user report + call printf ("writing [%s]: %s\n") + call pargstr (newfile) + call pargstr (IM_TITLE(imnew)) + call flush (STDOUT) + + call imunmap (imnew) + + # Update record counter + start_rec = start_rec + 1 +end diff --git a/noao/onedspec/irsiids/t_bswitch.x b/noao/onedspec/irsiids/t_bswitch.x new file mode 100644 index 00000000..0c7b71a5 --- /dev/null +++ b/noao/onedspec/irsiids/t_bswitch.x @@ -0,0 +1,924 @@ +include <error.h> +include <imhdr.h> +include <mach.h> +include <time.h> +include <smw.h> + +define MAX_NR_BEAMS 100 # Max number of instrument apertures +define MIN_RANGES 100 # Minimum spectra per beam if not given + +# T_BSWITCH -- Beam switch a series of spectra to produce a single +# sky subtracted spectrum. +# +# The spectra may be extinction corrected if not already done. +# +# The summation may include an optional statistical weighting +# based on the total countrate summed over a user definable +# piece of the spectrum. If the countrate is <= 0, the +# spectrum is given zero weight. +# +# The data may be organized as data from the IIDS/IRS are usually +# obtained - where the telescope is beam-switched so that the +# object is first in one aperture while sky is observed in the other, +# and then the process is reversed. +# +# If the instrument offers many apertures, "nebular" mode can be used +# to obtain the same effect. Here all apertures observe the object(s) +# at one time; then the telescope is moved so all apertures are observing +# sky. +# +# Both these methods are considered "idsmode". But if there are a different +# number of sky observations than object, an imbalance exists. +# To account for this possibility, all summations are performed by computing +# an average countrate over all observations. Sky countrates can then be +# subtracted from the object. Later the differential countrate is returned +# to an "equivalent" count by multiplying by the exposure time. +# +# Spectra must be dispersion corrected to employ either +# weighting or extinction correction. +# +# The series of spectra may be accumulated in subsets rather than +# over the entire series by specifying a subset rate. (E.g. for +# IIDS data a subset rate of 4 would produce a summed pair for +# every quadruple.) + +# Revisions made for WCS support and change from idsmtn.h structure to shdr.h +# structure. Because this program is an awful mess the changes were made a +# small as possible without altering the structure. (5/1/91, Valdes) + +procedure t_bswitch () + +char image[SZ_FNAME,MAX_NR_BEAMS+1] +char rec_numbers[SZ_LINE], title[SZ_LINE,MAX_NR_BEAMS] +char ofile[SZ_FNAME], stat_fil[SZ_FNAME] +int sfd, nrecsx +int i, infile, nrecs, def_beam, start_rec, nimage, sub_rate +int records[300], beam_stat[MAX_NR_BEAMS], ncols[MAX_NR_BEAMS] +bool idsmode, extinct, stat, weight, eof_test +pointer ids[MAX_NR_BEAMS+1] +pointer imnames[MAX_NR_BEAMS] # Hold pointers to pointers of image names +pointer imin, sp, obs + +# The following arrays are suffixed by either 'o' for object or 's' for sky + +int ico [MAX_NR_BEAMS], ics [MAX_NR_BEAMS] # nr obs in beam +real expo [MAX_NR_BEAMS], exps [MAX_NR_BEAMS] # exposure times +pointer accumo[MAX_NR_BEAMS+1], accums[MAX_NR_BEAMS+1] # beam accumulators +pointer counto[MAX_NR_BEAMS], counts[MAX_NR_BEAMS] # counts in each obs + +int clpopni(), clgeti(), get_next_image(), decode_ranges() +int open(), mod() +pointer immap() +bool clgetb(), streq() + +begin + call smark (sp) + call aclri (ids, MAX_NR_BEAMS+1) + + # Open input filename template + infile = clpopni ("input") + + # Get range specification + call clgstr ("records", rec_numbers, SZ_LINE) + if (decode_ranges (rec_numbers, records, 100, nrecs) == ERR) + call error (0, "Bad range specification") + + # If no ranges is given, filename expansion will occur, so + # we must will need some indication of the number of spectra. + if (nrecs == MAX_INT) + nrecsx = MIN_RANGES + else + nrecsx = nrecs + + # Get root name for new records and starting record number + call clgstr ("output", ofile, SZ_FNAME) + start_rec = clgeti ("start_rec") + + # Get filename for statistics + call clgstr ("stats", stat_fil, SZ_FNAME) + + # Assume spectra are in quadruples? + idsmode = clgetb ("ids_mode") + + # Perform de-extinction? + extinct = clgetb ("extinct") + + # Use weighting? + weight = clgetb ("weighting") + + # Accumulate by subsets? - A very large number implies no subsetting + sub_rate = clgeti ("subset") + + # Open statistics file if any + if (streq (stat_fil, "")) { + sfd = NULL + stat = false + } else { + stat = true + sfd = open (stat_fil, APPEND, TEXT_FILE) + } + + # Initialize beam-switch status + obs = NULL + call init_file (extinct, def_beam, ico, ics, beam_stat) + + + # Begin cycling through all images - accumulate if possible + # by beam number + + # Initialize range decoder + call reset_next_image () + + # Set up for subsets + nimage = 0 + eof_test = false + + repeat { + + while (get_next_image (infile, records, nrecs, image[1,def_beam], + SZ_FNAME) != EOF) { + + # Attempt to open image with extended header - + iferr (imin = immap (image[1,def_beam], READ_ONLY, 0)) { + call eprintf ("[%s]") + call pargstr (image[1,def_beam]) + call error (0, "Image not found or header info not available") + } + + # Add in to accumlators + call accum_image (imin, ids, accumo, accums, counto, counts, + ico, ics, expo, exps, image, beam_stat, idsmode, extinct, + weight, nrecsx, ncols, title, imnames, sfd, obs) + + call printf ("[%s] added\n") + call pargstr (image[1,def_beam]) + call flush (STDOUT) + + # Close current image + call imunmap (imin) + + # Test for subsets + nimage = nimage + 1 + if (mod (nimage, sub_rate) == 0) + go to 10 + } + + # Get here by running out of data + eof_test = true + + # Must be careful not to write out the last sums if subsets are + # in effect because the subset check would already have done so + # We can check because "nimage" will not have been bumped + # if EOF was encountered. + + if (mod (nimage, sub_rate) != 0) { + + # All data has been summed - generate spectra of the accumlations +10 call wrt_accum (ids, image, title, accumo, accums, ico, ics, + counto, counts, expo, exps, ncols, beam_stat, idsmode, weight, + extinct, ofile, start_rec, sub_rate) + + # Generate statistics output for this beam + if (stat) + call wrt_stats (sfd, accumo, accums, ico, ics, counto, counts, + expo, exps, image, beam_stat, title, imnames, weight) + + # Clear counters and accumulators + call reset_beams (accumo, accums, expo, exps, ico, ics, beam_stat, + ncols) + } + + } until (eof_test) + + # Put current record counter back into the parameter file for + # subsequent invocations + call clputi ("next_rec", start_rec) + + # Close out inputs, outputs, and space + do i = 1, MAX_NR_BEAMS+1 + call shdr_close (ids[i]) + if (obs != NULL) + call obsclose (obs) + call clpcls (infile) + call close (sfd) + call sfree (sp) +end + +# ACCUM_IMAGE -- Opens current pixel file, loads header elements, +# adds current spectrum to accumulator array(s), +# and updates the accumulator status array. +# If not in IDSMODE, then returns both object and +# sky sums for further consideration. +# IDSMODE requires an equal number of each, object and sky, in +# a sequence of OSSO-OSSO or OSSO-SOOS groups. + +procedure accum_image (imin, ids, accumo, accums, counto, counts, ico, ics, + expo, exps, image, beam_stat, idsmode, extinct, weight, nrecs, + ncols, title, imnames, sfd, obs) + +pointer imin, ids[ARB] +pointer imnames[ARB] # Saved image names for stat printout +pointer sfd # Statistics file +pointer obs # Observatory + +pointer accumo[ARB], accums[ARB] # Object and sky accumlators +pointer counto[ARB], counts[ARB] # counting stats +real expo [ARB], exps [ARB] # total exposure times +int ico [ARB], ics [ARB] # number of observations + +char image[SZ_FNAME, MAX_NR_BEAMS+1], title[SZ_LINE,MAX_NR_BEAMS] +char observatory[SZ_FNAME] +int beam_stat[ARB], ncols[ARB] +int dum_beam +bool idsmode, extinct, weight, exflag, newobs, obshead +real latitude + +int last_len[MAX_NR_BEAMS], name_nr[MAX_NR_BEAMS] +int ifile, nr_beams, i, j, def_beam, beam_nr +int nwaves, ic, nrecs +real airm, wave1, wave2, wt +pointer wave_tbl, extn_tbl, ipacc, ipc, mw + +real clgetr(), obsgetr() +pointer smw_openim() +errchk smw_openim, shdr_open, obsimopen + +begin + # Bump image file counter + ifile = ifile + 1 + + # Load header area + mw = smw_openim (imin) + call shdr_open (imin, mw, 1, 1, INDEFI, SHDATA, ids[def_beam]) + call smw_close (MW(ids[def_beam])) + + accumo[def_beam] = SY(ids[def_beam]) + + # Check for proper flags + call flag_chk (ids[def_beam], exflag) + + if (ifile == 1) { + + # Get region for statistics to operate over - + # Currently only one set of wavelengths is available, but + # at some point, it may be desirable to extend this to + # provide a start and ending wavelength for each aperture + # since an aperture must be considered as an independent + # instrument. + + # Insert defaults --> entire spectrum + # Now ask user for start and end - if =0.0, use defaults + wave1 = clgetr ("wave1") + wave2 = clgetr ("wave2") + + if (wave1 == 0.0) + wave1 = W0(ids[def_beam]) + if (wave2 == 0.0) + wave2 = W0(ids[def_beam]) + (IM_LEN(imin,1)-1) * + WP(ids[def_beam]) + + } + + # Determine beam number and add/sub in pixels + # Remember that IIDS/IRS "beams" are 0-indexed + + beam_nr = BEAM(ids[def_beam]) + 1 + if (beam_nr > MAX_NR_BEAMS || beam_nr < 1) + call error (0, "Illegal beam number") + + # Allocate space for this aperture if not already done + # Space must be allocated for 2 lines of spectra for + # each aperture - Line 1 is used to sum up the most + # recent object-sky spectra to maintain the local + # statistics. Line 2 is used for the net accumulation + # over the entire sequence. The statistics from Line 1 + # may be used to weigh the observations as they are + # added into the Line 2 accumulation. + # + # For non-IDSMODE the two lines are used for separate + # object and sky sums + + if (IS_INDEFI (beam_stat[beam_nr])) { + beam_stat[beam_nr] = 0 + + # Allocate space for the accumulators for this beam nr + call salloc (accumo[beam_nr], IM_LEN(imin,1), TY_REAL) + call salloc (accums[beam_nr], IM_LEN(imin,1), TY_REAL) + + # Zero object and sky accumulators + call amovkr (0.0, Memr[accumo[beam_nr]], IM_LEN(imin,1)) + call amovkr (0.0, Memr[accums[beam_nr]], IM_LEN(imin,1)) + + + # Allocate space for statistics array - For each beam, + # a series of up to 'nrecs' spectra may be read, and we + # want to keep track of the stats (=countrates) for each + # observation. For non-idsmode, need sky rates too. + call salloc (counto[beam_nr], nrecs, TY_REAL) + if (!idsmode) + call salloc (counts[beam_nr], nrecs, TY_REAL) + + # Allocate space for the image names + call salloc (imnames[beam_nr], nrecs, TY_INT) + name_nr[beam_nr] = 1 + do j = 1, nrecs + call salloc (Memi[imnames[beam_nr]+j-1], SZ_LINE, TY_CHAR) + + # Save number of points for checking purposes + last_len[beam_nr] = IM_LEN(imin,1) + ncols[beam_nr] = last_len[beam_nr] + + # Initialize exposure time + expo[beam_nr] = 0.0 + exps[beam_nr] = 0.0 + + nr_beams = nr_beams + 1 + } + + # If this is an object observation, save the image name + if (OFLAG(ids[def_beam]) == 1) { + call strcpy (image[1,def_beam], Memc[Memi[imnames[beam_nr]+ + name_nr[beam_nr]-1]], SZ_LINE) + name_nr[beam_nr] = name_nr[beam_nr] + 1 + } + + # If an object observation, save the header elements -- + # NOTE that if we get >1 objects before getting a sky, only + # the last observation header is saved! + + # The pixel data will be the sum of all objects until the + # |object-sky| count = 0 -- Thus, beam switching does not + # necessarily accumulate by pairs, but depends on how the + # sequence of observations are presented to the program. + + # The following test has been deleted so that headers + # will be saved for sky frames as well. This is necessary + # if BSWITCH is to perform the function of EXTINCTION + # only when sky frames are to be written out as well. + + if (OFLAG(ids[def_beam]) == 1 || !idsmode) { + # Save headers - could probably be done faster by AMOV + call shdr_copy (ids[def_beam], ids[beam_nr], NO) + + # Fix airmass if necessary + if (extinct && IS_INDEF (AM(ids[beam_nr]))) { + call clgstr ("observatory", observatory, SZ_FNAME) + call obsimopen (obs, imin, observatory, NO, newobs, obshead) + if (newobs) { + call obslog (obs, "BSWITCH", "latitude", STDOUT) + if (sfd != NULL) + call obslog (obs, "BSWITCH", "latitude", sfd) + } + latitude = obsgetr (obs, "latitude") + call get_airm (RA(ids[beam_nr]), DEC(ids[beam_nr]), + HA(ids[beam_nr]), ST(ids[beam_nr]), latitude, + AM(ids[beam_nr])) + } + + call strcpy (image[1,def_beam], image[1,beam_nr], SZ_FNAME) + + # Save length - Each beam may be independent sizes + ncols[beam_nr] = IM_LEN(imin,1) + + # Save title, too for same reason + call strcpy (IM_TITLE(imin), title[1,beam_nr], SZ_LINE) + } + + # Verify length + if (last_len[beam_nr] != ncols[beam_nr]) { + call eprintf ("[%s] -- Length not consistent %d\n") + call pargstr (image[1,beam_nr]) + call pargi (ncols[beam_nr]) + ncols[beam_nr] = min (ncols[beam_nr], last_len[beam_nr]) + } + last_len[beam_nr] = ncols[beam_nr] + + + # Check to see if a pair is obtained - then perform statistics + # and add into global accumulator + + if (idsmode) { + + # Add spectrum to local accumulation buffer --> Use SKY buffer + # At this point of deriving a sequentially local sum, weighting + # is not used. + + call add_spec (Memr[accumo[def_beam]], Memr[accums[beam_nr]], + beam_stat[beam_nr], OFLAG(ids[def_beam]), last_len[beam_nr]) + + # IDSMODE requires that every 2N observations produce an + # OBJECT-SKY pair + if (mod (ifile, 2*nr_beams) == 0) + + # Review all beams in use for non-zero pairings + do i = 1, MAX_NR_BEAMS + if (!IS_INDEFI (beam_stat[i]) && beam_stat[i] != 0) + call error (0, "Spectra are not in quadruples") + + + # Object and sky exposure times must be equal. + if (OFLAG(ids[def_beam]) == 1) { + expo[beam_nr] = expo[beam_nr] + IT(ids[def_beam]) + + # Increment number of object observations for this beam + ico[beam_nr] = ico[beam_nr] + 1 + } + + + if (beam_stat[beam_nr] == 0) { + # Add up all counts within a region for statistics of objects + # This must be kept separately for each beam number and for + # each observation + + # First convert to counts per second (CPS) + call adivkr (Memr[accums[beam_nr]], IT(ids[def_beam]), + Memr[accums[beam_nr]], last_len[beam_nr]) + + # Sum CPS in statistics region + call sum_spec (Memr[accums[beam_nr]], wave1, wave2, + W0(ids[def_beam]), WP(ids[def_beam]), Memr[counto[beam_nr]+ + ico[beam_nr]-1], last_len[beam_nr]) + + # De-extinct spectrum + if (extinct && !exflag) { + airm = AM(ids[beam_nr]) + call de_ext_spec (Memr[accums[beam_nr]], airm, + W0(ids[def_beam]), WP(ids[def_beam]), Memr[wave_tbl], + Memr[extn_tbl], nwaves, last_len[beam_nr]) + } + + # Add to global accumulator + # Use weights which are proportional to countrate, if desired + if (weight) { + wt = Memr[counto[beam_nr]+ico[beam_nr]-1] + call amulkr (Memr[accums[beam_nr]], wt, + Memr[accums[beam_nr]], last_len[beam_nr]) + } + + # And add into global sum + call aaddr (Memr[accums[beam_nr]], Memr[accumo[beam_nr]], + Memr[accumo[beam_nr]], last_len[beam_nr]) + } + + } else { + # Non IDSMODE -accumulate separate object and sky CPS sums + + # Set pointers and update obj-sky parameters + if (OFLAG(ids[def_beam]) == 1) { + beam_stat[beam_nr] = beam_stat[beam_nr] + 1 + ipacc = accumo[beam_nr] + ipc = counto[beam_nr] + ico[beam_nr] = ico[beam_nr] + 1 + ic = ico[beam_nr] + expo[beam_nr] = expo[beam_nr] + IT(ids[def_beam]) + } else { + beam_stat[beam_nr] = beam_stat[beam_nr] - 1 + ipacc = accums[beam_nr] + ipc = counts[beam_nr] + ics[beam_nr] = ics[beam_nr] + 1 + ic = ics[beam_nr] + exps[beam_nr] = exps[beam_nr] + IT(ids[def_beam]) + } + + # First convert to counts per second (CPS) + call adivkr (Memr[accumo[def_beam]], IT(ids[def_beam]), + Memr[accumo[def_beam]], last_len[beam_nr]) + + # Get counting stats + call sum_spec (Memr[accumo[def_beam]], wave1, wave2, + W0(ids[def_beam]), WP(ids[def_beam]), Memr[ipc+ic-1], + last_len[beam_nr]) + + # De-extinct spectrum + if (extinct && !exflag) { + airm = AM(ids[beam_nr]) + call de_ext_spec (Memr[accumo[def_beam]], airm, + W0(ids[def_beam]), WP(ids[def_beam]), Memr[wave_tbl], + Memr[extn_tbl], nwaves, last_len[beam_nr]) + } + + if (weight) { + wt = Memr[ipc+ic-1] + call amulkr (Memr[accumo[def_beam]], wt, Memr[accumo[def_beam]], + last_len[beam_nr]) + } + + # Add into appropriate accumulator + call aaddr (Memr[accumo[def_beam]], Memr[ipacc], Memr[ipacc], + last_len[beam_nr]) + } + + return + +# INIT_FILE -- Zero the file initializer, the beam counter, beam stats +# and read the extinction data if necessary + +entry init_file (extinct, dum_beam, ico, ics, beam_stat) + + ifile = 0 + nr_beams = 0 + def_beam = MAX_NR_BEAMS + 1 + dum_beam = def_beam + + do i = 1, MAX_NR_BEAMS { + beam_stat[i] = INDEFI + ico[i] = 0 + ics[i] = 0 + } + + # If extinction required, read in extinction file, and sensitivity file + if (extinct) + call get_extn (wave_tbl, extn_tbl, nwaves) + + return + +# INIT_NAME -- Reset name index counter for a beam number + +entry init_name (dum_beam) + + name_nr[dum_beam] = 1 + return +end + +# ACCUM_OUT -- Checks accumulator flags and writes out a new summed +# image if the count is zero + +procedure accum_out (accum, image, ncols, title, root, rec, beam_nr, + bsflag, itm, exflag) + +real accum[ARB], itm +char image[SZ_FNAME], title[SZ_LINE], root[SZ_FNAME] +int ncols, rec, beam_nr +int bsflag, exflag + +pointer imin, imout, spec +char bs_image[SZ_FNAME] + +pointer immap(), impl1r() + +begin + # Create new image with user area + # Use ROOT for spectrum name and increment starting record number + + call sprintf (bs_image, SZ_FNAME, "%s.%04d") + call pargstr (root) + call pargi (rec) + + rec = rec + 1 + + # Provide user info + call printf ("writing: [%s] %s\n") + call pargstr (bs_image) + call pargstr (title) + call flush (STDOUT) + + imin = immap (image, READ_ONLY, 0) + imout = immap (bs_image, NEW_COPY, imin) + + # Add standard image header + IM_NDIM(imout) = 1 + IM_LEN(imout,1) = ncols + IM_PIXTYPE(imout) = TY_REAL + call strcpy (title, IM_TITLE(imout), SZ_LINE) + + # Write out pixels + spec = impl1r (imout) + call amovr (accum, Memr[spec], ncols) + + # Update changed parameters + if(bsflag == 1) + call imaddi (imout, "BS-FLAG", bsflag) + call imaddr (imout, "EXPTIME", itm) + call imaddi (imout, "EX-FLAG", exflag) + + call imunmap (imin) + call imunmap (imout) + + # Store new image name back into image + call strcpy (bs_image, image, SZ_FNAME) +end + +# ACCUM_NORM - Normalize weighted rate and convert to counts + +procedure accum_norm (accum, nr, counts, exp, ncols, weight) + +real accum[ARB], counts[ARB], exp +int nr, ncols +bool weight + +real sum_wt +int i + +begin + # The accumulation is an array weighted by non-normalized weights + # Normalize to total weight to produce a true weighted average + # and multiply by the total exposure to produce + # an equivalent sum + + # Add up all weighting factors + if (weight) { + sum_wt = 0.0 + do i = 1, nr + sum_wt = sum_wt + counts[i] + } else + sum_wt = real (nr) + + if (sum_wt == 0.0) + sum_wt = 1.0 + + # Correct for exposure time + sum_wt = exp / sum_wt + + call amulkr (accum, sum_wt, accum, ncols) +end + +# WRT_ACCUM -- Write out accumulations as spectra + +procedure wrt_accum (ids, image, title, accumo, accums, ico, ics, + counto, counts, expo, exps, ncols, beam_stat, idsmode, weight, + extinct, ofile, start_rec, sub_rate) + +pointer ids[ARB] +char image[SZ_FNAME,MAX_NR_BEAMS+1], title[SZ_LINE,MAX_NR_BEAMS] + +pointer accumo[ARB], accums[ARB] +pointer counto[ARB], counts[ARB] +int ico [ARB], ics [ARB] +real expo [ARB], exps [ARB] +int ncols[ARB] +int beam_stat[ARB] +bool idsmode, weight, extinct +char ofile[SZ_FNAME] +int start_rec, sub_rate, bsflag + +int i, nr_beams +real exp_ratio + +begin + # First compute number of beams + nr_beams = 0 + do i = 1, MAX_NR_BEAMS + if (!IS_INDEFI (beam_stat[i]) && ((ico[i] > 0) || (ics[i] > 0))) + nr_beams = nr_beams + 1 + + # For all present apertures, write out a spectrum + do i = 1, MAX_NR_BEAMS { + + if (!IS_INDEFI (beam_stat[i]) && ((ico[i] > 0) || (ics[i] > 0))) { + if (beam_stat[i] != 0 && idsmode) { + call eprintf ("Non-equal number of obj-sky observations") + call eprintf (" beam: %d - residual: %d\n") + call pargi (i-1) + call pargi (beam_stat[i]) + + # Reset to 0 and force output + beam_stat[i] = 0 + } + + # The accumulator has a total CPS using non-normalized + # weights - apply normalization and exposure time to + # generate an equivalent COUNT sum. + call accum_norm (Memr[accumo[i]], ico[i], Memr[counto[i]], + expo[i], ncols[i], weight) + + if (!idsmode) { + # Separate object and sky sums require sky info + call accum_norm (Memr[accums[i]], ics[i], Memr[counts[i]], + exps[i], ncols[i], weight) + + # Then normalize sky exposure time to that of object + if (exps[i] != 0.0) + exp_ratio = expo[i]/exps[i] + else + exp_ratio = 1.0 + + # Check that some object observtion was made + # If not, then we only have sky data so multiply by -1 + # so that the subsequent subtraction will produce a + # positive sky + if (expo[i] == 0.0) + exp_ratio = -1.0 + + if (exp_ratio != 1.0) + call amulkr (Memr[accums[i]], exp_ratio, + Memr[accums[i]], ncols[i]) + + # Finally subtract sky from object equivalent counts + call asubr (Memr[accumo[i]], Memr[accums[i]], + Memr[accumo[i]], ncols[i]) + + } + # Set header flags + # BS flag is not set if the subset rate equals the + # number of apertures since each record in is copied out + if (sub_rate > nr_beams) + bsflag = 1 + else + bsflag = -1 + + if (OFLAG(ids[i]) == 1) + IT(ids[i]) = expo[i] + else + IT(ids[i]) = exps[i] + + if (extinct) + EC(ids[i]) = 0 + + # And write out spectrum, at last + call accum_out (Memr[accumo[i]], image[1,i], + ncols[i], title[1,i], ofile, start_rec, i, + bsflag, IT(ids[i]), EC(ids[i])) + + # Reset name entry counter + call init_name (i) + } + + } +end + +# RESET_BEAMS -- Zeroes the counters and accumulators for additional +# cases + +procedure reset_beams (accumo, accums, expo, exps, ico, ics, beam_stat, ncols) + +pointer accumo[ARB], accums[ARB] +real expo [ARB], exps [ARB] +int ico [ARB], ics [ARB] +int beam_stat[ARB], ncols[ARB] + +int i + +begin + do i = 1, MAX_NR_BEAMS + if (!IS_INDEFI (beam_stat[i])) { + + expo[i] = 0.0 + exps[i] = 0.0 + ico[i] = 0 + ics[i] = 0 + call amovkr (0.0, Memr[accumo[i]], ncols[i]) + call amovkr (0.0, Memr[accums[i]], ncols[i]) + } +end + +# WRT_STATS -- Write out statistics file + +procedure wrt_stats (fd, accumo, accums, ico, ics, counto, counts, + expo, exps, image, beam_stat, title, imnames, weight) + +int fd +pointer accumo[ARB], accums[ARB], counto[ARB], counts[ARB] +real expo[ARB], exps[ARB] +int ico[ARB], ics[ARB], beam_stat[ARB] +char image[SZ_FNAME,MAX_NR_BEAMS+1] +char title[ARB] +pointer imnames[ARB] +bool weight + +int i, j +real cmaxo, cmaxs +char ctime[SZ_TIME] + +long clktime() + +begin + # Issue time stamp + call cnvtime (clktime (long(0)), ctime, SZ_TIME) + call fprintf (fd, "%s\n\n") + call pargstr (ctime) + + # Issue message if weighted sums in effect + if (weight) + call fprintf (fd, "--> Using weighted averages <--\n\n") + + # Cycle over beams + do i = 1, MAX_NR_BEAMS { + if (!IS_INDEFI (beam_stat[i])) { + + # Write out Object stats if any + if (ico[i] > 0) { + call fprintf (fd, "Object statistics for beam %d -->[%s]\n") + call pargi (i-1) + call pargstr (image[1,i]) + call fprintf (fd, "Title: %s\n") + call pargstr (title) + + # Find maximum count value for this beam + cmaxo = Memr[counto[i]] + + do j = 1, ico[i] + cmaxo = max (cmaxo, Memr[counto[i]+j-1]) + + call fprintf (fd, + "Obs Relative CPS Image%12wPeak CPS = %10.3g\n") + call pargr (cmaxo) + + if (cmaxo == 0.0) + cmaxo = 1.0 + + do j = 1, ico[i] { + call fprintf (fd, "%3d %8.2f [%s]\n") + call pargi (j) + call pargr (Memr[counto[i]+j-1] / cmaxo) + call pargstr (Memc[Memi[imnames[i]+j-1]]) + } + } + + # Write out sky stats if any + if (ics[i] > 0) { + call fprintf (fd, "Sky statistics for beam %d\n") + call pargi (i-1) + + cmaxs = Memr[counts[i]] + + do j = 1, ics[i] + cmaxs = max (cmaxs, Memr[counts[i]+j-1]) + + call fprintf (fd, "Obs Relative CPS Peak CPS = %10.3g\n") + call pargr (cmaxs) + + if (cmaxs == 0.0) + cmaxs = 1.0 + + do j = 1, ics[i] { + call fprintf (fd, "%3d %8.2f\n") + call pargi (j) + call pargr (Memr[counts[i]+j-1] / cmaxs) + } + } + + call fprintf (fd, "\n\n") + } + } +end + + +# ADD_SPEC -- Accumulate spectrum into array - either add or subtract +# Returns status = net number of object - sky apectra +# = 0 for equal numbers to indicate further +# processing may take place + +procedure add_spec (inspec, accum, stat, flag, len) + +real inspec[ARB], accum[ARB] +int stat, flag, len + +int i, add_sub + +begin + add_sub = 0 + + # Is this an Object or Sky? + # If flag is neither 0 or 1, spectrum is ignored + if (flag == 0) + add_sub = -1 + if (flag == 1) + add_sub = +1 + + if (add_sub == 0) { + stat = INDEFI + return + } + + # Is accumulator to be cleared? + if (IS_INDEFI (stat) || stat == 0) { + call amulkr (inspec, real (add_sub), accum, len) + stat = add_sub + + } else { + # Add into accumulator + do i = 1, len + accum[i] = accum[i] + add_sub * inspec[i] + + stat = stat + add_sub + } +end + +# FLAG_CHK -- Check header flags prior to beam switching + +procedure flag_chk (ids, exflag) + +pointer ids +bool exflag + +int bsflag, imgeti() + +begin + # BS requires + # 1. dispersion corrected spectra + # 2. non-beam switched + # 3. may be either extinction corrected or not + + if (DC(ids) != DCLINEAR) + call error (0, "Spectrum not dispersion corrected") + + iferr (bsflag = imgeti (IM(ids), "BS-FLAG")) + bsflag = -1 + if (bsflag == 1) + call error (0, "Spectrum already beam-switched") + + if (EC(ids) == ECYES) + exflag = true + else + exflag = false +end diff --git a/noao/onedspec/irsiids/t_coefs.x b/noao/onedspec/irsiids/t_coefs.x new file mode 100644 index 00000000..656e777d --- /dev/null +++ b/noao/onedspec/irsiids/t_coefs.x @@ -0,0 +1,88 @@ +include <error.h> + +# COEFS -- Convert IIDS/IRS coeffients to IDENTIFY database entry. + +procedure t_coefs () + +int root # List of input root names +pointer database # Output database directory + +int i, nrecs, ncoefs +real coef +pointer sp, image, dtname, recs, im, dt + +real imgetr() +int clpopni(), imgeti(), get_next_image(), decode_ranges() +pointer immap(), dtmap1() +errchk imgetr, dtmap1 + +begin + call smark (sp) + call salloc (image, SZ_LINE, TY_CHAR) + call salloc (database, SZ_FNAME, TY_CHAR) + call salloc (dtname, SZ_FNAME, TY_CHAR) + call salloc (recs, 300, TY_INT) + + root = clpopni ("input") + call clgstr ("records", Memc[image], SZ_LINE) + call clgstr ("database", Memc[database], SZ_LINE) + + if (decode_ranges (Memc[image], Memi[recs], 100, nrecs) == ERR) + call error (0, "Bad range specification") + + # Loop over all input images - print name on STDOUT + call reset_next_image () + while (get_next_image (root, Memi[recs], nrecs, Memc[image], + SZ_LINE) != EOF) { + iferr (im = immap (Memc[image], READ_ONLY, 0)) { + call erract (EA_WARN) + next + } + + iferr (ncoefs = imgeti (im, "DF-FLAG")) + ncoefs = -1 + if (ncoefs > 1) { + call strcpy ("id", Memc[dtname], SZ_FNAME) + call imgcluster (Memc[image], Memc[dtname+2], SZ_FNAME) + dt = dtmap1 (Memc[database], Memc[dtname], APPEND) + + call dtptime (dt) + call dtput (dt, "begin\tidentify %s\n") + call pargstr (Memc[image]) + call dtput (dt, "\tid\t%s\n") + call pargstr (Memc[image]) + call dtput (dt, "\ttask\tidentify\n") + call dtput (dt, "\timage\t%s\n") + call pargstr (Memc[image]) + + # Convert coefficients + call dtput (dt, "\tcoefficients\t%d\n") + call pargi (ncoefs+4) + call dtput (dt, "\t\t2\n") + call dtput (dt, "\t\t%1d\n") + call pargi (ncoefs) + call dtput (dt, "\t\t1\n") + call dtput (dt, "\t\t1024\n") + + do i = 1, ncoefs { + call sprintf (Memc[dtname], SZ_FNAME, "DF%d") + call pargi (i) + coef = imgetr (im, Memc[dtname]) + call dtput (dt, "\t\t%10.4f\n") + call pargr (coef) + } + + call dtput (dt, "\n") + call dtunmap (dt) + } + + call printf ("[%s] %d coefficients written\n") + call pargstr (Memc[image]) + call pargi (max (0, ncoefs)) + call flush (STDOUT) + call imunmap (im) + } + + call clpcls (root) + call sfree (sp) +end diff --git a/noao/onedspec/irsiids/t_coincor.x b/noao/onedspec/irsiids/t_coincor.x new file mode 100644 index 00000000..ad2d8bd4 --- /dev/null +++ b/noao/onedspec/irsiids/t_coincor.x @@ -0,0 +1,102 @@ +include <error.h> +include <imhdr.h> + + +# T_COINCOR -- Apply coincidence corrections to spectra + +procedure t_coincor () + +int root, start_rec, ccmode, npts, nrecs, coflag +real dtime, power, expo +pointer sp, image, ofile, str, recs, imin, imout, pixin, pixout + +int clpopni(), clgeti(), clgwrd(), imgeti() +int get_next_image(), decode_ranges() +real clgetr(), imgetr() +pointer immap(), imgl1r(), impl1r() +errchk coincor + +begin + # Allocate memory + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (ofile, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (recs, 300, TY_INT) + + # Get parameters + root = clpopni ("input") + call clgstr ("output", Memc[ofile], SZ_FNAME) + if (Memc[ofile] != EOS) + start_rec = clgeti ("start_rec") + ccmode = clgwrd ("ccmode", Memc[str], SZ_LINE, ",photo,iids,") + dtime = clgetr ("deadtime") + power = clgetr ("power") + call clgstr ("records", Memc[str], SZ_LINE) + + # Initialize + if (decode_ranges (Memc[str], Memi[recs], 100, nrecs) == ERR) + call error (0, "Bad range specification") + call reset_next_image () + + # Loop over all input images by subsets + while (get_next_image (root, Memi[recs], nrecs, Memc[image], + SZ_FNAME) != EOF) { + + # Open input image and check coincidence flag + iferr (imin = immap (Memc[image], READ_WRITE, 0)) { + call erract (EA_WARN) + start_rec = start_rec + 1 + next + } + iferr (coflag = imgeti (imin, "CO-FLAG")) + coflag = -1 + if (coflag > 0) { + call printf ("[%s] already coincidence corrected\n") + call pargstr (IM_HDRFILE(imin)) + call flush (STDOUT) + call imunmap (imin) + next + } + + # Open output image + if (Memc[ofile] != EOS) { + call sprintf (Memc[str], SZ_LINE, "%s.%04d") + call pargstr (Memc[ofile]) + call pargi (start_rec) + start_rec = start_rec + 1 + + imout = immap (Memc[str], NEW_COPY, imin) + IM_PIXTYPE (imout) = TY_REAL + } else + imout = imin + + # Apply coincidence correction + pixin = imgl1r (imin) + pixout = impl1r (imout) + npts = IM_LEN (imin, 1) + iferr (expo = imgetr (imin, "EXPOSURE")) + iferr (expo = imgetr (imin, "ITIME")) + iferr (expo = imgetr (imin, "EXPTIME")) + expo = 1 + call coincor (Memr[pixin], Memr[pixout], npts, expo, coflag, + dtime, power, ccmode) + + # Update flag and write status + call imaddi (imout, "CO-FLAG", coflag) + call printf ("[%s] --> [%s] %s\n") + call pargstr (IM_HDRFILE(imin)) + call pargstr (IM_HDRFILE(imout)) + call pargstr (IM_TITLE(imout)) + call flush (STDOUT) + + # Close images + if (imout != imin) + call imunmap (imout) + call imunmap (imin) + } + + call clputi ("next_rec", start_rec) + call clpcls (root) + call sfree (sp) +end diff --git a/noao/onedspec/irsiids/t_flatdiv.x b/noao/onedspec/irsiids/t_flatdiv.x new file mode 100644 index 00000000..37186878 --- /dev/null +++ b/noao/onedspec/irsiids/t_flatdiv.x @@ -0,0 +1,276 @@ +include <imhdr.h> +include <error.h> + +define MAX_NR_BEAMS 100 # Max number of instrument apertures + +# T_FLATDIV -- Divide by a flat field spectrum. This is basically +# a simple division of two vectors but with the following +# additional functions: +# +# 1. Check the processing flag of the input spectra to avoid +# double processing, and set the flag if the processing is +# performed. +# 2. Trap division by zero errors +# 3. Optionally apply coincidence corrections + +procedure t_flatdiv () + +int root, start_rec +int nrecs +int len_flat +int ccmode, qd_flag +real dtime +real power +bool coincidence +pointer sp, image, str, ofile, flat, recs, bstat, flatsp, im + +int clpopni(), clgeti(), clgwrd(), imgeti() +int get_next_image(), decode_ranges() +real clgetr() +bool clgetb() +pointer immap() +errchk get_flatsp + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (ofile, SZ_FNAME, TY_CHAR) + call salloc (flat, SZ_FNAME, TY_CHAR) + call salloc (recs, 300, TY_INT) + call salloc (bstat, MAX_NR_BEAMS, TY_INT) + + # Open input file name template + root = clpopni ("input") + + # Get range specification if any + call clgstr ("records", Memc[str], SZ_LINE) + if (decode_ranges (Memc[str], Memi[recs], 100, nrecs) == ERR) + call error (0, "Bad range specification") + + # Get rootname for output files and starting record + # Subtract 1 from start_rec because 1 will be added later. + call clgstr ("output", Memc[ofile], SZ_FNAME) + start_rec = clgeti ("start_rec") - 1 + + # Get flat field spectrum root name + call clgstr ("flat_file", Memc[flat], SZ_FNAME) + + # Apply coincidence corrections? + coincidence = clgetb ("coincor") + if (coincidence) { + ccmode = clgwrd ("ccmode", Memc[str], SZ_LINE, ",photo,iids,") + dtime = clgetr ("deadtime") + power = clgetr ("power") + } + + # Initialize beam number status + call init_bs (Memi[bstat]) + + # Initialize range decoder + call reset_next_image () + + # Loop over all input images - divide and make new image. + # The output record number is incremented in all cases. + while (get_next_image (root, Memi[recs], nrecs, Memc[image], + SZ_FNAME) != EOF) { + start_rec = start_rec + 1 + + # Open image + iferr (im = immap (Memc[image], READ_ONLY, 0)) { + call erract (EA_WARN) + next + } + + # Get header + iferr (qd_flag = imgeti (im, "QD-FLAG")) + qd_flag = -1 + + # Verify divide flag + if (qd_flag != 0) { + + # Get flat field spectrum if needed + call get_flatsp (im, flatsp, Memc[flat], Memi[bstat], len_flat) + + # Calibrate the current spectrum and make a calibrated version + call divide (im, flatsp, len_flat, Memc[image], Memc[ofile], + start_rec, coincidence, ccmode, dtime, power) + + } else { + call eprintf ("[%s] already divided - ignored\n") + call pargstr (Memc[image]) + } + } + + # Update record number + call clputi ("next_rec", start_rec) + + # Free space + call sfree (sp) + call clpcls (root) +end + +# GET_FLATSP -- Load flat field spectrum for the current beam number + +procedure get_flatsp (im, sp, flat_file, beam_stat, len_flat) + +pointer im, sp +char flat_file[SZ_FNAME] +int beam_stat[ARB], len_flat + +int i +int beam, len[MAX_NR_BEAMS] +char sfname[SZ_FNAME] +pointer flatsp[MAX_NR_BEAMS], imflat + +int strlen(), imgeti() +pointer imgl1r(), immap() +errchk immap + +begin + # Determine beam number. + + iferr (beam = imgeti (im, "BEAM-NUM")) + beam = 0 + beam = beam + 1 + + # Validate beam number + if (beam < 1 || beam > MAX_NR_BEAMS) { + call eprintf (" Beam number out of range: %d - using 0\n") + call pargi (beam) + beam = 1 + } + + # Has this beam already been loaded? + if (IS_INDEFI (beam_stat[beam])) { + + # Create file name + call strcpy (flat_file, sfname, SZ_FNAME) + + # Flat field file names have beam number appended + call sprintf (sfname[strlen(sfname)+1], SZ_FNAME, ".%04d") + call pargi (beam-1) + + # Open spectrum + imflat = immap (sfname, READ_ONLY, 0) + + # Allocate space for this beam's sensitivity spectrum + call salloc (flatsp[beam], IM_LEN(imflat,1), TY_REAL) + + # Copy pixels into space + call amovr (Memr[imgl1r(imflat)], Memr[flatsp[beam]], + IM_LEN(imflat,1)) + + # Must be careful that no division by zero occurs. + do i = 1, IM_LEN(imflat,1) { + if (Memr[flatsp[beam]+i-1] == 0.0) + Memr[flatsp[beam]+i-1] = 1.0 + } + + # Mark this beam accounted for + beam_stat[beam] = 1 + len[beam] = IM_LEN(imflat,1) + + call imunmap (imflat) + } + + # Point to the spectrum + sp = flatsp[beam] + len_flat = len[beam] + +end + +# DIVIDE -- Perform the division and create new spectrum + +procedure divide (im, flat, len_flat, ifile, ofile, rec, + coincidence, ccmode, dtime, power) + +pointer im, flat +int len_flat, rec, ccmode +real dtime, power +char ifile[ARB], ofile[ARB] +bool coincidence + +real itm, imgetr() +int i, co_flag, imgeti() +int ncols, nlines +char calfname[SZ_FNAME], original[SZ_FNAME] +pointer imcal, rawpix, calpix + +pointer immap(), impl2r(), imgl2r() + +begin + # Find smallest length of the two possible spectra + ncols = min (IM_LEN (im, 1), len_flat) + + # Create new spectrum. Make up a name + call sprintf (calfname, SZ_FNAME, "%s.%04d") + call pargstr (ofile) + call pargi (rec) + + call xt_mkimtemp (ifile, calfname, original, SZ_FNAME) + imcal = immap (calfname, NEW_COPY, im) + + IM_NDIM(imcal) = IM_NDIM(im) + IM_LEN (imcal,1) = ncols + IM_PIXTYPE(imcal) = TY_REAL + + # Check for 2D spectrum + if (IM_NDIM(im) > 1) + nlines = IM_LEN(im,2) + else + nlines = 1 + + # Copy across the image title + call strcpy (IM_TITLE(im), IM_TITLE(imcal), SZ_LINE) + + # Operate on the pixels + do i = 1, nlines { + rawpix = imgl2r (im,i) + calpix = impl2r (imcal,i) + + # Apply coincidence correction if needed + co_flag = -1 + if (coincidence) { + iferr (co_flag = imgeti (im, "CO-FLAG")) + ; + if (co_flag < 1) { + iferr (itm = imgetr (im, "EXPOSURE")) + iferr (itm = imgetr (im, "ITIME")) + iferr (itm = imgetr (im, "EXPTIME")) + itm = 1. + call coincor (Memr[rawpix], Memr[rawpix], ncols, + itm, co_flag, dtime, power, ccmode) + } + } + + call adivr (Memr[rawpix], Memr[flat], Memr[calpix], ncols) + } + + call imaddi (imcal, "QD-FLAG", 0) + if (co_flag != -1) + call imaddi (imcal, "CO-FLAG", co_flag) + + # Send user report + call printf ("writing [%s]: %s\n") + call pargstr (original) + call pargstr (IM_TITLE(imcal)) + call flush (STDOUT) + + call imunmap (im) + call imunmap (imcal) + call xt_delimtemp (calfname, original) +end + +# INIT_BS -- Initialize beam status flags + +procedure init_bs (beam_stat) + +int beam_stat[ARB] + +int i + +begin + do i = 1, MAX_NR_BEAMS + beam_stat[i] = INDEFI +end diff --git a/noao/onedspec/irsiids/t_flatfit.x b/noao/onedspec/irsiids/t_flatfit.x new file mode 100644 index 00000000..c3558afc --- /dev/null +++ b/noao/onedspec/irsiids/t_flatfit.x @@ -0,0 +1,740 @@ +include <imhdr.h> +include <math/curfit.h> +include <gset.h> + +define MAX_NR_BEAMS 100 # Max number of instrument apertures + +define KEY "noao$lib/scr/flatfit.key" +define PROMPT "flatfit cursor options" + +# Definitions for Plotting modes +define PLT_FIT 1 # Plot the direct fit +define PLT_ERR 2 # Plot the errors in the fit +define PLT_LIN 3 # Plot the fit minus the linear part + +# T_FLATFIT -- Accumulate a series of flat field spectra to produce +# a grand sum and fit a function to the sum to produce a normalized +# flat containing the pixel-to-pixel variations. +# User interaction via the graphics cursor is provided. The following +# cursor commands are recognized: +# +# ? - Screen help +# / - Status line help +# e - Plot in residual error mode +# f - Plot in fit to the data mode +# o - Change order of fit +# l - Change lower rejection sigma +# u - Change upper rejection sigma +# r - Reset fit to include rejected pixels +# s - Change upper and lower sigmas to same value +# i - Iterate again +# n - Iterate N times +# q - Quit and accept current solution (also RETURN) +# + +procedure t_flatfit () + +pointer image # Image name to be fit +pointer images # Image name to be fit +pointer ofile # Output image file name +int function # Fitting function +int order # Order of fitting function +int recs # Spectral record numbers +int root, nrecs # CL and ranges flags +real expo # Exposure time +real dtime # Deadtime +real power # Power law coin. correction +real lower # Lower rejection sigma +real upper # Upper threshold sigma +int ngrow # Rejection radius +real div_min # Division min for option RESP +bool coincidence, all # Apply coincidence correction +bool interact # Interactive levels +pointer bstat # Status of each aperture +pointer npts # Length of spectrum +pointer esum # Accumulated exposure time +pointer accum # Pointers to beam accumulators +pointer title +int ccmode, beam +int niter + +int i +pointer sp, str, im + +int clgeti(), clgwrd(), clpopni(), imgeti() +int get_next_image(), decode_ranges() +real clgetr(), imgetr() +bool clgetb() +pointer immap() + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (images, MAX_NR_BEAMS, TY_POINTER) + call salloc (ofile, SZ_FNAME, TY_CHAR) + call salloc (recs, 300, TY_INT) + call salloc (bstat, MAX_NR_BEAMS, TY_INT) + call salloc (npts, MAX_NR_BEAMS, TY_INT) + call salloc (esum, MAX_NR_BEAMS, TY_REAL) + call salloc (accum, MAX_NR_BEAMS, TY_POINTER) + call salloc (title, MAX_NR_BEAMS, TY_POINTER) + call salloc (str, SZ_LINE, TY_CHAR) + call amovki (NULL, Memi[images], MAX_NR_BEAMS) + + # Get task parameters. + root = clpopni ("input") + + # Get input record numbers + call clgstr ("records", Memc[str], SZ_LINE) + if (decode_ranges (Memc[str], Memi[recs], 100, nrecs) == ERR) + call error (0, "Bad range specification") + + call clgstr ("output", Memc[ofile], SZ_LINE) + + call clgcurfit ("function", "order", function, order) + + lower = clgetr ("lower") + upper = clgetr ("upper") + ngrow = clgeti ("ngrow") + div_min = clgetr ("div_min") + + # Determine desired level of activity + interact = clgetb ("interact") + all = clgetb ("all_interact") + + niter = clgeti ("niter") + + # Is coincidence correction to be performed? + coincidence = clgetb ("coincor") + + if (coincidence) { + ccmode = clgwrd ("ccmode", Memc[str], SZ_LINE, ",photo,iids,") + dtime = clgetr ("deadtime") + power = clgetr ("power") + } + + call reset_next_image () + + # Clear all beam status flags + call amovki (INDEFI, Memi[bstat], MAX_NR_BEAMS) + call aclrr (Memr[esum], MAX_NR_BEAMS) + + call printf ("Accumulating spectra --\n") + call flush (STDOUT) + +10 while (get_next_image (root, Memi[recs], nrecs, Memc[image], + SZ_FNAME) != EOF) { + iferr (im = immap (Memc[image], READ_ONLY, 0)) { + call eprintf ("Header info not available for [%s]\n") + call pargstr (Memc[image]) + goto 10 + } + + iferr (beam = imgeti (im, "BEAM-NUM")) + beam = 0 + if (beam < 0 || beam > MAX_NR_BEAMS-1) + call error (0, "Invalid aperture number") + + iferr (expo = imgetr (im, "EXPOSURE")) + iferr (expo = imgetr (im, "ITIME")) + iferr (expo = imgetr (im, "EXPTIME")) + expo = 1 + + # Add spectrum into accumulator + if (IS_INDEFI (Memi[bstat+beam])) { + Memi[npts+beam] = IM_LEN (im,1) + call salloc (Memi[accum+beam], Memi[npts+beam], TY_REAL) + call aclrr (Memr[Memi[accum+beam]], Memi[npts+beam]) + Memi[bstat+beam] = 0 + + call salloc (Memi[title+beam], SZ_LINE, TY_CHAR) + call strcpy (IM_TITLE(im), Memc[Memi[title+beam]], SZ_LINE) + } + + call ff_accum_spec (im, Memi[npts], expo, Memi[bstat], beam+1, + Memi[accum], Memr[esum], coincidence, ccmode, dtime, power, + Memi[title]) + + call printf ("[%s] added to aperture %1d\n") + call pargstr (Memc[image]) + call pargi (beam) + call flush (STDOUT) + if (Memi[images+beam] == NULL) + call salloc (Memi[images+beam], SZ_FNAME, TY_CHAR) + call strcpy (Memc[image], Memc[Memi[images+beam]], SZ_FNAME) + + call imunmap (im) + } + + # Review all apertures containing data and perform fits. + # Act interactively if desired + do i = 0, MAX_NR_BEAMS-1 { + if (!IS_INDEFI (Memi[bstat+i])) { + call fit_spec (Memr[Memi[accum+i]], Memi[npts+i], Memr[esum+i], + interact, function, order, niter, lower, upper, ngrow, + div_min, i) + if (interact & !all) + interact = false + call wrt_fit_spec (Memc[Memi[images+i]], Memr[Memi[accum+i]], + Memr[esum+i], Memc[ofile], i, Memc[Memi[title+i]], + Memi[npts+i], order) + } + } + + call sfree (sp) + call clpcls (root) +end + +# ACCUM_SPEC -- Accumulate spectra by beams + +procedure ff_accum_spec (im, len, expo, beam_stat, beam, accum, expo_sum, + coincidence, ccmode, dtime, power, title) + +pointer im, accum[ARB], title[ARB] +real expo, expo_sum[ARB] +int beam_stat[ARB], beam, len[ARB] +bool coincidence +int ccmode +real dtime, power + +int npts, co_flag, imgeti() +pointer pix + +pointer imgl1r() + +begin + npts = IM_LEN (im, 1) + + # Map pixels and optionally correct for coincidence + pix = imgl1r (im) + if (coincidence) { + iferr (co_flag = imgeti (im, "CO-FLAG")) + co_flag = -1 + if (co_flag < 1) { + call coincor (Memr[pix], Memr[pix], npts, expo, co_flag, + dtime, power, ccmode) + } + } + + # Add in the current data + npts = min (npts, len[beam]) + + call aaddr (Memr[pix], Memr[accum[beam]], Memr[accum[beam]], npts) + + beam_stat[beam] = beam_stat[beam] + 1 + expo_sum [beam] = expo_sum [beam] + expo +end + +# WRT_FIT_SPEC -- Write out normalized spectrum + +procedure wrt_fit_spec (image, accum, expo_sum, ofile, beam, title, npts, order) + +char image[SZ_FNAME] +real accum[ARB], expo_sum +int beam, npts, order +char ofile[SZ_FNAME] +char title[SZ_LINE] + +char output[SZ_FNAME], temp[SZ_LINE] +pointer im, imnew, newpix + +pointer immap(), impl1r() +int strlen() + +begin + im = immap (image, READ_ONLY, 0) +10 call strcpy (ofile, output, SZ_FNAME) + call sprintf (output[strlen (output) + 1], SZ_FNAME, ".%04d") + call pargi (beam) + + # Create new image with a user area + # If an error occurs, ask user for another name to try + # since many open errors result from trying to overwrite an + # existing image. + + iferr (imnew = immap (output, NEW_COPY, im)) { + call eprintf ("Cannot create [%s] -- Already exists??\07\n") + call pargstr (output) + call clgstr ("output", ofile, SZ_FNAME) + go to 10 + } + + call strcpy ("Normalized flat:", temp, SZ_LINE) + call sprintf (temp[strlen (temp) + 1], SZ_LINE, "%s") + call pargstr (title) + call strcpy (temp, IM_TITLE (imnew), SZ_LINE) + IM_PIXTYPE (imnew) = TY_REAL + + newpix = impl1r (imnew) + call amovr (accum, Memr[newpix], npts) + + call imaddr (imnew, "EXPOSURE", expo_sum) + call imaddi (imnew, "QF-FLAG", order) + call imunmap (im) + call imunmap (imnew) + + call printf ("Fit for aperture %1d --> [%s]\n") + call pargi (beam) + call pargstr (output) + call flush (STDOUT) +end + +# FIT_SPEC -- Fit a line through the spectrum with user interaction + +procedure fit_spec (accum, npts, expo_sum, interact, function, + order, niter, lower, upper, ngrow, div_min, beam) + +real accum[ARB], expo_sum +bool interact +int function, order, niter, ngrow, npts, beam +real lower, upper, div_min + +int cc, key, gp, plt_mode +int i, initer, sum_niter, newgraph +real x1, y1, sigma, temp +pointer sp, wts, x, y, cv +bool first +char gtitle[SZ_LINE], command[SZ_FNAME] + +int clgcur(), clgeti() +pointer gopen() +real clgetr(), cveval() + +data plt_mode/PLT_FIT/ + +begin + # Perform initial fit + call smark (sp) + call salloc (wts, npts, TY_REAL) + call salloc (x , npts, TY_REAL) + call salloc (y , npts, TY_REAL) + + first = true + if (!interact) { + sum_niter = 0 + do i = 1, niter + call linefit (accum, npts, function, order, lower, upper, + ngrow, cv, first, Memr[wts], Memr[x]) + sum_niter = niter + + } else { + gp = gopen ("stdgraph", NEW_FILE, STDGRAPH) + call sprintf (gtitle, SZ_LINE, "Flat Field Sum - %f seconds ap:%1d") + call pargr (expo_sum) + call pargi (beam) + + key = 'r' + repeat { + switch (key) { + case 'e': # Plot errors + plt_mode = PLT_ERR + newgraph = YES + + case 'f': # Plot fit + plt_mode = PLT_FIT + newgraph = YES + + case 'o': # Change order + order = clgeti ("new_order") + # Reinstate all pixels + first = true + newgraph = YES + + case 'l': # Change lower sigma + lower = clgetr ("new_lower") + newgraph = YES + + case 'u': # Change upper sigma + upper = clgetr ("new_upper") + newgraph = YES + + case 'r': # Reset fit parameters + first = true + newgraph = YES + + case 's': # Change both rejection sigmas + lower = clgetr ("new_lower") + upper = lower + call clputr ("new_upper", upper) + newgraph = YES + + case 'i': # Iterate again - Drop thru + initer = 1 + newgraph = YES + + case 'n': # Iterate n times + initer = clgeti ("new_niter") + newgraph = YES + + case 'q': # Quit + break + + case '?': # Clear and help + call gpagefile (gp, KEY, PROMPT) + + case '/': # Status line help + call ff_sts_help + + case 'I': # Interrupt + call fatal (0, "Interrupt") + + default: + call printf ("\07\n") + } + + if (newgraph == YES) { + # Suppress an iteration if plot mode change requested + if (key != 'e' && key != 'f') { + if (first) { + sum_niter = 0 + initer = niter + call cvfree (cv) + } + do i = 1, initer + call linefit (accum, npts, function, order, lower, + upper, ngrow, cv, first, Memr[wts], Memr[x]) + sum_niter = sum_niter + initer + } + + switch (plt_mode) { + case PLT_FIT: + call plot_fit (gp, accum, cv, function, order, npts, + gtitle, Memr[wts], Memr[x], Memr[y], sigma) + case PLT_ERR: + call plot_fit_er (gp, accum, cv, function, order, npts, + gtitle, Memr[wts], Memr[x], Memr[y], sigma) + } + + newgraph = NO + } + } until (clgcur ("cursor",x1,y1,cc,key,command,SZ_FNAME) == EOF) + call gclose (gp) + } + + # Replace original data with the data/fit + do i = 1, npts { + temp = cveval (cv, real (i)) + if (temp == 0.0) + temp = max (temp, div_min) + accum[i] = accum[i] / temp + } + + call cvfree (cv) + call sfree (sp) + + # Save iteration count for next time + niter = sum_niter +end + +# LINEFIT -- Fit desired function thru data + +procedure linefit (pix, npts, function, order, lower, upper, ngrow, cv, + first, wts, x) + +real pix[ARB] # Data array to fit +int npts # Elements in array +int function # Type of fitting function +int order # Order of fitting function +real lower # Lower rejection threshold +real upper # Upper rejection threshold +int ngrow # Rejection growing radius +pointer cv +real wts[ARB] # Array weights +real x[ARB] +bool first + +int ier, i, nreject + +int reject() + +begin +10 if (first) { + do i = 1, npts { + x[i] = i + wts[i] = 1.0 + } + + # Initialize curve fitting. + call cvinit (cv, function, order, 1., real (npts)) + call cvfit (cv, x, pix, wts, npts, WTS_USER, ier) + nreject = 0 + first = false + } + + # Do pixel rejection if desired. + if ((lower > 0.) || (upper > 0.)) + nreject = reject (cv, x, pix, wts, npts, lower, upper, ngrow) + else + nreject = 0 + + if (nreject == ERR) { + call eprintf ("Cannot fit data -- too many points rejected??\n") + call cvfree (cv) + first = true + go to 10 + } +end + +# REJECT -- Reject points with large residuals from the fit. +# +# The sigma of the input to the fit is calculated. The rejection thresholds +# are set at -lower*sigma and upper*sigma. Points outside the rejection +# thresholds are rejected from the fit and flaged by setting their +# weights to zero. Finally, the remaining points are refit and a new +# fit line evaluated. The number of points rejected is returned. + +int procedure reject (cv, x, y, w, npoints, lower, upper, ngrow) + +pointer cv # Curve descriptor +real x[ARB] # Input ordinates +real y[ARB] # Input data values +real w[ARB] # Weights +int npoints # Number of input points +real lower # Lower rejection sigma +real upper # Upper rejection sigma +int ngrow # Rejection radius + +int i, j, n, i_min, i_max, nreject +real sigma, residual, resid_min, resid_max + +real cveval() + +begin + # Determine sigma of fit and set rejection limits. + sigma = 0. + n = 0 + do i = 1, npoints { + if (w[i] == 0.) + next + sigma = sigma + (y[i] - cveval (cv, x[i])) ** 2 + n = n + 1 + } + + sigma = sqrt (sigma / (n - 1)) + resid_min = -lower * sigma + resid_max = upper * sigma + + # Reject the residuals exceeding the rejection limits. + nreject = 0 + for (i = 1; i <= npoints; i = i + 1) { + if (w[i] == 0.) + next + residual = y[i] - cveval (cv, x[i]) + if ((residual < resid_min) || (residual > resid_max)) { + i_min = max (1, i - ngrow) + i_max = min (npoints, i + ngrow) + + # Reject points from the fit and flag them with zero weight. + do j = i_min, i_max { + call cvrject (cv, x[j], y[j], w[j]) + w[j] = 0. + nreject = nreject + 1 + } + i = i_max + } + } + + # Refit if points have been rejected. + if (nreject > 0) { + call cvsolve (cv, i) + if (i != OK) + return (ERR) + } + + return (nreject) +end + +# PLOT_FIT -- Plot the fit to the image line and data + +procedure plot_fit (gp, pix, cv, function, order, npts, gtitle, wts, xfit, + yfit, sigma) + +int gp, npts, function, order +real pix[ARB], wts[ARB], xfit[ARB], yfit[ARB] +pointer cv +real sigma +char gtitle[SZ_LINE] + +real x1, x2 +int i + +begin + # Set up plot + x1 = 1.0 + x2 = npts + + call gseti (gp, G_NMINOR, 0) + call gclear (gp) + call gsview (gp, 0.15, 0.95, 0.20, 0.9) + call gploto (gp, pix, npts, x1, x2, gtitle) + + # Now plot the fit + do i = 1, npts + xfit[i] = i + + call cvvector (cv, xfit, yfit, npts) + call gvline (gp, yfit, npts, x1, x2) + + # Compute sigma and write it out + call get_sigma (pix, yfit, wts, npts, sigma) + call show_status (function, order, sigma, npts, wts) +end + +# PLOT_FIT_ER -- Plot the error in the fit to the image line and data + +procedure plot_fit_er (gp, pix, cv, function, order, npts, gtitle, wts, xfit, + yfit, sigma) + +int gp, npts, function, order +real pix[ARB], wts[ARB], xfit[ARB], yfit[ARB] +pointer cv +real sigma +char gtitle[SZ_LINE] + +real x1, x2, y[2] +int i + +begin + # Set up plot + x1 = 1.0 + x2 = npts + y[1] = -0.0001 + y[2] = +0.0001 + + call cvvector (cv, xfit, yfit, npts) + + # Compute percentage errors + do i = 1, npts + if (pix[i] != 0.0) + yfit[i] = (pix[i] - yfit[i]) / pix[i] + else + yfit[i] = 0.0 + + call gseti (gp, G_NMINOR, 0) + call gclear (gp) + call gsview (gp, 0.15, 0.95, 0.20, 0.9) + + call gploto (gp, yfit, npts, x1, x2, + "Flat field fractional error in fit") + + # Draw a zero error line + call gline (gp, x1, y[1], x2, y[2]) + + # Compute sigma + call get_sigma0 (yfit, wts, npts, sigma) + call show_status (function, order, sigma, npts, wts) +end + +# SHOW_STATUS -- Show the fit status on status line + +procedure show_status (function, order, sigma, npts, wts) + +int function, order, npts +real sigma, wts[ARB] + +int i, nvals + +begin + # Count non-rejected points + nvals = 0 + do i = 1, npts + if (wts[i] != 0.0) + nvals = nvals + 1 + + call printf ("Fit type: %s order: %2d rms: %6.3f") + switch (function) { + case LEGENDRE: + call pargstr ("Legendre") + case CHEBYSHEV: + call pargstr ("Chebyshev") + case SPLINE3: + call pargstr ("Spline3") + case SPLINE1: + call pargstr ("Spline1") + default: + call pargstr ("???") + } + + call pargi (order) + call pargr (sigma) + + call printf (" points: %d out of %d") + call pargi (nvals) + call pargi (npts) + + call flush (STDOUT) +end + +# GET_SIGMA -- Compute rms error between two vectors whose average difference +# is zero. + +procedure get_sigma (y1, y2, wts, n, sigma) + +real y1[ARB], y2[ARB], wts[ARB], sigma +int n + +int i, nval +real sum + +begin + sum = 0.0 + nval = 0 + do i = 1, n + if (wts[i] != 0.0) { + sum = sum + (y1[i] - y2[i]) ** 2 + nval = nval + 1 + } + + sigma = sqrt (sum / (nval-1)) + return +end + +# GET_SIGMA0 -- Compute rms error of a vector + +procedure get_sigma0 (y1, wts, n, sigma) + +real y1[ARB], wts[ARB], sigma +int n + +int i, nval +real sum + +begin + sum = 0.0 + nval = 0 + do i = 1, n + if (wts[i] != 0.0) { + sum = sum + y1[i]**2 + nval = nval + 1 + } + + sigma = sqrt (sum / (nval-1)) + return +end + +# FF_STS_HELP -- Status line help for Flat Fit + +procedure ff_sts_help () + +int linenr, maxline + +data linenr/1/ +data maxline/2/ + +begin + switch (linenr) { + case 1: + call printf ("e=err plot f=data plot o=order l=lower sigma ") + call printf ("u=upper sigma s=both sigmas") + + case 2: + call printf ("r=incl reject i=iterate n=niterate q=quit ") + call printf ("?=help /=linehelp <CR>=quit") + } + + call flush (STDOUT) + + linenr = linenr + 1 + if (linenr > maxline) + linenr = 1 +end diff --git a/noao/onedspec/irsiids/t_slist1d.x b/noao/onedspec/irsiids/t_slist1d.x new file mode 100644 index 00000000..75837d50 --- /dev/null +++ b/noao/onedspec/irsiids/t_slist1d.x @@ -0,0 +1,163 @@ +include <error.h> +include <imhdr.h> +include <fset.h> +include <smw.h> + + +# SLIST1D -- Lists header information from IIDS/IRS format header +# This is the original T_SLIST. + +procedure t_slist1d () + +int root +int long_header +pointer sp, image, im, mw, sh, ptr +int i, nl, df, sm, qf, qd, bs, co + +int btoi(), imtgetim(), imgeti() +bool clgetb() +pointer imtopenp(), immap(), smw_openim() +errchk immap, smw_openim, shdr_open + +begin + call smark (sp) + call salloc (image, SZ_LINE, TY_CHAR) + + # Parameters + root = imtopenp ("input") + call clgstr ("records", Memc[image], SZ_LINE) + call odr_openp (root, Memc[image]) + long_header = btoi (clgetb ("long_header")) + + # Initialize + call fseti (STDOUT, F_FLUSHNL, YES) + + # Loop over all input images by subsets + while (imtgetim (root, Memc[image], SZ_FNAME) != EOF) { + + # Open image + iferr { + im = NULL + mw = NULL + ptr = immap (Memc[image], READ_ONLY, 0); im = ptr + ptr = smw_openim (im); mw = ptr + call shdr_open (im, mw, 1, 1, INDEFI, SHHDR, sh) + } then { + if (mw != NULL) { + call smw_close (mw) + if (sh != NULL) + MW(sh) = NULL + } + if (im != NULL) + call imunmap (im) + call erract (EA_WARN) + next + } + + nl = IM_LEN(im,2) + do i = 1, nl { + call shdr_open (im, mw, i, 1, INDEFI, SHHDR, sh) + + if (long_header == YES) { + call printf ("[%s] %4dpts %s\n") + call pargstr (IMNAME(sh)) + call pargi (SN(sh)) + call pargstr (TITLE(sh)) + + if (OFLAG(sh) == 1) { + call printf ("oflag = OBJECT, beam_number = %d") + call pargi (BEAM(sh)) + } else if (OFLAG (sh) == 0) { + call printf ("oflag = SKY, beam_number = %d") + call pargi (BEAM(sh)) + } + call printf (",\n") + + iferr (df = imgeti (im, "DF-FLAG")) + df = -1 + iferr (sm = imgeti (im, "SM-FLAG")) + sm = -1 + iferr (qf = imgeti (im, "QF-FLAG")) + qf = -1 + iferr (qd = imgeti (im, "QD-FLAG")) + qd = -1 + iferr (bs = imgeti (im, "BS-FLAG")) + bs = -1 + iferr (co = imgeti (im, "CO-FLAG")) + co = -1 + + # Airmass may not be in header. It could be computed if + # if the observatory latitude were available. + + call printf ("airmass = %5.3f,%25tW0 = %0.3f,") + call pargr (AM(sh)) + call pargr (W0(sh)) + call printf (" WPC = %0.5g, ITM = %.2f,\n") + call pargr (WP(sh)) + call pargr (IT(sh)) + call printf ("NP1 = %d, NP2 = %d,") + call pargi (NP1(sh)) + call pargi (NP2(sh)) + call printf (" UT = %0.1h, ST = %0.1h,\n") + call pargr (UT(sh)) + call pargr (ST(sh)) + call printf ("HA = %0.2h,") + call pargr (HA(sh)) + call printf (" RA = %0.2h, DEC = %0.1h,\n") + call pargr (RA(sh)) + call pargr (DEC(sh)) + call printf ( + "df = %d, sm = %d, qf = %d, dc = %d, qd = %d, ") + call pargi (df) + call pargi (sm) + call pargi (qf) + call pargi (DC(sh)) + call pargi (qd) + call printf ("ex = %d, bs = %d, ca = %d, co = %d") + call pargi (EC(sh)) + call pargi (bs) + call pargi (FC(sh)) + call pargi (co) + + call printf ("\n\n") + } else { + if (nl == 1) { + call printf ("[%s]:%s %4ds %4dpts %s\n") + call pargstr (IMNAME(sh)) + if (OFLAG(sh) == 1) + call pargstr ("o") + else + call pargstr ("s") + call pargr (IT(sh)) + call pargi (SN(sh)) + call pargstr (TITLE(sh)) + } else { + call printf ("[%s]:%s %6.2fs %4dpts %dspectra %s\n") + call pargstr (IMNAME(sh)) + if (OFLAG(sh) == 1) + call pargstr ("o") + else + call pargstr ("s") + call pargr (IT(sh)) + call pargi (SN(sh)) + call pargi (nl) + call pargstr (TITLE(sh)) + } + break + } + } + + call smw_close (mw) + if (sh != NULL) + MW(sh) = NULL + call imunmap (im) + } + + # Null out record string to avoid learn mode + call clpstr ("records", "") + + # Free space + call shdr_close (sh) + call imtclose (root) + call sfree (sp) +end diff --git a/noao/onedspec/irsiids/t_subsets.x b/noao/onedspec/irsiids/t_subsets.x new file mode 100644 index 00000000..6a2a61bf --- /dev/null +++ b/noao/onedspec/irsiids/t_subsets.x @@ -0,0 +1,121 @@ +include <error.h> +include <imhdr.h> + + +# T_SUBSETS -- Sub a series of spectra by pairs. A single spectrum +# is produced for every pair. +# + +procedure t_subsets () + +pointer image +pointer recstr, ofile +int root, start_rec, subset +int nrecs +int npts, nrem, ifile, tog +real expo, wtsum +pointer sp, recs, im[2], cur_pix, sp_sum + +real imgetr() +int clpopni(), clgeti() +int get_next_image(), decode_ranges() +pointer immap(), imgl1r() + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (ofile, SZ_FNAME, TY_CHAR) + call salloc (recstr, SZ_LINE, TY_CHAR) + call salloc (recs, 300, TY_INT) + + # Open input file name template + root = clpopni ("input") + + # Get range specification if any + call clgstr ("records", Memc[recstr], SZ_LINE) + if (decode_ranges (Memc[recstr], Memi[recs], 100, nrecs) == ERR) + call error (0, "Bad range specification") + + # Get rootname for output files and starting record + call clgstr ("output", Memc[ofile], SZ_FNAME) + start_rec = clgeti ("start_rec") + + # Initialize range decoder + call reset_next_image () + + #Initialize file counter + ifile = 0 + + # Set weighting value needed by spectrum writer + wtsum = 1.0 + + # Define subset of operation is a pair + subset = 2 + + # Loop over all input images by subsets + while (get_next_image (root, Memi[recs], nrecs, Memc[image], + SZ_FNAME) != EOF) { + + # Get toggle value + tog = mod (ifile, 2) + 1 + + # Open image + iferr (im[tog] = immap (Memc[image], READ_ONLY, 0)) { + call erract (EA_WARN) + next + } + + # Load data + cur_pix = imgl1r (im[tog]) + + # Allocate space for the sum + if (mod (ifile,2) == 0) { + npts = IM_LEN (im[tog],1) + call calloc (sp_sum, npts, TY_REAL) + + # Zero exposure counter + expo = 0.0 + + # Add first spectrum + call amovr (Memr[cur_pix], Memr[sp_sum], npts) + + iferr (expo = imgetr (im[tog], "EXPOSURE")) + iferr (expo = imgetr (im[tog], "ITIME")) + iferr (expo = imgetr (im[tog], "EXPTIME")) + expo = 1 + + call printf ("[%s] added\n") + call pargstr (Memc[image]) + call flush (STDOUT) + + } else { + # Subtract second spectrum + call asubr (Memr[sp_sum], Memr[cur_pix], Memr[sp_sum], + min (npts, IM_LEN(im[tog],1))) + call printf ("[%s] subtracted\n") + call pargstr (Memc[image]) + call flush (STDOUT) + call imunmap (im[2]) + + call wrt_set (Memr[sp_sum], subset, im[1], Memc[ofile], + start_rec, expo, wtsum, -1) + call mfree (sp_sum, TY_REAL) + } + + ifile = ifile + 1 + } + # Check that there are no remaining spectra in an unfulfilled subset + nrem = mod (ifile, 2) + if (nrem != 0) { + call mfree (sp_sum, TY_REAL) + + call eprintf ("Unfulfilled pair ignored\n") + } + + # Update record number + call clputi ("next_rec", start_rec) + + # Free space + call sfree (sp) + call clpcls (root) +end diff --git a/noao/onedspec/irsiids/t_sums.x b/noao/onedspec/irsiids/t_sums.x new file mode 100644 index 00000000..e28ebb35 --- /dev/null +++ b/noao/onedspec/irsiids/t_sums.x @@ -0,0 +1,239 @@ +include <error.h> +include <imhdr.h> + +define MAX_NR_BEAMS 100 # Max number of instrument apertures + +# T_SUMS -- Compute sums of strings of spectra according to +# Aperture number and object/sky flag. So for IIDS/IRS +# type spectra, 4 sums will be generated. +# In general, there will be 2N sums where N is the number +# apertures. + +procedure t_sums () + +pointer image # Image name to be added +pointer images # Image name to be added +pointer ofile # Output image file name +pointer recstr # Record number string +int recs # Spectral record numbers +int root, nrecs # CL and ranges flags +real expo # Exposure time +pointer bstat[2] # Status of each aperture +pointer npts[2] # Length of spectrum +pointer esum[2] # Accumulated exposure time +pointer accum[2] # Pointers to beam accumulators +pointer title[2] +int beam, object +int start_rec + +int i, j +pointer sp, work, im + +real imgetr() +int clgeti(), clpopni(), imgeti() +int get_next_image(), decode_ranges() +pointer immap() + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (images, MAX_NR_BEAMS, TY_POINTER) + call salloc (ofile, SZ_FNAME, TY_CHAR) + call salloc (recstr, SZ_LINE, TY_CHAR) + call salloc (recs, 300, TY_INT) + call salloc (accum, MAX_NR_BEAMS, TY_POINTER) + call salloc (title, MAX_NR_BEAMS, TY_POINTER) + call amovki (NULL, Memi[images], MAX_NR_BEAMS) + call salloc (work, 2*5*MAX_NR_BEAMS, TY_STRUCT) + bstat[1] = work + bstat[2] = work + MAX_NR_BEAMS + npts[1] = work + 2 * MAX_NR_BEAMS + npts[2] = work + 3 * MAX_NR_BEAMS + esum[1] = work + 4 * MAX_NR_BEAMS + esum[2] = work + 5 * MAX_NR_BEAMS + accum[1] = work + 6 * MAX_NR_BEAMS + accum[2] = work + 7 * MAX_NR_BEAMS + title[1] = work + 8 * MAX_NR_BEAMS + title[2] = work + 9 * MAX_NR_BEAMS + + # Get task parameters. + root = clpopni ("input") + + # Get input record numbers + call clgstr ("records", Memc[recstr], SZ_LINE) + if (decode_ranges (Memc[recstr], Memi[recs], 100, nrecs) == ERR) + call error (0, "Bad range specification") + + call clgstr ("output", Memc[ofile], SZ_LINE) + + start_rec = clgeti ("start_rec") + + call reset_next_image () + + # Clear all beam status flags + call amovki (INDEFI, Memi[bstat[1]], MAX_NR_BEAMS*2) + call aclrr (Memr[esum[1]], MAX_NR_BEAMS*2) + + call printf ("Accumulating spectra --\n") + call flush (STDOUT) + + while (get_next_image (root, Memi[recs], nrecs, Memc[image], + SZ_FNAME) != EOF) { + iferr (im = immap (Memc[image], READ_ONLY, 0)) { + call erract (EA_WARN) + next + } + + # Load header + iferr (beam = imgeti (im, "BEAM-NUM")) + beam = 0 + if (beam < 0 || beam > MAX_NR_BEAMS-1) + call error (0, "Invalid aperture number") + + # Select array: Object = array 2; sky = array 1 + iferr (object = imgeti (im, "OFLAG")) + object = 1 + if (object == 1) + object = 2 + else + object = 1 + + iferr (expo = imgetr (im, "EXPOSURE")) + iferr (expo = imgetr (im, "ITIME")) + iferr (expo = imgetr (im, "EXPTIME")) + expo = 1 + + # Add spectrum into accumulator + if (IS_INDEFI (Memi[bstat[object]+beam])) { + Memi[npts[object]+beam] = IM_LEN (im,1) + call salloc (Memi[accum[object]+beam], IM_LEN(im,1), TY_REAL) + call aclrr (Memr[Memi[accum[object]+beam]], IM_LEN(im,1)) + Memi[bstat[object]+beam] = 0 + + call salloc (Memi[title[object]+beam], SZ_LINE, TY_CHAR) + call strcpy (IM_TITLE(im), Memc[Memi[title[object]+beam]], + SZ_LINE) + } + + call su_accum_spec (im, Memi[npts[1]], expo, Memi[bstat[1]], + beam+1, Memi[accum[1]], Memr[esum[1]], Memi[title[1]], object) + + call printf ("[%s] %s spectrum added to aperture %1d\n") + call pargstr (Memc[image]) + if (object == 2) + call pargstr ("object") + else + call pargstr ("sky ") + call pargi (beam) + call flush (STDOUT) + + if (Memi[images+beam] == NULL) + call salloc (Memi[images+beam], SZ_FNAME, TY_CHAR) + call strcpy (Memc[image], Memc[Memi[images+beam]], SZ_FNAME) + call imunmap (im) + } + + # Review all apertures containing data and write sums + do i = 0, MAX_NR_BEAMS-1 + do j = 1, 2 + if (!IS_INDEFI (Memi[bstat[j]+i])) { + call wrt_spec (Memc[Memi[images+i]], Memr[Memi[accum[j]+i]], + Memr[esum[j]+i], Memc[ofile], start_rec, + Memc[Memi[title[j]+i]], Memi[npts[j]+i], i, j) + + start_rec = start_rec + 1 + } + + call clputi ("next_rec", start_rec) + call sfree (sp) + call clpcls (root) +end + +# ACCUM_SPEC -- Accumulate spectra by beams + +procedure su_accum_spec (im, len, expo, beam_stat, beam, accum, expo_sum, + title, object) + +pointer im, accum[MAX_NR_BEAMS,2], title[MAX_NR_BEAMS,2] +real expo, expo_sum[MAX_NR_BEAMS,2] +int beam_stat[MAX_NR_BEAMS,2], beam, len[MAX_NR_BEAMS,2] +int object + +int npts +pointer pix + +pointer imgl1r() + +begin + npts = IM_LEN (im, 1) + + # Map pixels and optionally correct for coincidence + pix = imgl1r (im) + + # Add in the current data + npts = min (npts, len[beam, object]) + + call aaddr (Memr[pix], Memr[accum[beam, object]], + Memr[accum[beam, object]], npts) + + beam_stat[beam, object] = beam_stat[beam, object] + 1 + expo_sum [beam, object] = expo_sum [beam, object] + expo +end + +# WRT_SPEC -- Write out normalized spectrum + +procedure wrt_spec (image, accum, expo_sum, ofile, start, title, npts, object, + beam) + +char image[SZ_FNAME] +real accum[ARB], expo_sum +int start, npts +char ofile[SZ_FNAME] +char title[SZ_LINE] +int object, beam + +char output[SZ_FNAME], temp[SZ_LINE] +pointer im, imnew, newpix + +pointer immap(), impl1r() +int strlen() + +begin + im = immap (image, READ_ONLY, 0) +10 call strcpy (ofile, output, SZ_FNAME) + call sprintf (output[strlen (output) + 1], SZ_FNAME, ".%04d") + call pargi (start) + + # Create new image with a user area + # If an error occurs, ask user for another name to try + # since many open errors result from trying to overwrite an + # existing image. + + iferr (imnew = immap (output, NEW_COPY, im)) { + call eprintf ("Cannot create [%s] -- Already exists??\07\n") + call pargstr (output) + call clgstr ("newoutput", ofile, SZ_FNAME) + go to 10 + } + + call strcpy ("Summation:", temp, SZ_LINE) + call sprintf (temp[strlen (temp) + 1], SZ_LINE, "%s") + call pargstr (title) + call strcpy (temp, IM_TITLE (imnew), SZ_LINE) + + newpix = impl1r (imnew) + call amovr (accum, Memr[newpix], npts) + + call imaddr (imnew, "EXPOSURE", expo_sum) + call imunmap (im) + call imunmap (imnew) + + call printf ("%s sum for aperture %1d --> [%s]\n") + if (object == 1) + call pargstr ("Object") + else + call pargstr ("Sky ") + call pargi (beam) + call pargstr (output) + call flush (STDOUT) +end diff --git a/noao/onedspec/irsiids/t_widstape.x b/noao/onedspec/irsiids/t_widstape.x new file mode 100644 index 00000000..1f96d146 --- /dev/null +++ b/noao/onedspec/irsiids/t_widstape.x @@ -0,0 +1,343 @@ +include <mach.h> +include <error.h> +include <imhdr.h> +include <smw.h> + +define SZ_IDSTITLE 64 # Length of IDSOUT title +define SZ_CARD 80 # Columns on a card + +# T_WIDSTAPE -- Convert each line of an IRAF image to IDSOUT text format. +# Each image line is treated as a one dimensional spectrum. +# A maximum IDSOUT length of 1024 points is enforced silently. +# +# There are two types of output: +# single -- All image lines are appended to a single IDSOUT file. +# multiple -- Each image line is appended to a different IDSOUT file. + +procedure t_widstape () + +pointer image # Image to be converted +pointer recs # Record numbers +pointer idsout # IDSOUT file or root name to be written +int block_size # Block size +bool ebcdic # ASCII or EBCDIC + +int i, mfd, root, nrecs +pointer sp, im, mw, sh, ptr + +int open(), mtopen(), clgeti(), clpopni() +int get_next_image(), decode_ranges(), mtfile(), mtneedfileno() +bool clgetb() +pointer immap(), smw_openim() +errchk immap, smw_openim, shdr_open, wrt_ids_rec + +begin + call smark (sp) + call salloc (image, SZ_LINE, TY_CHAR) + call salloc (idsout, SZ_FNAME, TY_CHAR) + call salloc (recs, 300, TY_INT) + + # Parameters + root = clpopni ("input") + call clgstr ("records", Memc[image], SZ_LINE) + call clgstr ("idsout", Memc[idsout], SZ_FNAME) + block_size = clgeti ("block_size") + ebcdic = clgetb ("ebcdic") + + # Set record numbers + if (decode_ranges (Memc[image], Memi[recs], 100, nrecs) == ERR) + call error (0, "Bad range specification") + + # Check that a realistic block size was requested + if (mod (block_size, SZ_CARD) == 0) + block_size = block_size / SZB_CHAR + else + call error (0, "Blocks not integral number of cards") + + # Open output tape file + # First determine if a file number was specified + if (mtfile (Memc[idsout]) == YES) { + + # If no file, check if new_tape was specified and if so, + # force file=1; otherwise force file=EOT + + if (mtneedfileno (Memc[idsout]) == YES) { + if (!clgetb("new_tape")) + call mtfname (Memc[idsout], EOT, Memc[idsout], SZ_FNAME) + + else + call mtfname (Memc[idsout], 1, Memc[idsout], SZ_FNAME) + } + mfd = mtopen (Memc[idsout], WRITE_ONLY, block_size) + } else + mfd = open (Memc[idsout], NEW_FILE, BINARY_FILE) + + # Loop over all files + call reset_next_image () + while (get_next_image (root, Memi[recs], nrecs, Memc[image], + SZ_LINE) != EOF) { + iferr { + im = NULL + mw = NULL + ptr = immap (Memc[image], READ_ONLY, 0); im = ptr + ptr = smw_openim (im); mw = ptr + + # Write out a spectrum for each line in the image + do i = 1, IM_LEN (im,2) { + call shdr_open (im, mw, i, 1, INDEFI, SHDATA, sh) + call wrt_ids_rec (mfd, sh, Memc[image], ebcdic) + } + + call printf ("copied - [%s]: %s\n") + call pargstr (IMNAME(sh)) + call pargstr (TITLE(sh)) + call flush (STDOUT) + } then + call erract (EA_WARN) + + if (mw != NULL) + call smw_close (mw) + if (im != NULL) + call imunmap (im) + } + + call shdr_close (sh) + call close (mfd) + call sfree (sp) +end + + +# WRT_IDS_REC -- Write one IIDS/IRS format record in IDSOUT form + +procedure wrt_ids_rec (mfd, sh, image, ebcdic) + +int mfd +pointer sh +char image[SZ_FNAME] +bool ebcdic + +# IDSOUT header parameters +char label[SZ_IDSTITLE] # Record label +int record # Record number +int uttime # UT time in seconds +int st # Siderial time in seconds +real ra # Right Ascension in hours +real dec # Declination in degrees +real ha # Hour angle in hours +real airmass # Air mass +int itime # Integration time +real wavelen1 # Wavelength of first pixel +real dispersion # Dispersion per pixel + +int i, rec_no, df, sm, qf, qd, bs, co +pointer sp, padline, bufline + +int strmatch(), imgeti() + +begin + call smark (sp) + call salloc (padline, SZ_LINE, TY_CHAR) + call salloc (bufline, SZ_LINE, TY_CHAR) + + # Fill in header parameters. + + call strcpy (TITLE(sh), label, SZ_IDSTITLE) + + # The following two calculations were causing floating overflows + # when the header values were indefinite. SEH 7-23-86 + if (IS_INDEF(UT(sh))) + uttime = INDEFI + else + uttime = UT(sh) * 3600. + + if (IS_INDEF(ST(sh))) + st = INDEFI + else + st = ST(sh) * 3600. + + ra = RA(sh) + dec = DEC(sh) + ha = HA(sh) + airmass = AM(sh) + itime = IT(sh) + wavelen1 = W0(sh) + dispersion = WP(sh) + + iferr (df = imgeti (IM(sh), "DF-FLAG")) + df = -1 + iferr (sm = imgeti (IM(sh), "SM-FLAG")) + sm = -1 + iferr (qf = imgeti (IM(sh), "QF-FLAG")) + qf = -1 + iferr (qd = imgeti (IM(sh), "QD-FLAG")) + qd = -1 + iferr (bs = imgeti (IM(sh), "BS-FLAG")) + bs = -1 + iferr (co = imgeti (IM(sh), "CO-FLAG")) + co = -1 + + # Create a padding line to fill the IDSOUT block to 1024 points. + + call sprintf (Memc[padline], SZ_LINE, + "%10.4e%10.4e%10.4e%10.4e%10.4e%10.4e%10.4e%10.4e\n") + do i = 1, 8 + call pargr (0.) + + # Line 1 -- Record number, etc. + rec_no = strmatch (image, ".") + call sscan (image[rec_no]) + call gargi (record) + + call sprintf (Memc[bufline], SZ_LINE, + "%5d%5d%15.7e%15.7e%5d%5d%5d%5d%5d%5d%10d") + call pargi (record) + call pargi (itime) + call pargr (wavelen1) + call pargr (dispersion) + call pargi (0) + call pargi (SN(sh)) + call pargi (BEAM(sh)) + call pargi (-1) + call pargi (-1) + call pargi (0) + call pargi (uttime) + + call putcard (mfd, Memc[bufline], ebcdic) + + # Line 2 -- Siderial time, RA, and Dec. + + call sprintf (Memc[bufline], SZ_LINE, + "%10d%15.7e%15.7e%5d%5d%5d%5d%5d%5d%5d%5d") + call pargi (st) + call pargr (ra) + call pargr (dec) + call pargi (0) + call pargi (df) + call pargi (sm) + call pargi (qf) + call pargi (DC(sh)) + call pargi (qd) + call pargi (EC(sh)) + call pargi (bs) + + call putcard (mfd, Memc[bufline], ebcdic) + + # Line 3 -- Hour angle, air mass, UT date, and exposure title. + + call sprintf (Memc[bufline], SZ_LINE, + "%5d%5d%2w%-3.3s%5d%15.7e%15.7e%27wEND") + call pargi (FC(sh)) + call pargi (co) + call pargstr ("IRF") + call pargi (OFLAG(sh)) + call pargr (ha) + call pargr (airmass) + + call putcard (mfd, Memc[bufline], ebcdic) + + # Line 4 -- Record label. + call sprintf (Memc[bufline], SZ_LINE, "%-77sEND") + call pargstr (TITLE(sh)) + + call putcard (mfd, Memc[bufline], ebcdic) + + # Lines 5 to 132 + + call putdata (mfd, Memr[SY(sh)], SN(sh), Memc[padline], + Memc[bufline], ebcdic) + + # Line 133 -- Blank line + + call sprintf (Memc[bufline], SZ_LINE, "%80w") + call putcard (mfd, Memc[bufline], ebcdic) +end + + +# PUTDATA -- Format and output extraction data to IDSOUT length of 1024 points. +# Special effort is made to make the zero padding efficient. + +procedure putdata (mfd, data, npts, padline, bufline, ebcdic) + +int mfd # IDSOUT file descriptor +real data[npts] # Data +int npts # Number of data points +char padline[ARB] # Padding string +char bufline[ARB] # Output buffer string +bool ebcdic # Convert to ebcdic + +int i, j, k, l, n +int index +double ddata + +int dtoc3() + +begin + j = min (1024, npts) # Maximum number of data points + k = j / 8 * 8 # Index of last data point in last complete line + if (k < j) + l = k + 8 # Index of last point in last line with data + else + l = k + + # Write all complete data lines. + + index = 1 + do i = 1, k { + ddata = double (data[i]) + n = dtoc3 (ddata, bufline[index], 10, 4, 'e', 10) + while (n < 10) { + bufline[index+n] = ' ' + n = n + 1 + } + index = index + 10 + if (mod (i, 8) == 0) { + call putcard (mfd, bufline, ebcdic) + index = 1 + } + } + + # Write partial data line. + + index = 1 + do i = k + 1, l { + if (i <= j) { + ddata = double (data[i]) + n = dtoc3 (ddata, bufline[index], 11, 5, 'e', 10) + } else + n = dtoc3 (0.D0, bufline[index], 11, 5, 'e', 10) + while (n < 10) { + bufline[index+n] = ' ' + n = n + 1 + } + index = index + 10 + if (mod (i, 8) == 0) { + call putcard (mfd, bufline, ebcdic) + index = 1 + } + } + + # Write remaining padding lines. + + do i = l + 1, 1024, 8 + call putcard (mfd, padline, ebcdic) +end + +# PUTCARD -- Convert to ebcdic if desired and write out card + +procedure putcard (mfd, bufline, ebcdic) + +int mfd +char bufline[ARB] +bool ebcdic + +char packline[SZ_LINE] + +begin + if (ebcdic) { + call ascii_to_ebcdic (bufline, packline, SZ_CARD) + call achtsb (packline, packline, SZ_CARD) + } else + call chrpak (bufline, 1, packline, 1, SZ_CARD) + + call write (mfd, packline, SZ_CARD/SZB_CHAR) +end diff --git a/noao/onedspec/irsiids/widstape.par b/noao/onedspec/irsiids/widstape.par new file mode 100644 index 00000000..33dee906 --- /dev/null +++ b/noao/onedspec/irsiids/widstape.par @@ -0,0 +1,8 @@ +# IDSOUT parameter file -- write a CYBER style IDSOUT tape + +idsout,s,a,,,,Output file or magtape +input,s,a,,,,Image root name to write +records,s,a,,,,Records to write +block_size,i,h,3200,80,10640,Tape block size in bytes +new_tape,b,h,no,,,Is this a new (blank) tape +ebcdic,b,h,no,,,Convert character code to ebcdic |