Source: 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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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