Skip to content

SHRIMP: Step 3

Jump to bottom Edit New page
sbodorkos edited this page Apr 18, 2017 · 7 revisions

2016.04.13

Step 3: Calculation of interpolated ratios of measured species

The next step is to calculate 'isotopic ratios', which in this context refers strictly to ratios in which both the numerator and the denominator correspond to species which have been directly measured as part of the SHRIMP data acquisition process. The following routine expects the index of the species constituting the numerator (denoted NUM) and the index of the species constituting the denominator (denoted DEN) to be defined. Indices from 1 to Nspecies are based on the order (within a scan) in which the species were acquired during analysis, which by convention (but not compulsorily) is in order of increasing mass within a scan. In the demo XML, species 196(Zr2O) has index 1 because it was acquired first, 204Pb = 2, background = 3, 206Pb = 4, 207Pb = 5, ..., 270(UO2) = 10.

The list of isotopic 'ratios of interest' is user-specified in SQUID, but can be hard-wired here for prototyping purposes. The SQUID processing of the demo XML used a 'recipe' in which 10 'ratios of interest' were defined (this has nothing to do with the fact that Nspecies = 10 in the demo XML: just a coincidence), as follows (see columns T–AM on the worksheet 'StandardData' in the SQUID-book 100142_G6147_original_frozen.xls):

  1. 204Pb/206Pb (i.e. NUM = 2, DEN = 4; an example where DEN > NUM: see below)
  2. 207Pb/206Pb (i.e. NUM = 5, DEN = 4; an example where NUM > DEN: see below)
  3. 208Pb/206Pb (i.e. NUM = 6, DEN = 4)
  4. 238U/196(Zr2O) (i.e. NUM = 7, DEN = 1)
  5. 206Pb/238U (i.e. NUM = 4, DEN = 7)
  6. 254(UO)/238U (i.e. NUM = 9, DEN = 7)
  7. 248(ThO)/254(UO) (i.e. NUM = 8, DEN = 9)
  8. 206Pb/270(UO2) (i.e. NUM = 4, DEN = 10)
  9. 270(UO2)/254(UO) (i.e. NUM = 10, DEN = 9)
  10. 206Pb/254(UO) (i.e. NUM = 4, DEN = 9)

Note that there is no significance in the sequence of the 'ratios of interest' specified above; it simply reflects a user-tendency to look at the lower-mass species at the top of the SHRIMP analytical run-table before looking at the higher-mass species at the base. Nor is it necessary to include all of the species measured in the set of 'ratios of interest': in fact, utilising background (index = 3 in the demo XML) as a numerator or denominator is not meaningful (though SQUID does not forbid it), as all the other species measured have already been background-corrected in Step 2.

For each analysis/fraction, we need three sets of data calculated in Step 2, and a further two harvested directly from the XML:

  1. NetPkCps from Step 2, as an array of size [Nscans, Nspecies] i.e. 6 x 10 in the demo XML.
  2. SBMcps from Step 2, as an array of size [Nscans, Nspecies] i.e. 6 x 10 in the demo XML.
  3. PkFerr from Step 2, as an array of size [Nscans, Nspecies] i.e. 6 x 10 in the demo XML.
  4. time_stamp_sec from XML, as an array of size [Nscans, Nspecies] i.e. 6 x 10 in the demo XML.
  5. count_time_sec from XML, as a vector of length [Nspecies] i.e. 10 in the demo XML.

This routine is most easily implemented by establishing at the outset which of NUM and DEN is measured first in the SHRIMP acquisition run-table, because this dictates the form of the arithmetic:

If DEN > NUM  --re-express indices so that (within a scan):   
  aOrd = NUM  --aOrd always represents whichever of NUM/DEN was measured first  
  bOrd = DEN  --bOrd always represents whichever of NUM/DEN was measured last    
Else  
  aOrd = DEN  
  bOrd = NUM  
End If

