An Interactive Least-Squares Fitting Package

  1. Introduction
  2. Program Commands
  3. Writing Multi-FRILLS Applications
  4. Genie Applications

[Next Page | Top]

Introduction

Multi-FRILLS is a subroutine library designed to be used in programs that perform interactive non-linear least-squares fitting of a parametrized function to one or more sets of data points. It provides a flexible and easy-to-use interface to a least-squares fitting routine with the minimum of programming on the user's part. The simplest use of Multi-FRILLS is to fit a function to a single set of data points, e.g. counts as a function of time in a single detector. More sophisticated programs can fit several data sets simultaneously, with some parameters global to all data sets and others specific to a particular data set.

The user is required to provide:

When subroutine MULTI_FRILLS has been called, the program operates in command mode. Following the prompt '#', the user types commands to

On completion, the fitted parameters are listed along with their estimated errors and a value for chi^2. The program returns to command mode so that the user can plot the fit and choose new initial conditions, e.g. by clearing fixed parameters, before trying again. The results of fits are recorded in an ASCII file and lists of parameters may be stored in files for subsequent printing or input as initial parameters to another Multi-FRILLS session. The plots may also be stored on disk in Postscript or LN03 files for output to a printer, or the calculated function may be written to ASCII files in tab-delimited form suitable for input to PC or Macintosh programs such as Excel, SigmaPlot, proFit or Kaleidagraph.

To optimize the function parameters, Multi-FRILLS uses either its own non-linear least-squares fitting routine which employs a standard quasi-Newton (inverse Hessian) algorithm, or a standard Levenberg-Marquadt algorithm. For a discussion of the general principles, see section 14.4 of "Numerical Recipes: The Art of Scientific Computing" (Cambridge University Press 1986). Using the calculated values passed back from the user-supplied subroutine, a chi^2 function is defined in the following way:

where v and s are the data values and errors, c are the calculated values, and nv and np are the number of data values and fitting parameters respectively. The best fit parameters are those which minimize chi^2. With this definition, the minimum of chi^2 should be close to

Much larger values of chi^2 mean either that the model function is not an accurate representation of the data or that the data errors are underestimated, whereas smaller values mean that the data errors are too large. The fitting method will find the local chi^2 minimum; there may be a global minimum elsewhere in parameter space. However, since the user can plot the calculated function and compare it to the data for each data set before attempting a fit, it should be possible to define reasonable starting parameters for most models. The criteria for chi^2 convergence may be adjusted along with other fitting constants such as the maximum number of iterations. See the description of the ALTER FIT command in the next section for more details. Estimates of the parameter errors are provided by the diagonal elements of the covariance matrix, which assumes that they are normally distributed i.e. uncorrelated. The covariance matrix can be inspected for evidence of such correlations.

Multi-FRILLS v1.3 is based on FRILLS which is described in the report

"FRILLS: An Interactive Least Squares Fitting Package"
R. Osborn, Rutherford Appleton Laboratory Report RAL 91-011 (1991)

Although the command line interface has changed little from FRILLS, there are some important differences in the functionality of Multi-FRILLS which are summarized here.

At present, Multi-FRILLS can only be installed on VMS machines on which Dec GKS has been installed. We are in the advanced stages of a project to convert the graphics to PGPLOT, which is a freely available library of Fortran graphics routines. It should then be possible to install Multi-FRILLS on platforms other than VMS.

Section 2 describes all the Multi-FRILLS commands in detail; if you are just using a Multi-FRILLS application and not writing your own, this is the only section you need to read. For programmers, a description of how to write Multi-FRILLS applications is given in Section 3 together with a simple example. A Multi-FRILLS application may launched from the data analysis package Open GENIE (or the older GENIEv2), with the data being transferred from one or more GENIE workspaces. A description of how to do this is given in Section 4.


[Previous Page | Next Page | Top]

Program Commands

An important feature of this fitting package is the flexibility that the user has in the way a fit is performed. Commands may be given in any order, with the minimum of constraints. For instance, parameters may be fixed before one fit and cleared before the next. At any stage, the function, defined by the current set of parameters, may be plotted on the terminal screen for any of the several data sets. This section describes all the currently available Multi-FRILLS commands for setting up the parameters, plotting the data and functions, and performing the least-squares fits.

The following conventions are followed by Multi-FRILLS:

In this manual, optional parameters are surrounded by square brackets. Values to be chosen by the user are surrounded by angle brackets.

Multi-FRILLS graphics work both on X-Windows screens and on terminals/emulators with both text and Tektronix 4014 graphics screens, e.g. VersaTerm Pro. On the latter, the terminal is automatically switched to the graphics screen after a plot command. However, the user must type [RETURN] after the plot is completed in order to return to text mode.

At any stage in a Multi-FRILLS session, help can be obtained using the HELP command:

Help
The format of a specified command is listed at the terminal. If no command is specified, all the available commands are listed.
Format: H [<command>]

[Previous Page | Next Page | Top]

Manipulating Parameter Values and Data

Use
Selects the data set to be used for future commands. All commands that manipulate or display parameter values apply only to those that are global or specific to the data set identified by the most recent USE command. Similarly, commands that manipulate, display or plot data will do so only for that data set. For example, the DISPLAY PARAMETERS command lists only global parameters and those associated with the data set chosen with the USE command. Similarly, the REMOVE command, which removes points in a specified x-axis range, applies only to that data set. The GO command, which invokes the minimization algorithm acts on all data sets at once although an individual data set may be eliminated from the fit by removing all its data points. The default on entry is to use the first data set.
Format: U <n>
Set
The user is prompted for values of the specified parameters, and, optionally, their upper and lower limits. If no parameter number is specified, all the global parameters and those parameters associated with the selected data set are listed and the user is prompted to give the new values of each one in turn. If the user types [RETURN], the parameter is left unchanged.
Format: S [<n1>,[<n2-n3>,[n4...]]]
Fix
The specified parameters are fixed in subsequent fits. They may be freed by the CLEAR PARAMETER command.
Format: F <n1>,[<n2-n3>,[n4...]]
Bind
The value of each specified parameter is bound to one other parameter. The user is prompted for the number of the other parameter and the ratio of the two parameters i.e. if p(i) is the specified parameter, and p(j) is the parameter to which it is bound, the required ratio is p(i) / p(j) . If no ratio is given, it is set to the ratio of the current values of the parameters. Parameter binding is cleared by the CLEAR PARAMETER command.
Format: B <n1>,[<n2-n3>,[n4...]]
Limit
The specified parameters are constrained in subsequent fits by upper and lower values, for which the user is prompted. If the limits are to be removed, set both values to 0. In the least-squares fit, constrained parameters are transformed into a sine function so that parameter derivatives do not diverge even at the limiting values. Nevertheless, limiting parameters should only be used when absolutely necessary.
Format: L <n1>,[<n2-n3>,[n4...]]
Input
The parameter values and constraints are read in from a disk file of the same format as produced by the OUTPUT command below.
Format: I <disk file name>
Output Parameters
The current parameters are output to a disk file. If none is specified, the default file name is SYS$SCRATCH:PAR.OUT. The parameters may be read subsequently from the disk file with the INPUT command.
Format: O P [<disk file name>]
Output Fit
The latest fitted parameters are output to a disk file along with their chi^2 value. If no disk file is specified, the default file name is SYS$SCRATCH:FIT.OUT. The parameters may be read subsequently from the disk file with the INPUT command.
Format: O F [<disk file name>]
Remove
Data points between the specified x-values are removed from the fit and will not appear in subsequent plots. They may be restored by the CLEAR DATA command.
Format: R <xmin> <xmax>
Modify
The data points between the specified x-values are modified for subsequent fits. The user is prompted for new y-values. They may be restored to their initial values by the CLEAR DATA command. This facility is used at the user's discretion.
Format: M <xmin> <xmax>
Clear Parameters
The specified parameters are freed to vary in subsequent fits. If no parameters are specified, all are freed.
Format: C P [<n1>,[<n2-n3>,[n4...]]]
Clear Data
The data between the specified limits are restored to their original values (if modified) and included in subsequent fits (if removed). If no limits are specified, all points are restored.
Format: C D [<xmin> <xmax>]

