Sq2.50 Procedural Framework: Part 3

Task parameters: Every Task has two fundamental binary switches that control its function and roles of its built-in expressions:

1. "Ratio of primary interest": The choices are 206Pb/238U or 208Pb/232Th, and exactly one must be chosen.
2. "Method for calculating secondary ratio": The choices are "Directly, via the reference zircon" or "Indirectly, via 232Th/238U", and exactly one must be chosen. SQUID 2.50 stores this as a Boolean named Switch.DirectAltPD.

The first of these is self-explanatory (and 99+% of the time, it is 206Pb/238U). Note that 'built-in' EqNum = -1 always refers to the ratio of primary interest, no matter which is chosen.

The secondary ratio is simply whichever of the two candidates was not chosen as the primary ratio, and there are two ways in which it can be calculated. The first (but more rare!) is to simply determine the secondary ratio directly, via calculation of calibration constants analogous to those used for the primary ratio. In the case of zircon, this would entail calculating a '208Pb/232Th calib. constant' alongside the 206Pb/238U calib. constant we are already familiar with. The problem with this approach is usually the reference zircon: it is unusual for any material to be perfectly fit for purpose in both the U-Pb and Th-Pb systems, as the geochemical behaviour of U and Th are quite different in most geological environments, and properly closed U-Th-Pb systems seem to be quite difficult to maintain in nature. If (and only if) the 'direct' method is employed, ' built-in' EqNum = -2 always refers to the ratio of secondary interest, no matter which is chosen. If the 'direct' method is not selected, EqNum = -2 is null, and EqNum = -3 will be non-null (see below).

The alternative (but more common!) method of calculating the secondary ratio is indirect, and relies on the appropriate rearrangement of the general equation:

(208Pb/232Th) = (206Pb/238U) * (238U/232Th) * (208Pb/206Pb)  

Note that one of the first two ratios is independently accounted for (as the primary ratio of interest, via EqNum = -1), and the other is the secondary ratio we are interested in evaluating. In general, the fourth ratio is relatively easily measured, which means we can determine our target ratio if we have a means of evaluating the third ratio. In the case of zircon, Williams et al. (1996) formulated an empirical relationship that works well enough; it takes the (familiar) form:

232Th/238U = (0.03446 * ["254/238"] + 0.868) * ["248/254"]  

where the ratios in square brackets represent isotopic "ratios of interest" measured directly via SHRIMP. If (and only if) the 'indirect' method is employed, ' built-in' EqNum = -3 always refers to the expression specifically governing 232Th/238U. If the 'indirect' method is not selected, EqNum = -3 is null, and EqNum = -2 will be non-null (see above).

