-
Notifications
You must be signed in to change notification settings - Fork 16
SHRIMP Sub: sqConcAge
sbodorkos edited this page Sep 6, 2019
·
4 revisions
This subroutine calculates a Concordia Age from the double-precision values in the selected 4-column Tera-Wasserburg range, where the uncertainties on the 238/206 and 207/206 values are 1sigma PERCENTAGE. It can operate either in Automatic mode (Automatic = TRUE, in which case input-rows are auto-rejected until either probability of fit or MSWD becomes acceptable, using a procedure analogous to the previously described FindCoherentGroup) or in Manual mode (Automatic = FALSE, in which case the calculation is attempted solely by reference to the pre-existing rejected/unrejected statuses of the input-rows).
sqConcAge(Automatic)
Automatic: Boolean dictating whether the subroutine should attempt auto-rejection of input-rows in order to meet threshold values for probability of fit or MSWD. Default value is FALSE.
Values of type Boolean
Automatic, b, Bad, NoXY
Values of type Integer
Area, Ccol, Col, df, i, MaxRind, N, n0, plHdrRw, PtNum, Rw
Values of type Double
Age, cProb, cMSWD, Emult, ErrX, ErrY, MaxRsd, MinCprob, MSWD, pdMinFract, pdMinProb, Prob, RhoXY, SigmaAge, SigmaXp, SigmaYp, SumsXY, Xbar, Ybar
Values of type String
s, tmp
Arrays of type Double
inpdat, Pts, X, Xerr, Y, Yerr
Start by setting some threshold values, and doing some error-trapping:
MinCprob = 0.1 --Threshold minimum probability of XY-equivalence, and concordance,
--to be satisfied in order to define success when Automatic = TRUE
plHdrRw = flHeaderRow(0) --Set plHdrRow to be column-headers on Grouped-sample sheet
FindStr "238/206*", , Ccol, plHdrRw, 1 --Find first column containing "238/206*" in
--plHdrRw, store its index-number as Ccol. As noted in Part 7c, Squid3 should find
--ALL occurrences of 238/206 in Tera-Wasserburg blocks (uncorrected, 204Pb-corrected,
--and 208Pb-corrected) and perform sqConcAge for each occurrence.
If (Ccol = 0)
Exit Sub
End If --(Ccol = 0)
If (Automatic = TRUE)
--Some error-trapping: not sure of purpose. Maybe to cover pbExtractAgeGroups
--pbExtractAgeGroups = FALSE? Because if TRUE, pdMinProb and pdMinFract are
--specified explicitly by user-operated 'spinners' in the 'Group Me' form.
If (pdMinProb = 0)
pdMinProb = 0.15 --I think 0.05 would be better, but leave alone at present
End If --(pdMinProb = 0)
If (pdMinFract = 0)
pdMinFract = 0.7 --Looks too high, would prefer 0.3, but leave alone for now
End If --(pdMinFract = 0)
--Now count the number of input-rows beneath the header-row:
Rw = plHdrRw
Do
Rw = Rw + 1
If (Rw > 999) --Too many?
Exit Sub
End If --(Rw > 999)
Loop Until Cells(Rw, Ccol).Borders(xlEdgeBottom).LineStyle = xlContinuous
--i.e. Ludwig stops counting rows when he finds the enclosing border at the base
--of the Tera-Wasserburg column-block! Very fragile, even with Automatic = TRUE!
--Now select the whole Tera-Wasserburg block, as (1 + plHdrRw) = first data-row,
--Ccol = leftmost T-W column, Rw = gross extent of input-rows, (3 + Ccol) = the
--3 columns to the right of Ccol, containing the relevant data:
frSr(1 + plHdrRw, Ccol, Rw, 3 + Ccol).Select --this becomes "Selection"
Selection.Font.Strikethrough = FALSE --Removes any Strikethrough (i.e. rejections)
--that might have been in place before this point; essentially ensures that when
--Automatic = TRUE, calculations start from "none rejected".
Else --i.e. (Automatic = FALSE), covering manual trigger (i.e. Redo) of ConcordiaAge
--In this scenario, a manual Selection will already exist (albeit potentially
--non-contiguous rows, as permitted by Excel using Ctrl-click). So the nested For
--below simply checks the selected Area(s) for "column geography", because these
--still need to precisely encompass Ccol to (3 + Ccol).
With Selection ' Can select ANY cells in the radiogenic T-W range.
For Area = 1 To .Areas.Count --number of non-contiguous blocks
For Col = .Areas(Area).Column To
( .Areas(Area).Column + .Areas(Area).Columns.Count - 1 )
If (Col < Ccol) OR (Col > (3 + Ccol))
s = "You must select one or more rows from the 4 radiogenic
Tera-Wasserburg concordia columns"
MsgBox s, , pscSq --SQUID-branded message-box with OK button
Exit Sub
End If --(Col < Ccol) OR (Col > (3 + Ccol))
Next Col
Next Area
End With
End If --(Automatic = TRUE)
Having confirmed the selection and column-geography, the next step is to inspect all the cells for Strikethrough (i.e. evidence of rejection), and build vectors X, Xerr, Y, Yerr. X and Y are input 238/206 and 207/206 from rows in which NO cells are rejected, and Xerr and Yerr are their (converted) 1sigma ABSOLUTE uncertainties (go figure!):
Do
N = 0
With Selection
For Area = 1 To .Areas.Count
--Check each cell for non-numeric, Strikethrough, or zero in key columns:
For Rw = .Areas(Area).Row To
( .Areas(Area).Row + .Areas(Area).Rows.Count - 1 )
b = TRUE
For Col = Ccol To (3 + Ccol)
Set tObj = Cells[Rw, Col] --makes each Cell a Range for checking
If (IsNumeric(tObj) = FALSE) OR (tObj.Font.Strikethrough = TRUE) OR
(tObj.Value = 0)
b = FALSE
Exit For
End If --(IsNumeric(tObj) = FALSE) OR (tObj.Font.Strikethrough = TRUE)
--OR (tObj.Value = 0)
Next Col
If (b = TRUE) --then the row-data are valid and unrejected
N = 1 + N
ReDim Preserve X[1 To N], Y[1 To N], Xerr[1 To N], Yerr[1 To N]
X[N] = Cells[Rw, Ccol] --238U/206Pb ratio-value
Y[N] = Cells[Rw, 2 + Ccol] --207Pb/206Pb ratio-value
--Now convert the respect PERCENTAGE input-errors to ABSOLUTE:
Xerr[N] = Cells[Rw, 1 + Ccol] * X[N] / 100 --238/206 error as ABS
Yerr[N] = Cells[Rw, 3 + Ccol] * Y[N] / 100 --207/206 error as ABS
End If --(b = TRUE)
Next Rw
Next Area
End With
If (N < 1) --no need to continue
Exit Sub
End If --(N < 1)
If (n0 = 0)
n0 = N --preserves "pre-WtdXYmean" value of N
End If --(n0 = 0)
--Now build the array-input Pts to WtdXYmean:
ReDim Pts[1 To N, 1 To 5]
For PtNum = 1 To N
Pts[PtNum, 1] = X[PtNum]
Pts[PtNum, 2] = Xerr[PtNum]
Pts[PtNum, 3] = Y[PtNum]
Pts[PtNum, 4] = Yerr[PtNum]
Pts[PtNum, 5] = 0 --Ludwig ignores error-correlations in Tera-Wasserburg data
Next PtNum
--Then do it! (Function/subroutine documented in LudwigLibrary)
Isoplot3.WtdXYmean Pts, N, Xbar, Ybar, SumsXY, ErrX, ErrY, RhoXY, Bad
df = 2 * (N - 1)
If (df > 0) --i.e. N > 1
MSWD = SumsXY / df
Prob = FDIST(MSWD, df, 1e9) --FDIST is an Excel function described previously
Else --i.e. N = 1
MSWD = 0
Prob = 1
End If --(df > 0)
--Is this result good enough? If Automatic = FALSE, then YES, because it is the
--result requested. Otherwise, we need to satisfy MinCprob or MSWD. To check:
If ( (Automatic = FALSE) OR ( (Prob >= MinCprob) OR (MSWD <= 1.4) ) )
--MSWD value seems arbitrary! Not completely convinced the constraint is needed
Exit Do --move on to next phase
End If --( (Automatic = FALSE) OR ( (Prob >= MinCprob) OR (MSWD <= 1.4) ) )
If (N = 1) --should never happen; should satisfy "Exit Do" conditions above)
Exit Sub
End If --(N = 1)
If ( (N - 1) / n0 ) < pdMinFract ) OR ( Bad = TRUE ) --then the 2D weighted
--mean calculation has failed, irrespective of Concordia considerations:
Selection.Font.Strikethrough = FALSE --removes all Strikethrough/reject marks
NoXY = TRUE
Exit Do --move on to next phase
End If --( (N - 1) / n0 ) < pdMinFract ) OR ( Bad = TRUE )
--But if we have not yet exited the Do Loop, then Automatic = TRUE, but MSWD and
--Prob do not yet satisfy the criteria. So, time for stepwise rejections. This is
--done by trimming the point with the highest weighted-residual. Recall that these
--are stored in a Range named Yf, generated by subroutine WtdXYmean!
MaxRsd = 0
For PtNum = 1 To N
If ( Yf.WtdResid[PtNum] = 0 ) OR ( Yf.WtdResid[PtNum] > MaxRsd )
--then we need to find out which row it is
With Selection
Col = .Column
For Rw = .Row To ( .Row + .Rows.Count - 1 )
--Following If (which checks that the 238/206 and 207/206 ratio-values
--are numbers) looks redundant based on previous error-trapping measures:
If ( IsNumeric(Cells[Rw, Col] = TRUE ) AND
( IsNumeric(Cells[Rw, 2 + Col] = TRUE )
If ( Cells[Rw, Col] = X[PtNum] ) AND
( Cells(Rw, 2 + Col] = Y[PtNum] )
--Check it hasn't already been rejected!
If ( frSr(Rw, Col, Rw, 3 + Col).Font.Strikethrough = FALSE )
MaxRsd = Yf.WtdResid[PtNum]
MaxRind = PtNum
End If --( frSr(Rw, Col, Rw, 3 + Col).Font.Strikethrough = FALSE )
Exit For --Rw = .Row To ( .Row + .Rows.Count - 1 )
Else
MaxRsd = Yf.WtdResid[PtNum]
MaxRind = PtNum
Exit For --Rw = .Row To ( .Row + .Rows.Count - 1 )
End If --( Cells[Rw, Col] = X[PtNum] ) AND
--( Cells(Rw, 2 + Col] = Y[PtNum] )
End If --( IsNumeric(Cells[Rw, Col] = TRUE ) AND
--( IsNumeric(Cells[Rw, 2 + Col] = TRUE )
Next Rw
End With
End If --( Yf.WtdResid[PtNum] = 0 ) OR ( Yf.WtdResid[PtNum] > MaxRsd )
Next PtNum
b = FALSE --Now find the new rejected point
With Selection
Col = .Column
For Rw = .Row To ( .Row + .Rows.Count - 1 )
--Following If (which checks that the 238/206 and 207/206 ratio-values
--are numbers) looks redundant based on previous error-trapping measures:
If ( IsNumeric(Cells[Rw, Col] = TRUE ) AND
( IsNumeric(Cells[Rw, 2 + Col] = TRUE )
If ( Cells[Rw, Col] = X[MaxRind] ) AND ( Cells[Rw, 2 + Col] = Y[MaxRind] )
If ( frSr(Rw, Col, Rw, 3 + Col).Font.Strikethrough = FALSE )
--then apply Strikethrough!
frSr(Rw, Col, Rw, 3 + Col).Font.Strikethrough = TRUE
b = TRUE
Exit For --Rw = .Row To ( .Row + .Rows.Count - 1 )
End If --( frSr(Rw, Col, Rw, 3 + Col).Font.Strikethrough = FALSE )
End If --( Cells[Rw, Col] = X[MaxRind] ) AND
--( Cells[Rw, 2 + Col] = Y[MaxRind] )
End If --( IsNumeric(Cells[Rw, Col] = TRUE ) AND
--( IsNumeric(Cells[Rw, 2 + Col] = TRUE )
Next Rw
End With
If (b = FALSE)
Exit Do --move on to next phase
End If --(b = FALSE)
Loop
This should conclude the calculation of the two-dimensional weighted mean. Rw is essentially a cursor control, moving below the input-rows in order to write results (or, alternatively, non-result strings that contain failure information):
Rw = Rw + 2 --SQUID 2.50 moves 2 rows beneath the final input-row, to prepare to
--write the output of the summary calculations
If ( (Bad = TRUE) AND (Automatic = FALSE) )
--SQUID-branded message-box with OK button
MsgBox "Could not calculate the X-Y wtd mean", , pscSq
Exit Sub
End If --( (Bad = TRUE) AND (Automatic = FALSE) )
If we have gotten this far, then there is an WtdXYmean value (even if it isn't a very good one). The last loop is to report the result, and augment the SigmaX value by the internal error of the Standard if applicable (this process is exactly analogous to the 1D process documented at the end of subroutine ExtractGroup). If applicable, the modified WtdXYmean data is then reassembled for calculation of a Concordia Age (in Tera-Wasserburg mode), if warranted.
Ludwig handles a lot of the outputs as strings, where he takes the numeric results, converts the numbers to strings (via his fsS function), and wraps explanatory text strings around the ex-numeric strings to make coherent composite strings. I hope it is reasonably self-explanatory (it is to me, because I know what the real output looks like in SQUID 2.50); if it poses problems, let me know and I will provide some examples of output.
If ( (Prob < MinCprob) AND (MSWD > 1.4) AND (NoXY = FALSE) )
--then there is an XYmean value, but it violates coherence constraints
If (Automatic = TRUE)
Exit Sub
Else --for Automatic = FALSE, make a failure advisory re probability of equiv.
If (Prob < 0.00001)
tmp = "~zero)" --string
Else
tmp = fsS( DRND( Prob, 2 ) ) & ")" --DRND rounds Prob to 2 sig. figures;
--fsS converts that result to a string
End If --(Prob < 0.00001)
--SQUID-branded message-box advising inadequate Pequiv, and its value as string
MsgBox "Probability of X-Y equivalence is too low (" & tmp, , pscSq
End If --(Automatic = TRUE)
Else --we can now start assembling WtdXYmean result for ConcordiaAge calculation
If ( (Prob < 0.15) AND (NoXY = FALSE) ) --then expand the errors by Student's t
--(at 68.3% confidence level) and SQRT(MSWD):
Emult = SQRT( MSWD ) * TINV( (1 - 0.683), df ) --TINV is Excel function
ErrX = ErrX * Emult --ErrX is ABSOLUTE error output from WtdXYmean
ErrY = ErrY * Emult --ErrY is ABSOLUTE error output from WtdXYmean
End If --( (Prob < 0.15) AND (NoXY = FALSE) )
SigmaXp = 100 * ErrX / Xbar --reconverts to PERCENTAGE!
SigmaYp = 100 * ErrY / Ybar --reconverts to PERCENTAGE!
If (N > 1) --should always be TRUE, based on foregoing error-trapping
--Perform final modifications to XYmean if needed, and write results
If (NoXY = TRUE)
s = "No coherent Concordia group"
Else
--Start writing
Fonts Rw, Ccol, "X-Y wtd mean (68%-conf. errs), incl. error from Standard:"
Rw = Rw + 1 --move down one row
Cells[Rw, Ccol] = Xbar --leftmost column, write WtdXYmean output Xbar
--Step to second-from-left column, and MODIFY SigmaXp by quadratic addition!
--The below formula adds the internal error of the standard to the WtdXYmean
--SigmaXp, exactly analogous to corresponding 1D calculation at the end of
--subroutine ExtractGroup (recall that DRND rounds SigmaXp to 3 significant
--figures, and fsS converts the DRND result to string):
Cells[Rw, 1 + Ccol] = " =SQRT( " & fsS( DRND( SigmaXp, 3 ) ) &
"^2 + ( StandardData!WtdMeanAPerr1 )^2 ) "
Cells[Rw, 2 + Ccol] = Ybar --step to third-from-left column, write Ybar
Cells[Rw, 3 + Ccol] = DRND( SigmaYp, 3 ) --step to rightmost column, write
--DRND-rounded SigmaYp
If (Automatic = TRUE) --then construct explanatory string for result
If (N = n0)
s = "on all" & fsS( N ) & " points."
ElseIf (N < n0) --probably no other possibility
s = "on " & fsS( N ) & " points of" & fsS( n0 ) & " (arbitrary rejections!)."
End If --(N = n0)
End If --(Automatic = TRUE)
End If --(NoXY = TRUE)
--Step down another row and write string s
Fonts (Rw + 1), Ccol, s
If (NoXY = FALSE) --prepare (modified) WtdXYmean output as ConcordiaAges input
ReDim Pts[1 To 1, 1 To 5]
Pts[1, 1] = Xbar
--Take care to pick up the quadratically modified Xbar if relevant:
Pts[1, 2] = ( Cells[Rw, 1 + Ccol].Value ) * Xbar / 100
Pts[1, 3] = Ybar
Pts[1, 4] = ErrY
Pts[1, 5] = 0
ReDim inpdat[1 To 1, 1 To 5] --can't see where this is used
For i = 1 To 5
inpdat[1, i] = Pts[1, i]
Next i
Rw = 2 + Rw --step down another 2 rows, and prepare the string for MSWD of
--equivalence and Prob-equivalence, reporting each to 2 significant figures:
s = "Probability of equivalence = 0" & fsS( DRND( Prob, 2 ) ) &
" (MSWD = " & fsS( Drnd( MSWD, 2 ) ) & ")"
Fonts Rw, Ccol, s
Rw = 1 + Rw --step down another row, in preparation for ConcordiaAges:
Inverse = TRUE --Isoplot needs this, to be working in Tera-Wasserburg mode,
--rather than conventional (Wetherill concordia) mode...
--Now attempt to calculate the ConcordiaAge via LudwigLibrary subroutine:
Isoplot3.ConcordiaAges Pts, 1, Bad, Age, SigmaAge, cMSWD, cProb
--Prepare reporting based on outputs:
If (Bad = TRUE) OR (cProb < 0.05) --then the result is discordant; bail out
--and report offending probability of concordance to 2 significant figures:
s = "DISCORDANT (probability of concordance =" & fsS( DRND( cProb, 2 ) ) & ")"
Fonts Rw, Ccol, s --write the string
Else --write a string containing the results! FIXED is an Excel function
--specifying rounding to specified number of DECIMAL PLACES
If (N = 1) --context not clear by reading; I think this N is the SECOND
--argument of Isoplot3.ConcordiaAges. Nothing else makes much sense.
tmp = "Concordia Age = " & fsS( FIXED( Age, 1 ) ) & " ± " &
fsS( FIXED( 1.96 * SigmaAge, 1 ) ) &
" (95% confidence, incl. error from Standard)"
Else
tmp = "Concordia Age = " & fsS( FIXED( Age, 1 ) ) & " ± " &
fsS( FIXED( 1.96 * SigmaAge, 1 ) ) & " (95% confidence)"
End If --(N = 1)
Fonts Rw, Ccol, tmp --write the string tmp, as alternative to s
Rw = 1 + Rw --step down 1 more row, to write probability of concordance:
tmp = "Probability of concordance = 0" & fsS( FIXED( cProb, 3 ) )
Fonts Rw, Ccol, tmp --write the string
End If --(Bad = TRUE) OR (cProb < 0.05)
End If --(NoXY = FALSE)
End If --(N > 1), really should be obsolete as noted above
End If --( (Prob < MinCprob) AND (MSWD > 1.4) AND (NoXY = FALSE) )
I think that's it!