Igor Code

The objective of this page is to serve a repository of useful Igor code that we keep needing for different applications. Feel free to add new code or to comment on existing code, but please don't delete old code w/o discussing with Jose.

There is a a "General Macros" igor procedure file that has a large suite of various Igor functions. It can be downloaded here. There is an Igor Help file (with clickable links) that helps to describe the various functions in this ipf that can be downloaded here. Some background: this ipf was created by Donna as a good-bye 'gift' to my NOAA friends. Two new versions of Igor have since been released, and several tools are now dated and/or obsolete. Still, I think it is a good reference for general programming techniques, and as a resource for a variety of tasks. - Donna

Doing a Robust Linear Regression

Linear regression gets pretty complex when one starts thinking about the details. There are 3 important decisions, which depend on the data which you are trying to fit:

• (A) Error in one or both coordinates:
  • (A1) Only one of the coordinates (Y) has error and the other does not (e.g. when fitting an aerosol measurement vs. time).
  • (A2) Or do both coordinates have errors of the same relative order-of-magnitude? (e.g. as when comparing two instruments that measure the same quantity)
• (B) Assumptions about the relative magnitudes of the errors in different points
  • (B1) All the points are assumed to have the same error
  • (B2) The error varies point-by-point.
• (C) Standard vs. Robust regression
  • (C1) The fit line is found by minimizing Chi^2, i.e. the sum of the squares of the deviations between each point and the regression line.
  • (C2) The fit line is found by minimizing the sum of the absolute values of the deviations between each point and the regression line.

This gives rise to 2 x 2 x 2 = 8 options to choose from when fitting your data. Which to use depends on the nature of the data and the use of the fit, and should be chosen intelligently. In practice if there is a good relationship and the points are scattered around a line with few outliers, all the methods will produce similar results. However when there are outliers or significant scatter, there can be significant differences and some of the methods are plainly incorrect for some datasets.

To do regressions in Igor in each case:

• A1-B1-C1: this is the default assumption of the linear fit
• A2-B1-C1: an approximation to this fit can be done with the /ODR=2 option in Igor.
• A1-B1-C2: this can be done in Igor, but you need to define a new custom fit function (which can be just a line) and then set V_FitOptions = 2 in the command line. (Going back to V_FitOptions = 0 gets you back to the regular A1-B1-C1 fits)
• A2-B1-C2: this can be accomplished approximately by doing A1-B1-C2 for Y vs X, then for X vs Y, and averaging the fitting parameters appropriately (and the fitting parameter uncertainties in quadrature). There are some papers that show that this is approximately correct.
• For all the B2 options, choose the fitting weights in the Igor Curve Fitting tab

Regridding a Wave into the Time Axis of Another Wave

To regrid a wave into another framework, use the following functions. For example you might have one curve, defined by a certain x- and y-wave (xwave1, ywave1). Let's say there are 10 points in these waves, and the maximum is 20 (and extreme case of length discrepancy). Another curve is defined by two other x- and y-waves, also with a maximum of 20, but with 20 points each (xwave2, ywave2). It is important that they have approximately the same bound values, or it will have a hard time interpolating between them.

If I want to look at the first curve in terms of the second x-wave, you have to regrid the first y-wave into the x-wave of the second. This can be done in two ways, and that is dependent on which inintial set of waves has more data points. If the wave you want to regrid into has more points than the original, then use the first of these functions (Interp_1). If the opposite is true (wave you want to regrid to has less points) use the second function (Interp_2). They can be run as is, after changing the wave names appropriately. The [p] refers to each point, and should stay with this syntax (replacing it with a number refers only to that specified cell).

Function Interp_1 ()

Wave xwave1, ywave1, xwave2 //ywave2 is not used

duplicate/o xwave2 ywave_regrid

ywave_regrid = ywave1 [binarysearchinterp(xwave1, xwave2[p])]

appendtotable ywave_regrid

End Function

Function Interp_2 ()

Wave xwave1, ywave2, xwave2

duplicate/o xwave1 ywave_regrid

ywave_regrid = ywave1 [fAverageXY(xwave2, xwave1, xwave2[p-1], xwave2[p])]

End Function

Compiled by Qi Zhang and Alex Huffman, 2004

March 2012 update

