A brief demonstration of the double spike toolbox

This is a quick guide to the main features of the double spike toolbox.

Contents

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
Welcome to the double spike toolbox.    John F. Rudge 2009.
Default parameters are stored in the global variable ISODATA.
Type "global ISODATA" to access and set these parameters.

ans = 

           isonum: [54 56 57 58]
             mass: [53.9396 55.9349 56.9354 57.9333]
         standard: [0.0585 0.9175 0.0212 0.0028]
         rawspike: [4x4 double]
          isoname: {'Fe54'  'Fe56'  'Fe57'  'Fe58'}
         isolabel: {'^{54}Fe'  '^{56}Fe'  '^{57}Fe'  '^{58}Fe'}
          element: 'Fe'
         isoindex: [1x58 double]
      revisoindex: [1x58 double]
            nisos: 4
          nratios: 3
          nspikes: 4
    rawspikelabel: {'spike 1'  'spike 2'  'spike 3'  'spike 4'}
       errormodel: [1x1 struct]

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,:)
ans =

    54    56    57    58


ans =

    0.0585    0.9175    0.0212    0.0028


ans =

   53.9396   55.9349   56.9354   57.9333


ans =

    0.0010    0.0697    0.9288    0.0005


ans =

    0.0039    0.1378    0.0125    0.8458

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])
ans =

    0.0488    0.9448    0.0007    0.0055    0.0001    0.0001


ans =

    0.0210    0.0002    0.0001    0.0007    0.0001    0.9778

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])
optspike =

    0.0025    0.1059    0.4413    0.4502


optprop =

    0.4574


opterr =

    0.0036


optisoinv =

    54    56    57    58


optspikeprop =

         0         0    0.4680    0.5320


optppmperamu =

   64.8550

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')
optspike =

         0    0.7728         0    0.2272
         0         0    0.4765    0.5235
         0    0.7531    0.2469         0
    0.7996         0         0    0.2004
    0.4487    0.5513         0         0
    0.7023         0    0.2977         0


optprop =

    0.5540
    0.4486
    0.6916
    0.2148
    0.3982
    0.2306


opterr =

    0.0032
    0.0035
    0.0045
    0.0093
    0.0109
    0.0263


optisoinv =

    54    56    57    58
    54    56    57    58
    54    56    57    58
    54    56    57    58
    54    56    57    58
    54    56    57    58


optspikeprop =

         0    0.7728         0    0.2272
         0         0    0.4765    0.5235
         0    0.7531    0.2469         0
    0.7996         0         0    0.2004
    0.4487    0.5513         0         0
    0.7023         0    0.2977         0


optppmperamu =

   57.1338
   62.4727
   80.1060
  165.5689
  193.3196
  468.5264

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')
optspike =

    0.0015    0.7957    0.0033    0.1996
    0.0025    0.1059    0.4413    0.4502
    0.0008    0.7654    0.2337    0.0002
    0.7863    0.0401    0.0031    0.1704
    0.4293    0.5701    0.0005    0.0001
    0.6906    0.0317    0.2775    0.0002


optprop =

    0.5586
    0.4574
    0.6961
    0.2121
    0.4014
    0.2328


opterr =

    0.0033
    0.0036
    0.0046
    0.0095
    0.0111
    0.0267


optisoinv =

    54    56    57    58
    54    56    57    58
    54    56    57    58
    54    56    57    58
    54    56    57    58
    54    56    57    58


optspikeprop =

         0    0.7641         0    0.2359
         0         0    0.4680    0.5320
         0    0.7487    0.2513         0
    0.7986         0         0    0.2014
    0.4360    0.5640         0         0
    0.7017         0    0.2983         0


optppmperamu =

   59.0686
   64.8550
   82.5825
  169.1806
  196.6900
  475.5390

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])
optspike =

    0.6344         0    0.3656         0
    0.7290    0.2710         0         0
         0    0.7183    0.2817         0
    0.2219         0         0    0.7781
         0         0    0.1093    0.8907
         0    0.3768         0    0.6232


optprop =

    0.5124
    0.5073
    0.5274
    0.5629
    0.6150
    0.5337


opterr =

    0.0031
    0.0032
    0.0034
    0.0045
    0.0048
    0.0056


optisoinv =

   204   206   207   208
   204   206   207   208
   204   206   207   208
   204   206   207   208
   204   206   207   208
   204   206   207   208


optspikeprop =

    0.6344         0    0.3656         0
    0.7290    0.2710         0         0
         0    0.7183    0.2817         0
    0.2219         0         0    0.7781
         0         0    0.1093    0.8907
         0    0.3768         0    0.6232

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,:)
ans =

         0    0.3964         0         0         0    0.6036
         0    0.4642         0         0    0.5358         0
         0         0    0.4353         0         0    0.5647
         0         0    0.4898         0    0.5102         0
         0         0         0    0.1605         0    0.8395
         0         0         0    0.1550         0    0.8450
         0         0         0    0.1545         0    0.8455
         0         0         0         0    0.4725    0.5275
         0         0         0    0.2298    0.7702         0
         0         0         0    0.2238    0.7762         0
         0         0         0    0.2162    0.7838         0
         0    0.4827    0.5173         0         0         0
         0    0.2126         0         0    0.7874         0
         0         0    0.3766         0    0.6234         0
         0    0.2933    0.7067         0         0         0
         0    0.0741         0    0.9259         0         0
         0         0    0.3098         0    0.6902         0
         0         0         0    0.4462    0.5538         0
         0    0.1221         0         0    0.8779         0
         0         0    0.1190         0         0    0.8810
         0    0.0429         0         0         0    0.9571
         0    0.0375         0         0         0    0.9625
         0         0    0.4774         0    0.5226         0
         0         0    0.5363         0    0.4637         0
         0         0         0         0    0.3084    0.6916
         0         0    0.3493         0         0    0.6507
         0         0         0    0.1584         0    0.8416
         0         0         0    0.1605         0    0.8395
         0         0         0    0.3211    0.6789         0
         0         0         0         0    0.3514    0.6486
         0         0    0.1725    0.8275         0         0

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,:)
ans =

    0.2657    4.4859    2.6304    2.6179
    0.2658    4.4858    2.6305    2.6180
    0.2657    4.4860    2.6302    2.6180
    0.2658    4.4860    2.6304    2.6180
    0.2657    4.4860    2.6303    2.6180
    0.2657    4.4860    2.6303    2.6179
    0.2657    4.4858    2.6302    2.6178
    0.2657    4.4861    2.6303    2.6180
    0.2658    4.4859    2.6302    2.6179
    0.2657    4.4861    2.6303    2.6181

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)
monteerror =

    0.0035


predictederror =

    0.0036

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
ans = 

    measured: [1x1 struct]
       spike: [1x1 struct]
    standard: [1x1 struct]


ans = 

    intensity: 10
            a: [2.0710e-10 2.0710e-10 2.0710e-10 2.0710e-10]
            b: [2.0027e-09 2.0027e-09 2.0027e-09 2.0027e-09]
            c: [0 0 0 0]

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(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(); % return error model to defaults
error1 =

    0.0036


error2 =

    0.0025


ans =

    0.6757