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

These functions 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. Adapted by Christoph Knote (NCAR) and edited by Jose-Luis Jimenez on Feb. 2013.