[Previous Page | Next Page | Top]

Displaying And Plotting Parameters And Data

Display Data
The data points of the current data set within the specified range are listed at the terminal. If no range is specified, all points are listed.
Format: D D [<xmin> <xmax>]
Display Parameters
The specified parameters are listed at the terminal. If none are specified, all the global parameters and those parameters associated with the current data set are listed. The first column is the unique parameter number and the second column is the data set with which it is associated. If this is 0, the parameter is global. The remaining columns list the parameter names, values and limits.
Format: D P [<n1>,[<n2-n3>,[n4...]]]
Display Calculation
The data points and the calculated values of the current data set within the specified range using the current parameters are listed at the terminal. If no range is specified, all points are listed.
Format: D C [<xmin> <xmax>]
Display Fit
The fitted parameters are listed at the terminal in a similar fashion to the DISPLAY PARAMETER command with, in addition, the chi^2 value and status of the fit.
Format: D F
Display Matrix
The correlation matrix of the fitted parameters is listed at the terminal.
Format: D M
Display License
The terms of the GNU General Public License are output to the terminal.
Format: D L
Plot Data
The data points of the current data set are plotted at the terminal with the option of setting the axis limits. If DATA is preceded by OVER, i.e. PLOT OVER DATA, the data are overlaid on the previous plot.
Format: P [O] D [<xmin> <xmax> [<ymin> <ymax>]]
Plot Calculation
The data points of the current data set and the calculated values using the current parameters are plotted at the terminal with the option of setting the axis limits. If CALCULATION is preceded by OVER, i.e. PLOT OVER CALCULATION, the calculated lines are overlaid on the previous plot and the data points are not replotted.
Format: P [O] C [<xmin> <xmax> [<ymin> <ymax>]]
Plot Fit
The data points of the current data set and their fitted values are plotted at the terminal with the option of setting the axis limits. If FIT is preceded by OVER, i.e. PLOT OVER FIT, the fitted lines are overlaid on the previous plot and the data points are not plotted.
Format: P [O] F [<xmin> <xmax> [<ymin> <ymax>]]
Plot Residuals
The difference between the data points of the current data set and their fitted values are plotted at the terminal with the option of setting the axis limits. If RESIDUALS is preceded by OVER, i.e. PLOT OVER RESIDUALS, the plotted points are overlaid on the previous plot.
Format: P [O] R [<xmin> <xmax> [<ymin> <ymax>]]
Plot <no. list>
The individual components of a fit are plotted separately. The numbers refer to the order in which the components are stored in the user-supplied function subroutine (see section 3, and common block /FR_CMPS/). If the list of component numbers is preceded by OVER, i.e. PLOT OVER , the lines are overlaid on the previous plot.
Format: P [O] [<n1>,[<n2-n3>,[n4...]]]
Output Calculation
The data points of the current data set and the calculated values of the fitting function and its components using the current parameters are written to an ASCII disk file in tab-delimited form. The first line contains the title, the second line contains column headings, and subsequent lines contain the data, function calculations and function components. Note that in general, there are more calculated points than data points, depending on the smoothing option selected, so there are two columns of x-values. Tab-delimited files are easily read by common PC and Macintosh programs such as Excel, Kaleidagraph and SigmaPlot. If no disk file is specified, the default file name is SYS$SCRATCH:CAL.OUT.
Format: O C [<disk file name>]
Keep
The plots on the screen are stored in a Postscript A4 (default), Postscript US letter or LN03 file, depending on the hard copy type selected with the ALTER HARDCOPY command. If no file name is specified, the default is SYS$SCRATCH:FRILLS.PS for Postscript files and SYS$SCRATCH:FRILLS.TEK for LN03 files. If a fit or calculation has been plotted, the parameters are included at the side of the plot even though they do not appear on the terminal or workstation screen.
Format: K [<disk file>]

[Previous Page | Next Page | Top]

Other Commands

Alter Fits
Adjust least-squares fitting constants. The user is prompted for the new value of each one in turn. Type [RETURN] to leave them unchanged. The parameters that can be changed are:
Format: A F

The last parameter illustrates where the quasi-Newton fitting algorithm used as the default in Multi-FRILLS parts company with other common non-linear least-squares methods. When the inverse Hessian method produces an increase in chi^2, Multi-FRILLS halves the most recently calculated parameter steps until a smaller value of chi^2 results. The minimization then continues using the new values. This simple procedure has proved to be very successful in the sorts of problems most frequently solved by Multi-FRILLS. However, the more commonly-used Levenberg-Marquadt algorithm is available as an option.

