SHRIMP: Step 4

2016.04.21

Step 4: Calculation of 'mean' ratio for each measured species for each analysis

Having derived the Ndod (= Nscans - 1) interpolated ratio values and errors in Step 3, the next step is to combine these values into a single isotopic ratio (per 'ratio of interest', per analysis), to match the values tabulated for the measured nuclide-ratios on the StandardData and SampleData sheets of the SQUID-workbook.

SQUID offers users a choice of calculation method: EITHER the error-weighted linear regression of interpolated ratios vs time, with the regression-error calculated at analysis midtime OR the error-weighted mean of interpolated ratios and its error, independent of time. (In fact, SQUID allows sub-choices, depending on whether each 'ratio of interest' involves nuclides of a single element or not, but we will ignore the finer distinctions for the moment, partly because Geoscience Australia doesn't use them anyway.)

This implies the existence of a Boolean 'UserLinFits' (= TRUE when the user preference is for mean ratios to be calculated via error-weighted linear regression vs time, FALSE when the user preference is for mean ratios to be calculated via error-weighted mean independent of time).

As stipulated in the early part of Step 3, single-scan data (i.e. Ndod = Rct = N = 0) will never be subjected to Dodson interpolation, and thus will never occur as input to WtdLinCorr. Similarly, at the end of Step 3, specific provision is made for two-scan data (i.e. Ndod = Rct = N = 1), because although Dodson interpolation has taken place, it has generated only a single value and error, so there is nothing to 'average'. Therefore, Step 4 is applicable only to nuclide-ratios collected over three or more scans (i.e. Ndod = Rct = N > 1; admittedly, this encompasses the vast majority of 'proper' SHRIMP data).

But in addition to this, Ludwig specified that linear regression to analysis midtime calculation is not permitted unless Ndod = Rct = N > 3, regardless of the specified UserLinFits (i.e. if UserLinFits = TRUE, but N = 2 or 3, then the calculation is performed as error-weighted average instead, without user consultation). These scenarios are best handled by logic preceding the call of the calculation routines (as per the code below), but it should be noted that the WtdLinCorr subroutine (invoked below) contains its own check, as a backup that is redundant in this particular (Step 4) application of WtdLinCorr subroutine.

The 'Step 4' code proceeds (noting that the case of interest for first-pass implementation and testing is UserLinFits = FALSE):

  If UserLinFits = True And Rct > 3
    --Calculate error-wtd linear regression of ratios vs time, with error calculated  
    --at mid burn-time, using index-adjacent ratio error-correlations
    WtdLinCorr 2, Rct, InterpRatVal, SigRho, MSWD, Probfit, 0, RatioInter, ...
      SigmaRatioInter, Bad, RatioSlope, SigmaRatioSlope, CovSlopeInter, RatioInterpTime

    MidTime = ( time_stamp_sec[Nscans, Nspecies] + time_stamp_sec[1, 1] ) / 2  
    RatioMean = (RatioSlope * MidTime) + RatioInter
    RatioMeanSig = sqrt( (MidTime * SigmaRatioSlope) ^ 2 + ...
      (SigmaRatioInter ^ 2) + (2 * MidTime * CovSlopeInter) )

  Else  
    --Calculate error-wtd mean, using index-adjacent ratio error-correlations
    WtdLinCorr 1, Rct, InterpRatVal, SigRho, MSWD, Probfit, 0, RatioMean, ...
      RatioMeanSig, Bad

  End If

The WtdLinCorr subroutine is documented separately () as it is a generalised function used by SQUID 2.50 in a variety of contexts. For the purpose of Step 4-specific orientation of the above code relative to WtdLinCorr:

• the first WtdLinCorr argument (i.e. 1 or 2) is the input integer value of Avg1LinRegr2.
• the second WtdLinCorr argument (i.e. Rct from Step 3) is the input integer value of N.
• the third WtdLinCorr argument (i.e. InterpRatVal from Step 3) is a vector of length N, as the input Y.
• the fourth WtdLinCorr argument (i.e. SigRho from Step 3) is an array of size N x N, as the input SigRho.
• the remaining arguments of WtdLinCorr are outputs of the subroutine.

Initial testing of Step 4 (i.e. UserLinFits = FALSE) requires only the WtdLinCorr, DeletePoint, and WtdAvCorr subroutines, all of which are now documented. Note that the subroutine WeightedLinearCorr (called by WtdLinCorr if UserLinFits = TRUE and N > 3) remains undocumented as of 21 June 2016: that case will be considered once the FALSE case is implemented.

Once the WtdLinCorr subroutine has returned values of RatioMean and RatioMeanSig:

If Bad = TRUE --error-trapping  
  RatioVal = [error value]
  RatioFractErr = [error value]
ElseIf RatioMean = 0  
  RatioVal = 1E-32
  RatioFractErr = 1
Else  
  RatioVal = RatioMean
  RatioFractErr = Max(1E-32, RatioMeanSig) / Abs(RatioVal)
End If  

The final step is to multiply each RatioFractErr by 100, to make it a percentage. This operation is part of the cell-filling routine in SQUID, which I don't think is particularly helpful.