procedure lickbase (input, output)
#########################################################################
# This is a procedure to baseline data taken at Mt. Hamilton            #
# using the default condition that the baseline is in the last column.  #
# The bias subtracted bias column is retained in the final image.       #
#                                                                       #
# original code was lickbl by Hank Donnelly long, long ago              #
# Modified from lickbl by Nan Ellman February 6, 1990 and called defbl  #
#                                                                       #
# Hacked ungloriously by Steve Allen in 1993 August to add all sorts    #
# of extra cruft as might be desired by someone running the task(s)     #
# in noao.imred.ccdred, especially ccdproc                              #
#                                                                       #
# The gain can be set to a different value if the image is binned       #
# by making use of the parameter gainb.  This is relevant for Lick      #
# images where the amplifier gain automatically reduces by a factor     #
# of 4 whenever any binning is done.  (Note that this means the IRAF    #
# GAIN parameter must increase by a factor of 4, because IRAF really    #
# the inverse gain in units of [e-/DN].)                                #
#                                                                       #
# If the parameters gainf, rdnoisef and/or gainbf are set then these    #
# values will be used if the CLKSPEED card has a value of 'Fast'.       #
#                                                                       #
# now can be run with output images replacing input images              #
#########################################################################

string  input                   {prompt = 'Input images'}
string  output                  {prompt = 'Output images'}
real    rdnoise=-1.             {prompt = "(Slow) Readout Noise (if >= 0)"}
real    rdnoisef=-1.            {prompt = "Fast Readout Noise   (if >= 0)"}
real    gain=-1.                {prompt = "(Slow) Gain          (if >  0)"}
real    gainf=-1.               {prompt = "Fast Gain            (if >  0)"}
real    gainb=-1.               {prompt = "(Slow) Binned Gain   (if >  0)"}
real    gainbf=-1.              {prompt = "Fast Binned Gain     (if >  0)"}
string  pixtrim=""              {prompt = "Trim pixels (physical pixels)"}
int     dispaxis=-1             {prompt = "Dispersion Axis (if >  0)"}
string  observat=""             {prompt = "Observatory database key"}

				# file structure pointers for reading lists
struct  *ilist                  {prompt = "Just leave this blank"}
struct  *olist                  {prompt = "Just leave this blank"}

begin
    string  l_input, l_output   # internal copies of input and output
    string  ifile, ofile        # temporary files containing expanded lists
    int     nin, nout           # number of entries in expanded lists
    string  in, out             # loop variables for stepping thru file lists
    string  biasframe           # temporary image file holding the bias column
    string  basesect            # the last column of the input image
    string  wholsec             # the whole input image, actually used to trim
    string  trimsec             # the trim section, not actually used to trim
    string  datasec             # the data section
    string  ccdsec              # readout window expressed a la IRAF
    string  lickpix             # next best thing to IRAF origsec
    string  ccdsum              # the on chip summation
    string  overstr             # value of the overscan card on output
    int     naxis1, naxis2      # size of the input image
    int     cdelt1, cdelt2      # binning or on chip summation
    int     crval1, crval2      # offset of readout window from chip origin
    string  s1,s2,s3,s4,s5,s6   # pieces of the date
    bool    binned              # is the image binned
    real    b1,b2               # binning factors
    int     ib1,ib2             # binning factors
    string  clkspeed            # fast or slow readout mode
    int     sp1,ep1,sp2,ep2     #
    bool    fastmode            # true if Lick fast readout mode

    print("Lickbase")
#   we make use of the noao.imred.bias task as a part of this script
    noao
	imred
	    bias

#   fix so script will work in the background
    cache ("sections")
    cache ("imgets")

#   Pick names for the temporary files
#   Note that your variable tmp$ had better point to a big enough place
    ifile = mktemp ("tmp$lickbl")
    ofile = mktemp ("tmp$lickbl")
#   biasframe might have a name collision because of the ".imh"
    biasframe = mktemp ("tmp$lickbl")

#   prompt the user for the query parameters
    l_input = input
    l_output = output

