Algorithm

NOTE!! This is now out of date - I haven't got around to writing a new one.

   subroutine lccorr

   foreach (srcTS, bkgTS) {

     # Do basic existence and sense checks on the datasets.

     RATE = COUNTS / bin duration
     ERROR = ERROR / bin duration

     # Retrieve region from DSS and make list of all pixels within it.

     # Calculate the total area of the region.

     if (use measured spectrum) {
       # Rebin spectrum if specified by user.
     } else {
       # Make a 'natural' or model spectrum.
     }

     # Obtain the effective area as a function of energy.

     # Calculate the correction for GTIs as a function of time.
     tdx = calcGtiCorrection()

     # Calculate the other time-dependent exposure factors (cosmic rays and
     # dead time).
     tdx = tdx * calcOtherTimeDepCorrns()

     # Now calculate spatial contributions.
     if (numCcdsInRegion == 1) {
       if (withspecset) {
         # Routine 00: 1 ccd/node, spectrum supplied.
         (numerator, denominator) = calcSpatial00(ccd, node, imageModel, etc)

       } else {
         # Routine 01: 1 ccd/node, model spectrum used.
         (numerator, denominator) = calcSpatial01(ccd, node, imageModel, etc)
       }
       FRACEXP = tdx * numerator / denominator

     } else {
       if (withspecset) {
         # Routine 10: many ccds/nodes, spectrum supplied.
         (numerator, denominator) = calcSpatial10(imageModel, etc)
         FRACEXP = tdx * numerator / denominator

       } else {
         # Routine 11: many ccds/nodes, model spectrum used.
         FRACEXP = calcFracexp11(tdx, imageModel, etc)
       }
     }

     # Correct for external-source OOTEs.

     # Correct for binning noise.

     if (applyCorrections) {
       foreach (timebin) {
         if (FRACEXP(timebin)) > 0) {
           RATE[timebin]  =  RATE[timebin] / FRACEXP(t)
           ERROR[timebin] = ERROR[timebin] / FRACEXP(t)
         } else {
           RATE[timebin]  = setToNull
           ERROR[timebin] = setToNull
         }
       }
     }
   } # end foreach (srcTS, bkgTS)

   # Normalise bkgTS to src collection area.

   if (subtractBkg) {
     foreach (timebin) {
       RATE[timebin] = RATE[timebin] - BACKV[timebin]

       ERROR[timebin] = sqrt(ERROR[timebin] * ERROR[timebin]
       + bkgERROR[timebin] * bkgERROR[timebin])
     }
   }

   # Delete bkg FRACEXP column from output file.

   # Release output file.

   end subroutine lccorr

The workings of the four `spatial' routines are very similar to each other. To serve as an example, the first is laid out below.

   subroutine calcSpatial00(ccd, node, imageModel, etc)

   # startEbin = first Ebin with centre energy greater than TS minimum E.
   # stopEbin = last Ebin with centre energy less than TS maximum E.

   numerator = 0.0
   denominator = 0.0
   if (imageModel eq 'point') { # eqn 33, 34 of SSC-LUX-TN-0053
     foreach Ebin (startEbin .. stopEbin) {
       call CAL_getPSFmap(PSFmap)
       expr = 0.0
       foreach pixel (listOfPixels) {
         call CAL_getFilterTransmission(filterTrans)

         totalQe = 0.0
         foreach patternID (0 .. 31) {
           qeOfPattern = CAL_getQuantumEfficiency(patternID etc)
           totalQe = totalQe + qeOfPattern
         } # next patternID

         expr = expr + filterTrans * totalQe * psf(srcX, srcY, pixelX, pixelY, PSFmap)
       } # next pixel
       call CAL_getVignettingFactor(vigFactor)
       expr = expr * pixelArea * vigFactor

       denominator = denominator + spectrum(Ebin) * dE / expr
       numerator   = numerator   + spectrum(Ebin) * dE
     } # next Ebin

   } elsif (imageModel eq 'uniform') { # eqn 35, 36 of SSC-LUX-TN-0053
     foreach Ebin (startEbin..stopEbin) {
       expr = 0.0
       foreach pixel (listOfPixels) {
         call CAL_getFilterTransmission(filterTrans)

         totalQe = 0.0
         foreach patternID (0 .. 31) {
           qeOfPattern = CAL_getQuantumEfficiency(patternID etc)
           totalQe = totalQe + qeOfPattern
         } # next patternID

         call CAL_getVignettingFactor(vigFactor)
         expr = expr + filterTrans * totalQe  * vigFactor
       } # next pixel

       denominator = denominator + spectrum(Ebin) * dE / expr
       numerator   = numerator   + spectrum(Ebin) * dE
     } # next Ebin
     denominator = denominator * regionArea
   }

   end subroutine calcSpatial00