So in general, one of EqNum = -2 and EqNum = -3 will always be null, and it is possible for both to be null, in the special case where only one system is being investigated (e.g. U-Pb geochronology where no proxy for Th is being measured (e.g. a 'normal' zircon list of mass-stations lacking the 248 (ThO) peak).

So, there are four permutations of the two binary switches, and they are all valid in an isotopic sense:

• Case 1: Primary ratio = 206Pb/238U, Switch.DirectAltPD = FALSE (this is the most common in SHRIMP geochronology)
• Case 2: Primary ratio = 206Pb/238U, Switch.DirectAltPD = TRUE
• Case 3: Primary ratio = 208Pb/232Th, Switch.DirectAltPD = FALSE
• Case 4: Primary ratio = 208Pb/232Th, Switch.DirectAltPD = TRUE

Iteration-specific parameters for a Task: The above relates to fundamental properties of an individual Task. There are a series of auxiliary settings, applied to a Task solely for a specific iteration of its use. In theory, these auxiliary settings should not interact with the fundamental Task properties... in practice, the way SQUID 2.50 is constructed, there is some overlap.

We have already dealt with two of these 'auxiliary' settings:

1. SBM-normalisation (on or off), and
2. Isotope-ratio calculation (weighted average or linear regression)

A third, which unfortunately does interact with the four Cases above, is the choice of 'index isotope'. The 'index isotope' refers to the identity of the isotope of Pb used to calculate the proportion of non-radiogenic Pb present in any given analysis, and perform the appropriate correction. The two most-frequently used options are 204Pb and 207Pb, as neither are used in the fundamental daughter-parent ratios that form the basis of U-Th-Pb geochronology. SQUID 2.50 does offer a third option: 208Pb. Note, however, that this option is not available to any Task involving direct calculation of a 208Pb/232Th calibration constant (irrespective of whether 208Pb/232Th is the primary or secondary ratio), and this restriction rules out Cases 2-4. As a result, when index isotope is added to the permutations, there are nine valid combinations requiring consideration:

• Case 1a: Primary ratio = 206Pb/238U, Switch.DirectAltPD = FALSE, index isotope = 204Pb
• Case 1b: Primary ratio = 206Pb/238U, Switch.DirectAltPD = FALSE, index isotope = 207Pb
• Case 1c: Primary ratio = 206Pb/238U, Switch.DirectAltPD = FALSE, index isotope = 208Pb
• Case 2a: Primary ratio = 206Pb/238U, Switch.DirectAltPD = TRUE, index isotope = 204Pb
• Case 2b: Primary ratio = 206Pb/238U, Switch.DirectAltPD = TRUE, index isotope = 207Pb
• Case 3a: Primary ratio = 208Pb/232Th, Switch.DirectAltPD = FALSE, index isotope = 204Pb
• Case 3b: Primary ratio = 208Pb/232Th, Switch.DirectAltPD = FALSE, index isotope = 207Pb
• Case 4a: Primary ratio = 208Pb/232Th, Switch.DirectAltPD = TRUE, index isotope = 204Pb
• Case 4b: Primary ratio = 208Pb/232Th, Switch.DirectAltPD = TRUE, index isotope = 207Pb

Traditionally, SQUID 2.50 requires the user to select a single index isotope, and the calculations for one of the nine permutations above are carried out. In general, I think it would be better if, once the Case number is established from the binary switch settings, the calculations were carried out for all valid index isotopes for that Case. This would save the user needing to test each index isotope via sequential iterations of data reduction (as SQUID 2.50 currently requires).

The superset of calculatable parameters, identified by column-header, is given below. In all cases, these are row-by-row calculations aimed at filling cells in an array where the rows represent individual analyses, as per the StandardData sheet produced by SQUID 2.50. Note that no single data-reduction will ever calculate all of these: every Case will exclude some of them. The parameters are grouped in the list below, to try and highlight some of the data structure.

• 4-corr206Pb/238Ucalibr.const

• 4-corr206Pb/238Ucalibr.const %err

• 4-corr206Pb/238U Age(Ma)

• 4-corr206Pb/238U Age(Ma) ±1s

• 7-corr206Pb/238Ucalibr.const

• 7-corr206Pb/238Ucalibr.const %err

• 7-corr206Pb/238U Age(Ma)

• 7-corr206Pb/238U Age(Ma) ±1s

• 8-corr206Pb/238Ucalibr.const

• 8-corr206Pb/238Ucalibr.const %err

• 8-corr206Pb/238U Age(Ma)

• 8-corr206Pb/238U Age(Ma) ±1s

• 4-corr208Pb/232Thcalibr.const

• 4-corr208Pb/232Thcalibr.const %err

• 4-corr208Pb/232Th Age(Ma)

• 4-corr208Pb/232Th Age(Ma) ±1s

• 7-corr208Pb/232Thcalibr.const

• 7-corr208Pb/232Thcalibr.const %err

• 7-corr208Pb/232Th Age(Ma)

• 7-corr208Pb/232Th Age(Ma) ±1s

• 4-corr%com206

• 7-corr%com206

• 8-corr%com206

• 4-corr%com208

• 7-corr%com208

• 4-corr208Pb*/206Pb*

• 4-corr208Pb*/206Pb* %err

• 7-corr208Pb*/206Pb*

• 7-corr208Pb*/206Pb* %err

• ppmU

• ppmTh

• 232Th/238U

• 232Th/238U %err

• 204overcts/sec(fr. 207)

• 204overcts/sec(fr. 208)

• 204/206(fr. 207)

• 204/206(fr. 207) %err

• 204/206(fr. 208)

• 7-corr206Pb/238Uconst.delta%

• 8-corr206Pb/238Uconst.delta%

• UncorrPb/Uconst

• UncorrPb/Uconst %err

• UncorrPb/Thconst

• UncorrPb/Thconst %err

• 4-corr207Pb/206Pbage

• 4-corr207Pb/206Pbage ±1s

• 4-corr207Pb/206Pb

• 4-corr207Pb/206Pb %err

• 4-corr207*/235

• 4-corr207*/235 %err

• 4-corr206*/238

• 4-corr206*/238 %err

• 4-corr errcorr

• 7-corr206*/238

• 7-corr206*/238 %err

• 7-corr errcorr

• 8-corr206*/238

• 8-corr206*/238 %err

• 8-corr errcorr

With this 'Task-anatomy' defined, we can rejoin the code. The next phase of Loop A is aimed at the evaluation of many of the columns listed above, with emphasis on the role of the Case-number in defining the exact expression used.

Loop A remains unfinished, but now that the first loop through the spot-groups has been completed (via Loop B and subsequent), it is time to reset the counters and prepare for the second loop through the spots. This second loop encompasses the calculation and placement, row-by-row, of the Daughter-Parent ratios/constants, using the built-in expressions.

Do --Start of Loop D ("second loop of spots": calculation and placement, row-by-row,  
   --of the Daughter-Parent ratios/constants)  

  Do --Start of Loop E (counting, and "basic" parsing of spot, as per Loop C)  
    piaSpotCt[-pbStd] = 1 + piaSpotCt[-pbStd]  
    piaSpotIndx[-pbStd] = 1 + piaSpotIndx[-pbStd]  
    piSpotNum = piaSpots[-pbStd, piaSpotIndx[-pbStd] ]  
  
    ParseRawData piSpotNum, FALSE, IgnoredChangedRunTable 
  Loop Until IgnoredChangedRunTable = FALSE --End of Loop E  

See separately documented subroutine ParseRawData for the argument list. Note that (1) FirstPass = FALSE here (as opposed to the prior call of ParseRawData in Sq2.50 Procedure Pt 2, and (2) the prior value of IgnoredChangedRunTable is not specified, but this is ultimately irrelevant as it is set to FALSE near the start of ParseRawData anyway. (I am not certain of the purpose of Loop E, given that every analysis that reaches this point has already been through Loop C.)

SQUID 2.50 proceeds by determining how many Daughter-Parent ratios require direct evaluation. If Switch.DirectAltPD = TRUE, then piNumDauPar = 2 (i.e. calibration constants must be calculated for both 206Pb/238U and 208Pb/232Th); otherwise, piNumDauPar = 1. SQUID 2.50 starts by eliminating the trivial case where the user has specified Switch.DirectAltPD = TRUE but provided no expression for evaluation of the calibration constant of the secondary ratio. The incomplete Loop D proceeds as follows:

  plSpotOutputRw = 1 + plSpotOutputRw  
  plOutputRw = plSpotOutputRw
  MaxDPnum = piNumDauPar  
  
  If MaxDPnum = 2 AND Switch.DirectAltPD = TRUE AND Task.Eqns[-2] = ""  
    MaxDPnum = 1  
  End If  

SQUID 2.50 then proceeds to evaluate as many calibration constants as are needed (at least one, via EqNum = -1, is mandatory), using the EqnInterp subroutine documented previously, and placing the result in the appropriate column. SQUID 2.50 does this slightly awkwardly, because it always stores the primary ratio calibration constant in the same place; it just applies different column-header names depending on whether the primary ratio is 206Pb/238U or 208Pb/232Th. I think it would probably be a better idea to simply have all the relevant columns permanently available and named, calculate whichever values are necessary/permitted, and then designated which of the calibration constants is the primary one.

  For DauParNum = 1 To MaxDPnum 
    EqnInterp Task.Eqns[-DauParNum], -DauParNum, EqnRes, EqnFerr, 1, TmpRej
  
  ValCol = IIf(pbStd, piaStdUnCorrAcol(DauParNum), IIf(pbCanDriftCorr, _
                piUnDriftCorrConstCol, piaAcol(DauParNum)))
  
  ErCol = IIf(pbStd, piaStdUnCorrAerCol(DauParNum), piaAeCol(DauParNum))

The next stage of the VBA code deals with 'looping' each spot through a set of procedures aimed at validating the analysis in itself, validating it against the governing Task, and evaluating Custom expressions and Built-in expression in a logical sequence, before proceeding to the Task-independent 'central processing of the U-Pb data.

Unfortunately, in the early sections of the code, Ludwig has made extensive use of the numeric values assigned to Booleans (in VBA, TRUE = -1 and FALSE = 0), as array-references (for a series of two-column, zero-addressed arrays). He sought to exploit the fact that, numerically, FALSE = -FALSE, and that -TRUE = 1. No doubt this is extremely 'clever', but it makes the code unnecessarily difficult to read and interpret. In this wiki, the original 'Booleans-as-indices' usage has been retained (to guard against the introduction of errors during translation), but wherever possible, plain-English comments have been added in order to clarify what is going on.

Some Definitions

Most of these appear in the VBA code, although they have been extended and reworded in places, for clarity.

• Integer piSpotNum (range from 1 to total number of spots analysed) is the index of the total vector of Spots, in time sequence.
• Integer vector piaNumSpots is a zero-addressed, two-element vector containing counts of the number of analyses by type. piaNumSpots[0] = number of Sample spots in the run (91 in demo XML); piaNumSpots[1] = number of Reference 206/238 material spots in the run (23 in demo XML).
• Integer array piaSpots is a two-row array in which the rows are zero-addressed (range 0 to 1) but the columns are not (range 1 to max(piaNumSpots) i.e. 91 in demo XML). The first row (i.e. piaSpots[0, x]) contains the index numbers (i.e piSpotNum values) for each Sample spot in the run; the second row (i.e. piaSpots[1, x]) contains the index numbers (i.e piSpotNum values) for each Reference 206/238 material spot in the run. Unused elements in the shorter of the two rows default to zero.
• Integer vector piaSpotIndx is a zero-addressed, two-element vector containing the index of piaSpots[w, x] of the current spot being processed. piaSpotIndx[0] = piaSpots[0, x] when a Sample spot is being processed, and piaSpotIndx[1] = piaSpots[1, x] when a Reference 206/238 material spot is being processed.
• Integer vector piaSpotCt is a zero-addressed, two-element vector containing counts of the number of analyses processed so far. piaSpotCt[0] = number of Sample spots processed so far; piaSpotCt[1] = number of Reference 206/238 material spots processed so far.

Paraphrased Code

pbStd = FALSE  

"Start of Reference Material-Sample Loop":

Do --Start of Loop A  
  pbStd = (pbStd = FALSE) Or (pbStdsOnly = TRUE)  
  --The first time through Loop A, this sets pbStd = TRUE,
  --in which case integer/index "[-pbStd]" = 1  
  
  pbDone = FALSE  
  
  Do --Start of Loop B ("first loop of spots")  

    Do --Start of Loop C (counting, and "basic" parsing of spot)  
      
      piaSpotCt[-pbStd] = 1 + piaSpotCt[-pbStd]  
      piaSpotIndx[-pbStd] = 1 + piaSpotIndx[-pbStd]  
      piSpotNum = piaSpots[-pbStd, piaSpotIndx[-pbStd] ]  
      
      ParseRawData piSpotNum, TRUE, IgnoredChangedRunTable, DateStr, TRUE, FALSE, TRUE  

See separately documented subroutine ParseRawData for the argument list. Note that (1) FirstPass = TRUE here (as opposed to the prior call of ParseRawData in GetConcStdData, documented in Sq2.50 Procedure Pt 1, and (2) the prior value of IgnoredChangedRunTable is not specified, but this is ultimately irrelevant as it is set to FALSE near the start of ParseRawData anyway. The incomplete Loop C proceeds as follows:

      If (IgnoredChangedRunTable = TRUE) And (SpotNscans > 1)  
        MsgBox("Run Table changes at spot number " & StR[piSpotNum] & " -- terminating.")
        CrashEnd  
      End If  
      
    Loop Until IgnoredChangedRunTable = FALSE --End of Loop C

Loop C usually only runs once because IgnoredChangedRunTable is usually FALSE throughout. When it does get set to TRUE (via failure of a ParseRawData test), it still only runs once, because CrashEnd is immediately triggered in the vast majority of cases. But note that there seems to be a gap in the VBA code logic here, in the case where an analysis with SpotNscans = 1 failed a ParseRawData test, resulting in IgnoredChangedRunTable = TRUE. It looks like such an analysis could never escape Loop C... but I don't have any 'real data' to test this.

The incomplete Loop B proceeds as follows:

    DateStr = psaSpotDateTime[piSpotNum]
    ParseTimedate DateStr, Seconds 

ParseTimeDate is a Ludwig function designed to convert a date-string (DateStr) into its representation as a number of seconds (Seconds) elapsed since the commencement of calendar 1990 (seems arbitrary; presumably chosen because SHRIMP output was not computerised prior to that). The incomplete Loop B proceeds as follows:

    If FirstSecond = 0
      FirstSecond = Seconds
    End If
    
    plSpotOutputRw = 1 + plSpotOutputRw

In the following, 'CFs' reflects a Ludwig function best considered as 'Cell Fill with String', and 'CF' reflects 'Cell Fill'. In each case, the function appears to have three arguments: the row the cursor is to be placed in, the column the cursor is to be placed in, and the value to be placed. The incomplete Loop B proceeds as follows:

    CFs plSpotOutputRw, 1, psSpotName
    CFs plSpotOutputRw, piDateTimeCol, DateStr
    CF plSpotOutputRw, piHoursCol, Excel.Fixed((Seconds - FirstSecond) / 3600, 3)

Excel.Fixed invokes the Excel spreadsheet function FIXED, which returns a text representation of a number rounded to a specified number (3 here) of decimal places (e.g. https://www.techonthenet.com/excel/formulas/fixed.php). The third of these three statements calculates the 'Hours' value, as shown in the third column of processed SQUID-books, as a double-precision number (to 3 decimal places) defining the time elapsed in the analytical session since the commencement of the first Reference 206/238 material analysis.

At the time the code was written, this was presumably supposed to rigorously reflect the first analysis collected, chronologically. However, at Geoscience Australia (and probably elsewhere), sessions rarely commence with an analysis of the Reference 206/238 material, with the result that a small handful of spots have negative Hours values. This is not an operational problem: the primary function of the Hours column is to provide easy chronological sorting of any set of analyses and to provide an easy 'axis' for time-dependent representations and calculations (e.g. X-Y plots where X is analytical session time, and calculation of any associated X-Y regressions).

The incomplete Loop B proceeds as follows:

    GetRatios Ratios(), RatioFractErrs(), False, BadSbm()  
    PlaceRawRatios plSpotOutputRw, Ratios, RatioFractErrs

This pair of statements 'gets and places' the 'isotope ratios' for this row (i.e. in the initial instance, the first (chronological) analysis of the reference 206Pb/238U material). In addition to the Spot Name, Spot Date/Time and Hours data placed previously by the CFs and CF statements above, this pair of statements results in the placement of:

• All spot-specific data harvested directly from the XML file (in SQUID 2.50, these are the QT1Y and QT1Z bit-values, the Stage X, Stage Y and Stage Z positions, and the primary beam intensity in nA).
• BackgroundCps and all TotCps values calculated at the end of SHRIMP: Step 2
• All RatioVal and RatioFractErr (the latter as percentage) values, as calculated for each 'ratio of interest' by SHRIMP: Step 4. Note that this does not include the results of any Task expressions, as these calculations have not yet been performed!

The procedure then turns firstly to the custom expressions of the governing Task. This first pass through the expressions is (I think) primarily aimed at finding which of the 'custom' expressions requires evaluation as input to any of the 'built-in' expressions. The incomplete Loop B proceeds as follows:

    For EqNum = 1 to Task.Neqns --Remember 'built-ins' have EqNum < 0

      If Switch.LA = FALSE And Switch.SC = FALSE
    
        If Switch.FO = TRUE Or Switch.AR = TRUE Or piaEqnRats[EqNum] = 0
          --Second condition is redundant, I think

          Formulae Task.EqnAsString[EqNum], EqNum, pbStd, piaEqCol[pbStd, EqNum]
          --to be documented/paraphrased separately
          
          If Switch.FO = FALSE
            EqnResu = Cells[SpotOutputRw, piaEqCol[pbStd, EqNum] ]
          End If
      
        Else --essentially equivalent to "If Switch.NU = TRUE"
      
          EqnInterp Task.EqnAsString[EqNum], EqNum, EqnRes, EqnFerr, 1, 0, TRUE  
          --See separate documentation of Sub EqnInterp for argument list 
      
          EqnResu = StR( EqnRes ) --double-to-string conversion 
          EqnFerro = StR( 100 * EqnFerr ) --formal construction of %err from Ferr
          CFs plSpotOutputRw, piaEqCol[pbStd, EqNum], EqnResu
          CFs plSpotOutputRw, piaEqECol[pbStd, EqNum], EqnFerro
          --Note "piaEqCol" vs "piaEqECol" in the last two lines
          
        End If
    
      End If
    
    Next EqNum
    
  Loop Until piaSpotIndx[-pbStd] = piaNumSpots[-pbStd] --End of Loop B

This marks the end of the first loop through a set of spots (i.e either reference 206Pb/238U materials or samples, depending on the stage of the processing. Note that Loop A is not yet closed, and the next stage of the processing revisits the Task's Custom expressions, this time looking at 'single-cell' (or 'single-array') results generated from columns of pre-existing data. The incomplete Loop A proceeds as follows:

  plHdrRw = flHeaderRow[pbStd] 
  --Fixed index-number of header-row: 6 for ref 206/238 material, 2 for samples
  Frw = FirstDatRw[-pbStd]
  --First data-row: essentially = 1 + plHdrRw 
  --(i.e. 7 for ref 206/238 material, 3 for samples)
  Lrw = LastDatRw[-pbStd]
  --Last data-row: essentially = piaNumSpots[-pbStd] + plHdrRw (i.e. for demo XML,
  --ref 206/238 material = 23 analyses, plHdrRw = 6, so Lrw = 29, and for
  --samples = 91 analyses, plHdrRw = 2, so Lrw = 93.

With Frw and Lrw defined, the next step is to place any 'single-cell' (or 'single-array') expressions for which Switch.LA (i.e. 'LAst') = FALSE. The incomplete Loop A proceeds as follows:

  For EqNum = 1 to Task.Neqns --Custom expressions only!
  
    If Switch.LA = FALSE And ( (Switch.AR = TRUE And Switch.ARrows = 1) Or 
      (Switch.SC = TRUE) ) --logic repaired from Ludwig; see below
      
      Formulae Task.EqnAsString[EqNum], EqNum, pbStd, Frw, piaEqCol[pbStd, EqNum]
      --to be documented/paraphrased separately
    End If
  
  Next EqNum

The above If statement incorporates repairs to the conditions originally imposed by Ludwig: his VBA code omitted "Switch.AR = TRUE And" under the assumption that this Boolean was implied by a non-zero value for Switch.ARrows. Unfortunately, a bug in the Task Editor permits the proliferation of non-zero values for Switch.ARrows and Switch.ARcols, even when Switch.AR = FALSE, and the real presence of this prohibited combination permitted expressions with Switch.AR = FALSE to satisfy the If condition (it is safe to say this was unintended!). Consequences included 're-evaluation' of NU-switched expressions as FO-switched (in the first data-row only), as observed by Jim Bowring in the 29 April 2017 iteration of the 100142_ShowcaseTask... output.

Loop A remains unfinished, but now that the first loop through the spot-groups has been completed (via Loop B and subsequent), it is time to reset the counters and prepare for the second loop through the spots.

plSpotOutputRw = plHdrRw
piaSpotIndx[-pbStd] = piaStartSpotIndx[-pbStd] - 1
piaSpotCt[-pbStd] = 0

The second loop through the spot-groups encompasses the calculation and placement, row-by-row, of the Daughter-Parent ratios/constants, using the built-in expressions. This will be documented in Procedural Framework: Part 3.