Alter Output
Controls the printing of output to the screen and to disk file during a fit. A list of options is provided if the output code is omitted. Type [RETURN] to leave the output code unchanged. The default is to print chi^2 to the screen for each iteration with no output to the disk file.
Format: A O [<ncode>]
The possible values of ncode are:
0	No output to terminal during fit.
1 Output chi^2 to terminal during fit.
2 Output chi^2 and params. to terminal during fit.
10 As for option 0, but, in addition, print fitted points in file.
11 As for option 1, but, in addition, print fitted points in file.
12 As for option 2, but, in addition, print fitted points in file.
Alter Markers
Changes the marker type, size and colour used when plotting data. If no parameter values are entered on the line the user is prompted for the new value of each one in turn. The meaning of the parameters are the same as the DEC GKS marker attributes. Type [RETURN] to leave the values unchanged.
Format: A M [<ntype>, [<size> ], <ncolour>]]]
Alter Lines
Changes the line-type, thickness and colour used when plotting calculations or fits. If no parameter values are entered on the command line the user is prompted for the new value of each one in turn. . The meaning of the parameters are the same as the DEC GKS line-type attributes. Type [RETURN] to leave the values unchanged.
Format: A L [<ntype>, [<thickness> ], ]]]
Alter Smoothing
Adjusts the smoothness of calculated and fitted functions when plotted to the screen. The user is prompted for the new value of two parameters:

Type [RETURN] to leave them unchanged. The function is calculated at N_PLT (default 1000) x-coordinates equally spaced over the range of the plot. N_PLT should be reduced if the function is computationally intensive to evaluate. If the user sets N_PLT=1, the function is calculated only at those x-coordinates at which there are data and which have not been REMOVEd. If the user sets SMOOTH=1, then from the N_PLT calculated data points a smooth curve is generated at 1000 equally spaced x-coordinates which passes through the calculated data points. This produces a smooth curve even if the function is evaluated at only a few points but may be an inaccurate representation of the true function if the number of nodes is small. By default SMOOTH=0 and the calculated values are joined by straight lines.

Format: A S

N.B. the smoothing options only affect the appearance of the calculated function in plots. chi^2 is only calculated using the data points themselves, so the fitting results are not affected by these options.

Alter Workstation
Alter the terminal emulation and workstation type. The user is prompted for new values. Type [RETURN] to leave the values unchanged. The default device emulation assumes either an X windows terminal or a terminal with Tektronix 4014 emulation. Alternatively, Tektronix 4014 emulation suitable for a PC running Kermit can be chosen. The default workstation type produces graphical output using X windows/Motif or Tektronix 4014 drivers; the program translates DECW$DISPLAY to determine which to use i.e. it tests to see if an X server has been specified either by the default terminal configuration or by giving a SET DISPLAY/CREATE command explicitly. The graphics driver can be changed by entering the corresponding GKS workstation type.
Format: A W
Alter Hardcopy
Alter the type of hardcopy output produced by the KEEP command. If no parameter value is given, the user is prompted for a new value. Type [RETURN] to leave the hardcopy type unchanged. Three types are available: Postscript A4 (default type), Postscript US letter, and LN03 file.
Format: A H <ntype>
Title
The graph title is changed. The user is prompted for the new title if is not entered on the command line. Type [RETURN] to leave the title unchanged.
Format: T <new-title>
Jump
A DCL sub-process is spawned. If a DCL command is given, it is executed and control returns to the Multi-FRILLS session. Otherwise, a DCL process is started from which the user must log out in order to return to Multi-FRILLS.
Format: J [<DCL command>]
@<user command>
The user may define commands to be invoked by preceding them with the '@' sign. The command line is passed to a subroutine called USER_COMMAND. A dummy subroutine of this name is supplied with the Multi-FRILLS library, but will be superceded if the user links another routine of the same name. This routine will need to parse the input line to obtain the required command and/or arguments. Look at FRILLS_SOURCES:USER_COMMAND.FOR for an example of using the integer*4 function GETWRD, supplied with Multi-FRILLS, to extract the first word after '@' from the command line.
Go
The free parameters are fitted to the data points. On completion, the fitted parameters are listed at the terminal along with a chi^2 value and the user is asked if the new parameters should be transferred to the initial parameter list. The fitted points may optionally be plotted using the command P F, or listed at the terminal using the command D F. The GO command minimizes chi^2 calculated over all the data sets.
Format: G
Exit
Control is passed back to the program calling MULTI_FRILLS.
Format: E

[Previous Page | Next Page | Top]

Writing Multi-FRILLS Applications

In this section all the essential initialization required by Multi-FRILLS is described, together with the input and output variables of the fitting function subroutine. An example application and an example session are also provided. The source code for the example program and the GENIE examples in Section 4, is provided with the Multi-FRILLS source code and installation package. Following installation of the package, two logical names should be defined:

If these have not been defined at the system level, please consult your system administrator, or whoever installed the Multi-FRILLS library to find out where the relevant subdirectories are located. You should then add the definitions to your own LOGIN.COM e.g.


$ DEFINE FRILLS_SOURCES DUA0:[SOME.SUBDIRECTORY.SOURCES]
$ DEFINE FRILLS_LIBRARY DUA0:[SOME.SUBDIRECTORY.LIBRARY]


[Previous Page | Next Page | Top]

Calling Programs

The user needs to provide the main program which should initialize the data points, axis labels and the parameter names before calling MULTI_FRILLS. It is also possible to preset initial parameter values and whether they are to be fitted or not. This is useful for those functions in which a number of parameters are never intended to be fitted, e.g. the sample temperature. Note however that these parameters can be cleared at a later stage by the CLEAR PARAMETER command, so the user should be careful in using this feature. The user-supplied subroutine that defines the fitting function, must be explicitly declared as EXTERNAL in the program that calls MULTI_FRILLS e.g.


      external GAUSS
The calling program should also contain the line

      include 'FRILLS.INC'
which gives access to all the common blocks used by Multi-FRILLS. The user needs to initialize the following variables:

NPTOT I*4 the total number of parameters
NAM(ip) C*20 the names of the parameters (ip=1:NPTOT)
INS(ip) I*4 the data sets with which each parameter is associated
NS I*4 the number of data sets

For each of the data sets (is=1, NS), the following values must be set:

ND(is) I*4 the number of points in each data set (is=1:NS)
X(i,is) R*8 the x-coordinates (i=1:ND(is), is=1:NS)
Y(i,is) R*8 the y-coordinates
YSIG(i,is) R*8 the error bars
XLAB(is), YLAB(is) C*5 names used for the x and y axes when DISPLAYing data
XCAP(is), YCAP(is) C*40 labels used for the x and y axes when PLOTting data
TITLE(is) C*80 title
SUBTITLE(3,is) C*60 optional subtitles to be placed under the plots

Optionally, the user can define each parameter as either global to all data sets or specific to one data set, by filling the array INS(ip) = 0 (global) or INS(ip) = is (specific to data set is). The default is for parameters to be global to all data sets. The array INS is used for display purposes only and is not used in the minimization algorithm.

The Multi-FRILLS command interface is accessed by a call to subroutine MULTI_FRILLS with the fitting function subroutine name as argument e.g.


      call multi_frills (GAUSS)
