-
Notifications
You must be signed in to change notification settings - Fork 16
SHRIMP: Sub SamRadiogenicCols
This subroutine (which is solely for SampleData) places, row-by-row, formulae to calculate radiogenic (i.e. corrected for common Pb) Pb/Pb, Pb/U, and Pb/Th ratios and ages in the SampleData sheet.
SamRadiogenicCols plaFirstDatRw, plaLastDatRw
plaFirstDatRw: Integer index number of the first row containing spot-by-spot data (for the Standard).
plaLastDatRw: Integer index number of the last row containing spot-by-spot data (for the Standard).
Values of type Boolean
pbTh, pbU
Values of type Integer
f, L, piNumDauPar
Values of type String
q, s, SA, sae, t, WtdMnA
Values of type Range
rw1, rw2, r4
The subroutine starts by calculating the Pb-isotope ratios used for the common-Pb correction. Note that these are not predicated on the user having specified "the right" numerator-denominator combinations as part of their Task definition: rather, SQUID 2.50 mandates that ratios ["204/206"], ["207/206"] and ["208/206"] are always calculated if the relevant mass-stations are present in the run-table, even if the user did not specifically request those ratios (i.e. they are evaluated 'in the background'). Remembering that SQUID 2.50 constructs its formulae as strings:
Alpha = " 1 / ["204/206"] --measured ["206/204"], i.e. total 206/204
Beta = " ["207/206"] / ["204/206"] " --measured ["207/204"], i.e. total 207/204
Gamma = " ["208/206"] / ["204/206"] " --measured ["208/204"], i.e. total 208/204
In each case, the measured (or total) Pb-isotope ratio can be considered as the sum of the radiogenic component and the common component. Using ["206/204"] as an example, R64 as shorthand for the ratio, and suffixes m, r, and c to denoted measured, common, and radiogenic components respectively: R64m = R64r + R64c. It is then possible to isolate the radiogenic component, by subtracting the common component from the measured value:
NetAlpha = " ( ["Alpha"] - sComm_64 ) " --i.e. R64r = R64m - R64c
NetBeta = " ( ["Beta"] - sComm_74 ) " --i.e. R74r = R74m - R74c
NetGamma = " ( ["Gamma"] - sComm_68 ) " --i.e. R84r = R84m - R84c
For the daughter isotopes of calibration-related interest (i.e 206Pb and 208Pb), it is possible to define proportionality factors describing the radiogenic component as a fraction of the total measured value:
radd6 = " ( ["NetAlpha"] / ["Alpha"] ) " --i.e. R64r / R64m
radd8 = " ( ["NetGamma"] / ["Gamma"] ) " --i.e. R84r / R84m
Start by calculating one or both of ["Total 206Pb/238U"] and/or ["Total 208Pb/232Th"], as required by Perms1-4. Note that as per previously-documented subroutine StdRadiogenicCols, this arithmetic is predicated on having only one calibration constant per daughter-parent pair (following the user-defined index isotope for the common Pb correction), whereas SQUID 3.0 will be calculating multiple calibration constants corresponding to all the candidate index-isotope values. The following expressions will therefore require generalisation:
For i = 1 To piNumDauPar -- recall piNumDauPar = 2 for Perm2 and Perm 4, else 1
If ( (i = 1) AND (pbU = TRUE) ) OR ( (i = 2) AND (pbTh = TRUE) )
-- first part of the If covers Perm1 and Perm2; the second part covers Perm4
m = piPb6U8_totCol --index number for ["Total 206Pb/238U"] column
--NOTE that SQUID 3.0 will need to key this to index-isotope, no matter how
--counter-intuitive that seems!
ElseIf ( (i = 1) AND (pbTh = TRUE) ) OR ( (i = 2) AND (pbU = TRUE) )
-- first part of the If covers Perm3 and Perm4; the second part covers Perm2
m = piPb8Th2_totCol --index number for ["Total 208Pb/232Th"] column
--NOTE that SQUID 3.0 will need to key this to index-isotope, no matter how
--counter-intuitive that seems!
Else
m = 0 --should never happen
End If
--Now place formulae for so-called ["Total 206Pb/238U"] and/or so-called
--["Total 208Pb/232Th"]. Note that in SQUID 3.0, because we are performing
--calculations for ALL candidate index-isotopes, these calculations need to
--be generalised to give results like ["Xcorr Tot 206Pb/238U"] (where X can
--be 4 or 7 (Perm1,2,4) or 8 (Perm1 only)), or ["Ycorr Total 208Pb/232Th"]
--(where Y can be 4 or 7 (Perm2,3,4)).
--We will need to evaluate the Xcorr... results carefully, because it seems
--to me that there should not be much variation, and it occurs to me that
--the calculation would be more efficently AND more correctly done by direct
--reference to our ["Uncorr Pb/U const."] and/or [Uncorr Pb/Th const."]
--columns that we calculated long ago (note that this approach would also
--yield the intuitively expected "single" ["Total 206Pb/238U"] and/or
--["Total 208Pb/232Th"]). But for the present, I will document Ludwig's
--calculations in SQUID 2.50.
--Calculate ["Total 206Pb/238U"] or ["Total 208Pb/232Th"]
If m = piPb6U8_totCol --Place row-by-row expression for ["Total 206Pb/238U"]
If i = 1 --Perm1 and Perm2
PlaceFormulae " = ["Xcorr 206Pb/238U calibr. const."] / WtdMeanA1 *
StdUPbRatio ", f, m, L --i.e. place expression in col m, rows f to L
PlaceFormulae " = SQRT( ["Xcorr 206Pb/238U calibr. const. %err"]^2
+ ExtPerrA1 ^ 2 ) ", f, m + 1, L
--i.e. place expression in col m + 1, rows f to L
--Recall that ExtPerrA is defined as max(R, S) where S is the ExtPerr
--value (at 1-sigma, as a percentage) returned by the relevant
--WtdMeanAcalc, and R is an arbitrary "minimum" value for S (R = 0.75%
--at Geoscience Australia, by convention). There is more detail on this
--at the very end of the documentation of subroutine WtdMeanAcalc: it
--is described under extBox.
ElseIf i = 2 --Perm4
PlaceFormulae " = ["Xcorr 206Pb/238U calibr. const."] / WtdMeanA2 *
StdUPbRatio ", f, m, L --i.e. place expression in col m, rows f to L
PlaceFormulae " = SQRT( ["Xcorr 206Pb/238U calibr. const. %err"]^2
+ ExtPerrA2 ^ 2 ) ", f, m + 1, L
--i.e. place expression in col m + 1, rows f to L
End If --i = 1 or 2
ElseIf m = piPb8Th2_totCol --Place expression for ["Total 208Pb/232Th"]
If i = 1 --Perm3 and Perm4
PlaceFormulae " = ["Xcorr 208Pb/232Th calibr. const."] / WtdMeanA1 *
StdThPbRatio ", f, m, L --i.e. place expression in col m, rows f to L
PlaceFormulae " = SQRT( ["Xcorr 208Pb/232Th calibr. const. %err"]^2
+ ExtPerrA1 ^ 2 ) ", f, m + 1, L
--i.e. place expression in col m + 1, rows f to L
ElseIf i = 2 --Perm2
PlaceFormulae " = ["Xcorr 208Pb/232Th calibr. const."] / WtdMeanA2 *
StdThPbRatio ", f, m, L --i.e. place expression in col m, rows f to L
PlaceFormulae " = SQRT( ["Xcorr 208Pb/232Th calibr. const. %err"]^2
+ ExtPerrA2 ^ 2 ) ", f, m + 1, L
--i.e. place expression in col m + 1, rows f to L
End If --i = 1 or 2
End If --m = pi...col
Next i