-
Notifications
You must be signed in to change notification settings - Fork 16
SHRIMP: Sub OverCtMeans
This subroutine (which is solely for the Standard) does a bit more than the name implies. Firstly it places, row-by-row, formulae to calculate the 204-corrected 207Pb/206Pb ratio and its 1sigma percentage uncertainty, as well as invoking (row-by-row) the relevant LudwigLibrary functions to calculate the associated 204-corrected 207Pb/206Pb date and its 1sigma absolute uncertainty.
Secondly, it identifies which columns can usefully have robust means calculated, performs those calculations (using LudwigLibrary function TukeysBiweight, with tuning 9) and places the output of the expression as a 3 x 1 array beneath the relevant column. The SQUID 2.50 subroutine requires the index number of the last row of analytical data as an input, so it can determine in which rows the "summary" results should be placed so that the calculated biweights appear directly beneath the input data.
OverCtMeans plaLastDatRw
plaLastDatRw: Integer index number of the last row containing spot-by-spot data (for the Standard).
Values of type Boolean
ShowOverCtCols
Values of type Integer
c, k, kk, plaFirstDatRw
Values of type String
t0, t1, t2, t3, t4, t5
The subroutine starts by calculating the 204Pb-corrected 207Pb/206Pb, its uncertainty, and the associated date and uncertainty.
If piaAgePb76_4Col[1] > 0 --i.e. if column ["4-corr207Pb/206Pb"] exists on StandardData:
t0 = "=( ["207/206"] / ["204/206"] - sComm_74 ) / ( 1 / ["204/206"] - sComm_64 )"
--Then place this formula in column ["4-corr207Pb/206Pb"] on StandardData:
PlaceFormulae t0, plaFirstDatRw, piStdPb76_4Col, plaLastDatRw
--Now calculate the associated %err:
t1 = "( ( ["207/206"] * ["207/206%err"] )^2 +
( ["204/206"] * ( ["4-corr207Pb/206Pb"] * sComm_64 - sComm_74) * ["204/206%err"] )^2 )"
t2 = "( ["207/206"] - ["204/206"] * sComm_74 )^2"
t3 = "=sqrt( t1 / t2 )"
--Then place this formula in column ["4-corr207Pb/206Pb%err"] on StandardData:
PlaceFormulae t3, plaFirstDatRw, piStdPb76_4eCol, plaLastDatRw
--Now invoke LudwigLibrary functions to calculate associated age and error:
t4 = "=AgePb76( ["4-corr207Pb/206Pb"] )"
--Then place this formula in column ["4-corr207Pb/206Pbage"] on StandardData:
PlaceFormulae t4, plaFirstDatRw, piaAgePb76_4Col, plaLastDatRw
--Convert ["4-corr207Pb/206Pb%err"] to absolute, to calculate age error:
t5 = "=AgeErPb76( ["4-corr207Pb/206Pb"] , ["4-corr207Pb/206Pb"] * ["4-corr207Pb/206Pb%err"] / 100 )"
--Then place this formula in column ["4-corr207Pb/206Pbage±1sigma"] on StandardData:
PlaceFormulae t5, plaFirstDatRw, piaAgePb76_4eCol, plaLastDatRw
End If
The second part of the subroutine calculates the various biweight means, for the range of columns for which they are relevant. In practice, the "entry" Boolean is always TRUE, because Ludwig set the nominally user-defined Boolean ShowOverCtCols to TRUE, and then removed the relevant check-box from the user form! Note that a possibly unintended consequence of "locking in" this Boolean is that it formally requires *all SQUID 2.50 U-Pb Geochronology Tasks to include 204Pb in the list of mass-stations.
__(There will come a time when it is necessary for us to formally assess the absolute minimum mass-stations in order to perform U-Pb (or Th-Pb) geochronology in SQUID 3.0. My feeling is that the absolute minimum lists are as follows:
- 206Pb/238U: 206Pb, 238U (or a proxy thereof), and ONE OF (204Pb or 207Pb)
- 208Pb/232Th: 208Pb, 232Th (or a proxy thereof), and ONE OF (204Pb or 207Pb)
At present, SQUID 2.50 enforces "BOTH OF". Certainly the "alternate" daughter-isotope (i.e. 208Pb for U-Pb, 206Pb for Th-Pb) should be optional, and probably Background should be optional too (even though including 204Pb in a run-table without including Background would be terrible practice. Something to revisit later, and the stakes are relatively low: it doesn't really matter exactly what the 'bare bones' list of mass-stations is, as long as we warn SQUID 3.0 users up-front... it would be an improvement on SQUID 2.50!)__
The subroutine proceeds:
If ShowOverCtCols = TRUE OR piaAgePb76_4Col[1] > 0 --i.e. if StandardData contains ["4-corr207Pb/206Pbage"]
If piaOverCts4Col[7] > 0 Or piaOverCts4Col[8] > 0 Or piaAgePb76_4Col[1] > 0
--i.e. if ANY of the columns ["204overcts/sec(fr. 207)"], ["204overcts/sec(fr. 208)"],
--["4-corr207Pb/206Pbage"] exist on the StandardData sheet, then define the extent of calculations
--to be performed:
If ShowOverCtCols = TRUE
kk = 1
Else
kk = 5
End If
For k = kk To 5
Select Case k --all column-indices refer to StandardData sheet:
Case 1: c = piaOverCts4Col[7] --index for ["204overcts/sec(fr. 207)"]
Case 2: c = piaOverCts4Col[8] --index for ["204overcts/sec(fr. 208)"]
Case 3: c = piacorrAdeltCol[7] --index for ["7-corr206Pb/238Uconst.delta%"]
Case 4: c = piacorrAdeltCol[8] --index for ["8-corr206Pb/238Uconst.delta%"]
Case 5: c = piaAgePb76_4Col[1] --index for ["4-corr207Pb/206Pbage"]
End Select
If c > 0
--Define Range bw as 3 x 1, with the first row immediately following analytical data:
Set bw = frSr(1 + plaLastDatRw, c, 3 + plaLastDatRw) --sets Range bw = 3 rows by 1 col
--Invoke LudwigLibrary function Biweight (tuning 9) on the data in the column above bw.
--Note that for k = 3 or 4, the values in this range have not yet been calculated!
bw.FormulaArray = "=Biweight(" & frSr(plaFirstDatRw, c, plaLastDatRw).Address & ",9)"
--Finally, add Names (and explanatory Notes) to result-cells. These are relevant because
--the Names can be used in Task-expressions, and the Notes matter because the columns are
--quite abstract concepts for non-SHRIMP geochronologists (as well as SHRIMP beginners):
--Note that I have recast the following If to make it longer but more transparent:
If k = 1
t1 = "StandardData!Pb204OverCts7corr"
AddName t1, TRUE, 1 + plaLastDatRw, c --Name for the Biweight Value element
t2 = "StandardData!Pb204OverCts7corrEr"
AddName t2, TRUE, 3 + plaLastDatRw, c --Name for the Biweight 95% conf. element
t3 = "Robust avg 204 overcts assuming 206Pb/238U-207Pb/235U age concordance"
ElseIf k = 2
t1 = "StandardData!Pb204OverCts8corr"
AddName t1, TRUE, 1 + plaLastDatRw, c --Name for the Biweight Value element
t2 = "StandardData!Pb204OverCts8corrEr"
AddName t2, TRUE, 3 + plaLastDatRw, c --Name for the Biweight 95% conf. element
t3 = "Robust avg 204 overcts assuming 206Pb/238U-208Pb/232Th age concordance"
ElseIf k = 3
t1 = "OverCtsDeltaP7corr"
AddName t1, TRUE, 1 + plaLastDatRw, c --Name for the Biweight Value element
t2 = "OverCtsDeltaP7corrEr"
AddName t2, TRUE, 3 + plaLastDatRw, c --Name for the Biweight 95% conf. element
t3 = "Robust avg of diff. between 207-corr. and 204-corr. calibr. const."
ElseIf k = 4
t1 = "OverCtsDeltaP8corr"
AddName t1, TRUE, 1 + plaLastDatRw, c --Name for the Biweight Value element
t2 = "OverCtsDeltaP8corrEr"
AddName t2, TRUE, 3 + plaLastDatRw, c --Name for the Biweight 95% conf. element
t3 = "Robust avg of diff. between 208-corr. and 204-corr. calibr. const."
Else --interestingly, SQUID 2.50 does not assign Names (t1, t2) to the Biweight
--elements calculated for the ["4-corr207Pb/206Pbage"] column. Probably an oversight!
t3 = "Robust average of 204-corrected 207/206 age"
End If --k = 1 to 5
--Now VBA function "Note" to add text t3 as Comment on the Biweight Value element:
Note 1 + plaLastDatRw, c, t3
--Corresponding Comment on the Biweight 95% conf. element is context-sensitive:
If k = 3 OR k = 4
t3 = "95%-conf. error in above difference"
Else
t3 = "95%-conf. error in above"
End If
Note 3 + plaLastDatRw, c, t3
End If --c > 0
Next k
End If --piaOverCts4Col[7] > 0 Or piaOverCts4Col[8] > 0 Or piaAgePb76_4Col[1] > 0
ClearObj bw
End If --ShowOverCtCols = TRUE OR piaAgePb76_4Col[1] > 0
End Sub