Home Forums Code snippets Data Reduction Scheme for LA Sr data (with CaAr and REE2+ corrections)

This topic contains 3 replies, has 2 voices, and was last updated by  Graham Hagen-Peter 1 year, 5 months ago.

  • Creator
    Topic
  • #8271

    Graham Hagen-Peter
    Participant

    Hello everyone,

    I have recently started doing Sr isotope measurements by laser ablation (mostly plagioclase) and quickly realized the challenges of this method due to numerous and complex interferences on the Sr (and Rb) masses of interest. I think that the downloadable Iolite package comes with two Sr data reduction schemes (“Sr_isotopes_CaAr” and “Sr_isotopes_REE”). I’ve included below a data reduction scheme that corrects for both CaAr and doubly charged REE interferences (in two different ways). I apologize in advance for the lengthy post, but I think it’s important to explain what the DRS is doing to eliminate any mystery. I will explain the calculations below, but please read several disclaimers first:

    1) The calculations and other code in the DRS are modified from code in previously existing Iolite DRS (“Sr_isotopes_CaAr” and “Sr_isotopes_REE”). Dr. Marillo-Sialer provided additional help with some code for importing half-masses. Therefore, the code is not original and is not entirely my own. I have simply modified it and/or combined different calculations.

    2) Much of the code is unchanged from existing Iolite code, but I’ve included the complete DRS to avoid problems with cutting and pasting.

    3) The code that imports the data and defines waves is set up for the particular collector configuration on the Nu Plasma II that we use in Aarhus. From low mass to high mass: L5-82; L4-83; L3-83.5; L2-84; L1-84.5; Ax-85; H1-85.5; H2-86; H3-86.5; H4-87; H5-87.5; H6-88; H7-89. If you use a different configuration, you may need to change the lines of code that look like this:
    Wave SrCaArYbLu88 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel(“88”, “wavename”, OptionalElementName = “H6″),name=”SrCaArYbLu88”)

    4) The calculations in this DRS have not been thoroughly checked over by anyone else, and I am not guaranteeing their accuracy. I strongly suggest that you check over the calculations (everything below the line: //CALCULATIONS BELOW HERE), as it’s possible that I’ve made errors. If you find and error, please announce it on the forum.

    5) This is DRS is not necessarily the best approach for you and your data. It simply addresses several different interferences that we have observed in our measurements here. One caveat is that this approach seems to decrease the internal precision of measurements by ~25%, relative to other approaches that I’ve implemented (Kr peak stripping, CaAr peak stripping without REE2+ correction). I think this is due to error propagation from using small signals on 82, 83, and half-masses for calculations, but I’m not exactly sure.

    6) I’ve successfully used this DRS with Iolite version 3.32 within Igor Pro 6 on a Windows 7 machine. It accurately reproduces the reference 87Sr/86Sr ratios of several in-house plagioclase standards that we use in Aarhus. However, corrected 87Sr/86Sr ratios for high Rb/Sr materials are routinely too high. This can potentially be mitigated by making adjustments to the Rb mass bias factor, but the DRS does not automatically do this. From my experience, this DRS works well for plagioclase and low Rb/Sr rock standards, but I’m not suggesting that it will be a universally effective approach to LA Sr data from all materials. You should apply this DRS to data from standards to see if it yields accurate and adequately precise results for your purposes.

    Here is an explanation of the DRS. This explanation is not exhaustive, so you may need to look into the calculations for more details:

    This DRS corrects for both CaAr and REE2+ in two different approaches. The first approach subtracts REE2+ from mass 82 and uses the residual signal on 82 to correct for CaAr. Doubly charged Dy and Er both interfere on mass 82. We monitor 167Er2+ on mass 83.5, but cannot measure an interference-free Dy2+ mass on our collector configuration. Therefore, this approach uses a canonical Dy/Er ratio to subtract the Dy interference from 82. The default ratio in the DRS is 1.8 (very loosely constrained by Dy/Er of several rock standards), but this can be changed in the DRS tab. The second approach subtracts Er2+ from mass 83 (Er is the only REE2+ interference on this mass) and uses the residual signal on mass 83 to correct for CaAr. The downside to this is 83(CaAr) is only about 1/5 as abundant as 82(CaAr). Once REE2+ have been subtracted from the 82 or 83 peaks, all masses of interest are corrected for CaAr and REE2+ interferences. For both approaches, Kr interferences are corrected using baseline subtraction, which should be okay as long as whichever Iolite baseline spline you use accurately reflects changes in the background. A preliminary mass bias factor is calculated using the baseline-subtracted (bot not interference-corrected) 88 and 86 signals (assuming these are dominated by Sr). A refined Sr mass bias factor is calculated after interferences have been stripped from 88Sr and 86Sr. This mass bias factor is applied throughout, but you can make adjustments to the Rb and CaAr miass bias factors in the DRS tab if you see fit. The “RbBeta_adjustment” and “CaArBeta_adjustment” factors simply add or subtract (if negative) from the Sr mass bias factor. The same could be done for the REE mass bias factor (it can even be directly calculated using 171Yb2+ and 173Yb2+ on the half masses), but just haven’t gotten around to implementing this. I doubt that it is important for low REE/Sr material (maybe important for apatite). Due to my lack of creativity, results calculated using the first CaAr correction method are named “result” (e.g. Sr87_86_Corr), and results calculated using the second method are named “result_b” (e.g. Sr87_86_Corr_b). Corrected 87Sr/86Sr ratios can also be normalized to a primary standard. The default is BHVO-2 glass, although I don’t recommend that you use this to normalize plagioclase data. Age-corrected 87Sr/86Sr are calculated using the present-day corrected (but not normalized) 87Sr/86Sr and 87Rb/86Sr and 87Rb decay constant from Villa et al. (2015; GCA). You can set the age of the sample in the DRS tab (in Ma). I suggest that you don’t change the “Rb87_85_reference” or “Sr88_86_reference” values unless you have a good reason to do so. I should probably take those options out of the DRS tab.

    
    #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= "2.0"  	//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 							="m88_H6"
    	GlobalString				ReferenceStandard 						="G_BHVO2G"
    	GlobalString				DefaultIntensityUnits					="V"
    	//**** 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 "//"
    	GlobalVariable			MaskThreshold 							=0.1
    	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			RbBeta_adjustment									= 0			//add or subtract this factor to SrBeta and apply to Rb
    	GlobalVariable			CaArBeta_adjustment									= 0			//add or subtract this factor to SrBeta and apply to CaAr
    	GlobalVariable			Rb87_85_reference 									= 0.385710	//Konter & Storm (2014)
    	GlobalVariable			Sr88_86_reference									= 8.37520938 	//Konter & Storm (2014)
    	GlobalVariable			Age												= 0
    	GlobalVariable			Dy_Er											= 1.8
    	GlobalString				PropagateSplineErrors					="Yes"
    	//**** If you'd like to set up some preferred settings for the report window you can set these here too
    	GlobalString				Report_DefaultChannel				="Sr87_86_Corr"
    	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)
    	Print "DRS initialised:  Sr laser ablation analysis module for the Nu plasma MC-ICP-MS, using CaAr interference corrections, \"" + nameofthisDRS + "\", Version " + DRS_Version_No + "\r"
    End //**end of initialisation routine
    
    //****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.
    	//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, PropagateSplineErrors
    	NVar MaskThreshold, MaskEdgeDiscardSeconds, BeamSecondsSensitivity, RbBeta_adjustment, CaArBeta_adjustment, Rb87_85_reference, Sr88_86_reference, Age	, Dy_Er
    	setdatafolder $currentdatafolder
    	
    	ProgressDialog()		//Open the progress dialog
    	
    	//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;..."
    	
    	//Do we have a baseline_1 spline for the index channel, as require this to proceed further?
    	DRSAbortIfNotSpline(IndexChannel, "Baseline_1")
    	SetProgress(10, "Subtracting baselines...")
    
    	//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) //Use total beam to get total no of points - it will incorporate all input waves, so is better than the index, which will only incorporate the points where it exists
    	wave IndexOut = $InterpOntoIndexTimeAndBLSub(IndexChannel,name="IndexChannel")	//Make an output wave for Index isotope (as baseline-subtracted intensity)
    	
    	Wave Yb86_5_in=$DRS_FindInputChannel("87", "full", OptionalElementName = "H3")
    	if (numpnts(Yb86_5_in) == 0)	//Check to see if half masses were measured. If no half masses were measured, print the following message in the history
    		print "No REE corrections are possible because the required half-mass beams were not collected when analysing. Please select another data reduction scheme"
    		return 0	//and stop the DRS, as REE calculations are not possible
    	endif
    
    	//baseline subtract all input channels - explicitly stating each channel is preferable for multi-collector data.
    	Wave Y89 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("89", "wavename"),name="Y89")
    	Wave SrCaArYbLu88 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("88", "wavename", OptionalElementName = "H6"),name="SrCaArYbLu88")
    	Wave Lu87_5 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("88", "wavename", OptionalElementName = "H5"),name="Lu87_5")	
    	Wave SrRbYb87 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("87", "wavename", OptionalElementName = "H4"),name="SrRbYb87")
    	Wave Yb86_5 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("87", "wavename", OptionalElementName = "H3"),name="Yb86_5")	
    	Wave SrCaArYb86 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("86", "wavename", OptionalElementName = "H2"),name="SrCaArYb86")
    	Wave RbYbEr85 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("85", "wavename", OptionalElementName = "Ax"),name="RbYbEr85")
    	//Wave Tm84_5 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("85", "wavename", OptionalElementName = "L1"),name="Tm84_5")	
    	Wave SrCaArErYb84 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("84", "wavename", OptionalElementName = "L2"),name="SrCaArErYb84")
    	Wave Er83_5 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("84", "wavename", OptionalElementName = "L3"),name="Er83_5")
    	Wave CaArEr83 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("83", "wavename", OptionalElementName = "L4"),name="CaArEr83")
    	Wave CaArErDy82 = $InterpOntoIndexTimeAndBLSub(DRS_FindInputChannel("82", "wavename", OptionalElementName = "L5"),name="CaArErDy82")
    	wave Raw83 = $MakeioliteWave("CurrentDRS","Raw83",n=NoOfPoints)
    	Wave Mx83_in=$DRS_FindInputChannel("83", "full", OptionalElementName = "L4")
    	Raw83 = Mx83_in //This one is not baseline subtracted as it is intended to be used a reference of the raw counts
    	
    	SetProgress(30, "Calculating intermediates...")
    	//now make some waves to be used in further calculations
    	wave PFract = $MakeioliteWave("CurrentDRS","PFract",n=NoOfPoints)
    	wave CaAr82 = $MakeioliteWave("CurrentDRS","CaAr82",n=NoOfPoints)
    	wave CaAr83 = $MakeioliteWave("CurrentDRS","CaAr83",n=NoOfPoints)
    	wave PSrCaArYb86 = $MakeioliteWave("CurrentDRS","PSrCaArYb86",n=NoOfPoints)
    	wave PSrCaArYbLu88 = $MakeioliteWave("CurrentDRS","PSrCaArYbLu88",n=NoOfPoints)
    	wave BetaSr = $MakeioliteWave("CurrentDRS","BetaSr",n=NoOfPoints)
    	wave BetaRb = $MakeioliteWave("CurrentDRS","BetaRb",n=NoOfPoints)
    	wave BetaCaAr = $MakeioliteWave("CurrentDRS","BetaCaAr",n=NoOfPoints)
    	wave Sr87 = $MakeioliteWave("CurrentDRS","Sr87",n=NoOfPoints)
    	wave CaArErYb84 = $MakeioliteWave("CurrentDRS","CaArErYb84",n=NoOfPoints)
    	wave Sr84 = $MakeioliteWave("CurrentDRS","Sr84",n=NoOfPoints)
    	wave Sr88 = $MakeioliteWave("CurrentDRS","Sr88",n=NoOfPoints)
    	wave Sr86 = $MakeioliteWave("CurrentDRS","Sr86",n=NoOfPoints)
    	wave Rb87 = $MakeioliteWave("CurrentDRS","Rb87",n=NoOfPoints)
    	wave Rb85 = $MakeioliteWave("CurrentDRS","Rb85",n=NoOfPoints)
    	wave Sr8786_Uncorr = $MakeioliteWave("CurrentDRS","Sr8786_Uncorr",n=NoOfPoints)
    	wave Sr87_86_Corr = $MakeioliteWave("CurrentDRS","Sr87_86_Corr",n=NoOfPoints)
    	wave Rb87_Sr86_Corr = $MakeioliteWave("CurrentDRS","Rb87_Sr86_Corr",n=NoOfPoints)
    	wave Sr8488_Uncorr = $MakeioliteWave("CurrentDRS","Sr8488_Uncorr",n=NoOfPoints)
    	wave Sr84_88_Corr = $MakeioliteWave("CurrentDRS","Sr84_88_Corr",n=NoOfPoints)
    	wave Rb87asPPM = $MakeioliteWave("CurrentDRS","Rb87asPPM",n=NoOfPoints)
    	wave CaArErYb84asPPM = $MakeioliteWave("CurrentDRS","CaArErYb84asPPM",n=NoOfPoints)
    	wave TotalSrBeam = $MakeioliteWave("CurrentDRS","TotalSrBeam",n=NoOfPoints)
    	
    	wave PSrCaArYb86_b = $MakeioliteWave("CurrentDRS","PSrCaArYb86_b",n=NoOfPoints)
    	wave PSrCaArYbLu88_b = $MakeioliteWave("CurrentDRS","PSrCaArYbLu88_b",n=NoOfPoints)
    	wave BetaSr_b = $MakeioliteWave("CurrentDRS","BetaSr_b",n=NoOfPoints)
    	wave BetaRb_b = $MakeioliteWave("CurrentDRS","BetaRb_b",n=NoOfPoints)
    	wave BetaCaAr_b = $MakeioliteWave("CurrentDRS","BetaCaAr_b",n=NoOfPoints)
    	wave Sr87_b = $MakeioliteWave("CurrentDRS","Sr87_b",n=NoOfPoints)
    	wave CaArErYb84_b = $MakeioliteWave("CurrentDRS","CaArErYb84_b",n=NoOfPoints)
    	wave Sr84_b = $MakeioliteWave("CurrentDRS","Sr84_b",n=NoOfPoints)
    	wave Sr88_b = $MakeioliteWave("CurrentDRS","Sr88_b",n=NoOfPoints)
    	wave Sr86_b = $MakeioliteWave("CurrentDRS","Sr86_b",n=NoOfPoints)
    	wave Rb87_b = $MakeioliteWave("CurrentDRS","Rb87_b",n=NoOfPoints)
    	wave Rb85_b = $MakeioliteWave("CurrentDRS","Rb85_b",n=NoOfPoints)
    	wave Sr8786_Uncorr_b = $MakeioliteWave("CurrentDRS","Sr8786_Uncorr_b",n=NoOfPoints)
    	wave Sr87_86_Corr_b = $MakeioliteWave("CurrentDRS","Sr87_86_Corr_b",n=NoOfPoints)
    	wave Rb87_Sr86_Corr_b = $MakeioliteWave("CurrentDRS","Rb87_Sr86_Corr_b",n=NoOfPoints)
    	wave Sr8488_Uncorr_b = $MakeioliteWave("CurrentDRS","Sr8488_Uncorr_b",n=NoOfPoints)
    	wave Sr84_88_Corr_b = $MakeioliteWave("CurrentDRS","Sr84_88_Corr_b",n=NoOfPoints)
    	wave Rb87asPPM_b = $MakeioliteWave("CurrentDRS","Rb87asPPM_b",n=NoOfPoints)
    	wave CaArErYb84asPPM_b = $MakeioliteWave("CurrentDRS","CaArErYb84asPPM_b",n=NoOfPoints)
    	wave TotalSrBeam_b = $MakeioliteWave("CurrentDRS","TotalSrBeam_b",n=NoOfPoints)
    	
    	
    	//CALCULATIONS BELOW HERE
    		
    	PFract = (Ln(Sr88_86_reference/(SrCaArYbLu88/SrCaArYb86)))/(Ln(87.9056/85.9093))										//calculate a preliminary fractionation factor
    	
    	//The following equations subtract CaAr using the signal on 82, itself corrected for REE2+ using the Er83_5 signal and canonical Dy/Er (because we don't measure Dy on a half mass).
    	CaAr82 = CaArErDy82 - (Er83_5 * 1.601 /22.869) / ((81.96460 / 83.46619) ^ PFract) - ((Er83_5 * 1.601 /22.869) / ((81.96460 / 83.46619) ^ PFract) * Dy_Er * 28.260)
    	
    	PSrCaArYb86 = (SrCaArYb86 - ((CaAr82 * .004 / .647) / ((85.9160721 / 81.9210049) ^ PFract)) - (Yb86_5 * 21.68 / 16.103) / ((85.968195 / 86.46911) ^ PFract))	//use this preliminary fract in stripping CaAr and REE from Sr 86
    	PSrCaArYbLu88 = (SrCaArYbLu88 - ((CaAr82 * .187 / .647) / ((87.9149151 / 81.9210049) ^ PFract)) - (Yb86_5 * 12.996 / 16.103) / ((87.97129 / 86.46911) ^ PFract) - (Lu87_5 * 2.599 /97.401) / ((87.971345 / 87.47039) ^ PFract))		//and Sr 88
    	BetaSr = (Ln(Sr88_86_reference/(PSrCaArYbLu88/PSrCaArYb86)))/(Ln(87.9056/85.9093))									//Use these CaAr and REE stripped Sr 86 and Sr 88 values to calculate a refined fractionation factor
    	BetaRb = BetaSr + RbBeta_adjustment					//allows for a modification of the Rb fractionation factor relative to the Sr fract. We have never used anything but 1 (i.e. no modification)																																
    	Rb87 = ((RbYbEr85-(Yb86_5 * 2.982 / 16.103) / ((84.967385 / 86.46911) ^ BetaSr) - (Er83_5 * 14.910 /22.869) / ((84.96774 / 83.46619) ^ BetaSr)) * Rb87_85_reference) / ((86.90918 / 84.9118) ^ BetaRb)				//use the Rb fract to calculate the amount of Rb on mass 87
    	Sr87 = (SrRbYb87 - Rb87 - (Yb86_5 * 32.026 /16.103) / ((86.96943 / 86.46911) ^ BetaSr))																															//subtract this amount of Rb 87 from the 87 beam
    	BetaCaAr = BetaSr + CaArBeta_adjustment			//allows for a modification of the CaAr fractionation factor relative to the Sr fract. We have never used anything but 1 (i.e. no modification)																																
    	CaArErYb84 = (CaAr82 * ((2.086/0.647) / ((83.917989 / 81.921122) ^ BetaCaAr)))+(Er83_5 * 26.978 /22.869)/ ((83.96619 / 83.46619) ^ BetaSr)+(Yb86_5 * 0.123 /16.103) / ((83.96694 / 86.46911) ^ BetaSr)								//calculate the amount of CaAr on mass 84
    	Sr84 = (SrCaArErYb84 - CaArErYb84)																							//subtract this amount from the 84 beam						
    	Sr88 = (SrCaArYbLu88 - ((CaAr82 * .187 / .647) / ((87.9149151 / 81.9210049) ^ BetaSr))-(Yb86_5 * 12.996 /16.103) / ((87.97129 / 86.46911) ^ BetaSr) - (Lu87_5 * 2.599 /97.401) / ((87.971345 / 87.47039) ^ BetaSr))					//calculate a final CaAr and REE corrected 88 beam
    	Sr86 = (SrCaArYb86 - ((CaAr82 * .004 / .647) / ((85.9160721 / 81.9210049) ^ BetaSr))-(Yb86_5 * 21.68 /16.103) / ((85.968195 / 86.46911) ^ BetaSr))					//calculate a final CaAr and REE corrected 86 beam
    	
    	//Create a mask wave, to be applied to the ratio-derived outputs . If a mask is not used the ratios become extremely noisy at low signals.
    	Wave Mask=$DRS_CreateMaskWave(IndexOut,MaskThreshold,MaskEdgeDiscardSeconds,"Export_Mask","StaticAbsolute") 
    
    	//Calculate the output waves
    	SetProgress(65, "Calculating output channels...")	
    	Sr8786_Uncorr = (SrRbYb87 / SrCaArYb86) * (86.9089 / 85.9093) ^ BetaSr										* Mask	
    	Sr87_86_Corr = (Sr87 / Sr86) * (86.9089 / 85.9093) ^ BetaSr													* Mask	
    	Rb87_Sr86_Corr = (Rb87 / Sr86) * (86.9089 / 85.9093) ^ BetaSr													* Mask
    	Sr8488_Uncorr = (SrCaArErYb84 / SrCaArYbLu88) * (83.9134 / 87.9056) ^ BetaSr									* Mask	
    	Sr84_88_Corr = (Sr84 / Sr88) * (83.9134 / 87.9056) ^ BetaSr													* Mask	
    	Rb87asPPM = (Rb87 / SrRbYb87) * 1000000																					* Mask	
    	CaArErYb84asPPM = (CaArErYb84 / SrCaArErYb84) * 100000																				* Mask					
    	TotalSrBeam = Sr88 + Sr84 + Sr86 + Sr87																					* Mask				
    	
    	
    	//The following equations subtract CaAr using the signal on 83, itself corrected for REE2+ using the Er83_5 signal.			
    	CaAr83 = CaArEr83 -  (Er83_5 * 33.503 /22.869) / ((82.96514975 / 83.46619) ^ PFract)
    	
    	PSrCaArYb86_b = (SrCaArYb86 - ((CaAr83 * .004 / .135) / ((85.9160721 / 82.921150) ^ PFract)) - (Yb86_5 * 21.68 / 16.103) / ((85.968195 / 86.46911) ^ PFract))	//use this preliminary fract in stripping CaAr and REE from Sr 86
    	PSrCaArYbLu88_b = (SrCaArYbLu88 - ((CaAr83 * .187 / .135) / ((87.9149151 / 82.921150) ^ PFract)) - (Yb86_5 * 12.996 / 16.103) / ((87.97129 / 86.46911) ^ PFract) - (Lu87_5 * 2.599 /97.401) / ((87.971345 / 87.47039) ^ PFract))		//and Sr 88
    	BetaSr_b = (Ln(Sr88_86_reference/(PSrCaArYbLu88_b/PSrCaArYb86_b)))/(Ln(87.9056/85.9093))									//Use these CaAr and REE stripped Sr 86 and Sr 88 values to calculate a refined fractionation factor
    	BetaRb_b = BetaSr_b + RbBeta_adjustment					//allows for a modification of the Rb fractionation factor relative to the Sr fract. We have never used anything but 1 (i.e. no modification)																																
    	Rb87_b = ((RbYbEr85-(Yb86_5 * 2.982 / 16.103) / ((84.967385 / 86.46911) ^ BetaSr_b) - (Er83_5 * 14.910 /22.869) / ((84.96774 / 83.46619) ^ BetaSr_b)) * Rb87_85_reference) / ((86.90918 / 84.9118) ^ BetaRb_b)				//use the Rb fract to calculate the amount of Rb on mass 87
    	Sr87_b = (SrRbYb87 - Rb87_b - (Yb86_5 * 32.026 /16.103) / ((86.96943 / 86.46911) ^ BetaSr_b))																															//subtract this amount of Rb 87 from the 87 beam
    	BetaCaAr_b = BetaSr_b + CaArBeta_adjustment			//allows for a modification of the CaAr fractionation factor relative to the Sr fract. We have never used anything but 1 (i.e. no modification)																																
    	CaArErYb84_b = (CaAr83 * ((2.086/0.135) / ((83.917989 / 82.921150) ^ BetaCaAr_b)))+(Er83_5 * 26.978 /22.869)/ ((83.96619 / 83.46619) ^ BetaSr_b)+(Yb86_5 * 0.123 /16.103) / ((83.96694 / 86.46911) ^ BetaSr_b)								//calculate the amount of CaAr and REE on mass 84
    	Sr84_b = (SrCaArErYb84 - CaArErYb84_b)																							//subtract this amount from the 84 beam						
    	Sr88_b = (SrCaArYbLu88 - ((CaAr83 * .187 / .135) / ((87.9149151 / 82.921150) ^ BetaSr_b))-(Yb86_5 * 12.996 /16.103) / ((87.97129 / 86.46911) ^ BetaSr_b) - (Lu87_5 * 2.599 /97.401) / ((87.971345 / 87.47039) ^ BetaSr_b))					//calculate a final CaAr and REE corrected 88 beam
    	Sr86_b = (SrCaArYb86 - ((CaAr83 * .004 / .135) / ((85.9160721 / 82.921150) ^ BetaSr_b))-(Yb86_5 * 21.68 /16.103) / ((85.968195 / 86.46911) ^ BetaSr_b))					//calculate a final CaAr and REE corrected 86 beam
    	
    	//Create a mask wave, to be applied to the ratio-derived outputs . If a mask is not used the ratios become extremely noisy at low signals.
    	Wave Mask=$DRS_CreateMaskWave(IndexOut,MaskThreshold,MaskEdgeDiscardSeconds,"Export_Mask","StaticAbsolute") 
    
    	//Calculate the output waves
    	SetProgress(65, "Calculating output channels...")	
    	Sr8786_Uncorr_b = (SrRbYb87 / SrCaArYb86) * (86.9089 / 85.9093) ^ BetaSr_b										* Mask	
    	Sr87_86_Corr_b = (Sr87_b / Sr86_b) * (86.9089 / 85.9093) ^ BetaSr_b												* Mask	
    	Rb87_Sr86_Corr_b = (Rb87_b / Sr86_b) * (86.9089 / 85.9093) ^ BetaSr_b													* Mask
    	Sr8488_Uncorr_b = (SrCaArErYb84 / SrCaArYbLu88) * (83.9134 / 87.9056) ^ BetaSr_b									* Mask	
    	Sr84_88_Corr_b = (Sr84_b / Sr88_b) * (83.9134 / 87.9056) ^ BetaSr_b													* Mask	
    	Rb87asPPM_b = (Rb87_b / SrRbYb87) * 1000000																					* Mask	
    	CaArErYb84asPPM_b = (CaArErYb84_b / SrCaArErYb84) * 100000																				* Mask					
    	TotalSrBeam_b = Sr88_b + Sr84_b + Sr86_b + Sr87_b																					* Mask				
    	
    
    	
    	
    	//Add all of the created outputs to the list of outputs (only waves in this list will be displayed or exported)
    	ListOfIntermediateChannels="Sr88;Sr86;Sr84;Rb85;Rb87asPPM;CaArEr84asPPM;Sr88_b;Sr86_b;Sr84_b;Rb85_b;Rb87asPPM_b;CaArEr84asPPM_b;"
    	ListOfOutputChannels="TotalSrBeam;Sr8786_Uncorr;Sr87_86_Corr;Rb87_Sr86_Corr;Sr8488_Uncorr;Sr84_88_Corr;BetaSr;TotalSrBeam_b;Sr8786_Uncorr_b;Sr87_86_Corr_b;Rb87_Sr86_Corr_b;Sr8488_Uncorr_b;Sr84_88_Corr_b;BetaSr_b;"
    
    	//Now check if there are integration periods for the selected reference standard, and if so, generate standard-normalised ratios
    	DRSAbortIfNotSpline("Sr87_86_Corr", ReferenceStandard)
    	SetProgress(70, "Calculating reference material corrected results...")
    	//first, map the spline onto index time
    	wave StdSpline_Sr87_86 = $InterpSplineOntoIndexTime("Sr87_86_Corr", ReferenceStandard)
    	//get the value from the standard
    	variable StdValue_Sr87_86 = GetValueFromStandard("87Sr_86Sr",ReferenceStandard)
    	If(numtype(StdValue_Sr87_86) == 2)		//There appears to be a couple of formats for 87Sr/86Sr in the std files. If we failed to get the value in the above line, try again
    		StdValue_Sr87_86 = GetValueFromStandard("87Sr/86Sr",ReferenceStandard)
    	Endif
    	//make a wave for the standard-corrected results
    	wave StdCorr_Sr87_86 =$MakeioliteWave("CurrentDRS","StdCorr_Sr87_86",n=NoOfPoints)
    	StdCorr_Sr87_86 = Sr87_86_Corr * StdValue_Sr87_86 / StdSpline_Sr87_86
    	ListOfOutputChannels+="StdCorr_Sr87_86;"
    		
    	//Now generate age-corrected values (using the observed Rb/Sr ratio)
    	Wave Sr87_86_AgeCorr=$MakeioliteWave("CurrentDRS","Sr87_86_AgeCorr",n=NoOfPoints)
    	Sr87_86_AgeCorr = Sr87_86_Corr - (Rb87_Sr86_Corr *(exp(0.000013972 * Age)-1))			
    	ListOfOutputChannels+="Sr87_86_AgeCorr;"
    	
    	//Now propagate errors (if required)
    	SetProgress(80, "Propagating errors...")
    	if(cmpstr(PropagateSplineErrors, "Yes") == 0)	//if the user wants to propagate spline errors
    		Propagate_Errors("All", "StdCorr_Sr87_86", "Sr87_86_Corr", ReferenceStandard)
    	endif
    	SetProgress(100, "Finished")
    	
    	//Now check if there are integration periods for the selected reference standard, and if so, generate standard-normalised ratios
    	DRSAbortIfNotSpline("Sr87_86_Corr_b", ReferenceStandard)
    	SetProgress(70, "Calculating reference material corrected results...")
    	//first, map the spline onto index time
    	wave StdSpline_Sr87_86_b = $InterpSplineOntoIndexTime("Sr87_86_Corr_b", ReferenceStandard)
    	//get the value from the standard
    	
    	//make a wave for the standard-corrected results
    	wave StdCorr_Sr87_86_b =$MakeioliteWave("CurrentDRS","StdCorr_Sr87_86_b",n=NoOfPoints)
    	StdCorr_Sr87_86_b = Sr87_86_Corr_b * StdValue_Sr87_86 / StdSpline_Sr87_86_b
    	ListOfOutputChannels+="StdCorr_Sr87_86_b;"
    		
    	//Now generate age-corrected values (using the observed Rb/Sr ratio)
    	Wave Sr87_86_AgeCorr_b=$MakeioliteWave("CurrentDRS","Sr87_86_AgeCorr_b",n=NoOfPoints)
    	Sr87_86_AgeCorr = Sr87_86_Corr_b - (Rb87_Sr86_Corr_b *(exp(0.000013972 * Age)-1))			
    	ListOfOutputChannels+="Sr87_86_AgeCorr_b;"
    	
    	//Now propagate errors (if required)
    	SetProgress(80, "Propagating errors...")
    	if(cmpstr(PropagateSplineErrors, "Yes") == 0)	//if the user wants to propagate spline errors
    		Propagate_Errors("All", "StdCorr_Sr87_86_b", "Sr87_86_Corr_b", ReferenceStandard)
    	endif
    	SetProgress(100, "Finished")
    	
    end   //****End of DRS function.  Write any required external sub-routines below this point****
    
    //****Start Export data function (optional).  If present in a DRS file, this function is called by the export Stats routine when it is about to save the export stats text matrix to disk.
    Function ExportFromActiveDRS(Stats,NameOfPathToDestinationFolder) //this line must be as written here
    	wave/T Stats //will be a wave reference to the stats text wave that is about to be saved
    	String NameOfPathToDestinationFolder //will be the name of the path to the destination folder for this export.
    	//This routine allows the segment/stats export routine to be intercepted if there are DRS-specific actions to be undertaken during data export.  For now, for U-Th there is just this warning:
    end	//end of DRS intercept of data export - export routine will now save the (~altered) stats wave in the folder it supplied.
    
    Function AutoBaselines(buttonstructure) //Build the main display and integration window --- This is based off a button, so has button structure for the next few lines
    	STRUCT WMButtonAction&buttonstructure
    	if( buttonstructure.eventCode != 2 )
    		return 0  // we only want to handle mouse up (i.e. a released click), so exit if this wasn't what caused it
    	endif  //otherwise, respond to the popup click
    	ClearAllTraces()
    	AutoTrace(0, "H6_m88", -.008, -.0052)	//see the autotrace function for what these mean.	set both max and min to zero for autoscale.
    	AutoTrace(1, "H4_m87", -0.01, -.0048, extraflag = "Primary")	//see the autotrace function for what these mean.
    	AutoTrace(2, "H2_m86", -.009, -.0029)	//see the autotrace function for what these mean.
    	AutoTrace(3, "Ax_m85", -0.01, -.0007)	//see the autotrace function for what these mean.
    	AutoTrace(4, "L2_m84", -.0075, -.0006, extraflag = "Right")	//see the autotrace function for what these mean.
    	AutoTrace(5, "L4_m83", -.007, -.0012)	//see the autotrace function for what these mean.
    	AutoTrace(6, "L5_m82", -.0067, -.00042, extraflag = "Hidden")	//see the autotrace function for what these mean.
    end
    
    Function AutoIntermediates(buttonstructure) //Build the main display and integration window --- This is based off a button, so has button structure for the next few lines
    	STRUCT WMButtonAction&buttonstructure
    	if( buttonstructure.eventCode != 2 )
    		return 0  // we only want to handle mouse up (i.e. a released click), so exit if this wasn't what caused it
    	endif  //otherwise, respond to the popup click
    	ClearAllTraces()
    	AutoTrace(0, "Sr87_86_Corr", 0, 0, extraflag = "Primary")	//see the autotrace function for what these mean.	set both max and min to zero for autoscale.
    	AutoTrace(2, "Sr8786_Uncorr", 0, 0, extraflag = "Hidden")	//see the autotrace function for what these mean.	set both max and min to zero for autoscale.
    	AutoTrace(2, "TotalSrBeam", 0, 0)	//see the autotrace function for what these mean.
    	AutoTrace(4, "Rb87asPPM", 0, 10000, extraflag = "Right")	//see the autotrace function for what these mean.
    	AutoTrace(8, "Sr84_88_Corr", 0.05, 0.08)	//see the autotrace function for what these mean.
    end
