-
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"] where calibration constants are available, 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
Now calculate the corresponding ["Total 238U/206Pb"], via simple inversion:
If piPb6U8_totCol > 0 --i.e. if ["Total 206Pb/238U"] exists
--Fill column ["Total 238U/206Pb"]:
PlaceFormulae " = 1 / ["Total 206Pb/238U"] ", f, piU8Pb6_totCol, L
--Fill column ["Total 238U/206Pb %err"]:
PlaceFormulae " = ["Total 238U/206Pb"] ", f, piU8Pb6_TotEcol, L
End If
Now calculate the "secondary" daughter-parent isotopic ratio, for the permutations where no calibration constant is available. Specifically, this means calculating ["Total 208Pb/232Th"] for Perm1, and ["Total 206Pb/238U"] for Perm3. Note that the previously documented analogue of this subroutine for Standard data (Sub StdRadiogenicCols in SQUID 2.50) does not contain this step, but Squid3 should probably be revised to include it.
If (pbU = TRUE) AND (Switch.DirectAltPD = FALSE) --i.e. Perm1 only
If piPb8Th2_totCol > 0 --i.e. if column ["Total 208Pb/232Th"] exists
--Fill column ["Total 208Pb/232Th"] on SampleData:
PlaceFormulae " = ["Total 206Pb/238U"] * ["208/206"] / ["232Th/238U"] ",
f, piPb8Th2_totCol, L
--Fill column ["Total 208Pb/232Th %err"] on SampleData:
PlaceFormulae " = SQRT( ["208/206 %err"]^2 + ["Total 206Pb/238U %err"]^2 +
["232Th/238U %err"]^2 ) ", f, piPb8Th2_totEcol, L
End If
ElseIf (pbTh = TRUE) AND (Switch.DirectAltPD = FALSE) AND (piPb6U8_totCol > 0)
--i.e. Perm3 only, and even then, only if column ["Total 206Pb/238U"] exists
--Fill column ["Total 238U/206Pb"] on SampleData:
PlaceFormulae " = ["Total 208Pb/232Th"] / ["208/206"] * ["232Th/238U"] ",
f, piPb6U8_totCol, L
--Fill column ["Total 238U/206Pb %err"] on SampleData:
PlaceFormulae " = SQRT( ["208/206 %err"]^2 + ["Total 208Pb/232Th %err"]^2 +
["232Th/238U %err"]^2 ) ", f, piPb6U8_totEcol, L
End If
Now, if a uranium concentration reference material has been defined (i.e. column ["ppmU"] contains data), it is time to calculate and populate the columns containing the ppm values of radiogenic Pb, according to the various index isotope. There are a total of five permutations, comprising ["X-corr ppm 206*"], where X = 4, 7 or 8, and [Y-corr ppm 208*"], where Y = 4 or 7.
Note that SQUID 2.50 uses the 'magic number' 0.859 hard-coded, and for initial testing, this should be replicated. In terms of models of physical constants, however, the value 0.859 is derived via:
0.859 = ( 1 - ( 1 / Present238U235U ) * 206 / 238
where 206 is the atomic mass of 206Pb, 238 is the atomic mass of 238U, and Ludwig assumes a value of Present238U235U of 137.88. More recently, other values of Present238U235U have been measured via EARTHTIME (e.g. 137.818 by Hiess et al., 2012), so the hard-coded value should ultimately be replaced by the expression containing Present238U235U above.
With this in mind, the SQUID 2.50 code populates each of these 5 column in turn:
If piaPpmUcol > 0 --i.e. ["ppmU"] exists on SampleData
For p = 1 To 5
Select Case p
Case 1 --populate column ["4-corr ppm 206*"]
PlaceFormulae " = ["Total 206Pb/238U"] * ["ppmU"] * 0.859 *
( 1 - ["204/206"] * sComm_64 ) ", f, piaRadDauCol_46, L
Case 2 --populate column ["7-corr ppm 206*"]
PlaceFormulae " = ["Total 206Pb/238U"] * ["ppmU"] * 0.859 *
( 1 - ["204/206 fr. 207"] * sComm_64 ) ", f, piaRadDauCol_76, L
Case 3 --populate column ["8-corr ppm 206*"]
PlaceFormulae " = ["Total 206Pb/238U"] * ["ppmU"] * 0.859 *
( 1 - ["204/206 fr. 208"] * sComm_64 ) ", f, piaRadDauCol_86, L
Case 4 --populate column [4-corr ppm 208*"], noting that this calculation draws
--on the row-by-row results determined in Case 1
PlaceFormulae " = ["4-corr ppm 206*"] * ["4-corr 208Pb*/206Pb*"] * 208 / 206 ",
f, piaRadDauCol_48, L
--where 208 is the atomic mass of 208Pb and 206 is the atomic mass of 206Pb.
Case 5 --populate column [7-corr ppm 208*"], noting that this calculation draws
--on the row-by-row results determined in Case 2
PlaceFormulae " = ["7-corr ppm 206*"] * ["7-corr 208Pb*/206Pb*"] * 208 / 206 ",
f, piaRadDauCol_78, L
--where 208 is the atomic mass of 208Pb and 206 is the atomic mass of 206Pb.
End Select
Next p
End If
If column ["4corr 206*/238"] exists, calculate 204Pb-corrected 206Pb*/238U values (and their uncertainties), and place them in columns ["4corr 206*/238"] and ["4corr 206*/238 %err"], along with their corresponding dates, which are placed in columns ["204corr 206Pb/238U Age"] and ["204corr 206Pb/238U Age 1serr"]. Strictly, these calculations are predicated on the presence of the column ["204/206"], but this code permits the calculation to proceed even if ["204/206"] does not exist (albeit by instead calculating an ["Uncorr 206Pb/238U Age"]):
If piPb6U8_4col > 0 --i.e. if column ["4corr 206*/238"] exists
--Populate column ["4corr 206*/238"]:
PlaceFormulae = " = ["Total 206Pb/238U"] * ["radd6"] ", f, piPb6U8_4col, L
--Now populate column ["4corr 206*/238 %err"]. Ludwig has commented the expression:
--Var(Rad6/8) = Var(Tot6/8) + (sComm_64/(Alpha-sComm_64))^2 * Var(Alpha)
--and this seems to be correct in detail, although easier to read if rearranged:
--Var(Rad6/8) = Var(Tot6/8) + Var(Alpha) * ( (sComm_64/(Alpha-sComm_64))^2 )
--The expression for t1 (i.e. sqrt[Var(Rad6/8)]) is a bit easier to read:
t1 = " = SQRT( ["Total 206Pb/238U %err"]^2 + ( sComm_64 * ["204/206 %err"] /
(1 / ["204/206"] - sComm_64) )^2 )" _
PlaceFormulae t1, f, piPb6U8_4ecol, L
If piU8Pb6_4col > 0 Then ' 238U/206Pb* PlaceFormulae "=1/ Pb6U8_4 ", f, piU8Pb6_4col, L PlaceFormulae "= Pb6U8_4e ", f, piU8Pb6_4ecol, L End If
' 204-corr 206238 age PlaceFormulae "=LN(1+ Pb6U8_4 )/" & pscLm8, f, piAgePb6U8_4col, L ' 4-corr206/238 age-err t1 = fsRa(f, piAgePb6U8_4col) d1 = "(" & NetAlpha & ".01* Pb6U8_totE )^2" d3 = "(.01* Pb46eCol " & psaC64(0) & ")^2" d4 = "( Pb6U8_tot * Pb46 )^2" d5 = "(1/" & pscLm8 & "/Exp(" & pscLm8 & " AgePb6U8_4 ))^2" tmp = d5 & "" & d4 & "(" & d1 & "+" & d3 & ")" tmp = "=SQRT(" & tmp & ")" PlaceFormulae tmp, f, piAgePb6U8_4ecol, L
ElseIf piPb46col = 0 And piAgePb6U8_4col > 0 Then ' age from ucorr 206/238 PlaceFormulae "=LN(1+ Pb6U8_tot )/" & pscLm8, f, piAgePb6U8_4col, L CFs plHdrRw, piAgePb6U8_4col, "total|206Pb|/238U|age", True PlaceFormulae "= Pb6U8_tot /" & pscLm8 & "/(1+ Pb6U8_tot )* Pb6U8_totE /100", _ f, piAgePb6U8_4ecol, L End If ' If piPb6U8_4col > 0