-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathVizualAgeDRS.ipf
executable file
·2605 lines (2421 loc) · 171 KB
/
VizualAgeDRS.ipf
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
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#pragma rtGlobals=1 // Use modern global access method. - Leave this line as is, as 1st line!
#pragma ModuleName= Iolite_ActiveDRS //Leave this line as is, as 2nd line!
StrConstant DRS_Version_No= "2015.06" //Leave this line as is, as 3rd line!
//****End of Header Lines - do not disturb anything above this line!****
//****The global strings (SVar) and variables (NVar) below must always be present. Do not alter their names, alter only text to the right of the "=" on each line.**** (It is important that this line is left unaltered)
GlobalString IndexChannel ="U238"
GlobalString ReferenceStandard ="Z_91500"
GlobalString DefaultIntensityUnits ="CPS"
//**** Below are some optional global strings and variables with pre-determined behaviour. If you wish to include these simply remove the two "//" at the beginning of the line. Similarly, if you wish to omit them, simply comment them using "//"
GlobalString BeamSecondsMethod = "Rate of Change"
GlobalVariable BeamSecondsSensitivity =1
GlobalString CurveFitType ="Exponential"
GlobalVariable MaskThreshold =1000
GlobalVariable MaskEdgeDiscardSeconds =1
//**** Any global strings or variables you wish to use in addition to those above can be placed here. You may name these how you wish, and have as many or as few as you like**** (It is important that this line is left unaltered)
GlobalVariable Sample238_235Ratio =137.88
GlobalVariable DefaultStartMask =0
GlobalVariable DefaultEndMask =0
GlobalVariable MaxBeamDuration =60
GlobalVariable FitOutlierTolerance =1.5
GlobalString Ignore_235U ="Yes"
//**** If you'd like to set up some preferred settings for the report window you can set these here too
GlobalString Report_DefaultChannel ="Final206_238"
GlobalString Report_AverageMethod ="weighted (2 S.E.)"
GlobalString Report_UncertaintyMethod ="2 S.E. (absolute)"
GlobalString Report_OutlierMethod ="None"
//**** End of optional global strings and variables**** (It is important that this line is left unaltered)
//certain optional globals are built in, and have pre-determined lists. these are currently: "StandardisationMethod", "OutputUnits"
//Note that the above values will be the initial values every time the DRS is opened for the first time, but can then be overwritten for that experiment after that point via the button "Edit DRS Variables". This means that the settings for the above will effectively be stored within the experiment, and not this DRS (this is a good thing)
//DO NOT EDIT THE ABOVE IF YOU WISH TO EDIT THESE VALUES WITHIN A PARTICULAR EXPERIMENT. THESE ARE THE STARTING VALUES ONLY. THEY WILL ONLY BE LOOKED AT ONCE WHEN THE DRS IS OPENED FOR THE FIRST TIME (FOR A GIVEN EXPERIMENT).
//**** Initialisation routine for this DRS. Will be called each time this DRS is selected in the "Select DRS" popup menu (i.e. usually only once).
Function InitialiseActiveDRS() //If init func is required, this line must be exactly as written. If init function is not required it may be deleted completely and a default message will print instead at initialisation.
SVAR nameofthisDRS=$ioliteDFpath("Output","S_currentDRS") //get name of this DRS (which should have been already stored by now)
//###########################################################
// JAP
Print "DRS initialised: VizualAge U(Th)Pb zircon laser ablation module for ICP-MS, \"" + nameofthisDRS + "\", Version " + DRS_Version_No + "\r"
// Initialize VisualAge
If (strlen(FunctionInfo("VAInit")) == 0)
printabort("VisualAge isn't loaded!")
Else
Execute "VAInit()"
EndIf
// !JAP
//###########################################################
End //**end of initialisation routine
//###########################################################
// JAP
// Function used to make a 204Hg wave if all we have is a 204Pb wave
Function MakeHgWave()
Wave Pb204 = $ioliteDFPath("input", "Pb204")
Wave Pb204_time = $ioliteDFPath("input", "Pb204_time")
Variable NoOfPoints = numpnts(Pb204)
Wave Hg204 = $makeioliteWave("input", "Hg204", n=NoOfPoints)//$IoliteDFPath("input", "Hg204")
Wave Hg204_time = $makeioliteWave("input", "Hg204_time", n = NoOfPoints)
Redimension/D Hg204
Redimension/D Hg204_time
Hg204_time = Pb204_time
Hg204 = Pb204
SVAR ListOfInputChannels=$ioliteDFpath("input","GlobalListOfInputChannels") //Get reference to "GlobalListOfInputChannels", in the Input folder, and is a list of the form "ChannelName1;ChannelName2;..."
ListOfInputChannels = RemoveFromList("Hg204", ListOfInputChannels)
ListOfInputChannels += "Hg204;"
End
Function CalculateDose()
Wave Uppm = $ioliteDFpath("CurrentDRS", "Approx_U_PPM")
Wave Thppm = $ioliteDFpath("CurrentDRS", "Approx_Th_PPM")
Wave FinalAge206_238 = $iolitedfpath("CurrentDRS", "FinalAge206_238")
Wave FinalAge207_206 = $iolitedfpath("CurrentDRS", "FinalAge207_206")
Variable Npts = numpnts(Uppm)
Wave Dose = $makeiolitewave("CurrentDRS", "Dose", n = Npts)
// Loop through each time slice:
Variable N238, N235, N232
Variable l235 = 9.8485e-10
Variable l238 = 1.55125e-10
Variable l232 = 0.49475e-10
Variable k = 137.88
Variable ct = 0
Variable Navo = 6.0221413E+23
Variable i
For( i = 0; i < Npts; i = i + 1 )
N238 = 0.001*Uppm[i]*1e-6*Navo/238
N235 = 0.001*(Uppm[i]/k)*1e-6*Navo/235
N232 = 0.001*Thppm[i]*1e-6*Navo/232
ct = FinalAge206_238[i]
If (ct > 2000)
ct = FinalAge207_206[i]
EndIf
Dose[i] = 8*N238*(exp(l238*ct*1e6)-1) + 7*N235*(exp(l235*ct*1e6)-1) + 6*N232*(exp(l232*ct*1e6)-1)
EndFor
SVAR ListOfOutputChannels=$ioliteDFpath("Output","ListOfOutputChannels")
ListOfOutputChannels = RemoveFromList("Dose", ListOfOutputChannels)
ListOfOutputChannels += "Dose;"
End
// !JAP
//###########################################################
//****Start of actual Data Reduction Scheme. This is run every time raw data is added or the user presses the "crunch data" button. Try to keep it to no more than a few seconds run-time!
Function RunActiveDRS() //The DRS function name must be exactly as written here. Enter the function body code below.
ProgressDialog() //Start progress indicator
//the next 5 lines reference all of the global strings and variables in the header of this file for use in the main code of the DRS that follows.
string currentdatafolder = GetDataFolder(1)
setdatafolder $ioliteDFpath("DRSGlobals","")
SVar IndexChannel, ReferenceStandard, DefaultIntensityUnits, UseOutlierRejection, BeamSecondsMethod, CurveFitType, Ignore_235U
NVar MaskThreshold, MaskEdgeDiscardSeconds, BeamSecondsSensitivity, MaxBeamDuration, DefaultStartMask, DefaultEndMask, FitOutlierTolerance, Sample238_235Ratio
setdatafolder $currentdatafolder
//convert the long names of CurveFitType in the user interface into short labels
string ShortCurveFitType
string UserInterfaceList = "Exponential plus optional linear;Linear;Exponential;Double exponential;Smoothed cubic spline;Running median"
string ShortLabelsList = "LinExp;Lin;Exp;DblExp;Spline;RunMed"
ShortCurveFitType = StringFromList(WhichListItem(CurveFitType, UserInterfaceList, ";", 0, 0), ShortLabelsList, ";") //this line extracts the short label corresponding to the user interface label in the above string.
if(cmpstr(ShortCurveFitType, "") == 0) //if for some reason the above substitution didn't work, then need to throw an error, as that will have to be fixed
printabort("Sorry, the DRS failed to recognise the down-hole fractionation model you chose")
endif
//Do we have a baseline_1 spline for the index channel, as require this to proceed further?
DRSabortIfNotWave(ioliteDFpath("Splines", IndexChannel+"_Baseline_1")) //Abort if [index]_Baseline_1 is not in the Splines folder, otherwise proceed with DRS code below..
SetProgress(5, "Starting baseline subtraction...")
//Next, create a reference to the Global list of Output channel names, which must contain the names of all outputs produced by this routine, and to the inputs
SVAR ListOfOutputChannels=$ioliteDFpath("Output","ListOfOutputChannels") //"ListOfOutputChannels" is already in the Output folder, and will be empty ("") prior to this function being called.
SVAR ListOfIntermediateChannels=$ioliteDFpath("Output","ListOfIntermediateChannels")
SVAR ListOfInputChannels=$ioliteDFpath("input","GlobalListOfInputChannels") //Get reference to "GlobalListOfInputChannels", in the Input folder, and is a list of the form "ChannelName1;ChannelName2;..."
//Now create the global time wave for intermediate and output waves, based on the index isotope time wave ***This MUST be called "index_time" as some or all export routines require it, and main window will look for it
wave Index_Time = $MakeIndexTimeWave() //create the index time wave using the external function - it tries to use the index channel, and failing that, uses total beam
variable NoOfPoints=numpnts(Index_Time) //Make a variable to store the total number of time slices for the output waves
//THIS DRS IS A SPECIAL CASE, and has been built to allow a 'partial' data crunch, beginning after the downhole correction of ratios
NVar OptionalPartialCrunch = $ioliteDFpath("CurrentDRS","OptionalPartialCrunch")
if(NVar_Exists(OptionalPartialCrunch)!=1) //if the OptionalPartialCrunch NVar doesn't exist yet then make it here. this will only happen once, the first time the DRS is crunched
variable/g $ioliteDFpath("CurrentDRS","OptionalPartialCrunch") = 0 //so make the global variable and set it to 0
NVar OptionalPartialCrunch = $ioliteDFpath("CurrentDRS","OptionalPartialCrunch") //and reference it
endif
//The below Svar is used throughout the DRS, so place it outside the below if command
String/g $ioliteDFpath("CurrentDRS","Measured_UPb_Inputs")
SVar Measured_UPb_Inputs = $ioliteDFpath("CurrentDRS","Measured_UPb_Inputs")
if(OptionalPartialCrunch!=1) //if this is a normal crunch data then do all of this stuff, otherwise skip to the 'else' after DownHoleCurveFit()
wave IndexOut = $InterpOntoIndexTimeAndBLSub(IndexChannel) //Make an output wave for Index isotope (as baseline-subtracted intensity)
//baseline subtract all input channels. will sieve out the U Pb ones specifically afterwards
variable CurrentChannelNo
CurrentChannelNo = 0
variable NoOfChannels
NoOfChannels = itemsinlist(ListOfInputChannels) //Create local variables to hold the current input channel number and the total number of input channels
String NameOfCurrentChannel
String CurrentElement //Create a local string to contain the name of the current channel, and its corresponding element
Do //Start to loop through the available channels
NameOfCurrentChannel=StringFromList(CurrentChannelNo,ListOfInputChannels) //Get the name of the nth channel from the input list
//Can no longer use the below test, as some inputs from multicollectors are too complex and will not be recognised as elements
//CurrentElement=GetElementFromIsotope(NameOfCurrentChannel) //get name of the element
if(cmpstr(NameOfCurrentChannel, IndexChannel)!=0) //if this element is not "null" (i.e. is an element), and it is not the index isotope, then..
wave ThisChannelBLsub = $InterpOntoIndexTimeAndBLSub(NameOfCurrentChannel) //use this external function to interpolate the input onto index_time then subtract it's baseline
ListOfIntermediateChannels+=NameOfCurrentChannel+"_" + DefaultIntensityUnits +";" //Add the name of this new output channel to the list of outputs
endif //Have now created a (baseline-subtracted channel) output wave for the current input channel, unless it was TotalBeam or index
SetProgress(5+((CurrentChannelNo+1)/NoOfChannels)*10,"Processing baselines") //Update progress for each channel
CurrentChannelNo+=1 //So move the counter on to the next channel..
While(CurrentChannelNo<NoOfChannels) //..and continue to loop until the last channel has been processed.
ListOfIntermediateChannels+=IndexChannel+"_"+DefaultIntensityUnits+";" //Add the name of this new output channel to the list of outputs
//Now all baseline subtracted waves have been created.
//###########################################################
// JAP
// Do 204Pb = 204Total - F*202Hg
// where F is the ratio of 204Hg/202Hg determined from the baseline
Print "Checking if Hg correction should be applied..."
If (FindListItem("Hg204", ListOfInputChannels) != -1 && FindListItem("Pb204", ListOfInputChannels) != -1)
Print "Doing Hg correction..."
Wave Hg204 = $IoliteDFPath("input", "Hg204")
Wave Hg202 = $IoliteDFPath("input", "Hg202")
Wave Pb204 = $ioliteDFPath("input", "Pb204")
Wave Hg204_Spline = $InterpSplineOntoIndexTime("Hg204", "Baseline_1")
Wave Hg202_Spline = $InterpSplineOntoIndexTime("Hg202", "Baseline_1")
Wave HgRatio = $MakeIoliteWave("CurrentDRS", "HgRatio", n=NoOfPoints)
HgRatio = Hg204_Spline/Hg202_Spline
// Use this line to use the determined HgRatio spline:
Pb204 = Hg204 - HgRatio*Hg202
// Use this line to use the expected Hg ratio:
//Pb204 = Hg204 - 0.22987*Hg202
EndIf
// !JAP
//###########################################################
//make a mask for ratios, don't put it on baseline subtracted intermediates, as the full range is useful on these.
Wave MaskLowCPSBeam=$DRS_CreateMaskWave(IndexOut,MaskThreshold,MaskEdgeDiscardSeconds,"MaskLowCPSBeam","StaticAbsolute") //This mask currently removes all datapoints below 1000 CPS on U238, with a sideways effect of 1 second.
//The below function is called to detect which inputs are present - the format of the inputs can vary depending on the machine used to acquire the data.
//the function returns a list of the inputs in ascending mass order, with 204 at the very end if present. The list can then be used below to reference the waves in this function
//As part of the below function that detects a variety of input channel names (they vary depending on the mass spec used), make a global string to use as a reference of which Hg, Pb, Th, U isotopes have been measured (using a key=value; system)
Measured_UPb_Inputs = "200=no;202=no;204=no;206=no;207=no;208=no;232=no;235=no;238=no;"
//the "no" values in this string will be replaced by the name of the input channel for each isotope that was measured
VAGenerateUPbInputsList(ListOfIntermediateChannels)
//Now have a key=value string storing either "no" if a channel wasn't measured, or the name of the channel if it was, e.g. "200=no;202=no;204=no;206=Pb206;207=Pb207;208=Pb208;232=Th232;235=no;238=U238;"
//can now use this string to reference the relevant baseline-subtracted waves
//At its most basic level this DRS will expect at least Pb 206 and U238. The below lines check if these two are present, and report a failure if they're not
if(cmpstr(StringByKey("206", Measured_UPb_Inputs, "=", ";", 0), "no") ==0 || cmpstr(StringByKey("238", Measured_UPb_Inputs, "=", ";", 0), "no") ==0)
printabort("It appears that 206Pb or U238 were not measured. The DRS requires that as a minimum these two isotopes were measured.")
endif
//In addition to the key=value string used above, want to make flags for which ratios can be calculated - these can then be used throughout the rest of the DRS (note that 204 uses a separate flag)
variable/G $ioliteDFpath("CurrentDRS","Calculate_206_238")
variable/G $ioliteDFpath("CurrentDRS","Calculate_207_235")
variable/G $ioliteDFpath("CurrentDRS","Calculate_208_232")
variable/G $ioliteDFpath("CurrentDRS","Calculate_207_206")
variable/G $ioliteDFpath("CurrentDRS","Calculate_206_208")
NVar Calculate_206_238 = $ioliteDFpath("CurrentDRS","Calculate_206_238")
NVar Calculate_207_235 = $ioliteDFpath("CurrentDRS","Calculate_207_235")
NVar Calculate_208_232 = $ioliteDFpath("CurrentDRS","Calculate_208_232")
NVar Calculate_207_206 = $ioliteDFpath("CurrentDRS","Calculate_207_206")
NVar Calculate_206_208 = $ioliteDFpath("CurrentDRS","Calculate_206_208")
//Now set each one, depending on whether the required waves are present (for both UPb ratios check using 238, with the assumption that it will always be available and ok to use, even if 235 wasn't measured)
Calculate_206_238 = (cmpstr(StringByKey("206", Measured_UPb_Inputs, "=", ";", 0), "no") !=0 && cmpstr(StringByKey("238", Measured_UPb_Inputs, "=", ";", 0), "no") !=0)? 1 : 0
Calculate_207_235 = (cmpstr(StringByKey("207", Measured_UPb_Inputs, "=", ";", 0), "no") !=0 && cmpstr(StringByKey("238", Measured_UPb_Inputs, "=", ";", 0), "no") !=0)? 1 : 0
Calculate_208_232 = (cmpstr(StringByKey("208", Measured_UPb_Inputs, "=", ";", 0), "no") !=0 && cmpstr(StringByKey("232", Measured_UPb_Inputs, "=", ";", 0), "no") !=0)? 1 : 0
Calculate_207_206 = (cmpstr(StringByKey("207", Measured_UPb_Inputs, "=", ";", 0), "no") !=0 && cmpstr(StringByKey("206", Measured_UPb_Inputs, "=", ";", 0), "no") !=0)? 1 : 0
Calculate_206_208 = (cmpstr(StringByKey("206", Measured_UPb_Inputs, "=", ";", 0), "no") !=0 && cmpstr(StringByKey("208", Measured_UPb_Inputs, "=", ";", 0), "no") !=0)? 1 : 0
//Now start referencing the isotopes used in ratio calculation
string ThisChannelName //(re-use this string for each of the channels below)
//Hg200
ThisChannelName = StringByKey("200", Measured_UPb_Inputs, "=", ";", 0)
if(cmpstr(ThisChannelName, "no") != 0)
Wave Hg200_Beam = $ioliteDFpath("CurrentDRS", ThisChannelName)
endif
//Hg202
ThisChannelName = StringByKey("202", Measured_UPb_Inputs, "=", ";", 0)
if(cmpstr(ThisChannelName, "no") != 0)
Wave Hg202_Beam = $ioliteDFpath("CurrentDRS", ThisChannelName)
endif
//Pb204
ThisChannelName = StringByKey("204", Measured_UPb_Inputs, "=", ";", 0)
if(cmpstr(ThisChannelName, "no") != 0)
Wave Pb204_Beam = $ioliteDFpath("CurrentDRS", ThisChannelName)
endif
//Pb206
ThisChannelName = StringByKey("206", Measured_UPb_Inputs, "=", ";", 0)
if(cmpstr(ThisChannelName, "no") != 0)
Wave Pb206_Beam = $ioliteDFpath("CurrentDRS", ThisChannelName)
endif
//Pb207
ThisChannelName = StringByKey("207", Measured_UPb_Inputs, "=", ";", 0)
if(cmpstr(ThisChannelName, "no") != 0)
Wave Pb207_Beam = $ioliteDFpath("CurrentDRS", ThisChannelName)
endif
//Pb208
ThisChannelName = StringByKey("208", Measured_UPb_Inputs, "=", ";", 0)
if(cmpstr(ThisChannelName, "no") != 0)
Wave Pb208_Beam = $ioliteDFpath("CurrentDRS", ThisChannelName)
endif
//Th232
ThisChannelName = StringByKey("232", Measured_UPb_Inputs, "=", ";", 0)
if(cmpstr(ThisChannelName, "no") != 0)
Wave Th232_Beam = $ioliteDFpath("CurrentDRS", ThisChannelName)
endif
//U235
ThisChannelName = StringByKey("235", Measured_UPb_Inputs, "=", ";", 0)
if(cmpstr(ThisChannelName, "no") != 0)
Wave U235_Beam = $ioliteDFpath("CurrentDRS", ThisChannelName)
endif
//U238
ThisChannelName = StringByKey("238", Measured_UPb_Inputs, "=", ";", 0)
if(cmpstr(ThisChannelName, "no") != 0)
Wave U238_Beam = $ioliteDFpath("CurrentDRS", ThisChannelName)
endif
//have now referenced all relevant channels that have been measured
//Now as a last check, just confirm that 206 and 238 do not have zero point waves
if((numpnts(Pb206_Beam)==0)||(numpnts(U238_Beam)==0))
abort "One of the Pb, Th, or U channels are empty or missing, things are going to end badly..."
endif
//now check if a 204 beam has been measured and set a flag appropriately so that it can be used elsewhere
variable Was204Measured
ThisChannelName = StringByKey("204", Measured_UPb_Inputs, "=", ";", 0)
if(cmpstr(ThisChannelName, "no") != 0) //if 204 was measured
//set the flag to 1
Was204Measured = 1
//and reference the relevant wave
wave Pb204_Beam = $ioliteDFpath("CurrentDRS", ThisChannelName)
else //otherwise set the flag to 0
Was204Measured = 0
endif
//now, if 204 has been measured, the flag will be set to 1 and the wave has been referenced, otherwise it will be set to 0
SetProgress(20,"Calculating raw ratios...") //Update progress for each channel
if(Calculate_206_238 == 1)
Wave Raw_206_238=$MakeioliteWave("CurrentDRS","Raw_206_238",n=NoOfPoints)
Raw_206_238 = Pb206_Beam/U238_Beam * MaskLowCPSBeam
Wave Raw_Age_206_238=$MakeioliteWave("CurrentDRS","Raw_Age_206_238",n=NoOfPoints)
Raw_Age_206_238 = Ln(Raw_206_238 + 1) / 0.000155125
ListOfIntermediateChannels+="Raw_206_238;Raw_Age_206_238;"
endif
if(Calculate_207_235 == 1)
Wave Raw_207_235=$MakeioliteWave("CurrentDRS","Raw_207_235",n=NoOfPoints)
if(waveexists(U235_Beam) == 1 && cmpstr(Ignore_235U, "No") == 0)
Raw_207_235 = Pb207_Beam/U235_Beam * MaskLowCPSBeam
else
Raw_207_235 = Pb207_Beam/U238_Beam * Sample238_235Ratio * MaskLowCPSBeam
endif
Wave Raw_Age_207_235=$MakeioliteWave("CurrentDRS","Raw_Age_207_235",n=NoOfPoints)
Raw_Age_207_235 = Ln((Raw_207_235) + 1) / 0.00098485
ListOfIntermediateChannels+="Raw_207_235;Raw_Age_207_235;"
endif
if(Calculate_208_232 == 1)
Wave Raw_208_232=$MakeioliteWave("CurrentDRS","Raw_208_232",n=NoOfPoints)
Raw_208_232 = Pb208_Beam/Th232_Beam * MaskLowCPSBeam
Wave Raw_Age_208_232=$MakeioliteWave("CurrentDRS","Raw_Age_208_232",n=NoOfPoints)
Raw_Age_208_232 = Ln(Raw_208_232 + 1) / 0.000049475
ListOfIntermediateChannels+="Raw_208_232;Raw_Age_208_232;"
endif
if(Calculate_207_206 == 1)
//Call the function that will generate a lookup table to be used in calculating 207/206 ages
VAGenerate207206LookupTable()
Wave Raw_207_206=$MakeioliteWave("CurrentDRS","Raw_207_206",n=NoOfPoints)
Raw_207_206 = Pb207_Beam/Pb206_Beam * MaskLowCPSBeam
Wave Raw_Age_207_206=$MakeioliteWave("CurrentDRS","Raw_Age_207_206",n=NoOfPoints)
wave LookupTable_76 = $ioliteDFpath("CurrentDRS","LookupTable_76")
wave LookupTable_age = $ioliteDFpath("CurrentDRS","LookupTable_age")
Raw_Age_207_206 = interp(Raw_207_206, LookupTable_76, LookupTable_age)
ListOfIntermediateChannels+="Raw_207_206;Raw_Age_207_206;"
endif
if(Calculate_206_208 == 1)
Wave Raw_206_208=$MakeioliteWave("CurrentDRS","Raw_206_208",n=NoOfPoints)
Raw_206_208 = Pb206_Beam/Pb208_Beam * MaskLowCPSBeam
ListOfIntermediateChannels+="Raw_206_208;"
endif
//Now deal with 204 ratios if available
if(Was204Measured == 1)
Wave Raw_206_204=$MakeioliteWave("CurrentDRS","Raw_206_204",n=NoOfPoints)
Raw_206_204 = Pb206_Beam/Pb204_Beam * MaskLowCPSBeam
ListOfIntermediateChannels+="Raw_206_204;"
if(waveexists(Pb207_Beam) == 1)
Wave Raw_207_204=$MakeioliteWave("CurrentDRS","Raw_207_204",n=NoOfPoints)
Raw_207_204 = Pb207_Beam/Pb204_Beam * MaskLowCPSBeam
ListOfIntermediateChannels+="Raw_207_204;"
endif
if(waveexists(Pb208_Beam) == 1)
Wave Raw_208_204=$MakeioliteWave("CurrentDRS","Raw_208_204",n=NoOfPoints)
Raw_208_204 = Pb208_Beam/Pb204_Beam * MaskLowCPSBeam
ListOfIntermediateChannels+="Raw_208_204;"
endif
endif
//now want to add in a channel for U/Th ratio if both were measured (already know U238 was measured, so only need to check 232)
if(waveexists(Th232_Beam) == 1)
Wave Raw_U_Th_Ratio=$MakeioliteWave("CurrentDRS","Raw_U_Th_Ratio",n=NoOfPoints)
Raw_U_Th_Ratio = (U238_Beam/1.0000)/(Th232_Beam/1.0000) * MaskLowCPSBeam //currently using simple isotopic ratio here, can convert to elemental using 0.99275 (238U) and XXXXX (232Th)
ListOfIntermediateChannels+="Raw_U_Th_Ratio;"
endif
//Now make the BeamSeconds wave (used as a proxy for hole depth during downhole fractionation correction)
wave BeamSeconds=$DRS_MakeBeamSecondsWave(IndexOut,BeamSecondsSensitivity, BeamSecondsMethod) //This is determined by an external function which can be fine-tuned using the single sensitivity parameter. Let me know if it fails!
ListOfIntermediateChannels+="Beam_Seconds;"
//up to this point no standard is required. Need to choose at least one standard integration at this point. (e.g. Z_91500)
//Next, are we ready to proceed to producing the remaining outputs?
DRSAbortIfNotSpline(StringFromList(0,ListOfIntermediateChannels), ReferenceStandard)
// //Have now checked that at least one Z_91500 has been selected, can proceed with the following, which is for down-hole fractionation correction
SetProgress(30,"Starting down-hole curve fit...") //Update progress for each channel
//the following lines are for the down-hole correction of ratios
// JAP
NVAR HoldDHC = root:Packages:VisualAge:Options:HoldDHC
if(Calculate_208_232 == 1 && HoldDHC == 0)
VADownHoleCurveFit("Raw_208_232", OptionalWindowNumber = 0) //the optional window number can be set, if it is it allows the function to stagger the windows so that they don't all overlap completely. if it's missing it defaults to 1
endif
if(Calculate_207_235 == 1 && HoldDHC == 0)
VADownHoleCurveFit("Raw_207_235", OptionalWindowNumber = 1) //the optional window number can be set, if it is it allows the function to stagger the windows so that they don't all overlap completely. if it's missing it defaults to 1
endif
if(Calculate_206_238 == 1 && HoldDHC == 0)
VADownHoleCurveFit("Raw_206_238", OptionalWindowNumber = 2) //the optional window number can be set, if it is it allows the function to stagger the windows so that they don't all overlap completely. if it's missing it defaults to 1
endif
// !JAP
//Note that the reverse order of the ratios here just means that the topmost graph is the most commonly used (i.e. 6/38 ratio)
//NOTE: Although it may be confusing and is not a particularly nice solution, it is necessary to add waves to the list of intermediates and outputs here so that they won't be duplicated unnecessarily during a partial data crunch
//This is an unnecessarily long list, but is done this way in order to produce the desired order for the output channels
if(Calculate_207_235 == 1)
ListOfIntermediateChannels+="DC207_235;"
ListOfOutputChannels+="Final207_235;"
endif
if(Calculate_206_238 == 1)
ListOfIntermediateChannels+="DC206_238;"
ListOfOutputChannels+="Final206_238;"
endif
if(Calculate_207_206 == 1)
ListOfIntermediateChannels+="DC207_206;"
ListOfOutputChannels+="Final207_206;"
endif
if(Calculate_208_232 == 1)
ListOfIntermediateChannels+="DC208_232;"
ListOfOutputChannels+="Final208_232;"
endif
if(Calculate_206_208 == 1)
ListOfIntermediateChannels+="DC206_208;"
ListOfOutputChannels+="Final206_208;"
endif
if(Calculate_207_235 == 1)
ListOfIntermediateChannels+="DCAge207_235;"
ListOfOutputChannels+="FinalAge207_235;"
endif
if(Calculate_206_238 == 1)
ListOfIntermediateChannels+="DCAge206_238;"
ListOfOutputChannels+="FinalAge206_238;"
endif
if(Calculate_208_232 == 1)
ListOfIntermediateChannels+="DCAge208_232;"
ListOfOutputChannels+="FinalAge208_232;"
endif
if(Calculate_207_206 == 1)
ListOfIntermediateChannels+="DCAge207_206;"
ListOfOutputChannels+="FinalAge207_206;"
endif
if(Was204Measured == 1) //if 204 was measured
ListOfIntermediateChannels+="DC206_204;"
ListOfOutputChannels+="Final206_204;"
if(waveexists(Pb207_Beam) == 1)
ListOfIntermediateChannels+="DC207_204;"
ListOfOutputChannels+="Final207_204;"
endif
if(waveexists(Pb208_Beam) == 1)
ListOfIntermediateChannels+="DC208_204;"
ListOfOutputChannels+="Final208_204;"
endif
endif
//now want to add in channels for U, Th, Pb abundances
if(cmpstr(StringByKey("238", Measured_UPb_Inputs, "=", ";", 0), "no") != 0)
ListOfOutputChannels+="Approx_U_PPM;"
endif
if(cmpstr(StringByKey("232", Measured_UPb_Inputs, "=", ";", 0), "no") != 0)
ListOfOutputChannels+="Approx_Th_PPM;"
endif
if(cmpstr(StringByKey("208", Measured_UPb_Inputs, "=", ";", 0), "no") != 0)
ListOfOutputChannels+="Approx_Pb_PPM;"
endif
if(cmpstr(StringByKey("238", Measured_UPb_Inputs, "=", ";", 0), "no") != 0 && cmpstr(StringByKey("232", Measured_UPb_Inputs, "=", ";", 0), "no") != 0)
ListOfOutputChannels+="FInal_U_Th_Ratio;"
endif
//THIS IS A BIG ELSE. The following occurs if a 'partial' crunch data has been chosen. Any waves used below need to be referenced here
else
//Allow all waves to be referenced, even if some don't exist
Wave Raw_206_238=$ioliteDFpath("CurrentDRS","Raw_206_238")
Wave Raw_207_235=$ioliteDFpath("CurrentDRS","Raw_207_235")
Wave Raw_208_232=$ioliteDFpath("CurrentDRS","Raw_208_232")
Wave Raw_207_206=$ioliteDFpath("CurrentDRS","Raw_207_206")
Wave Raw_206_208=$ioliteDFpath("CurrentDRS","Raw_206_208")
wave BeamSeconds=$ioliteDFpath("CurrentDRS", "Beam_Seconds")
endif //everything after here will be executed during both the 'partial' and normal data crunches
OptionalPartialCrunch = 0 //Important: the first thing is to set the optional crunch back to the default of a full data crunch.
if(Calculate_206_238 == 1)
Wave DC206_238=$MakeioliteWave("CurrentDRS","DC206_238",n=NoOfPoints)
Wave DCAge206_238=$MakeioliteWave("CurrentDRS","DCAge206_238",n=NoOfPoints)
endif
if(Calculate_207_235 == 1)
Wave DC207_235=$MakeioliteWave("CurrentDRS","DC207_235",n=NoOfPoints)
Wave DCAge207_235=$MakeioliteWave("CurrentDRS","DCAge207_235",n=NoOfPoints)
endif
if(Calculate_208_232 == 1)
Wave DC208_232=$MakeioliteWave("CurrentDRS","DC208_232",n=NoOfPoints)
Wave DCAge208_232=$MakeioliteWave("CurrentDRS","DCAge208_232",n=NoOfPoints)
endif
string CoefficientWaveName, ratio, SmoothWaveName, SplineWaveName, AverageBeamSecsName //various strings required by the different fit types below
strswitch(ShortCurveFitType)
case "LinExp":
if(Calculate_206_238 == 1)
ratio = "Raw_206_238"
CoefficientWaveName = "LECoeff_" + ratio
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
NVar Variable_b = $ioliteDFpath("CurrentDRS","LEVarB_"+ratio) //variable b is the linear component of the equation
DC206_238 = Raw_206_238 / (1 + (Variable_b/Coefficients[0])*BeamSeconds + (Coefficients[1]/Coefficients[0])*Exp(-Coefficients[2]*BeamSeconds))
endif
if(Calculate_207_235 == 1)
ratio = "Raw_207_235"
CoefficientWaveName = "LECoeff_" + ratio
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
NVar Variable_b = $ioliteDFpath("CurrentDRS","LEVarB_"+ratio) //variable b is the linear component of the equation
DC207_235 = Raw_207_235 / (1 + (Variable_b/Coefficients[0])*BeamSeconds + (Coefficients[1]/Coefficients[0])*Exp(-Coefficients[2]*BeamSeconds))
endif
if(Calculate_208_232 == 1)
ratio = "Raw_208_232"
CoefficientWaveName = "LECoeff_" + ratio
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
NVar Variable_b = $ioliteDFpath("CurrentDRS","LEVarB_"+ratio) //variable b is the linear component of the equation
DC208_232 = Raw_208_232 / (1 + (Variable_b/Coefficients[0])*BeamSeconds + (Coefficients[1]/Coefficients[0])*Exp(-Coefficients[2]*BeamSeconds))
endif
break
case "Exp":
if(Calculate_206_238 == 1)
CoefficientWaveName = "ExpCoeff_" + "Raw_206_238"
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
DC206_238 = Raw_206_238 / (1+(Coefficients[1]/Coefficients[0])*Exp(-Coefficients[2]*BeamSeconds))//this equation is trying to change the magnitude of the original std wave to equal 1 at beamseconds = infinity, otherwise it will alter the ratios depending on the value obtained for the standard (could exploit this by factoring in a simultaneous drift correction?)
endif
if(Calculate_207_235 == 1)
CoefficientWaveName = "ExpCoeff_" + "Raw_207_235"
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
DC207_235 = Raw_207_235 / (1+(Coefficients[1]/Coefficients[0])*Exp(-Coefficients[2]*BeamSeconds))
endif
if(Calculate_208_232 == 1)
CoefficientWaveName = "ExpCoeff_" + "Raw_208_232"
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
DC208_232 = Raw_208_232 / (1+(Coefficients[1]/Coefficients[0])*Exp(-Coefficients[2]*BeamSeconds))
endif
break
case "DblExp":
if(Calculate_206_238 == 1)
CoefficientWaveName = "DblExpCoeff_" + "Raw_206_238"
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
DC206_238 = Raw_206_238 / (1 + (Coefficients[1]/Coefficients[0])*Exp(-Coefficients[2]*BeamSeconds) + (Coefficients[3]/Coefficients[0])*Exp(-Coefficients[4]*BeamSeconds))//this equation (y = K0+K1*exp(-K2*x)+K3*exp(-K4*x)) is trying to change the magnitude of the original std wave to equal 1 at beamseconds = infinity
endif
if(Calculate_207_235 == 1)
CoefficientWaveName = "DblExpCoeff_" + "Raw_207_235"
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
DC207_235 = Raw_207_235 / (1 + (Coefficients[1]/Coefficients[0])*Exp(-Coefficients[2]*BeamSeconds) + (Coefficients[3]/Coefficients[0])*Exp(-Coefficients[4]*BeamSeconds))//this equation (y = K0+K1*exp(-K2*x)+K3*exp(-K4*x)) is trying to change the magnitude of the original std wave to equal 1 at beamseconds = infinity
endif
if(Calculate_208_232 == 1)
CoefficientWaveName = "DblExpCoeff_" + "Raw_208_232"
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
DC208_232 = Raw_208_232 / (1 + (Coefficients[1]/Coefficients[0])*Exp(-Coefficients[2]*BeamSeconds) + (Coefficients[3]/Coefficients[0])*Exp(-Coefficients[4]*BeamSeconds))//this equation (y = K0+K1*exp(-K2*x)+K3*exp(-K4*x)) is trying to change the magnitude of the original std wave to equal 1 at beamseconds = infinity
endif
break
case "Lin":
if(Calculate_206_238 == 1)
CoefficientWaveName = "LinCoeff_" + "Raw_206_238"
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
DC206_238 = Raw_206_238 / (1+(Coefficients[1]/Coefficients[0])*BeamSeconds)//this equation is trying to change the magnitude of the original std wave to equal 1 at beamseconds = 0, otherwise it will alter the ratios depending on the value obtained for the standard (could exploit this by factoring in a simultaneous drift correction?)
endif
if(Calculate_207_235 == 1)
CoefficientWaveName = "LinCoeff_" + "Raw_207_235"
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
DC207_235 = Raw_207_235 / (1+(Coefficients[1]/Coefficients[0])*BeamSeconds)
endif
if(Calculate_208_232 == 1)
CoefficientWaveName = "LinCoeff_" + "Raw_208_232"
wave Coefficients = $ioliteDFpath("CurrentDRS",CoefficientWaveName)
DC208_232 = Raw_208_232 / (1+(Coefficients[1]/Coefficients[0])*BeamSeconds)
endif
break
case "RunMed":
if(Calculate_206_238 == 1)
ratio = "Raw_206_238"
SmoothWaveName = "SmthFitCurve_"+Ratio
wave SmoothedWave = $ioliteDFpath("CurrentDRS",SmoothWaveName)
wave AverageBeamSeconds = $ioliteDFpath("CurrentDRS","AverageBeamSecs_"+ratio)
DC206_238 = Raw_206_238 / ForInterp(Beamseconds, AverageBeamSeconds, SmoothedWave) * SmoothedWave[0]
endif
if(Calculate_207_235 == 1)
ratio = "Raw_207_235"
SmoothWaveName = "SmthFitCurve_"+Ratio
wave SmoothedWave = $ioliteDFpath("CurrentDRS",SmoothWaveName)
wave AverageBeamSeconds = $ioliteDFpath("CurrentDRS","AverageBeamSecs_"+ratio)
DC207_235 = Raw_207_235 / ForInterp(Beamseconds, AverageBeamSeconds, SmoothedWave) * SmoothedWave[0]
endif
if(Calculate_208_232 == 1)
ratio = "Raw_208_232"
SmoothWaveName = "SmthFitCurve_"+Ratio
wave SmoothedWave = $ioliteDFpath("CurrentDRS",SmoothWaveName)
wave AverageBeamSeconds = $ioliteDFpath("CurrentDRS","AverageBeamSecs_"+ratio)
DC208_232 = Raw_208_232 / ForInterp(Beamseconds, AverageBeamSeconds, SmoothedWave) * SmoothedWave[0]
endif
break
case "Spline":
if(Calculate_206_238 == 1)
ratio = "Raw_206_238"
SplineWaveName = "SplineCurve_"+Ratio
wave SplineWave = $ioliteDFpath("CurrentDRS",SplineWaveName)
wave SplineBeamSeconds = $ioliteDFpath("CurrentDRS","SplineBeamSecs_"+ratio) //note that the spline method creates its own beamseconds wave that extends beyond averagebeamseconds
DC206_238 = Raw_206_238 / ForInterp(Beamseconds, SplineBeamSeconds, SplineWave) * SplineWave[0]
endif
if(Calculate_207_235 == 1)
ratio = "Raw_207_235"
SplineWaveName = "SplineCurve_"+Ratio
wave SplineWave = $ioliteDFpath("CurrentDRS",SplineWaveName)
wave SplineBeamSeconds = $ioliteDFpath("CurrentDRS","SplineBeamSecs_"+ratio)
DC207_235 = Raw_207_235 / ForInterp(Beamseconds, SplineBeamSeconds, SplineWave) * SplineWave[0]
endif
if(Calculate_208_232 == 1)
ratio = "Raw_208_232"
SplineWaveName = "SplineCurve_"+Ratio
wave SplineWave = $ioliteDFpath("CurrentDRS",SplineWaveName)
wave SplineBeamSeconds = $ioliteDFpath("CurrentDRS","SplineBeamSecs_"+ratio)
DC208_232 = Raw_208_232 / ForInterp(Beamseconds, SplineBeamSeconds, SplineWave) * SplineWave[0]
endif
break
endswitch
if(Calculate_206_238 == 1)
DCAge206_238 = Ln(DC206_238 + 1) / 0.000155125
Wave Final206_238=$MakeioliteWave("CurrentDRS","Final206_238",n=NoOfPoints)
Wave FinalAge206_238=$MakeioliteWave("CurrentDRS","FinalAge206_238",n=NoOfPoints)
endif
if(Calculate_207_235 == 1)
DCAge207_235 = Ln((DC207_235) + 1) / 0.00098485
Wave Final207_235=$MakeioliteWave("CurrentDRS","Final207_235",n=NoOfPoints)
Wave FinalAge207_235=$MakeioliteWave("CurrentDRS","FinalAge207_235",n=NoOfPoints)
endif
if(Calculate_208_232 == 1)
DCAge208_232 = Ln(DC208_232 + 1) / 0.000049475
Wave Final208_232=$MakeioliteWave("CurrentDRS","Final208_232",n=NoOfPoints)
Wave FinalAge208_232=$MakeioliteWave("CurrentDRS","FinalAge208_232",n=NoOfPoints)
endif
//at the moment I don't think Pb-Pb ratios need any treatment, so they are left as they were...
if(Calculate_207_206 == 1)
Wave DC207_206=$MakeioliteWave("CurrentDRS","DC207_206",n=NoOfPoints)
Wave DCAge207_206=$MakeioliteWave("CurrentDRS","DCAge207_206",n=NoOfPoints)
DC207_206 = Raw_207_206
wave LookupTable_76 = $ioliteDFpath("CurrentDRS","LookupTable_76")
wave LookupTable_age = $ioliteDFpath("CurrentDRS","LookupTable_age")
DCAge207_206 = interp(DC207_206, LookupTable_76, LookupTable_age)
Wave Final207_206=$MakeioliteWave("CurrentDRS","Final207_206",n=NoOfPoints)
Wave FinalAge207_206=$MakeioliteWave("CurrentDRS","FinalAge207_206",n=NoOfPoints)
endif
if(Calculate_206_208 == 1)
Wave DC206_208=$MakeioliteWave("CurrentDRS","DC206_208",n=NoOfPoints)
DC206_208 = Raw_206_208
Wave Final206_208=$MakeioliteWave("CurrentDRS","Final206_208",n=NoOfPoints)
endif
if(Was204Measured == 1) //if 204 was measured
Wave DC206_204=$MakeioliteWave("CurrentDRS","DC206_204",n=NoOfPoints)
DC206_204 = Raw_206_204
Wave Final206_204=$MakeioliteWave("CurrentDRS","Final206_204",n=NoOfPoints)
if(waveexists(Pb207_Beam) == 1)
Wave DC207_204=$MakeioliteWave("CurrentDRS","DC207_204",n=NoOfPoints)
DC207_204 = Raw_207_204
Wave Final207_204=$MakeioliteWave("CurrentDRS","Final207_204",n=NoOfPoints)
endif
if(waveexists(Pb208_Beam) == 1)
Wave DC208_204=$MakeioliteWave("CurrentDRS","DC208_204",n=NoOfPoints)
DC208_204 = Raw_208_204
Wave Final208_204=$MakeioliteWave("CurrentDRS","Final208_204",n=NoOfPoints)
endif
endif
SetProgress(40,"Calculating final ratios...") //Update progress for each channel
//so, have done down-hole correction, now need to do drift correction (note that there is often a substantial offset in raw and down hole corr. values from true values)
//(relevant waves were already made above in the if statements)
//and make some waves for the approximated elemental concentrations
//now want to add in a channel for U, Th, Pb abundances
if(cmpstr(StringByKey("238", Measured_UPb_Inputs, "=", ";", 0), "no") != 0)
Wave Approx_U_PPM=$MakeioliteWave("CurrentDRS","Approx_U_PPM",n=NoOfPoints)
endif
if(cmpstr(StringByKey("232", Measured_UPb_Inputs, "=", ";", 0), "no") != 0)
Wave Approx_Th_PPM=$MakeioliteWave("CurrentDRS","Approx_Th_PPM",n=NoOfPoints)
endif
if(cmpstr(StringByKey("208", Measured_UPb_Inputs, "=", ";", 0), "no") != 0)
Wave Approx_Pb_PPM=$MakeioliteWave("CurrentDRS","Approx_Pb_PPM",n=NoOfPoints)
endif
if(cmpstr(StringByKey("238", Measured_UPb_Inputs, "=", ";", 0), "no") != 0 && cmpstr(StringByKey("232", Measured_UPb_Inputs, "=", ";", 0), "no") != 0)
Wave FInal_U_Th_Ratio=$MakeioliteWave("CurrentDRS","FInal_U_Th_Ratio",n=NoOfPoints)
endif
//now call an external function to do the actual drift correction. The reason for using an external function is that it's also called during export for error propagation.
VADriftCorrectRatios() //this function has optional range parameters, but they can be left blank and the function will operate on the entire wave
//replacing the original error propagation with the below generic version (no blossoming of errors with low N integrations)
//propagate each ratio separately, using the relevant down-hole corrected ratio
//############################################################
// JAP
// Get rid of old waves:
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalDiscPercent")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAnd207_235")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAnd206_238")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAnd208_232")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAgeAnd207_206")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalFracCommonPb")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalFracLostPb")
KillWaves/Z $ioliteDFpath("CurrentDRS", "AndersenDeltaAge")
KillWaves/Z $ioliteDFpath("CurrentDRS", "AndersenSolution")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAgeAnd206_238")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAgeAnd207_235")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAgeAnd208_232")
KillWaves/Z $ioliteDFpath("CurrentDRS", "Final238_206")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAnd238_206")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAnd207_206")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalPbC206_238")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalPbC207_235")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalPbC208_232")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalPbC238_206")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalPbC207_206")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAgePbC206_238")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAgePbC207_235")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAgePbC208_232")
KillWaves/Z $ioliteDFpath("CurrentDRS", "FinalAgePbC207_206")
// Remove everything from the output list (they'll be added back later):
ListOfOutputChannels = RemoveFromList("Final238_206;FinalAnd238_206;FinalAnd207_206;", ListOfOutputChannels)
ListOfOutputChannels = RemoveFromList("FinalAnd207_235;FinalAnd206_238;", ListOfOutputChannels)
ListOfOutputChannels = RemoveFromList("FinalAnd208_232;FinalAgeAnd207_206;FinalAgeAnd206_238;", ListOfOutputChannels)
ListOfOutputChannels = RemoveFromList("FinalAgeAnd207_235;FinalAgeAnd208_232;FinalFracCommonPb;", ListOfOutputChannels)
ListOfOutputChannels = RemoveFromList("FinalFracLostPb;AndersenSolution;FinalDiscPercent;", ListOfOutputChannels)
ListOfOutputChannels = RemoveFromList("FinalPbC206_238;FinalPbC207_235;FinalPbC208_232;FinalPbC207_206;FinalPbC238_206;", ListOfOutputChannels)
ListOfOutputChannels = RemoveFromList("FinalAgePbC206_238;FinalAgePbC207_235;FinalAgePbC208_232;FinalAgePbC207_206;", ListOfOutputChannels)
NVAR doPbPb = root:Packages:VisualAge:Options:PbPbOption_Calculate
NVAR doAndersen = root:Packages:VisualAge:Options:AndersenOption_Calculate
NVAR do204 = root:Packages:VisualAge:Options:PbCOption_Calculate
NVAR Have204 = root:Packages:VisualAge:Was204Measured
Have204 = Was204Measured
// Make the required waves:
Wave Final238_206=$MakeioliteWave("CurrentDRS", "Final238_206", n=NoOfPoints)
Wave FinalDiscPercent=$MakeioliteWave("CurrentDRS", "FinalDiscPercent", n = NoOfPoints)
If (doPbPb)
Killwaves/Z $ioliteDFpath("CurrentDRS", "FinalAge207_206")
ListOfOutputChannels = RemoveFromList("FinalAge207_206;", ListOfOutputChannels)
Wave FinalAge207_206=$MakeioliteWave("CurrentDRS", "FinalAge207_206", n = NoOfPoints)
EndIf
If (doAndersen)
Wave FinalAnd207_235=$MakeioliteWave("CurrentDRS", "FinalAnd207_235", n = NoOfPoints)
Wave FinalAnd206_238=$MakeioliteWave("CurrentDRS", "FinalAnd206_238", n = NoOfPoints)
Wave FinalAnd208_232=$MakeioliteWave("CurrentDRS", "FinalAnd208_232", n = NoOfPoints)
Wave FinalAnd238_206=$MakeioliteWave("CurrentDRS", "FinalAnd238_206", n=NoOfPoints)
Wave FinalAnd207_206=$MakeioliteWave("CurrentDRS", "FinalAnd207_206", n=NoOfPoints)
Wave FinalAgeAnd207_206=$MakeioliteWave("CurrentDRS", "FinalAgeAnd207_206", n = NoOfPoints)
Wave FinalAgeAnd206_238=$MakeioliteWave("CurrentDRS", "FinalAgeAnd206_238", n = NoOfPoints)
Wave FinalAgeAnd207_235=$MakeioliteWave("CurrentDRS", "FinalAgeAnd207_235", n = NoOfPoints)
Wave FinalAgeAnd208_232=$MakeioliteWave("CurrentDRS", "FinalAgeAnd208_232", n = NoOfPoints)
Wave FinalFracCommonPb=$MakeioliteWave("CurrentDRS", "FinalFracCommonPb", n = NoOfPoints)
Wave FinalFracLostPb=$MakeioliteWave("CurrentDRS", "FinalFracLostPb", n = NoOfPoints)
Wave AndersenDeltaAge=$MakeioliteWave("CurrentDRS", "AndersenDeltaAge", n = NoOfPoints)
Wave AndersenSolution=$MakeioliteWave("CurrentDRS", "AndersenSolution", n = NoOfPoints)
EndIf
If (do204 && Was204Measured)
Wave FinalPbC206_238=$MakeioliteWave("CurrentDRS", "FinalPbC206_238", n = NoOfPoints)
Wave FinalPbC207_235=$MakeioliteWave("CurrentDRS", "FinalPbC207_235", n = NoOfPoints)
Wave FinalPbC208_232=$MakeioliteWave("CurrentDRS", "FinalPbC208_232", n = NoOfPoints)
Wave FinalPbC238_206=$MakeioliteWave("CurrentDRS", "FinalPbC238_206", n = NoOfPoints)
Wave FinalPbC207_206=$MakeioliteWave("CurrentDRS", "FinalPbC207_206", n = NoOfPoints)
Wave FinalAgePbC206_238=$MakeioliteWave("CurrentDRS", "FinalAgePbC206_238", n = NoOfPoints)
Wave FinalAgePbC207_235=$MakeioliteWave("CurrentDRS", "FinalAgePbC207_235", n = NoOfPoints)
Wave FinalAgePbC208_232=$MakeioliteWave("CurrentDRS", "FinalAgePbC208_232", n = NoOfPoints)
Wave FinalAgePbC207_206=$MakeioliteWave("CurrentDRS", "FinalAgePbC207_206", n = NoOfPoints)
EndIf
// Calculate the 38/06 ratio:
ListOfOutputChannels += "Final238_206;"
Final238_206 = 1/Final206_238
// Do 7/6 age calculation if desired:
If (doPbPb)
SetProgress(50, "Calculating 207Pb/206Pb ages")
ListOfOutputChannels += "FinalAge207_206;"
PbPbAges()
EndIf
// Calculate a rough measure of discordance:
ListOfOutputChannels += "FinalDiscPercent;"
SetProgress(60, "Calculating degree of discordance")
CalculateDisc()
CalculateDose()
// Do Andersen correction if desired:
If (doAndersen)
ListOfOutputChannels += "FinalAnd207_235;FinalAnd206_238;FinalAnd208_232;FinalAnd238_206;FinalAnd207_206;"
ListOfOutputChannels += "FinalAgeAnd207_206;FinalAgeAnd207_235;FinalAgeAnd206_238;FinalAgeAnd208_232;"
ListOfOutputChannels += "FinalFracCommonPb;FinalFracLostPb;AndersenSolution;"
// Then iterate Andersen routine until ages don't change:
NVAR reAndersen = root:Packages:VisualAge:RecrunchAndersen
NVAR maxAndersenItrs = root:Packages:VisualAge:Options:AndersenOption_MaxRecalc
SetProgress(70, "Doing Andersen's common-Pb correction")
Print "[VisualAge] Andersen routine will execute " + num2str(maxAndersenItrs) + " times or less."
Variable numAndersenItrs = 0
AndersenDeltaAge = -1
AndersenSolution = 0
Do
Print "[VisualAge] Starting Andersen iteration " + num2str (numAndersenItrs + 1) + "."
Andersen()
numAndersenItrs = numAndersenItrs + 1
While ( reAndersen && numAndersenItrs < maxAndersenItrs )
FinalAnd238_206 = 1/FinalAnd206_238
Variable j
For (j = 0; j < NoOfPoints; j = j + 1)
FinalAnd207_206[j] = Ratio7_6(FinalAgeAnd207_206[j]*1e6)
EndFor
EndIf
// Do 204Pb common Pb correction if desired:
If (do204 && Was204Measured)
ListOfOutputChannels += "FinalPbC207_235;FinalPbC206_238;FinalPbC208_232;FinalPbC207_206;FinalPbC238_206;"
ListOfOutputChannels += "FinalAgePbC207_206;FinalAgePbC207_235;FinalAgePbC206_238;FinalAgePbC208_232;"
SetProgress(80, "Doing 204Pb common-Pb correction")
Calculate204PbCorrection()
EndIf
// !JAP
//############################################################
SetProgress(90,"Propagating errors...") //Update progress for each channel
string ListOfOutputsToPropagate
if(Calculate_206_238 == 1)
ListOfOutputsToPropagate = "Final206_238;FinalAge206_238"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC206_238", ReferenceStandard)
endif
if(Calculate_207_235 == 1)
ListOfOutputsToPropagate = "Final207_235;FinalAge207_235"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC207_235", ReferenceStandard)
endif
if(Calculate_208_232 == 1)
ListOfOutputsToPropagate = "Final208_232;FinalAge208_232"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC208_232", ReferenceStandard)
endif
if(Calculate_207_206 == 1)
ListOfOutputsToPropagate = "Final207_206;FinalAge207_206"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC207_206", ReferenceStandard)
endif
if(Calculate_206_208 == 1)
ListOfOutputsToPropagate = "Final206_208"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC206_208", ReferenceStandard)
endif
If (doAndersen)
if (Calculate_206_238 == 1)
ListOfOutputsToPropagate = "FinalAnd206_238;FinalAgeAnd206_238"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC206_238", ReferenceStandard)
endif
if (Calculate_207_235 == 1)
ListOfOutputsToPropagate = "FinalAnd207_235;FinalAgeAnd207_235"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC207_235", ReferenceStandard)
endif
if (Calculate_208_232 == 1)
ListOfOutputsToPropagate = "FinalAnd208_232;FinalAgeAnd208_232"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC208_232", ReferenceStandard)
endif
if (Calculate_207_206 == 1)
ListOfOutputsToPropagate = "FinalAnd207_206;FinalAgeAnd207_206"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC207_206", ReferenceStandard)
endif
EndIf
if (do204 && Was204Measured)
if (Calculate_206_238 == 1)
ListOfOutputsToPropagate = "FinalPbC206_238;FinalAgePbC206_238"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC206_238", ReferenceStandard)
endif
if (Calculate_207_235 == 1)
ListOfOutputsToPropagate = "FinalPbC207_235;FinalAgePbC207_235"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC207_235", ReferenceStandard)
endif
if (Calculate_208_232 == 1)
ListOfOutputsToPropagate = "FinalPbC208_232;FinalAgePbC208_232"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC208_232", ReferenceStandard)
endif
if (Calculate_207_206 == 1)
ListOfOutputsToPropagate = "FinalPbC207_206;FinalAgePbC207_206"
Propagate_Errors("All", ListOfOutputsToPropagate, "DC207_206", ReferenceStandard)
endif
endif
SetProgress(100,"Finished DRS.") //Update progress for each channel
end //****End of DRS function. Write any required external sub-routines below this point****
//############################################################
// JAP
//------------------------------------------------------------------------
// Calculates discordance of each point
//------------------------------------------------------------------------
Function CalculateDisc()
Wave FinalDiscPercent=$ioliteDFpath("CurrentDRS", "FinalDiscPercent")
Wave Final207_235=$ioliteDFpath("CurrentDRS", "Final207_235")
Wave Final206_238=$ioliteDFpath("CurrentDRS", "Final206_238")
Wave FinalAge206_238=$ioliteDFpath("CurrentDRS", "FinalAge206_238")
Wave FinalAge207_206=$ioliteDFpath("CurrentDRS", "FinalAge207_206")
Variable Npts = numpnts(FinalDiscPercent)
Variable i
For( i = 0; i < Npts; i = i + 1 )
//FinalDiscPercent[i] = DiscPercent(Final207_235[i], Final206_238[i])
FinalDiscPercent[i] = DiscPercent2(FinalAge206_238[i], FinalAge207_206[i])
EndFor
End
//------------------------------------------------------------------------
// Calculates the 207/206 age of each point
//------------------------------------------------------------------------
Function PbPbAges()
Print "[VisualAge] Starting calculation of 207Pb/206Pb ages."
// Get time for calculation start:
Variable calcStartTime = DateTime
// Get required waves from iolite:
Wave Final207_206=$ioliteDFpath("CurrentDRS", "Final207_206")
Wave FinalAge207_235=$ioliteDFpath("CurrentDRS", "FinalAge207_235")
Wave FinalAge206_238=$ioliteDFpath("CurrentDRS", "FinalAge206_238")
Wave FinalAge207_206=$ioliteDFpath("CurrentDRS", "FinalAge207_206")
// Get number of wave points:
Variable Npts = numpnts(FinalAge207_206)
SVAR wavesForGuess = root:Packages:VisualAge:Options:PbPbOption_WavesForGuess
// Loop through each point, and call the PbPb routine for each:
Variable i
For ( i = 1; i < Npts; i = i + 1 )
// Get current 7/6 ratio:
Variable m = (Final207_206[i] + Final207_206[i-1] + Final207_206[i+1])/3
// Calculate a reasonable guess at the age:
Variable guess = 1e6*(AgeFromList(wavesForGuess, i))
// If the ratio or age seem unreasonable set age to NaN and skip
If (numtype(m) == 2 || guess <= 1 || guess > 5e9 || numtype(guess) == 2)
FinalAge207_206[i] = NaN
Continue
EndIf
// Call Newton's method PbPb function:
FinalAge207_206[i] = CalculatePbPbAge(m, guess)
EndFor
// Get time of completion:
Variable calcStopTime = DateTime
// Spit out some info:
Print "[VisualAge] ...Done. Calculation duration: " + num2str(calcStopTime - calcStartTime) + " s."
End
//------------------------------------------------------------------------
// Calculates 204Pb corrections for each point
//------------------------------------------------------------------------
Function Calculate204PbCorrection()
Variable calcStartTime = DateTime
Print "[VisualAge] Starting 204Pb correction."
Wave Final207_204 = $ioliteDFpath("CurrentDRS", "Final207_204")
Wave Final206_204 = $ioliteDFpath("CurrentDRS", "Final206_204")
Wave Final208_204 = $ioliteDFpath("CurrentDRS", "Final208_204")
Wave Final207_206 = $ioliteDFpath("CurrentDRS", "Final207_206")
Wave FinalAge206_238 = $ioliteDFpath("CurrentDRS", "FinalAge206_238")
Wave FinalAge207_206 = $ioliteDFpath("CurrentDRS", "FinalAge207_206")
Wave FinalPbC206_238 = $ioliteDFpath("CurrentDRS", "FinalPbC206_238")
Wave FinalPbC238_206 = $ioliteDFpath("CurrentDRS", "FinalPbC238_206")
Wave FinalPbC207_235 = $ioliteDFpath("CurrentDRS", "FinalPbC207_235")
Wave FinalPbC208_232 = $ioliteDFpath("CurrentDRS", "FinalPbC208_232")
Wave FinalPbC207_206 = $ioliteDFpath("CurrentDRS", "FinalPbC207_206")
Wave Final207_235 = $ioliteDFpath("CurrentDRS", "Final207_235")
Wave Final206_238 = $ioliteDFpath("CurrentDRS", "Final206_238")
Wave Final208_232 = $ioliteDFpath("CurrentDRS", "Final208_232")
Wave FinalAgePbC207_235 = $ioliteDFpath("CurrentDRS", "FinalAgePbC207_235")
Wave FinalAgePbC206_238 = $ioliteDFpath("CurrentDRS", "FinalAgePbC206_238")
Wave FinalAgePbC208_232 = $ioliteDFpath("CurrentDRS", "FinalAgePbC208_232")
Wave FinalAgePbC207_206 = $ioliteDFpath("CurrentDRS", "FinalAgePbC207_206")
SVAR WavesForGuess = root:Packages:VisualAge:Options:PbCOption_WavesForGuess
NVAR UsePbComp = root:Packages:VisualAge:Options:Option_UsePbComp
NVAR Pb64 = root:Packages:VisualAge:Options:Option_Common64
NVAR Pb74 = root:Packages:VisualAge:Options:Option_Common74
NVAR Pb84 = root:Packages:VisualAge:Options:Option_Common84