-
Notifications
You must be signed in to change notification settings - Fork 1
/
example.m
179 lines (158 loc) · 8.82 KB
/
example.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
%% A brief demonstration of the double spike toolbox
% This is a quick guide to the main features of the double spike toolbox.
%
%% Startup
% The function _dsstartup_ is used to initialise a number of parameters to default values. These
% are stored in a global variable called ISODATA. This is a structure with fieldnames corresponding
% to the different elements. Information includes the isotope numbers, the atomic masses, the
% standard compositions, and the default error model coefficients. Some of these values are
% loaded in from the file 'maininput.csv' which can modified by the user as necessary.
% _dsstartup_ is called automatically the first time one of the double spike toolbox functions
% is used. The information stored about Fe is shown below.
dsstartup;
global ISODATA;
ISODATA.Fe
%% Default values
% The default values can be accessed and set through the ISODATA global variable. For example,
% shown below for Fe are the isotope numbers, standard value,
% the atomic masses, and the third and fourth Oak Ridge National Labs (ORNL) spike
% compositions (corresponding to the 57Fe and 58Fe spikes). Note that
% all values are expressed as proportions of each isotope, rather than as isotopic ratios.
ISODATA.Fe.isonum
ISODATA.Fe.standard
ISODATA.Fe.mass
ISODATA.Fe.rawspike(3,:)
ISODATA.Fe.rawspike(4,:)
%% 2D Error surfaces 1
% The function _errorcurve2d_ plots a 2D error surface as a contour plot, showing how
% the error varies with both double spike to sample mixture proportions and the proportions in which
% the two single spikes are mixed to make the double spike. The example below is for a 57Fe-58Fe
% double spike, using pure spikes.
figure(1);
errorcurve2d('Fe','pure',[57 58]);
%% 2D Error surfaces 2
% Calculations can also be performed using real spikes rather than hypothetical pure spikes. Below
% is an example using the third and fourth ORNL spikes (corresponding to 57Fe and 58Fe spikes).
% This produces Figure 2 in the manuscript. Spike purity usually has only a little effect.
figure(2);
errorcurve2d('Fe','real',[3 4]);
%% 2D Error surfaces 3
% Some isotopes have more than four isotopes, and in these cases the isotopes to use in the
% inversion must be specified. Here is an example for Ca, with a 42Ca-48Ca double spike
% and 40Ca, 42Ca, 44Ca, 48Ca used in inversion. The second and sixth ORNL spikes are the 42Ca and 48Ca
% spikes. This produces Figure 5 in the manuscript.
figure(3);
ISODATA.Ca.rawspike(2,:)
ISODATA.Ca.rawspike(6,:)
errorcurve2d('Ca','real',[2 6],[40 42 44 48])
%% Error curves 1
% The function _errorcurve_ plots the error in either the fractionation factor alpha
% or a chosen ratio as a function of the double spike to sample proportion. Note again that all
% compositions are specified by proportions of each isotope rather than by ratios. Here we look
% at the error curve for a double spike which is 50% 57Fe and 50% 58Fe.
figure(4);
spike=[0 0 0.5 0.5]; % define a spike composition as 50% 57Fe and 50% 58Fe
errorcurve('Fe',spike);
%% Error curves 2
% The function _errorcurve2_ plots the error in either the fractionation factor alpha
% or a chosen ratio as a function of the proportion of the two single spikes that make the double spike.
% The proportion in which double spike and sample are mixed must be specified, as must the single spikes
% to use. We give an example here for a pure 57Fe-58Fe double spike with 50% double spike to 50% sample
figure(5);
errorcurve2('Fe','pure',0.5,[57 58])
%% Optimal spikes 1
% The function _optimalspike_ finds double spike compositions which minimise the error on alpha
% or a chosen ratio. This can be done either for pure spikes, or with the real spikes available
% from Oak Ridge National Labs.
% The following example finds the best 57Fe-58Fe spike using the real spikes available from ORNL.
% The 3rd and 4th spikes correspond to the 57Fe and 58Fe spikes. The optimal double spike turns
% out to be quite close to a 50-50 mix of the available spikes (optspikeprop). The actual double
% spike compositions are in optspike, the optimal double spike-sample mixing proportions in optprop,
% the error estimates in opterr, rescaled error estimates in optppmperamu, and the isotopes that
% were used in the inversion in optisoinv.
[optspike,optprop,opterr,optisoinv,optspikeprop,optppmperamu]=optimalspike('Fe','real',[3 4])
%% Optimal spikes 2
% If the isotopes to spike are not specified, the _optimalspike_ function checks all possible combinations.
% An example for Fe pure spikes is shown below. This produces Table 1 in the manuscript.
[optspike,optprop,opterr,optisoinv,optspikeprop,optppmperamu]=optimalspike('Fe','pure')
%% Optimal spikes 3
% An example for Fe ORNL spikes is shown below. This produces Table 2 in the manuscript.
[optspike,optprop,opterr,optisoinv,optspikeprop,optppmperamu]=optimalspike('Fe','real')
%% Optimal spikes 4
% By default, _optimalspike_ minimises the error on alpha, but for radiogenic work we often
% wish to minimise the error on a particular ratio. An example of this is Pb. Shown below is the
% result of minimising the error on 206Pb/204Pb. This produces part of Table 3 in the manuscript.
[optspike,optprop,opterr,optisoinv,optspikeprop]=optimalspike('Pb','pure',[],[],[206 204])
%% Optimal spikes 5
% For elements with more than four isotopes, such as Ca, _optimalspike_
% tries all combinations of four isotopes in the inversion.
% This produces Table 4 in the manuscript. We only display the first 31 rows of the optimal spike
% composition.
[optspike,optprop,opterr,optisoinv,optspikeprop,optppmperamu]=optimalspike('Ca','pure');
optspike(1:31,:)
%% Optimal error curves 1
% The function _errorcurveoptimalspike_ calculates the _optimalspike_, and then plots the _errorcurve_.
figure(6);
errorcurveoptimalspike('Fe','real',[3 4]);
%% Optimal error curves 2
% If the isotopes to spike are not specified, all possible combinations are shown.
% This produces Figure 3 in the manuscript.
figure(7);
errorcurveoptimalspike('Fe','real');
%% Monte Carlo fake mass spec runs
% _monterun_ performs Monte Carlo fake mass spec runs.
% The example below uses a 50-50 spike-sample mix, the pure 50-50 57Fe-58Fe spike used earlier, a
% natural fractionation of -0.2, an instrumental fractionation of 1.8, with 1000 Monte Carlo samples.
% The first 10 mixture measurements are shown.
measured=monterun('Fe',0.5,spike,-0.2,1.8,1000);
measured(1:10,:)
%% Double spike inversions
% _dsinversion_ performs the double spike inversion. Here we run
% the double spike inversion on the Monte-Carlo generated data with the chosen spike. We then
% produce a figure showing how the value of alpha varies over the run.
out=dsinversion('Fe',measured,spike);
figure(8);
plot(out.alpha);
ylabel('\alpha');
%% Error estimates
% The _errorestimate_ routine estimates the errors by linear error propagation.
% Here we compare the error got from the Monte-Carlo simulation with that
% predicted by linear error propagation. The fact that these are close is a good
% validation of the linear error propagation method.
monteerror=std(out.alpha)
predictederror=errorestimate('Fe',0.5,spike,[],[],-0.2,1.8)
%% Error model 1
% The coefficients of the error model are in contained in the global variable ISODATA.
% The coefficients can be specified for the measured, standard (unspiked run), and double
% spike compositions. See Appendix C of the manuscript for their definition.
ISODATA.Fe.errormodel
ISODATA.Fe.errormodel.measured
%% Error model 2
% The function _seterrormodel_ can be used to simply change the intensity and integration time
% for the default error model. In the example below, a doubling of the intensity decreases the error
% by roughly a factor of 1/sqrt(2).
error1=errorestimate('Fe',0.5,spike,[],[],-0.2,1.8)
seterrormodel('Fe',20,8); % set a 20 V total beam with 8 second integrations
error2=errorestimate('Fe',0.5,spike,[],[],-0.2,1.8)
error2/error1
seterrormodel('Fe'); % return error model to defaults
%% Error model 3
% The default error model of the double spike toolbox fixes the total voltage of the beams for the overall
% mixture at a given level. When sample-limited it may be more appropriate to consider an error model
% where the voltage for the sample is fixed (see John 2012, J. Anal. At. Spectrom.).
% In the toolbox this can be achieved by changing the error model type from 'fixed-total' to 'fixed-sample'.
% The code below recreates Figure 8 of Klaver and Coath 2019, Geostandards and Geoanalytical Research.
isoinv = [58 60 61 62];
spike6062 = [0.0132 0.3295 0.0014 0.6547 0.0012];
spike6162 = [0.0109 0.0081 0.4496 0.5297 0.0017];
seterrormodel('Ni')
% fix the sample intensity at 0.5 V for all runs
ISODATA.Ni.errormodel.measured.intensity = 0.5;
ISODATA.Ni.errormodel.measured.type = 'fixed-sample';
figure(9);
errorcurve('Ni', spike6062, isoinv);
hold on
errorcurve('Ni', spike6162, isoinv);
ylim([0 0.06]);
legend('^{60}Ni-^{62}Ni spike', '^{61}Ni-^{62}Ni spike')
seterrormodel('Ni'); % return error model to defaults