The include file FRILLS_SOURCES:FRILLS.INC gives access to many more variables which allow the ambitious programmer more control over the initial state of Multi-FRILLS. Details are given in the next section. Whilst the use of include files simplifies the writing of Multi-FRILLS applications, it opens up the danger that the user will inadvertently use variable names that are also contained in the Multi-FRILLS common blocks. The simplest way to reduce this possibility is for the user to avoid the IMPLICIT statement and explicitly declare all variables in the calling program. The DEC FORTRAN compiler will print a warning message if any variable has been declared twice and the compiler option FORTRAN/WARNINGS=DECLARATIONS will print a warning if there are any undeclared variables. This option is used in the command files provided with the library. The common blocks are included in the appendix to this manual for reference.

Diagnostic output from the fit is written to unit 41. If this is not explicitly opened by the calling program, it will be sent to SYS$SCRATCH:FRILLS.TXT.


[Previous Page | Next Page | Top]

Fitting Subroutines

The subroutine defining the fitting function, in following example GAUSS, is expected to calculate the fitting function for only one data set at each call. It has no arguments passed to it directly, but it must contain the statement


      include 'FRILLS_SOURCES:FUNCTION.INC'
      
which gives access to the following input variables:

P(ip) R*8 the current values of the parameters.
IS I*4 the index of the data set for which calculated values are required.
NV I*4 the number of points in the data set.
X(iv) R*8 the x-coordinates of the points (iv=1:NV).
X(iv) R*8 the x-coordinates of the points (iv=1:NV).

and the variables to contain the output from the fitting function:

YCAL(iv) R*8 the calculated function at each of the x values (iv=1:NV)
NC I*4 the number of separate function components that make up YCAL
YC(iv, ic) R*8 the calculated function components (iv=1:NV; ic=1:NC)

Filling NC and YC is optional. They permit the individual components of a multi-component function to be plotted separately with the PLOT command. Although the data arrays in FUNCTION.INC only contain one data set, the parameter array contains all the parameters. The subroutine defining the function must distinguish which parameters are required by the data set being processed. We provide an example of how this may be done.


[Previous Page | Next Page | Top]

Linking Applications

Multi-FRILLS applications can be compiled and linked either on VAX/VMS or Alpha/VMS computers on which DEC GKS has been installed. To compile and link a program to Multi-FRILLS, use the following command:


$ @FRILLS_LIBRARY:FRILLS_LINK <program name(s)> [<library name(s)>]
This compiles and links the program both to the Multi-FRILLS object library FRILLS_LIBRARY:MULTI_FRILLS.OLB, and to the DEC Fortran Binding GKS library SYS$LIBRARY:GKSFORBND.OLB. Other object files or libraries can be attached to <program name> by adding the second parameter to the command line. These must explicitly include the /LIB or /OPT qualifiers as they would appear in a normal link statement. Alternatively, the user can edit FRILLS_LIBRARY:FRILLS_LINK.COM to include these libraries automatically. To allow for reasonable numbers of parameters, data sets and the lengths of data sets Multi-FRILLS requires a large amount of computer memory. The program may fail to compile as it might exceed the user's page file quota. Please see your system administrator if this happens.


[Previous Page | Next Page | Top]

Array Sizes

The Multi-FRILLS library currently has the following defaults:

The number of parameters that are allowed to vary in a fit should be much less than 150 for acceptable running time.

The limits are contained in the file FRILLS_SOURCES:FRILLS.PAR which may be edited. However, the whole Multi-FRILLS library would need to be recompiled before linking the new program. If Multi-FRILLS is installed centrally at your site, consult the system administrator about changing the options. Alternatively, you could download the Multi-FRILLS library in order to have your own customized copy.


[Previous Page | Next Page | Top]

Example Program

The following program is an example of a Multi-FRILLS application to fit a single Gaussian function to a number of data sets. Although the parameters can be different for each data set, Multi-FRILLS allows the user to bind parameters together. For example, if the user believes that the centres of the peak should be the same in each data set, they can all be bound to a common value. In the calling program, two data sets are read in from FRILLS_SOURCES:GAUSS_FIT.DAT, though the rest of the program would accommodate more data sets. The source code is contained in FRILLS_SOURCES:GAUSS_FIT.F. The data sets provided with the library have been produced by defining a Gaussian on top of a sloping backgroun. Statistical fluctuations have been simulated by assuming a Gaussian distribution of errors with variances equal to the data values. For comparison with the fitted parameters, the Gaussian and background parameters were

Parameter Dataset 1 Dataset 2
Height 200 100
FWHM 1.5 3
Centre 5 5
Background Constant 10 10
Background Slope 1 1


FRILLS_SOURCES:GAUSS_FIT.F
!------------------------------------------------------------------------------------
! Multi-FRILLS  - An Interactive Least Squares Fitting Package for Multiple Datasets
! 
! Copyright (C) 1996 R. Osborn, T. G. Perring (See MULTI_FRILLS.FOR for details)
!
!     Version   Date    Name       Comments
!     =======   ====    ====       ========
!       1.0   96/08/23  RO/TGP     First release of multiple dataset version of FRILLS
!
!------------------------------------------------------------------------------------
!     program GAUSS_FIT
!
!     Program to fit single Gaussian lineshape to more than one dataset
!
!     Supplied as an example of a simple Multi-FRILLS application
!
      program GAUSS_FIT

      include 'FRILLS.INC'
      external GAUSS
      integer*4 i, j, k

!Read in the data
      open (unit = 10, file = 'FRILLS_SOURCES:GAUSS_FIT.DAT', status = 'OLD')
      ns = 2                                !No. of datasets
      nd(1) = 41                            !No. of points in dataset 1
      nd(2) = 41                            !No. of points in dataset 2
      do i = 1,41
            read (10, *) x(i,1), y(i,1), ysig(i,1), y(i,2), ysig(i,2)
            x(i,2) = x(i,1)
      end do
      close (10)
!Define axis labels and titles
      xcap(1) = 'Energy Transfer [meV]'     !Extended X-caption for graphical output
      xcap(2) = xcap(1)
      xlab(1) = 'En'                        !Abbreviated X-label for data columns
      xlab(2) = xlab(1)
      ycap(1) = 'Neutron Counts'            !Extended Y-caption for graphical output
      ycap(2) = ycap(1)
      ylab(1) = 'Cts'                       !Abbreviated Y-label for data columns
      ylab(2) = ylab(1)
      title(1) = 'Sample Dataset 1'         !Title for dataset 1
      title(2) = 'Sample Dataset 2'         !Title for dataset 2
