-
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
pbCalc8corrConcPlotRats, pbTh, pbU, Switch.DirectAltPD
Values of type Integer
f, L, m, piNumDauPar
Values of type String
d2, t5
Values of type Double
Alpha, Beta, Gamma, d1, d3, d4, d5, NetAlpha, NetBeta, NetGamma, radd6, radd8, sComm_68, t1, t3, TotPb6U8AbsErr, TotPb76AbsErr, TotPb8Th2AbsErr
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'). Define row-by-row values of Alpha, Beta and Gamma:
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 each measured value to generate row-by-row values of NetAlpha, NetBeta and NetGamma:
NetAlpha = Alpha - sComm_64 --i.e. R64r = R64m - R64c
NetBeta = Beta - sComm_74 --i.e. R74r = R74m - R74c
NetGamma = Gamma - sComm_84 --i.e. R84r = R84m - R84c
For the daughter isotopes of calibration-related interest (i.e 206Pb and 208Pb), it is possible to define, row-by-row, proportionality factors describing the radiogenic component as a fraction of the total measured value. The calculated values of radd6 and radd8 are used later in this subroutine:
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 efficiently 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
--Place row-by-row expression for ["Total 206Pb/238U"]
PlaceFormulae " = ["Uncorr Pb/U const"] / WtdMeanA1 * StdUPbRatio ", f, m, L
--i.e. place expression in col m, rows f to L
--Place row-by-row expression for ["Total 206Pb/238U %err"]
PlaceFormulae " = SQRT( ["Uncorr Pb/U 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
--Place row-by-row expression for ["Total 206Pb/238U"]
PlaceFormulae " = ["Uncorr Pb/U const"] / WtdMeanA2 * StdUPbRatio ", f, m, L
--i.e. place expression in col m, rows f to L
--Place row-by-row expression for ["Total 206Pb/238U %err"]
PlaceFormulae " = SQRT( ["Uncorr Pb/U 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
--Place row-by-row expression for ["Total 208Pb/232Th"]
PlaceFormulae " = ["Uncorr Pb/Th const"] / WtdMeanA1 * StdThPbRatio ", f, m, L
--i.e. place expression in col m, rows f to L
--Place row-by-row expression for ["Total 208Pb/232Th %err"]
PlaceFormulae " = SQRT( ["Uncorr Pb/Th 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
--Place row-by-row expression for ["Total 208Pb/232Th"]
PlaceFormulae " = ["Uncorr Pb/Th const"] / WtdMeanA2 * StdThPbRatio ", f, m, L
--i.e. place expression in col m, rows f to L
--Place row-by-row expression for ["Total 208Pb/232Th %err"]
PlaceFormulae " = SQRT( ["Uncorr Pb/Th 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 = piPb8Th2_totCol
Next i
The next step in the SQUID 2.50 involves the calculation of ["Total 238U/206Pb"] by simple inversion. However, that calculation is premature, as its precursor ["Total 206Pb/238U"] has not yet been calculated for Perm3. Therefore, I have switched the order of those successive If loops, so the calculation of the "secondary" daughter-parent isotopic ratio (for the permutations where no calibration constant is available) happens first, and the inversion happens second:
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 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 206Pb/238U"] on SampleData:
PlaceFormulae " = ["Total 208Pb/232Th"] / ["208/206"] * ["232Th/238U"] ",
f, piPb6U8_totCol, L
--Fill column ["Total 206Pb/238U %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 calculate ["Total 238U/206Pb"], via simple inversion, everywhere ["Total 206Pb/238U"] exists:
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 206Pb/238U %err"] ", f, piU8Pb6_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 - ["7-corr 204Pb/206Pb"] * sComm_64 ) ", f, piaRadDauCol_76, L
Case 3 --populate column ["8-corr ppm 206*"]
PlaceFormulae " = ["Total 206Pb/238U"] * ["ppmU"] * 0.859 *
( 1 - ["8-corr 204Pb/206Pb"] * 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 ["Total 206Pb/238U Age"]). Squid3 could (and probably should) calculate BOTH ["204corr 206Pb/238U Age"] and ["Total 206Pb/238U Age"] as a matter of course, just because we can:
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 commented the expression:
--Var(Rad6/8) = Var(Tot6/8) + (sComm_64/(Alpha-sComm_64))^2 * Var(Alpha)
--which is correct in detail, although easier to read if rearranged:
--Var(Rad6/8) = Var(Tot6/8) + Var(Alpha) * ( (sComm_64/(Alpha-sComm_64))^2 )
--Defining ["4corr 206*/238 %err"] = sqrt[Var(Rad6/8)]:
PlaceFormulae " = SQRT( ["Total 206Pb/238U %err"]^2 + ( sComm_64 *
["204/206 %err"] / ( 1 / ["204/206"] - sComm_64) )^2 ) ", f, piPb6U8_4ecol, L
If piU8Pb6_4col > 0 --i.e. if column ["4corr 238/206*"] exists
--Populate column ["4corr 238/206*"] via simple inversion:
PlaceFormulae " = 1 / ["4corr 206*/238"] ", f, piU8Pb6_4col, L
--Populate column ["4corr 238/206* %err"]:
PlaceFormulae "= ["4corr 206*/238 %err"] ", f, piU8Pb6_4ecol, L
End If
--Now populate column ["204corr 206Pb/238U Age"]:
PlaceFormulae " = LN( 1 + ["4corr 206*/238"] ) / Lambda238Ma ",
f, piAgePb6U8_4col, L
--where Lambda238Ma has units "Ma^-1", so Age will be in Ma.
--Now calculate and populate column ["204corr 206Pb/238U Age 1serr"]:
d1 = ( NetAlpha * ["Total 206Pb/238U %err"] / 100 )^2
d3 = ( ["204/206 %err"] * sComm_64 / 100 )^2
d4 = ( ["Total 206Pb/238U"] * ["204/206"] )^2
d5 = ( 1 / Lambda238Ma /
EXP( Lambda238Ma * ["204corr 206Pb/238U Age"] ) )^2
--where Lambda238Ma has units "Ma^-1".
PlaceFormulae " = SQRT( d5 * d4 * (d1 + d3) ) ", f, piAgePb6U8_4ecol, L
ElseIf (piPb46col = 0) AND (piAgePb6U8_4col > 0) --this curious ElseIf
--covers the possibility that a column ["204corr 206Pb/238U Age"] might
--exist even in the absence of column ["204/206"], which would
--ordinarily be a prerequisite for the calculation. In this scenario,
--Ludwig simply calculates the Age from the **uncorrected** 206/238,
--and modifies the relevant column-header accordingly.
--Populate the column **designated** ["204corr 206Pb/238U Age"] with
--an expression **reflecting** ["Total 206Pb/238U Age"]:
PlaceFormulae " = LN( 1 + ["Total 206Pb/238U"] ) / Lambda238Ma,
f, piAgePb6U8_4col, L
--where Lambda238Ma has units "Ma^-1", so Age will be in Ma.
--...and then change its header to ["Total 206Pb/238U Age"]!!
CFs plHdrRw, piAgePb6U8_4col, "total|206Pb|/238U|age", TRUE
--Now populate column **designated** ["204corr 206Pb/238U Age 1serr"]
--with an expression **reflecting** ["Total 206Pb/238U Age 1serr"]:
PlaceFormulae " = ["Total 206Pb/238U"] / Lambda238Ma /
(1 + ["Total 206Pb/238U"] ) * ["Total 206Pb/238U %err"] / 100 ",
f, piAgePb6U8_4ecol, L
--where Lambda238Ma has units "Ma^-1", so Age 1serr will be in Ma.
End If --piPb6U8_4col > 0 --i.e. if column ["4corr 206*/238"] exists.
The next step is to perform the superset of calculations enabled by the measurement of ["207/206"]. These start with 4corr 207*/206* analogous to the 206*/238 calculations described above, but also encompass 4corr 207*/235, plus 207Pb-corrected 206*/238. The code continues:
If piPb76col > 0 --i.e. if column ["207/206"] exists
If piU8Pb6_totCol > 0 --i.e. if column ["Total 238U/206Pb"] exists
--Place ["Total 207Pb/206Pb"]; used for Tera-Wasserburg plots and calculations
PlaceFormulae "= ["207/206"] ", f, piPb76_totCol, L
--Place ["Total 207Pb/206Pb %err"]
PlaceFormulae " = ["207/206 %err"] ", f, piPb76_totEcol, L
End If --column ["Total 238U/206Pb"] exists
If (piPb46col > 0) AND (piU8Pb6_4col > 0) --i.e. if BOTH columns
--["204/206"] and ["4corr 238/206*"] exist
--Populate column ["4corr 207*/206*"]
PlaceFormulae " ABS( NetBeta / NetAlpha ") ",
f, piPb76_4col, L
--Populate column ["4corr 207*/206* %err"]
--Complicated expression; note that it uses ["4corr 207*/206*"]:
t1 = ( ( ["207/206"] - ["4corr 207*/206*"] ) * ["204/206 %err"]
/ 100 / ["204/206"] )^2
t3 = ( ["207/206 %err"] / ["204/206"] / 100 * ["207/206"] )^2
PlaceFormulae " = ABS( SQRT( t1 + t3 ) / NetAlpha * 100 / ["4corr 207*/206*"] ) ",
f, piPb76_4eCol, L
If piaAgePb76_4Col > 0 --i.e. if column ["204corr 207Pb/206Pb Age"] exists,
--the form of the next calculations depend on whether column ["204/206"]
--exists, analogous to the previous calculation of 206*/238 (i.e. if **NO**
--["204/206"] exists, a ["Total 207Pb/206Pb Age"] will be calculated; else
--["204corr 207Pb/206Pb Age"] will be calculated as would be expected.
If piPb46col > 0 --i.e. if column ["204/206"] exists
--Place LudwigLibrary expression for ["204corr 207Pb/206Pb Age"]:
PlaceFormulae " = AgePb76( ["4corr 207*/206*"] ) ",
f, piaAgePb76_4Col, L
--Place LudwigLibrary expression for ["204corr 207Pb/206Pb Age 1serr"]:
PlaceFormulae " = AgeErPb76( ["4corr 207*/206*"],["4corr 207*/206* %err"]
,,,, TRUE ) ", f, piaAgePb76_4eCol, L
Else --i.e. if column ["204/206"] does *NOT* exist. Note that this Else
--appears to be redundant, as this If loop cannot be traversed unless BOTH
--["204/206"] and ["4corr 238/206*"] exist, as per If clause near the top
--of this block:
--Place LudwigLibrary expression for ["Total 207Pb/206Pb Age"]
--in the column **designated** ["204corr 207Pb/206Pb Age"]...
PlaceFormulae " = AgePb76( ["Total 207/206"] ) ",
f, piaAgePb76_4Col, L
--...and then change the column-header to reflect the contents!
CFs plHdrRw, piaAgePb76_4Col, ["Total 207Pb/206Pb Age"], TRUE
--Place LudwigLibrary expression for ["Total 207Pb/206Pb Age 1serr"]
--in the column **designated** ["204corr 207Pb/206Pb Age 1serr"]:
PlaceFormulae " = AgeErPb76( ["Total 207*/206*"],["Total 207*/206* %err"]
,,,,TRUE ) ", f, piaAgePb76_4eCol, L
End If --piPb46col > 0 --i.e. if column ["204/206"] exists
End If --piaAgePb76_4Col > 0 --i.e. if ["204corr 207Pb/206Pb Age"] exists
--Now place expressions for ["4corr 207*/235"] and ["4corr 207*/235 %err"]:
If (piPb7U5_4col > 0) AND (piPb6U8_4col > 0) --i.e. if BOTH columns
--["4corr 207*/235"] and ["4corr 206*/238"] exist
--Calculate and place ["4corr 207*/235"]:
PlaceFormulae " = ["4corr 207*/206*"] * ["4corr 206*/238"] * Present238U235U ",
f, piPb7U5_4col, L
--Calculate and place ["4corr 207*/235 %err"]:
PlaceFormulae " = SQRT( ["4corr 207*/206* %err"]^2 +
["4corr 206*/238 %err"]^2 ), f, piPb7U5_4ecol, L
--Calculate and place rho (error correlation) value ["4corr err corr"]:
PlaceFormulae " = ["4corr 206*/238 %err"] / ["4corr 207*/235 %err"] ",
f, piPb7U5Pb6U8_4rhoCol, L
End If --BOTH columns ["4corr 207*/235"] and ["4corr 206*/238"] exist
End If --BOTH columns ["204/206"] and ["4corr 238/206*"] exist
--Now calculate and place ["204corr Discordance"], noting that SQUID 2.50
--column-headers do not specify the index isotope prefix:
If (piDiscordCol > 0) AND (piPb46col > 0) --i.e. if BOTH columns
--["204corr Discordance"] and ["204/206"] exist
--SQUID 2.50 defines Discordance as 100 * ( 1 - R68m / R68i ) where
--R68i is the concordant 206*/238 RATIO (corresponding to the measured
--207*/206* RATIO) and R68m is the measured 206*/238 RATIO. First calculate
--the row-by-row values of (204corr) R68i:
R68i = EXP( Lambda238Ma * ["204corr 207Pb/206Pb Age"] ) - 1
--Then place the expression for ["204corr Discordance"]:
PlaceFormulae " = 100 * ( 1 - ["4corr 206*/238"] / R68i ) ",
f, piDiscordCol, L
End If --BOTH columns ["204corr Discordance"] and ["204/206"] exist
--Now calculate and place ["207corr 206Pb/238U Age"] and its uncertainty.
--Note that this code does *NOT* encompass the ["7corr 206*/238"] RATIO, but it
--SHOULD; the relevant code is misplaced a little further along. Below I have
--marked the spot at which ["7corr 206*/238"] and ["7corr 206*/238 %err"]
--SHOULD be evaluated!
If piAgePb6U8_7col > 0 --i.e. if column ["207corr 206Pb/238U Age"] exists
--Place LudwigLibrary expression for ["207corr 206Pb/238U Age"]:
PlaceFormulae " = Age7corr( ["Total 206Pb/238U"], ["207/206"], sComm_76 ) ",
f, piAgePb6U8_7col, L
--Place LudwigLibrary expression for ["207corr 206Pb/238U Age 1serr"], noting
--that Ludwig function AgeEr7corr has a total of seven arguments. Two of these
--arguments are *ABSOLUTE* uncertainties, which I have evaluated in advance,
--for clarity. Calculate TotPb6U8AbsErr and TotPb76AbsErr as:
TotPb6U8AbsErr = ["Total 206Pb/238U %err"] / 100 * ["Total 206Pb/238U"]
TotPb76AbsErr = ["Total 207Pb/206Pb %err"] / 100 * ["Total 207Pb/206Pb"]
PlaceFormulae " = AgeEr 7corr( ["207corr 206Pb/238U Age"],
["Total 206Pb/238U"], TotPb6U8AbsErr, ["207/206"], TotPb76AbsErr,
sComm_76, 0 ) ", f, piAgePb6U8_7ecol, L
--This is the spot at which ["7corr 206*/238"] and ["7corr 206*/238 %err"]
--SHOULD be evaluated! The relevant lines are marked in the next code-block.
End If --piAgePb6U8_7col > 0 --i.e. if ["207corr 206Pb/238U Age"] exists
End If --piPb76col > 0 --i.e. if column ["207/206"] exists
The next step is to calculate Xcorr 208Pb/232Th ratios, dates and uncertainties (where X = 4 or 7). The 204Pb-corrected calculations are analogous to those for 206*/238 and 207*/206*. Note also that the following code contains an out-of-place block relating to ["7corr 206*/238"] and ["7corr 206*/238 %err"]. This is undesirable at two levels: (1) it should have been evaluated in the previous code-block, subject to the If conditions that are relevant to the calculation, and (2) The If conditions employed below are not appropriate to the calculation: they could result in the evaluation being skipped for no good reason. The code continues:
If (piPb8Th2_totCol > 0) AND (piPb86col > 0) --i.e. if BOTH columns
--["Total 208Pb/232Th"] and ["208/206"] exist
If piPb46col > 0 --i.e. if column ["204/206"] exists
--Place the expression for ["4corr 208*/232"]:
PlaceFormulae " = ["Total 208Pb/232Th"] * radd8 ", f, piPb8Th2_4col, L
--Place the expression for ["4corr 208*/232 %err"]. Ludwig notes that this
--expression "neglects the 208/206 error":
PlaceFormulae " = SQRT( ["Total 208Pb/232Th %err"]^2 +
( sComm_84 / NetGamma )^2 * ["204/206 %err"]^2 )", f, piPb8Th2_4eCol, L
End If
If piPb46col > 0 --i.e. if column ["204/206"] exists
--Place expression for ["204corr 208Pb/232Th Age"]:
PlaceFormulae " = LN( 1 + ["4corr 208*/232"] ) / Lambda232Ma ",
f, piAgePb8Th2_4col, L
--where Lambda232Ma has units "Ma^-1", so Age will be in Ma.
--Place expression for ["204corr 208Pb/232Th Age 1serr"]:
PlaceFormulae " = ["4corr 208*/232"] / Lambda232Ma /
(1 + ["4corr 208*/232"] ) * ["4corr 208*/232 %err"] / 100 ",
f, piAgePb8Th2_4eCol, L
--where Lambda232Ma has units "Ma^-1", so Age 1serr will be in Ma.
Else --i.e. if column ["204/206"] does *NOT* exist
--Place expression for ["Total 208Pb/232Th Age"] in the column
--**designated** ["204corr 208Pb/232Th Age"]...
PlaceFormulae " = LN( 1 + ["Total 208Pb/232Th"] ) / Lambda232Ma ",
f, piAgePb8Th2_4col, L
--...and then change the column-header to reflect the contents!
CFs plHdrRw, piAgePb8Th2_4col, ["Total 208Pb/232Th Age"], TRUE
--Place expression for ["Total 208Pb/232Th Age 1serr"] in the column
--**designated** ["204corr 208Pb/232Th Age 1serr"]:
PlaceFormulae " = ["Total 208Pb/232Th"] ) / Lambda232Ma /
(1 + ["Total 208Pb/232Th"] ) * ["Total 208Pb/232Th %err"] / 100 ",
f, piAgePb8Th2_4eCol, L
End If --piPb46col > 0 --i.e. if column ["204/206"] exists
--Ludwig's next If block is superfluous: it repeats the evaluation of
--["4corr 208*/232 %err"], adding an unnecessary Boolean, and an absolute
--value of NetGamma, which is moot as the NetGamma-bearing term is squared
--anyway. So I have commented-out the following lines:
--
--If (pbTh = TRUE) AND (piPb46col > 0) --i.e. if (Perm3 or Perm4) AND column
----["204/206"] exists:
--
-- PlaceFormulae " = SQRT( ["Total 208Pb/232Th %err"]^2 +
-- ( sComm_84 / ABS( NetGamma ) )^2 * ["204/206 %err"]^2 )",
-- f, piPb8Th2_4eCol, L
--
--End If
--Ludwig's next If block is OUT OF PLACE; it should have appeared at the base
--of the previous code-block, in the location I have marked there. I have
--included the code below, but I have inserted some comments of my own.
--**NOTE ALSO** that Ludwig's code below also contained an arithmetical bug:
--he used the wrong Lambda at one point, which results the value being too small
--by a factor of about 3!!
If piAgePb6U8_7col > 0 --i.e. if column ["207corr 206Pb/238U Age"] exists. This
--criterion on its own is appropriate for the following calculation, and it is
--also consistent with the If criterion already in place at the proper
--destination (shown in the previous code-block) for the following expressions.
--But note that this criterion is **NOT** appropriate in combination with the
--2-column If criteria applied at the top of the current code-block. As written,
--the following calculation would not be performed if ["208/206"] was not
--measured, which is ridiculous.
--CODE TO BE TRANSPLANTED TO PREVIOUS CODE-BLOCK STARTS **HERE**
--Place expression for ["7corr 206*/238"]:
PlaceFormulae " = EXP ( Lambda238Ma * ["207corr 206Pb/238U Age"] ) - 1 ",
f, piPb6U8_7col, L
--where Lambda238Ma has units "Ma^-1".
--Place expression for ["7corr 206*/238 %err"]. Ludwig's code contains a bug
--in the following line (commented out), where his expression incorrectly
--commences with Lambda232Ma:
--PlaceFormulae " = Lambda232Ma * EXP( Lambda238Ma * ["207corr 206Pb/238U Age"] ) *
-- ["207corr 206Pb/238U Age 1serr"] / ["7corr 206*/238"] * 100 ",
-- f, piPb6U8_7ecol, L
--The corrected expression commences with Lambda238Ma:
PlaceFormulae " = Lambda238Ma * EXP( Lambda238Ma * ["207corr 206Pb/238U Age"] ) *
["207corr 206Pb/238U Age 1serr"] / ["7corr 206*/238"] * 100 ", f, piPb6U8_7ecol, L
--CODE TO BE TRANSPLANTED TO PREVIOUS CODE-BLOCK ENDS **HERE**
End If --piAgePb6U8_7col > 0 i.e. if column ["207corr 206Pb/238U Age"] exists.
--Code transplantation should result in abolition of the previous If.
--Now perform calculations related to ["7corr 208*/232"]:
If piAgePb8Th2_7col > 0 --i.e. if column ["207corr 208Pb/232Th Age"] exists:
--Place LudwigLibrary expression for ["207corr 208Pb/232Th Age"], which has a
--total of 7 arguments:
PlaceFormulae " = Age7corrPb8Th2( ["Total 206Pb/238"], ["Total 208Pb/232Th"],
["208/206"], ["207/206"], sComm_64, sComm_76, sComm_86 ) ",
f, piAgePb8Th2_7col, L
--Place LudwigLibrary expression for ["207corr 208Pb/232Th Age 1serr"], which
--has a total of 11 arguments (i.e. includes %err values for the first four
--arguments of the previous function). Note also that the input order of
--["207/206"] and ["208/206"] is *REVERSED* in this second function!
PlaceFormulae " = AgeErr7corrPb8Th2( ["Total 206Pb/238"],
["Total 206Pb/238 %err"],["Total 208Pb/232Th"], ["Total 208Pb/232Th %err"],
["207/206"], ["207/206 %err"], ["208/206"], ["208/206 %err"], sComm_64,
sComm_76, sComm_86 ) ", f, piAgePb8Th2_7ecol, L
--Place expression for ["7corr 208*/232"]:
PlaceFormulae " = EXP ( Lambda232Ma * ["207corr 208Pb/232Th Age"] ) - 1 ",
f, piPb8Th2_7col, L
--where Lambda232Ma has units "Ma^-1".
--Place expression for ["7corr 208*/232 %err"]:
PlaceFormulae " = Lambda232Ma * EXP( Lambda232Ma * ["207corr 208Pb/232Th Age"] ) *
["207corr 208Pb/232Th Age 1serr"] / ["7corr 208*/232"] * 100 ",
f, piPb8Th2_7ecol, L
End If --piAgePb8Th2_7col > 0 i.e. if column ["207corr 208Pb/232Th Age"] exists
The next thing SQUID 2.50 does is calculate ["208corr 206Pb/238U Age"] and its uncertainty. As for the 207corr equivalent, this is a hangover from SQUID 1, where these two ages and their errors were presented without the actual source ratios! SQID 2.50 rectifies that, but only "in reverse" (i.e. it still calculates the Age first, and then derives the ratios from the Ages). Note that the following If condition seems irrelevant at best (incorrect at worst), as neither ["204/206"] nor any of its dependents are required in order to perform the LudwigLibrary calculations 'Age8Corr' or 'AgeEr8Corr'. I think it is possible that the condition should read "If piPb86col > 0" (i.e. if column ["208/206"] exists), because that would make more sense, but note that piPb86col > 0 is already a precondition of the entire loop, and so would be redundant here.
After some deliberation, I have left the If condition (piPb4col > 0) as it was specified by Ludwig. One possibility is that the If condition should be removed entirely; perhaps these particular calculations should be mandatory if the foregoing If conditions have been satisfied.
For use as arguments in ensuing LudwigLibrary functions, first define sComm_68 as the inverse of sComm_86, and TotPb8Th2AbsErr as the Th-Pb equivalent of the previously defined TotPb6U8AbsErr, and then proceed:
sComm_68 = 1 / sComm_86
TotPb8Th2AbsErr = ["Total 208Pb/232Th %err"] / 100 * ["Total 208Pb/232Th"]
--If piPb46col > 0 --i.e. if ["204/206"] exists, although this doesn't make sense!
--Place LudwigLibrary expression for ["208corr 206Pb/238U Age"]:
PlaceFormulae " = Age8Corr( ["Total 206Pb/238U"], ["Total 208Pb/232Th"],
["232Th/238U"], sComm_68 ) ", f, piAgePb6U8_8col, L
--Place LudwigLibrary expression for ["208corr 206Pb/238U Age 1serr"], noting
--that the function has a total of 9 arguments:
PlaceFormulae " = AgeEr8Corr( ["208corr 206Pb/238U Age"], ["Total 206Pb/238U"],
TotPb6U8AbsErr, ["Total 208Pb/232Th"], TotPb8Th2AbsErr, ["232Th/238U"], 0,
sComm_68, 0 ) ", f, piAgePb6U8_8ecol, L
End If --piPb46col > 0 --i.e. if ["204/206"] exists.
SQUID 2.50 has a user-controlled Boolean (pbCalc8corrConcPlotRats) which controls whether the 208corr isotopic ratios required for Wetherill and Tera-Wasserburg concordia plots are calculated (possibly because they frequently give results that are ridiculous or impossible in the case of high-Th analyses). But in Squid3, I would think that errors in calculations would be handled separately, and that this Boolean ought to at least be TRUE by default, if not hard-coded TRUE. Then, if column ["208corr 206Pb/238U Age"] exists, the code continues:
If (pbCalc8corrConcPlotRats = TRUE) AND (piAgePb6U8_8col > 0)
--i.e. if user has requested 208corr Concordia ratios AND the column
--["208corr 206Pb/238U Age"] exists:
If piU8Pb6_8col > 0 --i.e. if column ["8corr 238/206*"] exists
--PLace LudwigLibrary expression for ["8corr 206*/238"] first:
PlaceFormulae " = Pb206U238rad( ["208corr 206Pb/238U Age"] ) ",
f, piPb6U8_8col, L
--Then place expression for ["8corr 206*/238 %err"]:
PlaceFormulae " = Lambda238Ma * ( 1 + ["8corr 206*/238"] ) *
["208corr 206Pb/238U Age 1serr"] * 100 / ["8corr 206*/238"] ",
f, piPb6U8_8ecol, L
--Then place expression for ["8corr 238/206*"] by simple inversion:
PlaceFormulae " = 1 / ["8corr 206*/238"] ", f, piU8Pb6_8col, L
--Finally, ["8corr 238/206* %err"] is equivalent to ["8corr 206*/238 %err"]:
PlaceFormulae " = ["8corr 206*/238 %err"] ", f, piU8Pb6_8ecol, L
End If --piU8Pb6_8col > 0 --i.e. if column ["8corr 238/206*"] exists
If piPb7U5_8col > 0 --i.e. if column ["8corr 207*/235"] exists
--Place LudwigLibrary expression for ["8corr 207*/235"]:
PlaceFormulae " = Rad8corPb7U5( ["208corr 206Pb/238U Age"], ["Total 206Pb/238U"],
["207/206"], sComm_76 ) ", f, piPb7U5_8col, L
--Define row-by-row TotPb7U5 as follows (noting that SQUID 2.50 hard-codes a
--numeric value of 137.88 in lieu of physical constant Present238U235U):
TotPb7U5 = ["Total 206Pb/238U"] * ["207/206"] / Present238U235U
--Now place LudwigLibrary expression for [8corr 207*/235 %err"], noting that
--the function has a total of 12 arguments:
PlaceFormulae " = Rad8corPb7U5perr( ["Total 206Pb/238U"],
["Total 206Pb/238U %err"], ["8corr 206*/238"], TotPb7U5, ["232Th/238U"],
["232Th/238U %err"], ["207/206"], ["207/206 %err"], ["208/206"],
["208/206 %err"], sComm_76, sComm_86 ) ", f, piPb7U5_8ecol, L
--Now place LudwigLibrary expression for [8corr err corr"], noting that the
--function has a total of 11 arguments (as per previous, without TotPb7U5):
PlaceFormulae " = Rad8corConcRho( ["Total 206Pb/238U"], ["Total 206Pb/238U %err"],
["8corr 206*/238"], ["232Th/238U"], ["232Th/238U %err"], ["207/206"],
["207/206 %err"], ["208/206"], ["208/206 %err"], sComm_76, sComm_86 ) ",
f, piPb7U5Pb6U8_8rhoCol, L
End If --piPb7U5_8col > 0 --i.e. if column ["8corr 207*/235"] exists
The final step in this loop is to evaluate ["8corr 207*/206*"] and its uncertainty, and then use those to generate ["208corr 207Pb/206Pb Age"] and its error:
If piPb76_8col > 0 --i.e. if column ["8corr 207*/206*"] exists
--Place expression for ["8corr 207*/206*"], noting that once again SQUID 2.50
--uses the hard-coded numeric value 137.88 in lieu of Present238U235U::
PlaceFormulae " = ["8corr 207*/235"] / ["8corr 206*/238"] / Present238U235U ",
f, piPb76_8col, L
--Place expression for ["8corr 207*/206* %err"]:
PlaceFormulae " = SQRT( ["8corr 207*/235 %err"]^2 + ["8corr 206*/238 %err"]^2 -
2 * ["8corr 207*/235 %err"] * ["8corr 206*/238 %err"] * ["8corr err corr"] ) ",
f, piPb76_8ecol, L
End If --piPb76_8col > 0 --i.e. if column ["8corr 207*/206*"] exists
If piAgePb76_8col > 0 --i.e. if column ["208corr 207Pb/206Pb Age"] exists
--Place LudwigLibrary expression for ["208corr 207Pb/206Pb Age"]:
PlaceFormulae " = AgePb76( ["8corr 207*/206*"] ) ", f, piAgePb76_8col, L
--Place LudwigLibrary expression for ["208corr 207Pb/206Pb Age 1serr"]:
PlaceFormulae " = AgeErPb76( ["8corr 207*/206*"], ["8corr 207*/206* %err"]
,,,, TRUE ) ", f, piAgePb76_8ecol, L
--PLACE EXPRESSION FOR 208CORR DISCORDANCE (SQUID 2.50 DOES NOT EVALUATE IT)
--INSERT IDENTICAL EXPRESSION TO 204CORR DISCORDANCE, BUIT WITH 8CORR RATIOS
--AND AGES SUBSTITUTED!
End If --piAgePb76_8col > 0 --i.e. column ["208corr 207Pb/206Pb Age"] exists
End If --(pbCalc8corrConcPlotRats = TRUE) AND (piAgePb6U8_8col > 0)
--i.e. if user has requested 208corr Concordia ratios AND the column
--["208corr 206Pb/238U Age"] exists
End If --If (piPb8Th2_totCol > 0) AND (piPb86col > 0) --i.e. if BOTH columns
--["Total 208Pb/232Th"] and ["208/206"] exist
--The subroutine finishes with some formatting to colour the rows corresponding to
--the Uranium concentration reference material pale yellow.
End Sub