#   expand the input and output image templates
#   (i.e., the given names may have contained wildcard characters)
    sections (l_input, > ifile)
    nin = sections.nimages
    sections (l_output, > ofile)
    nout = sections.nimages

    if (nin != nout) {
#       this is an error condition, we need number in equal to number out
#       as always, take care to delete our temporary files
	delete (ifile // "," // ofile, ver-, >& "dev$null")
	error (1, "Input and output lists don't match")
    }


#   set file pointers to look at our input and output lists, and loop
    ilist = ifile; olist = ofile
    while (fscan (ilist, in) != EOF && fscan (olist, out) != EOF) {

#       be sure that the baseline correction has not been done before now
	s1 = ""
	hselect( in, "OVERSCAN", yes) | scan (s1)
	if (s1 != "") {
	    print( in//":  already baselined?  No action taken.")
#           go on to the next image in the while loop
	    next
	}

#       read the header cards to get the image size
	imgets(in, "i_naxis1", value="")
	naxis1=int(imgets.value)
	imgets(in, "i_naxis2", value="")
	naxis2=int(imgets.value)
#       read the header cards to get the binning (or on chip summation)
	imgets(in, "cdelt1", value="")
	cdelt1=int(imgets.value)
	imgets(in, "cdelt2", value="")
	cdelt2=int(imgets.value)
#       read the header cards to get the window offset from the CCD origin
	imgets(in, "crval1", value="")
	crval1=int(imgets.value)
	imgets(in, "crval2", value="")
	crval2=int(imgets.value)

#       create the strings describing the image sections
#       In this script we will not trim the image, ...
	wholsec='[1:'+str(naxis1)+',1:'+str(naxis2)+']'
#       ... but we will write out the trim section for ccdproc
	# simple trimsec which is as big as the actual imaging area
	# this is the default which may be replaced if pixtrim is consistent
	trimsec='[1:'+str(naxis1-1)+',1:'+str(naxis2)+']'
	# now proceed to see if we change based on pixtrim
	if (pixtrim != "") {
	    # convert from 1-based physical pixels on the original detector to
	    # IRAF-style pixels matching the binning and readout window
	    parsec(pixtrim)
	    #sp1 = max(1       , (parsec.sp[1]-crval1+cdelt1-1)/cdelt1 + 1 )
	    #ep1 = min(naxis1-1, (parsec.ep[1]-crval1+1)       /cdelt1      )
	    #sp2 = max(1       , (parsec.sp[2]-crval2+cdelt2-1)/cdelt2 + 1 )
	    #ep2 = min(naxis2  , (parsec.ep[2]-crval2+1)       /cdelt2      )
#       print("debug crvaldelt ",crval1,cdelt1,crval2,cdelt2)
#       print("debug parsec ",parsec.sp[1],parsec.ep[1],parsec.sp[2],parsec.ep[2])
	    sp1 = max(1       , (parsec.sp[1]-crval1+cdelt1-2)/cdelt1 + 1 )
	    ep1 = min(naxis1-1, (parsec.ep[1]-crval1  )       /cdelt1     )
	    sp2 = max(1       , (parsec.sp[2]-crval2+cdelt2-2)/cdelt2 + 1 )
	    ep2 = min(naxis2  , (parsec.ep[2]-crval2  )       /cdelt2     )
#       print("debug sep1, sep2 ", sp1, ep1, sp2, ep2   )
	    if ( ep1 >= sp1 && ep2 >= sp2 ) {
		# actually reset the trimsec
		trimsec='['+str(sp1)+':'+str(ep1)+','+str(sp2)+':'+str(ep2)+']'
	    } else {
		# do not reset trimsec and warn user
		print( in//":  pixtrim results in null trimsec, no trimming done.")
	    }
	}
#       this would and could be called biassect if we were from Tucson
	basesect= '['+str(naxis1)+':'+str(naxis1)+',1:'+str(naxis2)+']'
#       readout window as would have been given by a NOAO data system
	datasec='['+str(1)+':'+str(naxis1-1)+','+str(1)+':'+str(naxis2)+']'
#       WARNING:  It is not obvious what CCDSEC should be when binned.
#       We adopt here the assumption that a binned detector is an entirely
#       different detector which has fewer pixels.  (This seems to be most
#       consistent with the IRAF way of viewing the world.)
#       If in the Lick header ( CRVALn % CDELTn == 0 ), i.e. there
#       is an integral number of binned pixels preceding the readout
#       window, then this definition of CCDSEC makes sense.  If the
#       above is not true, then all bets are off because the result
#       must depend on which detector, how many amplifiers, and other
#       factors which have never been specified for Lick CCD readouts.
	ccdsec='['+str(1+(crval1/cdelt1))+':'+str(naxis1+(crval1/cdelt1)-1)+','+str(1+(crval2/cdelt2))+':'+str(naxis2+(crval2/cdelt2))+']'
#       binning or on chip summation parameters
#       Note that as far as CCDPROC is concerned, a chip using binning
#       is a different "instrument" than the same chip not using binning.
	ccdsum=str(cdelt1)+' '+str(cdelt2)
#       finally, set the OVERSCAN value to indicate who we are
	date | scan(s1,s2,s3,s4,s5,s6)
	overstr='LICKBASE '+s1+' '+s6+' '+s2+' '+s3+' '+s4+' '+s5
#       It would also be nice if we could set the ORIGSEC card in the
#       fashion of a NOAO site.  However, this keyword is not used
#       anywhere in CCDPROC and there is no way we could infer it from
#       a Lick-generated header.  (well, ok, knowing that the origin pixel
#       of a Lick detector is (0,0) we could tell if the readout window did
#       not include the origin pixel.  But we could not tell how much bigger
#       than the readout window is the actual full imaging area of the CCD.)
#       So we punt on trying to generate ORIGSEC as such ...
#
#       However, it may still be a valuable thing to have an IRAF-style
#       indication of what the original image pixel limits were.  The
#       full frame of the CCD may well be bigger, and this may be inferable
#       in the manner described above, but we shall ignore that here.
#       LICKPIX is a 1-based (not 0-based) rectangle indicating which
#       imaging pixels were actually part of the readout window.
#       LICKPIX should not be modified by ccdproc when it trims things.
#       (Note that DATASEC and TRIMSEC cards are destroyed and CCDSEC is
#       modified by ccdproc when it trims the image.)
	lickpix='['+str(crval1+1)+':'+str(crval1+((naxis1-1)*cdelt1))+','+str(crval2+1)+':'+str(crval2+(naxis2*cdelt2))+']'

#       Regenerate a bias frame consisting of a column-replicated version
#       of the bias which was applied by the Lick data acquisition system.
#       First, grab the last column (in[naxis1,*]) and stuff it into biasframe
	imstack (in // basesect, biasframe)
#       imtrans is needed since the "column" taken via imstack
#       was stored as a "row" (biasframe[nrows,1]) instead of
#       what we would have expected (biasframe[1,nrows]).
#       (Recall that images are usually denoted [naxis1(cols),naxis2(rows)])
	imtrans (biasframe, biasframe)
#       block replicate the image up to the size of the input frame
	blkrep (biasframe, biasframe, naxis1, 1)

#       use noao.imred.bias.colbias to perform any necessary refinement
#       of the bias which was originally applied when this frame was
#       taken by the Lick data aquisition system
	colbias (in, out, trim=wholsec, bias=basesect)
#       Add in all the header cards for the sake of ccdproc
	hedit(out, 'biassec' , basesect, add=yes, show-, ver-)
	hedit(out, 'trimsec' , trimsec , add=yes, show-, ver-)
	hedit(out, 'datasec' , datasec , add=yes, show-, ver-)
	hedit(out, 'ccdsec'  , ccdsec  , add=yes, show-, ver-)
	hedit(out, 'lickpix' , lickpix , add=yes, show-, ver-)
	hedit(out, 'ccdsum'  , ccdsum  , add=yes, show-, ver-)
	hedit(out, 'overscan', overstr , add=yes, show-, ver-)
	########
	if ( observat != "" )
	hedit(out, 'observat', observat, add=yes, show-, ver-)
	if ( dispaxis > 0 )
	hedit(out, 'dispaxis', dispaxis, add=yes, show-, ver-)
	########
	# determine if binned
	binned = ( (cdelt1 > 1) || (cdelt2 > 1) )
	########
	# determine if fast mode readout
	clkspeed=""
	hselect( out, 'clkspeed', yes) | scan(clkspeed)
	fastmode = (clkspeed == "Fast" )
	########
	if ( fastmode && rdnoisef >= 0. ) {
	    hedit(out, 'rdnoise' , rdnoisef , add=yes, show-, ver-)
	} else if ( rdnoise >= 0. ) {
	    hedit(out, 'rdnoise' , rdnoise , add=yes, show-, ver-)
	}
	########
	if (fastmode && binned && gainbf > 0. ) {
	    hedit(out, 'gain'    , gainbf  , add=yes, show-, ver-)
	} else if (fastmode && gainf > 0.) {
	    hedit(out, 'gain'    , gainf   , add=yes, show-, ver-)
	} else if (binned && gainb > 0.) {
	    hedit(out, 'gain'    , gainb   , add=yes, show-, ver-)
	} else if ( gain > 0. ) {
	    hedit(out, 'gain'    , gain    , add=yes, show-, ver-)
	}
	########

#       Finish the job by adding back in the bias which was originally
#       taken out by the Lick data acquisition system.
	imarith (out, "+", biasframe, out)

#       Get rid of the temporary biasframe
	imdelete (biasframe, ver-, >& "dev$null")
    }
#   clear these two lists so we don't remember their content in the par file
    ilist = ""; olist = ""

#   get rid of the two temporary files we created as workspace
    delete (ifile // "," // ofile, ver-)
end