The routine Interp_3 below is a modification of interp_1, which can deal with gaps on the data series. Adapted by Anna Ripoll and Jose Jimenez in March 2012.

Function Interp_3 (xwave1, ywave1, xwave2, maxdur_secs) // Regridding a Wave with gaps into the Time Axis of Another Wave without gaps

Variable maxdur_secs // We need to create a variable which define the maximum duration of the gap in seconds

Wave xwave1, ywave1, xwave2 // xwave1 is the date/time of the initial wave, ywave1 is the data of the initial wave, xwave2 is the date/time of the final wave and ywave2 is not used

variable i1,i2 // i1 is the point from xwave1 and i2 is the point from xwave2

duplicate/o xwave2, ywave_regrid2 // ywave_regrid2 is the new data wave

ywave_regrid2 = NaN // it just cleans the ywave_regrid2 in case there is anyone already

ywave_regrid2 = ywave1 [binarysearchinterp(xwave1, xwave2[p])] // it regrid the ywave1 into the xwave2

// Now we want to stop the regridding process when there is a gap, so we need to filter out

for(i1=0; i1<NumPnts(xwave1)-1; i1 += 1) // For all points of the xwave1

If ((xwave1[i1+1]-xwave1[i1])>=maxdur_secs) // If the xwave1 points are further from each other than the maximum duration gap we created

for(i2=0; i2<NumPnts(xwave2)-1; i2 += 1) // For all the xwave2 points which are inside the gap

If ((xwave2[i2]>(xwave1[i1]+maxdur_secs/2)) && (xwave2[i2]<(xwave1[i1+1]-maxdur_secs/2))) // If they are far from the bound (you can choose how far, in this case maxdur_secs/2)

ywave_regrid2[i2]=NaN // then we decide don't have that data

Endif

EndFor

EndIf

endfor

appendtotable ywave_regrid2

End Function

Looking for an Object in Experiments in a Directory

Function findThisObjectInExpInPath(ObjectType, ObjectName, PathName, ThisSubFolder)

//imu 3 Sept 2009

// This function will help you look for lost objects in your experiments.

// ObjectType must be "string", "variable", or "wave".

// ObjectName must be the exact name of the object, and can be semicolon-separated list.

// Before running this function you must create an Igor Path to the folder you want to search in.

// Command line: NewPath ThisPathName "C:Program Files:WaveMetrics:Igor Pro Folder:"

// or Misc -> New Path

// Send the pathName as a string to the function, e.g, "ThisPathName"

// ThisSubFolder refers to a subfolder to check in the experiments you're looking in.

// To look in root: use "".

// To look in all folders, use "allFolders". This is a cheat for selecting the /R flag for LoadData.

// Note that if you look in allFolders, the complete directory structure of the opened file will be

// created in the current experiment, even if the object is not found.

//

// You can further customize this function for overwriting old data and storing loaded data in a new dataFolder

// by including the /O and /T flags, respectively, with LoadData

string ObjectName, ObjectType, PathName, ThisSubFolder

variable idex, LoadFlag

// check that ObjectType is valid and set LoadFlag for LoadData

// Note that LoadData could look for multiple objectTypes at once, but I'm not using it that way

if(stringmatch(ObjectType, "string"))

LoadFlag = 4

elseif(stringmatch(ObjectType, "variable"))

LoadFlag = 2

elseif(stringmatch(ObjectType, "wave"))

LoadFlag = 1

else

abort "ObjectType must be \"string\", \"variable\", or \"wave\"."

endif

// get ; -separated list of file names in PathName

string ExperimentsInFolderList = indexedFile($PathName, -1, ".pxp")

For(idex = 0; idex < ItemsInList(ExperimentsInFolderList); idex += 1)

// Searching in a Subfolder or in root:

if(strlen(ThisSubFolder) == 0) // root

LoadData/J=objectName /L=(loadFlag)/P=$pathName StringFromList(idex,ExperimentsInFolderList)

elseif(stringmatch(ThisSubFolder, "allFolders")) // recursive look in all subfolders

LoadData/J=objectName/L=(loadFlag)/P=$pathName/R StringFromList(idex,ExperimentsInFolderList)

else // subfolder provided

LoadData/J=objectName /L=(loadFlag)/P=$pathName/S=ThisSubFolder StringFromList(idex,ExperimentsInFolderList)

endif

endfor

