Source: src/idl_cvs/maken.pro
src/idl_cvs/maken.pro
NAME: MAKEN PURPOSE: Make an array of N values, linear between two given limits. CATEGORY: CALLING SEQUENCE: x = makex( first, last, num) INPUTS: first, last = array start and end values. in num = number of values from first to last. in KEYWORD PARAMETERS: OUTPUTS: x = array of values. out COMMON BLOCKS: NOTES: MODIFICATION HISTORY: Ray Sterner, 26 Sep, 1984. Johns Hopkins University Applied Physics Laboratory. Copyright (C) 1984, Johns Hopkins University/Applied Physics Laboratory This software may be used, copied, or redistributed as long as it is not sold and this copyright notice is reproduced on each copy made. This routine is provided as is without any express or implied warranties whatsoever. Other limitations apply as described in the file disclaimer.txt.
Source: src/idl_cvs/makenxy.pro
src/idl_cvs/makenxy.pro
NAME: MAKENXY PURPOSE: Make 2-d x and y coordinate arrays of specified dimensions. CATEGORY: CALLING SEQUENCE: makenxy, x1, x2, nx, y1, y2, ny, xarray, yarray INPUTS: x1 = min x coordinate in output rectangular array. in x2 = max x coordinate in output rectangular array. in nx = Number of steps in x. in y1 = min y coordinate in output rectangular array. in y2 = max y coordinate in output rectangular array. in ny = Number of steps in y. in KEYWORD PARAMETERS: OUTPUTS: xarray, yarray = resulting rectangular arrays. out COMMON BLOCKS: NOTES: MODIFICATION HISTORY: R. Sterner, 1996 Jul 11 Copyright (C) 1996, Johns Hopkins University/Applied Physics Laboratory This software may be used, copied, or redistributed as long as it is not sold and this copyright notice is reproduced on each copy made. This routine is provided as is without any express or implied warranties whatsoever. Other limitations apply as described in the file disclaimer.txt.
Source: src/idl_cvs/idlv4_to_v5.pro
src/idl_cvs/idlv4_to_v5.pro
NAME: MATCH PURPOSE: Routine to match values in two vectors. CALLING SEQUENCE: match, a, b, suba, subb, [ COUNT =, /SORT ] INPUTS: a,b - two vectors to match elements, numeric or string datatype OUTPUTS: suba - subscripts of elements in vector a with a match in vector b subb - subscripts of the positions of the elements in vector b with matchs in vector a. suba and subb are ordered such that a[suba] equals b[subb] OPTIONAL INPUT KEYWORD: /SORT - By default, MATCH uses two different algorithm: (1) the /REVERSE_INDICES keyword to HISTOGRAM is used for integer data, while a sorting algorithm is used for non-integer data. The histogram algorithm is usually faster, except when the input vectors are sparse and contain very large numbers, possibly causing memory problems. Use the /SORT keyword to always use the sort algorithm. OPTIONAL KEYWORD OUTPUT: COUNT - set to the number of matches, integer scalar SIDE EFFECTS: The obsolete system variable !ERR is set to the number of matches; however, the use !ERR is deprecated in favor of the COUNT keyword RESTRICTIONS: The vectors a and b should not have duplicate values within them. You can use rem_dup function to remove duplicate values in a vector EXAMPLE: If a = [3,5,7,9,11] & b = [5,6,7,8,9,10] then IDL> match, a, b, suba, subb, COUNT = count will give suba = [1,2,3], subb = [0,2,4], COUNT = 3 and suba[a] = subb[b] = [5,7,9] METHOD: For non-integer data types, the two input vectors are combined and sorted and the consecutive equal elements are identified. For integer data types, the /REVERSE_INDICES keyword to HISTOGRAM of each array is used to identify where the two arrays have elements in common. HISTORY: D. Lindler Mar. 1986. Fixed "indgen" call for very large arrays W. Landsman Sep 1991 Added COUNT keyword W. Landsman Sep. 1992 Fixed case where single element array supplied W. Landsman Aug 95 Converted to IDL V5.0 W. Landsman September 1997 Use a HISTOGRAM algorithm for integer vector inputs for improved performance W. Landsman March 2000 Work again for strings W. Landsman April 2000
Source: src/idl_cvs/maxstruct.pro
src/idl_cvs/maxstruct.pro
NAME: maxstruct PURPOSE: find the maximum of a structure (ignores structures within the input structure) CATEGORY: Programming, Math EXAMPLE: IDL> structure = { a: 2.3, b:10, c:[2.1,22.3], d:'Hi Mom' } IDL> print, maxstruct(structure) 22.3000 IDL> structure2 = { a: 2.3, b:10, c:[2.1,99], d:'Hi Mom' } IDL> structure = [structure, structure2] IDL> print, maxstruct(structure) 99.0000 IDL> inStruc = {a:77., b:88., c:101.} IDL> structure = { a: 2.3, b:10, c:[2.1,22.3], d:'Hi Mom', e:inStruc } IDL> print, totstruct(structure) 101.000 NOTES: will not consider values in sub-structures MODIFICATIONS 25-Mar-2008 handles arrays of structures, and go into nested structures [BD] 01-Dec-2006 Ignore non-numeric structure elements [BD]
Source: src/idl_cvs/mdir.pro
src/idl_cvs/mdir.pro
NAME: mdir PURPOSE: Directory listing of MDSplus data files with wild cards CATEGORY: MDSplus, files CALLING SEQUENCE: IDL> list = mdir(shot, tree=tree) INPUTS: shot - shot number (defaults to 0). Can be a character string with wildcards. tree - an MDSplus e.g. 'operations' (defaults to 'wf') KEYWORD PARAMETERS: inputs: tree - an MDSplus e.g. 'operations' (defaults to 'wf') dataroot - file system root for MDSplus data (default = '/nstxdata') extension - extension of file wanted, default = 'datafile' server - MDSplus server for data, if not set will use denied - if set, will not return permission-denied messages nDigits - when trying to extract a shot number from a filename, assume this number of digits. verbose - if set, will print many informational messages debug - if set, debug output will be printed outputs: allShots - tries to extract shot number from all files. date - dates of files found ctimes - creation times of files found files - actual file names found OUTPUTS: list - list of data files found EXAMPLE: IDL> print, mdir( '13000*', tree='wf' ) IDL> list = mdir( '555555*', tree='engineering', date=date) IDL> result = mdir( '9??????', tree='mptscalib', date=date, allshots=shots, $ nDigits=7, files=files, /verbose, /debug ) (the above took 18 minutes!) NOTES: if on a system with a /nstxdata directory, will spawn "find" otherwise, you will have to enter your password to ssh into the server. You may also have to respond "yes" to a question about connecting. You must have the same username on the server as the machine you are on. You may, instead, use tree_exists.pro, which is slower. If you get tired of entering your password when looking for data on a remote host, see http://w3.pppl.gov/~bdavis/swdoc/SSH_login_without_pw.txt MODIFICATION HISTORY: 03-Mar-2010 parse server name from environmental variable for remote access 19-Nov-2009 Tidied up, added keywords, let shotstring work with wildcards. 27-Mar-2009 Added ctime output for Creation time, e.g. "Mar 19 14:37" or "May 19 2005" 18-Mar-2009 Added date keyword [Bd] 15-Nov-2008 Written by Bill Davis, PPPL
Source: src/idl_cvs/mdsdat.pro
src/idl_cvs/mdsdat.pro
NAME: MDSDAT PURPOSE: "replacement" for gadat.pro. Get MDSplus signal data without having to open a tree (much faster, if a tree is specified) CATEGORY: MDSplus CALLING SEQUENCE: mdsdat, X, Y, signame, SHOT [,err=ERR] [,debug=DEBUG] $ [,xmin=XMIN] [,xmax=XMAX] $ [,npts=NPTS] KEYWORDS: (all are optional) INPUT/OUTPUT keywords: NPTS: (INPUT) Number of samples to return. Defaults to 64K. Cannot exceed 64K (if it does, it is set to 64K). (OUTPUT) Number of samples actually returned. INPUT KEYWORDS: DEBUG: If set (ie. called as gadat,/debug), debug printouts on the progress of gadat are reported XMIN: specifies the "start" time for the time window within which to return the data. Specify the time in milliseconds. If not set, defaults to -10000 ms. XMAX: specifies the "end" time for the time window within which to return the data. Specify the time in milliseconds. If not set, defaults to 100000 ms. XUNITS - if non-zero, return units of x-axis out YUNITS - if non-zero, return units of y-axis out ALWAYSOPEN - always open the tree (for when other opens or MDSconnects are done) OUTPUT KEYWORDS: ERR: PTDATA/MDSplus error code. 0 = success, anything else = failure PLABEL: plot label for pointname OUTPUTS: X: The time base of the requested pointname in seconds Y: Data from the requested signal name. See keywords NPTS and A_NPTS for details on the number of samples returned. EXAMPLE: IDL> mdsdat, iptime, ip, 'wf::IP', 120200 IDL> mdsdat, efittime, ipmeas, 'efit01::IPMEAS', 120200 IDL> COMMON BLOCKS: MDSDAT_LOCAL SIDE EFFECTS: RESTRICTIONS: MODIFICATION HISTORY: 07-May-2008 added /AlwaysOpen keyword for when other opens or MDSconnects are done outside of mdsdat. This will make it slower. 04-Feb-2008 Written by Bill Davis
Source: src/idl_cvs/mdsdatavsshot.pro
src/idl_cvs/mdsdatavsshot.pro
NAME: MDSdataVsShot PURPOSE: Print single values of MDSplus signals vs. shot number. The value can be the max or min, or at a specific time, including time of Max Ip, Max Neutrons, etc. The shot list may be reduced by certain qualifiers, such as to shots where the max of an arbitrary signal is > a given value. Output can be in a format suitable for importing into Kaleidograph, Excel, etc. CATEGORY: MDSplus, Output CALLING SEQUENCE: MDSdataVsShot, shot, nshots, signal INPUTS: MDStree - e.g. 'operations'. Default is 'NSTX' shot - MDSplus shot number, i.e, ID. signal - a single MDSplus signal name, or array of names. (ignored if Keyword "Scope" is present.) OUTPUTS: none returned KEYWORDS: (all Optional) format - output format. Defaults to '(100(g15.6, a1))' (a tab or space is written between the columns) tabs - if set, columns will be separated by tabs instead of spaces headings - if set, the signal labels are at the top of each column in row 1. scope - if present, gets the signal names and column titles from this scope file. (if present, the signal input is ignored). THIS SHOULD BE THE FULL PATHNAME OF THE FILE, or it will default to general areas. status - 1=OK, 0 means shot, signal, or scope file weren't found. units - if set, will append units to headings WhatToList - 'MaxValues', 'MinValues', (time of the following:) 'MaxNeutrons', 'MaxIP', 'MaxBetaN', 'MaxDENS', 'MaxWMHD' EXAMPLES: LIMITATION: HISTORY: 11-Feb-05 bug fix for when 1 dimension of 2-d array smaller than nsmooth 04-Jun-04 Bug fix when lots of bad shots 04-Nov-02 default to skipping a shot if any requested signal missing Written 17-Jul-02 by Bill Davis at the request of Steve Paul
Source: src/idl_cvs/mdsdeclareevent.pro
src/idl_cvs/mdsdeclareevent.pro
NAME: MDSDeclareEvent PURPOSE: set an MDS event with shot number associated with the event CATEGORY: MDSplus, Events CALLING SEQUENCE: IDL> MDSDeclareEvent, eventName [, shot_num] INPUTS: eventName - string of MDS event name, e.g., 'NE_DENSITY_CALC' shot_num - NSTX shot number (OPTIONAL -- will default to current shot) KEYWORD PARAMETERS: Optional: VERBOSE - If set, print debugging information NOCONNECT - if set, will not try to connect to an MDSplus server OUTPUTS: none -- just sets the event COMMON BLOCKS: NONE EXAMPLE: NOTES: Works from Unix or VMS MODIFICATION HISTORY: 11-Feb-04 use SETEVENT in MDSVALUE on Unix. Add connect option on unix. 20-Dec-99 Written [bd]
Source: src/idl_cvs/mdseditbranch.pro
src/idl_cvs/mdseditbranch.pro
NAME: MDSeditBranch PURPOSE: Read a branch of MDSplus values, allow editing of scalar values in a widget and write back out. CATEGORY: MDSplus CALLING SEQUENCE: MDSeditBranch, shot=shot, tree=tree, branch=branch INPUTS: KEYWORD PARAMETERS: Inputs: shot - shot number to edit tree - MDSplus treename to edit (DEF='cameras') branch - branch of tree to edit values in (DEF=".FC_1.SETTINGS") cautions - array, where, if=1, will warn user about changing (not implemented) order - if =0 will order inputs by Integers, reals & then strings comments - tips, to go after fields (not implemented) ignore - array, where, if=1, will not present that element for editing ncols - # of columns in display title - title of widget maxPerCol - maxFieldSize - if not present, takes reasonable defaults nInFirstCol - debug - if set, prints debug info and stops in a few places verbose - if set, prints information messages OUTPUTS: NONE unless user chooses to write new values to a tree (user must have priviledges to write to tree) EXAMPLE setenv cameras_path /u/bdavis/TCLcode idl IDL> MDSeditBranch, 142000, tree='cameras', branch='.FC_1.SETTINGS', $ ignore='BACKIMAGE' NOTES: For text nodes, you can enter \n for a new line MODIFICATION HISTORY: 27-Mar-2010 Written by Bill Davis, PPPL
Source: src/idl_cvs/mdseventwait.pro
src/idl_cvs/mdseventwait.pro
NAME: mdseventwait PURPOSE: Widget to Wait for an MDSplus event, or the user clicks the CANCEL button. CATEGORY: MDSplus; Event Widgets CALLING SEQUENCE: IDL> shot = mdseventwait( someEvent, cancel=cancel ) INPUTS: someEvent - (string) an MDSplus event (defualts to 'nstx_acq_done') KEYWORD PARAMETERS: Optional: cancel - if user hits cancel button, cancel = 1 (& shot=-1), else 0 group_leader - higher level widget to place this on top of OUTPUTS: shot - shot # (data) set like: IDL> setmdsshotevent,'dum',999999 COMMON BLOCKS: NONE EXAMPLE: RELATED ROUTINE: NOTES: Only works on X. MODIFICATION HISTORY: 29-Aug-00 Written by Bill Davis for David Swain
Source: src/idl_cvs/mdsgetsig.pro
src/idl_cvs/mdsgetsig.pro
NAME: MDSGETSIG PURPOSE: Get a signal, units and axes values easily from MDSplus CATEGORY: MDSplus CALLING SEQUENCE: d = MDSGETSIG( signal, SIGUNITS=fUnits, XAXIS=time, XUNITS=timeUnits) INPUTS: signal - MDSplus Signal name KEYWORD PARAMETERS: Keywords: SIGUNITS - if non-zero, return units of signal out XAXIS - if non-zero, return array of x-axis out XUNITS - if non-zero, return units of x-axis out YAXIS - if non-zero, return array of y-axis out YUNITS - if non-zero, return units of y-axis out ZAXIS - if non-zero, return array of Z-axis out ZUNITS - if non-zero, return units of Z-axis out status - if non-zero, return status (1=OK) out OUTPUTS: d - 1- or 2-D data out COMMON BLOCKS: NONE EXAMPLE MDSCONNECT, 'yourhost.pppl.gov:8501' MDSOPEN, 'testtree', 10 sig = 'onedim' data = MDSGETSIG( signal,SIGUNITS=fUnits, XAXIS=time,XUNITS=timeUnits )' PLOT, time, data, XTITLE = timeUnits, YTITLE = fUnits, $ TITLE = signal + ' - shot ' + STRTRIM( STRING(shot), 2 ) NOTES: If something like x units were NOT stored for the signal, and they are asked for, IDL will complain about an access violation, but continue. When signals are TDI combinations of others, use the units and xaxis from the first signal (could be erroneous) MODIFICATION HISTORY: 10-Aug-2009 Check for data being a number (not '*') when taking min & max 23-Jun-2008 when finding Xunits of a signal like dim_of(\...), use it 13-May-08 handle \rf::refl_avgdens[*,33] 26-Oct-07 handle 2-D square bracket clauses, as in an array subscript 31-Aug-07 made xAxis work when SigMath TDI used 06-Apr-07 made work for 'maxval(\activespec::chers_best:vt,1)' 03-Jun-05 Added Z data 01-May-01 When signals are TDI combinations of others, use the units and xaxis from the first signal (could be erroneous) 11-Nov-99 test for input signal name not defined 10-Aug-99 Use Arg_Present so don't have to set keywords before call 26-Jan-99 Trim trailing blanks from unit strings [BD] 13-Jan-99 Written by Bill Davis, PPPL
Source: src/idl_cvs/mdslargestnode.pro
src/idl_cvs/mdslargestnode.pro
NAME: mdsLargestNode PURPOSE: find largest node in a tree CATEGORY: MDSplus CALLING SEQUENCE: IDL> size = mdsLargestNode( tree, searchStr, shot=shot, $ biggestSig=biggestSig ) INPUTS: tree - MDSplus tree. Defaults to 'PCS' optional: searchStr - what to use to search for nodes. Defaults to '*' KEYWORD PARAMETERS: Inputs: shot = NSTX shot number (defaults to latest shot) usage - Defaults to 'SIGNAL'. Can be "ANY", "TEXT" or "NUMERIC" verbose - if set, will print many informational messages debug - if set, debug output will be printed Outputs: status - if odd, then success biggestSig - name of biggest signal found OUTPUTS: size - size of contents of node in bytes EXAMPLES: IDL> size =mdsLargestNode( 'WF', shot=136000, biggestSig=biggestSig ) IDL> size=mdsLargestNode( 'operations', 'PCS', shot=136000, $ biggestSig=biggestSig ) IDL> help BIGGESTSIG STRING = '\OPERATIONS::TOP.PCS.MAGCONTROL.CURRENT:IPLAS_PF1B_1' SIZE LONG = 90 MODIFICATION HISTORY: 22-Aug-2013 Written by Bill Davis for Keith Erickson
Source: src/idl_cvs/mdslist.pro
src/idl_cvs/mdslist.pro
NAME: mdslist PURPOSE: Print ASCII file containing data for MDSplus signal(s). The "dimension" of one of the signals (typically time) is in the first column, followed by data values at that time. Output may be in the format suitable for importing into Kaleidograph, Excel, etc. CATEGORY: MDSplus, Output CALLING SEQUENCE: mdslist, MDStree, shot, signames mdslist, MDStree, shot, signames, $ t1=t1, t2=t2, npts=npts, format=format, tabs=tabs, $ headings=headings, SigForX=SigForX, scope=scope INPUTS: MDStree - e.g. 'operations'. Default is 'NSTX' shot - MDSplus shot number, i.e, ID. signames - a single MDSplus signal name, or array of names. (ignored if Keyword "Scope" is present.) OUTPUTS: none returned KEYWORDS: (all Optional) t1 - Starting x-value desired (typically time, in seconds). Defaults to beginning of data for 1st signal. t2 - Ending x-value desired. Defaults to end of data npts - number of points desired. Defaults to all points. nth - just list every nth point format - output format. Defaults to '(100(g15.6, a1))' (a tab or space is written between the columns) tabs - if set, columns will be separated by tabs instead of spaces headings - if set, the signal labels are at the top of each column in row 1. SigForX - The signal from which to extract the output x-values. Especially relevant when signals of different timebase are requested (defaults to that of the first signal). scope - if present, gets the signal names and column titles from this scope file. (if present, the signames input is ignored). THIS SHOULD BE THE FULL PATHNAME OF THE FILE, or it will default to general areas. outFile - name of output file. Default is 'mdslist.txt' status - 1=OK, 0 means shot, signames, or scope file weren't found. units - if set, will append units to headings EXAMPLES: To write the first 100 points of time and ip pairs to file mdsoutput_106138.dat: IDL> mdslist,'wf',106138,'\ip',npts=100 To write 3 columns, of time, Ip, and Itf, from 0.1-0.2 sec to file mdsoutput_106138.dat: IDL> mdslist,'wf',106138,['\ip','\itf'],t1=.1,t2=.2 To make a tab-delimited file of signals from an MDSplus Scope, with column headings, on efit times: IDL> mdslist,'',106138,scope='/u/Old/bdavis/bv.scope', $ SigForX='EFIT', /headings, /tabs To write time and stored energy data at 1 KHz (like all signals in WF tree): IDL> mdslist,'',106138,'\EFIT01::WMHD/1000', SigForX='\WF::IP' LIMITATION: Would be faster, if it just opened trees that were needed. HISTORY: 01-May-2012 start time different than data wasn't working 06-Apr-2011 added nth keyword for Hiro Takahashi 05-Jan-2010 handle missing signal. Print good error messages and warning on web page and in file. 11-Jul-2008 add /XML output option 24-Mar-08 puts tabs between columns if tabs='on' (as well as if tabs=1) 07-Jan-08 Close file before e-mailing, because not all of file begin sent 26-Sep-07 Add message on how to get file from Linux Cluster nodes 09-Aug-07 bug fix with Yunits not defined 29-Mar-07 changed openMDSshot to mdsopen, so SFLIP and such would work 13-Sep-06 Don't worry if time and data arrays different length 15-Feb-05 Handle when time not first dimension in 2-d array 20-Sep-01 converted from mdsmkfile 05-Sep-01 Added multi-column and interpolation options. Written 13-jul-00 by Bill Davis at the request of Charlie Neumeyer
Source: src/idl_cvs/mdsloadscalar.pro
src/idl_cvs/mdsloadscalar.pro
NAME: MDSloadScalar PURPOSE: Load MDSplus scalar with one simple call (after opening shot) CATEGORY: MDSplus CALLING SEQUENCE: MDSloadScalar, sig, data, SIGUNITS = sigUnits, XAXIS = xAxis, $ XUNITS = xUnits, YAXIS = yAxis, YUNITS = yUnits INPUTS: sig - MDSplus Signal name data - 1- or 2-D data array to load KEYWORD PARAMETERS: Keywords: stat=val - if non-zero, return status (1=OK) out OUTPUTS: NONE COMMON BLOCKS: NONE EXAMPLE stat=openMDSshot('engineering', 100700) sig = '\Some_sig' MDSloadScalar, sig, 3.14159 MODIFICATION HISTORY: 02-Nov-99 Written by Bill Davis, PPPL
Source: src/idl_cvs/mdsloadsig.pro
src/idl_cvs/mdsloadsig.pro
NAME: MDSLOADSIG PURPOSE: Load MDSplus signals with one simple call (after opening shot) Optionally add node and/or tag to a shot tree. CATEGORY: MDSplus CALLING SEQUENCE: MDSLOADSIG, sig, data, SIGUNITS=sigUnits, XAXIS=xAxis, $ XUNITS=xUnits, YAXIS=yAxis, YUNITS=yUnits INPUTS: sigName - MDSplus Signal name data - 1- or 2-D data array to load (can be 3-D, but must have all keywords) KEYWORD PARAMETERS: Keywords: SIGUNITS=val - units of signal in XAXIS=val - array of x-axis in XUNITS=val - units of x-axis in YAXIS=val - array of y-axis in YUNITS=val - units of y-axis in ZAXIS=val - array of z-axis in ZUNITS=val - units of z-axis in FullX - if set, store full x-axis, even if dx constant (1-D only) SkipCheck - if set, do not check time base for constancy (assume) AddNode - if set, will try and add the node to the tree TagToAdd - Tag to add shot - needed if adding node or tag tree - needed if adding node or tag JustAdd - if set, just add node and/or tag, but don't load data stat=val - bad if even, OK if odd out quiet - passed to MDSplus routines (default=1) OUTPUTS: NONE COMMON BLOCKS: NONE EXAMPLE mdsopen,'wf', 99999, stat=stat sigName = '\ip' x = findgen(10) y = sin( x/n_elements( x ) * 3 * 2 * !PI ) MDSLOADSIG, sigName, y, SIGUNITS='Y Units', XAXIS=x, XUNITS='X Units' or, to just add a signal and tag: MDSLOADSIG, sigName, TagToAdd=tag, SHOT=shot, TREE=tree, /JUSTADD MODIFICATION HISTORY: 23-Jul-2013 added noclose and forceAdd keywords 05-Mar-2009 check output keyword on dir in mdstcl (needed for Linux) 09-Jul-2008 use dt_nicenumber, instead of nicenumber, for digitizing rates of 4, etc. 03-Jul-2008 Added dt keyword for Stefan Gerhardt, with ~4e-6 dt. 24-Apr-06 change mds$tcl to mdstcl for Unix. 22-Jun-04 allow 3-d, but must have all parameters 27-Jan-03 Just store x0, dx, y0, and dy for 2-D signals, if xaxis input 23-May-01 Make, as a default, signals with constant timebase store with t0 & delta t, rather than storing whole time array. 16-Nov-00 Added the AddNode & TagToAdd keywords 28-Aug-00 Don't Use Arg_Present 07-Jan-99 Written by Bill Davis, PPPL
Source: src/idl_cvs/mdsmkfile.pro
src/idl_cvs/mdsmkfile.pro
NAME: mdsmkfile PURPOSE: Write an ascii file containing time and data for MDSplus signal(s). Time is in the first column, followed by data values at that time. Output may be in the format suitable for importing into Kaleidograph, Excel, etc. CATEGORY: MDSplus, Output, File I/O CALLING SEQUENCE: mdsmkfile, MDStree, shot, signals mdsmkfile, MDStree, shot, signals, filename=filename, $ t1=t1, t2=t2, npts=npts, format=format, tabs=tabs, $ headings=headings, timeSig=timeSig, scope=scope INPUTS: MDStree - e.g. 'operations'. Default is 'NSTX' shot - MDSplus shot number, i.e, ID. signals - a single MDSplus signal name, or array of names. (ignored if Keyword "Scope" is present.) OUTPUTS: none returned file named mdsoutput_nnnnnn.dat is written, where nnnnnn is the shot number, or, if a scope is used, the scope name is in the place of mdsoutput. If Keyword "Filename" is present, that name is used. KEYWORDS: (all Optional) filename - output filename (default is mdsoutput_nnnnnn.dat, where nnnnnn is the shot number, or, if a scope is used, the scope name -- without extension -- is in the place of mdsoutput). DeltaTime - time step for output - Default is that of the first signal specified. If = 0, then do not list time. t1 - Starting time desired, in seconds. Defaults to beginning of data for 1st signal. t2 - Ending time desired. Defaults to end of data npts - number of points desired. Defaults to all points. format - output format. Defaults to '(100(g15.6, a1))' (a tab or space is written between the columns) tabs - if set, columns will be separated by tabs instead of spaces headings - if set, the signal labels are at the top of each column in row 1. timeSig - The signal from which to extract the output timebase. Especially relevant when signals of different timebase are requested (defaults to that of the first signal). If DeltaTime is specified, this keyword is overridden. scope - if present, gets the signal names and column titles from this scope file. (if present, the Signals input is ignored). THIS SHOULD BE THE FULL PATHNAME OF THE FILE, or it will default to general areas. status - 1=OK, 0 means shot, signals, or scope file weren't found. EXAMPLES: To write the first 100 points of time and ip pairs to file mdsoutput_106138.dat: IDL> mdsmkfile,'wf',106138,'\ip',npts=100 To write 3 columns, of time, Ip, and Itf, from 0.1-0.2 sec to file mdsoutput_106138.dat: IDL> mdsmkfile,'wf',106138,['\ip','\itf'],t1=.1,t2=.2 To make a tab-delimited file of signals from an MDSplus Scope, with column headings, on efit times: IDL> mdsmkfile,'',106138,scope='/u/bdavis/bv.scope', $ timeSig='EFIT', /headings, /tabs To write time and stored energy data at 1 KHz (like all signals in WF tree): IDL> mdsmkfile,'',106138,'\EFIT01::WMHD/1000', timeSig='\WF::IP' LIMITATION: HISTORY: 12-Oct-01 Added deltaTime keyword [BD] 05-Sep-01 Added multi-column and interpolation options. Written 13-jul-00 by Bill Davis at the request of Charlie Neumeyer
Source: src/idl_cvs/mdsnodechanges.pro
src/idl_cvs/mdsnodechanges.pro
NAME: MDSnodeChanges PURPOSE: Print times when the last change was made to an MDSplus node. CATEGORY: MDSplus, Tree Status, dates CALLING SEQUENCE: MDSnodeChanges, node, shot1=shot1, shot2=shot2, skip=skip, tree=tree INPUTS: node = MDSplus node name, e.g. \OPERATIONS::TOP.RAWDATA:LM_H908_01 KEYWORDS: shot1 = starting shot number to process shot2 = last shot number to process (else shot1+400) skip = skip factor tree = MDSplus tree (defaults to nstx -- will be faster if specific tree specified) BEFORE - show only nodes written before this date (defaults to '1-jan-1971'), so can see what never got written. AFTER - show only nodes written after this date (defaults to '1-jan-1971'), so can see only nodes that HAVE been written. OUTFILE - a file name to contain a script related all shots for which a valid node was found in MDSplus PREFIX - string to prefix the shot # in the script. DEFAULT is '/bin/rm -R -f ' MAX - if set, will print max of signal VERBOSE - if set, will tell you what shots could not be opened. EXAMPLE: IDL> tree='microwave' IDL> node='\top.electron_den:line_density' IDL> MDSnodeChanges, node, shot1=101300, shot2=101330, tree=tree IDL> MDSnodeChanges, '\cameras2::FastSoftXray:frames',shot1=124459, $ shot2=125000, /after to produce a script to delete shot directories that have been put in MDSplus: IDL> MDSnodeChanges, node, shot1=shot1, shot2=shot2, /after, outfile='shots.txt' COMMON BLOCKS: NOTES: see MDStreeChanges for who last changed a tree. Note that "OWNER" returns 0 from getnci call. This runs quickly. LIMITATIONS: MODIFICATION HISTORY: 20-Aug-2009 added Max keyword 30-Apr-2008 added outfile and prefix keywords, so directories can be deleted after verifying they are in MDSplus 09-Jul-2007 added AFTER keyword, so you can see what nodes have been written 30-Apr-2007 add BEFORE keyword, so you can see what nodes never got written 26-Apr-2007 "Who" stopped working (returning numberical value) 09-May-2006 Allow for numerical or '*' for who. 20-Dec-99 Written for Rajesh Maingi.
Source: src/idl_cvs/mdsscalars.pro
src/idl_cvs/mdsscalars.pro
NAME: mdsscalars PURPOSE: Print ASCII file containing values of MDSplus scalars(s) for a list of shots. CATEGORY: MDSplus, Output CALLING SEQUENCE: mdsscalars, MDStree, shot, signames INPUTS: MDStree - e.g. 'operations'. Default is 'NSTX' shot - MDSplus shot number, i.e, ID. Can be an array, or range of shots. signames - a single MDSplus signal name, or array of names. OUTPUTS: none returned KEYWORDS: (all Optional) format - output format. Defaults to '(100(g15.6, a1))' (a tab or space is written between the columns) tabs - if set, columns will be separated by tabs instead of spaces headings - if set, the signal labels are at the top of each column in row 1. status - 1=OK, 0 means shot, signames, or scope file weren't found. EXAMPLES: IDL> mdsscalars,'',118956,'\ENGINEERING::TOP.EPICS.GAS.PARAMETERS:INJ1_READY' LIMITATION: Would be faster, if it just opened trees that were needed. HISTORY: Written 20-Aug-2006 by Bill Davis
Source: src/idl_cvs/mdssetevent.pro
src/idl_cvs/mdssetevent.pro
NAME: MDSSETEVENT PURPOSE: Generates an MDSplus event CATEGORY: MDSplus, Events CALLING SEQUENCE: MDSSETEVENT,event[,/QUIET][,STATUS=STATUS][,DATA=data][,/INFO] INPUT PARAMETERS: event = Name of MDSplus event to generate. Keywords: QUIET = prevents IDL error if TCL command fails STATUS = return status, low bit set = success DATA = bytarr(12) to send with event INFO - if this keyword is set the VMS error will be displayed but control will be returned to the caller instead of doing a longjump back to the IDL prompt. OUTPUTS: istat = return status, low bit set = success COMMON BLOCKS: None. SIDE EFFECTS: None. RESTRICTIONS: None. PROCEDURE: Straightforward. Makes a call to MDSplus shared image library procedure MDS$OPEN and checks for errors. MODIFICATION HISTORY: renamed from mds$setevent & made to work on UNIX, too [BD] VERSION 1.0, CREATED BY T.W. Fredian, April 22,1991
Source: src/idl_cvs/mdsshotdate.pro
src/idl_cvs/mdsshotdate.pro
NAME: MDSshotDate PURPOSE: Print time and date of an NSTX shot. CATEGORY: Dates, MDSplus, Tree Status CALLING SEQUENCE: IDL> date = MDSshotDate( shot, time=time ) INPUTS: shot = MDSplus shot number KEYWORDS: optional output: time = time of shot, e.g., 12:25:02 longDate - returned date in form of yyyymmdd, e.g., 20030613 optonally returned: year month day OUTPUT: date - e.g., '11-FEB-2003 15:09:37.60' EXAMPLE: IDL> date=mdsshotdate(110109,time=time) IDL> print, date, ' ', time 11-FEB-2003 15:09:37.60 COMMON BLOCKS: none NOTES: This runs fairly quickly, but SQL access is quicker. Also see SURVEYDB to quickly see first & last shot of day. LIMITATIONS: MODIFICATION HISTORY: 22-Apr=2010 added noOpen keyword 06-Apr-2009 fixed bug with new date format coming from MDSplus & introduced slash keyword, to get date in the form 03-Apr-2009 allow for missing value 14-Jan-2009 Changed node read to be \NSTX::TOP.ACQ_INFO.STATISTICS:STORE_START 04-Sep-2008 replace openMDSshot with mdsopen 09-Feb-2004 Added status (shots before around 104611 don't work) 03-Jun-02 use \WF::ACQ_START 03-Oct-00 Written by Bill Davis
Source: src/idl_cvs/mdssiginterp.pro
src/idl_cvs/mdssiginterp.pro
NAME: mdsSigInterp PURPOSE: call TDI to interpolate signals, but handle expressions in which multiple signals are included. (assumes all same timebase within a given TDI expression) CATEGORY: MDSplus, Timing CALLING SEQUENCE: IDL> interpData = mdsSigInterp(inSig, outXsig_in, status=status) INPUTS: inSig - character string of signal to interpolate outXsig_in - signal to get timebase for output of inSig OUTPUTS: returns array on desired timebase KEYWORDS: (all Optional) inXsig_in - signal to find timebase for. Default is dim_of(inSig) minTime - if set, just returns time range common to both signals outTimes - time base of returned signal status - 1=OK, 0 means problem EXAMPLE: stat=openMDSshot('nstx', 109070) interpData = mdsSigInterp( '\wf::ip', '\engineering::ip1', outTimes=outTimes, /mintime ) plot,outtimes,interpdata HISTORY: 10-Apr-03 outTimes and minTime keywords added 05-Nov-01 Written by Bill Davis
Source: src/idl_cvs/mdstime2dbtime.pro
src/idl_cvs/mdstime2dbtime.pro
NAME: mdsTime2dbTime PURPOSE: Convert date & time from MDSplus to database format CATEGORY: Dates, times CALLING SEQUENCE: sqltime = mdsTime2dbTime( time ) INPUTS: time - from MDSplus, e.g. '10-MAR-2008 09:50:46.36' KEYWORDS OUTPUT sqltime - in SQL format, e.g., Mar 10 2008 09:50AM EXAMPLE: IDL> print,mdsTime2dbTime( '10-MAR-2008 09:50:46.36') MAR 10 2008 09:50AM MODIFICATION HISTORY: WRITTEN by 10-Mar-2008 Bill Davis
Source: src/idl_cvs/mdstreechanges.pro
src/idl_cvs/mdstreechanges.pro
NAME: MDStreeChanges PURPOSE: Print times when the last change was made to an MDSplus tree, and who did it. CATEGORY: MDSplus, Tree Status, dates CALLING SEQUENCE: MDStreeChanges, tree, shot1=shot1, shot2=shot2, skip=skip INPUTS: tree = MDSplus tree name, e.g. 'Microwave' or 'NSTX' KEYWORDS: shot1 = starting shot number to process shot2 = last shot number to process skip = skip factor (default to no skip since = only list changes for a shot if since this date e.g. - '20-sep-2001' COMMON BLOCKS: none EXAMPLE: IDL> mdstreechanges, 'microwave', shot1=101300, shot2=101301 Latest changes from shot 101300 to 101301 for the microwave tree Shot Who Date Node ---- --- ---- ---- 101300 [BDAVIS] 4-APR-2000 .ELECTRON_DEN:LINE_DENS_DF 101301 [BDAVIS] 4-APR-2000 .ELECTRON_DEN:LINE_DENS_DF MODIFICATION HISTORY: 01-Feb-02 added text for repeated lines 28-Jan-02 added SINCE keyword [BD] 20-Dec-99 Written by Bill Davis for Rajesh Maingi.
Source: src/idl_cvs/mdsunits.pro
src/idl_cvs/mdsunits.pro
NAME: MDSunits PURPOSE: return units of an MDS tag, or the time, or time units CATEGORY: MDSplus CALLING SEQUENCE: units = MDSunits( tag ) INPUTS: tag - MDSplus tag or node name in KEYWORD PARAMETERS Input: TIME - If set, Time returned (unless UNITS also set) UNITS - If set, Units of tag returned (unless TIME also set) Returned: STATUS- Status of call (can be used in IDL logical expressions) OUTPUTS: units - units EXAMPLES: MDSOPEN, 'tftr', 89725 tag = '.waveforms:mb_ip_sl' f = MDSVALUE( tag ) time = MDSUNITS( tag, /TIME ) flabel = MDSUNITS( tag, /UNITS ) timelabel = MDSUNITS( tag, /TIME, /UNITS ) plot, time, f/1000., xtitle = timelabel, ytitle = 'milli'+flabel COMMON BLOCKS: NOTES: MODIFICATION HISTORY: 29-Mar-99 Return status from MDSVALUE [BD] 21-Dec-98 Written by Bill Davis
Source: src/idl_cvs/mdsw.pro
src/idl_cvs/mdsw.pro
NAME: mdsw PURPOSE: Interactively plot MDS signals with Crosshairs and many options CATEGORY: MDSplus, Plotting, Widget Example, Printing Plots, Crosshairs CALLING SEQUENCE: mdsw mdsw, xsize=800, ysize=600 KEWORDS: XSIZE - initial horizontal size of graphics window YSIZE - initial vertical size of graphics window shot - MDSplus shot number, i.e, ID tree - tree to open. Defaults to NSTX COMMON BLOCKS: none EXAMPLE: IDL> mdsw,tree='engineering',shot=130000,charsize=2, Az=140, $ signal='\engineering::ppsum' NOTES: Your display may have to be set to 256 colors to see the crosshairs. LIMITATIONS: Some of the MDSplus signals in the list may not have any data. You need IDL version 5, or greater to run this application. MODIFICATION HISTORY: 16-Mar-2009 Option to not display scalars in Signal list 19-Jan-2008 Added option to change Az & Ax on surface plots 16-Jan-2009 added zAxisW calls for RGA data 17-Oct-05 Add Save to Tek option. 08-Mar-02 Fixed screen dump cabability for 24-bit color, and some minor things. 05-Oct-01 Improved restoring settings. Allowed for "Next Shot" to be missing. 18-May-01 Use Vcomp for "visual compression" if more points than resolution 07-May-01 Added Reload color table. Better behavior when copying plots that don't fill the screen. 08-Mar-01 Bug fixed on Reference Shots with signals of different length. 02-Feb-01 Use Label sub-node as title, if available. Access to nodes, even if not tags. Display text nodes. Display Signal List initially. Plot 2-D data. 30-Nov-00 Allow "Save Settings" to include Signal List & !X, !Y, !P. (old settings files won't work). Axis selection with middle mouse button just redraws last plot. 20-Nov-00 Rearranged menus (added Appearance Menu); better plot spacing. 13-Nov-00 Added "Assume same x-label" & "Assume same x-axis" options 30-Oct-00 Added widget for !P settings, so can specify symbols, etc. 26-Oct-00 Fixed "Copy Plot to Tek" to return to X color table afterwards. 18-Oct-00 Added "Edit Selected Signal" option in Signal Menu. TDI expressions included in name will be displayed. Added button to show treenames in list. 10-Oct-00 Added "Label Every Other tick options" & better tick labels Made postscript output a little better for grid of plots 03-Oct-00 Add "Always Include 0 on Y-axis" (not default) and "Draw dashes at y=0" (now default) options. 28-Sep-00 Allow yaxisw to control Y-axis range, in some situations, at least. 20-Sep-00 Just label "No Data" if data not found. Put < & > on Signal Menu 07-Oct-99 Redraw whole screen when mouse used. Alphabetize signal list. Default to "Only zoom X" 19-Aug-99 Allow overlays of data with different dimensions 23-Jun-99 Added Reference shot. Added "Plot all on one page" option. 12-Apr-99 Added Tek output option [BD] 29-Mar-99 Added X & Y Style editing. [BD] 15-Mar-99 Add "Get Next Good Shot" buttons [BD] 08-Mar-99 Add overlays and legends. [BD] 26-Feb-99 Multiplot option and x-only zooming [BD] 11-Feb-99 Added a pick list for MDS signals 26-Jan-99 Modified for using MDS at PPPL by Bill Davis GA Crosshair routines originally written by John Ferron.
Source: src/idl_cvs/mdsw_noch.pro
src/idl_cvs/mdsw_noch.pro
NAME: mdsw_noch PURPOSE: Plot one MDS signal simply (without Crosshairs) CATEGORY: Example CALLING SEQUENCE: mdsw_noch COMMON BLOCKS: mdsw_noch NOTES: Your display may have to be set to 256 colors to see the crosshairs. LIMITATIONS: The printers and help file used here only work from UNIX. Only one version of this program may be run from a given IDL session. MODIFICATION HISTORY: 11-Feb-99 Using UNSETUP_X rather than SETUP_X 28-Jan-99 Took out cross-hair routines 26-Jan-99 Modified for using MDS at PPPL by Bill Davis 8/29/97 CH_example written by Bill Davis,
Source: src/idl_cvs/medsmooth.pro
src/idl_cvs/medsmooth.pro
NAME: MEDSMOOTH PURPOSE: Median smoothing of a vector, including point near it's ends. CATEGORY: Smoothing, Math CALLING SEQUENCE: SMOOTHED = MEDSMOOTH( VECTOR, WINDOW_WIDTH ) INPUTS: VECTOR = The vector to be smoothed WINDOW = The full width of the window over which the median is determined for each point. OUTPUT: Function returns the smoothed vector SUBROUTINES CALLED: MEDIAN, to find the median PROCEDURE: Each point is replaced by the median of the nearest WINDOW of points. The width of the window shrinks towards the ends of the vector, so that only the first and last points are not filtered. These points are replaced by forecasting from smoothed interior points. REVISION HISTORY: Written, H. Freudenreich, STX, 12/89 H.Freudenreich, 8/90: took care of end-points by shrinking window.
Source: src/idl_cvs/medsmoothedge.pro
src/idl_cvs/medsmoothedge.pro
NAME: medsmoothedge PURPOSE: Median smoothing of an array, including handling edges. CATEGORY: Smoothing CALLING SEQUENCE: SMarray - medsmoothedge( ARRAY, WINDOW_WIDTH ) INPUTS: ARRAY - The array to be smoothed WINDOW_WIDTH - The full width of the window over which the median is determined for each point. OUTPUT: SMarray - smoothed array SUBROUTINES CALLED: MEDIAN, to find the median, and smoothedge PROCEDURE: Each point is replaced by the median of the nearest WINDOW of points. the rows aroung the edge, which are not smoothed by MEDIAN, are replaced by those from calling smoothedge.pro (so are not median-smoothed) REVISION HISTORY: 28-Aug-2012 Written by Bill Davis.
Source: src/idl_cvs/minmax.pro
src/idl_cvs/minmax.pro
NAME: minmax PURPOSE: print the min and the max of an array CATEGORY: programming CALLING SEQUENCE: minmax, a INPUTS: a- an array in KEYWORD PARAMETERS OUTPUTS: none (just printing to IDL session) EXAMPLE: COMMON BLOCKS: NOTES: MODIFICATION HISTORY: 1998 Written by Bill Davis
Source: src/idl_cvs/mkdatedir.pro
src/idl_cvs/mkdatedir.pro
NAME: mkdatedir PURPOSE: Make a directory of the form YYYYMMDD (if it doesnt't exist) and move files there if the files were last modified on that date. CATEGORY: Files, Dates CALLING SEQUENCE: IDL> mkdatedir, path=path, prefix=prefix, ext=ext, date=date INPUTS: none KEYWORD PARAMETERS: PATH - directory PREFIX - file prefix for deletion EXT - file extension for deletion; file search will be path+prefix+'*.'+ext DATE - date for search, of form yyyymmdd, e.g., 20050617 (default=TODAY) ALL - if set, all files will be moved to the right directory (overrides DATE) VERBOSE - if set, will print informational output OUTPUTS: none LIMITATIONS Only works on Unix/Linux MODIFICATION HISTORY: 07-Mar-2007 Added keyword all 22-Aug-2005 Written by Bill Davis, PPPL
Source: src/idl_cvs/mk_bitarray.pro
src/idl_cvs/mk_bitarray.pro
NAME: mk_bitarray PURPOSE: Create an array of 1's & 0's corresponding to input bits set (works for negative numbers, too, unlike similar routines) CATEGORY: Bits CALLING SEQUENCE: IDL> array= mk_bitarray( input ) IDL> array= mk_bitarray( input, nbits=5 ) INPUTS: input = whatever; might be something like !X.STYLE KEYWORD PARAMETERS: NBITS=nbits - length of returned array (default to input type) PRINT - if set, will print bits in groups of fours. OUTPUTS: Byte array containing 1's and 0's out COMMON BLOCKS: NONE EXAMPLE: IDL> print,mk_bitarray(3, nbits=8) 1 1 0 0 0 0 0 0 IDL> dum = mk_bitarray( 1025, /print )' 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 LIMITATION: Only works on a scalar. MODIFICATION HISTORY: 05-Jun-00 default nbits to input type. add print keyword. 30-Mar-99 Written by Bill Davis, PPPL
Source: src/idl_cvs/mk_color.pro
src/idl_cvs/mk_color.pro
NAME: mk_color PURPOSE: Create color tables with 12 fixed colors at the top. These colors can be referenced by name. CATEGORY: Colors, Graphics CALLING SEQUENCE: index = mk_color( color ) index = mk_color( color, TABLE=5, MAX_COLORS_TO_USE = 128 ) OPTIONAL INPUTS: COLOR: A string with the name of the color. 'WHITE','BLACK','YELLOW', 'RED','GREEN','BLUE','MAGENTA','LTBLUE', and (if device supports more than 8 colors) 'GREY' are allowed. Alternately, COLOR may be an integer from 0 - 8. If the COLOR input is anything else, or not present, an index for yellow is returned. KEYWORD PARAMETERS: LOAD: If set, return a structure containing all 8 colors. e.g. IDL> c=MK_COLOR(/load) IDL> plot,x,y,color=c.red MAX_COLORS_TO_USE: The maximum number of colors to use for all of IDL . The default is 64 so IDL doesn't grab all the colors on a 256-color monitor. If an IDL window, or a call to LOADCT, was made before the first call to mk_color, the number of colors will already have been allocated. TABLE: The number of the pre-defined color table to load, from 0 to 40. The same as for LOADCT. The default is 0. FILE: If present, will load color table from this file. These colors will be squeezed into the number of colors used for the palette (N_NONFIXED colors). N_NONFIXED (returned) number of colors other than the nine fixed ones. e.g., MAX_COLORS_TO_USE-12. SILENT: If this keyword is set, the Color Table message is suppressed. SEARCH - if set, will do a least squares fit to color table to find closest index to color -- good if color table was set outside of mk_color GAcolors - if set, create a color table with colors beginning at zero, in the same order as GA color_setup routine. COMMON BLOCKS: mk_color_local SIDE EFFECTS: Limits the maximum number of colors for this IDL session. (The default is 64 colors). RESTRICTIONS: Reserving the right number of colors only works if this is called before any IDL window is drawn, and before LOADCT is called. Good Luck if you try and run with 24-bit color. If a device allows less than 12 colors (like Tek), the colors beyond the number allowed will not be available. EXAMPLES: Simply: plot,indgen(100),color=mk_color('red'),background=mk_color('white') or: plot,[0,1],[0,10]-1 for i=0,8 do oplot,[i,i],color=mk_color(i) ; to see all colors x-y overlays: PLOT,INDGEN(100),/NODATA, COLOR=MK_COLOR('blue') ; axes in blue PLOT,INDGEN(100), COLOR=MK_COLOR('red') ; data in red OPLOT,INDGEN(100)/2, COLOR=MK_COLOR('yellow') ; overlay yellow Using color contours AND x-y plots: array_2D = DIST(150) ; generate a test array yellow_index = mk_color( TABLE=5, MAX_COLORS_TO_USE = 128, $ N_NONFIXED = n_nonfixed) WINDOW, XSIZE=400, YSIZE=250 ; make wide window TV, BYTSCL( array_2D, TOP = n_nonfixed-1 ), 25, 50 ; plot scaled 2-D image PLOT, HISTOGRAM(array_2D), COLOR=mk_color('LtBlue'), $ ; lt. blue axes POSITION=[0.55, 0.15, 0.95, 0.95], /NOERASE, /NODATA OPLOT, HISTOGRAM(array_2D) > 10, COLOR=mk_color('Magenta') ; data in magenta whatever = mk_color(TABLE=39) ; change colors of contours, but not lines NOTES: Tek Colors are 0=black, 1=white, 2=red, 3=green, 4=blue, 5=cyan, 6=magenta, 7=yellow, 8=orange. default Green is [0,192,0] so it won't be so light. MODIFICATION HISTORY: 29-Sep-2009 made dark green actually forest green, so shows up better. make green [0,192,0] rather than [0,255,0], so not so bright. 24-May-04 moved creation of X-window when Z-buffer the device (necessary when mk_color first called with PS or Z.) 22-Apr-2004 fixed bug when Z-buffer set on first call. 14-Nov-03 - fixed new bug for Tek colors. 13-Mar-03 - if requested color not found, try to find black (instead of yellow) if ask for 'FOREGROUND' or 'BACKGROUND' return !p values. GAcolor Keyword for doing colors like GA (same order and position) 20-Aug-02 - set !p.color & !p.background to previous colors (closest r,g,b value), when called the first time, or when changing color tables. Add FILE keyword 07-Mar-02 - Change !d.n_colors to !d.table_size 08-May-01 - When in Tek, make !p.color=black (was going to index 15) 23-Apr-01 - added three more colors (orange, purple, darkgreen) 26-Jan-01 - limit structure returned to !d.n_colors (e.g., when=2). 28-Sep-00 - Added bottom keyword (passes straight through to loadct) 24-Jan-00 - If color found at multiple indices, return highest valid one 02-Jan-00 - Both /LOAD and individual colors return valid indices in 24-bit mode. 05-Oct-99 - BD Synonyms & Structure return as in new D. Fanning routine. 09-Jun-99 - BD If desired color not found, return index of nearest color. 11-May-99 - BD make return valid values for 24-bit (decomposed) color 28-Jan-99 - BD handle device with less than 9 colors. 14-Dec-98 - BD allow COLOR input to be an integer from 0-8. Return index of 0 when color not found in table. 04-Dec-98 - Bill Davis changed GETCOLORS to load color table, return color index, etc. Colors were added. GETCOLORS Written by: David Fanning, 10 February 96.
Source: src/idl_cvs/mk_filename.pro
src/idl_cvs/mk_filename.pro
NAME: mk_filename PURPOSE: File names with system independent symbolic directories. CATEGORY: Files, OS Specific, File I/O CALLING SEQUENCE: f = mk_filename(symdir, name) INPUTS: symdir = symbolic directory name. in name = file name. in KEYWORD PARAMETERS: Keywords: /NOSYM means directory given is not a symbolic name. OPSYS=os specifiy operating system to over-ride actual. Use VMS, WINDOWS, or MACOS. UNIX is default. SUBDIR=s Subdirectory below the symdir, e.g., SUBDIR='dir1/dir2' DELIM=c delimiter character between directory and file name. This is returned. OUTPUTS: f = file name including directory. out COMMON BLOCKS: NOTES: Notes: symdir is a logical name for VMS and an environmental variable for UNIX and WINDOWS. Ex: DEFINE IDL_IDLUSR d0:[publib.idl] for VMS set IDL_IDLUSR=c:\IDL\LIB\IDLUSR for WINDOWS. setenv IDL_IDLUSR /usr/pub/idl for UNIX. Then in IDL: f=mk_filename('IDL_IDLUSR','tmp.tmp') will be the name of the file tmp.tmp in IDL_IDLUSR. MODIFICATION HISTORY: R. Sterner, 4 Feb, 1991 R. Sterner, 27 Mar, 1991 --- added /NOSYM R. Sterner, 21 May, 1991 --- If not a listed opsys assume unix. R. Sterner, 4 Jun, 1991 --- added DOS. R. Sterner, 2 Jul, 1991 --- added DELIM. R. Sterner, 17 Jan, 1992 --- added OPSYS= and avoided double // R. Sterner, 1994 Feb 7 --- Added MacOS. R. Sterner, 1994 Feb 14 --- Changed DOS to windows. B. Davis - nenamed from filename.pro. Added SUBDIR keyword Copyright (C) 1991, Johns Hopkins University/Applied Physics Laboratory This software may be used, copied, or redistributed as long as it is not sold and this copyright notice is reproduced on each copy made. This routine is provided as is without any express or implied warranties whatsoever. Other limitations apply as described in the file disclaimer.txt.
Source: src/idl_cvs/mk_html_help.pro
src/idl_cvs/mk_html_help.pro
NAME: MK_HTML_HELP PURPOSE: Given a list of IDL procedure files (.PRO), VMS text library files (.TLB), or directories that contain such files, this procedure generates a file in the HTML format that contains the documentation for those routines that contain a DOC_LIBRARY style documentation template. The output file is compatible with World Wide Web browsers. CATEGORY: Help, documentation. CALLING SEQUENCE: MK_HTML_HELP, Sources, Outfile INPUTS: Sources: A string or string array containing the name(s) of the .pro or .tlb files (or the names of directories containing such files) for which help is desired. If a source file is a VMS text library, it must include the .TLB file extension. If a source file is an IDL procedure, it must include the .PRO file extension. All other source files are assumed to be directories. Outfile: The name of the output file which will be generated. KEYWORDS: BY_CATEGORY: if set, will group output by their first Category INCLUDE: set to 'CATEGORY' or 'PURPOSE' to include at top of the file TITLE: If present, a string which supplies the name that should appear as the Document Title for the help. VERBOSE: Normally, MK_HTML_HELP does its work silently. Setting this keyword to a non-zero value causes the procedure to issue informational messages that indicate what it is currently doing. !QUIET must be 0 for these messages to appear. STRICT: If this keyword is set to a non-zero value, MK_HTML_HELP will adhere strictly to the HTML format by scanning the the document headers for characters that are reserved in HTML (<,>,&,"). These are then converted to the appropriate HTML syntax in the output file. By default, this keyword is set to zero (to allow for faster processing). EXTENSION: if want to do for files other than .pro files STARTHEAD: characters at the beginning of a line which indicate the header start (default = ';+') STOPHEAD: characters at the beginning of a line which indicate the header end (default = ';-') caveat: string that will be displayed above the links, like a copyright notice. EXAMPLE: if on Unix: % cd directory-where-web-pages-will-be (e.g., /w/nstx.pppl.gov/htdocs/nstx/Software/Programming) % rm pro_sources.txt % /bin/ls -1 the-directory-with-IDL-code/*.pro > pro_sources.txt (e.g., ./src/idl_cvs) % idl IDL> sources = READ_LIST('pro_sources.txt') IDL> MK_HTML_HELP, sources, 'idl_routines.html', include='purpose',/by_category IDL> exit This just gives a relative path, but the web link will find it. When /BY_CATEGORY is specified, an additional page, e.g., idl_routines_alphabet.html, will be produced as well. COMMON BLOCKS: None. SIDE EFFECTS: A help file with the name given by the Outfile argument is created. RESTRICTIONS: The following rules must be followed in formatting the .pro files that are to be searched. (a) The first line of the documentation block contains only the characters ";+", starting in column 1. (b) There must be a line which contains the string "NAME:", which is immediately followed (same or next line) by the name of the procedure or function being described in that documentation block. If this NAME field is not present, the name of the source file will be used. (c) Likewise, for organizing by category, after the purpose, The "CATEGORY:" line should follow. If fitting, use one of existing categories as the first category (the only one sorted on). If less than two routines are in a category, they are listed in the Misc category. (d) The last line of the documentation block contains the characters ";-", starting in column 1. (e) Every other line in the documentation block contains a ";" in column 1. Note that a single .pro file can contain multiple procedures and/or functions, each with their own documentation blocks. If it is desired to have "invisible" routines in a file, i.e. routines which are only for internal use and should not appear in the help file, simply leave out the ";+" and ";-" lines in the documentation block for those routines. MODIFICATION HISTORY: 04-May-01 Add /BY_CATEGORY switch [BD] 31-Oct-00 Add link to actual code [BD] 07-Jun-99 Added Extension keyword [BD] 01-Apr-99 Allow CATEGORY or PURPOSE to be included in the top of of the HTML file [Bill Davis] July 17, 1995, DD, RSI. Added code to alphabetize the subjects; At the end of each description block in the HTML file, added a reference to the source .pro file. July 13, 1995, Mark Rivers, University of Chicago. Added support for multiple source directories and multiple documentation headers per .pro file. July 5, 1995, DD, RSI. Original version.
Source: src/idl_cvs/mk_jpeg.pro
src/idl_cvs/mk_jpeg.pro
NAME: MK_JPEG PURPOSE: make a bit-map color jpeg file of an IDL window for (viewable) insertion into MS Word, or Powerpoint documents. CATEGORY: Graphics, Publication CALLING SEQUENCE: IDL> mk_jpeg, FILENAME=filename, WINDOW=window, /VERBOSE INPUTS: None required KEYWORD PARAMETERS: Optional Inputs: filename - Name of output file; Default='idljpeg.jpg' Window - window number to dump to file VERBOSE - If set, print debugging information OUTPUTS: Just the file out COMMON BLOCKS: NONE EXAMPLE: IDL> LoadCT, 5 IDL> TH_IMAGE_CONT, dist(400) ; make a color contour with color bar IDL> mk_jpeg, filename='myjpeg.jpg' NOTES: do not call mk_color, or other things that limit the number of colors, before running this, so you get 256 colors. You must have priviledges to write the file. MODIFICATION HISTORY: 28-Jul-2010 add status to detect stale NFS handles 09-Aug-2009 Allow filename to be first argument 30-Aug-2004 Added test for truecolor (24-bit color), which seems to be needed for Windows 22-Feb-00 Written by Bill Davis, PPPL; with thanks to Dave Fanning
Source: src/idl_cvs/mk_mpeg.pro
src/idl_cvs/mk_mpeg.pro
NAME: mk_mpeg PURPOSE: Write a sequence of images from a 3-D array, or a series of tiff files, as an mpeg movie. CATEGORY: animation CALLING SEQUENCE: mk_mpeg, 'movie.mpg' ,ims or mk_mpeg, 'movie.mpg', files=files INPUTS: ims: sequence of images as a 3D array with dimensions [sx, sy, nims] where sx = xsize of images sy = ysize of images nims = number of images OPTIONAL INPUTS: KEYWORD PARAMETERS: delafter: if set delete temporary array after movie was created you should actually always do it otherwise you get problems with permissions on multiuser machines (since /tmp normally has the sticky bit set) rep: if given means repeat every image 'rep' times (as a workaround to modify replay speed). i.e., if = 2, make 2 copies of each frame. files: file list. If just one value, needs to include a wild card justone: just plot D-alpha despeckle : if set call despeckle routine (slow, but less "intrusive" than median) OUTPUTS: None OPTIONAL OUTPUTS: COMMON BLOCKS: SIDE EFFECTS: creates some files in TMPDIR which are only removed when the delafter keyword is used RESTRICTIONS: depends on the program mpeg_encode from University of California, Berkeley, which must be installed in /usr/local/bin PROCEDURE: writes a parameter file based on the dimensions of the image array + the sequence of images in ppm format into a temporary directory; finally spawns mpeg_encode to build the movie EXAMPLE: IDL> shot = 111840 IDL> f = file_search( GETENV("NSTXUSR")+'/divcam/'+strtrim(shot,2)+'/*.tif') IDL> mk_mpeg, strtrim(shot,2)+'_divcam'+'.mpg',files=f,shot=shot,/fliphoriz LIMITATIONS: mpeg_encode must be in your path (not currently on PPPL Linux) MODIFICATION HISTORY: 22-Apr-04 lots of kludges to make nice movies of camera data from Hiroshima University on NSTX. [BD] 29-Apr-02 Modified write_mpeg to accept a filelist or wildcard spec [BD]. Mon Nov 18 13:13:53 1996, Christian Soeller grabbed original from the net and made slight modifications
Source: src/idl_cvs/mk_pdmenu.pro
src/idl_cvs/mk_pdmenu.pro
NAME: mk_pdmenu PURPOSE: make a structure for a pull down menu for CW_PDMENU routine CATEGORY: GUI, Programming CALLING SEQUENCE: menu = mk_pdmenu( list ) INPUTS: list - string array in KEYWORD PARAMETERS: maxPerSubMenu = max # of menu items per sub-menu. If not specified, will try to limit max number to fit on screen. MenuTitle = menu title (Default=Select) OUTPUTS: menu - structure which can be input to CW_PDMenu EXAMPLE: pro test_pdmenu base = WIDGET_BASE() colID=WIDGET_BASE(base,/col) Printer_List = list_printer_unix( ) Printer_List = Printer_List( sort( Printer_List )) printer_pdmenu = MK_PDMENU( Printer_List, maxPerSubMenu=20 ) menu = CW_PDMENU(colID, printer_pdmenu, /RETURN_FULL_NAME, /MBAR) butID = WIDGET_BUTTON(colID, VALUE='Quit', UVALUE='Quit') WIDGET_CONTROL, /REALIZE, base repeat begin ev = WIDGET_EVENT(base) PRINT, ev.value end until ev.value eq 'Quit' WIDGET_CONTROL, /DESTROY, base stop end COMMON BLOCKS: NOTES: MODIFICATION HISTORY: 11-Apr-03 make default maxPerSubMenu fit on screen, if not given 27-Jan-03 Automatically make submenus when list of items is large 26-Mar-99 by Bill Davis
Source: src/idl_cvs/mk_png.pro
src/idl_cvs/mk_png.pro
NAME: MK_png PURPOSE: make a bit-map color png file of an IDL window for (viewable) insertion into MS Word, or Powerpoint documents. CATEGORY: Graphics, Publication CALLING SEQUENCE: IDL> mk_png, FILENAME=filename, WINDOW=window, /VERBOSE INPUTS: None required KEYWORD PARAMETERS: Optional Inputs: filename - Name of output file; Default='idlpng.jpg' Window - window number to dump to file VERBOSE - If set, print debugging information OUTPUTS: Just the file out COMMON BLOCKS: NONE EXAMPLE: IDL> LoadCT, 5 IDL> TH_IMAGE_CONT, dist(400) ; make a color contour with color bar IDL> mk_png, filename='mypng.jpg' NOTES: do not call mk_color, or other things that limit the number of colors, before running this, so you get 256 colors. You must have priviledges to write the file. MODIFICATION HISTORY: 16-Jun-2011 Written by Bill Davis, PPPL for Charles Skinner
Source: src/idl_cvs/mk_scope_gif.pro
src/idl_cvs/mk_scope_gif.pro
NAME: mk_scope_gif PURPOSE: make gif files from scope-like output, using a scope file as input CATEGORY: MDSplus, Plot files KEYWORD INPUTS: shots - a string indicating shots, e.g., 107694 108305 108330-108350 or 108100+20 (this returns shots 108100, 108101,...,108120) scopeFile - if not specified will prompt rightSpace - fractioanal space to right of last plot (0-1) leftSpace - fractioanal space to left of first plot (0-1) colSpace - fractioanal space between plot columns (0-1) vertSpace - fractioanal space between plot rows (0-1) t1 - start time of plots (sec) t2 - end time of plots (sec) nRows - # of rows of plots nCols - # of columns of plots pixmap - if set, will not display plots on screen (much faster) EXAMPLE: (on Unix:) mk_scope_gif,shots='109071-109075',scope='wf.scope', /pixmap or mk_scope_gif,shots='109071+5' ; (will be prompted for file) or mk_scope_gif, shots='109070+2',scope='/u/kaye/Scope/sk.scope', $ leftspace=.07, rightspace=0.005, $ tweenspace=0.008, vertSpace = 0.03, $ nrows=8, ncols=4, charsize=1.3 HISTORY: Written by Bill Davis for Stan Kaye, Dec. 2004
Source: src/idl_cvs/mk_shotlist.pro
src/idl_cvs/mk_shotlist.pro
NAME: mk_shotlist PURPOSE: Parse an input string for shots into an array of shots CATEGORY: Input, MDSplus CALLING SEQUENCE: IDL> shotlist = mk_shotlist(input) INPUTS: input =e.g., 107694 108305 108330-108350 (108100+20 will search 108100-108120) if a shot number is given as 0, use latest MDSplus shot KEYWORD PARAMETERS: dupsOK - if=0, will remove duplicate shot numbers OUTPUTS: opt = returned value out COMMON BLOCKS: NONE EXAMPLES: IDL> print,mk_shotlist('112300+3') 112300 112301 112302 112303 IDL> print,mk_shotlist('112300-4') 112300 112299 112298 112297 112296 IDL> print,mk_shotlist('0+3') 113862 113863 113864 113865 IDL> print,mk_shotlist('0-3') 113862 113861 113860 113859 IDL> print,mk_shotlist('112300-112302 112305+2') 112300 112301 112302 112305 112306 112307 IDL> print,mk_shotlist('107694 108305 108330-108332') 107694 108305 108330 108331 108332 IDL> print,mk_shotlist('4(138846) 4(138847)') IDL> print,mk_shotlist('138846-2 4(138847)') NOTES: When there are duplicate shots, the list is ordered, and the duplicates removed. MODIFICATION HISTORY: 22-Jul-2013 allow inputs with multipliers, like 4(138846) *** made the default dupsOK=1 23-Nov-2009 fixed bug with test for negative shot number 08-Aug-07 added dupsOK keyword 21-Aug-06 Allow for blanks around "-" and "+" 14-Jul-04 Written by Bill Davis, PPPL
Source: src/idl_cvs/monevents.pro
src/idl_cvs/monevents.pro
NAME: monevents PURPOSE: Widget to monitor occurrence of MDSplus events & optional Shot #, if passed with the event declaration. The time of the event can also be displayed. EXPLANATION: Lists various events related to an NSTX shot cycle (as a default). As the events are declared, the boxes turn color. When the first event in the list is recognized, the other boxes are set to yellow (unless /noClear was set). Shot numbers declared with the events may optionally be displayed. CATEGORY: Events, MDSplus, Monitoring Widget CALLING SEQUENCE: IDL> monevents to get a small # of key events: IDL> monevents, EVENTS_WANTED=['NSTX_SOC', $ 'NSTX_SOD', 'PC_908C16N08_ACQ', $ ; Digitizer that EFIT needs 'MPTS_FIT_DONE', 'NSTX_ACQ_DONE' , $ 'phoenixDone_p1' ], $ /SHOT, /NoClear to get Diagnostic events and labels: IDL> monevents, file_input='/u/bdavis/cvs/idl_cvs/eventlabelsdiags.txt', $ tabs=0, maxper=12, /shot, /noclear, /WebGifs to get all MPTS-relevant events: IDL> monevents,/shot,add=['NSTX_ABC', 'MPTS_FORCEFIT', $ 'NSTX_ABS','NSTX_KLC','TS_RAW_READY'] to see all MEMS at once: IDL> monevents,/MEMS,tabs=0, maxper=22, /shot, clear_event='NSTX_SOC' ,/noclear for MEMS testing: IDL> monevents,/shot,add=['sflip_files','test2','qcs_test'],tabs=0,/nocle, $ logfile='$HOME/ShotCycleEvents.log' for watching TF timing: IDL> monevents,/shot,tabs=0,/nocle, $ EVENTS_WANTED=['NSTX_SOD','PC_908C36N02_ACQ','PC_908C36N06_ACQ', $ 'PC_908C36N10_ACQ','PC_908C36N14_ACQ','TFMON_ACQ_DONE']; for monitoring TF joint events, with a logfile in your home directory: for watching EFIT timing: IDL> monevents,/shot,/noclear,tabs=0, $ EVENTS_WANTED=['NSTX_SOD','FDIA_CALC_DONE', $ 'OP_H908_01_ACQUIRED', 'OP_H908_02_ACQUIRED', 'OP_H908_04_ACQUIRED', $ 'OP_H908_05_ACQUIRED', 'OP_H908_07_ACQUIRED', 'OP_H908_08_ACQUIRED', $ 'OP_H908_09_ACQUIRED', 'OP_H908_11_ACQUIRED', 'OP_H908_13_ACQUIRED', $ 'PHOENIXGO_P1', 'MPTS_FIT_DONE', 'PHOENIXGO_P2V1', $ 'phoenixDone_p1', $ 'NSTX_ACQ_DONE', $ 'phoenixDone_p2' ], $ logfile='$HOME/EfitEvents.log' INPUTS: none. KEYWORD PARAMETERS: (optional) EVENTS_WANTED - String array of events to watch for LABELS - strings array to display instead of events FILE_INPUT - a text file of events to watch for (labels can follow each, in quotes) format like: NSTX_SOC "Start of Cycle" BA_L8212_01_ACQ "Bolometers" GS_908C28N14_ACQ "Gas Injection" ... CLEAR_EVENT - event to clear green lights on (default to EVENTS_WANTED[0] ADD_EVENTS - these events will be added to default events. MEMS - if set, will read the file for MEMS events, and use all there. LogFile - will log events to this file. LogEvents - if set, and LogFile not set, will log events to MonEvens.log nCols - number of columns wanted. Will default to 15 per column maxPerCol - another way to determine # of columns. Default=15. tabs - will default to tabs if more than one column. if=0, all on 1. NOCLEAR - if set, will not clear text boxes when first event is recognized PRINT - if set, will print messages when events are recognized SHOWSHOTS - if set, will display shot numbers declared with the events TIME - if set, will display time-of-day when event was declared BIG - if set, fonts are big WALL - if set, events of interest for display wall, and big font GIFseconds - seconds between writes of GIF file of window contents (good for making web-accessible displays). If not present, don't make files WebGifs - if set, will create a gif (to be read from the web) whenever widget is updated (DEFAULT=0) BD - display events favored by a certain programmer. GROUP_LEADER - Group_Leader ID VERBOSE - if set, lots of information is output OUTPUTS: none. TESTING IDL> monevents, EVENTS_WANTED=['MY_TEST_EVENT_1','MY_TEST_event_2'], $ /Shots and from a VMS computer: IDL> SetMDSShotEvent,'MY_TEST_EVENT_1', 999999 COMMON BLOCKS: none ROUTINES USED: MK_COLOR, usingXwindows, MDSEVENT EXAMPLE: IDL> monevents IDL> monevents, num=450, prefix='gt', time=0, tab=0, npercol=30, show=0, shot=0 LIMITATION: Shot numbers are only returned when running on VMS currently. NOTES: The time of day will be printed for the first event, and the other times will be cleared. Then the time difference will be printed as mm:ss, so :12 would mean it came 12 seconds after the first event. The status fields are initially the background color. A real program might want to save !D.WINDOW before the TV commands and restore it immediately afterwards. MODIFICATION HISTORY: 12-Apr-2010 add num and prefix keywords so can generate event names 01-Oct-2008 Removed mdsconnect [BD] 18-Aug-2008 Added labels to be displayed instead of events, if desired. added WebGifs & GifSeconds keywords. Add 10-Jul-2008 made logging of events clearer 21-Feb-07 added /logevents and logfile as keywords 08-Feb-06 Added clear_event keyword, and made it show up 1st in list 09-Nov-05 added tabs and QCS keyword 03-Feb-04 added wall keyword 10-Feb-03 added keyword add_event [BD] 14-Jun-01 Default to showing time 20-Jul-00 Written by Bill Davis, PPPL
Source: src/idl_cvs/monevents_noxmanager.pro
src/idl_cvs/monevents_noxmanager.pro
NAME: monevents version to run in batch mode. PURPOSE: Widget to monitor occurrence of MDSplus events & optional Shot #, if passed with the event declaration. The time of the event can also be displayed. EXPLANATION: Lists various events related to an NSTX shot cycle (as a default). As the events are declared, the boxes turn color. When the first event in the list is recognized, the other boxes are set to yellow (unless /noClear was set). Shot numbers declared with the events may optionally be displayed. CATEGORY: Events, MDSplus, Monitoring Widget CALLING SEQUENCE: IDL> monevents to get a small # of key events: IDL> monevents, EVENTS_WANTED=['NSTX_SOC', $ 'NSTX_TMINUS60', 'NSTX_SOD', $ 'NSTX_ACQ_STARTED' ,'PC_908C16N08_ACQ', $ ; Digitizer that EFIT needs 'MPTS_FIT_DONE', 'NSTX_ACQ_DONE' , $ 'SFLIP_IN_TREE' ], $ /SHOT, /NoClear to get Diagnostic events and labels: IDL> monevents_noxmanager, file_input='/u/bdavis/cvs/idl_cvs/eventlabelsdiags.txt', $ tabs=0, maxper=12, /shot, /noclear, /WebGifs, $ outFile='NSTX_Diagnostic_Status.gif' IDL> monevents_noxmanager, file_input='/u/bdavis/cvs/idl_cvs/eventlabelssmall.txt', $ tabs=0, maxper=12, /shot, /noclear, /WebGifs, $ outFile='NSTX_Diagnostic_Status_Example.gif' to get all MPTS-relevant events: IDL> monevents,/shot,add=['NSTX_ABC', 'MPTS_FORCEFIT', $ 'NSTX_ABS','NSTX_KLC','TS_RAW_READY'] to see all QCS at once: IDL> monevents,/qcs,tabs=0,maxper=22, /shot, clear_event='NSTX_SOC' ,/noclear for QCS testing: IDL> monevents,/shot,add=['sflip_files','test2','qcs_test'],tabs=0,/nocle, $ logfile='$HOME/ShotCycleEvents.log' for watching TF timing: IDL> monevents,/shot,tabs=0,/nocle, $ EVENTS_WANTED=['NSTX_SOD','PC_908C36N02_ACQ','PC_908C36N06_ACQ', $ 'PC_908C36N10_ACQ','PC_908C36N14_ACQ','TFMON_ACQ_DONE']; for monitoring TF joint events, with a logfile in your home directory: for watching EFIT timing: IDL> monevents,/shot,/noclear,tabs=0, $ EVENTS_WANTED=['NSTX_SOD','PC_908C16N08_ACQ', 'FDIA_CALC_DONE', $ 'OP_H908_01_ACQUIRED', 'OP_H908_02_ACQUIRED', 'OP_H908_04_ACQUIRED', $ 'OP_H908_05_ACQUIRED', 'OP_H908_07_ACQUIRED', 'OP_H908_08_ACQUIRED', $ 'OP_H908_09_ACQUIRED', 'OP_H908_11_ACQUIRED', 'OP_H908_13_ACQUIRED', $ 'PHOENIXGO_P1', 'MPTS_FIT_DONE', 'PHOENIXGO_P2V1', $ 'phoenixDone_p1', $ 'NSTX_ACQ_DONE', $ 'phoenixDone_p2' ], $ logfile='$HOME/EfitEvents.log' INPUTS: none. KEYWORD PARAMETERS: (optional) EVENTS_WANTED - String array of events to watch for LABELS - strings array to display instead of events FILE_INPUT - a text file of events to watch for (labels can follow each, in quotes) format like: NSTX_SOC "Start of Cycle" BA_L8212_01_ACQ "Bolometers" GS_908C28N14_ACQ "Gas Injection" ... CLEAR_EVENT - event to clear green lights on (default to EVENTS_WANTED[0] ADD_EVENTS - these events will be added to default events. QCS - if set, will read the file for QCS events, and use all there. LogFile - will log events to this file. LogEvents - if set, and LogFile not set, will log events to MonEvens.log nCols - number of columns wanted. Will default to 15 per column maxPerCol - another way to determine # of columns. Default=15. tabs - will default to tabs if more than one column. if=0, all on 1. NOCLEAR - if set, will not clear text boxes when first event is recognized PRINT - if set, will print messages when events are recognized SHOWSHOTS - if set, will display shot numbers declared with the events TIME - if set, will display time-of-day when event was declared BIG - if set, fonts are big WALL - if set, events of interest for display wall, and big font GIFseconds - seconds between writes of GIF file of window contents (good for making web-accessible displays). If not present, don't make files WebGifs - if set, will create a gif (to be read from the web) whenever widget is updated (DEFAULT=0) GROUP_LEADER - Group_Leader ID VERBOSE - if set, lots of information is output OUTPUTS: none. TESTING IDL> monevents, EVENTS_WANTED=['MY_TEST_EVENT_1','MY_TEST_event_2'], $ /Shots and from a VMS computer: IDL> SetMDSShotEvent,'MY_TEST_EVENT_1', 999999 COMMON BLOCKS: none ROUTINES USED: MK_COLOR, usingXwindows, MDSEVENT EXAMPLE: IDL> monevents LIMITATION: Shot numbers are only returned when running on VMS currently. NOTES: The time of day will be printed for the first event, and the other times will be cleared. Then the time difference will be printed as mm:ss, so :12 would mean it came 12 seconds after the first event. The status fields are initially the background color. A real program might want to save !D.WINDOW before the TV commands and restore it immediately afterwards. MODIFICATION HISTORY: 01-Oct-2008 Removed mdsconnect [BD] 18-Aug-2008 Added labels to be displayed instead of events, if desired. added WebGifs & GifSeconds keywords. Add 10-Jul-2008 made logging of events clearer 21-Feb-07 added /logevents and logfile as keywords 08-Feb-06 Added clear_event keyword, and made it show up 1st in list 09-Nov-05 added tabs and QCS keyword 03-Feb-04 added wall keyword 10-Feb-03 added keyword add_event [BD] 14-Jun-01 Default to showing time 20-Jul-00 Written by Bill Davis, PPPL
Source: src/idl_cvs/monitorjob.pro
src/idl_cvs/monitorjob.pro
NAME: monitorjob PURPOSE: Monitor an IDL job, and restart it if not present CATEGORY: Utility CALLING SEQUENCE: IDL> monitorjob, delay=delay, job=job INPUTS: none required (will be prompted for job name, if not a keyword) KEYWORD PARAMETERS: delay - delay in seconds between checks (default to 60) job - name of IDL batchfile for spawning job to monitor. For the job file, just include the call to the procedure e.g., /u/bdavis/cvs/idl_cvs/batchmdsw.pro OUTPUTS: NONE EXAMPLE: (from the directory with batchmdsw.pro:) IDL> monitorjob, job='batchmdsw' MODIFICATION HISTORY: 09-Jun-2006 Written by Bill Davis, PPPL
Source: src/idl_cvs/monthnames.pro
src/idl_cvs/monthnames.pro
NAME: MONTHNAMES PURPOSE: Returns a string array of month names. CATEGORY: Dates CALLING SEQUENCE: mnam = monthnames() INPUTS: monthNum - if present, just the name for this month returned KEYWORD PARAMETERS: nchars - max number of characters in month names OUTPUTS: mnam = string array of 13 items: out ['Error','January',...'December'] COMMON BLOCKS: NOTES: MODIFICATION HISTORY: 28-Jan-2009 Added nchars keyword 22-Feb-2006 Added monthNum optional input [BD] R. Sterner, 18 Sep, 1989 Copyright (C) 1989, Johns Hopkins University/Applied Physics Laboratory This software may be used, copied, or redistributed as long as it is not sold and this copyright notice is reproduced on each copy made. This routine is provided as is without any express or implied warranties whatsoever. Other limitations apply as described in the file disclaimer.txt.
Source: src/idl_cvs/monthnumber.pro
src/idl_cvs/monthnumber.pro
NAME: monthnumber PURPOSE: Returns a month number given a month name. CATEGORY: Dates CALLING SEQUENCE: monthNum = monthnumber(monthStr) INPUTS: monthStr - e.g., 'Dec', or 'MAR' KEYWORD PARAMETERS: OUTPUTS: monthNum = number of month, 1-12, or -1 if not found out COMMON BLOCKS: NOTES: MODIFICATION HISTORY: 09-Jul-2013 made to work with an array 29-Jan-00 Written by Bill Davis
Source: src/idl_cvs/mpcurvefit.pro
src/idl_cvs/mpcurvefit.pro
NAME: MPCURVEFIT AUTHOR: Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 craigm@lheamail.gsfc.nasa.gov UPDATED VERSIONs can be found on my WEB PAGE: http://cow.physics.wisc.edu/~craigm/idl/idl.html PURPOSE: Perform Levenberg-Marquardt least-squares fit (replaces CURVEFIT) MAJOR TOPICS: Curve and Surface Fitting CALLING SEQUENCE: YFIT = MPCURVEFIT(X, Y, WEIGHTS, P, [SIGMA,] FUNCTION_NAME=FUNC, ITER=iter, ITMAX=itmax, CHISQ=chisq, NFREE=nfree, DOF=dof, NFEV=nfev, COVAR=covar, [/NOCOVAR, ] [/NODERIVATIVE, ] FUNCTARGS=functargs, PARINFO=parinfo, FTOL=ftol, XTOL=xtol, GTOL=gtol, TOL=tol, ITERPROC=iterproc, ITERARGS=iterargs, NPRINT=nprint, QUIET=quiet, ERRMSG=errmsg, STATUS=status) DESCRIPTION: MPCURVEFIT fits a user-supplied model -- in the form of an IDL function -- to a set of user-supplied data. MPCURVEFIT calls MPFIT, the MINPACK-1 least-squares minimizer, to do the main work. Given the data and their uncertainties, MPCURVEFIT finds the best set of model parameters which match the data (in a least-squares sense) and returns them in the parameter P. MPCURVEFIT returns the best fit function. The user must supply the following items: - An array of independent variable values ("X"). - An array of "measured" *dependent* variable values ("Y"). - An array of weighting values ("WEIGHTS"). - The name of an IDL function which computes Y given X ("FUNC"). - Starting guesses for all of the parameters ("P"). There are very few restrictions placed on X, Y or FUNCT. Simply put, FUNCT must map the "X" values into "Y" values given the model parameters. The "X" values may represent any independent variable (not just Cartesian X), and indeed may be multidimensional themselves. For example, in the application of image fitting, X may be a 2xN array of image positions. MPCURVEFIT carefully avoids passing large arrays where possible to improve performance. See below for an example of usage. USER FUNCTION The user must define a function which returns the model value. For applications which use finite-difference derivatives -- the default -- the user function should be declared in the following way: PRO MYFUNCT, X, P, YMOD ; The independent variable is X ; Parameter values are passed in "P" YMOD = ... computed model values at X ... END The returned array YMOD must have the same dimensions and type as the "measured" Y values. User functions may also indicate a fatal error condition using the ERROR_CODE common block variable, as described below under the MPFIT_ERROR common block definition. See the discussion under "ANALYTIC DERIVATIVES" and AUTODERIVATIVE in MPFIT.PRO if you wish to compute the derivatives for yourself. AUTODERIVATIVE is accepted and passed directly to MPFIT. The user function must accept one additional parameter, DP, which contains the derivative of the user function with respect to each parameter at each data point, as described in MPFIT.PRO. CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD The behavior of MPFIT can be modified with respect to each parameter to be fitted. A parameter value can be fixed; simple boundary constraints can be imposed; limitations on the parameter changes can be imposed; properties of the automatic derivative can be modified; and parameters can be tied to one another. These properties are governed by the PARINFO structure, which is passed as a keyword parameter to MPFIT. PARINFO should be an array of structures, one for each parameter. Each parameter is associated with one element of the array, in numerical order. The structure can have the following entries (none are required): .VALUE - the starting parameter value (but see the START_PARAMS parameter for more information). .FIXED - a boolean value, whether the parameter is to be held fixed or not. Fixed parameters are not varied by MPFIT, but are passed on to MYFUNCT for evaluation. .LIMITED - a two-element boolean array. If the first/second element is set, then the parameter is bounded on the lower/upper side. A parameter can be bounded on both sides. Both LIMITED and LIMITS must be given together. .LIMITS - a two-element float or double array. Gives the parameter limits on the lower and upper sides, respectively. Zero, one or two of these values can be set, depending on the values of LIMITED. Both LIMITED and LIMITS must be given together. .PARNAME - a string, giving the name of the parameter. The fitting code of MPFIT does not use this tag in any way. However, the default ITERPROC will print the parameter name if available. .STEP - the step size to be used in calculating the numerical derivatives. If set to zero, then the step size is computed automatically. Ignored when AUTODERIVATIVE=0. This value is superceded by the RELSTEP value. .RELSTEP - the *relative* step size to be used in calculating the numerical derivatives. This number is the fractional size of the step, compared to the parameter value. This value supercedes the STEP setting. If the parameter is zero, then a default step size is chosen. .MPSIDE - the sidedness of the finite difference when computing numerical derivatives. This field can take four values: 0 - one-sided derivative computed automatically 1 - one-sided derivative (f(x+h) - f(x) )/h -1 - one-sided derivative (f(x) - f(x-h))/h 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h) Where H is the STEP parameter described above. The "automatic" one-sided derivative method will chose a direction for the finite difference which does not violate any constraints. The other methods do not perform this check. The two-sided method is in principle more precise, but requires twice as many function evaluations. Default: 0. .MPMAXSTEP - the maximum change to be made in the parameter value. During the fitting process, the parameter will never be changed by more than this value in one iteration. A value of 0 indicates no maximum. Default: 0. .TIED - a string expression which "ties" the parameter to other free or fixed parameters. Any expression involving constants and the parameter array P are permitted. Example: if parameter 2 is always to be twice parameter 1 then use the following: parinfo(2).tied = '2 * P(1)'. Since they are totally constrained, tied parameters are considered to be fixed; no errors are computed for them. [ NOTE: the PARNAME can't be used in expressions. ] .MPPRINT - if set to 1, then the default ITERPROC will print the parameter value. If set to 0, the parameter value will not be printed. This tag can be used to selectively print only a few parameter values out of many. Default: 1 (all parameters printed) Future modifications to the PARINFO structure, if any, will involve adding structure tags beginning with the two letters "MP". Therefore programmers are urged to avoid using tags starting with the same letters; otherwise they are free to include their own fields within the PARINFO structure, and they will be ignored. PARINFO Example: parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $ limits:[0.D,0]}, 5) parinfo(0).fixed = 1 parinfo(4).limited(0) = 1 parinfo(4).limits(0) = 50.D parinfo(*).value = [5.7D, 2.2, 500., 1.5, 2000.] A total of 5 parameters, with starting values of 5.7, 2.2, 500, 1.5, and 2000 are given. The first parameter is fixed at a value of 5.7, and the last parameter is constrained to be above 50. INPUTS: X - Array of independent variable values. Y - Array of "measured" dependent variable values. Y should have the same data type as X. The function FUNCT should map X->Y. WEIGHTS - Array of weights to be used in calculating the chi-squared value. If WEIGHTS is specified then the ERR parameter is ignored. The chi-squared value is computed as follows: CHISQ = TOTAL( (Y-FUNCT(X,P))^2 * ABS(WEIGHTS) ) Here are common values of WEIGHTS: 1D/ERR^2 - Normal weighting (ERR is the measurement error) 1D/Y - Poisson weighting (counting statistics) 1D - Unweighted P - An array of starting values for each of the parameters of the model. The number of parameters should be fewer than the number of measurements. Also, the parameters should have the same data type as the measurements (double is preferred). Upon successful completion the new parameter values are returned in P. If both START_PARAMS and PARINFO are passed, then the starting *value* is taken from START_PARAMS, but the *constraints* are taken from PARINFO. SIGMA - The formal 1-sigma errors in each parameter, computed from the covariance matrix. If a parameter is held fixed, or if it touches a boundary, then the error is reported as zero. If the fit is unweighted (i.e. no errors were given, or the weights were uniformly set to unity), then SIGMA will probably not represent the true parameter uncertainties. *If* you can assume that the true reduced chi-squared value is unity -- meaning that the fit is implicitly assumed to be of good quality -- then the estimated parameter uncertainties can be computed by scaling SIGMA by the measured chi-squared value. DOF = N_ELEMENTS(X) - N_ELEMENTS(P) ; deg of freedom CSIGMA = SIGMA * SQRT(CHISQ / DOF) ; scaled uncertainties RETURNS: Returns the array containing the best-fitting function. KEYWORD PARAMETERS: CHISQ - the value of the summed, squared, weighted residuals for the returned parameter values, i.e. the chi-square value. COVAR - the covariance matrix for the set of parameters returned by MPFIT. The matrix is NxN where N is the number of parameters. The square root of the diagonal elements gives the formal 1-sigma statistical errors on the parameters IF errors were treated "properly" in MYFUNC. Parameter errors are also returned in PERROR. To compute the correlation matrix, PCOR, use this: IDL> PCOR = COV * 0 IDL> FOR i = 0, n-1 DO FOR j = 0, n-1 DO $ PCOR(i,j) = COV(i,j)/sqrt(COV(i,i)*COV(j,j)) If NOCOVAR is set or MPFIT terminated abnormally, then COVAR is set to a scalar with value !VALUES.D_NAN. DOF - number of degrees of freedom, computed as DOF = N_ELEMENTS(DEVIATES) - NFREE Note that this doesn't account for pegged parameters (see NPEGGED). ERRMSG - a string error or warning message is returned. FTOL - a nonnegative input variable. Termination occurs when both the actual and predicted relative reductions in the sum of squares are at most FTOL (and STATUS is accordingly set to 1 or 3). Therefore, FTOL measures the relative error desired in the sum of squares. Default: 1D-10 FUNCTION_NAME - a scalar string containing the name of an IDL procedure to compute the user model values, as described above in the "USER MODEL" section. FUNCTARGS - A structure which contains the parameters to be passed to the user-supplied function specified by FUNCT via the _EXTRA mechanism. This is the way you can pass additional data to your user-supplied function without using common blocks. By default, no extra parameters are passed to the user-supplied function. GTOL - a nonnegative input variable. Termination occurs when the cosine of the angle between fvec and any column of the jacobian is at most GTOL in absolute value (and STATUS is accordingly set to 4). Therefore, GTOL measures the orthogonality desired between the function vector and the columns of the jacobian. Default: 1D-10 ITER - the number of iterations completed. ITERARGS - The keyword arguments to be passed to ITERPROC via the _EXTRA mechanism. This should be a structure, and is similar in operation to FUNCTARGS. Default: no arguments are passed. ITERPROC - The name of a procedure to be called upon each NPRINT iteration of the MPFIT routine. It should be declared in the following way: PRO ITERPROC, FUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $ PARINFO=parinfo, QUIET=quiet, ... ; perform custom iteration update END ITERPROC must either accept all three keyword parameters (FUNCTARGS, PARINFO and QUIET), or at least accept them via the _EXTRA keyword. FUNCT is the user-supplied function to be minimized, P is the current set of model parameters, ITER is the iteration number, and FUNCTARGS are the arguments to be passed to FUNCT. FNORM should be the chi-squared value. QUIET is set when no textual output should be printed. See below for documentation of PARINFO. In implementation, ITERPROC can perform updates to the terminal or graphical user interface, to provide feedback while the fit proceeds. If the fit is to be stopped for any reason, then ITERPROC should set the common block variable ERROR_CODE to negative value (see MPFIT_ERROR common block below). In principle, ITERPROC should probably not modify the parameter values, because it may interfere with the algorithm's stability. In practice it is allowed. Default: an internal routine is used to print the parameter values. ITMAX - The maximum number of iterations to perform. If the number is exceeded, then the STATUS value is set to 5 and MPFIT returns. Default: 200 iterations NFEV - the number of FUNCT function evaluations performed. NFREE - the number of free parameters in the fit. This includes parameters which are not FIXED and not TIED, but it does include parameters which are pegged at LIMITS. NOCOVAR - set this keyword to prevent the calculation of the covariance matrix before returning (see COVAR) NODERIVATIVE - if set, then the user function will not be queried for analytical derivatives, and instead the derivatives will be computed by finite differences (and according to the PARINFO derivative settings; see above for a description). NPRINT - The frequency with which ITERPROC is called. A value of 1 indicates that ITERPROC is called with every iteration, while 2 indicates every other iteration, etc. Note that several Levenberg-Marquardt attempts can be made in a single iteration. Default value: 1 PARINFO - Provides a mechanism for more sophisticated constraints to be placed on parameter values. When PARINFO is not passed, then it is assumed that all parameters are free and unconstrained. Values in PARINFO are never modified during a call to MPFIT. See description above for the structure of PARINFO. Default value: all parameters are free and unconstrained. QUIET - set this keyword when no textual output should be printed by MPFIT STATUS - an integer status code is returned. All values other than zero can represent success. It can have one of the following values: 0 improper input parameters. 1 both actual and predicted relative reductions in the sum of squares are at most FTOL. 2 relative error between two consecutive iterates is at most XTOL 3 conditions for STATUS = 1 and STATUS = 2 both hold. 4 the cosine of the angle between fvec and any column of the jacobian is at most GTOL in absolute value. 5 the maximum number of iterations has been reached 6 FTOL is too small. no further reduction in the sum of squares is possible. 7 XTOL is too small. no further improvement in the approximate solution x is possible. 8 GTOL is too small. fvec is orthogonal to the columns of the jacobian to machine precision. TOL - synonym for FTOL. Use FTOL instead. XTOL - a nonnegative input variable. Termination occurs when the relative error between two consecutive iterates is at most XTOL (and STATUS is accordingly set to 2 or 3). Therefore, XTOL measures the relative error desired in the approximate solution. Default: 1D-10 YERROR - upon return, the root-mean-square variance of the residuals. EXAMPLE: ; First, generate some synthetic data npts = 200 x = dindgen(npts) * 0.1 - 10. ; Independent variable yi = gauss1(x, [2.2D, 1.4, 3000.]) ; "Ideal" Y variable y = yi + randomn(seed, npts) * sqrt(1000. + yi); Measured, w/ noise sy = sqrt(1000.D + y) ; Poisson errors ; Now fit a Gaussian to see how well we can recover p0 = [1.D, 1., 1000.] ; Initial guess yfit = mpcurvefit(x, y, 1/sy^2, p0, $ ; Fit a function FUNCTION_NAME='GAUSS1P',/autoderivative) print, p Generates a synthetic data set with a Gaussian peak, and Poisson statistical uncertainty. Then the same function is fitted to the data to see how close we can get. GAUSS1 and GAUSS1P are available from the same web page. COMMON BLOCKS: COMMON MPFIT_ERROR, ERROR_CODE User routines may stop the fitting process at any time by setting an error condition. This condition may be set in either the user's model computation routine (MYFUNCT), or in the iteration procedure (ITERPROC). To stop the fitting, the above common block must be declared, and ERROR_CODE must be set to a negative number. After the user procedure or function returns, MPFIT checks the value of this common block variable and exits immediately if the error condition has been set. By default the value of ERROR_CODE is zero, indicating a successful function/procedure call. REFERENCES: MINPACK-1, Jorge More', available from netlib (www.netlib.org). "Optimization Software Guide," Jorge More' and Stephen Wright, SIAM, *Frontiers in Applied Mathematics*, Number 14. MODIFICATION HISTORY: Translated from MPFITFUN, 25 Sep 1999, CM Alphabetized documented keywords, 02 Oct 1999, CM Added QUERY keyword and query checking of MPFIT, 29 Oct 1999, CM Check to be sure that X and Y are present, 02 Nov 1999, CM Documented SIGMA for unweighted fits, 03 Nov 1999, CM Changed to ERROR_CODE for error condition, 28 Jan 2000, CM Copying permission terms have been liberalized, 26 Mar 2000, CM Propagated improvements from MPFIT, 17 Dec 2000, CM Corrected behavior of NODERIVATIVE, 13 May 2002, CM Documented RELSTEP field of PARINFO (!!), CM, 25 Oct 2002 Make more consistent with comparable IDL routines, 30 Jun 2003, CM Minor documentation adjustment, 03 Feb 2004, CM Fix error in documentation, 26 Aug 2005, CM Convert to IDL 5 array syntax (!), 16 Jul 2006, CM $Id: mpcurvefit.pro,v 2.3 2009/04/15 14:43:47 bdavis Exp $
Source: src/idl_cvs/mpeg_from_screen.pro
src/idl_cvs/mpeg_from_screen.pro
NAME: mpeg_from_screen PURPOSE: Example of creating a color Mpeg animation with text overlayed. Works on both 8-bit and 24-bit monitors. Can read from files, but those returning 24-bit data will not have the text overlaying option. CATEGORY: Animation CALLING SEQUENCE: IDL> mpeg_from_screen or IDL> mpeg_from_screen, data3d=data3D, mpeg_filename='mymovie.mpg' or IDL> mpeg_from_screen, ndups=5, $ files=['out1.jpg','out2.jpg','out3.jpg','out4.jpg','out5.jpg'] INPUTS: none required KEYWORD PARAMETERS: All optional inputs: data3D - optional 3-D array ( nx, ny, nFrames ) nFrames - # of frames when using dummy data. files - a string array of filenames. Can be of type 'jpg','tif','tif','bmp','jpeg,'png', or 'gif' mpeg_filename - output filename for mpeg movie ndups - number of times each frame of movie is duplicated (default=3) showFrame - if =0, will not show frame around image (irrelevant for rgb files) useScreen - if=0, don't bother using screen for output (i.e., no text to overlay) (irrelevant for rgb files) border - # of pixels around data (default=25) ctb - color table to use (default=39) charsize - character size (default=1.5) charthick - character thickness (default=2) pixMap - if = 0, will write and read each frame from screen (much slower) OUTPUTS: mpeg file MODIFICATION HISTORY: 20-Oct-2004 Optionally operate off a list of files 15-Oct-2004 Written by Bill Davis
Source: src/idl_cvs/mpfit.pro
src/idl_cvs/mpfit.pro
NAME: MPFIT AUTHOR: Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 craigm@lheamail.gsfc.nasa.gov UPDATED VERSIONs can be found on my WEB PAGE: http://cow.physics.wisc.edu/~craigm/idl/idl.html PURPOSE: Perform Levenberg-Marquardt least-squares minimization (MINPACK-1) MAJOR TOPICS: Curve and Surface Fitting CALLING SEQUENCE: parms = MPFIT(MYFUNCT, start_parms, FUNCTARGS=fcnargs, NFEV=nfev, MAXITER=maxiter, ERRMSG=errmsg, NPRINT=nprint, QUIET=quiet, FTOL=ftol, XTOL=xtol, GTOL=gtol, NITER=niter, STATUS=status, ITERPROC=iterproc, ITERARGS=iterargs, COVAR=covar, PERROR=perror, BESTNORM=bestnorm, PARINFO=parinfo) DESCRIPTION: MPFIT uses the Levenberg-Marquardt technique to solve the least-squares problem. In its typical use, MPFIT will be used to fit a user-supplied function (the "model") to user-supplied data points (the "data") by adjusting a set of parameters. MPFIT is based upon MINPACK-1 (LMDIF.F) by More' and collaborators. For example, a researcher may think that a set of observed data points is best modelled with a Gaussian curve. A Gaussian curve is parameterized by its mean, standard deviation and normalization. MPFIT will, within certain constraints, find the set of parameters which best fits the data. The fit is "best" in the least-squares sense; that is, the sum of the weighted squared differences between the model and data is minimized. The Levenberg-Marquardt technique is a particular strategy for iteratively searching for the best fit. This particular implementation is drawn from MINPACK-1 (see NETLIB), and seems to be more robust than routines provided with IDL. This version allows upper and lower bounding constraints to be placed on each parameter, or the parameter can be held fixed. The IDL user-supplied function should return an array of weighted deviations between model and data. In a typical scientific problem the residuals should be weighted so that each deviate has a gaussian sigma of 1.0. If X represents values of the independent variable, Y represents a measurement for each value of X, and ERR represents the error in the measurements, then the deviates could be calculated as follows: DEVIATES = (Y - F(X)) / ERR where F is the analytical function representing the model. You are recommended to use the convenience functions MPFITFUN and MPFITEXPR, which are driver functions that calculate the deviates for you. If ERR are the 1-sigma uncertainties in Y, then TOTAL( DEVIATES^2 ) will be the total chi-squared value. MPFIT will minimize the chi-square value. The values of X, Y and ERR are passed through MPFIT to the user-supplied function via the FUNCTARGS keyword. Simple constraints can be placed on parameter values by using the PARINFO keyword to MPFIT. See below for a description of this keyword. MPFIT does not perform more general optimization tasks. See TNMIN instead. MPFIT is customized, based on MINPACK-1, to the least-squares minimization problem. USER FUNCTION The user must define a function which returns the appropriate values as specified above. The function should return the weighted deviations between the model and the data. For applications which use finite-difference derivatives -- the default -- the user function should be declared in the following way: FUNCTION MYFUNCT, p, X=x, Y=y, ERR=err ; Parameter values are passed in "p" model = F(x, p) return, (y-model)/err END See below for applications with analytical derivatives. The keyword parameters X, Y, and ERR in the example above are suggestive but not required. Any parameters can be passed to MYFUNCT by using the FUNCTARGS keyword to MPFIT. Use MPFITFUN and MPFITEXPR if you need ideas on how to do that. The function *must* accept a parameter list, P. In general there are no restrictions on the number of dimensions in X, Y or ERR. However the deviates *must* be returned in a one-dimensional array, and must have the same type (float or double) as the input arrays. User functions may also indicate a fatal error condition using the ERROR_CODE common block variable, as described below under the MPFIT_ERROR common block definition (by setting ERROR_CODE to a number between -15 and -1). ANALYTIC DERIVATIVES In the search for the best-fit solution, MPFIT by default calculates derivatives numerically via a finite difference approximation. The user-supplied function need not calculate the derivatives explicitly. However, if you desire to compute them analytically, then the AUTODERIVATIVE=0 keyword must be passed. As a practical matter, it is often sufficient and even faster to allow MPFIT to calculate the derivatives numerically, and so AUTODERIVATIVE=0 is not necessary. Also, the user function must be declared with one additional parameter, as follows: FUNCTION MYFUNCT, p, dp, X=x, Y=y, ERR=err model = F(x, p) if n_params() GT 1 then begin ; Compute derivatives dp = make_array(n_elements(x), n_elements(p), value=x[0]*0) for i = 0, n_elements(p)-1 do $ dp(*,i) = FGRAD(x, p, i) endif return, (y-model)/err END where FGRAD(x, p, i) is a user function which must compute the derivative of the model with respect to parameter P(i) at X. When finite differencing is used for computing derivatives (ie, when AUTODERIVATIVE=1), the parameter DP is not passed. Therefore functions can use N_PARAMS() to indicate whether they must compute the derivatives or not. Derivatives should be returned in the DP array. DP should be an m x n array, where m is the number of data points and n is the number of parameters. dp[i,j] is the derivative at the ith point with respect to the jth parameter. The derivatives with respect to fixed parameters are ignored; zero is an appropriate value to insert for those derivatives. Upon input to the user function, DP is set to a vector with the same length as P, with a value of 1 for a parameter which is free, and a value of zero for a parameter which is fixed (and hence no derivative needs to be calculated). This input vector may be overwritten as needed. If the data is higher than one dimensional, then the *last* dimension should be the parameter dimension. Example: fitting a 50x50 image, "dp" should be 50x50xNPAR. CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD The behavior of MPFIT can be modified with respect to each parameter to be fitted. A parameter value can be fixed; simple boundary constraints can be imposed; limitations on the parameter changes can be imposed; properties of the automatic derivative can be modified; and parameters can be tied to one another. These properties are governed by the PARINFO structure, which is passed as a keyword parameter to MPFIT. PARINFO should be an array of structures, one for each parameter. Each parameter is associated with one element of the array, in numerical order. The structure can have the following entries (none are required): .VALUE - the starting parameter value (but see the START_PARAMS parameter for more information). .FIXED - a boolean value, whether the parameter is to be held fixed or not. Fixed parameters are not varied by MPFIT, but are passed on to MYFUNCT for evaluation. .LIMITED - a two-element boolean array. If the first/second element is set, then the parameter is bounded on the lower/upper side. A parameter can be bounded on both sides. Both LIMITED and LIMITS must be given together. .LIMITS - a two-element float or double array. Gives the parameter limits on the lower and upper sides, respectively. Zero, one or two of these values can be set, depending on the values of LIMITED. Both LIMITED and LIMITS must be given together. .PARNAME - a string, giving the name of the parameter. The fitting code of MPFIT does not use this tag in any way. However, the default ITERPROC will print the parameter name if available. .STEP - the step size to be used in calculating the numerical derivatives. If set to zero, then the step size is computed automatically. Ignored when AUTODERIVATIVE=0. This value is superceded by the RELSTEP value. .RELSTEP - the *relative* step size to be used in calculating the numerical derivatives. This number is the fractional size of the step, compared to the parameter value. This value supercedes the STEP setting. If the parameter is zero, then a default step size is chosen. .MPSIDE - the sidedness of the finite difference when computing numerical derivatives. This field can take four values: 0 - one-sided derivative computed automatically 1 - one-sided derivative (f(x+h) - f(x) )/h -1 - one-sided derivative (f(x) - f(x-h))/h 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h) Where H is the STEP parameter described above. The "automatic" one-sided derivative method will chose a direction for the finite difference which does not violate any constraints. The other methods do not perform this check. The two-sided method is in principle more precise, but requires twice as many function evaluations. Default: 0. .MPMAXSTEP - the maximum change to be made in the parameter value. During the fitting process, the parameter will never be changed by more than this value in one iteration. A value of 0 indicates no maximum. Default: 0. .TIED - a string expression which "ties" the parameter to other free or fixed parameters as an equality constraint. Any expression involving constants and the parameter array P are permitted. Example: if parameter 2 is always to be twice parameter 1 then use the following: parinfo[2].tied = '2 * P[1]'. Since they are totally constrained, tied parameters are considered to be fixed; no errors are computed for them. [ NOTE: the PARNAME can't be used in expressions. ] .MPPRINT - if set to 1, then the default ITERPROC will print the parameter value. If set to 0, the parameter value will not be printed. This tag can be used to selectively print only a few parameter values out of many. Default: 1 (all parameters printed) .MPFORMAT - IDL format string to print the parameter within ITERPROC. Default: '(G20.6)' An empty string will also use the default. Future modifications to the PARINFO structure, if any, will involve adding structure tags beginning with the two letters "MP". Therefore programmers are urged to avoid using tags starting with the same letters; otherwise they are free to include their own fields within the PARINFO structure, and they will be ignored. PARINFO Example: parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $ limits:[0.D,0]}, 5) parinfo[0].fixed = 1 parinfo[4].limited(0) = 1 parinfo[4].limits(0) = 50.D parinfo[*].value = [5.7D, 2.2, 500., 1.5, 2000.] A total of 5 parameters, with starting values of 5.7, 2.2, 500, 1.5, and 2000 are given. The first parameter is fixed at a value of 5.7, and the last parameter is constrained to be above 50. HARD-TO-COMPUTE FUNCTIONS: "EXTERNAL" EVALUATION The normal mode of operation for MPFIT is for the user to pass a function name, and MPFIT will call the user function multiple times as it iterates toward a solution. Some user functions are particularly hard to compute using the standard model of MPFIT. Usually these are functions that depend on a large amount of external data, and so it is not feasible, or at least highly impractical, to have MPFIT call it. In those cases it may be possible to use the "(EXTERNAL)" evaluation option. In this case the user is responsible for making all function *and derivative* evaluations. The function and Jacobian data are passed in through the EXTERNAL_FVEC and EXTERNAL_FJAC keywords, respectively. The user indicates the selection of this option by specifying a function name (MYFUNCT) of "(EXTERNAL)". No user-function calls are made when EXTERNAL evaluation is being used. At the end of each iteration, control returns to the user, who must reevaluate the function at its new parameter values. Users should check the return value of the STATUS keyword, where a value of 9 indicates the user should supply more data for the next iteration, and re-call MPFIT. The user may refrain from calling MPFIT further; as usual, STATUS will indicate when the solution has converged and no more iterations are required. Because MPFIT must maintain its own data structures between calls, the user must also pass a named variable to the EXTERNAL_STATE keyword. This variable must be maintained by the user, but not changed, throughout the fitting process. When no more iterations are desired, the named variable may be discarded. INPUTS: MYFUNCT - a string variable containing the name of the function to be minimized. The function should return the weighted deviations between the model and the data, as described above. For EXTERNAL evaluation of functions, this parameter should be set to a value of "(EXTERNAL)". START_PARAMS - An array of starting values for each of the parameters of the model. The number of parameters should be fewer than the number of measurements. Also, the parameters should have the same data type as the measurements (double is preferred). This parameter is optional if the PARINFO keyword is used (but see PARINFO). The PARINFO keyword provides a mechanism to fix or constrain individual parameters. If both START_PARAMS and PARINFO are passed, then the starting *value* is taken from START_PARAMS, but the *constraints* are taken from PARINFO. RETURNS: Returns the array of best-fit parameters. KEYWORD PARAMETERS: AUTODERIVATIVE - If this is set, derivatives of the function will be computed automatically via a finite differencing procedure. If not set, then MYFUNCT must provide the (analytical) derivatives. Default: set (=1) NOTE: to supply your own analytical derivatives, explicitly pass AUTODERIVATIVE=0 BESTNORM - the value of the summed squared weighted residuals for the returned parameter values, i.e. TOTAL(DEVIATES^2). COVAR - the covariance matrix for the set of parameters returned by MPFIT. The matrix is NxN where N is the number of parameters. The square root of the diagonal elements gives the formal 1-sigma statistical errors on the parameters IF errors were treated "properly" in MYFUNC. Parameter errors are also returned in PERROR. To compute the correlation matrix, PCOR, use this example: IDL> PCOR = COV * 0 IDL> FOR i = 0, n-1 DO FOR j = 0, n-1 DO $ PCOR[i,j] = COV[i,j]/sqrt(COV[i,i]*COV[j,j]) If NOCOVAR is set or MPFIT terminated abnormally, then COVAR is set to a scalar with value !VALUES.D_NAN. DOF - number of degrees of freedom, computed as DOF = N_ELEMENTS(DEVIATES) - NFREE Note that this doesn't account for pegged parameters (see NPEGGED). ERRMSG - a string error or warning message is returned. EXTERNAL_FVEC - upon input, the function values, evaluated at START_PARAMS. This should be an M-vector, where M is the number of data points. EXTERNAL_FJAC - upon input, the Jacobian array of partial derivative values. This should be a M x N array, where M is the number of data points and N is the number of parameters. NOTE: that all FIXED or TIED parameters must *not* be included in this array. EXTERNAL_STATE - a named variable to store MPFIT-related state information between iterations (used in input and output to MPFIT). The user must not manipulate or discard this data until the final iteration is performed. FASTNORM - set this keyword to select a faster algorithm to compute sum-of-square values internally. For systems with large numbers of data points, the standard algorithm can become prohibitively slow because it cannot be vectorized well. By setting this keyword, MPFIT will run faster, but it will be more prone to floating point overflows and underflows. Thus, setting this keyword may sacrifice some stability in the fitting process. FTOL - a nonnegative input variable. Termination occurs when both the actual and predicted relative reductions in the sum of squares are at most FTOL (and STATUS is accordingly set to 1 or 3). Therefore, FTOL measures the relative error desired in the sum of squares. Default: 1D-10 FUNCTARGS - A structure which contains the parameters to be passed to the user-supplied function specified by MYFUNCT via the _EXTRA mechanism. This is the way you can pass additional data to your user-supplied function without using common blocks. Consider the following example: if FUNCTARGS = { XVAL:[1.D,2,3], YVAL:[1.D,4,9], ERRVAL:[1.D,1,1] } then the user supplied function should be declared like this: FUNCTION MYFUNCT, P, XVAL=x, YVAL=y, ERRVAL=err By default, no extra parameters are passed to the user-supplied function, but your function should accept *at least* one keyword parameter. [ This is to accomodate a limitation in IDL's _EXTRA parameter-passing mechanism. ] GTOL - a nonnegative input variable. Termination occurs when the cosine of the angle between fvec and any column of the jacobian is at most GTOL in absolute value (and STATUS is accordingly set to 4). Therefore, GTOL measures the orthogonality desired between the function vector and the columns of the jacobian. Default: 1D-10 ITERARGS - The keyword arguments to be passed to ITERPROC via the _EXTRA mechanism. This should be a structure, and is similar in operation to FUNCTARGS. Default: no arguments are passed. ITERPRINT - The name of an IDL procedure, equivalent to PRINT, that ITERPROC will use to render output. ITERPRINT should be able to accept at least four positional arguments. In addition, it should be able to accept the standard FORMAT keyword for output formatting; and the UNIT keyword, to redirect output to a logical file unit (default should be UNIT=1, standard output). These keywords are passed using the ITERARGS keyword above. The ITERPRINT procedure must accept the _EXTRA keyword. NOTE: that much formatting can be handled with the MPPRINT and MPFORMAT tags. Default: 'MPFIT_DEFPRINT' (default internal formatter) ITERPROC - The name of a procedure to be called upon each NPRINT iteration of the MPFIT routine. ITERPROC is always called in the final iteration. It should be declared in the following way: PRO ITERPROC, MYFUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $ PARINFO=parinfo, QUIET=quiet, DOF=dof, PFORMAT=pformat, $ UNIT=unit, ... ; perform custom iteration update END ITERPROC must either accept all three keyword parameters (FUNCTARGS, PARINFO and QUIET), or at least accept them via the _EXTRA keyword. MYFUNCT is the user-supplied function to be minimized, P is the current set of model parameters, ITER is the iteration number, and FUNCTARGS are the arguments to be passed to MYFUNCT. FNORM should be the chi-squared value. QUIET is set when no textual output should be printed. DOF is the number of degrees of freedom, normally the number of points less the number of free parameters. See below for documentation of PARINFO. PFORMAT is the default parameter value format. UNIT is passed on to the ITERPRINT procedure, and should indicate the file unit where log output will be sent (default: standard output). In implementation, ITERPROC can perform updates to the terminal or graphical user interface, to provide feedback while the fit proceeds. If the fit is to be stopped for any reason, then ITERPROC should set the common block variable ERROR_CODE to negative value between -15 and -1 (see MPFIT_ERROR common block below). In principle, ITERPROC should probably not modify the parameter values, because it may interfere with the algorithm's stability. In practice it is allowed. Default: an internal routine is used to print the parameter values. ITERSTOP - Set this keyword if you wish to be able to stop the fitting by hitting the predefined ITERKEYSTOP key on the keyboard. This only works if you use the default ITERPROC. ITERKEYSTOP - A keyboard key which will halt the fit (and if ITERSTOP is set and the default ITERPROC is used). ITERSTOPKEY may either be a one-character string with the desired key, or a scalar integer giving the ASCII code of the desired key. Default: 7b (control-g) NOTE: the default value of ASCI 7 (control-G) cannot be read in some windowing environments, so you must change to a printable character like 'q'. MAXITER - The maximum number of iterations to perform. If the number is exceeded, then the STATUS value is set to 5 and MPFIT returns. Default: 200 iterations NFEV - the number of MYFUNCT function evaluations performed. NFREE - the number of free parameters in the fit. This includes parameters which are not FIXED and not TIED, but it does include parameters which are pegged at LIMITS. NITER - the number of iterations completed. NOCOVAR - set this keyword to prevent the calculation of the covariance matrix before returning (see COVAR) NPEGGED - the number of free parameters which are pegged at a LIMIT. NPRINT - The frequency with which ITERPROC is called. A value of 1 indicates that ITERPROC is called with every iteration, while 2 indicates every other iteration, etc. Be aware that several Levenberg-Marquardt attempts can be made in a single iteration. Also, the ITERPROC is *always* called for the final iteration, regardless of the iteration number. Default value: 1 PARINFO - Provides a mechanism for more sophisticated constraints to be placed on parameter values. When PARINFO is not passed, then it is assumed that all parameters are free and unconstrained. Values in PARINFO are never modified during a call to MPFIT. See description above for the structure of PARINFO. Default value: all parameters are free and unconstrained. PERROR - The formal 1-sigma errors in each parameter, computed from the covariance matrix. If a parameter is held fixed, or if it touches a boundary, then the error is reported as zero. If the fit is unweighted (i.e. no errors were given, or the weights were uniformly set to unity), then PERROR will probably not represent the true parameter uncertainties. *If* you can assume that the true reduced chi-squared value is unity -- meaning that the fit is implicitly assumed to be of good quality -- then the estimated parameter uncertainties can be computed by scaling PERROR by the measured chi-squared value. DOF = N_ELEMENTS(X) - N_ELEMENTS(PARMS) ; deg of freedom PCERROR = PERROR * SQRT(BESTNORM / DOF) ; scaled uncertainties QUIET - set this keyword when no textual output should be printed by MPFIT RESDAMP - a scalar number, indicating the cut-off value of residuals where "damping" will occur. Residuals with magnitudes greater than this number will be replaced by their logarithm. This partially mitigates the so-called large residual problem inherent in least-squares solvers (as for the test problem CURVI, http://www.maxthis.com/- curviex.htm). A value of 0 indicates no damping. Default: 0 Note: RESDAMP doesn't work with AUTODERIV=0 STATUS - an integer status code is returned. All values greater than zero can represent success (however STATUS EQ 5 may indicate failure to converge). It can have one of the following values: -16 a parameter or function value has become infinite or an undefined number. This is usually a consequence of numerical overflow in the user's model function, which must be avoided. -15 to -1 these are error codes that either MYFUNCT or ITERPROC may return to terminate the fitting process (see description of MPFIT_ERROR common below). If either MYFUNCT or ITERPROC set ERROR_CODE to a negative number, then that number is returned in STATUS. Values from -15 to -1 are reserved for the user functions and will not clash with MPFIT. 0 improper input parameters. 1 both actual and predicted relative reductions in the sum of squares are at most FTOL. 2 relative error between two consecutive iterates is at most XTOL 3 conditions for STATUS = 1 and STATUS = 2 both hold. 4 the cosine of the angle between fvec and any column of the jacobian is at most GTOL in absolute value. 5 the maximum number of iterations has been reached 6 FTOL is too small. no further reduction in the sum of squares is possible. 7 XTOL is too small. no further improvement in the approximate solution x is possible. 8 GTOL is too small. fvec is orthogonal to the columns of the jacobian to machine precision. 9 A successful single iteration has been completed, and the user must supply another "EXTERNAL" evaluation of the function and its derivatives. This status indicator is neither an error nor a convergence indicator. XTOL - a nonnegative input variable. Termination occurs when the relative error between two consecutive iterates is at most XTOL (and STATUS is accordingly set to 2 or 3). Therefore, XTOL measures the relative error desired in the approximate solution. Default: 1D-10 EXAMPLE: p0 = [5.7D, 2.2, 500., 1.5, 2000.] fa = {X:x, Y:y, ERR:err} p = mpfit('MYFUNCT', p0, functargs=fa) Minimizes sum of squares of MYFUNCT. MYFUNCT is called with the X, Y, and ERR keyword parameters that are given by FUNCTARGS. The resulting parameter values are returned in p. COMMON BLOCKS: COMMON MPFIT_ERROR, ERROR_CODE User routines may stop the fitting process at any time by setting an error condition. This condition may be set in either the user's model computation routine (MYFUNCT), or in the iteration procedure (ITERPROC). To stop the fitting, the above common block must be declared, and ERROR_CODE must be set to a negative number. After the user procedure or function returns, MPFIT checks the value of this common block variable and exits immediately if the error condition has been set. This value is also returned in the STATUS keyword: values of -1 through -15 are reserved error codes for the user routines. By default the value of ERROR_CODE is zero, indicating a successful function/procedure call. COMMON MPFIT_PROFILE COMMON MPFIT_MACHAR COMMON MPFIT_CONFIG These are undocumented common blocks are used internally by MPFIT and may change in future implementations. THEORY OF OPERATION: There are many specific strategies for function minimization. One very popular technique is to use function gradient information to realize the local structure of the function. Near a local minimum the function value can be taylor expanded about x0 as follows: f(x) = f(x0) + f'(x0) . (x-x0) + (1/2) (x-x0) . f''(x0) . (x-x0) ----- --------------- ------------------------------- (1) Order 0th 1st 2nd Here f'(x) is the gradient vector of f at x, and f''(x) is the Hessian matrix of second derivatives of f at x. The vector x is the set of function parameters, not the measured data vector. One can find the minimum of f, f(xm) using Newton's method, and arrives at the following linear equation: f''(x0) . (xm-x0) = - f'(x0) (2) If an inverse can be found for f''(x0) then one can solve for (xm-x0), the step vector from the current position x0 to the new projected minimum. Here the problem has been linearized (ie, the gradient information is known to first order). f''(x0) is symmetric n x n matrix, and should be positive definite. The Levenberg - Marquardt technique is a variation on this theme. It adds an additional diagonal term to the equation which may aid the convergence properties: (f''(x0) + nu I) . (xm-x0) = -f'(x0) (2a) where I is the identity matrix. When nu is large, the overall matrix is diagonally dominant, and the iterations follow steepest descent. When nu is small, the iterations are quadratically convergent. In principle, if f''(x0) and f'(x0) are known then xm-x0 can be determined. However the Hessian matrix is often difficult or impossible to compute. The gradient f'(x0) may be easier to compute, if even by finite difference techniques. So-called quasi-Newton techniques attempt to successively estimate f''(x0) by building up gradient information as the iterations proceed. In the least squares problem there are further simplifications which assist in solving eqn (2). The function to be minimized is a sum of squares: f = Sum(hi^2) (3) where hi is the ith residual out of m residuals as described above. This can be substituted back into eqn (2) after computing the derivatives: f' = 2 Sum(hi hi') f'' = 2 Sum(hi' hj') + 2 Sum(hi hi'') (4) If one assumes that the parameters are already close enough to a minimum, then one typically finds that the second term in f'' is negligible [or, in any case, is too difficult to compute]. Thus, equation (2) can be solved, at least approximately, using only gradient information. In matrix notation, the combination of eqns (2) and (4) becomes: hT' . h' . dx = - hT' . h (5) Where h is the residual vector (length m), hT is its transpose, h' is the Jacobian matrix (dimensions n x m), and dx is (xm-x0). The user function supplies the residual vector h, and in some cases h' when it is not found by finite differences (see MPFIT_FDJAC2, which finds h and hT'). Even if dx is not the best absolute step to take, it does provide a good estimate of the best *direction*, so often a line minimization will occur along the dx vector direction. The method of solution employed by MINPACK is to form the Q . R factorization of h', where Q is an orthogonal matrix such that QT . Q = I, and R is upper right triangular. Using h' = Q . R and the ortogonality of Q, eqn (5) becomes (RT . QT) . (Q . R) . dx = - (RT . QT) . h RT . R . dx = - RT . QT . h (6) R . dx = - QT . h where the last statement follows because R is upper triangular. Here, R, QT and h are known so this is a matter of solving for dx. The routine MPFIT_QRFAC provides the QR factorization of h, with pivoting, and MPFIT_QRSOLV provides the solution for dx. REFERENCES: MINPACK-1, Jorge More', available from netlib (www.netlib.org). "Optimization Software Guide," Jorge More' and Stephen Wright, SIAM, *Frontiers in Applied Mathematics*, Number 14. More', Jorge J., "The Levenberg-Marquardt Algorithm: Implementation and Theory," in *Numerical Analysis*, ed. Watson, G. A., Lecture Notes in Mathematics 630, Springer-Verlag, 1977. MODIFICATION HISTORY: Translated from MINPACK-1 in FORTRAN, Apr-Jul 1998, CM Fixed bug in parameter limits (x vs xnew), 04 Aug 1998, CM Added PERROR keyword, 04 Aug 1998, CM Added COVAR keyword, 20 Aug 1998, CM Added NITER output keyword, 05 Oct 1998 D.L Windt, Bell Labs, windt@bell-labs.com; Made each PARINFO component optional, 05 Oct 1998 CM Analytical derivatives allowed via AUTODERIVATIVE keyword, 09 Nov 1998 Parameter values can be tied to others, 09 Nov 1998 Fixed small bugs (Wayne Landsman), 24 Nov 1998 Added better exception error reporting, 24 Nov 1998 CM Cosmetic documentation changes, 02 Jan 1999 CM Changed definition of ITERPROC to be consistent with TNMIN, 19 Jan 1999 CM Fixed bug when AUTDERIVATIVE=0. Incorrect sign, 02 Feb 1999 CM Added keyboard stop to MPFIT_DEFITER, 28 Feb 1999 CM Cosmetic documentation changes, 14 May 1999 CM IDL optimizations for speed & FASTNORM keyword, 15 May 1999 CM Tried a faster version of mpfit_enorm, 30 May 1999 CM Changed web address to cow.physics.wisc.edu, 14 Jun 1999 CM Found malformation of FDJAC in MPFIT for 1 parm, 03 Aug 1999 CM Factored out user-function call into MPFIT_CALL. It is possible, but currently disabled, to call procedures. The calling format is similar to CURVEFIT, 25 Sep 1999, CM Slightly changed mpfit_tie to be less intrusive, 25 Sep 1999, CM Fixed some bugs associated with tied parameters in mpfit_fdjac, 25 Sep 1999, CM Reordered documentation; now alphabetical, 02 Oct 1999, CM Added QUERY keyword for more robust error detection in drivers, 29 Oct 1999, CM Documented PERROR for unweighted fits, 03 Nov 1999, CM Split out MPFIT_RESETPROF to aid in profiling, 03 Nov 1999, CM Some profiling and speed optimization, 03 Nov 1999, CM Worst offenders, in order: fdjac2, qrfac, qrsolv, enorm. fdjac2 depends on user function, qrfac and enorm seem to be fully optimized. qrsolv probably could be tweaked a little, but is still <10% of total compute time. Made sure that !err was set to 0 in MPFIT_DEFITER, 10 Jan 2000, CM Fixed small inconsistency in setting of QANYLIM, 28 Jan 2000, CM Added PARINFO field RELSTEP, 28 Jan 2000, CM Converted to MPFIT_ERROR common block for indicating error conditions, 28 Jan 2000, CM Corrected scope of MPFIT_ERROR common block, CM, 07 Mar 2000 Minor speed improvement in MPFIT_ENORM, CM 26 Mar 2000 Corrected case where ITERPROC changed parameter values and parameter values were TIED, CM 26 Mar 2000 Changed MPFIT_CALL to modify NFEV automatically, and to support user procedures more, CM 26 Mar 2000 Copying permission terms have been liberalized, 26 Mar 2000, CM Catch zero value of zero a(j,lj) in MPFIT_QRFAC, 20 Jul 2000, CM (thanks to David Schlegel ) MPFIT_SETMACHAR is called only once at init; only one common block is created (MPFIT_MACHAR); it is now a structure; removed almost all CHECK_MATH calls for compatibility with IDL5 and !EXCEPT; profiling data is now in a structure too; noted some mathematical discrepancies in Linux IDL5.0, 17 Nov 2000, CM Some significant changes. New PARINFO fields: MPSIDE, MPMINSTEP, MPMAXSTEP. Improved documentation. Now PTIED constraints are maintained in the MPCONFIG common block. A new procedure to parse PARINFO fields. FDJAC2 now computes a larger variety of one-sided and two-sided finite difference derivatives. NFEV is stored in the MPCONFIG common now. 17 Dec 2000, CM Added check that PARINFO and XALL have same size, 29 Dec 2000 CM Don't call function in TERMINATE when there is an error, 05 Jan 2000 Check for float vs. double discrepancies; corrected implementation of MIN/MAXSTEP, which I still am not sure of, but now at least the correct behavior occurs *without* it, CM 08 Jan 2001 Added SCALE_FCN keyword, to allow for scaling, as for the CASH statistic; added documentation about the theory of operation, and under the QR factorization; slowly I'm beginning to understand the bowels of this algorithm, CM 10 Jan 2001 Remove MPMINSTEP field of PARINFO, for now at least, CM 11 Jan 2001 Added RESDAMP keyword, CM, 14 Jan 2001 Tried to improve the DAMP handling a little, CM, 13 Mar 2001 Corrected .PARNAME behavior in _DEFITER, CM, 19 Mar 2001 Added checks for parameter and function overflow; a new STATUS value to reflect this; STATUS values of -15 to -1 are reserved for user function errors, CM, 03 Apr 2001 DAMP keyword is now a TANH, CM, 03 Apr 2001 Added more error checking of float vs. double, CM, 07 Apr 2001 Fixed bug in handling of parameter lower limits; moved overflow checking to end of loop, CM, 20 Apr 2001 Failure using GOTO, TERMINATE more graceful if FNORM1 not defined, CM, 13 Aug 2001 Add MPPRINT tag to PARINFO, CM, 19 Nov 2001 Add DOF keyword to DEFITER procedure, and print degrees of freedom, CM, 28 Nov 2001 Add check to be sure MYFUNCT is a scalar string, CM, 14 Jan 2002 Addition of EXTERNAL_FJAC, EXTERNAL_FVEC keywords; ability to save fitter's state from one call to the next; allow '(EXTERNAL)' function name, which implies that user will supply function and Jacobian at each iteration, CM, 10 Mar 2002 Documented EXTERNAL evaluation code, CM, 10 Mar 2002 Corrected signficant bug in the way that the STEP parameter, and FIXED parameters interacted (Thanks Andrew Steffl), CM, 02 Apr 2002 Allow COVAR and PERROR keywords to be computed, even in case of '(EXTERNAL)' function, 26 May 2002 Add NFREE and NPEGGED keywords; compute NPEGGED; compute DOF using NFREE instead of n_elements(X), thanks to Kristian Kjaer, CM 11 Sep 2002 Hopefully PERROR is all positive now, CM 13 Sep 2002 Documented RELSTEP field of PARINFO (!!), CM, 25 Oct 2002 Error checking to detect missing start pars, CM 12 Apr 2003 Add DOF keyword to return degrees of freedom, CM, 30 June 2003 Always call ITERPROC in the final iteration; add ITERKEYSTOP keyword, CM, 30 June 2003 Correct bug in MPFIT_LMPAR of singularity handling, which might likely be fatal for one-parameter fits, CM, 21 Nov 2003 (with thanks to Peter Tuthill for the proper test case) Minor documentation adjustment, 03 Feb 2004, CM Correct small error in QR factorization when pivoting; document the return values of QRFAC when pivoting, 21 May 2004, CM Add MPFORMAT field to PARINFO, and correct behavior of interaction between MPPRINT and PARNAME in MPFIT_DEFITERPROC (thanks to Tim Robishaw), 23 May 2004, CM Add the ITERPRINT keyword to allow redirecting output, 26 Sep 2004, CM Correct MAXSTEP behavior in case of a negative parameter, 26 Sep 2004, CM Fix bug in the parsing of MINSTEP/MAXSTEP, 10 Apr 2005, CM Fix bug in the handling of upper/lower limits when the limit was negative (the fitting code would never "stick" to the lower limit), 29 Jun 2005, CM Small documentation update for the TIED field, 05 Sep 2005, CM Convert to IDL 5 array syntax (!), 16 Jul 2006, CM If MAXITER equals zero, then do the basic parameter checking and uncertainty analysis, but do not adjust the parameters, 15 Aug 2006, CM Added documentation, 18 Sep 2006, CM $Id: mpfit.pro,v 2.4 2009/04/15 14:43:47 bdavis Exp $
Source: src/idl_cvs/mpfit2dfun.pro
src/idl_cvs/mpfit2dfun.pro
NAME: MPFIT2DFUN AUTHOR: Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 craigm@lheamail.gsfc.nasa.gov UPDATED VERSIONs can be found on my WEB PAGE: http://cow.physics.wisc.edu/~craigm/idl/idl.html PURPOSE: Perform Levenberg-Marquardt least-squares fit to a 2-D IDL function MAJOR TOPICS: Curve and Surface Fitting CALLING SEQUENCE: parms = MPFIT2DFUN(MYFUNCT, X, Y, Z, ERR, start_parms, ...) DESCRIPTION: MPFIT2DFUN fits a user-supplied model -- in the form of an IDL function -- to a set of user-supplied data. MPFIT2DFUN calls MPFIT, the MINPACK-1 least-squares minimizer, to do the main work. MPFIT2DFUN is a specialized version for two-dimensional data. Given the data and their uncertainties, MPFIT2DFUN finds the best set of model parameters which match the data (in a least-squares sense) and returns them in an array. The user must supply the following items: - Two arrays of independent variable values ("X", "Y"). - An array of "measured" *dependent* variable values ("Z"). - An array of "measured" 1-sigma uncertainty values ("ERR"). - The name of an IDL function which computes Z given (X,Y) ("MYFUNCT"). - Starting guesses for all of the parameters ("START_PARAMS"). There are very few restrictions placed on X, Y, Z, or MYFUNCT. Simply put, MYFUNCT must map the (X,Y) values into Z values given the model parameters. The (X,Y) values are usually the independent X and Y coordinate positions in the two dimensional plane, but need not be. MPFIT2DFUN carefully avoids passing large arrays where possible to improve performance. See below for an example of usage. USER FUNCTION The user must define a function which returns the model value. For applications which use finite-difference derivatives -- the default -- the user function should be declared in the following way: FUNCTION MYFUNCT, X, Y, P ; The independent variables are X and Y ; Parameter values are passed in "P" ZMOD = ... computed model values at (X,Y) ... return, ZMOD END The returned array YMOD must have the same dimensions and type as the "measured" Z values. User functions may also indicate a fatal error condition using the ERROR_CODE common block variable, as described below under the MPFIT_ERROR common block definition. See the discussion under "ANALYTIC DERIVATIVES" and AUTODERIVATIVE in MPFIT.PRO if you wish to compute the derivatives for yourself. AUTODERIVATIVE is accepted and passed directly to MPFIT. The user function must accept one additional parameter, DP, which contains the derivative of the user function with respect to each parameter at each data point, as described in MPFIT.PRO. CREATING APPROPRIATELY DIMENSIONED INDEPENDENT VARIABLES The user must supply appropriate independent variables to MPFIT2DFUN. For image fitting applications, this variable should be two-dimensional *arrays* describing the X and Y positions of every *pixel*. [ Thus any two dimensional sampling is permitted, including irregular sampling. ] If the sampling is regular, then the x coordinates are the same for each row, and the y coordinates are the same for each column. Call the x-row and y-column coordinates XR and YC respectively. You can then compute X and Y as follows: X = XR # (YC*0 + 1) eqn. 1 Y = (XR*0 + 1) # YC eqn. 2 For example, if XR and YC have the following values: XR = [ 1, 2, 3, 4, 5,] ;; X positions of one row of pixels YC = [ 15,16,17 ] ;; Y positions of one column of pixels Then using equations 1 and 2 above will give these values to X and Y: X : 1 2 3 4 5 ;; X positions of all pixels 1 2 3 4 5 1 2 3 4 5 Y : 15 15 15 15 15 ;; Y positions of all pixels 16 16 16 16 16 17 17 17 17 17 Using the above technique is suggested, but *not* required. You can do anything you wish with the X and Y values. This technique only makes it easier to compute your model function values. CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD The behavior of MPFIT can be modified with respect to each parameter to be fitted. A parameter value can be fixed; simple boundary constraints can be imposed; limitations on the parameter changes can be imposed; properties of the automatic derivative can be modified; and parameters can be tied to one another. These properties are governed by the PARINFO structure, which is passed as a keyword parameter to MPFIT. PARINFO should be an array of structures, one for each parameter. Each parameter is associated with one element of the array, in numerical order. The structure can have the following entries (none are required): .VALUE - the starting parameter value (but see the START_PARAMS parameter for more information). .FIXED - a boolean value, whether the parameter is to be held fixed or not. Fixed parameters are not varied by MPFIT, but are passed on to MYFUNCT for evaluation. .LIMITED - a two-element boolean array. If the first/second element is set, then the parameter is bounded on the lower/upper side. A parameter can be bounded on both sides. Both LIMITED and LIMITS must be given together. .LIMITS - a two-element float or double array. Gives the parameter limits on the lower and upper sides, respectively. Zero, one or two of these values can be set, depending on the values of LIMITED. Both LIMITED and LIMITS must be given together. .PARNAME - a string, giving the name of the parameter. The fitting code of MPFIT does not use this tag in any way. However, the default ITERPROC will print the parameter name if available. .STEP - the step size to be used in calculating the numerical derivatives. If set to zero, then the step size is computed automatically. Ignored when AUTODERIVATIVE=0. This value is superceded by the RELSTEP value. .RELSTEP - the *relative* step size to be used in calculating the numerical derivatives. This number is the fractional size of the step, compared to the parameter value. This value supercedes the STEP setting. If the parameter is zero, then a default step size is chosen. .MPSIDE - the sidedness of the finite difference when computing numerical derivatives. This field can take four values: 0 - one-sided derivative computed automatically 1 - one-sided derivative (f(x+h) - f(x) )/h -1 - one-sided derivative (f(x) - f(x-h))/h 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h) Where H is the STEP parameter described above. The "automatic" one-sided derivative method will chose a direction for the finite difference which does not violate any constraints. The other methods do not perform this check. The two-sided method is in principle more precise, but requires twice as many function evaluations. Default: 0. .MPMINSTEP - the minimum change to be made in the parameter value. During the fitting process, the parameter will be changed by multiples of this value. The actual step is computed as: DELTA1 = MPMINSTEP*ROUND(DELTA0/MPMINSTEP) where DELTA0 and DELTA1 are the estimated parameter changes before and after this constraint is applied. Note that this constraint should be used with care since it may cause non-converging, oscillating solutions. A value of 0 indicates no minimum. Default: 0. .MPMAXSTEP - the maximum change to be made in the parameter value. During the fitting process, the parameter will never be changed by more than this value. A value of 0 indicates no maximum. Default: 0. .TIED - a string expression which "ties" the parameter to other free or fixed parameters. Any expression involving constants and the parameter array P are permitted. Example: if parameter 2 is always to be twice parameter 1 then use the following: parinfo[2].tied = '2 * P[1]'. Since they are totally constrained, tied parameters are considered to be fixed; no errors are computed for them. [ NOTE: the PARNAME can't be used in expressions. ] Future modifications to the PARINFO structure, if any, will involve adding structure tags beginning with the two letters "MP". Therefore programmers are urged to avoid using tags starting with the same letters; otherwise they are free to include their own fields within the PARINFO structure, and they will be ignored. PARINFO Example: parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $ limits:[0.D,0]}, 5) parinfo[0].fixed = 1 parinfo[4].limited(0) = 1 parinfo[4].limits(0) = 50.D parinfo[*].value = [5.7D, 2.2, 500., 1.5, 2000.] A total of 5 parameters, with starting values of 5.7, 2.2, 500, 1.5, and 2000 are given. The first parameter is fixed at a value of 5.7, and the last parameter is constrained to be above 50. INPUTS: MYFUNCT - a string variable containing the name of an IDL function. This function computes the "model" Z values given the X,Y values and model parameters, as described above. X - Array of "X" independent variable values, as described above. These values are passed directly to the fitting function unmodified. Y - Array of "Y" independent variable values, as described above. X and Y should have the same data type. Z - Array of "measured" dependent variable values. Z should have the same data type as X and Y. The function MYFUNCT should map (X,Y)->Z. ERR - Array of "measured" 1-sigma uncertainties. ERR should have the same data type as Z. ERR is ignored if the WEIGHTS keyword is specified. START_PARAMS - An array of starting values for each of the parameters of the model. The number of parameters should be fewer than the number of measurements. Also, the parameters should have the same data type as the measurements (double is preferred). This parameter is optional if the PARINFO keyword is used (see MPFIT). The PARINFO keyword provides a mechanism to fix or constrain individual parameters. If both START_PARAMS and PARINFO are passed, then the starting *value* is taken from START_PARAMS, but the *constraints* are taken from PARINFO. RETURNS: Returns the array of best-fit parameters. KEYWORD PARAMETERS: BESTNORM - the value of the summed, squared, weighted residuals for the returned parameter values, i.e. the chi-square value. COVAR - the covariance matrix for the set of parameters returned by MPFIT. The matrix is NxN where N is the number of parameters. The square root of the diagonal elements gives the formal 1-sigma statistical errors on the parameters IF errors were treated "properly" in MYFUNC. Parameter errors are also returned in PERROR. To compute the correlation matrix, PCOR, use this: IDL> PCOR = COV * 0 IDL> FOR i = 0, n-1 DO FOR j = 0, n-1 DO $ PCOR[i,j] = COV[i,j]/sqrt(COV[i,i]*COV[j,j]) If NOCOVAR is set or MPFIT terminated abnormally, then COVAR is set to a scalar with value !VALUES.D_NAN. DOF - number of degrees of freedom, computed as DOF = N_ELEMENTS(DEVIATES) - NFREE Note that this doesn't account for pegged parameters (see NPEGGED). ERRMSG - a string error or warning message is returned. FTOL - a nonnegative input variable. Termination occurs when both the actual and predicted relative reductions in the sum of squares are at most FTOL (and STATUS is accordingly set to 1 or 3). Therefore, FTOL measures the relative error desired in the sum of squares. Default: 1D-10 FUNCTARGS - A structure which contains the parameters to be passed to the user-supplied function specified by MYFUNCT via the _EXTRA mechanism. This is the way you can pass additional data to your user-supplied function without using common blocks. By default, no extra parameters are passed to the user-supplied function. GTOL - a nonnegative input variable. Termination occurs when the cosine of the angle between fvec and any column of the jacobian is at most GTOL in absolute value (and STATUS is accordingly set to 4). Therefore, GTOL measures the orthogonality desired between the function vector and the columns of the jacobian. Default: 1D-10 ITERARGS - The keyword arguments to be passed to ITERPROC via the _EXTRA mechanism. This should be a structure, and is similar in operation to FUNCTARGS. Default: no arguments are passed. ITERPROC - The name of a procedure to be called upon each NPRINT iteration of the MPFIT routine. It should be declared in the following way: PRO ITERPROC, MYFUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $ PARINFO=parinfo, QUIET=quiet, ... ; perform custom iteration update END ITERPROC must either accept all three keyword parameters (FUNCTARGS, PARINFO and QUIET), or at least accept them via the _EXTRA keyword. MYFUNCT is the user-supplied function to be minimized, P is the current set of model parameters, ITER is the iteration number, and FUNCTARGS are the arguments to be passed to MYFUNCT. FNORM should be the chi-squared value. QUIET is set when no textual output should be printed. See below for documentation of PARINFO. In implementation, ITERPROC can perform updates to the terminal or graphical user interface, to provide feedback while the fit proceeds. If the fit is to be stopped for any reason, then ITERPROC should set the common block variable ERROR_CODE to negative value (see MPFIT_ERROR common block below). In principle, ITERPROC should probably not modify the parameter values, because it may interfere with the algorithm's stability. In practice it is allowed. Default: an internal routine is used to print the parameter values. MAXITER - The maximum number of iterations to perform. If the number is exceeded, then the STATUS value is set to 5 and MPFIT returns. Default: 200 iterations NFEV - the number of MYFUNCT function evaluations performed. NITER - the number of iterations completed. NOCOVAR - set this keyword to prevent the calculation of the covariance matrix before returning (see COVAR) NPRINT - The frequency with which ITERPROC is called. A value of 1 indicates that ITERPROC is called with every iteration, while 2 indicates every other iteration, etc. Note that several Levenberg-Marquardt attempts can be made in a single iteration. Default value: 1 PARINFO - Provides a mechanism for more sophisticated constraints to be placed on parameter values. When PARINFO is not passed, then it is assumed that all parameters are free and unconstrained. Values in PARINFO are never modified during a call to MPFIT. See description above for the structure of PARINFO. Default value: all parameters are free and unconstrained. PERROR - The formal 1-sigma errors in each parameter, computed from the covariance matrix. If a parameter is held fixed, or if it touches a boundary, then the error is reported as zero. If the fit is unweighted (i.e. no errors were given, or the weights were uniformly set to unity), then PERROR will probably not represent the true parameter uncertainties. *If* you can assume that the true reduced chi-squared value is unity -- meaning that the fit is implicitly assumed to be of good quality -- then the estimated parameter uncertainties can be computed by scaling PERROR by the measured chi-squared value. DOF = N_ELEMENTS(Z) - N_ELEMENTS(PARMS) ; deg of freedom PCERROR = PERROR * SQRT(BESTNORM / DOF) ; scaled uncertainties QUIET - set this keyword when no textual output should be printed by MPFIT STATUS - an integer status code is returned. All values other than zero can represent success. It can have one of the following values: 0 improper input parameters. 1 both actual and predicted relative reductions in the sum of squares are at most FTOL. 2 relative error between two consecutive iterates is at most XTOL 3 conditions for STATUS = 1 and STATUS = 2 both hold. 4 the cosine of the angle between fvec and any column of the jacobian is at most GTOL in absolute value. 5 the maximum number of iterations has been reached 6 FTOL is too small. no further reduction in the sum of squares is possible. 7 XTOL is too small. no further improvement in the approximate solution x is possible. 8 GTOL is too small. fvec is orthogonal to the columns of the jacobian to machine precision. WEIGHTS - Array of weights to be used in calculating the chi-squared value. If WEIGHTS is specified then the ERR parameter is ignored. The chi-squared value is computed as follows: CHISQ = TOTAL( (Z-MYFUNCT(X,Y,P))^2 * ABS(WEIGHTS) ) Here are common values of WEIGHTS: 1D/ERR^2 - Normal weighting (ERR is the measurement error) 1D/Z - Poisson weighting (counting statistics) 1D - Unweighted XTOL - a nonnegative input variable. Termination occurs when the relative error between two consecutive iterates is at most XTOL (and STATUS is accordingly set to 2 or 3). Therefore, XTOL measures the relative error desired in the approximate solution. Default: 1D-10 YFIT - the best-fit model function, as returned by MYFUNCT. EXAMPLE: p = [2.2D, -0.7D, 1.4D, 3000.D] x = (dindgen(200)*0.1 - 10.) # (dblarr(200) + 1) y = (dblarr(200) + 1) # (dindgen(200)*0.1 - 10.) zi = gauss2(x, y, p) sz = sqrt(zi>1) z = zi + randomn(seed, 200, 200) * sz p0 = [0D, 0D, 1D, 10D] p = mpfit2dfun('GAUSS2', x, y, z, sz, p0) Generates a synthetic data set with a Gaussian peak, and Poisson statistical uncertainty. Then the same function (but different starting parameters) is fitted to the data to see how close we can get. It is especially worthy to notice that the X and Y values are created as full images, so that a coordinate is attached to each pixel independently. This is the format that GAUSS2 accepts, and the easiest for you to use in your own functions. COMMON BLOCKS: COMMON MPFIT_ERROR, ERROR_CODE User routines may stop the fitting process at any time by setting an error condition. This condition may be set in either the user's model computation routine (MYFUNCT), or in the iteration procedure (ITERPROC). To stop the fitting, the above common block must be declared, and ERROR_CODE must be set to a negative number. After the user procedure or function returns, MPFIT checks the value of this common block variable and exits immediately if the error condition has been set. By default the value of ERROR_CODE is zero, indicating a successful function/procedure call. REFERENCES: MINPACK-1, Jorge More', available from netlib (www.netlib.org). "Optimization Software Guide," Jorge More' and Stephen Wright, SIAM, *Frontiers in Applied Mathematics*, Number 14. MODIFICATION HISTORY: Written, transformed from MPFITFUN, 26 Sep 1999, CM Alphabetized documented keywords, 02 Oct 1999, CM Added example, 02 Oct 1999, CM Tried to clarify definitions of X and Y, 29 Oct 1999, CM Added QUERY keyword and query checking of MPFIT, 29 Oct 1999, CM Check to be sure that X, Y and Z are present, 02 Nov 1999, CM Documented PERROR for unweighted fits, 03 Nov 1999, CM Changed to ERROR_CODE for error condition, 28 Jan 2000, CM Copying permission terms have been liberalized, 26 Mar 2000, CM Propagated improvements from MPFIT, 17 Dec 2000, CM Documented RELSTEP field of PARINFO (!!), CM, 25 Oct 2002 Add DOF keyword to return degrees of freedom, CM, 23 June 2003 Minor documentation adjustment, 03 Feb 2004, CM Fix the example to prevent zero errorbars, 28 Mar 2005, CM Defend against users supplying strangely dimensioned X and Y, 29 Jun 2005, CM Convert to IDL 5 array syntax (!), 16 Jul 2006, CM $Id: mpfit2dfun.pro,v 2.2 2009/04/15 14:43:47 bdavis Exp $
Source: src/idl_cvs/mpfit2dpeak.pro
src/idl_cvs/mpfit2dpeak.pro
NAME: MPFIT2DPEAK AUTHOR: Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 craigm@lheamail.gsfc.nasa.gov UPDATED VERSIONs can be found on my WEB PAGE: http://cow.physics.wisc.edu/~craigm/idl/idl.html PURPOSE: Fit a gaussian, lorentzian or Moffat model to data MAJOR TOPICS: Curve and Surface Fitting CALLING SEQUENCE: yfit = MPFIT2DPEAK(Z, A [, X, Y, /TILT ...] ) DESCRIPTION: MPFIT2DPEAK fits a gaussian, lorentzian or Moffat model using the non-linear least squares fitter MPFIT. MPFIT2DPEAK is meant to be a drop-in replacement for IDL's GAUSS2DFIT function (and requires MPFIT and MPFIT2DFUN). The choice of the fitting function is determined by the keywords GAUSSIAN, LORENTZIAN and MOFFAT. By default the gaussian model function is used. [ The Moffat function is a modified Lorentzian with variable power law index. ] The two-dimensional peak has independent semimajor and semiminor axes, with an optional rotation term activated by setting the TILT keyword. The baseline is assumed to be a constant. GAUSSIAN A[0] + A[1]*exp(-0.5*u) LORENTZIAN A[0] + A[1]/(u + 1) MOFFAT A[0] + A[1]/(u + 1)^A[7] u = ( (x-A[4])/A[2] )^2 + ( (y-A[5])/A[3] )^2 where x and y are cartesian coordinates in rotated coordinate system if TILT keyword is set. The returned parameter array elements have the following meanings: A[0] Constant baseline level A[1] Peak value A[2] Peak half-width (x) -- gaussian sigma or half-width at half-max A[3] Peak half-width (y) -- gaussian sigma or half-width at half-max A[4] Peak centroid (x) A[5] Peak centroid (y) A[6] Rotation angle (radians) if TILT keyword set A[7] Moffat power law index if MOFFAT keyword set By default the initial starting values for the parameters A are estimated from the data. However, explicit starting values can be supplied using the ESTIMATES keyword. Also, error or weighting values can optionally be provided; otherwise the fit is unweighted. RESTRICTIONS: If no starting parameter ESTIMATES are provided, then MPFIT2DPEAK attempts to estimate them from the data. This is not a perfect science; however, the author believes that the technique implemented here is more robust than the one used in IDL's GAUSS2DFIT. The author has tested cases of strong peaks, noisy peaks and broad peaks, all with success. INPUTS: Z - Two dimensional array of "measured" dependent variable values. Z should be of the same type and dimension as (X # Y). X - Optional vector of x positions for a single row of Z. X[i] should provide the x position of Z[i,*] Default: X values are integer increments from 0 to NX-1 Y - Optional vector of y positions for a single column of Z. Y[j] should provide the y position of Z[*,j] Default: Y values are integer increments from 0 to NY-1 OUTPUTS: A - Upon return, an array of best fit parameter values. See the table above for the meanings of each parameter element. RETURNS: Returns the best fitting model function as a 2D array. KEYWORDS: ** NOTE ** Additional keywords such as PARINFO, BESTNORM, and STATUS are accepted by MPFIT2DPEAK but not documented here. Please see the documentation for MPFIT for the description of these advanced options. CHISQ - the value of the summed squared residuals for the returned parameter values. CIRCULAR - if set, then the peak profile is assumed to be azimuthally symmetric. When set, the parameters A[2) and A[3) will be identical and the TILT keyword will have no effect. DOF - number of degrees of freedom, computed as DOF = N_ELEMENTS(DEVIATES) - NFREE Note that this doesn't account for pegged parameters (see NPEGGED). ERROR - upon input, the measured 1-sigma uncertainties in the "Z" values. If no ERROR or WEIGHTS are given, then the fit is unweighted. ESTIMATES - Array of starting values for each parameter of the model. Default: parameter values are estimated from data. GAUSSIAN - if set, fit a gaussian model function. The Default. LORENTZIAN - if set, fit a lorentzian model function. MOFFAT - if set, fit a Moffat model function. MEASURE_ERRORS - synonym for ERRORS, for consistency with built-in IDL fitting routines. NEGATIVE - if set, and ESTIMATES is not provided, then MPFIT2DPEAK will assume that a negative peak is present -- a valley. Specifying this keyword is not normally required, since MPFIT2DPEAK can determine this automatically. NFREE - the number of free parameters in the fit. This includes parameters which are not FIXED and not TIED, but it does include parameters which are pegged at LIMITS. PERROR - upon return, the 1-sigma uncertainties of the parameter values A. These values are only meaningful if the ERRORS or WEIGHTS keywords are specified properly. If the fit is unweighted (i.e. no errors were given, or the weights were uniformly set to unity), then PERROR will probably not represent the true parameter uncertainties. *If* you can assume that the true reduced chi-squared value is unity -- meaning that the fit is implicitly assumed to be of good quality -- then the estimated parameter uncertainties can be computed by scaling PERROR by the measured chi-squared value. DOF = N_ELEMENTS(Z) - N_ELEMENTS(A) ; deg of freedom PCERROR = PERROR * SQRT(BESTNORM / DOF) ; scaled uncertainties QUIET - if set then diagnostic fitting messages are suppressed. Default: QUIET=1 (i.e., no diagnostics) SIGMA - synonym for PERROR (1-sigma parameter uncertainties), for compatibility with GAUSSFIT. Do not confuse this with the Gaussian "sigma" width parameter. TILT - if set, then the major and minor axes of the peak profile are allowed to rotate with respect to the image axes. Parameter A[6] will be set to the clockwise rotation angle of the A[2] axis in radians. WEIGHTS - Array of weights to be used in calculating the chi-squared value. If WEIGHTS is specified then the ERR parameter is ignored. The chi-squared value is computed as follows: CHISQ = TOTAL( (Z-MYFUNCT(X,Y,P))^2 * ABS(WEIGHTS) ) Here are common values of WEIGHTS: 1D/ERR^2 - Normal weighting (ERR is the measurement error) 1D/Y - Poisson weighting (counting statistics) 1D - Unweighted The ERROR keyword takes precedence over any WEIGHTS keyword values. If no ERROR or WEIGHTS are given, then the fit is unweighted. EXAMPLE: ; Construct a sample gaussian surface in range [-5,5] centered at [2,-3] x = findgen(100)*0.1 - 5. & y = x xx = x # (y*0 + 1) yy = (x*0 + 1) # y rr = sqrt((xx-2.)^2 + (yy+3.)^2) ; Gaussian surface with sigma=0.5, peak value of 3, noise with sigma=0.2 z = 3.*exp(-(rr/0.5)^2) + randomn(seed,100,100)*.2 ; Fit gaussian parameters A zfit = mpfit2dpeak(z, a, x, y) REFERENCES: MINPACK-1, Jorge More', available from netlib (www.netlib.org). "Optimization Software Guide," Jorge More' and Stephen Wright, SIAM, *Frontiers in Applied Mathematics*, Number 14. MODIFICATION HISTORY: New algorithm for estimating starting values, CM, 31 Oct 1999 Documented, 02 Nov 1999 Small documentation fixes, 02 Nov 1999 Documented PERROR for unweighted fits, 03 Nov 1999, CM Copying permission terms have been liberalized, 26 Mar 2000, CM Small cosmetic changes, 21 Sep 2000, CM Corrected bug introduced by cosmetic changes, 11 Oct 2000, CM :-) Added POSITIVE keyword, 17 Nov 2000, CM Removed TILT in common, in favor of FUNCTARGS approach, 23 Nov 2000, CM Added SYMMETRIC keyword, documentation for TILT, and an example, 24 Nov 2000, CM Changed SYMMETRIC to CIRCULAR, 17 Dec 2000, CM Really change SYMMETRIC to CIRCULAR!, 13 Sep 2002, CM Add error messages for SYMMETRIC and CIRCLE, 08 Nov 2002, CM Make more consistent with comparable IDL routines, 30 Jun 2003, CM Defend against users supplying strangely dimensioned X and Y, 29 Jun 2005, CM Convert to IDL 5 array syntax (!), 16 Jul 2006, CM $Id: mpfit2dpeak.pro,v 2.2 2009/04/15 14:43:47 bdavis Exp $
Source: src/idl_cvs/mpfitellipse.pro
src/idl_cvs/mpfitellipse.pro
NAME: MPFITELLIPSE AUTHOR: Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 craigm@lheamail.gsfc.nasa.gov UPDATED VERSIONs can be found on my WEB PAGE: http://cow.physics.wisc.edu/~craigm/idl/idl.html PURPOSE: Approximate fit to points forming an ellipse MAJOR TOPICS: Curve and Surface Fitting CALLING SEQUENCE: parms = MPFITELLIPSE(X, Y, start_parms, [/TILT, WEIGHTS=wts, ...]) DESCRIPTION: MPFITELLIPSE fits a closed elliptical or circular curve to a two dimensional set of data points. The user specifies the X and Y positions of the points, and an optional set of weights. The ellipse may also be tilted at an arbitrary angle. IMPORTANT NOTE: this fitting program performs simple ellipse fitting. It will not work well for ellipse data with high eccentricity. More robust answers can usually be obtained with "orthogonal distance regression." (See FORTRAN package ODRPACK on netlib.org for more information). The best fitting ellipse parameters are returned from by MPFITELLIPSE as a vector, whose values are: P[0] Ellipse semi axis 1 P[1] Ellipse semi axis 2 ( = P[0] if CIRCLE keyword set) P[2] Ellipse center - x value P[3] Ellipse center - y value P[4] Ellipse rotation angle (radians) if TILT keyword set If the TILT keyword is set, then the P[0] is meant to be the semi-major axis, and P[1] is the semi-minor axis, and P[4] represents the tilt of the semi-major axis with respect to the X axis. If the TILT keyword is not set, the P[0] and P[1] represent the ellipse semi-axes in the X and Y directions, respectively. The returned semi-axis lengths should always be positive. The user may specify an initial set of trial parameters, but by default MPFITELLIPSE will estimate the parameters automatically. Users should be aware that in the presence of large amounts of noise, namely when the measurement error becomes significant compared to the ellipse axis length, then the estimated parameters become unreliable. Generally speaking the computed axes will overestimate the true axes. For example when (SIGMA_R/R) becomes 0.5, the radius of the ellipse is overestimated by about 40%. This unreliability is also pronounced if the ellipse has high eccentricity, as noted above. Users can weight their data as they see appropriate. However, the following prescription for the weighting may serve as a good starting point, and appeared to produce results comparable to the typical chi-squared value. WEIGHTS = 0.75/(SIGMA_X^2 + SIGMA_Y^2) where SIGMA_X and SIGMA_Y are the measurement error vectors in the X and Y directions respectively. However, this has not been robustly tested, and it should be pointed out that this weighting may only be appropriate for a set of points whose measurement errors are comparable. If a more robust estimation of the parameter values is needed, the so-called orthogonal distance regression package should be used (ODRPACK, available in FORTRAN at www.netlib.org). INPUTS: X - measured X positions of the points in the ellipse. Y - measured Y positions of the points in the ellipse. START_PARAMS - an array of starting values for the ellipse parameters, as described above. This parameter is optional; if not specified by the user, then the ellipse parameters are estimated automatically from the properties of the data. RETURNS: Returns the best fitting model ellipse parameters. Returned values are undefined if STATUS indicates an error condition. KEYWORDS: ** NOTE ** Additional keywords such as PARINFO, BESTNORM, and STATUS are accepted by MPFITELLIPSE but not documented here. Please see the documentation for MPFIT for the description of these advanced options. CIRCULAR - if set, then the curve is assumed to be a circle instead of ellipse. When set, the parameters P[0] and P[1] will be identical and the TILT keyword will have no effect. PERROR - upon return, the 1-sigma uncertainties of the returned ellipse parameter values. These values are only meaningful if the WEIGHTS keyword is specified properly. If the fit is unweighted (i.e. no errors were given, or the weights were uniformly set to unity), then PERROR will probably not represent the true parameter uncertainties. If STATUS indicates an error condition, then PERROR is undefined. QUIET - if set then diagnostic fitting messages are suppressed. Default: QUIET=1 (i.e., no diagnostics] STATUS - an integer status code is returned. All values greater than zero can represent success (however STATUS EQ 5 may indicate failure to converge). Please see MPFIT for the definitions of status codes. TILT - if set, then the major and minor axes of the ellipse are allowed to rotate with respect to the data axes. Parameter P[4] will be set to the clockwise rotation angle of the P[0] axis in radians, as measured from the +X axis. P[4] should be in the range 0 to !dpi. WEIGHTS - Array of weights to be used in calculating the chi-squared value. The chi-squared value is computed as follows: CHISQ = TOTAL( (Z-MYFUNCT(X,Y,P))^2 * ABS(WEIGHTS)^2 ) Users may wish to follow the guidelines for WEIGHTS described above. EXAMPLE: ; Construct a set of points on an ellipse, with some noise npts = 50. ph0 = 2*!pi*randomu(seed,npts) x = 50. + 32.*cos(ph0) + 4.0*randomn(seed, npts) y = -75. + 65.*sin(ph0) + 0.1*randomn(seed, npts) ; Compute weights function weights = 1./npts weights = sqrt(1./npts) weights = 0.75/(4.0^2 + 0.1^2) ; Fit ellipse and plot result p = mpfitellipse(x, y) phi = dindgen(101)*2D*!dpi/100 plot, x, y, psym=1 oplot, p[2]+p[0]*cos(phi), p[3]+p[1]*sin(phi), color='ff'xl ; Fit ellipse and plot result - WITH TILT p = mpfitellipse(x, y, /tilt, BESTNORM=chiSq, WEIGHTS=weights) phi = dindgen(101)*2D*!dpi/100 ; New parameter P[4] gives tilt of ellipse w.r.t. coordinate axes ; We must rotate a standard ellipse to this new orientation xm = p[2] + p[0]*cos(phi)*cos(p[4]) + p[1]*sin(phi)*sin(p[4]) ym = p[3] - p[0]*cos(phi)*sin(p[4]) + p[1]*sin(phi)*cos(p[4]) plot, x, y, psym=1 oplot, xm, ym, color='ff'xl REFERENCES: MINPACK-1, Jorge More', available from netlib (www.netlib.org). "Optimization Software Guide," Jorge More' and Stephen Wright, SIAM, *Frontiers in Applied Mathematics*, Number 14. MODIFICATION HISTORY: Ported from MPFIT2DPEAK, 17 Dec 2000, CM More documentation, 11 Jan 2001, CM Example corrected, 18 Nov 2001, CM Change CIRCLE keyword to the correct CIRCULAR keyword, 13 Sep 2002, CM Add error messages for SYMMETRIC and CIRCLE, 08 Nov 2002, CM Found small error in computation of _EVAL (when CIRCULAR) was set; sanity check when CIRCULAR is set, 21 Jan 2003, CM Convert to IDL 5 array syntax (!), 16 Jul 2006, CM Move STRICTARR compile option inside each function/procedure, 9 Oct 2006 Add disclaimer about the suitability of this program for fitting ellipses, 17 Sep 2007, CM Clarify documentation of TILT angle; make sure output contains semi-major axis first, followed by semi-minor; make sure that semi-axes are always positive (and can handle negative inputs) 17 Sep 2007, CM Output tilt angle is now in range 0 to !DPI, 20 Sep 2007, CM Some documentation clarifications, including to remove reference to the "ERR" keyword, which does not exist, 17 Jan 2008, CM Swapping of P[0] and P[1] only occurs if /TILT is set, 06 Nov 2009, CM Document an example of how to plot a tilted ellipse, 09 Nov 2009, CM Check for MPFIT error conditions and return immediately, 23 Jan 2010, CM $Id: mpfitellipse.pro,v 2.4 2011/01/18 21:15:25 bdavis Exp $
Source: src/idl_cvs/mpfitexpr.pro
src/idl_cvs/mpfitexpr.pro
NAME: MPFITEXPR AUTHOR: Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 craigm@lheamail.gsfc.nasa.gov UPDATED VERSIONs can be found on my WEB PAGE: http://cow.physics.wisc.edu/~craigm/idl/idl.html PURPOSE: Perform Levenberg-Marquardt least-squares fit to arbitrary expression MAJOR TOPICS: Curve and Surface Fitting CALLING SEQUENCE: MYFUNCT = 'X*(1-X)+3' parms = MPFITEXPR(MYFUNCT, XVAL, YVAL, ERR, start_parms, ...) DESCRIPTION: MPFITEXPR fits a user-supplied model -- in the form of an arbitrary IDL expression -- to a set of user-supplied data. MPFITEXPR calls MPFIT, the MINPACK-1 least-squares minimizer, to do the main work. Given the data and their uncertainties, MPFITEXPR finds the best set of model parameters which match the data (in a least-squares sense) and returns them in an array. The user must supply the following items: - An array of independent variable values ("X"). - An array of "measured" *dependent* variable values ("Y"). - An array of "measured" 1-sigma uncertainty values ("ERR"). - A text IDL expression which computes Y given X. - Starting guesses for all of the parameters ("START_PARAMS"). There are very few restrictions placed on X, Y or the expression of the model. Simply put, the expression must map the "X" values into "Y" values given the model parameters. The "X" values may represent any independent variable (not just Cartesian X), and indeed may be multidimensional themselves. For example, in the application of image fitting, X may be a 2xN array of image positions. Some rules must be obeyed in constructing the expression. First, the independent variable name *MUST* be "X" in the expression, regardless of the name of the variable being passed to MPFITEXPR. This is demonstrated in the above calling sequence, where the X variable passed in is called "XVAL" but the expression still refers to "X". Second, parameter values must be referred to as an array named "P". If you do not pass in starting values for the model parameters, MPFITEXPR will attempt to determine the number of parameters you intend to have (it does this by looking for references to the array variable named "P"). When no starting values are passed in, the values are assumed to start at zero. MPFITEXPR carefully avoids passing large arrays where possible to improve performance. See below for an example of usage. EVALUATING EXPRESSIONS This source module also provides a function called MPEVALEXPR. You can use this function to evaluate your expression, given a list of parameters. This is one of the easier ways to compute the model once the best-fit parameters have been found. Here is an example: YMOD = MPEVALEXPR(MYFUNCT, XVAL, PARMS) where MYFUNCT is the expression (see MYFUNCT below), XVAL is the list of "X" values, and PARMS is an array of parameters. The returned array YMOD contains the expression MYFUNCT evaluated at each point in XVAL. PASSING PRIVATE DATA TO AN EXPRESSION The most complicated optimization problems typically involve other external parameters, in addition to the fitted parameters. While it is extremely easy to rewrite an expression dynamically, in case one of the external parameters changes, the other possibility is to use the PRIVATE data structure. The user merely passes a structure to the FUNCTARGS keyword. The user expression receives this value as the variable PRIVATE. MPFITEXPR nevers accesses this variable so it can contain any desired values. Usually it would be an IDL structure so that any named external parameters can be passed to the expression. CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD The behavior of MPFIT can be modified with respect to each parameter to be fitted. A parameter value can be fixed; simple boundary constraints can be imposed; limitations on the parameter changes can be imposed; properties of the automatic derivative can be modified; and parameters can be tied to one another. These properties are governed by the PARINFO structure, which is passed as a keyword parameter to MPFIT. PARINFO should be an array of structures, one for each parameter. Each parameter is associated with one element of the array, in numerical order. The structure can have the following entries (none are required): .VALUE - the starting parameter value (but see the START_PARAMS parameter for more information). .FIXED - a boolean value, whether the parameter is to be held fixed or not. Fixed parameters are not varied by MPFIT, but are passed on to MYFUNCT for evaluation. .LIMITED - a two-element boolean array. If the first/second element is set, then the parameter is bounded on the lower/upper side. A parameter can be bounded on both sides. Both LIMITED and LIMITS must be given together. .LIMITS - a two-element float or double array. Gives the parameter limits on the lower and upper sides, respectively. Zero, one or two of these values can be set, depending on the values of LIMITED. Both LIMITED and LIMITS must be given together. .PARNAME - a string, giving the name of the parameter. The fitting code of MPFIT does not use this tag in any way. However, the default ITERPROC will print the parameter name if available. .STEP - the step size to be used in calculating the numerical derivatives. If set to zero, then the step size is computed automatically. Ignored when AUTODERIVATIVE=0. This value is superceded by the RELSTEP value. .RELSTEP - the *relative* step size to be used in calculating the numerical derivatives. This number is the fractional size of the step, compared to the parameter value. This value supercedes the STEP setting. If the parameter is zero, then a default step size is chosen. .MPSIDE - the sidedness of the finite difference when computing numerical derivatives. This field can take four values: 0 - one-sided derivative computed automatically 1 - one-sided derivative (f(x+h) - f(x) )/h -1 - one-sided derivative (f(x) - f(x-h))/h 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h) Where H is the STEP parameter described above. The "automatic" one-sided derivative method will chose a direction for the finite difference which does not violate any constraints. The other methods do not perform this check. The two-sided method is in principle more precise, but requires twice as many function evaluations. Default: 0. .MPMINSTEP - the minimum change to be made in the parameter value. During the fitting process, the parameter will be changed by multiples of this value. The actual step is computed as: DELTA1 = MPMINSTEP*ROUND(DELTA0/MPMINSTEP) where DELTA0 and DELTA1 are the estimated parameter changes before and after this constraint is applied. Note that this constraint should be used with care since it may cause non-converging, oscillating solutions. A value of 0 indicates no minimum. Default: 0. .MPMAXSTEP - the maximum change to be made in the parameter value. During the fitting process, the parameter will never be changed by more than this value. A value of 0 indicates no maximum. Default: 0. .TIED - a string expression which "ties" the parameter to other free or fixed parameters. Any expression involving constants and the parameter array P are permitted. Example: if parameter 2 is always to be twice parameter 1 then use the following: parinfo[2].tied = '2 * P[1]'. Since they are totally constrained, tied parameters are considered to be fixed; no errors are computed for them. [ NOTE: the PARNAME can't be used in expressions. ] Future modifications to the PARINFO structure, if any, will involve adding structure tags beginning with the two letters "MP". Therefore programmers are urged to avoid using tags starting with the same letters; otherwise they are free to include their own fields within the PARINFO structure, and they will be ignored. PARINFO Example: parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $ limits:[0.D,0]}, 5) parinfo[0].fixed = 1 parinfo[4].limited[0] = 1 parinfo[4].limits[0] = 50.D parinfo[*].value = [5.7D, 2.2, 500., 1.5, 2000.] A total of 5 parameters, with starting values of 5.7, 2.2, 500, 1.5, and 2000 are given. The first parameter is fixed at a value of 5.7, and the last parameter is constrained to be above 50. INPUTS: MYFUNCT - a string variable containing an IDL expression. The only restriction is that the independent variable *must* be referred to as "X" and model parameters *must* be referred to as an array called "P". Do not use symbol names beginning with the underscore, "_". The expression should calculate "model" Y values given the X values and model parameters. Using the vector notation of IDL, this can be quite easy to do. If your expression gets complicated, you may wish to make an IDL function which will improve performance and readibility. The resulting array should be of the same size and dimensions as the "measured" Y values. X - Array of independent variable values. Y - Array of "measured" dependent variable values. Y should have the same data type as X. The function MYFUNCT should map X->Y. ERR - Array of "measured" 1-sigma uncertainties. ERR should have the same data type as Y. ERR is ignored if the WEIGHTS keyword is specified. START_PARAMS - An array of starting values for each of the parameters of the model. The number of parameters should be fewer than the number of measurements. Also, the parameters should have the same data type as the measurements (double is preferred). This parameter is optional if the PARINFO keyword is used (see MPFIT). The PARINFO keyword provides a mechanism to fix or constrain individual parameters. If both START_PARAMS and PARINFO are passed, then the starting *value* is taken from START_PARAMS, but the *constraints* are taken from PARINFO. If no parameters are given, then MPFITEXPR attempts to determine the number of parameters by scanning the expression. Parameters determined this way are initialized to zero. This technique is not fully reliable, so users are advised to pass explicit parameter starting values. RETURNS: Returns the array of best-fit parameters. KEYWORD PARAMETERS: BESTNORM - the value of the summed, squared, weighted residuals for the returned parameter values, i.e. the chi-square value. COVAR - the covariance matrix for the set of parameters returned by MPFIT. The matrix is NxN where N is the number of parameters. The square root of the diagonal elements gives the formal 1-sigma statistical errors on the parameters IF errors were treated "properly" in MYFUNC. Parameter errors are also returned in PERROR. To compute the correlation matrix, PCOR, use this: IDL> PCOR = COV * 0 IDL> FOR i = 0, n-1 DO FOR j = 0, n-1 DO $ PCOR[i,j] = COV[i,j]/sqrt(COV[i,i]*COV[j,j]) If NOCOVAR is set or MPFIT terminated abnormally, then COVAR is set to a scalar with value !VALUES.D_NAN. DOF - number of degrees of freedom, computed as DOF = N_ELEMENTS(DEVIATES) - NFREE Note that this doesn't account for pegged parameters (see NPEGGED). ERRMSG - a string error or warning message is returned. FTOL - a nonnegative input variable. Termination occurs when both the actual and predicted relative reductions in the sum of squares are at most FTOL (and STATUS is accordingly set to 1 or 3). Therefore, FTOL measures the relative error desired in the sum of squares. Default: 1D-10 FUNCTARGS - passed directly to the expression as the variable PRIVATE. Any user-private data can be communicated to the user expression using this keyword. Default: PRIVATE is undefined in user expression GTOL - a nonnegative input variable. Termination occurs when the cosine of the angle between fvec and any column of the jacobian is at most GTOL in absolute value (and STATUS is accordingly set to 4). Therefore, GTOL measures the orthogonality desired between the function vector and the columns of the jacobian. Default: 1D-10 ITERARGS - The keyword arguments to be passed to ITERPROC via the _EXTRA mechanism. This should be a structure, and is similar in operation to FUNCTARGS. Default: no arguments are passed. ITERPROC - The name of a procedure to be called upon each NPRINT iteration of the MPFIT routine. It should be declared in the following way: PRO ITERPROC, MYFUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $ PARINFO=parinfo, QUIET=quiet, ... ; perform custom iteration update END ITERPROC must either accept all three keyword parameters (FUNCTARGS, PARINFO and QUIET), or at least accept them via the _EXTRA keyword. MYFUNCT is the user-supplied function to be minimized, P is the current set of model parameters, ITER is the iteration number, and FUNCTARGS are the arguments to be passed to MYFUNCT. FNORM should be the chi-squared value. QUIET is set when no textual output should be printed. See below for documentation of PARINFO. In implementation, ITERPROC can perform updates to the terminal or graphical user interface, to provide feedback while the fit proceeds. If the fit is to be stopped for any reason, then ITERPROC should set the common block variable ERROR_CODE to negative value (see MPFIT_ERROR common block below). In principle, ITERPROC should probably not modify the parameter values, because it may interfere with the algorithm's stability. In practice it is allowed. Default: an internal routine is used to print the parameter values. MAXITER - The maximum number of iterations to perform. If the number is exceeded, then the STATUS value is set to 5 and MPFIT returns. Default: 200 iterations NFEV - the number of MYFUNCT function evaluations performed. NFREE - the number of free parameters in the fit. This includes parameters which are not FIXED and not TIED, but it does include parameters which are pegged at LIMITS. NITER - the number of iterations completed. NOCOVAR - set this keyword to prevent the calculation of the covariance matrix before returning (see COVAR) NPEGGED - the number of free parameters which are pegged at a LIMIT. NPRINT - The frequency with which ITERPROC is called. A value of 1 indicates that ITERPROC is called with every iteration, while 2 indicates every other iteration, etc. Note that several Levenberg-Marquardt attempts can be made in a single iteration. Default value: 1 PARINFO - Provides a mechanism for more sophisticated constraints to be placed on parameter values. When PARINFO is not passed, then it is assumed that all parameters are free and unconstrained. Values in PARINFO are never modified during a call to MPFIT. See description above for the structure of PARINFO. Default value: all parameters are free and unconstrained. PERROR - The formal 1-sigma errors in each parameter, computed from the covariance matrix. If a parameter is held fixed, or if it touches a boundary, then the error is reported as zero. If the fit is unweighted (i.e. no errors were given, or the weights were uniformly set to unity), then PERROR will probably not represent the true parameter uncertainties. *If* you can assume that the true reduced chi-squared value is unity -- meaning that the fit is implicitly assumed to be of good quality -- then the estimated parameter uncertainties can be computed by scaling PERROR by the measured chi-squared value. DOF = N_ELEMENTS(X) - N_ELEMENTS(PARMS) ; deg of freedom PCERROR = PERROR * SQRT(BESTNORM / DOF) ; scaled uncertainties QUIET - set this keyword when no textual output should be printed by MPFIT STATUS - an integer status code is returned. All values other than zero can represent success. It can have one of the following values: 0 improper input parameters. 1 both actual and predicted relative reductions in the sum of squares are at most FTOL. 2 relative error between two consecutive iterates is at most XTOL 3 conditions for STATUS = 1 and STATUS = 2 both hold. 4 the cosine of the angle between fvec and any column of the jacobian is at most GTOL in absolute value. 5 the maximum number of iterations has been reached 6 FTOL is too small. no further reduction in the sum of squares is possible. 7 XTOL is too small. no further improvement in the approximate solution x is possible. 8 GTOL is too small. fvec is orthogonal to the columns of the jacobian to machine precision. WEIGHTS - Array of weights to be used in calculating the chi-squared value. If WEIGHTS is specified then the ERR parameter is ignored. The chi-squared value is computed as follows: CHISQ = TOTAL( (Y-MYFUNCT(X,P))^2 * ABS(WEIGHTS) ) Here are common values of WEIGHTS: 1D/ERR^2 - Normal weighting (ERR is the measurement error) 1D/Y - Poisson weighting (counting statistics) 1D - Unweighted XTOL - a nonnegative input variable. Termination occurs when the relative error between two consecutive iterates is at most XTOL (and STATUS is accordingly set to 2 or 3). Therefore, XTOL measures the relative error desired in the approximate solution. Default: 1D-10 YFIT - the best-fit model function, as returned by MYFUNCT. EXAMPLE: ; First, generate some synthetic data x = dindgen(200) * 0.1 - 10. ; Independent variable yi = gauss1(x, [2.2D, 1.4, 3000.]) + 1000 ; "Ideal" Y variable y = yi + randomn(seed, 200) * sqrt(yi) ; Measured, w/ noise sy = sqrt(y) ; Poisson errors ; Now fit a Gaussian to see how well we can recover expr = 'P[0] + GAUSS1(X, P[1:3])' ; fitting function p0 = [800.D, 1.D, 1., 500.] ; Initial guess p = mpfitexpr(expr, x, y, sy, p0) ; Fit the expression print, p plot, x, y ; Plot data oplot, x, mpevalexpr(expr, x, p) ; Plot model Generates a synthetic data set with a Gaussian peak, and Poisson statistical uncertainty. Then a model consisting of a constant plus Gaussian is fit to the data. COMMON BLOCKS: COMMON MPFIT_ERROR, ERROR_CODE User routines may stop the fitting process at any time by setting an error condition. This condition may be set in either the user's model computation routine (MYFUNCT), or in the iteration procedure (ITERPROC). To stop the fitting, the above common block must be declared, and ERROR_CODE must be set to a negative number. After the user procedure or function returns, MPFIT checks the value of this common block variable and exits immediately if the error condition has been set. By default the value of ERROR_CODE is zero, indicating a successful function/procedure call. REFERENCES: MINPACK-1, Jorge More', available from netlib (www.netlib.org). "Optimization Software Guide," Jorge More' and Stephen Wright, SIAM, *Frontiers in Applied Mathematics*, Number 14. MODIFICATION HISTORY: Written, Apr-Jul 1998, CM Added PERROR keyword, 04 Aug 1998, CM Added COVAR keyword, 20 Aug 1998, CM Added NITER output keyword, 05 Oct 1998 D.L Windt, Bell Labs, windt@bell-labs.com; Added ability to return model function in YFIT, 09 Nov 1998 Parameter values can be tied to others, 09 Nov 1998 Added MPEVALEXPR utility function, 09 Dec 1998 Cosmetic documentation updates, 16 Apr 1999, CM More cosmetic documentation updates, 14 May 1999, CM Made sure to update STATUS value, 25 Sep 1999, CM Added WEIGHTS keyword, 25 Sep 1999, CM Changed from handles to common blocks, 25 Sep 1999, CM - commons seem much cleaner and more logical in this case. Alphabetized documented keywords, 02 Oct 1999, CM Added QUERY keyword and query checking of MPFIT, 29 Oct 1999, CM Check to be sure that X and Y are present, 02 Nov 1999, CM Documented PERROR for unweighted fits, 03 Nov 1999, CM Removed ITERPROC='' when quiet EQ 1, 10 Jan 2000, CM Changed to ERROR_CODE for error condition, 28 Jan 2000, CM Updated the EXAMPLE, 26 Mar 2000, CM Copying permission terms have been liberalized, 26 Mar 2000, CM Propagated improvements from MPFIT, 17 Dec 2000, CM Correct reference to _WTS in MPFITEXPR_EVAL, 25 May 2001, CM (thanks to Mark Fardal) Added useful FUNCTARGS behavior (as yet undocumented), 04 Jul 2002, CM Documented FUNCTARGS/PRIVATE behavior, 22 Jul 2002, CM Added NFREE and NPEGGED keywords, 13 Sep 2002, CM Documented RELSTEP field of PARINFO (!!), CM, 25 Oct 2002 Add DOF keyword, CM, 31 Jul 2003 Add FUNCTARGS keyword to MPEVALEXPR, CM 08 Aug 2003 Minor documentation adjustment, 03 Feb 2004, CM $Id: mpfitexpr.pro,v 2.2 2009/04/15 14:43:47 bdavis Exp $
Source: src/idl_cvs/mpfitfun.pro
src/idl_cvs/mpfitfun.pro
NAME: MPFITFUN AUTHOR: Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 craigm@lheamail.gsfc.nasa.gov UPDATED VERSIONs can be found on my WEB PAGE: http://cow.physics.wisc.edu/~craigm/idl/idl.html PURPOSE: Perform Levenberg-Marquardt least-squares fit to IDL function MAJOR TOPICS: Curve and Surface Fitting CALLING SEQUENCE: parms = MPFITFUN(MYFUNCT, X, Y, ERR, start_params, ...) DESCRIPTION: MPFITFUN fits a user-supplied model -- in the form of an IDL function -- to a set of user-supplied data. MPFITFUN calls MPFIT, the MINPACK-1 least-squares minimizer, to do the main work. Given the data and their uncertainties, MPFITFUN finds the best set of model parameters which match the data (in a least-squares sense) and returns them in an array. The user must supply the following items: - An array of independent variable values ("X"). - An array of "measured" *dependent* variable values ("Y"). - An array of "measured" 1-sigma uncertainty values ("ERR"). - The name of an IDL function which computes Y given X ("MYFUNCT"). - Starting guesses for all of the parameters ("START_PARAMS"). There are very few restrictions placed on X, Y or MYFUNCT. Simply put, MYFUNCT must map the "X" values into "Y" values given the model parameters. The "X" values may represent any independent variable (not just Cartesian X), and indeed may be multidimensional themselves. For example, in the application of image fitting, X may be a 2xN array of image positions. MPFITFUN carefully avoids passing large arrays where possible to improve performance. See below for an example of usage. USER FUNCTION The user must define a function which returns the model value. For applications which use finite-difference derivatives -- the default -- the user function should be declared in the following way: FUNCTION MYFUNCT, X, P ; The independent variable is X ; Parameter values are passed in "P" YMOD = ... computed model values at X ... return, YMOD END The returned array YMOD must have the same dimensions and type as the "measured" Y values. User functions may also indicate a fatal error condition using the ERROR_CODE common block variable, as described below under the MPFIT_ERROR common block definition. See the discussion under "ANALYTIC DERIVATIVES" and AUTODERIVATIVE in MPFIT.PRO if you wish to compute the derivatives for yourself. AUTODERIVATIVE is accepted by MPFITFUN and passed directly to MPFIT. The user function must accept one additional parameter, DP, which contains the derivative of the user function with respect to each parameter at each data point, as described in MPFIT.PRO. CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD The behavior of MPFIT can be modified with respect to each parameter to be fitted. A parameter value can be fixed; simple boundary constraints can be imposed; limitations on the parameter changes can be imposed; properties of the automatic derivative can be modified; and parameters can be tied to one another. These properties are governed by the PARINFO structure, which is passed as a keyword parameter to MPFIT. PARINFO should be an array of structures, one for each parameter. Each parameter is associated with one element of the array, in numerical order. The structure can have the following entries (none are required): .VALUE - the starting parameter value (but see the START_PARAMS parameter for more information). .FIXED - a boolean value, whether the parameter is to be held fixed or not. Fixed parameters are not varied by MPFIT, but are passed on to MYFUNCT for evaluation. .LIMITED - a two-element boolean array. If the first/second element is set, then the parameter is bounded on the lower/upper side. A parameter can be bounded on both sides. Both LIMITED and LIMITS must be given together. .LIMITS - a two-element float or double array. Gives the parameter limits on the lower and upper sides, respectively. Zero, one or two of these values can be set, depending on the values of LIMITED. Both LIMITED and LIMITS must be given together. .PARNAME - a string, giving the name of the parameter. The fitting code of MPFIT does not use this tag in any way. However, the default ITERPROC will print the parameter name if available. .STEP - the step size to be used in calculating the numerical derivatives. If set to zero, then the step size is computed automatically. Ignored when AUTODERIVATIVE=0. This value is superceded by the RELSTEP value. .RELSTEP - the *relative* step size to be used in calculating the numerical derivatives. This number is the fractional size of the step, compared to the parameter value. This value supercedes the STEP setting. If the parameter is zero, then a default step size is chosen. .MPSIDE - the sidedness of the finite difference when computing numerical derivatives. This field can take four values: 0 - one-sided derivative computed automatically 1 - one-sided derivative (f(x+h) - f(x) )/h -1 - one-sided derivative (f(x) - f(x-h))/h 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h) Where H is the STEP parameter described above. The "automatic" one-sided derivative method will chose a direction for the finite difference which does not violate any constraints. The other methods do not perform this check. The two-sided method is in principle more precise, but requires twice as many function evaluations. Default: 0. .MPMINSTEP - the minimum change to be made in the parameter value. During the fitting process, the parameter will be changed by multiples of this value. The actual step is computed as: DELTA1 = MPMINSTEP*ROUND(DELTA0/MPMINSTEP) where DELTA0 and DELTA1 are the estimated parameter changes before and after this constraint is applied. Note that this constraint should be used with care since it may cause non-converging, oscillating solutions. A value of 0 indicates no minimum. Default: 0. .MPMAXSTEP - the maximum change to be made in the parameter value. During the fitting process, the parameter will never be changed by more than this value. A value of 0 indicates no maximum. Default: 0. .TIED - a string expression which "ties" the parameter to other free or fixed parameters. Any expression involving constants and the parameter array P are permitted. Example: if parameter 2 is always to be twice parameter 1 then use the following: parinfo[2].tied = '2 * P[1]'. Since they are totally constrained, tied parameters are considered to be fixed; no errors are computed for them. [ NOTE: the PARNAME can't be used in expressions. ] Future modifications to the PARINFO structure, if any, will involve adding structure tags beginning with the two letters "MP". Therefore programmers are urged to avoid using tags starting with the same letters; otherwise they are free to include their own fields within the PARINFO structure, and they will be ignored. PARINFO Example: parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $ limits:[0.D,0]}, 5) parinfo[0].fixed = 1 parinfo[4].limited[0] = 1 parinfo[4].limits[0] = 50.D parinfo[*].value = [5.7D, 2.2, 500., 1.5, 2000.] A total of 5 parameters, with starting values of 5.7, 2.2, 500, 1.5, and 2000 are given. The first parameter is fixed at a value of 5.7, and the last parameter is constrained to be above 50. INPUTS: MYFUNCT - a string variable containing the name of an IDL function. This function computes the "model" Y values given the X values and model parameters, as desribed above. X - Array of independent variable values. Y - Array of "measured" dependent variable values. Y should have the same data type as X. The function MYFUNCT should map X->Y. ERR - Array of "measured" 1-sigma uncertainties. ERR should have the same data type as Y. ERR is ignored if the WEIGHTS keyword is specified. START_PARAMS - An array of starting values for each of the parameters of the model. The number of parameters should be fewer than the number of measurements. Also, the parameters should have the same data type as the measurements (double is preferred). This parameter is optional if the PARINFO keyword is used (see MPFIT). The PARINFO keyword provides a mechanism to fix or constrain individual parameters. If both START_PARAMS and PARINFO are passed, then the starting *value* is taken from START_PARAMS, but the *constraints* are taken from PARINFO. RETURNS: Returns the array of best-fit parameters. KEYWORD PARAMETERS: BESTNORM - the value of the summed squared residuals for the returned parameter values. COVAR - the covariance matrix for the set of parameters returned by MPFIT. The matrix is NxN where N is the number of parameters. The square root of the diagonal elements gives the formal 1-sigma statistical errors on the parameters IF errors were treated "properly" in MYFUNC. Parameter errors are also returned in PERROR. To compute the correlation matrix, PCOR, use this: IDL> PCOR = COV * 0 IDL> FOR i = 0, n-1 DO FOR j = 0, n-1 DO $ PCOR[i,j] = COV[i,j]/sqrt(COV[i,i]*COV[j,j]) If NOCOVAR is set or MPFIT terminated abnormally, then COVAR is set to a scalar with value !VALUES.D_NAN. CASH - when set, the fit statistic is changed to a derivative of the CASH statistic. The model function must be strictly positive. DOF - number of degrees of freedom, computed as DOF = N_ELEMENTS(DEVIATES) - NFREE Note that this doesn't account for pegged parameters (see NPEGGED). ERRMSG - a string error or warning message is returned. FTOL - a nonnegative input variable. Termination occurs when both the actual and predicted relative reductions in the sum of squares are at most FTOL (and STATUS is accordingly set to 1 or 3). Therefore, FTOL measures the relative error desired in the sum of squares. Default: 1D-10 FUNCTARGS - A structure which contains the parameters to be passed to the user-supplied function specified by MYFUNCT via the _EXTRA mechanism. This is the way you can pass additional data to your user-supplied function without using common blocks. By default, no extra parameters are passed to the user-supplied function. GTOL - a nonnegative input variable. Termination occurs when the cosine of the angle between fvec and any column of the jacobian is at most GTOL in absolute value (and STATUS is accordingly set to 4). Therefore, GTOL measures the orthogonality desired between the function vector and the columns of the jacobian. Default: 1D-10 ITERARGS - The keyword arguments to be passed to ITERPROC via the _EXTRA mechanism. This should be a structure, and is similar in operation to FUNCTARGS. Default: no arguments are passed. ITERPROC - The name of a procedure to be called upon each NPRINT iteration of the MPFIT routine. It should be declared in the following way: PRO ITERPROC, MYFUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $ PARINFO=parinfo, QUIET=quiet, ... ; perform custom iteration update END ITERPROC must either accept all three keyword parameters (FUNCTARGS, PARINFO and QUIET), or at least accept them via the _EXTRA keyword. MYFUNCT is the user-supplied function to be minimized, P is the current set of model parameters, ITER is the iteration number, and FUNCTARGS are the arguments to be passed to MYFUNCT. FNORM should be the chi-squared value. QUIET is set when no textual output should be printed. See below for documentation of PARINFO. In implementation, ITERPROC can perform updates to the terminal or graphical user interface, to provide feedback while the fit proceeds. If the fit is to be stopped for any reason, then ITERPROC should set the common block variable ERROR_CODE to negative value (see MPFIT_ERROR common block below). In principle, ITERPROC should probably not modify the parameter values, because it may interfere with the algorithm's stability. In practice it is allowed. Default: an internal routine is used to print the parameter values. MAXITER - The maximum number of iterations to perform. If the number is exceeded, then the STATUS value is set to 5 and MPFIT returns. Default: 200 iterations NFEV - the number of MYFUNCT function evaluations performed. NFREE - the number of free parameters in the fit. This includes parameters which are not FIXED and not TIED, but it does include parameters which are pegged at LIMITS. NITER - the number of iterations completed. NOCOVAR - set this keyword to prevent the calculation of the covariance matrix before returning (see COVAR) NPEGGED - the number of free parameters which are pegged at a LIMIT. NPRINT - The frequency with which ITERPROC is called. A value of 1 indicates that ITERPROC is called with every iteration, while 2 indicates every other iteration, etc. Note that several Levenberg-Marquardt attempts can be made in a single iteration. Default value: 1 PARINFO - Provides a mechanism for more sophisticated constraints to be placed on parameter values. When PARINFO is not passed, then it is assumed that all parameters are free and unconstrained. Values in PARINFO are never modified during a call to MPFIT. See description above for the structure of PARINFO. Default value: all parameters are free and unconstrained. PERROR - The formal 1-sigma errors in each parameter, computed from the covariance matrix. If a parameter is held fixed, or if it touches a boundary, then the error is reported as zero. If the fit is unweighted (i.e. no errors were given, or the weights were uniformly set to unity), then PERROR will probably not represent the true parameter uncertainties. *If* you can assume that the true reduced chi-squared value is unity -- meaning that the fit is implicitly assumed to be of good quality -- then the estimated parameter uncertainties can be computed by scaling PERROR by the measured chi-squared value. DOF = N_ELEMENTS(X) - N_ELEMENTS(PARMS) ; deg of freedom PCERROR = PERROR * SQRT(BESTNORM / DOF) ; scaled uncertainties QUIET - set this keyword when no textual output should be printed by MPFIT STATUS - an integer status code is returned. All values other than zero can represent success. It can have one of the following values: 0 improper input parameters. 1 both actual and predicted relative reductions in the sum of squares are at most FTOL. 2 relative error between two consecutive iterates is at most XTOL 3 conditions for STATUS = 1 and STATUS = 2 both hold. 4 the cosine of the angle between fvec and any column of the jacobian is at most GTOL in absolute value. 5 the maximum number of iterations has been reached 6 FTOL is too small. no further reduction in the sum of squares is possible. 7 XTOL is too small. no further improvement in the approximate solution x is possible. 8 GTOL is too small. fvec is orthogonal to the columns of the jacobian to machine precision. WEIGHTS - Array of weights to be used in calculating the chi-squared value. If WEIGHTS is specified then the ERR parameter is ignored. The chi-squared value is computed as follows: CHISQ = TOTAL( (Y-MYFUNCT(X,P))^2 * ABS(WEIGHTS) ) Here are common values of WEIGHTS: 1D/ERR^2 - Normal weighting (ERR is the measurement error) 1D/Y - Poisson weighting (counting statistics) 1D - Unweighted XTOL - a nonnegative input variable. Termination occurs when the relative error between two consecutive iterates is at most XTOL (and STATUS is accordingly set to 2 or 3). Therefore, XTOL measures the relative error desired in the approximate solution. Default: 1D-10 YFIT - the best-fit model function, as returned by MYFUNCT. EXAMPLE: ; First, generate some synthetic data npts = 200 x = dindgen(npts) * 0.1 - 10. ; Independent variable yi = gauss1(x, [2.2D, 1.4, 3000.]) ; "Ideal" Y variable y = yi + randomn(seed, npts) * sqrt(1000. + yi); Measured, w/ noise sy = sqrt(1000.D + y) ; Poisson errors ; Now fit a Gaussian to see how well we can recover p0 = [1.D, 1., 1000.] ; Initial guess (cent, width, area) p = mpfitfun('GAUSS1', x, y, sy, p0) ; Fit a function print, p Generates a synthetic data set with a Gaussian peak, and Poisson statistical uncertainty. Then the same function is fitted to the data (with different starting parameters) to see how close we can get. COMMON BLOCKS: COMMON MPFIT_ERROR, ERROR_CODE User routines may stop the fitting process at any time by setting an error condition. This condition may be set in either the user's model computation routine (MYFUNCT), or in the iteration procedure (ITERPROC). To stop the fitting, the above common block must be declared, and ERROR_CODE must be set to a negative number. After the user procedure or function returns, MPFIT checks the value of this common block variable and exits immediately if the error condition has been set. By default the value of ERROR_CODE is zero, indicating a successful function/procedure call. REFERENCES: MINPACK-1, Jorge More', available from netlib (www.netlib.org). "Optimization Software Guide," Jorge More' and Stephen Wright, SIAM, *Frontiers in Applied Mathematics*, Number 14. MODIFICATION HISTORY: Written, Apr-Jul 1998, CM Added PERROR keyword, 04 Aug 1998, CM Added COVAR keyword, 20 Aug 1998, CM Added ITER output keyword, 05 Oct 1998 D.L Windt, Bell Labs, windt@bell-labs.com; Added ability to return model function in YFIT, 09 Nov 1998 Analytical derivatives allowed via AUTODERIVATIVE keyword, 09 Nov 1998 Parameter values can be tied to others, 09 Nov 1998 Cosmetic documentation updates, 16 Apr 1999, CM More cosmetic documentation updates, 14 May 1999, CM Made sure to update STATUS, 25 Sep 1999, CM Added WEIGHTS keyword, 25 Sep 1999, CM Changed from handles to common blocks, 25 Sep 1999, CM - commons seem much cleaner and more logical in this case. Alphabetized documented keywords, 02 Oct 1999, CM Added QUERY keyword and query checking of MPFIT, 29 Oct 1999, CM Corrected EXAMPLE (offset of 1000), 30 Oct 1999, CM Check to be sure that X and Y are present, 02 Nov 1999, CM Documented PERROR for unweighted fits, 03 Nov 1999, CM Changed to ERROR_CODE for error condition, 28 Jan 2000, CM Corrected errors in EXAMPLE, 26 Mar 2000, CM Copying permission terms have been liberalized, 26 Mar 2000, CM Propagated improvements from MPFIT, 17 Dec 2000, CM Added CASH statistic, 10 Jan 2001 Added NFREE and NPEGGED keywords, 11 Sep 2002, CM Documented RELSTEP field of PARINFO (!!), CM, 25 Oct 2002 Add DOF keyword to return degrees of freedom, CM, 23 June 2003 Convert to IDL 5 array syntax (!), 16 Jul 2006, CM $Id: mpfitfun.pro,v 2.3 2009/04/15 14:43:47 bdavis Exp $
Source: src/idl_cvs/mpfitpeak.pro
src/idl_cvs/mpfitpeak.pro
NAME: MPFITPEAK AUTHOR: Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 craigm@lheamail.gsfc.nasa.gov UPDATED VERSIONs can be found on my WEB PAGE: http://cow.physics.wisc.edu/~craigm/idl/idl.html PURPOSE: Fit a gaussian, lorentzian or Moffat model to data MAJOR TOPICS: Curve and Surface Fitting CALLING SEQUENCE: yfit = MPFITPEAK(X, Y, A, NTERMS=nterms, ...) DESCRIPTION: MPFITPEAK fits a gaussian, lorentzian or Moffat model using the non-linear least squares fitter MPFIT. MPFITPEAK is meant to be a drop-in replacement for IDL's GAUSSFIT function (and requires MPFIT and MPFITFUN). The choice of the fitting function is determined by the keywords GAUSSIAN, LORENTZIAN and MOFFAT. By default the gaussian model function is used. [ The Moffat function is a modified Lorentzian with variable power law index. (Moffat, A. F. J. 1969, Astronomy & Astrophysics, v. 3, p. 455-461) ] The functional form of the baseline is determined by NTERMS and the function to be fitted. NTERMS represents the total number of parameters, A, to be fitted. The functional forms and the meanings of the parameters are described in this table: GAUSSIAN# Lorentzian# Moffat# Model A[0]*exp(-0.5*u^2) A[0]/(u^2 + 1) A[0]/(u^2 + 1)^A[3] A[0] Peak Value Peak Value Peak Value A[1] Peak Centroid Peak Centroid Peak Centroid A[2] Gaussian Sigma HWHM% HWHM% A[3] + A[3] * + A[3] * Moffat Index A[4] + A[4]*x * + A[4]*x * + A[4] * A[5] + A[5]*x * Notes: # u = (x - A[1])/A[2] % Half-width at half maximum * Optional depending on NTERMS By default the initial starting values for the parameters A are estimated from the data. However, explicit starting values can be supplied using the ESTIMATES keyword. Also, error or weighting values can optionally be provided; otherwise the fit is unweighted. MPFITPEAK fits the peak value of the curve. The area under a gaussian peak is A[0]*A[2]*SQRT(2*!DPI); the area under a lorentzian peak is A[0]*A[2]*!DPI. RESTRICTIONS: If no starting parameter ESTIMATES are provided, then MPFITPEAK attempts to estimate them from the data. This is not a perfect science; however, the author believes that the technique implemented here is more robust than the one used in IDL's GAUSSFIT. The author has tested cases of strong peaks, noisy peaks and broad peaks, all with success. Users should be aware that if the baseline term contains a strong linear component then the automatic estimation may fail. For automatic estimation to work the peak amplitude should dominate over the the maximum baseline. INPUTS: X - Array of independent variable values, whose values should monotonically increase. Y - Array of "measured" dependent variable values. Y should have the same data type and dimension as X. OUTPUTS: A - Upon return, an array of NTERMS best fit parameter values. See the table above for the meanings of each parameter element. RETURNS: Returns the best fitting model function. KEYWORDS: ** NOTE ** Additional keywords such as PARINFO, BESTNORM, and STATUS are accepted by MPFITPEAK but not documented here. Please see the documentation for MPFIT for the description of these advanced options. CHISQ - the value of the summed squared residuals for the returned parameter values. DOF - number of degrees of freedom, computed as DOF = N_ELEMENTS(DEVIATES) - NFREE Note that this doesn't account for pegged parameters (see NPEGGED). ERROR - upon input, the measured 1-sigma uncertainties in the "Y" values. If no ERROR or WEIGHTS are given, then the fit is unweighted. ESTIMATES - Array of starting values for each parameter of the model. The number of parameters should at least be three (four for Moffat), and if less than NTERMS, will be extended with zeroes. Default: parameter values are estimated from data. GAUSSIAN - if set, fit a gaussian model function. The Default. LORENTZIAN - if set, fit a lorentzian model function. MOFFAT - if set, fit a Moffat model function. MEASURE_ERRORS - synonym for ERRORS, for consistency with built-in IDL fitting routines. NEGATIVE / POSITIVE - if set, and ESTIMATES is not provided, then MPFITPEAK will assume that a negative/positive peak is present. Default: determined automatically NFREE - the number of free parameters in the fit. This includes parameters which are not FIXED and not TIED, but it does include parameters which are pegged at LIMITS. NTERMS - An integer describing the number of fitting terms. NTERMS must have a minimum value, but can optionally be larger depending on the desired baseline. For gaussian and lorentzian models, NTERMS must be three (zero baseline), four (constant baseline) or five (linear baseline). Default: 4 For the Moffat model, NTERMS must be four (zero baseline), five (constant baseline), or six (linear baseline). Default: 5 PERROR - upon return, the 1-sigma uncertainties of the parameter values A. These values are only meaningful if the ERRORS or WEIGHTS keywords are specified properly. If the fit is unweighted (i.e. no errors were given, or the weights were uniformly set to unity), then PERROR will probably not represent the true parameter uncertainties. *If* you can assume that the true reduced chi-squared value is unity -- meaning that the fit is implicitly assumed to be of good quality -- then the estimated parameter uncertainties can be computed by scaling PERROR by the measured chi-squared value. DOF = N_ELEMENTS(X) - N_ELEMENTS(PARMS) ; deg of freedom PCERROR = PERROR * SQRT(BESTNORM / DOF) ; scaled uncertainties QUIET - if set then diagnostic fitting messages are suppressed. Default: QUIET=1 (i.e., no diagnostics) SIGMA - synonym for PERROR (1-sigma parameter uncertainties), for compatibility with GAUSSFIT. Do not confuse this with the Gaussian "sigma" width parameter. WEIGHTS - Array of weights to be used in calculating the chi-squared value. If WEIGHTS is specified then the ERROR keyword is ignored. The chi-squared value is computed as follows: CHISQ = TOTAL( (Y-MYFUNCT(X,P))^2 * ABS(WEIGHTS) ) Here are common values of WEIGHTS: 1D/ERR^2 - Normal weighting (ERR is the measurement error) 1D/Y - Poisson weighting (counting statistics) 1D - Unweighted The ERROR keyword takes precedence over any WEIGHTS keyword values. If no ERROR or WEIGHTS are given, then the fit is unweighted. YERROR - upon return, the root-mean-square variance of the residuals. EXAMPLE: ; First, generate some synthetic data npts = 200 x = dindgen(npts) * 0.1 - 10. ; Independent variable yi = gauss1(x, [2.2D, 1.4, 3000.]) + 1000 ; "Ideal" Y variable y = yi + randomn(seed, npts) * sqrt(1000. + yi); Measured, w/ noise sy = sqrt(1000.D + y) ; Poisson errors ; Now fit a Gaussian to see how well we can recover the original yfit = mpfitpeak(x, y, a, error=sy) print, p Generates a synthetic data set with a Gaussian peak, and Poisson statistical uncertainty. Then the same function is fitted to the data. REFERENCES: MINPACK-1, Jorge More', available from netlib (www.netlib.org). "Optimization Software Guide," Jorge More' and Stephen Wright, SIAM, *Frontiers in Applied Mathematics*, Number 14. MODIFICATION HISTORY: New algorithm for estimating starting values, CM, 31 Oct 1999 Documented, 02 Nov 1999 Small documentation fixes, 02 Nov 1999 Slight correction to calculation of dx, CM, 02 Nov 1999 Documented PERROR for unweighted fits, 03 Nov 1999, CM Copying permission terms have been liberalized, 26 Mar 2000, CM Change requirements on # elements in X and Y, 20 Jul 2000, CM (thanks to David Schlegel ) Added documentation on area under curve, 29 Aug 2000, CM Added POSITIVE and NEGATIVE keywords, 17 Nov 2000, CM Added reference to Moffat paper, 10 Jan 2001, CM Added usage message, 26 Jul 2001, CM Documentation clarification, 05 Sep 2001, CM Make more consistent with comparable IDL routines, 30 Jun 2003, CM Assumption of sorted data was removed, CM, 06 Sep 2003, CM Add some defensive code against divide by zero, 30 Nov 2005, CM Add some defensive code against all Y values equal to each other, 17 Apr 2005, CM Convert to IDL 5 array syntax (!), 16 Jul 2006, CM $Id: mpfitpeak.pro,v 2.2 2009/04/15 14:43:47 bdavis Exp $
Source: src/idl_cvs/mplottitle.pro
src/idl_cvs/mplottitle.pro
NAME: mplottitle PURPOSE: write a plot title with NSTX label and a list of shots on left CATEGORY: Plotting CALLING SEQUENCE: IDL> mplottitle, screenTitle, shots, HLP=hlp, $ COLORS=colors, LOGO=logo, $ CHARMAG=charMag, TopMargin=TopMargin, RIGHT=right, $ INSIDE=inside, XPOS=xPos INPUTS: screenTitle = whatever you want as the main plot title. shots = array of shot numbers. KEYWORD PARAMETERS: HLP - When set, help information is printed. colors - color indices of each shot number (default = !p.color) logo - Characters to draw above shot list. Default = '=NSTX= TopMargin - in characters, extra space to shift down the characters right - if set, labels put on left inside - if set, lists shots inside plot box xPos - xPos, in Device Coordinates, for left side of shot list (overrides /RIGHT & /INSIDE) OUTPUTS: none COMMON BLOCKS: NONE EXAMPLES: IDL> !Y.OMARGIN=[0,1] IDL> plot,indgen(100) IDL> mplottitle, 'This is my title', [100523,100544] IDL> plot,indgen(100) IDL> cs=mk_color(/load, /quiet) IDL> colors=[cs.black, cs.red, cs.green, cs.blue, cs.magenta] IDL> mplottitle,' ', shots, colors=colors, logo='Shots:', $ charMag=1.35, topMargin=0.5 NOTES: You'll probably want a little more room at the top of your plots, e.g. IDL> !Y.OMARGIN=[0,1] MODIFICATION HISTORY: 19-Jan-06 Added psym keyword 01-Oct-01 Added keywords for inside & right [BD] 21-Sep-99 Written by Bill Davis, PPPL
Source: src/idl_cvs/mptscam.pro
src/idl_cvs/mptscam.pro
NAME: mptscam PURPOSE: Animate several plots, including MPTS variables CATEGORY: Animation, MPTS CALLING SEQUENCE: IDL> mptscam, shot, xSize=xSize, ySize=ySize, sideSigs=sideSigs, $ MAX_COLORS_TO_USE = max_colors_to_use, dirs=dirs, $ noSides=noSides, colortable=colortable, $ mpegFile=mpegFile, rf=rf, nb=nb, pixMap=pixMap, $ gaps=gaps INPUTS: shot - NSTX shot # to display KEYWORD PARAMETERS: Optional Inputs: ... dirs - directories to search. If not set, typical NSTX camera is set to 'none', no camera file will be used directories will be searched. defPath - path to find most recent "Miro*.cin" file for default selection first - if set, will just return the first match (for speed) pick - if set, and more than one match found, will pop a dialog box for user to pick one. bigger - DEFAULT=1, so will reduce other graphs on left, so camera image is bigger. nSmooth - # of pixels to median smooth over (3 is great for bit noise) verbose - if set, will print many informational messages debug - if set, debug output will be printed OUTPUTS: VCR-like widget; can create mpegs EXAMPLE: for output example, see http://nstx.pppl.gov/nstx/Software/Diagnostics/MPTS/mptsdata.html IDL> mptscam, 130376, dir='/p/nstxcam-archive/Miro2-7988/2008/' IDL> mptscam, 140300, dir='/p/nstxcam/Phantom73-8032/2010/' IDL> f=findcamfiles(141833, /pick) ; select camera from those available IDL> fdecomp, f, disk, dir, name IDL> mptscam, 141833, dir=dir IDL> mptscam, 141833, dir='/p/nstxcam/Phantom4-6878/2010/' NOTES: This routine has not been tested for many combinations of signals Logic for making movies with less than the full time needs improving. MODIFICATION HISTORY: 12-Oct-2010 cleaned up for new camera searching, and added /bigger keyword 18-Aug-2009 limit radial plots from 18-160 cm per Adam McLean. 14-Aug-2009 major changes for new camera data locations 16-Nov-01 Fine tuned for xsize=600, ysize=400 (for publication) Written by Bill Davis
Source: src/idl_cvs/mptsout.pro
src/idl_cvs/mptsout.pro
NAME: mptsout PURPOSE: write a table of MPTS data and errors from the edge along with GPI average CATEGORY: MPTS, MDSplus CALLING SEQUENCE: mptsout, shot, r1=r1, r2=r2, t1=t1, t2=t2, outfile=outfile INPUTS: shot - MDSplus shot number, i.e, ID. KEYWORDS: (all Optional) r1 - starting radius (cm) for printed values (default=120 cm) r2 - starting radius (cm) for printed values (default=156 cm) t1 - Starting x-value desired (typically time, in seconds). Defaults to 0.3 s t2 - Ending x-value desired. Defaults to 0.9 s outFile - name defaults to 'MPTS_'+shot+'.txt' EXAMPLE: IDL> for shot=138844,138847 do mptsout, shot HISTORY: Written 04-Jun-2013 by Bill Davis for Stewart Zweben
Source: src/idl_cvs/mptssurface.pro
src/idl_cvs/mptssurface.pro
NAME: mptssurface PURPOSE: A graphical interface for plotting shaded surfaces of NSTX Multi-point Thomson Scattering data (MPTS) from MDSplus. Summary plots of Ip, D-alpha, Neutral Beam and RF input power, and stored energy (W-mhd from EFIT). Times of NB and RF heating are indicated on the surfaces by darker shading. Different controls are provided to change the viewing angle and other plot parameters. Projection of data maxima vs. radius may be projected on the "back wall". (For critical accuracy of values, especially on the edges, please check with the Thomson diagnosticians personally.) CATEGORY: 3-D Plotting, MDSplus, NSTX, Widgets CALLING SEQUENCE: mptssurface, shot INPUT PARAMETERS: shot - MDSplus shot # KEYWORD PARAMETERS: autoRefresh - if set, will automatically refresh when MDSplus event MPTS_REFRESH is declared refreshEvent - default is 'MPTS_REFRESH' warning - if=0, no widgets will be popped for warnings charsize - character size of plots (DEFAULT=2) density - if set, will display Density instead of Temperature initially contour - if=0, will not overplot contours on the surface project - if=1, will project maxima of variable vs. radius on "back wall" GROUP: The widget ID of the widget that calls mptssurface. When this keyword is specified, the death of the caller results in the death of mptssurface. BLOCK: Set this keyword to have XMANAGER block when this application is registered. By default the Xmanager keyword NO_BLOCK is set to 1 to provide access to the command line if active command line processing is available. Note that setting BLOCK for this application will cause all widget applications to block, not only this application. For more information see the NO_BLOCK keyword to XMANAGER. LIMITATIONS: Some liberties have been taken "massaging" the data from MDSplus to make nice looking surfaces, such as median smoothing with more on the edges, and interpolating between time points and radii of the raw data. mptssurface does not accept all of the keywords that the IDL command SURFACE does. EXAMPLE: IDL> mptssurface, 142000, /project MODIFICATION HISTORY: 06-Jun-2013 made shades work when zooming in 01-Aug-2011 Added Kelvin keyword for presentation to non-physicists 30-Nov-2010 Added color shading for other variables onto surface plot. 11-Nov-2010 Added contours and projection of maxima on back wall 28-Oct-2010 improved data filtering per Ben LeBlanc and Ahmed Diallo 22-Oct-2010 converted from MPTSsurface for Joel Hosea 20-May-04 Better message when no analysis data, or times out-of-range 10-Feb-03 Made shaded surface the default [BD] 25-Feb-02 List shot date and time on plot [BD] 01-Feb-01 Added "Copy Plot to Tek Window" option 14-Sep-00 Added min & max R options & better shot buttons. [BD] 29-Aug-00 Made from RSI xsurface [BD] Created from a template written by: Steve Richards, January, 1991.
Source: src/idl_cvs/msg_bell.pro
src/idl_cvs/msg_bell.pro
NAME: Msg_Bell PURPOSE: Print a message to the IDL session and ring a bell CATEGORY: Programming. Terminal Outpout CALLING SEQUENCE: Msg_Bell, msg, nrings INPUTS: msg - string to print to the IDL session nrings - (OPTIONAL) # of times to ring KEYWORD PARAMETERS: OUTPUTS: EXAMPLE: COMMON BLOCKS: NOTES: MODIFICATION HISTORY: 1997 Written by Bill Davis
Source: src/idl_cvs/multiplot.pro
src/idl_cvs/multiplot.pro
NAME: MULTIPLOT Purpose: Create multiple plots with shared axes. Category: Plotting Explanation: This procedure makes a matrix of plots with *SHARED AXES*, either using parameters passed to multiplot or !p.multi in a non-standard way. It is good for data with one or two shared axes and retains all the versatility of the plot commands (e.g. all keywords and log scaling). The plots are connected with the shared axes, which saves space by omitting redundant ticklabels and titles. Multiplot does this by setting !p.position, !x.tickname and !y.tickname automatically. A call (multiplot,/reset) restores original values. Note: This method may be superseded by future improvements in !p.multi by RSI. For now, it's a good way to gang plots together. CALLING SEQUENCE: multiplot[pmulti][,/help][,/initialize][,/reset][,/rowmajor] Examples: multiplot,/help ; print this header. ; Then copy & paste, from your xterm, the following lines to test: x = findgen(100) ; MULTIPLOT t=exp(-(x-50)^2/300) ; ------------------------- erase ; | | | u=exp(-x/30) ; | | | y = sin(x) ; | UL plot | UR plot | r = reverse(y*u) ; | | | !p.multi=[0,2,2,0,0] ; | | | multiplot ; y------------------------- plot,x,y*u,title='MULTIPLOT' ; l| | | multiplot & plot,x,r ; a| | | multiplot ; b| LL plot | LR plot | plot,x,y*t,ytit='ylabels' ; e| | | multiplot ; l| | | plot,x,y*t,xtit='xlabels' ; s------------------------- multiplot,/reset ; xlabels wait,2 & erase ; TEST multiplot,[1,3] ; H------------------------ plot,x,y*u,title='TEST' ; E| plot #1 | multiplot ; I------------------------ plot,x,y*t,ytit='HEIGHT' ; G| plot #2 | multiplot ; H------------------------ plot,x,r,xtit='PHASE' ; T| plot #3 | multiplot,/reset ; ------------------------ ; PHASE multiplot,[1,1],/init,/verbose ; one way to return to single plot % MULTIPLOT: Initialized for 1x1, plotted across then down (column major). Optional Inputs: pmulti = 2-element or 5-element vector giving number of plots, e.g., multiplot,[1,6] ; 6 plots vertically multiplot,[0,4,2,0,0] ; 4 plots along x and 2 along y multiplot,[0,4,2,0,1] ; ditto, except rowmajor (down 1st) multiplot,[4,2],/rowmajor ; identical to previous line Optional Keywords: help = flag to print header initialize = flag to begin only---no plotting, just setup, e.g., multiplot,[4,2],/init,/verbose & multiplot & plot,x,y reset = flag to reset system variables to values prior to /init default = flag to restore IDL's default value for system variables rowmajor = flag to number plots down column first (D=columnmajor) verbose = flag to output informational messages Outputs: !p.position = 4-element vector to place a plot !x.tickname = either '' or else 30 ' ' to suppress ticknames !y.tickname = either '' or else 30 ' ' to suppress ticknames !p.noerase = 1 Common blocks: multiplot---to hold saved variables and plot counter. See code. Side Effects: Multiplot sets a number of system variables: !p.position, !p.multi, !x.tickname, !y.tickname, !P.noerase---but all can be reset with the call: multiplot,/reset Restrictions: 1. If you use !p.multi as the method of telling how many plots are present, you have to set !p.multi at the beginning each time you use multiplot or call multiplot with the /reset keyword. 2. There's no way to make an xtitle or ytitle span more than one plot, except by adding spaces to shift it or to add it manually with xyouts. 3. There is no way to make plots of different sizes; each plot covers the same area on the screen or paper. Procedure: This routine makes a matrix of plots with common axes, as opposed to the method of !p.multi where axes are separated to allow labels. Here the plots are joined and labels are suppressed, except at the left edge and the bottom. You tell multiplot how many plots to make using either !p.multi (which is then reset) or the parameter pmulti. However, multiplot keeps track of the position by itself because !p.multi interacts poorly with !p.position. Modification history: 14-Dec-00 Account for !x.omargin and !x.ymargin [Bill Davis, PPPL] write, 21-23 Mar 94, Fred Knight (knight@ll.mit.edu) alter plot command that sets !x.window, etc. per suggestion of Mark Hadfield (hadfield@storm.greta.cri.nz), 7 Apr 94, FKK add a /default keyword restore IDL's default values of system vars, 7 Apr 94, FKK modify two more sys vars !x(y).tickformat to suppress user-formatted ticknames, per suggestion of Mark Hadfield (qv), 8 Apr 94, FKK Converted to IDL V5.0 W. Landsman September 1997
Source: src/idl_cvs/mvshotrangedir.pro
src/idl_cvs/mvshotrangedir.pro
NAME: mvshotrangedir PURPOSE: Put files fitting a certain partern into a sub-directory of the form nnnn00s, where nnnn00 represents a shot range, e.g., 114100s whould contain shots 114100-114199. CATEGORY: Files, Dates CALLING SEQUENCE: IDL> mvshotrangedir, path=path, prefix=prefix, ext=ext, date=date INPUTS: none KEYWORD PARAMETERS: path - directory prefix - file prefix for deletion ext - file extension for deletion; file search will be path+prefix+'*.'+ext verbose - if set, will print informational output OUTPUTS: none COMMON BLOCKS: NONE EXAMPLE: IDL> mvshotrangedir,/verb IDL> for i=1408, 1420 do mvshotrangedir, fourDigits=strtrim(i), /verb NOTES: LIMITATIONS Only works on Unix/Linux MODIFICATION HISTORY: 21-Oct-2010 made work 23-May-2006 do all files in dir 15-May-2006 Written by Bill Davis, PPPL