-
Notifications
You must be signed in to change notification settings - Fork 16
SHRIMP: Sub WtdMeanAcalc
This subroutine evaluates, for the Standard, the weighted mean (and associated parameters) of each relevant set of common-Pb corrected calibration constant values. This mean is the value to which all spot-by-spot calibration constants determined for the unknowns will be calibrated, and all spot-by-spot, 'directly-calculated daughter/parent dates' (i.e. 206Pb/238U and/or 208Pb/232Th as appropriate) are calculated from there, because SHRIMP (and secondary ion mass spectrometry in general) is an indirect dating technique. Thus this is a critical step.
WtdMeanAcalc BadSbm(), Adrift(), AdriftErr()
BadSbm: Array containing index numbers of spots for which 'bad' SBM values exist (e.g values lower than SBMzero).
Adrift: Array of double-precision values arising from 'correction' of the measured calibration-constant values for secular drift associated with the duration of the analytical session (i.e. "Hours"). Currently out of scope.
AdriftErr: Array of double-precision uncertainties associated with values arising from 'correction' of the measured calibration-constant values for secular drift associated with the duration of the analytical session (i.e. "Hours"). Currently out of scope.
Values of type Boolean
NoReject, NoUPbConstAutoReject, pbCanDriftCorr, pbU, pbTh, [Task]Switch.DirectAltPD
Values of type Integer
c,DauParNum, ErN, h, i, j, LargeErRegN, Npts, Nrej, NtotRej, NumDauPar, r, rw1, rwn
Values of type Double
ExtPerr, ExtPtSigma, Lambda, MedianEr, MSWD, Nmadd, Prob, pscLm2, pscLm8, StdThPbRatio, StdUPbRatio, WtdMean, WtdMeanErr
Values of type String
s, Ele, t1, t2, t4, t5
Vectors of type Integer
LargeErRej, Rejected
Vectors of type Double
Acol, Aecol, Arange, Aerange, ErrVals
Arrays of type Double
Adatrange
Arrays of type Variant/Mixed
wW
The function commences with some setting up. NoUPbConstAutoReject is a user-defined Boolean that dictates whether SQUID 2.50 will be permitted to reject 'aberrant' calibration-constant values (according to its own conventions, and subject to its own rules) in order to calculate the 'best' weighted mean calibration-constant value. In practice, it is useful to allow SQUID 2.50 to identify and exclude its own rejections (i.e. user usually sets NoUPbConstAutoReject = FALSE), just to get an idea of which values it considers aberrant. Once the calculation is complete, the user has the authority to 'unreject' SQUID's outliers, and that option is frequently exercised. After all, the reference material is supposed to be isotopically homogeneous (as usually demonstrated by ID-TIMS data) and if the SHRIMP-derived population of calibration constants displays significant dispersion, then this is likely to be an instrumental effect which must be faithfully propagated to analyses of the unknowns. In some cases, rejections might be allowed to stand if there is independent evidence of (hopefully transient!) sub-optimal analytical conditions, such as instability in the primary beam or secondary beam, or poor focus of the beam on the sample.
pbCanDriftCorr is another user-defined Boolean that (if set to TRUE) instructs SQUID 2.50 to evaluate the spot-by-spot measured values of the calibration constants in terms of when they were acquired (during the course of the analytical session, typically described as "secular drift"). In practice, this can be quite a useful thing to be able to do, but at this stage of SQUID 3.0 development, it should be considered out-of-scope, and a working value of pbCanDriftCorr = FALSE should be assumed.
rw1 = plaFirstDatRw[-pbStd] --first row of StandardData
rwn = plaLastDatRw[-pbStd] --last row of StandardData
If NoUPbConstAutoReject = TRUE AND pbCanDriftCorr = FALSE
NoReject = TRUE
Else
NoReject = FALSE
End If
As before, catch non-meaningful cases where NumDauPar has been set to 2, without the proper specification of a direct calibration-expression for the 'second' DauPar, via Task.saEqns[-2]:
If NumDauPar = 2 AND Switch.DirectAltPD = TRUE AND Task.saEqns[-2] = "" --i.e. "null"
NumDauPar = 1
End If
Now a long 'For' loop, which runs once for each calibration-constant mean that needs to be calculated:
For DauParNum = 1 to NumDauPar
Acol = piaSacol[DauParNum] --select appropriate column of calib. constant values
Aecol = piaSaEcol[DauParNum] --appropriate column of calib. constant uncertainties
s = fsS(DauParNum) --string representation of DauParNum, used for labelling purposes
If (pbU = TRUE AND NumDauPar = 1) OR (pbTh = TRUE AND NumDauPar = 2)
--we are dealing specifically with a 206Pb/238U calib. constant:
Lambda = pscLm8 --i.e. 238U decay constant in units of "Ma^-1"
Ele = "U"
Else --we are dealing specifically with a 208Pb/232Th calib. constant:
Lambda = pscLm2 --i.e. 232Th decay constant in units of "Ma^-1"
Ele = "Th"
End If
Next, Name and define some Ranges: SQUID 2.50 uses the Ludwig custom functions "AddName" and "frSr" for this purpose.
"AddName" has 6 arguments (first 4 mandatory, last 2 optional) corresponding to (1) the string Name being defined, (2) a Boolean defining the sheet the Name applies to (TRUE = StandardData, FALSE = SampleData), (3) the address of the top row, (4) the address of the left column, (5) the address of the bottom row, and (6) the address of the right column, for a contiguous cell-range.
"frSr" has only the final 4 arguments of AddName, and can accept Names as input. The For loop continues:
AddName "Arr_" & s, TRUE, rw1, Acol, rwn, Acol
--Name is thus "Arr_1" or "Arr_2", for single Acol alone
AddName "Aer_" & s, TRUE, rw1, Aecol, rwn, Aecol
--Name is thus "Aer_1" or "Aer_2", for single Aecol alone
AddName "Adat_" & s, TRUE, rw1, Acol, rwn, Aecol
--Name is thus "Adat_1" or "Adat_2", for both columns
Set Arange = Range("Arr_" & s)
Set AerRange = Range("Aer_" & s)
Set AdatRange = Range ("Adat_" & s)
Npts = Count( AdatRange ) / 2 --seems indirect!
ErN = Count( AerRange )
ReDim ErrVals[1 to ErN] --apparently not zero-addressed
Now clean the uncertainty values in the Range, and the ErrVals calculation-input. Then calculate the median (and median absolute deviation) of the ErrVals, as inputs for (potential) filtering of calibration-constant value-uncertainty pairs from the eventual calculation of the weighted mean:
For i = 1 to ErN
If (fbIsNum( AdatRange[i, 2] ) = TRUE) AND (IsEmpty( AdatRange[i, 2] ) = FALSE)
--function fbIsNum is in the Ludwig Library
ErrVals[i] = AdatRange[i, 2]
Else
ErrVals[i] = 0
AdatRange[i, 2] = ""
End If
Next i
MedianEr = Median( AerRange ) --why clean col2 of AdatRange if not to use it here?
GetMAD ErrVals, ErN, MedianEr, 0, 0 --GetMAD is a Ludwig Library function (I hope)
--It is not obvious what (if anything) is achieved by this invocation of
--subroutine GetMAD, given that function fdNmad is invoked in the next line:
Nmadd = fdNmad( ErrVals ) --function fdNmad defined separately
LargeErRegN = 0
If the user has specified that they would like SQUID to attempt 'auto-rejection' of suspect calibration-constant measurements (i.e. NoReject = FALSE; unfortunately, very many of Ludwig's logic tests are double-negatives!), the next step is to assess the uncertainties of the calibration-constants, in a search for anomalous-looking values. Two criteria are used for ErrVal-based rejection:
- ErrVals residuals (relative to MedianEr) that exceed Nmadd by more than an order of magnitude
- ErrVals of 0
The unfinished 'For DauParNum...' clause continues:
If NoReject = FALSE
ReDim LargeErRej[1 to 99] --won't have more than 99 rejections!
For i = 1 to ErN
If ( ABS( ErrVals[i] - MedianEr > 10 * Nmadd ) OR ( ErrVals[i] = 0 )
--The following 'For' applies Strikethrough font to both the offending
--ErrVals[i] value AND the associated calibration-constant *VALUE*
--in the column immediately to its left. This is significant because
--when SQUID invokes functions like WtdAv in a spreadsheet environment,
--rows with Strikethrough are *EXCLUDED* from the calculation.
For j = 0 to 1
AerRange[i, j].Font.Strikethrough = TRUE
Next j
LargeErRegN = 1 + LargeErRegN
LargeErRej[LargeErRegN] = i
End If
Next i
If LargeErRegN > 0
ReDim Preserve LargeErRej[1 to LargeErRegN]
End If
End If --NoReject = FALSE
SQUID 2.50 then nominates a place for the weighted mean result to be places (the two-column output array will appear directly beneath the relevant calibration-constant value-uncertainty pair of columns, with a 3-row gap in between:
h = 3 + rwn
AddName "WtdMeanA" & s, TRUE, h, Acol --names a single cell
AddName "WtdMeanAPerr" & s, TRUE, 1 + h, Acol
--names the cell in the row directly BENEATH WtdMeanA
Now perform the actual weighted mean calculation (which might involve a detailed review of the operation of subroutine WtdAv). Remember that for the initial implementation of SQUID 3.0, we assign pbCanDriftCorr = FALSE. The unfinished 'For DauParNum...' clause continues:
If pbCanDriftCorr = TRUE
[out-of-scope stuff]
Else
wW = WtdAv Adatrange, TRUE, TRUE, 1, (NOT NoReject), TRUE, 1
The output wW of WtdAv is an array with 2 columns and up to 7 rows. The left-hand column contains mostly numeric data (values, errors, MSWD, probability of fit, etc.), and the right-hand column contains text labels for the data in the left-hand column. Remember that analyses for which Font.Strikethrough was applied (above) are excluded from the calculation in advance, and so are not assessed as part of WtdAv's "CanReject" provisions - "CanReject" data-points are separate and additional, and the two rejection-sources are amalgamated into a single list of 'rejected points' in the code below.
Before that, however, Ludwig out elements of the array wW and gives them variable-names. Unfortunately, this simple-looking process is complicated by the fact the the output of WtdAv is conditional: the third row of the 7 x 2 array does not appear if it is not applicable (which is an awful 'feature'!), so this complicates addressing of the ensuing values in what might be a 6 x 2 array. The unfinished Else continues:
WtdMean = wW[1, 1]
WtdMeanErr = wW[2, 1]
Assigning the remaining values depends on whether the data are so dispersed (MSWD >> 1) that the subroutine WtdAv calculated a 'constant external error' (i.e. a constant additional uncertainty which, when added in quadrature to EACH of the data-point uncertainties, would yield a weighted mean value with an MSWD of ~1). The unfinished Else continues:
If wW[3, 2] = "MSWD" --then data were NOT dispersed, and NO constant
--external error was calculated:
ExtPtSigma = 0
MSWD = wW[3, 1]
Prob = wW[5, 1]
Else --data WERE dispersed, a constant external error WAS calculated,
--and the addresses of subsequent array-elements is different:
ExtPtSigma = wW[3, 1]
MSWD = wW[4, 1]
Prob = wW[6, 1]
End If
In SQUID 2.50, when WtdAv is invoked with Boolean input CanReject = TRUE (as is the case here), the final row of the output array is a space-delimited list of the index-numbers of analyses rejected (remembering that, for the invocation being documented here, analyses identified for rejection by the "LargeErRej"-related criteria documented above were never assigned an index-number by WtdAv: those analyses were ignored before WtdAv even started applying its "CanReject" criteria).
The SQUID 2.50 code continues:
ParseLine wW[7, 1], Rejected(), Nrej, " " --subroutine not worth documenting
Subroutine ParseLine is basically a generalised text-editing subroutine, which transforms a delimited string into an array of the values (numeric, string, or mixed) separated by the delimiting character (specified above as a single space " "). In terms of Java implementation, all that is required here is:
- generation of a vector array (Rejected) whose elements are the analysis index-numbers specified in wW[7, 1]
- definition of integer Nrej = length(Rejected)
In SQUID 2.50, the upshot of all this is that when after the subroutine WtdAv has been invoked by subroutine WtdMeanAcalc, there exists up to two sets of rejected analyses: (1) those potentially generated by "LargeErRej" provisions above, and (2) those potentially generated by CanReject = TRUE within the WtdAv routine. With respect to the overall set of calibration constant value-uncertainty pairs, we need the union/superset of these two sets. The unfinished Else continues:
If LargeErRegN > 0
For i = 1 To Nrej
For j = 1 To LargeErRegN
--If LargeErRej[j] < Rejected[i] --Ludwig original, which is incorrect.
--The following line (Bodorkos 2018-02-25) is the correct replacement:
If LargeErRej[j] < Rejected[i] OR LargeErRej[j] = Rejected[i]
Rejected[i] = 1 + Rejected[i]
End If
Next j
Next i
End If --LargeErRegN > 0
NtotRej = Nrej + LargeErRegN
If NtotRej > 0
ReDim Preserve Rejected[1 To NtotRej]
For i = (Nrej + 1) To NtotRej
Rejected[i] = LargeErRej[i - Nrej]
Next i
BubbleSort Rejected --subroutine Bubblesort is, according to Ludwig,
--a "quick and dirty string-sorter" which can also be used for numbers.
--I have not documented it; does not seem worthwhile. All that is
--required here is that the elements of the vector Rejected are
--rearranged into ascending order.
Nrej = NtotRej
t1 = Rejected[1] --Ludwig now constructs the *COMMA*-delimited string
--of index-numbers of rejected analyses, for eventual placement at the
--base of the "WtdMeanA" summary on StandardData sheet.
For i = 2 To Nrej
t1 = t1 & ", " & Rejected[i] --separate index nos with "comma-space"
Next i
wW[7, 1] = t1
Else --i.e. NtotRej = 0
t1 = "none"
End If --NtotRej > 0
--Finally, Ludwig assembles the various elements in their 'final' form,
--as shown beneath the spot-rows on StandardData sheet:
Cells[h, Acol] = WtdMean
Cells[h + 1, Acol] = WtdMeanErr
--Assignment of Cells[h + 2, Acol] is separate, outside this Else statement
Cells[h + 3, Acol] = MSWD
Cells[h + 4, Acol] = Prob
Cells[h + 5, Acol] = wW[7, 1]
--Finally, derive an element-dependent label-string:
t2 = "Wtd Mean of Std Pb/" & Ele & " calibr."
End If --(pbCanDriftCorr = TRUE)
Set ExtPerr = Cells[h + 2, Acol] --identifies the cell to be known as ExtPerr
ExtPerr.Name = "ExtPerr" & s --names it ExtPerr1 or ExtPerr2, based on DauParNum
ExtPerr = ExtPtSigma --assigns the double-precision value calculated by WtdAv
If pbCanDriftCorr = TRUE
[out-of-scope stuff]
End If
If MSWD >= 100
MSWD = Drnd(MSWD, 3)
--Drnd is a function (Isoplot3 - Pub) that takes a double-precision number
--(the first argument), and rounds it to a set number of significant figures
--(the second argument). I don't understand the need for this If statement.
End If
If pbCanDriftCorr = TRUE
[out-of-scope stuff]
ElseIf Cells(h + 5, Acol).Value > 0
--i.e. if there ARE some rejected spots,then show the rejected calib.-
--constant values with strikethrough, and colour them yellow:
For i = 1 To Nrej
j = Adatrange[Rejected[i], 1].Row --selects the row
Fonts j, Acol, , 3 + Acol, StrikeThrough=True
--applies StrikeThrough to four columns in that row, the leftmost of
--which is the calibration-constant value column
IntClr vbYellow, j, Acol, , 3 + Acol
--applies interior colour yellow to the same range of cells
Next i
End If
The next thing SQUID 2.50 does is define the cell-ranges which will contain the "Age(Ma)" data, the "Age(Ma)±1sigma" and the "Age(Ma)±2sigma". In all cases, these three columns are stored immediately to the right of their source calibration-constant value-uncertainty pairs, as defined below (recalling that rw1 and rwn are the first and last data-rows respectively on sheet StandardData):
AddName "Aadat" & s, TRUE, rw1, 1 + Aecol, rwn, 1 + Aecol
AddName "Aaerdat1_" & s, TRUE, rw1, 2 + Aecol, rwn, 2 + Aecol
AddName "Aaerdat2_" & s, TRUE, rw1, 3 + Aecol, rwn, 3 + Aecol
Next, assemble formulae for age-calculation, which naturally depend on the identity of the daughter-parent (i.e. 206Pb/238U or 208Pb/232Th), but also the identity of the index isotope (204-corr, 207-corr, etc.) under our "calculate everything" philosophy. So below, I have used 'X' as shorthand identifier for all 'permissible' index isotopes in a given permutation (recall the Cases outlined in the Introduction of Procedural Framework Part 3).
If (pbU = TRUE AND NumDauPar = 1) OR (pbTh = TRUE AND NumDauPar = 2)
--we are dealing specifically with a 206Pb/238U calib. constant:
t1 = "=LN( 1 + ["X-corr206Pb/238Ucalibr.const"] / ("WtdMeanA" & s) * StdUPbRatio ) / Lambda"
--Now fill column ["X-corr206Pb/238U Age(Ma)"]:
For i = rw1 To rwn
PlaceFormulae t1, i, piaSAgeCol[DauParNum]
Next i
t4 = "EXP( Lambda * ["X-corr206Pb/238U Age(Ma)"] )"
t5 = ["X-corr206Pb/238Ucalibr.const %err"] / 100 * ( t4 - 1 ) / Lambda / t4
--Now fill columns ["X-corr206Pb/238U Age(Ma)±1sigma"] and
--["X-corr206Pb/238U Age(Ma)±2sigma"]
For i = rw1 To rwn
PlaceFormulae t5, i, piaSAgeECol[DauParNum] --±1sigma
PlaceFormulae "2 * t5", i, 1 + piaSAgeECol[DauParNum] --±2sigma
Next i
Else
--we are dealing specifically with a 208Pb/232Th calib. constant:
t1 = "=LN( 1 + ["X-corr208Pb/232Thcalibr.const"] / ("WtdMeanA" & s) * StdThPbRatio ) / Lambda"
--Now fill column ["X-corr208Pb/232Th Age(Ma)"]:
For i = rw1 To rwn
PlaceFormulae t1, i, piaSAgeCol[DauParNum]
Next i
t4 = "EXP( Lambda * ["X-corr208Pb/232Th Age(Ma)"] )"
t5 = ["X-corr208Pb/232Thcalibr.const %err"] / 100 * ( t4 - 1 ) / Lambda / t4
--Now fill columns ["X-corr208Pb/232Th Age(Ma)±1sigma"] and
--["X-corr208Pb/232Th Age(Ma)±2sigma"]
For i = rw1 To rwn
PlaceFormulae t5, i, piaSAgeECol[DauParNum] --±1sigma
PlaceFormulae "2 * t5", i, 1 + piaSAgeECol[DauParNum] --±2sigma
Next i
End If
The SQUID 2.50 code now embarks on a massive round of plotting and formatting, most of which seems irrelevant to the arithmetic. Buried inside it, however, is a call to the complex subroutine "ExtractGroup", which I will document separately and in detail, as it is also invoked when assessing populations of sample/unknown analyses. Its task is to find and extract the largest "statistically coherent" subgroup of analyses (i.e. with a weighted mean exceeding a specified probability-of-fit threshold), and this functionality is very useful to analysts trying to extract some sort of "rock age" from a high-n dataset complicated by lots of inherited zircon, Pb loss, etc.
If pbCanDriftCorr = TRUE
[out-of-scope stuff]
Else --i.e. pbCanDriftCorr = FALSE
ExtractGroup TRUE, 0, UPbConst, FALSE, TRUE, 0, [NULL], DauParNum, [NULL]
--subroutine documented separately
End If
Finally, I also found the code that manually over-rides a calculated external error that is "too small" (in the subjective opinion of the analyst) with a user-specified static minimum value. I have included some of the addressing code, because like everything else so far, this manipulation is theoretically applicable to every set of calibration-constant calculations:
If pbCanDriftCorr = FALSE
With Range("WtdMeanA" & s)
r = Range.Row
c = Range.Column
End With
End If
Set extBox = frSr(r, c)
With extBox
extBox.Name = "ExtPerrA" & s --note that extBox.Name is the identity
--of the external uncertainty that is actually propagated to
--measurements of the unknowns
extBox.Formula = "=MAX( "User-specified minimum external 1sigma %err", ExtPerr )
--user-specified value is set in Preferences: GA uses 0.75% based on
--long-term observations
End With
Next DauParNum
End Sub