end

Other Tricks

• Making a radical in an axis label: \B\S\F'\F'Wingdings'l

RK4 and Euler's Methods

The functions below are adapted from Numerical Recipes and allow the integration of a system of ordinary differential equations in time using the Euler method and the 4th Order Runge-Kutta Method. Courtesy of Christoph Knote (NCAR).

// *******************************************************************************************************************

// HW4.4 // Integration of ordinary differential equations using RK4 and Euler // adapted by Christoph Knote on Jan 2013 from: Press et al., Numerical Recipes, 1st edition, SUBROUTINE RK4, p. 553 // extended by returning rates for different processes / species // Further adapted by J.L. Jimenez on 13-Feb-2013 to conform to naming conventions of our courses

CONSTANT MOLEC_CM3 = 2e19 // approx number of molecules per cm3 CONSTANT iNO2 = 0 // Constants for the position of each species on the Y and other arrays // to make code more understandable CONSTANT iNO = 1 CONSTANT iO3 = 2 CONSTANT NSPEC = 3 // Number of species CONSTANT NRATE = 1 // Number of rates

// Master Function to integrate an example problem in time -- can be modified to tackle other problems

Function Integrate_RK4(NO2_0_ppb, NO_0_ppb, O3_0_ppb, StepSize_s, EndTime_s, UseRK4, MakeNewPlot)

Variable NO2_0_ppb, NO_0_ppb, O3_0_ppb Variable StepSize_s, EndTime_s, UseRK4, MakeNewPlot

Variable NPointsTime = EndTime_s / StepSize_s Variable iTimeStep

Make/O/N=(NPointsTime) NO, NO2, O3, wTime NO2[0] = NO2_0_ppb * 1e-9 * MOLEC_CM3 NO[0] = NO_0_ppb * 1e-9 * MOLEC_CM3 O3[0] = O3_0_ppb * 1e-9 * MOLEC_CM3

Make/O/N=(NSPEC) wSpec, w_dSpecdt, wSpecNext Make/O/N=(NRATE) wRates

wSpec[iNO2] = NO2[0] wSpec[iNO] = NO[0] wSpec[iO3] = O3[0]

For (iTimeStep = 0 ; iTimeStep < NPointsTime ; iTimeStep += 1)

// Fill the wY array with the current species concentrations wSpec[iNO2] = NO2[iTimeStep] wSpec[iNO] = NO[iTimeStep] wSpec[iO3] = O3[iTimeStep]

// calculate the time and save it into a wave (for Derivs, and also for plotting purposes) wTime[iTimeStep] = iTimeStep * StepSize_s

// Calculate the rate of change of each species (total rate) with the derivs function Derivs(wTime[iTimeStep], wSpec, w_dSpecdt, wRates)

If (UseRK4 == 1)

// Step forward 1 time step using RK4 as a subroutine RK4(wSpec, w_dSpecdt, NSPEC, wTime[iTimeStep], StepSize_s, wSpecNext, wRates, NRATE)

// Save the output NO2[iTimeStep + 1] = wSpecNext[iNO2] NO[iTimeStep + 1] = wSpecNext[iNO] O3[iTimeStep + 1] = wSpecNext[iO3]

Else

// Use Euler's method NO2[iTimeStep + 1] = NO2[iTimeStep] + w_dSpecdt[iNO2] * StepSize_s NO[iTimeStep + 1] = NO[iTimeStep] + w_dSpecdt[iNO] * StepSize_s O3[iTimeStep + 1] = O3[iTimeStep] + w_dSpecdt[iO3] * StepSize_s

EndIf

EndFor

If (MakeNewPlot == 1) Display NO2 vs wTime AppendToGraph NO vs wTime AppendToGraph O3 vs wTime Legend EndIf

End Function

// *******************************************************************************************************************

// Function RK4 // purpose: integrate Y vs Time using 4th order Runge-Kutta // // inputs: // wY, wave of length Nspecies, values of the variables Y at one point in time // w_dydt, wave of length Nspecies, derivatives of Y at one point in time (calculated from the Function Derivs) // N_species, scalar, number of variables in Y // TimeNow, scalar, current time // dt, scalar, integration time step // // output: wYout, wave of length Nspecies, values of Y at Time + dt // // for rates of individual reactions // input and output: // rates, wave of length nrates, (production/loss/net) rates of selected reactions at x (calculated from DERIVS) // nrates, scalar, number of variables in rates