!Loop over no. of datasets defining parameter names and keep count of parameter nos.
      i = 0
      do j = 1,ns
            write (nam(i+1), '(a,i2,a)') 'Set ', j, ' Height'
            write (nam(i+2), '(a,i2,a)') 'Set ', j, ' FWHM'
            write (nam(i+3), '(a,i2,a)') 'Set ', j, ' Centre'
            write (nam(i+4), '(a,i2,a)') 'Set ', j, ' Bkgd Constant'
            write (nam(i+5), '(a,i2,a)') 'Set ', j, ' Bkgd Slope'
            do k = 1,5
                  ins(i+k) = j             !Parameters are associated with dataset k
                  inp(i+k) = i + k         !All parameters are normally varied
            end do
            i = i + 5
      end do
      nptot = i                            !Total number of parameters
      np = i                               !Number of parameters to be varied
!Open output file to contain diagnostic details of fit 
      open (unit = 41, file = 'SYS$SCRATCH:GAUSS_FIT.OUT', status = 'NEW')
!Pass control to Multi-FRILLS program      
      call multi_frills (GAUSS)
!Close output file
      close (unit=41)         
 
      stop
      end

!------------------------------------------------------------------------------------
!     subroutine GAUSS
!
!     Subroutine defining the Gaussian function to be fitted to the data
!     Only one dataset is calculated in each subroutine call identified by 'is'
!
!     Function parameters and dependent variables are passed through FUNCTION.INC
 
      subroutine GAUSS

      include 'FUNCTION.INC'

      integer*4 i, j
      real*8 h, w, c, bc, bs, s, const
      parameter (const = 2D0*sqrt(2D0*log(2D0)))

!Definition of parameters
      i = 5 * (is - 1)
      h  = p(i+1)
      w  = p(i+2)
      c  = p(i+3)
      bc = p(i+4)
      bs = p(i+5)
!Two components to the fit are defined - the peak and the background respectively
      nc = 2
!Loop over data points
      s = w / const                        !Convert FWHM to standard deviation
      do i = 1,nv
            yc(i,1) = h * exp (-(x(i)-c)**2/(2D0*s**2))
            yc(i,2) = bc + bs * x(i)
            ycal(i) = yc(i,1) + yc(i,2)
      end do

      return
      end 

[Previous Page | Next Page | Top]

Example Session

We reproduce here an annotated "typical" Multi-FRILLS session using the GAUSS-FIT application shown in the previous section (cf the fitted parameters with those used in producing the simulated data). During the Multi-FRILLS session, the user is prompted for commands with '#'.

  1. The application is compiled and linked to the Multi-FRILLS library.
    
    $ @frills_library:frills_link gauss_fit
    
  2. To start the program, type
    
    $ run gauss_fit
    Multi-FRILLS version 1.3, Copyright (C) 1996, 1997 R. Osborn, T. G. Perring
    

    The program reads in the data and then enters the Multi-FRILLS interface, indicated by the command prompt #.

  3. To plot the data from the first data set, type
    
    # plot data
    
    All Multi-FRILLS command names and keywords are uniquely identified by their first letter. In this instance it is sufficient to type just P D. The plotting range can be controlled by typing additional parameters on the same line. Try:
    
    # p d 0 10 0 250
    
  4. To plot the data from the second data set, you need the USE command
    
    # u 2
    # p d
    
  5. To set the parameters for the first peak and display what we have set,
    
    # u 1	(Reset to first data set)
    # s
    Dataset 1
    ==========
    Give the parameter value (with min and max if required)
        1  Set  1 Height       : 200
        2  Set  1 FWHM         : 4
        3  Set  1 Centre       : 4
        4  Set  1 Bkgd Constant:
        5  Set  1 Bkgd Slope   :
    # d p
               Dataset 1
               ==========
               Parameter set             P             Pmin         Pmax
        1   1  Set  1 Height          200.00
        2   1  Set  1 FWHM            4.0000
        3   1  Set  1 Centre          4.0000
        4   1  Set  1 Bkgd Constant  0.00000E+00
        5   1  Set  1 Bkgd Slope     0.00000E+00
    

    Note that only the first five parameters are displayed. This is because they are the only ones associated with the first data set. In the output to the DISPLAY PARAMETER command, the first column is the parameter number and the second column is the data set with which it is associated. If any global parameters were defined in this application, these would have 0 in the second column and would be displayed with all the data sets.

  6. The function with these parameters can be viewed with PLOT CALCULATION
    
    # p c
    
  7. It is clear that the width and centre need to be adjusted. Parameters can be set individually.
    
    # s 2,3
    Give the parameter value (with min and max if required)
        2  Set  1 FWHM         : 2
        3  Set  1 Centre       : 4.5
    
  8. Now set the parameters for the second data set.
    
    # u 2
    # s
    Dataset 2
    ==========
    Give the parameter value (with min and max if required)
        6  Set  2 Height       : 100
        7  Set  2 FWHM         : 3
        8  Set  2 Centre       : 4
        9  Set  2 Bkgd Constant:
       10  Set  2 Bkgd Slope   :
    
  9. We want to constrain the centre of the second data set peak to be the same as the first. Therefore we bind parameter 8 to be equal to parameter 3.
    
    # b 8
    Give the binding parameter no. and the ratio
    Set  2 Centre       : 3 1
    # d p
               Dataset 2
               ==========
               Parameter set             P             Pmin         Pmax
        6   2  Set  2 Height          100.00
        7   2  Set  2 FWHM            3.0000
        8   2  Set  2 Centre          4.5000          Bound to parameter 3
        9   2  Set  2 Bkgd Constant  0.00000E+00
       10   2  Set  2 Bkgd Slope     0.00000E+00
    
  10. We may also decide to fix the background slope to be 0 by fixing parameters 5 and 10.
    
    # f 5,10
    
  11. We are ready to try a least-squares fit to the data by giving the command to GO.
    
    # g
    
      Iteration    1
      Old Chi**2 =       0.23542033D+02          New Chi**2 =       0.29901154D+01
    
      Iteration    2
      Old Chi**2 =       0.29901154D+01          New Chi**2 =       0.12794530D+01
    
      Iteration    3
      Old Chi**2 =       0.12794530D+01          New Chi**2 =       0.11186471D+01
    
      Iteration    4
      Old Chi**2 =       0.11186471D+01          New Chi**2 =       0.11171174D+01
    
      Iteration    5
      Old Chi**2 =       0.11171174D+01          New Chi**2 =       0.11171038D+01
    
               Job completed: Agreement factor converged
    
                           Chi ** 2   1.1171
    
               Parameter set             P            Sigma
        1   1  Set  1 Height          190.94         7.7412
        2   1  Set  1 FWHM            1.5542        0.47325E-01
        3   1  Set  1 Centre          5.0163        0.20700E-01
        4   1  Set  1 Bkgd Constant   13.875        0.78519
        5   1  Set  1 Bkgd Slope     0.00000E+00      Fixed
        6   2  Set  2 Height          103.57         4.2283
        7   2  Set  2 FWHM            2.9573        0.11575
        8   2  Set  2 Centre          5.0163          Bound to parameter 3
        9   2  Set  2 Bkgd Constant   15.225         1.1257
       10   2  Set  2 Bkgd Slope     0.00000E+00      Fixed
    
    Do you wish to restore old parameters [N]?
    
  12. We can now plot the results of the fit with PLOT FIT.
    
    # p f
    
  13. To see if varying the background slope improves the fit, we can clear parameters 5 and 10 and redo the fit. Note that we do not want to clear all parameters or we would lose the binding of the two centres.
    
    # c p 5,10
    # g
    
      Iteration    1
      Old Chi**2 =       0.11477094D+01          New Chi**2 =       0.77440763D+00
    
      Iteration    2
      Old Chi**2 =       0.77440763D+00          New Chi**2 =       0.77435676D+00
    
               Job completed: Agreement factor converged
    
                           Chi ** 2  0.77436
    
               Parameter set             P            Sigma
        1   1  Set  1 Height          191.10         6.4791
        2   1  Set  1 FWHM            1.5393        0.39569E-01
        3   1  Set  1 Centre          4.9978        0.17568E-01
        4   1  Set  1 Bkgd Constant   10.640        0.99233
        5   1  Set  1 Bkgd Slope     0.79141        0.18196
        6   2  Set  2 Height          102.95         3.5557
        7   2  Set  2 FWHM            2.9139        0.96785E-01
        8   2  Set  2 Centre          4.9978          Bound to parameter 3
        9   2  Set  2 Bkgd Constant   12.206         1.1879
       10   2  Set  2 Bkgd Slope     0.81031        0.19721
    
    Do you wish to restore old parameters [N]?
    
  14. The fit has improved so we are now ready to plot the results. In the following, we plot the fit to dataset 1 and then alter the line-type to give dashed and dotted lines (2 and 3 respectively) before over-plotting the two components to the fit; the peak and the background. The resulting plot is stored in as a postscript file in GAUSS1.PS.
    
    # u 1
    # p f
    # a l 2
    # p o 1
    # a l 3
    # p o 2
    # k gauss1.ps
    Plot stored in UD1:[OSBORN]GAUSS1.PS;1
    
  15. This process is repeated for dataset 2 before OUTPUTting the resulting parameters to an ASCII file. This file may be read in during another Multi-FRILLS session using the INPUT command. Finally, we exit from the program.
    
    # u 2
    # a l 1
    # p f
    # a l 2
    # p o 1
    # a l 3
    # p o 2
    # k gauss2.ps
    Plot stored in UD1:[OSBORN]GAUSS2.PS;2
    # o f gauss.fit
    Fitted parameters stored in UD1:[OSBORN]GAUSS.FIT;1
    # e
    