Viewing 3 replies - 1 through 3 (of 3 total)
  • Author
    Replies
  • #8277

    Noreen Evans
    Participant

    Graham, this is awesome! We too are just starting down the Sr isotope path with our NPII and this post has saved me hours and hours of coding and troubleshooting. I’ll let you know how we go but I have managed to export the DRS and compile it in Iolite 3.4 (Igor 6.37). Once we get all the masses on the correct detectors, we’ll give it a shot.
    Thanks again for sharing your hard work with others and thanks too to Tephy for all her help with my own feeble attempt at amending the existing Sr isotope REE DRS. Cheers, Noreen.

    #8278

    Noreen Evans
    Participant

    Graham, one thing I noted. You list 85.5 as sitting on H1 but you don’t define a wave for 85.5. So you just collect on 85.5 in case you want to utilize the 171Yb2+ in future? Maybe I missed something here…

    #8360

    Graham Hagen-Peter
    Participant

    Hey Noreen,

    I’m glad that this is useful for you. You’re right about the wave for 85.5. The reason that the DRS doesn’t import this channel is that mass 86.5 (173Yb2+) is used for the Yb correction. It might be useful to add this channel because it would then be possible to determine a mass bias factor (173Yb/171Yb) for the REE. I will probably try this once I find some time. Thanks for pointing that out.

    One thing that I should add is that I have decided that it’s best to apply only the CaAr correction for plagioclase data (at least for the data we’ve collected here in Aarhus). This is because there is no signal on the half-masses for most plagioclase, and I think this is responsible for the increased internal errors. However, if there is some signal on the half-masses, the REE2+ correction is probably important. It seems to work well for relatively low REE/Sr and Rb/Sr rock standards (BIR and BHVO-2).

    -Graham

Viewing 3 replies - 1 through 3 (of 3 total)

You must be logged in to reply to this topic.