// requires: // Function DERIVS(TimeNow, wY, w_dydt, rates) to calculate dydx, rates //

Function RK4(wSpec, w_dSpecdt, Nspecies, TimeNow, dt, wSpecOut, wRates, NRates) // main input / output variables as in Numerical Recipes // Note that you need to call Derivs outside of this function once and pass the variables to RK4 // All times are in s, concentrations in molec cm-3, and rates in molec cm-3 s-1

Wave wSpec, w_dSpecdt Variable Nspecies, TimeNow, dt Wave wSpecOut

// Variables to keep track of the invidiual rates Wave wRates Variable NRates

// variables used internally Variable half_dt = 0.5 * dt Variable Sixth_dt = dt / 6.0 Variable iSpec // index for the number of species Variable iRates // index for the rates

Duplicate/O wSpec, wSpec_Temp1, w_dSpec_Temp1, w_dSpec_Temp2, wRates, wRates_Temp1, wRates_Temp2

// sub-step 1 -- Euler 1/2 step For (iSpec = 0 ; iSpec < NSPEC ; iSpec += 1)

wSpec_Temp1[iSpec] = wSpec[iSpec] + half_dt * w_dSpecdt[iSpec]

EndFor

// sub-step 2 -- another Euler 1/2 step using the info from step 1 Derivs(TimeNow + half_dt, wSpec_Temp1, w_dSpec_Temp1, wRates_Temp1)

For (iSpec = 0 ; iSpec < NSPEC ; iSpec += 1) wSpec_Temp1[iSpec] = wSpec[iSpec] + half_dt * w_dSpec_Temp1[iSpec] EndFor

// sub-step 3 Derivs(TimeNow + half_dt, wSpec_Temp1, w_dSpec_Temp2, wRates_Temp2)

For (iSpec = 0 ; iSpec < NSPEC ; iSpec += 1) wSpec_Temp1[iSpec] = wSpec[iSpec] + dt * w_dSpec_Temp2[iSpec] w_dSpec_Temp2[iSpec] = w_dSpec_Temp1[iSpec] + w_dSpec_Temp2[iSpec] EndFor

For (iRates = 0 ; iRates < NRATE ; iRates += 1) wRates_Temp2[iRates] = wRates_Temp1[iRates] + wRates_Temp2[iRates] EndFor // <-

// sub-step 4 DERIVS(TimeNow + dt, wSpec_Temp1, w_dSpec_Temp1, wRates_Temp1)

// Final step -- accumulate increments with proper weights For (iSpec = 0 ; iSpec < NSPEC ; iSpec += 1) wSpecOut[iSpec] = wSpec[iSpec] + Sixth_dt * (w_dSpecdt[iSpec] + w_dSpec_Temp1[iSpec] + 2 * w_dSpec_Temp2[iSpec]) EndFor

for (iRates = 0 ; iRates < NRATE ; iRates += 1) wRates[iRates] = 1/6 * (wRates[iRates] + wRates_Temp1[iRates] + 2 * wRates_Temp2[iRates]) endfor

// remove temporary waves KillWaves wSpec_Temp1, w_dSpec_Temp1, w_dSpec_Temp2, wRates_Temp1, wRates_Temp2

End

// *******************************************************************************************************************

// Function to calculate the derivative of a variable vs time at a given point in time // This function is called several times from RK4 Function Derivs(TimeNow, wSpec, w_dSpecdt, wRates) Variable TimeNow Wave wSpec, w_dSpecdt, wRates

Variable T = 300

// Example rates -- note that you will need to choose the rates appropriate for your problem

// NO + O3 Variable k_NO_O3 = 3.0e-12*exp(-1500/T)

// O3 photolysis Variable JO3_s1 = 0.0001

w_dSpecdt[iNO2] = + k_NO_O3 * wSpec[iNO] * wSpec[iO3] w_dSpecdt[iNO] = - k_NO_O3 * wSpec[iNO] * wSpec[iO3] w_dSpecdt[iO3] = - k_NO_O3 * wSpec[iNO] * wSpec[iO3] - JO3_s1 * wSpec[iO3]

// wRates used to keep track of individual components of the rates -- not used for the time being

End Function