[Previous Page | Next Page | Top]

Genie Applications

Open GENIE is an interactive data analysis package developed at the ISIS Pulsed Neutron Facility, Rutherford Appleton Laboratory . It is in widespread use at ISIS, and other neutron facilities, as a means of rapid visualization and simple analysis of neutron scattering data. The neutron data are grouped into workspaces which contain the x and y-values, errors, units, titles etc of an individual spectrum. We have written a number of subroutines to ease the transfer of these workspaces into the Multi-FRILLS applications, automatically copying the data, graph labels and titles into the Multi-FRILLS common blocks.

Open GENIE supercedes GENIE v2, which would only operate on VMS machines and had a more limited command language. However, the old version of GENIE is still supported by Multi-FRILLS and will also be described here. Source code for both versions are stored in the examples subdirectory i.e. if FRILLS_LIBRARY points to DUA0:[MULTI_FRILLS], the examples will be in either DUA0:[MULTI_FRILLS.EXAMPLES.OPENGENIE] or DUA0:[MULTI_FRILLS.EXAMPLES.GENIE].

In the following sections, we discuss how to run, write and link GENIE applications. Finally, we give an example of a program that fits several GENIE workspaces at once. The source code is in the example file MULTI_RESFIT.F. The program models the peak shape in all workspaces as the convolution of a Gaussian with an exponential decay with a sloping background. A number of the functions used to define the final peak shape are included with other common peak functions in the file FRILLS_SOURCES:PEAK_FUNCTIONS.F which is also included with the Multi-FRILLS package. Other examples of GENIE applications also included are:


[Previous Page | Next Page | Top]

Running GENIE Applications

Open GENIE

To run a Multi-FRILLS application from within Open GENIE, the user should first load the GCL command file:


>> load "frills_library:frills_genie.gcl"

In order to run the program, the user types


>> <result> = frills_genie("<program name>", <workspace>)
where

The operation of the Multi-FRILLS application is identical to a standalone application once it has been started.

GENIE v2

To run a GENIE application from within GENIE v2, the user should run a GENIE command file:


>> @frills_library:frills_genie
which will prompt the user for

Generally, Multi-FRILLS programs for GENIE can act on any number of workspaces up to the maximum number for which the program was designed. The operation of the Multi-FRILLS application is identical to a standalone application once it has been started.

It is possible to write GENIE applications that explicitly act on one workspace only. These should be invoked directly by the GENIE TRANSFORM command:


>> TRANSFORM W<m> <program-name> W<n>


[Previous Page | Next Page | Top]

Writing GENIE Applications

Open GENIE

Open GENIE applications are written as subroutines which are then linked as shareable objects and loaded into an Open GENIE session. Several general purpose subroutines have been provided so that an Open GENIE application is as simple to write as a stand-alone Multi-FRILLS application. In fact, it is normally simpler, because the data, graph captions and titles are read in automatically from the GENIE workspaces.

In the following, we shall use MULTI_RESFIT.F as an example. The first lines of the application should contain the following :


      subroutine MULTI_RESFIT (pars_get, pars_put)
      include 'FRILLS.INC'
      include 'GENIE_MODULES.INC'
      external RESFUN, pars_get, pars_put
PARS_GET and PARS_PUT are used by Open GENIE to pass the workspace information to the Multi-FRILLS application. RESFUN is the name of the subroutine containing the fitting function.

It is advisable to check that the application is consistent with the version of Open GENIE being run.


      if (module_version_ok (genie_major, genie_minor) .eq. 0) then
         call module_error ("MULTI_RESFIT", "Error: Open GENIE Version mismatch", 
     +                      "Please recompile this module")   
            return
      end if