At the start, background-corrected total counts are computed for each (using NetPkCps), and the number of scans counted, to determine whether the default calculation method (double-interpolation, following Dodson (1978): http://dx.doi.org/10.1088/0022-3735/11/4/004) is applicable to each specific combination of analysis/fraction NetPkCps data and desired 'ratio of interest':

TotCtsNUM = 0  
TotCtsDEN = 0  
  
For j = 1 to Nscans  
  TotCtsNUM = TotCtsNUM + ( NetPkCps[j, NUM] * count_time_sec[NUM] )  
  TotCtsDEN = TotCtsDEN + ( NetPkCps[j, DEN] * count_time_sec[DEN] )  
Next j  

Ndod = Nscans - 1

There are two situations in which double-interpolation will not be attempted:

  1. Single-scan data (i.e. Ndod = 0)
  2. At least one of the numerator or the denominator is characterised by very low counts (an arbitrary cut-off of 32 is used for total counts).

In these situations, interpolated times, values and fractional errors of ratios are either assumed or calculated arithmetically as follows:

If (TotCtsNUM < 32 Or TotCtsDEN < 32 Or Ndod = 0)  
  
  If TotCtsNUM = 0  
    RatioVal = 1E-32 --declared as SQUID's constant "tiny" value  
    RatioFractErr = 1  
  Elseif TotCtsDEN = 0  
    RatioVal = 1E+16  
    RatioFractErr = 1  
  Else  
    RatioVal = ( TotCtsNUM / count_time_sec[NUM] ) /  
               ( TotCtsDEN / count_time_sec[DEN] )  
    RatioFractErr = sqrt( 1 / Abs(TotCtsNUM) + 1 / Abs(TotCtsDEN) )
  End If
   
  --Define RatioInterpTime, InterpRatVal, RatValFerr as vectors of length 1  
   
  RatioInterpTime[1] = 0.5 * ( min( time_stamp_sec[1, NUM], time_stamp_sec[1, DEN] )  
                 + max( time_stamp_sec[Nscans, NUM], time_stamp_sec[Nscans, DEN] ) )  
  InterpRatVal[1] = RatioVal  
  RatValFerr[1] = RatioFractErr  

  --Define RatEqTime, RatEqVal, RatEqErr as vectors of length 1  
  --and write them to the Within-Spot Ratios sheet:  
      
      RatEqTime[1] = RatioInterpTime[1]  
      RatEqVal[1] = InterpRatVal[1]  
      RatEqErr[1] = Abs( RatValFerr[1] * InterpRatVal[1] )  
      
      --Bodorkos 2017-04-04: next 2 lines inserted to force VBA-Java match  
      RatioVal[1] = SigDigStrDbl( RatioVal[1], 12 ) --12 sig. digits  
      RatioFractErr[1] = SigDigStrDbl( RatioFractErr[1], 12 ) --12 sig. digits  

      Exit Sub --[i.e. go to 'sanity check' at base of wiki page,  
      --and exit this routine: bypass the rest of Step 3]  
        
End If

In most 'normal' geochronology datasets, the only species with a count-rate routinely low enough to meet the above criterion is 204Pb, so the above code is usually only invoked to calculate ratios involving 204Pb. In the demo XML, the only 'ratio of interest' involving 204Pb is 204/206... and even then, it is possible that some analyses/fractions will contain sufficient 204Pb (TotCtsNUM > 32) to skip this section of the code, and proceed to the main line of code implementing double-interpolation (below).

The rest of Step 3 pertains only to data satisfying (TotCtsNUM >= 32 And TotCtsDEN >= 32 And Ndod > 0).

PkF is defined as a vector of length Ndod. It is populated with the time elapsed (within scan j) between the first- and last-measured peaks, as a proportion of the inter-scan time of the first measured peak (i.e. from scan j to scan j+1), as follows:

For j = 1 to Ndod    
    PkF[j] = ( time_stamp_sec[j, bOrd] - time_stamp_sec[j, aOrd] ) /  
             ( time_stamp_sec[j + 1, aOrd] - time_stamp_sec[j, aOrd] )
Next j

Some mean values are calculated for later use:

AvPkF = average(PkF)  
f1 = ( 1 - AvPkF ) / 2  
f2 = ( 1 + AvPkF ) / 2  
RhoIJ = ( 1 - AvPkF ^ 2 ) / ( 1 + AvPkF ^ 2 ) / 2

Redefine each of RatioInterpTime, InterpRatVal and RatValFerr as vectors of length Ndod. Then:

Rct = 0  

For Snum = 1 to Ndod  
  Sn1 = Snum + 1  
  TotT = time_stamp_sec[Snum, aOrd] + time_stamp_sec[Snum, bOrd] +
         time_stamp_sec[Sn1, aOrd] + time_stamp_sec[Sn1, bOrd]
  MeanT = TotT / 4
  RatioInterpTime[Snum] = MeanT --'time_stamps' for interpolations
  
  ZerPkCt(Snum) = False --scan-by-scan error-trapping  
  ZerPkCt(sn1) = False  --mostly aimed at picking up  
  HasZerPk = False      --internal zeroes and error-values  

  For NumDenom = 1 To 2 --total counts for peaks
    k = Snum + NumDenom - 1
    aNetCPS = NetPkCps[k, aOrd]  
    bNetCPS = NetPkCps[k, bOrd]  

    If (aNetCPS = [error value] Or bNetCPS = [error value])  
      HasZerPk = True  
      ZerPkCt[k] = True  
      GoTo NextScanNum --trivial error-trapping  
    End If  

    aPkCts[NumDenom] = aNetCPS * count_time_sec[aOrd]
    bPkCts[NumDenom] = bNetCPS * count_time_sec[bOrd]

    If UseSBM --user-defined Boolean indicating whether SBM-normalisation  
              --is employed or not.  
      If ( SBMcps[k, aOrd] <= 0 Or SBMcps[k, aOrd] = [error value] Or  
           SBMcps[k, bOrd] <= 0 Or SBMcps[k, aOrd] = [error value] )  
        ZerPkCt(k) = True  
        GoTo NextScanNum --trivial error-trapping  
      End If
    
    End If

  Next NumDenom

  For k = 1 to 2  
    If k = 1  
      NumDenom = 2  
    Else  
      NumDenom = 1  
    End If  
  
    a = aPkCts[k]  
    b = aPkCts[NumDenom]  
    If ( a <= 0 And b > 16 )  
      ZerPkCt[Snum + k - 1] = True  
    End If  
  
    a = bPkCts[k]  
    b = bPkCts[NumDenom]  
    If ( a <= 0 And b > 16 )  
      ZerPkCt[Snum + k - 1] = True  
    End If  
  
  Next k  
  
  If ( ZerPkCt[Snum] = True Or ZerPkCt[Sn1] = True )  
    GoTo NextScanNum --not sure if this error-trapping is trivial or not  
  End If  
  
  aPk1 = NetPkCps[Snum, aOrd] --earlier scan, first-measured peak  
  bPk1 = NetPkCps[Snum, bOrd] --earlier scan, last-measured peak  
  aPk2 = NetPkCps[Sn1, aOrd]  --later scan, first-measured peak  
  bPk2 = NetPkCps[Sn1, bOrd]  --later scan, last-measured peak  
  
  If UseSBM --implementation of SBM-normalisation after error-trapping  
    aPk1 = aPk1 / SBMcps[Snum, aOrd]  
    bPk1 = bPk1 / SBMcps[Snum, bOrd]  
    aPk2 = aPk2 / SBMcps[Sn1, aOrd]  
    bPk2 = bPk2 / SBMcps[Sn1, bOrd]  
  End If  
  
  --Now revisit PkF  
  ScanDeltaT = time_stamp_sec[Sn1, aOrd] - time_stamp_sec[Snum, aOrd]  
  bTfract = time_stamp_sec[Snum, bOrd] - time_stamp_sec[Snum, aOrd]  
  PkF[Snum] = bTfract / ScanDeltaT  
  ff1 = ( 1 - PkF[Snum] ) / 2  
  ff2 = ( 1 + PkF[Snum] ) / 2  
  aInterp = ( ff1 * aPk1 ) + ( ff2 * aPk2 )  
  bInterp = ( ff2 * bPk1 ) + ( ff1 * bPk2 )  
  
  --Redefine numerator and denominator for ratio calculation:  
  If NUM < DEN  
    Rnum = aInterp  
    Rden = bInterp  
  Else  
    Rnum = bInterp  
    Rden = aInterp  
  End If  
  
  If Rden <> 0
    Rct = Rct + 1
    InterpRatVal[Rct] = Rnum / Rden
    --Approximate internal ratio variance as sum of variances from first peaks
    a1PkSig = PkFerr[Snum, aOrd] * aPk1  
    a2PkSig = PkFerr[Sn1, aOrd] * aPk2  
    b1PkSig = PkFerr[Snum, bOrd] * bPk1  
    b2PkSig = PkFerr[Sn1, bOrd] * bPk2  
    
    If UseSBM  
      a1PkSig = sqrt( a1PkSig ^ 2 +  
                ( aPk1 ^ 2 / SBMcps[Snum, aOrd] / count_time_sec[aOrd] ) )  
      a2PkSig = sqrt( a2PkSig ^ 2 +  
                ( aPk2 ^ 2 / SBMcps[Sn1, aOrd] / count_time_sec[aOrd] ) )  
      b1PkSig = sqrt( b1PkSig ^ 2 +  
                ( bPk1 ^ 2 / SBMcps[Snum, bOrd] / count_time_sec[bOrd] ) )  
      b2PkSig = sqrt( b2PkSig ^ 2 +  
                ( bPk2 ^ 2 / SBMcps[Sn1, bOrd] / count_time_sec[bOrd] ) )  
    End If  
    
    If aInterp = 0 Or bInterp = 0  
      RatValFerr[Rct] = 1  
      RatValSig[Rct] = 1E-32 --"tiny"  
      SigRho[Rct, Rct] = 1E-32 --"tiny"  
    Else  
      Term1 = ( ( f1 * a1PkSig ) ^ 2 + ( f2 * a2PkSig ) ^ 2 )  
      Term2 = ( ( f2 * b1PkSig ) ^ 2 + ( f1 * b2PkSig ) ^ 2 )  
      RatValFvar = ( Term1 / aInterp ^ 2 ) + ( Term2 / bInterp ^ 2 )  
      RatValVar = RatValFvar * ( InterpRatVal[Rct] ^ 2 )  
      RatValFerr[Rct] = sqrt( RatValFvar )  
      RatValSig[Rct] =  max( 1E-10, sqrt( RatValVar ) )  
      SigRho[Rct, Rct] = RatValSig[Rct]  
  
      If Rct > 1        
        If ZerPkCt[Snum - 1] = True  
          RhoIJ = 0  
        Else  
          RhoIJ = ( 1 - PkF[Snum] ^ 2 ) / ( 1 + PkF[Snum] ^ 2 ) / 2  
        End If  
    
        SigRho[Rct, Rct - 1] = RhoIJ  
        SigRho[Rct - 1, Rct] = RhoIJ  
      End If  
  
    End If  
    
  End If  

NextScanNum  
Next Snum

If Rct = 0  
  RatioVal = [error value]  
  RatioFractErr = [error value]  
  --[exit this routine: bypass the rest of Step 3]  
 
ElseIf Rct = 1  
  RatioVal = InterpRatVal[1]

  If RatioVal = 0  
    RatioVal = 1E-32  
    RatioFractErr = 1  
  Else  
    RatioFractErr = RatValFerr[1]  
  End If  
            
  --Bodorkos 2017-04-04: next 2 lines inserted to force VBA-Java match  
  RatioVal[1] = SigDigStrDbl( RatioVal[1], 12 ) --12 sig. digits  
  RatioFractErr[1] = SigDigStrDbl( RatioFractErr[1], 12 ) --12 sig. digits  

Else
  --Redefine each of RatioInterpTime, InterpRatVal and RatValFerr  
  --as vectors of length Rct, preserving their contents
  
  --Define RatEqTime, RatEqVal, RatEqErr as vectors of length Rct  
  --and write them to the Within-Spot Ratios sheet:  
      
  For j = 1 to Rct
    RatEqTime[j] = RatioInterpTime[j]  
    RatEqVal[j] = InterpRatVal[j]  
    RatEqErr[j] = Abs( RatValFerr[j] * InterpRatVal[j] )  
  Next j  
    
  --**INSERT "STEP 4" HERE, UP TO THE END OF THE FINAL CODE-BLOCK**
End If

Note that when the double-interpolation routine is invoked, the values for RatEqTime, RatEqVal and RatEqErr are only written out when Rct >= 2. I've triple-checked the code and this is definitely the SQUID implementation. I don't think it's right — writing nothing at all to the Within-Spot Ratios sheet when Rct = 0 or 1 doesn't seem correct — but I can only assume that this situation arises so rarely that 'normal' SQUID-users are unaffected.

Sanity check: For each analysis/fraction, there will be Ndod (i.e. Nscans - 1) rows, so for the demo XML, there will be 114 analyses x (Ndod = 5) rows per analysis = 570 rows total.

I think we will need four columns at left, with the first three being spot name, date/time, and 'scan number' (i.e. Ndod) as per the Step 1 sanity-check. The fourth column ('analysis type') should identify whether the analysis has been identified as being of the 206Pb/238U reference zircon ('standard'; corresponding to prefix 'T' in the demo XML), or not ('unknown', all other prefixes in the XML).

Columns 5–34 inclusive will constitute three sets of values for each of the 10 'ratios of interest'. Each set of values is a vector of length <= Ndod, as follows:

  1. RatEqTime (label "Time")
  2. RatEqVal (label "[name of NUM]/[name of DEN]")
  3. RatEqErr (label "+/- 1sig absolute")

Finally, we'll need some exotic sorting, to facilitate comparison with the SQUID-book:

  1. 'analysis type' alphabetical ascending (or whatever is needed to ensure that all 'standards' precede all 'unknowns')
  2. 'date/time' ascending (so within the 'analysis type' subdivision, analyses appear in chronological order, as per Step 1 sanity-check)
  3. 'scan (Ndod) number' ascending (so within each analysis, interpolation-rows appear in chronological order)

When the rows are sorted, the contents of the 30 columns (5–34) will be closely comparable to the 30 columns (B–AE) in worksheet 'Within-Spot Ratios' in SQUID-book 100142_G6147_original_frozen.xls.

Note on 'Mean times' and assumed perfect periodicity of the Dodson double-interpolation as implemented by SQUID
As documented above, the mean time of each interpolated is defined by the arithmetic mean of the four measurement times contributing to the ratio:

For Snum = 1 to Ndod  
  Sn1 = Snum + 1  
  TotT = time_stamp_sec[Snum, aOrd] + time_stamp_sec[Snum, bOrd] +
         time_stamp_sec[Sn1, aOrd] + time_stamp_sec[Sn1, bOrd]
  MeanT = TotT / 4
  RatioInterpTime[Snum] = MeanT  
Next Snum

where Snum is the scan number, aOrd is the index of whichever of the numerator and denominator is measured first within a scan, and bOrd is the index of whichever of the numerator and denominator is measured last within a scan.

This calculation creates the impression that the target interpolated values of aOrd and bOrd can be obtained in a simple geometric way, by linear interpolation of the scan-by-scan values to the target 'mean time'. But this is not the case, because the double-interpolation procedure illustrated in Ludwig (2009: http://dx.doi.org/10.1016/j.chemgeo.2009.07.004) and implemented in SQUID 2.50 uses ScanDeltaT and bTfract, which are defined as:

  ScanDeltaT = time_stamp_sec[Sn1, aOrd] - time_stamp_sec[Snum, aOrd]  
  bTfract = time_stamp_sec[Snum, bOrd] - time_stamp_sec[Snum, aOrd]

An important consequence of this is that the true time value of the peak not included in these two parameters (i.e time_stamp_sec[Sn1, bOrd]) is not used in the interpolation procedure. Instead, it is directly replaced with time_stamp_sec[Snum, bOrd] + ScanDeltaT (i.e. time_stamp_sec[Snum, bOrd] + time_stamp_sec[Sn1, aOrd] – time_stamp_sec[Snum, aOrd]) using the explicit assumption that the ScanDeltaT values for the first-measured and last-measured peaks are identical (i.e. perfect periodicity).

This approach is probably convenient in terms of arithmetic implementation, but it is not a particularly accurate assumption for real SHRIMP data. The experience at Geoscience Australia is that 1-second offsets in ScanDeltaT values are common, presumably arising from delays in processing instructions between the SHRIMP and the control software. But in addition, the double-interpolation procedure implemented in SQUID 2.50 takes no account of non-trivial peak-time offsets, such as that potentially arising from re-measurement of a peak where the original measurement was a 'stuff-up'. Such peak measurements are labelled 'S' by the SHRIMP control software in the acquired total-count data (independent of SQUID 2.50), indicating the failure of an F-test and/or linear interpolation test applied to the 10 constituent integrations at the peak. In this situation, the SHRIMP control software will try at least once to re-centre and re-measure the peak, thereby extending the inter-scan times and overall duration of the analysis.

The impact of these inconsistencies is likely to be minor for isotopic ratios calculated (in Step 4) as 'spot averages', because the interpolated values will only be affected by the distance between the assumed and actual mean times for each interpolation. However, it will be greater for isotopic ratios calculated (in Step 4) via linear regression to analysis mid-time, as the assigned 'mean times' for the interpolated ratios (and their differences from the 'true times' to which the ratio values and their uncertainties are actually interpolated) will affect the fit of the regression.

Next: Step 4: Calculation of 'spot-mean' for each 'ratio of interest' for each analysis

Home - SHRIMP Development - U Series Development

Clone this wiki locally