Upgrading to GitLab 11.11.0.

Commit ffebf678 authored by Ole Streicher's avatar Ole Streicher

Initial upstream version 1.82

parents
This diff is collapsed.
MARKWARDT IDL PROGRAMS
INSTALLATION INSTRUCTIONS
Craig Markwardt
craigm@cow.physics.wisc.edu
31 Jan 2000
The following instructions apply to all functions and procedures
included in the Markwardt IDL library.
DOWNLOADING
Download new versions of from Craig Markwardt's web page:
http://cow.physics.wisc.edu/~craigm/idl/idl.html
Program modification dates appear on the web page, which you can
compare agains your own copy. You can also check the modification
history of the file itself to see how recent it is.
COMPATIBILITY
All programs are intended to be compatible with both Unix, Windows and
MacIntosh versions of IDL, unless otherwise noted. If you are
downloading archive, then the most suitable archive file type can be
found from this chart, depending on your operating system:
Windows ZIP archive (must have unzipping utility)
MacIntosh ZIP archive (must have unzipping utility)
Unix TAR or ZIP archive (must have gunzip or unzip utility)
INSTALLATION
Program files should be placed in your IDL PATH. You may either copy
them to an existing directory in your path, or else create a new
directory for them to reside in. Please see your IDL documentation
for instructions on how to add directories to your IDL path.
You may read the individual program's documentation, stored with the
program itself, for more information on usage.
If the files you downloaded came in an archive, you must unarchive
them. It is recommended that you unarchive them into a separate
directory for convenience's sake.
For Windows and MacIntosh systems, you should use an appropriate
unzipping utility to expand the ZIP archive you downloaded.
For Unix systems, you should use one of the following commands,
whichever is appropriate:
gzip -dc <archive>.tar.gz | tar xvf - (for TAR file; or...)
unzip <archive>.zip (for ZIP file)
where <archive> is the name of the archive you downloaded.
IDL is a registered trademark of RSI
MPFIT PACKAGE
MARKWARDT IDL PROGRAMS
Craig Markwardt
Craig.Markwardt@gmail.com
22 Nov 2009
The following instructions apply to the MPFIT and TNMIN packages of
functions for curve fitting under IDL, available from the Markwardt
IDL Library.
MPFIT is a set of routines for robust least-squares minimization
(curve fitting), using arbitrary user written IDL functions or
procedures. MPFIT is based on the well-known and tested MINPACK-1
FORTRAN package of routines available at www.netlib.org. The relevant
sections of code have been translated almost directly from the FORTRAN
equivalent.
MPFIT functions are designed for consistent usage and have some
special features not found in the standard IDL routines. MPFIT
functions permit you to fix any function parameters, as well as to set
simple upper and lower parameter bounds. See the documentation under
PARINFO for instructions on how to use this facility.
When referring to MPFIT in your scientific papers, please cite
the following paper. You can use the NASA ADS link below to get
a citation in the preferred format for your journal.
Markwardt, C. B. 2008, "Non-Linear Least Squares Fitting in IDL
with MPFIT," in proc. Astronomical Data Analysis Software and
Systems XVIII, Quebec, Canada, ASP Conference Series, Vol. XXX, eds.
D. Bohlender, P. Dowler & D. Durand (Astronomical Society of the
Pacific: San Francisco), p. 251-254 (ISBN: 978-1-58381-702-5)
http://arxiv.org/abs/0902.2850
Link to NASA ADS: http://adsabs.harvard.edu/abs/2009ASPC..411..251M
Link to ASP: http://aspbooks.org/a/volumes/table_of_contents/411
Refer to the MPFIT website as:
http://purl.com/net/mpfit
TNMIN
TNMIN solves general minimization problems (generally not curve
fitting), and has similar features to the MPFIT family of functions.
As of this writing, TNMIN requires the user function to supply
derivatives.
DOWNLOADING
Download new versions of from Craig Markwardt's web page:
http://purl.com/net/mpfit
Program modification dates appear on the web page, which you can
compare agains your own copy. You can also check the modification
history of the file itself to see how recent it is.
Please see the file INSTALL for installation instructions.
MANIFEST
The following functions are included:
INSTALL - installation instructions
MPREADME - this file
MPFIT - main fitting engine, required for other driver functions
MPFITFUN - driver function for 1D function fitting
MPFIT2DFUN - driver function for 2D function fitting (images)
MPCURVEFIT - drop-in replacement for IDL's CURVEFIT, requires MPFIT
MPFITEXPR - driver function for fitting expressions interactively
MPFITPEAK - driver function for fitting Gaussian, Lorentzian or Moffat peaks
MPFIT2DPEAK - driver function for fitting 2D peaks
GAUSS1 - example 1D gaussian function
GAUSS1P - example 2D gaussian *procedure*
fakedata.sav - example 1D gaussian data
USAGE
The general theory of curve fitting is beyond the scope of this
document. However, I can offer you a few suggestions. First, read
the fitting tutorial found on my web page here:
http://purl.com/net/mpfit
(and click through to the Fitting section). This should get you
started on the basics of 1D fitting, and in fact 2D fitting too since
the principles are almost the same.
Second, read the documentation! Each of the program files is
extensively self-documented in their comment headers. You can either
read the files directly or download the documentation from my web page
in the documentation section.
Finally, don't be afraid to experiment.
RECOMMENDATIONS
There are a lot of fitting functions available, each one optimized for
a specific task. Allow me to suggest which one to use:
* For 1D curve fitting, use MPFITFUN;
* For 2D surface fitting, use MPFIT2DFUN;
* For existing programs which already use CURVEFIT, use MPCURVEFIT;
* For general non-linear minimization problems, use TNMIN
The main engine, MPFIT, is required in all cases, since the driver
functions call it.
IDL is a registered trademark of RSI
File added
This diff is collapsed.
;+
; NAME:
; GAUSS1
;
; AUTHOR:
; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
; craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
; Compute Gaussian curve given the mean, sigma and area.
;
; MAJOR TOPICS:
; Curve and Surface Fitting
;
; CALLING SEQUENCE:
; YVALS = GAUSS1(XVALS, [MEAN, SIGMA, AREA], SKEW=skew)
;
; DESCRIPTION:
;
; This routine computes the values of a Gaussian function whose
; X-values, mean, sigma, and total area are given. It is meant to be
; a demonstration for curve-fitting.
;
; XVALS can be an array of X-values, in which case the returned
; Y-values are an array as well. The second parameter to GAUSS1
; should be an array containing the MEAN, SIGMA, and total AREA, in
; that order.
;
; INPUTS:
; X - Array of X-values.
;
; [MEAN, SIGMA, AREA] - the mean, sigma and total area of the
; desired Gaussian curve.
;
; INPUT KEYWORD PARAMETERS:
;
; SKEW - You may specify a skew value. Default is no skew.
;
; PEAK - if set then AREA is interpreted as the peak value rather
; than the area under the peak.
;
; RETURNS:
;
; Returns the array of Y-values.
;
; EXAMPLE:
;
; p = [2.2D, 1.4D, 3000.D]
; x = dindgen(200)*0.1 - 10.
; y = gauss1(x, p)
;
; Computes the values of the Gaussian at equispaced intervals
; (spacing is 0.1). The gaussian has a mean of 2.2, standard
; deviation of 1.4, and total area of 3000.
;
; REFERENCES:
;
; MODIFICATION HISTORY:
; Written, Jul 1998, CM
; Correct bug in normalization, CM, 01 Nov 1999
; Optimized for speed, CM, 02 Nov 1999
; Added copyright notice, 25 Mar 2001, CM
; Added PEAK keyword, 30 Sep 2001, CM
;
; $Id: gauss1.pro,v 1.4 2001/10/13 17:41:48 craigm Exp $
;
;-
; Copyright (C) 1998,1999,2001, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
function gauss1, x, p, skew=skew, peak=peak, _EXTRA=extra
sz = size(x)
if sz(sz(0)+1) EQ 5 then smax = 26D else smax = 13.
if n_elements(p) GE 3 then norm = p(2) else norm = x(0)*0 + 1
u = ((x-p(0))/(abs(p(1)) > 1e-20))^2
mask = u LT (smax^2)
if NOT keyword_set(peak) then norm = norm / (sqrt(2.D * !dpi)*p(1))
f = norm * mask * exp(-0.5*temporary(u) * mask)
mask = 0
if n_elements(skew) GT 0 then $
f = (1.D + skew * (x-p(0))/p(1))*f
return, f
end
;+
; NAME:
; GAUSS1P
;
; AUTHOR:
; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
; craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
; Compute Gaussian curve given the mean, sigma and area (procedure).
;
; MAJOR TOPICS:
; Curve and Surface Fitting
;
; CALLING SEQUENCE:
; GAUSS1, XVALS, [MEAN, SIGMA, AREA], YVALS, SKEW=skew
;
; DESCRIPTION:
;
; This routine computes the values of a Gaussian function whose
; X-values, mean, sigma, and total area are given. It is meant to be
; a demonstration for curve-fitting.
;
; XVALS can be an array of X-values, in which case the returned
; Y-values are an array as well. The second parameter to GAUSS1
; should be an array containing the MEAN, SIGMA, and total AREA, in
; that order.
;
; INPUTS:
; X - Array of X-values.
;
; [MEAN, SIGMA, AREA] - the mean, sigma and total area of the
; desired Gaussian curve.
;
; YVALS - returns the array of Y-values.
;
;
; KEYWORD PARAMETERS:
;
; SKEW - You may specify a skew value. Default is no skew.
;
; EXAMPLE:
;
; p = [2.2D, 1.4D, 3000.D]
; x = dindgen(200)*0.1 - 10.
; gauss1p, x, p, y
;
; Computes the values of the Gaussian at equispaced intervals
; (spacing is 0.1). The gaussian has a mean of 2.2, standard
; deviation of 1.4, and total area of 3000.
;
; REFERENCES:
;
; MODIFICATION HISTORY:
; Transcribed from GAUSS1, 13 Dec 1999, CM
; Added copyright notice, 25 Mar 2001, CM
;
; $Id: gauss1p.pro,v 1.2 2001/03/25 18:55:12 craigm Exp $
;
;-
; Copyright (C) 1999,2001, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
pro gauss1p, x, p, f, skew=skew, _EXTRA=extra
sz = size(x)
if sz(sz(0)+1) EQ 5 then smax = 26D else smax = 13.
if n_elements(p) GE 3 then norm = p(2) else norm = x(0)*0 + 1
u = ((x-p(0))/(abs(p(1)) > 1e-20))^2
mask = u LT (smax^2) ;; Prevent floating underflow
f = norm * mask * exp(-0.5*temporary(u) * mask) / (sqrt(2.D * !dpi)*p(1))
mask = 0
if n_elements(skew) GT 0 then $
f = (1.D + skew * (x-p(0))/p(1))*f
return
end
;+
; NAME:
; GAUSS2
;
; AUTHOR:
; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
; craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
; Compute Gaussian curve given the mean, sigma and area.
;
; MAJOR TOPICS:
; Curve and Surface Fitting
;
; CALLING SEQUENCE:
; YVALS = GAUSS2(X, Y, [XCENT, YCENT, SIGMA, PEAK])
;
; DESCRIPTION:
;
; This routine computes the values of a Gaussian function whose
; X-values, mean, sigma, and total area are given. It is meant to be
; a demonstration for curve-fitting.
;
; XVALS can be an array of X-values, in which case the returned
; Y-values are an array as well. The second parameter to GAUSS1
; should be an array containing the MEAN, SIGMA, and total AREA, in
; that order.
;
; INPUTS:
; X - 2-dimensional array of "X"-values.
; Y - 2-dimensional array of "Y"-values.
;
; XCENT - X-position of gaussian centroid.
; YCENT - Y-position of gaussian centroid.
;
; SIGMA - sigma of the curve (X and Y widths are the same).
;
; PEAK - the peak value of the gaussian function.
;
; RETURNS:
;
; Returns the array of Y-values.
;
; 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.)
; z = gauss2(x, y, p)
;
; Computes the values of the Gaussian at equispaced intervals in X
; and Y (spacing is 0.1). The gaussian has a centroid position of
; (2.2, -0.7), standard deviation of 1.4, and peak value of 3000.
;
; REFERENCES:
;
; MODIFICATION HISTORY:
; Written, 02 Oct 1999, CM
; Added copyright notice, 25 Mar 2001, CM
;
; $Id: gauss2.pro,v 1.2 2001/03/25 18:55:13 craigm Exp $
;
;-
; Copyright (C) 1999,2001, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
function gauss2, x, y, p, _EXTRA=extra
u = ((x-p(0))/p(2))^2 + ((y-p(1))/p(2))^2
mask = u LT 100
f = p(3) * mask * exp(-0.5D * temporary(u) * mask)
mask = 0
return, f
end
;+
; NAME:
; LINFITEX
;
; 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:
; Model function for fitting line with errors in X and Y
;
; MAJOR TOPICS:
; Curve and Surface Fitting
;
; CALLING SEQUENCE:
; parms = MPFIT('LINFITEX', start_parms, $
; FUNCTARGS={X: X, Y: Y, SIGMA_X: SIGMA_X, SIGMA_Y: SIGMA_Y}, $
; ...)
;
; DESCRIPTION:
;
; LINFITEX is a model function to be used with MPFIT in order to
; fit a line to data with errors in both "X" and "Y" directions.
; LINFITEX follows the methodology of Numerical Recipes, in the
; section entitled, "Straight-Line Data with Errors in Both
; Coordinates."
;
; The user is not meant to call LINFITEX directly. Rather, the
; should pass LINFITEX as a user function to MPFIT, and MPFIT will in
; turn call LINFITEX.
;
; Each data point will have an X and Y position, as well as an error
; in X and Y, denoted SIGMA_X and SIGMA_Y. The user should pass
; these values using the FUNCTARGS convention, as shown above. I.e.
; the FUNCTARGS keyword should be set to a single structure
; containing the fields "X", "Y", "SIGMA_X" and "SIGMA_Y". Each
; field should have a vector of the same length.
;
; Upon return from MPFIT, the best fit parameters will be,
; P[0] - Y-intercept of line on the X=0 axis.
; P[1] - slope of the line
;
; NOTE that LINFITEX requires that AUTODERIVATIVE=1, i.e. MPFIT
; should compute the derivatives associated with each parameter
; numerically.
;
; INPUTS:
; P - parameters of the linear model, as described above.
;
; KEYWORD INPUTS:
; (as described above, these quantities should be placed in
; a FUNCTARGS structure)
; X - vector, X position of each data point
; Y - vector, Y position of each data point
; SIGMA_X - vector, X uncertainty of each data point
; SIGMA_Y - vector, Y uncertainty of each data point
;
; RETURNS:
; Returns a vector of residuals, of the same size as X.
;
; EXAMPLE:
;
; ; X and Y values
; XS = [2.9359964E-01,1.0125043E+00,2.5900450E+00,2.6647639E+00,3.7756164E+00,4.0297413E+00,4.9227958E+00,6.4959011E+00]
; YS = [6.0932738E-01,1.3339731E+00,1.3525699E+00,1.4060204E+00,2.8321848E+00,2.7798350E+00,2.0494456E+00,3.3113062E+00]
;
; ; X and Y errors
; XE = [1.8218818E-01,3.3440986E-01,3.7536234E-01,4.5585755E-01,7.3387712E-01,8.0054945E-01,6.2370265E-01,6.7048335E-01]
; YE = [8.9751285E-01,6.4095122E-01,1.1858428E+00,1.4673588E+00,1.0045623E+00,7.8527629E-01,1.2574003E+00,1.0080348E+00]
;
; ; Best fit line
; p = mpfit('LINFITEX', [1d, 1d], $
; FUNCTARGS={X: XS, Y: YS, SIGMA_X: XE, SIGMA_Y: YE}, $
; perror=dp, bestnorm=chi2)
; yfit = p[0] + p[1]*XS
;
;
; REFERENCES:
;
; Press, W. H. 1992, *Numerical Recipes in C*, 2nd Ed., Cambridge
; University Press
;
; MODIFICATION HISTORY:
; Written, Feb 2009
; Documented, 14 Apr 2009, CM
; $Id: linfitex.pro,v 1.3 2009/04/15 04:17:52 craigm Exp $
;
;-
; Copyright (C) 2009, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
;
function linfitex, p, $
x=x, y=y, sigma_x=sigma_x, sigma_y=sigma_y, $
_EXTRA=extra
a = p[0] ;; Intercept
b = p[1] ;; Slope
f = a + b*x
resid = (y - f)/sqrt(sigma_y^2 + (b*sigma_x)^2)
return, resid
end
This diff is collapsed.
;+
; NAME:
; MPCHITEST
;
; 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:
; Compute the probability of a given chi-squared value
;
; MAJOR TOPICS:
; Curve and Surface Fitting, Statistics
;
; CALLING SEQUENCE:
; PROB = MPCHITEST(CHI, DOF, [/SIGMA, /CLEVEL, /SLEVEL ])
;
; DESCRIPTION:
;
; The function MPCHITEST() computes the probability for a value drawn
; from the chi-square distribution to equal or exceed the given value
; CHI. This can be used for confidence testing of a measured value
; obeying the chi-squared distribution.
;
; P_CHI(X > CHI; DOF) = PROB
;
; In specifying the returned probability level the user has three
; choices:
;
; * return the confidence level when the /CLEVEL keyword is passed;
; OR
;
; * return the significance level (i.e., 1 - confidence level) when
; the /SLEVEL keyword is passed (default); OR
;
; * return the "sigma" of the probability (i.e., compute the
; probability based on the normal distribution) when the /SIGMA
; keyword is passed.
;
; Note that /SLEVEL, /CLEVEL and /SIGMA are mutually exclusive.
;
; INPUTS:
;
; CHI - chi-squared value to be tested.
;
; DOF - scalar or vector number, giving the number of degrees of
; freedom in the chi-square distribution.
;
; RETURNS:
;
; Returns a scalar or vector of probabilities, as described above,
; and according to the /SLEVEL, /CLEVEL and /SIGMA keywords.
;
; KEYWORD PARAMETERS:
;
; SLEVEL - if set, then PROB describes the significance level
; (default).
;
; CLEVEL - if set, then PROB describes the confidence level.
;
; SIGMA - if set, then PROB is the number of "sigma" away from the
; mean in the normal distribution.
;
; EXAMPLES:
;
; print, mpchitest(1300d,1252d)
;
; Print the probability for a chi-squared value with 1252 degrees of
; freedom to exceed a value of 1300, as a confidence level.
;
; REFERENCES:
;
; Algorithms taken from CEPHES special function library, by Stephen
; Moshier. (http://www.netlib.org/cephes/)
;
; MODIFICATION HISTORY:
; Completed, 1999, CM
; Documented, 16 Nov 2001, CM
; Reduced obtrusiveness of common block and math error handling, 18
; Nov 2001, CM
; Convert to IDL 5 array syntax (!), 16 Jul 2006, CM
; Move STRICTARR compile option inside each function/procedure, 9 Oct 2006
; Add usage message, 24 Nov 2006, CM
; Really add usage message, with /CONTINUE, 23 Sep 2009, CM
;
; $Id: mpchitest.pro,v 1.10 2009/10/05 16:22:44 craigm Exp $
;-
; Copyright (C) 1997-2001, 2006, 2009, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
forward_function cephes_igamc, cephes_igam
;; Set machine constants, once for this session. Double precision
;; only.
pro cephes_setmachar
COMPILE_OPT strictarr
common cephes_machar, cephes_machar_vals
if n_elements(cephes_machar_vals) GT 0 then return
if (!version.release) LT 5 then dummy = check_math(1, 1)
mch = machar(/double)
machep = mch.eps
maxnum = mch.xmax
minnum = mch.xmin
maxlog = alog(mch.xmax)
minlog = alog(mch.xmin)
maxgam = 171.624376956302725D
cephes_machar_vals = {machep: machep, maxnum: maxnum, minnum: minnum, $
maxlog: maxlog, minlog: minlog, maxgam: maxgam}
if (!version.release) LT 5 then dummy = check_math(0, 0)
return
end
function cephes_igam, a, x
;
; Incomplete gamma integral
;
;
;
; SYNOPSIS:
;
; double a, x, y, igam();
;
; y = igam( a, x );
;
; DESCRIPTION:
;
; The function is defined by
;
; x
; -
; 1 | | -t a-1
; igam(a,x) = ----- | e t dt.
; - | |
; | (a) -
; 0
;
;
; In this implementation both arguments must be positive.
; The integral is evaluated by either a power series or
; continued fraction expansion, depending on the relative
; values of a and x.
;
; ACCURACY:
;
; Relative error:
; arithmetic domain # trials peak rms
; IEEE 0,30 200000 3.6e-14 2.9e-15
; IEEE 0,100 300000 9.9e-14 1.5e-14
COMPILE_OPT strictarr
common cephes_machar, machvals
MAXLOG = machvals.maxlog
MACHEP = machvals.machep
if x LE 0 OR a LE 0 then return, 0.D
if x GT 1. AND x GT a then return, 1.D - cephes_igamc(a, x)
ax = a * alog(x) - x - lngamma(a)
if ax LT -MAXLOG then begin
; message, 'WARNING: underflow', /info