The simplest way of inputting GENIE workspaces into the application is to replace the call to MULTI_FRILLS with one to MULTI_FRILLS_OPENGENIE


      call multi_frills_opengenie (RESFUN, pars_get, pars_put)
This transfers data (assumed to be in histogram mode), including titles and graph captions, from PARS_GET to the Multi-FRILLS common blocks, and then calls MULTI_FRILLS itself. After completion of the Multi-FRILLS session it transfers the function, calculated with the latest set of parameters, into GENIE workspaces. The user need only provide the names and number of parameters, i.e.

NPTOT the total number of parameters
NAM(ip) the names of the parameters (ip=1:NPTOT)

and optionally fill the array elements INS(ip) = 0 (global) or INS(ip) = is (specific to data set is). The subroutine that calculates the fitting function should be written exactly as for standalone Multi-FRILLS applications. For those users who wish to manipulate or inspect the data before MULTI_FRILLS is called, two routines, MULTI_FRILLS_OPENGENIE_IN and MULTI_FRILLS_OPENGENIE_OUT, are provided. In this case the calling sequence would be


      call multi_frills_opengenie_in (pars_get)
                  .
                  .
      call multi_frills (FUN)
                  .
                  .
      call multi_frills_opengenie_out (FUN, pars_put)
where FUN is the name of the fitting function subroutine. One reason for using this sequence is if the number of workspaces being fitted is not known in advance. The call to MULTI_FRILLS_OPENGENIE_IN determines the value of ns after which the parameters can be initialized. This is illustrated in the example program. Otherwise the GENIE program is written exactly as with MULTI_FRILLS_OPENGENIE.

Since Open GENIE modules are subroutines rather than standalone programs, they should end with


      return
      end
rather than a STOP statement.

IMPORTANT: Don't ever put a STOP statement in an Open GENIE application. It will terminate the whole Open GENIE session, not just the Multi-FRILLS application. One unfortunate side effect of this method of invoking external procedures is that a fatal error in the procedure could prove fatal to the whole Open GENIE session. Extra care must therefore be taken in writing the fitting functions so that they do not crash. If in doubt, save essential data into an Open GENIE file before running FRILLS_GENIE.

Genie v2

GENIE v2 applications are very similar to regular Multi-FRILLS programs. Firstly, they are written as programs, rather than subroutines, which are spawned from a GENIE v2 session. In other words, the first lines would be :


      program MULTI_RESFIT
 
      include 'FRILLS.INC'
      external RESFUN
where RESFUN is the subroutine containing the fitting function. The simplest way of inputting GENIE workspaces into the application is to replace the call to MULTI_FRILLS with one to MULTI_FRILLS_GENIE

      call multi_frills_genie (RESFUN)
This reads data (assumed to be in histogram mode), including titles and graph captions, from a file written when the command file FRILLS_GENIE is invoked within GENIE. The data is read into the MULTI_FRILLS common blocks before MULTI_FRILLS itself is called within MULTI_FRILLS_GENIE. After completion of the Multi-FRILLS session it transfers the function, calculated with the latest set of parameters, into GENIE workspaces by FRILLS_GENIE. The user need only provide the names and number of parameters, i.e.

NPTOT the total number of parameters
NAM(ip) the names of the parameters (ip=1:NPTOT)

and optionally fill the array elements INS(ip) = 0 (global) or INS(ip) = is (specific to data set is). The subroutine that calculates the fitting function should be written exactly as for standalone Multi-FRILLS applications. For those users who wish to manipulate or inspect the data before MULTI_FRILLS is called, two routines, MULTI_FRILLS_GENIE_IN and MULTI_FRILLS_GENIE_OUT, are provided. In this case the calling sequence would be


      call multi_frills_genie_in
                  .
                  .
      call multi_frills (FUN)
                  .
                  .
      call multi_frills_genie_out (FUN)
where FUN is the name of the fitting function subroutine. One reason for using this sequence is if the number of workspaces being fitted is not known in advance. The call to MULTI_FRILLS_GENIE_IN determines the value of ns after which the parameters can be initialized. Otherwise the GENIE program is written exactly as with MULTI_FRILLS_GENIE.

Multi-FRILLS applications that act on one workspace only, i.e. which are invoked with the GENIE TRANSFORM command rather than the call to FRILLS_GENIE, can be written with SINGLE_FRILLS_GENIE replacing the call to MULTI_FRILLS. Alternatively, if the user wishes to manipulate the data before entering Multi-FRILLS, the following calling sequence should be used:


      call single_frills_genie_in
                  .
                  .
      call multi_frills (fun)
                  .
                  .
      call single_frills_genie_out (fun)


[Previous Page | Next Page | Top]

Linking GENIE Applications

Open GENIE

An Open GENIE application is written as a subroutine that needs to be linked as a shareable object module to Open GENIE. This is achieved by compiling and linking the source code with the command file OPENGENIE_LINK.COM e.g.


$ @FRILLS_LIBRARY:OPENGENIE_LINK <program name> [<library names>].
The module can then be loaded with the MODULE/LOAD command. However, this is done automatically by FRILLS_GENIE.GCL, the procedure invoked when running the Multi-FRILLS application.

GENIE v2

The linking procedure for GENIE v2 applications is similar to that for standalone Multi-FRILLS programs except that the GENIE library needs to be included by adding GENIE_LIB/LIB after the Multi-FRILLS library. The command file, GENIE_LINK.COM, has been written to compile and link GENIE applications.


$ @FRILLS_LIBRARY:GENIE_LINK <program name> [<library name(s)>]
As with FRILLS_LINK.COM, GENIE_LINK.COM may be edited to add more object libraries to the link statement.


[Previous Page | Top]

Example Program

The following program fits a form of the resolution function of pulsed neutron chopper spectrometers that works reasonably well in many cases. It has been written to accommodate multiple data sets although there are no global parameters. One use of such a multiple data set version would be to fit the elastic peaks from a number of different vanadium runs with the resolution parameters bound to a common value. The program is invoked by invoking the procedure FRILLS_LIBRARY:FRILLS_GENIE.GCL from within an open GENIE session. The main program calls MULTI_FRILLS_OPENGENIE_IN first in order to determine the number of data sets for which parameters need to be initialized. It then calls MULTI_FRILLS itself to transfer control to the Multi-FRILLS application. Finally, MULTI_FRILLS_OPENGENIE_OUT is called to transfer the results of the fit back to GENIE before exiting.


 MULTI_RESFIT.F
!------------------------------------------------------------------------------------
! Multi-FRILLS  - An Interactive Least Squares Fitting Package for Multiple Datasets
! 
! Copyright (C) 1996 R. Osborn, T. G. Perring (See MULTI_FRILLS.FOR for details)
!
!     Version   Date    Name         Comments
!     =======   ====    ====         ========
!       1.0   96/08/23  RO/TGP       First release of multiple dataset version of FRILLS
!       1.1   97/08/04   RO          Removed VMS logical name from include statement
!       1.2   97/09/08   TGP         Addition of Levenberg-Marquadt option
!       1.3   99/01/06  RO/TGP       First version compatible with Open GENIE
!
!------------------------------------------------------------------------------------
!     subroutine MULTI_RESFIT
!
!     Fitting of chopper spectrometer resolution function to multiple GENIE workspaces
!
      subroutine MULTI_RESFIT (pars_get, pars_put)

      include 'FRILLS.INC'
      include 'GENIE_MODULES.INC'
      external RESFUN, pars_get, pars_put
      integer*4 i, j, k

!Check that version of open Genie is OK
      if (module_version_ok (genie_major, genie_minor) .eq. 0) then
            call module_error ("MULTI_RESFIT", "Error: Open GENIE Version mismatch", 
     +                         "Please recompile this module")   
            return
      end if
! Give a welcome message
      write (*,*) '#'
      write (*,*) '# MULTI_RESFIT - Chopper Spectrometer Resolution Fitting Program'
      write (*,*) '#'
      write (*,*) '# Type H for a list of the available commands'
      write (*,*) '#'
      write (*,*) '# Details of fits are stored in SYS$SCRATCH:MULTI_RESFIT.TXT'
      write (*,*) '#'
! Read in GENIE workspaces
      call multi_frills_opengenie_in (pars_get)
! Initialize fitting parameters for each dataset
      i = 0
      do j = 1,ns
            write (nam(i+1), '(a,i2,a)') 'Set ', j, ' Bkgd Constant'
            write (nam(i+2), '(a,i2,a)') 'Set ', j, ' Bkgd Slope'
            write (nam(i+3), '(a,i2,a)') 'Set ', j, ' Intensity'
            write (nam(i+4), '(a,i2,a)') 'Set ', j, ' Gamma'
            write (nam(i+5), '(a,i2,a)') 'Set ', j, ' Sigma'
            write (nam(i+6), '(a,i2,a)') 'Set ', j, ' Centre'
            do k = i+1, i+6
                  ins(k) = j
            end do
            i = i + 6
      end do 
      nptot = i
! Output is directed to SYS$SCRATCH:MULTI_RESFIT.TXT
      open (unit = 41, file = 'SYS$SCRATCH:MULTI_RESFIT.TXT', status = 'NEW')
! Call MULTI_FRILLS_GENIE to initialize data using GENIE workspace and call Multi-FRILLS
      call multi_frills (RESFUN)
! Write out GENIE workspaces
      call multi_frills_opengenie_out(RESFUN, pars_put)
! Close the file containing the fitting details
      close (unit=41)

      return
      end

!------------------------------------------------------------------------------------
!     RESFUN
!
!     Subroutine to calculate fitting function
!
!     The function contains two components
!     1) Decaying exponential convolved with a Gaussian defining the peak shape
!     2) Sloping background
!
      subroutine RESFUN

      include 'FUNCTION.INC'
      real*8 bc, bs, int, gamma, sigma, centre, EXP_DECAY
      integer*4 i

! Definition of parameters
      i = (is-1) * 6
      bc = p(i+1)            !Background constant
      bs = p(i+2)            !Background slope
      int = p(i+3)           !Peak intensity
      gamma = p(i+4)         !Peak decay constant
      sigma = p(i+5)         !Peak Gaussian sigma
      centre = p(i+6)        !Peak centre (1st moment)
!Loop over data points defining peak function and two components
      nc = 2
      do i = 1,nv
            yc(i,1) = int * EXP_DECAY (x(i)-centre, gamma, sigma)
            yc(i,2) = bc + bs*x(i)
            ycal(i) = yc(i,1) + yc(i,2)
      end do

      return
      end 
!------------------------------------------------------------------------------------
!     EXP_DECAY
!
!     Returns value of decaying exponential convolved with a Gaussian
!
      real*8 function EXP_DECAY (eps, gamma, sigma)

      real*8 eps, gamma, sigma, GAUSS, x, ERFC, pi
      parameter (pi = 3.1415926535897932)

      if (sigma .le. 0.0 .and. gamma .le. 0.0) then
            EXP_DECAY = 0.0
      else if (gamma .le. 0.0) then
            EXP_DECAY = GAUSS (eps, sigma)
      else if (sigma .le. 0.0) then
            x = eps - gamma
            if (x .ge. 0.0) then
                  EXP_DECAY = 0.0
            else if (x .lt. 50.0D0*gamma) then
                  EXP_DECAY = exp (x / gamma) / gamma
            else
                  EXP_DECAY = 0.0
            end if
      else
            x = eps - gamma
            if (x .lt. 5.0*sigma) then
                  EXP_DECAY = exp(0.5D0*sigma**2/gamma**2 + x/gamma) * ERFC((sigma/gamma + x/sigma) / sqrt(2D0)) / (2D0*gamma)
            else
                  EXP_DECAY = 0.0
            end if
      end if

      return
      end 
!------------------------------------------------------------------------------------
!     GAUSS
!
!     Returns value of normalized Gaussian centred at 0 with standard deviation SIGMA
!
      real*8 function GAUSS (x, sigma)

      real*8 x, sigma, pi
      parameter (pi = 3.1415926535897932)

      GAUSS = exp(-x**2/(2D0*sigma**2)) / (sigma * sqrt(2D0*pi))

      return
      end
!------------------------------------------------------------------------------------
!     ERFC
!
!     Complementary Error Function
!     Algorithm from Abramovitz and Stegun p.299
!     Accurate to 3E-07 for both +ve and -ve x.
!
      real*8 function ERFC (x)

      real*8 a1, a2, a3, a4, a5, a6, x, z
      data a1, a2, a3 /0.0705230784, 0.0422820123, 0.0092705272/
      data a4, a5, a6 /0.0001520143, 0.0002765672, 0.0000430638/

      z = abs(x)
      ERFC = 1.0 + a1*z + a2*z**2 + a3*z**3 + a4*z**4 + a5*z**5 + a6*z**6

      if (ERFC .ge. 10.0) then
            ERFC = 0.0
      else
            ERFC = 1.0 / ERFC**16
      end if
      if (x .lt. 0.0) ERFC = 2.0 - ERFC

      return
      end

[Previous Page | Top]

Comments to: Ray Osborn <ROsborn@anl.gov>
Revised: Wednesday, October 27, 1999

Copyright © 1999, Ray Osborn and Toby Perring. All rights reserved. [Legal Disclaimer]