xgaufit.src

        Copyright (c) Kapteyn Laboratorium Groningen 2010
        All Rights Reserved.

Name:         xgaufit.src
Creator:      terlouw
Host:         azarquiel
Date:         Jun 28, 2010
Contents:     xgaufit.make xgaufit.dc1 xgaufit.h gui.h dofit.h
              utils.h fitopt.h import.h filtopt.h param.h
              parcol.h outputresults.h xgaufit.c gui.c dofit.c
              utils.c allocs.c fitopt.c import.c filtopt.c
              param.c parcol.c outputresults.c

#>            xgaufit.make
#----------------------------------------------------------------------
# Makefile for XGAUFIT, jun 27, 1999
#
#
# Pack sources into render.src with:     > make -f xgaufit.make pack
# Clean executable and objects           > make -f xgaufit.make clean
# Unpack source files with:              > $gip_exe/xfile xgaufit.src
# Install source file with:              > p -reserve xgaufit.src
#                                        > p -install xgaufit.src
#
# Local compilation for testing:         > p xgaufit.make
#
# The CC_OPTS etc. are set by 'compile'
#----------------------------------------------------------------------

SHELL = /bin/sh

OPTS  = $(CC_OPTS)
COMP  = $(CC_COMP)
LIBS  = $(CC_LIBS)

.SUFFIXES: .c .o

DOCUMENTS = xgaufit.dc1

INCLUDES  = xgaufit.h gui.h dofit.h utils.h  fitopt.h import.h filtopt.h param.h parcol.h outputresults.h

SOURCES   = xgaufit.c gui.c dofit.c utils.c allocs.c fitopt.c import.c filtopt.c param.c parcol.c outputresults.c

OBJECTS   = xgaufit.o gui.o dofit.o utils.o allocs.o fitopt.o import.o filtopt.o param.o parcol.o outputresults.o


default:: xgaufit

pack::
	$${gip_sys}/pack.csh xgaufit.src xgaufit.make \
        $(DOCUMENTS) $(INCLUDES) $(SOURCES) $(PGINCLUDE)

clean::
	rm -f xgaufit $(OBJECTS)

xgaufit.o: xgaufit.h

dofit.o:    xgaufit.h

allocs.o:   xgaufit.h

gui.o: gui.h


.c.o:  
	@echo "$(CC_COMP)  -c $(CC_OPTS) $(X11_OPTS) $< ";\
	$(CC_COMP) -c $(CC_OPTS) $(X11_OPTS) $<

xgaufit: $(OBJECTS)
	@echo "$(CC_COMP) -o xgaufit $(OBJECTS) $(X11_LIBS) $(CC_LIBS) $(XT_LIBS)";\
	$(CC_COMP) -o xgaufit $(OBJECTS) $(X11_LIBS) $(CC_LIBS) $(XT_LIBS)

#<

#>            xgaufit.dc1
                            COPYRIGHT (c) 1999
                      Kapteyn Astronomical Institute
                University of Groningen, The Netherlands
                           All Rights Reserved.    

Program:       XGAUFIT

Purpose:       XGAUFIT estimates parameters for multi component one
               dimensional gauss functions using data in selected
               profiles. It can use these estimates to fit the
               parameters to the data in a least squares fit.
               Output is send to subsets in a set.

Category:      ANALYSIS, PROFILES, FITTING

File:          xgaufit.src

Author:        M.G.R. Vogelaar
               (GUI: J.P. Terlouw)

Keywords:      This is a list of keywords from an example run of the 
               application. It should be obvious to which gui elements 
               these keywords correspond.
                
               AMPPERC= 50
               AXUNITS= KM/S
               BOXAXES= RA DEC
               CRITAMP= 0
               CRITDISP= 0
               DATAUNITS= W.U.
               FUNCTION= 0
               GAUSET=fittedprofiles
               GRIDRANGE= -30 30
               LAB= 0.01
               MAXITS= 200
               NCOMP= 1
               OUTSET=fittedparams
               PARAM= 
               PHYSRANGE= 552.668 1548.98
               PROFILE= freq
               PROFINDX= 
               Q= 1 2 3 4 5
               RMS= 1
               SETNAME= real
               TOLERANCE= 0.001
               VSYS= 1050
               Z0_EST= 0
               Z0_FIX= YES
               Z1_EST= 0
               Z1_FIX= YES
               Z2_EST= 0
               Z2_FIX= YES
               
               
Buttons etc.: 
               
               Notes about the Graphical User Interface (GUI):
               
               User input in this application is coupled to the Hermes keyword 
               mechanism. Keywords can still be specified on the Hermes command 
               line and also the syntax for input of numbers and text is the 
               same as for keywords. Also, Emacs-like text editing in input 
               fields is possible and we added mouse functionality. 
               E.g. the right mouse button will erase text in an input field. 
               With the middle button you can paste text stored in your mouse 
               buffer.

               <Esc>, <right mouse button> and ctrl-U
                                           Wipe current text from input field
               ctrl-A                      Go to beginning of field
               ctrl-E                      Go to end of field
               ctrl-K                      Wipe from cursor position to end of 
                                           line (ctrl-Y does NOT yank back)
               <TAB>                       Go to next input field.


               Control keys in the Help window:
               ctrl-S                      Search forward
               ctrl-R                      Search backwards
               ctrl-V                      Scroll text one page down
                                           (same as Page-Up)
               ctrl-P                      Scroll text one line up
               ctrl-Z                      Scroll text one line down
               

               Notes:

               1) Do not close pop-up windows using the cross symbol 
               on the upper windows bar. This will crash your application. 
               Instead use the close button or repeat pressing the button
               in the main gui with which you opened your pop-up window. 

               2) On some systems, a menu item is un-selectable if the 
               Caps-Lock is on!
               


               MAIN WINDOW:



               Setname: input field.
               
               Enter name of set which contains your profile data.

               
               Profile axis: Menu
               
               List of possible profile directions.
               

               Grid range: Input field.
               
               Default is the entire profile length in grids. 
               If you change one of the values, the corresponding 
               physical range is altered too.
               

               Phys range: Input field.
               
               Default is the entire profile length in physical coordinates.
               If you change one of the values, the corresponding 
               grid range is altered too.
               

               Box: Input field.

               The limits in your set between which profiles are examined.
               A box consist of two vectors. The first gives the low limits and
               the second gives the high limits. The dimension of
               a vector is the set dimension - 1. If you press <enter>
               you will get a default box which includes all profiles in 
               the set.
               

               Vsys: Input field

               Value used in sorting when you want to sort in distance
               to a position along the profile axis. Usually this will be
               the systemic velocity of a galaxy.
               
               
               Parameter output: Input field
               
               Name of the output set on disk which holds all the fitted
               parameters in subsets. Size of the output set is determined
               by the program. 
               
               
               Fitted data output: Input field
    
               Name of the output set on disk which holds the fitted data.
               The output set is copied from the input set. Initially all
               data is set to blank. You can get a residual cube by subtracting
               this set from the input set e.g. with task SUB.
                              
              
               GO/STOP: Button
               
               Start calculating initial estimates and fitting the 
               parameters to the data in the profiles. If you want to
               abort this process, press this button again and wait a while 
               for the program to respond.
                                                       

               BUTTONS in button bar:
               
               FILE
               ====
               
               Menu options:
               1) Get estimates from these subsets instead of calculating them.
               You must give the name of the set and a subset specification.
               E.g. myestset param 1 2 3.
               You need as much subsets as components x 3.
               2) exit program


               FILTERS
               =======
               
               Amp. range: Input field
               
               No values means no filtering. One value filters out all 
               amplitudes below this value and two values filter out all 
               amplitudes below the first and all amplitudes above the 
               second value.
               
               Vel. range: Input field 
               Disp.range: Input field 

               See above. These last input fields need values in physical 
               coordinates. 
                                         
               
               
               FIT OPTIONS
               ===========
               
               Tolerance: Input field 
               
               Stop criterion for least squares fit (see description below).
               
               Lambda: Input field
               
               Mixing parameter (see description below)
               
               Max. its: Input field
               
               If the stop criterion is not reached after this number 
               of least squares iterations, then abort the fitting of the
               current profile. Write blanks in the output set in the 
               corresponding position.
               
               RMS: Input field

               Set the value with which the value of chi-squared is scaled to 
               get a value for the reduced chi-squared that is related to
               the noise in your profile. The default is 1.
 
               Min. ampl: Input field 
               
               Discard initial estimates with an amplitude lower than this
               value.
               
               Min. disp: Input field 
               
               Reject initial estimates with a dispersion lower than this
               value.
               
               Q: Input field
               
               Smoothing factors (see description below)
               

               PARAM
               =====
               
               Set fit parameters to fixed or free in the fit. If you want
               to set a parameter to fixed, you have to enter a value
               for this parameter.
               For the background terms, an initial estimate must be 
               given in the input fields. The defaults for all the 
               background terms are zero.

              
               FUNCTION
               ========
               
               Menu with the following functions (see description below)
               Gauss
               Gauss Hermite (h3)
               Gauss Hermite (h3,h4)               
               Voigt
               


               NCOMP
               =====
               
               Menu with possible number of components per function.
               

               SORT
               ====
               
               Menu with sort options:
               None
               Peak,
               Disp,
               Vsys
               
               Select sort option if you have selected more than one 
               component to fit and want to sort these components. 
               See also see description below.


               HELP
               ====
               
               Pop up this dc1 document




Description:   BASIC INPUT

               The decomposition of profiles into Gaussian components
               can be a powerful method in the study of the kinematics
               of neutral hydrogen. The name of the set for which you
               want to study your profiles is entered in the main gui
               of the program. The input field is called 'Setname'.
               
               If a set is known, then the profile direction is selected
               from a menu labeled 'Profile axis'. Usually this will 
               be the FREQ or VEL axis of your data set. Note that you can
               select a profile direction as long as your input set has
               more than one axis so that input is not restricted to
               a three dimensional data cube. The length of your 
               profile axis must be greater than 1 and less than 4096.
               
               The length of a profile is set in the input fields
               'Grid range' and 'Phys. range'. A default is displayed
               after the name of a set and the profile axis are known.
               You can change these values either in grids or physical
               coordinates. If you change one field you see that the 
               values in the other field will change also.
               
HINT:          Inspect with XGAUPROF whether you can 
               exclude unwanted data by decreasing the grid range. However 
               if your grid range is too small you can run into trouble 
               while fitting base lines.
               
               The 'Box' field is blank. Pressing <enter> results in
               a default box which is the entire data range in grids.
               Note that the box can be given in any coordinates as you
               were used to in the non-gui GIPSY tasks.
                                             
               In general, the program will work on profiles in any
               direction and instead of velocity you could think of any
               kind of physical coordinate. The program automatically
               inserts the correct units in the header of the
               output set.

               
               PREPARING THE FIT
               
               Estimates calculated by the program
               -----------------------------------
               
               You can fit the parameters of four different functions
               to the profile data. The results are written to a
               GIPSY output set on disk, each parameter in a different
               subset per profile position. But the least squares fit
               needs a good starting point for these parameters and 
               an important part of this program is dedicated to finding 
               good estimates  by fitting second degree polynomials to the 
               profiles (See: Schwarz, 1968, Bull.Astr.Inst. Netherlands,
               Volume 19, 405-413). An important role here play the so 
               called smoothing parameter 'Q'. The initial estimate
               routine uses 2Q+1 points around a maximum to fit a second
               order polynomial. The coefficient of the second-order term
               of the polynomial is an approximation of the second derivative
               of the observed profile. With a moments method
               the center and dispersion are estimated. The maximum amplitude
               is derived from the observed profile.

               In practice there
               is not one smoothing factor that gives the best (in terms 
               of chi-squared) estimates for all profiles. Therefore you
               can enter in the 'FIT OPTIONS' menu a series of smoothing
               factors ('Q'), which are all used in all profiles to get the 
               best estimates. You can exclude estimates by entering
               values for a minimum amplitude ('Min.ampl') and the
               minimum dispersion ('Min.disp'), but usually you leave
               these fields blank and set a filter on the least squares fit 
               results.

               Note that there is one special case build in the program.
               If a profile contains blanks or zero-level values except
               for one, then the estimate routine fails to return an
               estimate. One can treat this peak as a peak with
               known amplitude and center, and a full width at half maximum
               of one pixel (which is automatically converted to a dispersion).
               However, the least squares fit routine will not fit such 
               profiles because it converges to zero dispersion.
                                                    


               Estimates from an input set
               ---------------------------

               Alternatively it is also possible to read estimates from 
               a set. This set could be an output set of a previous run of
               XGAUFIT. If not, then you must be sure that the box of
               your estimates set matches the box of your input set
               and that the data has the same physical coordinates as
               the set with the profiles. 
               Not only the set name must be entered, but also the subsets,
               e.g.:  myestimates param 1 2 3
               
               With this option it is possible to separate the data with which 
               you create estimates, and create fits. For instance you can 
               create estimates with clipped or conditional transferred
               data and use these estimates to fit parameters using the
               original profile data in a second step. Such methods avoid
               the noisy parts of your data and result in 'cleaner' output.


               THE FIT PROCEDURE
               -----------------
               
               The actual fitting is a least-squares fit of a function
               to a set of data points. The method used is described
               in: Marquardt, J.Soc.Ind.Appl.Math. 11, 431 (1963).
               It is a mixture of the steepest descent method and the
               Taylor method.  

               
               Weights
               -------
               
               The fit routine uses equal weights for all data that
               is not blank. Blank values get weight zero.
               

               Tolerance
               ---------

               The fit is controlled by three parameters which are set
               in the 'FIT OPTIONS' window. Most important parameter is
               the stopping criterion 'Tolerance'.               
               Fitting is stopped when successive iterations fail
               to produce a decrement in reduced chi-squared less
               than TOLERANCE=  The value cannot be less than a
               certain minimum as set by the system. This means that
               maximum accuracy can be obtained with TOLERANCE=0.0

               Usually only signal is present in a small region in the
               profile while the chi square is calculated over the entire
               baseline. In this case the Chi square primarily represents
               the fitting of the baseline instead of the gauss and then
               a very tiny decrease of the Chi square might very well be
               significant for the fitting of the gauss to the signal.
               So one should be very careful not setting TOLERANCE= too
               high. If TOLERANCE=0 then the minimization is stopped at
               machine precision. Increasing the tolerance results in 
               a shorter process time. If the tolerance cannot be reached
               within the maximum number of iterations, the fit failed and
               blanks are written in the output set.
HINT:          In most cases the default tolerance 0.001 gives fast 
               performance and good fit results.
              
              
               Mixing factor
               ------------

               Lambda in the 'FIT OPTIONS' window is the mixing parameter in 
               the least squares fit routine.
               Lambda determines the initial weight of steepest descent
               method relative to the Taylor method. It should be
               a small value (i.e. 0.001).
               

               Maximum number of iterations
               ----------------------------               

               Fitting is also controlled by the maximum number of
               iterations 'Max. Its.' in the 'FIT OPTIONS' window. If the
               maximum number of iterations is exceeded, then the fit
               failed and blanks are written in the output set.


               Fixed or free parameters               
               ------------------------
               
               In the 'PARAM' window you can alter the status of the 
               fit parameters. Default, all function related parameters
               are 'free' parameters (i.e. to be fitted) and the 
               background parameters (a constant, a linear and a quadratic
               term) are fixed at a value 0.0.
                              
               
               Error on fitted parameters
               --------------------------
               
               The errors on the fitted parameters are calculated as
               real errors by means of the covariance matrix in the fit
               and the average deviation between data and fit over the
               entire profile. However, if the region with signal is
               short compared to the full length of the profile, the
               average deviation between data and fit is primarily
               determined along the part of the baseline where no signal
               is present and then it becomes almost equal to the noise
               in the profile and the real errors determined in this way
               virtually become formal errors.


               Note that if Hanning smoothing has been applied, the subsets 
               are two by two correlated. The errors on the parameters
               are then a factor sqrt(2) higher.



               MULTI COMPONENTS AND SORTING
               ----------------------------

               If you fitted parameters of multicomponent functions, 
               then it could be useful to sort the fitted values. Usually 
               you will sort components in INCREASING order of distance
               to a user given value on the profile axis (often the
               systemic velocity) as given in Vsys in the main window.
               The sort menu is found in the upper right area of the 
               main gui. Other sorting options are: sort components in 
               DECREASING order of peak value and sort components in 
               DECREASING order of dispersion


               FILTERS
               -------
               
               The program has three filters for the output. You can 
               enter ranges for amplitude, velocity (or whatever is
               on the profile axis) and dispersion. Fitted parameters 
               outside these ranges are filtered out, i.e. they are set to 
               blank in the output. If no values are entered then nothing
               will be filtered. If only one value is entered for one
               of the filters, then everything BELOW this value is 
               filtered out. In the gui log you can monitor how many profile 
               fits are rejected by these filters.
               
Warning:       The velocity filter is straightforward for standard Gauss 
               functions. If one fits parameters of the Gauss-Hermite
               series, then the filter applies to parameter b which is 
               only the mean (velocity) of the profile distribution if the 
               parameters h3 and h4 are zero. 
               
               
              
               FUNCTIONS
               ---------

               The function menu provides the choice of four functions
               of which you can fit the parameters to the data in your
               profiles.               

               1) A standard Gaussian. Parameters are amplitude, center,
                  and dispersion,
               2) Gauss-Hermite polynomial (h3). Parameters are a,b,c,
                  (which are NOT the same as amplitude, center and
                  dispersion) and h3. The parameters are translated to
                  the more familiar amplitude, center, dispersion and
                  skewness. All parameters mentioned above are saved
                  in the output set.
               3) Gauss-Hermite polynomial (h3, h4). Same as above, but 
                  an extra parameter h4 is included. The value for this
                  h4 parameter is translated to kurtosis.
                  Note that if one of h3 and h4 is not zero, the mean of 
                  the distribution is not the position of the maximum. 
                  The program calculates this parameter separately.

                  (Reference; Marel, P. van der, Franx, M., A new method 
                  for the identification of non-gaussian line profiles 
                  in elliptical galaxies. A.J., 407 525-539, 1993 April 20).

               4) Voigt. Fit Lorentzian and Doppler widths together 
                  with center and amplitude. 
                  
                  The mathematics behind all this is worked out
                  in a separate document:
                  (http://www.astro.rug.nl/~gipsy/xgauprof/index.html)
              

               OUTPUT
               ------
               
               A name is sufficient to specify the output set.
               If the set exists, you are asked to confirm to write
               it over. Before writing any data, an existing set will 
               always be deleted first to avoid complaints about
               incompatible sizes if you change boxes or grid ranges and 
               do not alter the name of the output set.

               The program creates the new set on disk and copies relevant 
               header items. It replaces the profile axis
               by a new axis called PARAM-XGAUFIT.  Its length
               depends on the maximum number of components to be fitted and
               the selected function. It is also possible to save the 
               estimates only (Est.only button in 'FIT OPTIONS' window).
               The Hermes log displays the structure of your output.
               This information, together with relevant units, is also 
               written to the comment fields header of the output
               set so that you can always retrieve subset information
               e.g.:
               HISTORY ITEM=C MODE=R INSET=myout
               or:
               FIXHED INSET=myout ITEM=
               which shows header items at set level or 
               FIXHED INSET=myout p ITEM=
               to show units , minimum, maximum and number of blanks at
               subset level.

               Not only the fitted parameters are saved,
               but also the errors on these parameters, the background 
               parameters, total flux, number of iterations, reduced
               Chi-squared of the fit, the standard deviation of the residuals
               and the smoothing factor 'Q' used by the estimate routine.

               After finishing the fitting process, the program writes
               a description of the subsets to the comment fields of the
               header of the output set. The units, minimum and maximum values
               and the number of blanks are written on subset level in 
               the FITS keywords BUNIT, DATAMAX, DATAMIN and NBLANK.


               Output of the fitted profiles
               -----------------------------
               
               There is an option to reconstruct your data using the fitted
               parameters. These 'fitted profiles' are written to a set 
               which is a copy of the input set, i.e. the header is a copy
               and the data is set to 'blank'. The size of the output is
               limited by the values for the box and the profile range
               (Grid range or Phys. range). The header keywords for the
               minimum and maximum data value and the number of blanks are
               NOT updated. This update can be done with task MNMX.
                              
               Each fitted profile is written to this output set. If the fit
               was not successful a blank appears in all the profile subsets
               of the output set.
               
               A message is appended to the HISTORY of the output set.
               The message contains the task name and time/date.

              
Updates:       Jul 20, 1999: VOG, Document created.
               Jul 27, 2000: VOG, g->Yfit was not initialized to NULL.
                                  This seemed only to cause problems on PC's(?)
               Nov  6, 2001: JPT, Wait for processing finished before quitting.
               Apr 15, 2009: VOG, Removed macro NINT
               Jun 28, 2010: JPT, Changed log font.
#<

#>            xgaufit.h
/*
                           COPYRIGHT (c) 1999
                      Kapteyn Astronomical Institute
                University of Groningen, The Netherlands
                           All Rights Reserved.
*/


/*  xgaufit.c: include files     */

#include    "stdio.h"        /* Defines ANSI C input and output utilities */
#include    "stdlib.h"       /* Defines the ANSI C functions for number */
                             /* conversion, storage allocation, and similar tasks.*/
#include    "string.h"       /* Declares the ANSI C string functions*/
                             /* like:strcpy, strcat etc.*/
#include    "math.h"         /* Declares the mathematical functions and macros.*/
#include    "gipsyc.h"       /* Defines the ANSI-F77 types for Fortran to C intface */
                             /* including def. of char2str,str2char,tofchar,zadd */
                             /* and macros tobool and toflog */
#include    "float.h"        /* Definition of FLT_MAX etc.*/
#include    "ctype.h"        /* Declares ANSI C functions for testing characters */
                             /* like: isalpha, isdigit etc. also tolower, toupper.*/


/* Common includes */

#include    "init.h"         /* Declare task running to HERMES and initialize.*/
#include    "finis.h"        /* Informs HERMES that servant quits and cleans up the mess.*/
#include    "anyout.h"       /* General character output routine for GIPSY programs.*/
#include    "setfblank.h"    /* Subroutine to set a data value to the universal BLANK.*/
#include    "error.h"        /* User error handling routine. */
#include    "myname.h"       /* Obtain the name under which a GIPSY task is being run.*/
#include    "nelc.h"         /* Characters in F-string discarding trailing blanks.*/
#include    "minmax1.h"      /* Min, max of an array */
#include    "axunit.h"       /* Returns the units of the physical coordinates of an axis in a set */
#include    "cotrans.h"      /* Transformation from grid coordinates to physical coordinates vv. */

/* User input routines */

#include    "userfio.h"      /* Easy-C companions for user interface routines.*/
#include    "userint.h"      /* User input interface routines.*/
#include    "userlog.h"
#include    "userreal.h"
#include    "userdble.h"
#include    "usertext.h"
#include    "userchar.h"
#include    "usercharu.h"
#include    "userfio.h"      /* Easy-C companions of userxxx_c */
#include    "reject.h"       /* Reject user input.*/
#include    "cancel.h"       /* Remove user input from table maintained by HERMES.*/
#include    "keyevents.h"    /* Definitions for event-driven operation */
#include    "events.h"       /* Describes routines for event-driven programs.*/
#include    "eventmonitor.h" 
#include    "wkey.h"
#include    "deputy.h"
#include    "subst.h"
#include    "status.h"
#include    "setdblank.h"


/* Input of sets */

#include    "gdsinp.h"       /* Input of set, subsets, return # subsets.*/
#include    "gdsout.h"       /* Set for output */
#include    "gdspos.h"       /* Define a position in a subset.*/
#include    "gdsbox.h"       /* Define a box inside/around a subset.*/
#include    "gdsasn.h"       /* Copy input set */
#include    "gdscpa.h"       /* Change properties of an axis in output set */
#include    "gdscss.h"       /* Change size of output before gdsout call */
#include    "gds_exist.h"
#include    "gds_close.h"    /* Close set. */
#include    "gdsc_range.h"   /* Return lower left and upper right corner of a subset.*/
#include    "gdsc_word.h"
#include    "gdsc_ndims.h"   /* Return the dimensionality of a coordinate word.*/
#include    "gdsc_grid.h"    /* Extract grid value.*/
#include    "gdsc_fill.h"    /* return coordinate word filled with a grid */
                             /* value for each axis.*/
#include    "gdsi_read.h"    /* Reads data from (part of) a set.*/
#include    "gdsi_write.h"
#include    "gdsc_name.h"
#include    "gdsd_rchar.h"
#include    "gds_delete.h"


#include    "gui.h" 
#include    "ggi.h"
#include    "utils.h"


#define MYMAX(a,b)     ( (a) > (b) ? (a) : (b) )
#define MYMIN(a,b)     ( (a) > (b) ? (b) : (a) )
#define BETWEEN(a,b,c) ( (a) < (b) ? (b) : ((a) > (c) ? (c) : (a) ) )
#define ABS(a)         ( (a) < 0 ? (-(a)) : (a) )
#define PI             3.141592653589793
#define RAD(a)         ( (a) * 0.017453292519943295769237 )
#define DEG(a)         ( (a) * 57.295779513082320876798155 )

#define RELEASE        "2.1"      /* Version number */
#define MAXAXES        10         /* Max. axes in a set */
#define MAXSUBSETS     4096       /* Max. allowed subsets */
#define MAXSUBOUT      40
#define MAXBUF         32*1024    /* Buffer size for I/O */
#define STRLEN         256        /* Max length of strings */
#define FILENAMELEN    256        /* Max length of file names */
#define FITSLEN        20         /* Max length of header items etc.*/
#define MAXPROFILES    10000      /* Max number of profiles to process at once */
#define NOERRORLOOP    100        /* Special code for gdsinp */
#define NONE           0          /* Default levels in userxxx routines */
#define REQUEST        1
#define HIDDEN         2
#define EXACT          4
#define YES            1          /* C versions of .TRUE. and .FALSE. */
#define NO             0
#define FIRSTLINE      0
#define SECONDLINE     1
#define DEBUG          16
#define MAXFUNCTIONS   4          /* Gauss, Skewed Gauss, Voigt */
#define FIEID          100        /* Code to identify one of the functions */ 

#define INIT           1          /* Flag for gui_log initialization */
#define UPDATE         0

#define MAXQ           16
#define MAXGAUSS       8          /* Max. gaussians in estimates routine */
#define MAXPAR         5          /* Max. parameters in one function */
#define MAXPARALL      MAXPAR*MAXGAUSS+3  /* Max. parameters in all gaussians */
                                  /* + 3 for background Z0+Z1x+Z2x^2 */
                                  
#define MAXLOGLINES    21                                 
#define GAUSS          0
#define GAUSSHERMITE1  1
#define GAUSSHERMITE2  2
#define VOIGT          3

#define MAXITERS       200
#define TOLERANCE      0.001
#define LAMBDA         0.01

#define NEWAXISNAME   "PARAM-XGAUFIT"   /* Profile axis is replaced by this in the output */


typedef struct                    /* Parameters for the gui log */
{
   int      noestimates;
   int      iterstoosmall;
   int      nofit;
   int      ampfiltered;
   int      velfiltered;
   int      dispfiltered;      
   double   rchimax;
}
logtype;


typedef struct                    /* Parameters for the user filters */
{
   double  amp[2];
   double  disp[2];
   double  vel[2];
}
filtertype;



typedef struct
{
   int      setexist;
   int      boxexist;
   fchar    Setin;                         /* Name of 2-or 3 dim. input set */
   char     setbuf[STRLEN];
   char     setname[STRLEN];
   int      setnameexist;
   fint     subin[MAXSUBSETS];             /* Subset coordinate words */
   fint     axnum[MAXAXES];
   fint     subdim;                        /* Dimensionality of the subsets */
   fint     axesoutsidesub;
   fchar    Setout;
   char     setoutbuf[STRLEN];
   fint     subout[MAXSUBSETS]; 
   fint     axnumout[MAXAXES];
   int      outsetexist;
   
   int      gausetexist;                   /* The fiited data, not parameters */
   fchar    Setgau;
   char     setgaubuf[STRLEN];
   fint     subgau[MAXSUBSETS];       
    
   fchar    Setest;
   char     setestbuf[STRLEN];
   fint     subest[MAXSUBSETS]; 
   int      altset;
   int      naltsubs;
      
   
   fint     profilelen;
   fint     numprofiles;
   double   **profiles;
   double   **results;
   double   **estfromset;

   int      profileexist;
   char     profaxname[FITSLEN];
   fint     profcw;                        /* Coordinate word of profile */
   fint     profaxnum;                     /* Axis number of profile in set */
   int      currentprofilenumber;   

   fint     blo[MAXAXES];                  /* Low  edge of box in grids */
   fint     bhi[MAXAXES];                  /* High edge of box in grids */
   fint     flo[MAXAXES];                  /* Low  edge of frame in grids */
   fint     fhi[MAXAXES];                  /* High edge of frame in grids */   
   fint     modvec[MAXAXES];               /* Index to vector help variables */

   float    datamin;
   float    datamax;
   fint     imagesize;
   double   *xgrids;
   double   *xphys;
   double   *spline;
   double   *median;
   double   *xphysscale;
   double   *Yres;                         /* Profile residuals after fit */
   float    *Yfit;                         /* Reconstructed (fitted) profile */
   double   gridspac;
   fint     plohi[2];                      /* Range of profile in grids */
   fint     profgridlo;
   fint     profgridhi;   
   char     dataunits[FITSLEN+1];
   char     axunits[FITSLEN+1];
   int      tophys;                        /* Is conversion grid -> phys. possible? */
   double   xgridminmax[2];
   char     axisnames[MAXAXES*FITSLEN];   
   int      inactive;                      /* Global flag to return immediately */
   int      rangeexist;
   int      count;                         /* Global profile counter */

   int      nvar;                          /* Variable parameters per function */
   double   tolerance;
   double   lab;
   int      maxits;
   fint     smoothingpar[MAXQ];
   fint     numq;                          /* Number of smoothing parameters */
   fint     iters;
   double   chi2;
   double   rchi2;
   double   sigmaresidual;
   double   rms;
   double   vsys;
   double   critamp;
   double   critdisp;
   double   lorentz;
   double   h3;
   double   h4;
   int      weights;
   int      ncomp;
   int      sortopt;  
   fint     gaussiansfound;
   float    usedQ;
   double   estimates[MAXPARALL];          /* MAXGAUSS*MAXPAR + zero level */
   double   fixedvalues[MAXPARALL];        /* These are used as fixed estimates */
   double   errorpar[MAXPARALL];           
   fint     fixedfree[MAXPARALL];   
   bool     estonly;
   double   *wdat;                         /* Weights for lsqfit function */
   double   *work;                         /* Work array for sorted profile */
   double   *filtered;
   bool     estfilter;
   double   ampperc;                       /* Calculate width at this percentage */
   int      function;                      /* Which function must be fitted? */
   fint     noutsubs;                      /* # subsets in output set */
   int      fitpars;                       /* Number of output parameters */
   int      errpars;                       /* And their errors */
   
   logtype      log;
   filtertype   filter;
} globaltype;



extern int  allocs( globaltype *, int );
extern void dminmax( double *, fint *, double *, double * );
extern void helpmenu( char** names );
#<

#>            gui.h
#if !defined(_gui_h_)
#define _gui_h_
void gui(void);                             /* activate GUI */
void gui_status(char *message, int lineno); /* display status message */
void gui_log(char *message, int lineno);    /* display log message */
void gui_error(char *message);              /* display error message */
void gui_axnames(char **names);             /* change axis menu entries */
#endif
#<

#>            dofit.h
/*
                           COPYRIGHT (c) 1999
                      Kapteyn Astronomical Institute
                University of Groningen, The Netherlands
                           All Rights Reserved.
*/


extern void dofit( globaltype * );
extern void setuserguess( globaltype * );
#<

#>            utils.h
/*
                           COPYRIGHT (c) 1999
                      Kapteyn Astronomical Institute
                University of Groningen, The Netherlands
                           All Rights Reserved.
*/


extern void clearstr( fchar );
extern void hermesstatus( char * );
extern void error_status( char *, char *, char * );
extern void clearstatus( void );
extern void dminmax( double *, fint *, double *, double * );
extern int grid2phys2( fchar, fint, fint, int, double * );
extern void scalearray( double  *, double  *, int ); 
extern double **dmatrix( int, int, int, int );  
extern void freedmatrix( double **, int, int );
extern float **fmatrix( int, int, int, int );  
extern void freefmatrix( float **, int, int );
extern bool UserLog(char *);
extern void getdble( ident, char  *, int, void  *);  
extern float tofloat( double );
#<

#>            fitopt.h
#if !defined(_fitopt_h_)
#define _fitopt_h_
ident Fitopt(char*);
#endif

#<

#>            import.h
#if !defined(_import_h_)
#define _import_h_
void Import(char *key);
#endif
#<

#>            filtopt.h
#if !defined(_filtopt_h_)
#define _filtopt_h_
ident Filtopt(char*);
#endif

#<

#>            param.h
#if !defined(_param_h_)
#define _param_h_
void Param(char *key);
#endif
#<

#>            parcol.h
#if !defined(_parcol_h_)
#define _parcol_h_
ident BuildColumn(int nrows, char *title, char *parname);
#endif
#<

#>            outputresults.h
#if !defined(_outputres_h_)
#define _outputres_h_
void outputresults( globaltype * );
void w( double, double, double *, double * );
double funcd_c( double *, double *, fint *, fint * );
void disp2fwhm( double * );
void fwhm2disp( double * );
#endif
#<

#>            xgaufit.c
/* xgaufit.c -XT
                            COPYRIGHT (c) 1999
                      Kapteyn Astronomical Institute
                University of Groningen, The Netherlands
                           All Rights Reserved.

*/

#include    "xgaufit.h"
#include    "cmain.h"        /* Defines the main body of a C program with */
                             /* MAIN_PROGRAM_ENTRY and IDENTIFICATION */
#include    "gui.h"
#include    "dofit.h"
#include    "utils.h"
#include    "time.h"
#include    "timer.h"
#include    "gdsd_writec.h"
#include    "gdsd_delete.h"
#include    "gdsd_wchar.h"
#include    "gdsd_wvar.h"
#include    "minmax3.h"
#include    "wminmax.h"
#include    "listctrl.h"


#define      KEY_OUTSET     tofchar("OUTSET=")
#define      KEY_GAUSET     tofchar("GAUSET=")



fint         setlevel;       /* To get header items at set level. */
float        fblank;         /* Global value for BLANK. */
double       dblank;         /* Double precision version */
double       maxarg;
char         taskname[20];




static void getint( ident  id,
                    char  *key,
                    int    code,
                    void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: Keyword handler for (f)ints.                      */
/*------------------------------------------------------------*/
{
    fint       *userval = (fint  *) arg;
    fint       nitems  = 1;
    fint       dfault  = HIDDEN;

    (void) userint_c( userval,
                      &nitems,
                      &dfault,
                      tofchar(key),
                      tofchar(" ") );
}



static void getampfilter( ident  id,
                          char  *key,
                          int    code,
                          void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: AMPRANGE= keyword handler                         */
/* These values serve as a filter for the output. Amplitudes  */
/* outside this range are filtered. Restrict the minimum      */
/* amplitude value to 0.                                      */
/*------------------------------------------------------------*/
{
   globaltype *g=(globaltype*)arg;
   double     userval[2];
   fint       nitems  = 2;
   fint       dfault  = HIDDEN;
   fint       r;

   r = userdble_c( userval,
                   &nitems,
                   &dfault,
                   tofchar(key),
                   tofchar(" ") );

   if (r == 0)
      g->filter.amp[0] = g->filter.amp[1] = dblank;
   if (r == 1)
   {
      g->filter.amp[0] = MYMAX( 0.0, userval[0] );
      g->filter.amp[1] = dblank;
   }
   if (r == 2)
   {
      g->filter.amp[0] = MYMAX( 0.0, userval[0] );
      g->filter.amp[1] = MYMAX( 0.0, userval[1] );
   }
}




static void getdispfilter( ident  id,
                           char  *key,
                           int    code,
                           void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: DISPRANGE= keyword handler                        */
/* These values serve as a filter for the output. Dispersions */
/* outside this range are filtered. Restrict the minimum      */
/* dispersion value to 0.                                     */
/*------------------------------------------------------------*/
{
   globaltype *g=(globaltype*)arg;
   double     userval[2];
   fint       nitems  = 2;
   fint       dfault  = HIDDEN;
   fint       r;

   r = userdble_c( userval,
                   &nitems,
                   &dfault,
                   tofchar(key),
                   tofchar(" ") );

   if (r == 0)
      g->filter.disp[0] = g->filter.disp[1] = dblank;
   if (r == 1)
   {
      g->filter.disp[0] = MYMAX( 0.0, userval[0] );
      g->filter.disp[1] = dblank;
   }
   if (r == 2)
   {
      g->filter.disp[0] = MYMAX( 0.0, userval[0] );
      g->filter.disp[1] = MYMAX( 0.0, userval[1] );
   }
}




static void getvelfilter( ident  id,
                          char  *key,
                          int    code,
                          void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: VELRANGE= keyword handler                         */
/* These values serve as a filter for the output. Velocities  */
/* (or the mean of the profile distribution) outside this     */
/* range are filtered.                                        */
/*------------------------------------------------------------*/
{
   globaltype *g=(globaltype*)arg;
   double     userval[2];
   fint       nitems  = 2;
   fint       dfault  = HIDDEN;
   fint       r;

   r = userdble_c( userval,
                   &nitems,
                   &dfault,
                   tofchar(key),
                   tofchar(" ") );

   if (r == 0)
      g->filter.vel[0] = g->filter.vel[1] = dblank;
   if (r == 1)
   {
      g->filter.vel[0] = userval[0];
      g->filter.vel[1] = dblank;
   }
   if (r == 2)
   {
      g->filter.vel[0] = userval[0];
      g->filter.vel[1] = userval[1];
   }
}




static void estimatetime( int      count,
                          int      maxcount,
                          double  *rt,
                          double  *ct )
/*------------------------------------------------------------*/
/* PURPOSE: Give an estimate how long the process needs to    */
/*          finish.                                           */
/*------------------------------------------------------------*/
{
   static double    realtime, cputime;
   fint             start = 0;
   fint             elapse = 1;


   if (count == 0)
   {
      timer_c( &cputime, &realtime, &start );
   }
   else
   {
      double remainingtime;
      int    minutes;
      double seconds;
      double average;
      char   mes[128];

      *rt = realtime;
      *ct = cputime;
      timer_c( ct, rt, &elapse );
      average = *rt / count;
      remainingtime = (double) (maxcount-(count+1)) * average;
      minutes = (int) remainingtime/60;
      seconds = remainingtime - minutes*60;
      sprintf( mes, "Estimated time to finish fitting: %d min %d sec",
               minutes, (int) seconds );
      gui_status( mes, FIRSTLINE );
      /* Update the progressbar */
      wkeyf("PROGRESS=%f", (double) (count+1) / (double) maxcount );
   }
}



static void resetlog( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Reset log file variables.                         */
/*------------------------------------------------------------*/
{
   int   i;

   g->log.noestimates = 0;
   g->log.nofit = 0;
   g->log.iterstoosmall = 0;
   g->log.ampfiltered = 0;
   g->log.velfiltered = 0;
   g->log.dispfiltered = 0;

   for (i = 0; i <= MAXLOGLINES; i++)
      gui_log( " ", i );
}



static void updatelog( globaltype *g,
                       int         flag,
                       int         count )
/*------------------------------------------------------------*/
/* PURPOSE: Update log file.                                  */
/*------------------------------------------------------------*/
{
   char           mes[128];
   int            line = 1;

   if (flag == INIT)
   {
      g->log.rchimax = 0.0;
      return;
   }
   sprintf( mes, "Processed profiles:                     %d (%d)",
            count, g->numprofiles );
   gui_log( mes, line++ );
   sprintf( mes, "Profiles without initial estimates:     %d",
            g->log.noestimates );
   gui_log( mes, line++ );
   sprintf( mes, "Profiles without least-squares fit:     %d",
            g->log.nofit );
   gui_log( mes, line++ );
   sprintf( mes, "No fit because # iterations > maxiters: %d",
            g->log.iterstoosmall );
   gui_log( mes, line++ );
   sprintf( mes, "Maximum red. chi^2:                     %6g",
            g->log.rchimax );
   gui_log( mes, line++ );
   sprintf( mes, "Rejected by amplitude filter:           %d",
            g->log.ampfiltered );
   gui_log( mes, line++ );
   sprintf( mes, "Rejected by velocity filter:            %d",
            g->log.velfiltered );
   gui_log( mes, line++ );
   sprintf( mes, "Rejected by dispersion filter:          %d",
            g->log.dispfiltered );
   gui_log( mes, line++ );
}




static int  setoutputsubsetlen( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Set the number of output subsets.                 */
/*------------------------------------------------------------*/
{
   int   result;
   int   estipars = 3;

   /*--------------------------------------------------*/
   /* Start with the number of estimates and standard  */
   /* parameters per component. There are always three */
   /* estimates per component. The number of standard  */
   /* parameters is entered by 'g->nvars'. Each funct- */
   /* ion can have its additional parameters per       */
   /* component.                                       */
   /*--------------------------------------------------*/

   if (g->estonly)
      return( 3 * g->ncomp );

   if (g->function == GAUSS)
   {
      /* Fit Amp, Cent, Disp. Extra: Area  */
      g->errpars = g->fitpars = 4;
   }
   else if (g->function == GAUSSHERMITE1)
   {
      /* Fit a, b, c, h3. */
      /* Extra: Area, Mean, Disp, Skewness, Real max, Pos of max.*/
      g->errpars = g->fitpars = 10;
   }
   else if (g->function == GAUSSHERMITE2)
   {
      /* Fit a, b, c, h3, h4. */
      /* Extra: Area, Mean, Disp, Skewness, Kurtosis, Real max, Pos of max.*/
      g->errpars = g->fitpars = 12;
   }
   else if (g->function == VOIGT)
   {
      /* Fit Amp, Cent, Disp. Lorentz. Extra: Area */
      g->errpars = g->fitpars = 5;
   }

   result = (estipars + g->fitpars + g->errpars) * g->ncomp;

   /*--------------------------------------------------*/
   /* Add three subsets for the backgound terms and    */
   /* three for the errors.                            */
   /*--------------------------------------------------*/
   result += 6;

   /*--------------------------------------------------*/
   /* Add subsets for the total flux, number of        */
   /* iterations, reduced chi-squared , standard       */
   /* deviation of the residuals and Q.                */
   /*--------------------------------------------------*/
   result += 5;

   return( result );
}



static void writeheaderitem( globaltype  *g,
                             char        *contents,
                             char        *units,
                             int          nr,
                             float       *minval,
                             float       *maxval,
                             fint        *nblanks )
/*------------------------------------------------------------*/
/* PURPOSE: Write subset id to set header and units to subset */
/*          header.                                           */
/*------------------------------------------------------------*/
{
   fint    r = 0;
   char    mes[STRLEN];


   sprintf( mes, 
           "%2d: %-35s %-10s %-10.3g %-10.3g  %-d", 
            nr, contents, units, minval[nr-1], maxval[nr-1], nblanks[nr-1] );
   anyoutf( 1, mes );  
   r = 0;
   gdsd_wvar_c( g->Setout,
                tofchar("COMMENT"),
                &setlevel,
                tofchar(mes),
                &r );
                 
   r = 0;
   gdsd_wchar_c( g->Setout,
                 tofchar("BUNIT"),
                 &(g->subout[nr-1]),            /* Index is subset number - 1 */
                 tofchar(units), &r );
}



static void writeheader( globaltype *g,
                         float      *minval,
                         float      *maxval,
                         fint       *nblanks )
/*------------------------------------------------------------*/
/* PURPOSE: Update output set header with subset info.        */
/*------------------------------------------------------------*/
{
   int      subnr;
   int      i;
   char     key[FITSLEN];
   char     contents[80];
   char     areaunits[80];
   char    *ptr;
   time_t   now;
   fint     r = 0;



   now = time( NULL );
   ptr = ctime( &now );

   /* sprintf( areaunits, "%s.%s", g->dataunits, g->axunits ); */
   strcpy( areaunits, "area" );
   anyoutf( 3, " " );
   anyoutf( 3, "==================  XGAUFIT =================" );
   anyoutf( 3, "Time & date           : %s", ptr );
   {
      char fname[STRLEN];
      if (g->function == GAUSS)         strcpy( fname, "Gauss function" );
      if (g->function == GAUSSHERMITE1) strcpy( fname, "Gauss-Hermite series with only h3!");
      if (g->function == GAUSSHERMITE2) strcpy( fname, "Gauss-Hermite series with h3 and h4" );
      if (g->function == VOIGT)         strcpy( fname, "Voigt line-shape" );
      anyoutf( 3, "Fitted function       : %s", fname );
   }
   anyoutf( 3, "Profile grid range    : %d %d", g->plohi[0], g->plohi[1] );
   anyoutf( 3, "Results set name      : %s", g->Setout.a );
   anyoutf( 3, "Profile direction     : %s", g->profaxname );
   anyoutf( 3, "Rms used for r.chi^2  : %g", g->rms );
   anyoutf( 3, " " );

   r = 0;
   gdsd_delete_c( g->Setout, tofchar("SUBTIME"), &setlevel, &r );
   sprintf( contents, "%s", ptr );
   r = 0;
   gdsd_wchar_c( g->Setout, tofchar("SUBTIME"), &setlevel, tofchar(contents), &r );

   /* If any of the SUBSETn keywords exist, delete them first */
   for (i = 0; i < g->noutsubs; i++)
   {
      sprintf( key, "SUBSET%d", i+1 );
      r = 0;
      gdsd_delete_c( g->Setout, tofchar(key), &setlevel, &r );
   }
   

   /* Write table header in comment field of set */
   {
      char   mes[STRLEN];
      fint   r = 0;
      int    l;
      
      l = sprintf( mes, "%3s %-35s %-10s %-10s %-10s %-10s",
                  "sub", "description", "units", "min.", "max.", "blanks" );
      anyoutf( 1, mes );
      gdsd_wvar_c( g->Setout,
                   tofchar("COMMENT"),
                   &setlevel,
                   tofchar(mes),
                   &r );
      memset( mes, '=', l ); mes[l] = '\0';
      anyoutf( 1, mes );
      r = 0;
      gdsd_wvar_c( g->Setout,
                   tofchar("COMMENT"),
                   &setlevel,
                   tofchar(mes),
                   &r );                           
   }
  

   /* Write names of subsets in header */
   subnr = 1;
   for (i = 1; i <= g->ncomp; i++)
   {
      sprintf( contents, "Estimate max amplitude" );
      writeheaderitem( g, contents, g->dataunits, subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Estimate center" );
      writeheaderitem( g, contents, g->axunits, subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Estimate dispersion" );
      writeheaderitem( g, contents, g->axunits, subnr++,
                       minval, maxval, nblanks );
   }
   if (!g->estonly)
   {
      for (i = 1; i <= g->ncomp; i++)
      {
         if (g->function == GAUSS || g->function == VOIGT)
         {
            sprintf( contents, "Fitted area" );
            writeheaderitem( g, contents, areaunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted area" );
            writeheaderitem( g, contents, areaunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Fitted center" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted center" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Fitted dispersion" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted dispersion" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            if (g->function == VOIGT)
            {
               sprintf( contents, "Fitted Lorentz factor", i );
               writeheaderitem( g, contents, g->axunits, subnr++,
                                minval, maxval, nblanks );

               sprintf( contents, "Error in fitted Lorentz factor");
               writeheaderitem( g, contents, g->axunits, subnr++,
                                minval, maxval, nblanks );
            }
            sprintf( contents, "Fitted max. amplitude" );
            writeheaderitem( g, contents, g->dataunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted max. amplitude", i );
            writeheaderitem( g, contents, g->dataunits, subnr++,
                             minval, maxval, nblanks );
         }
         if (g->function == GAUSSHERMITE1 || g->function == GAUSSHERMITE2)
         {
            sprintf( contents, "Fitted parameter a" );
            writeheaderitem( g, contents, g->dataunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted parameter a" );
            writeheaderitem( g, contents, g->dataunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Fitted parameter b" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted parameter b" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Fitted parameter c" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted parameter c" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Fitted parameter h3" );
            writeheaderitem( g, contents, "NUMBER", subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted parameter h3" );
            writeheaderitem( g, contents, "NUMBER", subnr++,
                             minval, maxval, nblanks );

            if (g->function == GAUSSHERMITE2)
            {
               sprintf( contents, "Fitted parameter h4" );
               writeheaderitem( g, contents, "NUMBER", subnr++,
                                minval, maxval, nblanks );

               sprintf( contents, "Error in fitted parameter h4" );
               writeheaderitem( g, contents, "NUMBER", subnr++,
                                minval, maxval, nblanks );
            }

            sprintf( contents, "Fitted area" );
            writeheaderitem( g, contents, areaunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted area" );
            writeheaderitem( g, contents, areaunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Fitted mean" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted mean" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Fitted dispersion" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted dispersion" );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Fitted skewness" );
            writeheaderitem( g, contents, "NUMBER", subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted skewness" );
            writeheaderitem( g, contents, "NUMBER", subnr++,
                             minval, maxval, nblanks );

            if (g->function == GAUSSHERMITE2)
            {
               sprintf( contents, "Fitted kurtosis" );
               writeheaderitem( g, contents, "NUMBER", subnr++,
                                minval, maxval, nblanks );

               sprintf( contents, "Error in fitted kurtosis" );
               writeheaderitem( g, contents, "NUMBER", subnr++,
                                minval, maxval, nblanks );
            }

            sprintf( contents, "Fitted real maximum" );
            writeheaderitem( g, contents, g->dataunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted real maximum" );
            writeheaderitem( g, contents, g->dataunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Fitted position of real max." );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );

            sprintf( contents, "Error in fitted position of real max." );
            writeheaderitem( g, contents, g->axunits, subnr++,
                             minval, maxval, nblanks );
         }
      }

      /* Common to all functions */

      sprintf( contents, "Background constant term Z0");
      writeheaderitem( g, contents, g->dataunits, subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Error in Z0");
      writeheaderitem( g, contents, g->dataunits, subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Background linear term Z1");
      writeheaderitem( g, contents, "NUMBER", subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Error in  Z1");
      writeheaderitem( g, contents, "NUMBER", subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Background quadratic term Z2");
      writeheaderitem( g, contents, "NUMBER", subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Error in Z2");
      writeheaderitem( g, contents, "NUMBER", subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Total flux");
      writeheaderitem( g, contents, areaunits, subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Number of lsqfit iterations");
      writeheaderitem( g, contents, "NUMBER", subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Reduced chi-squared using an rms=%g", g->rms);
      writeheaderitem( g, contents, "NUMBER", subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Standard deviation of residuals");
      writeheaderitem( g, contents, "NUMBER", subnr++,
                       minval, maxval, nblanks );

      sprintf( contents, "Smoothing factor Q");
      writeheaderitem( g, contents, "NUMBER", subnr++,
                       minval, maxval, nblanks );
   }
   anyoutf( 3, " " );
}



static void createoutputsetondisk( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Set the number of output subsets and create set.  */
/*------------------------------------------------------------*/
{
   fint     dfault = HIDDEN+NOERRORLOOP;   /* Return immediately when error occurs */
   fint     class = 2;
   fint     pmask;
   double   outcdelt = 1.0;
   double   outcrota = 0.0;
   double   outcrpix = 0.0;
   double   outcrval = 0.0;
   fint     showdev = DEBUG;
   fint     axcountout;
   fint     maxaxes = MAXAXES;          /* Max num. of axes the program can deal with.*/
   fint     nsubsout;
   fint     ax;
   fint     r;


   /*--------------------------------------------------*/
   /* Always start with a new output set! Then the user*/
   /* can change box and grid range and keep the output*/
   /* set name the same, without troubles with resizing*/
   /* the output set that already existed.             */
   /*--------------------------------------------------*/
   r = 0;
   if (gds_exist_c(tofchar(g->setoutbuf), &r))   
   {
      r = 0;
      gds_delete_c( tofchar(g->setoutbuf), &r ); 
   }
               

   /*--------------------------------------------------*/
   /* There is a difference between the numbering of   */
   /* the axes for the input and the output. For the   */
   /* output we already had the axes permutations and  */
   /* the axis that we want to change is now always    */
   /* the last axis.                                   */
   /*--------------------------------------------------*/
   ax = g->subdim+1;

   /* Create an output set */
   gdsasn_c( tofchar("INSET="), KEY_OUTSET, &class );
   gdscss_c( KEY_OUTSET, g->blo, g->bhi );
   /*--------------------------------------------------*/
   /* Create subsets for estimates and fits for all    */
   /* components. The number of estimates is 3 per     */
   /* component. The number of fit parameters depends  */
   /* in the selected function. Extra parameters are   */
   /* those that are the same for all components       */
   /* e.g. chi2 etc.                                   */
   /*--------------------------------------------------*/


   g->noutsubs = setoutputsubsetlen( g );

   /* pmask=14=(2 + 4 + 8) means CTYPE,CRPIX & CRVAL are defined.*/
   pmask = 14;
   gdscpa_c(
         KEY_OUTSET,              /* Keyword associated with a GDSOUT call */
         &ax,                     /* The (output) number of the axis to be changed.*/
         &(g->noutsubs),          /* Size of the axis.*/
         &outcdelt,               /* Increment in physical units along axis.*/
         &outcrota,               /* Rotation angle of axis.*/
         &outcrpix,               /* Reference pixel of axis.*/
         &outcrval,               /* Physical reference value at reference pixel.*/
         tofchar(NEWAXISNAME),    /* New axis name.*/
         tofchar(" "),            /* Physical units of axis.*/
         &pmask );                /* Bit mask denoting which of the six  values above are defined.*/


   dfault = HIDDEN + NOERRORLOOP;
   strcat( g->Setout.a, " " );
   strcat( g->Setout.a, NEWAXISNAME );
   nsubsout = gdsout_c(
                     g->Setout,
                     g->subout,
                     &(g->noutsubs),
                     &dfault,
                     KEY_OUTSET,
                     tofchar(""),
                     &showdev,
                     g->axnumout,
                     &axcountout,
                     &maxaxes );

   g->Setout.a[nelc_c(g->Setout)] = '\0';
}




static void creategausscube( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Copy input set to create an output set for the    */
/* fitted data (not the parameters). A class 2 input          */
/* set is entered, therefore we cannot use gdsasn etc to copy */
/* the input because then we lose the last axis. Program COPY */
/* can do the task better. We copy the n-dim input set with   */
/* the n-dim new box and reread this copied set with gdsinp   */
/* without subsets, e.g. as a set.                            */
/*------------------------------------------------------------*/
{
   fint     r;
   

   /*--------------------------------------------------*/
   /* Always start with a new output set! Then the user*/
   /* can change box and grid range and keep the output*/
   /* set name the same, without troubles with resizing*/
   /* the output set that already existed.             */
   /*--------------------------------------------------*/
   r = 0;
   if (gds_exist_c(tofchar(g->setgaubuf), &r))   
   {
      r = 0;
      gds_delete_c( tofchar(g->setgaubuf), &r ); 
   }
   /* Prepare the call to deputy */               
   {
      char   mes[512];
      char   grid[32];      
      int    i;
      fint   r;
      fint   status;
      fint   hermesold, hermesnew; 
      
      strcpy( mes, "" );
      for (i = 0; i < g->subdim; i++)
      {
         sprintf( grid, "%d ", g->blo[i] );
         strcat( mes, grid );
      }
      sprintf( grid, "%d ", g->plohi[0] );
      strcat( mes, grid );  
      for (i = 0; i < g->subdim; i++)
      {
         char   grid[32];
         sprintf( grid, "%d ", g->bhi[i] );
         strcat( mes, grid );
      }
      sprintf( grid, "%d ", g->plohi[1] );      
      strcat( mes, grid );  
      anyoutf( DEBUG, mes );
      r = 0; subst_c( tofchar("BOX=R_BOX="), &r );
      wkeyf( "R_BOX=%s", mes );
      r = 0; subst_c( tofchar("INSET=R_INSET="), &r );  
      wkeyf( "R_INSET=%s", g->Setin.a );
      r = 0; subst_c( tofchar("OUTSET=R_OUTSET="), &r );
      wkeyf( "R_OUTSET=%s", g->Setgau.a );   
      wkeyf( "MAKEBLANK=Y" );
      hermesnew = 0;        /* No output for deputy in log file */
      hermesold = listctrl_c( &hermesnew ); 
      deputy_c( tofchar("COPY"), &status ); 
      (void) listctrl_c( &hermesold ); 
   }      
}




static int allocatebuffers( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Allocate memory for data reading and writing.     */
/*------------------------------------------------------------*/
{
   char    mes[STRLEN];
   /*--------------------------------------------------*/
   /* First the profiles for estimates and fit.        */
   /*--------------------------------------------------*/

   sprintf( mes, "Allocated %d Kbytes for input data",
            (sizeof(double) * MAXBUF * g->profilelen)/1024 );
   gui_log( mes, 10 );
   g->profiles = dmatrix( 0, 0, g->profilelen-1, MAXBUF-1 );
   if (g->profiles == NULL)
   {
      return( NO );
   }
   /*--------------------------------------------------*/
   /* Index notation f/dmatrix is M[y][x] if fmatrix   */
   /* called as fmatrix( xlo, ylo, xhi, yhi ). So      */
   /* change this order to get notation M[x][y].       */
   /*--------------------------------------------------*/


   /*--------------------------------------------------*/
   /* Try to allocate memory for the output.           */
   /*--------------------------------------------------*/
   sprintf( mes, "Allocated %d Kbytes for output",
            (sizeof(double) * MAXBUF * g->noutsubs)/1024 );
   gui_log( mes, 11 );

   g->results = dmatrix( 0, 0,  MAXBUF-1, g->noutsubs-1 );
   if (g->results == NULL)
   {
      return( NO );
   }


   /*--------------------------------------------------*/
   /* Try to allocate memory for a buffer for          */
   /* alternative profiles for the estimates routine.  */
   /*--------------------------------------------------*/
   if (g->altset)
   {
      int   numestsubsets = 3 * g->ncomp;
      sprintf( mes, "Allocated %d Kbytes for alternative set",
               (sizeof(double) * MAXBUF * numestsubsets)/1024 );
      gui_log( mes, 12 );

      g->estfromset = dmatrix( 0, 0, numestsubsets-1, MAXBUF-1 );
      if (g->estfromset == NULL)
      {
         return( NO );
      }
   }

   return( YES );                /* All ok... */
}




static void setsort( ident  id,
                     char  *key,
                     int    code,
                     void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: Set sorting on or off.                            */
/*          Number of components must be > 1                  */
/*------------------------------------------------------------*/
{
   globaltype *g = (globaltype*) arg;
   fint       userval = 1;
   fint       nitems  = 1;
   fint       dfault  = HIDDEN;


   (void) userint_c( &userval,
                     &nitems,
                     &dfault,
                     tofchar(key),
                     tofchar(" ") );

   if (userval > 0 && userval <= 3 && g->ncomp > 1)
      g->sortopt = userval;
   else
      g->sortopt= 0;
}




static void processdata( ident  id,
                         char  *key,
                         int    code,
                         void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: GO= keyword handler. Start loop over all          */
/*          profiles in current box.                          */
/*------------------------------------------------------------*/
{
   globaltype *g=(globaltype*)arg;
   bool        value;
   float       perc = 0.0;
   int         aborted = NO;
   char        mes[128];
   int         i;
   fint        tidI[MAXSUBSETS];
   fint        tidO[MAXSUBOUT];
   fint        tidE[3*MAXGAUSS];    /* Transfer id's for estimate set */
   int         rest;
   double      rt, ct;
   float       *minval = NULL, *maxval = NULL;
   fint        *nblanks = NULL, *mcount = NULL;


   value = UserLog( key );
   if (!value)
      return;

   clearstatus();
   if (!g->boxexist)
   {
      error_status( "No box defined!", "Enter a valid box first!", "" );
      wkeyf( key );
      return;
   }
   if (!g->outsetexist)
   {
      error_status( "No output set entered!",
                    "Enter a valid output set name first!", "" );
      wkeyf( key );
      return;
   }

   /*--------------------------------------------------*/
   /* If there is an alternative set for the estimates */
   /* then check whether the number of input estimate  */
   /* subsets equals the number of selected components */
   /* times 3 estimate parameters.                     */
   /*--------------------------------------------------*/
   if (g->altset)
   {
      if (g->naltsubs != 3*g->ncomp)
      {
         char     mes1[STRLEN];
         char     mes2[STRLEN];
         sprintf( mes1, "Exactly %d estimate subsets required!", 3*g->ncomp );
         sprintf( mes2, "You entered however %d estimate subsets", g->naltsubs );
         error_status( mes1, mes2, "Check estimate input or change number of components" );
         wkeyf( key );
         return;
      }
   }

   g->estonly = UserLog( "ESTONLY=" );           /* Calculate estimates only? */
   setsort( NULL, "SORT=", KEYCHANGE, g );              /* Is sorting wanted? */
   createoutputsetondisk( g );
   
   /* Create a template set for the output of the fitted cube */

   if (g->gausetexist)                           
      creategausscube( g );
      
   resetlog( g );
   /*--------------------------------------------------*/
   /* Get initial guesses for fixed parameters.        */
   /*--------------------------------------------------*/
   setuserguess( g );


   /*--------------------------------------------------*/
   /* Allocate memory for the IO buffers.              */
   /*--------------------------------------------------*/
   if (!allocatebuffers( g ))
   {
      error_status( "Cannot allocate memory for IO buffers!",
                    "Release some memory and start again",
                    "or use a machine with more memory." );
      wkeyf( key );
      return;
   }
   minval  = (float *) calloc( g->noutsubs, sizeof(float) );
   maxval  = (float *) calloc( g->noutsubs, sizeof(float) );
   nblanks = (fint *)  calloc( g->noutsubs, sizeof(fint) );
   mcount  = (fint *)  calloc( g->noutsubs, sizeof(fint) );
   if (minval == NULL || maxval == NULL || nblanks == NULL || mcount == NULL)
   {
      error_status( "Cannot allocate memory for minmax vars.",
                    "Release some memory and start again",
                    "or use a machine with more memory." );
      wkeyf( key );
      return;                                     
   }


   /*--------------------------------------------------*/
   /* Next part provides a loop over the data using a  */
   /* buffer for MAXBUF profiles of arbitrary length.  */
   /*--------------------------------------------------*/
   g->count = 0;
   rest = g->numprofiles;
   /* Initialize the transfer id's */
   for (i = 0; i < g->profilelen; i++)
   {
      tidI[i] = 0;
   }
   for (i = 0; i < g->noutsubs; i++)
   {
      tidO[i] = 0;
      mcount[i] = 0;
   }
   if (g->altset)
   {
      for (i = 0; i < g->ncomp * 3; i++)
      {
         tidE[i] = 0;
      }
   }

   /* Start reading, processing and writing */

   updatelog( g, INIT, 0 );
   do
   {
      int     s;
      float   buf[MAXBUF];
      fint    bufsize = MYMIN( MAXBUF, rest );
      fint    bufnumprofiles;
      fint    pixelsread;



      /*--------------------------------------------------*/
      /* First read a buffer with profiles to fit. The    */
      /* data can also be used to calculate estimates.    */
      /* The estimates can also be imported from a set.   */
      /*--------------------------------------------------*/
      for (s = 0; s < g->profilelen; s++)
      {
         fint    cwlo, cwhi;

         cwlo   = gdsc_fill_c( g->Setin, &(g->subin[s]), g->blo );
         cwhi   = gdsc_fill_c( g->Setin, &(g->subin[s]), g->bhi );
         gdsi_read_c( g->Setin,
                      &cwlo, &cwhi,
                      buf,
                      &bufsize,
                      &pixelsread,
                      &tidI[s] );
         for (i = 0; i < pixelsread; i++)
         {
            /* Transfer the subset data to the profiles buffers */
            float y = buf[i];
            if (y != fblank)
               g->profiles[i][s] = (double) y;
            else
               g->profiles[i][s] = dblank;
         }
      }
      /* Now we have MAXBUF (or less) profiles and we can process them */

      bufnumprofiles = pixelsread;

      /*--------------------------------------------------*/
      /* If the user wants estimates from a set, then     */
      /* read now them and convert to doubles.            */
      /*--------------------------------------------------*/
      if (g->altset)
      {
         for (s = 0; s < g->ncomp * 3; s++)
         {
            fint    cwlo, cwhi;

            cwlo   = gdsc_fill_c( g->Setest, &(g->subest[s]), g->blo );
            cwhi   = gdsc_fill_c( g->Setest, &(g->subest[s]), g->bhi );
            gdsi_read_c( g->Setest,
                         &cwlo, &cwhi,
                         buf,
                         &bufnumprofiles,
                         &pixelsread,
                         &tidE[s] );

            for (i = 0; i < bufnumprofiles; i++)
            {
               /* Transfer the subset estimates to the estimate buffer */
               float y = buf[i];
               if (y != fblank)
                  g->estfromset[i][s] = (double) y;
               else
                  g->estfromset[i][s] = dblank;
            }
         }
      }

      /*--------------------------------------------------*/
      /* Do the fitting for all the profiles in the       */
      /* input buffer. Update a log every 'blocksize'     */
      /* profiles.                                        */
      /*--------------------------------------------------*/
      g->currentprofilenumber = 0;      /* Reset this counter after reading a buffer */
      for (i = 0; i < bufnumprofiles; i++, g->currentprofilenumber++, g->count++)
      {
         int   maxcount = g->numprofiles;
         int   blocksize = 200;

         if (!(g->count%blocksize))
         {
            estimatetime( g->count, maxcount, &rt, &ct );
            if ( !UserLog(key) )
            {
               aborted = YES;
               perc = 100.0 * (float) i / (float) maxcount;
               break;
            }
            updatelog( g, UPDATE, g->count );
            GgiHandleEvents();
         }
         dofit( g );
      }

      /*--------------------------------------------------*/
      /* Convert the results to floats and store them in  */
      /* the output set.                                  */
      /*--------------------------------------------------*/
      for (s = 0; s < g->noutsubs; s++)
      {
         fint    cwlo, cwhi;
         fint    pixelsdone;

         cwlo   = gdsc_fill_c( g->Setout, &(g->subout[s]), g->blo );
         cwhi   = gdsc_fill_c( g->Setout, &(g->subout[s]), g->bhi );
         for (i = 0; i < bufnumprofiles; i++)
         {
            buf[i] = tofloat( g->results[s][i] );
         }
         gdsi_write_c( g->Setout,
                       &cwlo, &cwhi,
                       buf,
                       &bufnumprofiles,
                       &pixelsdone,
                       &tidO[s] );
                       
         minmax3_c( buf, 
                    &bufnumprofiles, 
                    &minval[s],
                    &maxval[s],
                    &nblanks[s],
                    &mcount[s] );                    
      }
      rest -= pixelsread;
   }
   while (tidI[0] != 0);         /* Repeat actions if there are profiles left */


   /*----------------------------------------*/
   /* Update header of new set.              */         
   /*----------------------------------------*/
   {
      fint  rmove = YES;  
      fint  nsub = g->noutsubs;
      wminmax_c( g->Setout,
                 g->subout,
                 minval, maxval,
                 nblanks,
                 &nsub,
                 &rmove );
   }

   /*----------------------------------------*/
   /* Update the gui log etc.                */
   /*----------------------------------------*/
   wkeyf("PROGRESS=1.0");
   updatelog( g, UPDATE, g->numprofiles );
   sprintf( mes, "Processed %d profiles in %d sec (%d cpu sec)",
            g->numprofiles, (int) rt, (int) ct );
   gui_log( mes, MAXLOGLINES );
   GgiHandleEvents();

   clearstatus();
   if (aborted)
   {
      char   mes[128];
      sprintf( mes, "Aborted after processing %.1f percent of the data", perc );
      gui_status( mes, FIRSTLINE );
   }
   else
   {
      char   mes[128];
      fint   r = 0;
      
      sprintf( mes, "Ready writing parameters to set: %s", g->Setout.a );
      gui_status( mes, FIRSTLINE );
      if (g->gausetexist) 
      {
         char    *ptr;
         time_t   now;
              
         now = time( NULL );
         ptr = ctime( &now );
         sprintf( mes, "Ready writing data to set: %s", g->Setgau.a );  
         gui_status( mes, SECONDLINE );
         sprintf( mes, "Model set created by XGAUFIT at: %s", ptr ); 
         gdsd_wvar_c( g->Setgau, 
                      tofchar("HISTORY"),
                      &setlevel,
                      tofchar(mes),  
                      &r ); 
      }
      wkeyf( key );
   }

   /* Write description, min, max, blanks and units to header */
   writeheader( g, minval, maxval, nblanks ); 


   free( minval );
   free( maxval );
   free( nblanks );
   free( mcount );
  
   if (g->altset)
   {
      freedmatrix( g->estfromset, 0, 0 );
   }
   freedmatrix( g->results, 0, 0 );
   freedmatrix( g->profiles, 0, 0 );

   /* Close output set */
   {
      fint r = 0;
      gds_close_c( g->Setout, &r );
   }
   wkeyf("PROGRESS=0.0" );
}




static void availableaxes( fchar Set )
/*------------------------------------------------------------*/
/* PURPOSE: Fill gui menu with axis names of current set      */
/*------------------------------------------------------------*/
{
   static char  axnames[MAXAXES][FITSLEN];
   static char  *ax[MAXAXES];
   int          i;
   fint         r;
   fint         dim;


  /* ax = axnames;*/
   dim = gdsc_ndims_c( Set, &setlevel );
   for (i = 0; i < dim; i++)
   {
      fchar  Name;
      fint   axnum = (fint) i+1;                 /* Axis numbers start with 1 */
      Name.a = axnames[i];
      Name.l = FITSLEN;
      r = 0;
      gdsc_name_c( Name, Set, &axnum, &r );
      ax[i] = Name.a;
      ax[i][nelc_c(Name)] = '\0';
      strtok( ax[i], " -" );
   }
   ax[i] = (char *) NULL;     /* Terminate array */
   gui_axnames( ax );
}



static void getbox( ident id,
                    char  *key,
                    int   code,
                    void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: BOX= Process user input for a box.                */
/*------------------------------------------------------------*/
{
   fint        boxopt = 0;
   fint        dfault = HIDDEN;
   fint        showdev = DEBUG;
   fint        r;
   char        dummytxt[STRLEN+1];
   fchar       Dummytxt;
   globaltype *g=(globaltype*)arg;


   if (!g->setexist)
   {
      rejectf( "BOX=", "You cannot define a box without a set!" );
      g->imagesize = 0;
      g->boxexist = NO;
      return;                     /* Nothing changed so return old image data */
   }

   if (g->inactive)
   {
      g->boxexist = NO;
      g->inactive = NO;
      cancel("CURPOS=");
      return;
   }

   /* Check on empty strings */
   Dummytxt.l = STRLEN;
   Dummytxt.a = dummytxt;
   r = usertext_c( Dummytxt, &dfault, tofchar(key), tofchar(" ") );
   if (!r)
   {
      char   mes1[128];
      char   mes2[128];
      int    i;

      /* Prepare the default for 'gdsbox'. With option 16 */
      /* it does not accept an empty string as input.     */

      strcpy( mes1, "BOX=" );
      for (i = 0; i < g->subdim; i++)
      {
         sprintf( mes2, " %d", g->flo[i] );
         strcat( mes1, mes2 );
      }
      for (i = 0; i < g->subdim; i++)
      {
         sprintf( mes2, " %d", g->fhi[i] );
         strcat( mes1, mes2 );
      }
      wkeyf( mes1 );
      return;
   }

   dfault = REQUEST;
   boxopt = 16;
   gdsbox_c( g->blo, g->bhi,
             g->Setin, g->subin,
             &dfault,
             tofchar(key),
             tofchar(" "),
             &showdev,
             &boxopt );

   if (boxopt == -1)
   {
      /* 'gdsbox' detected an error */
      fchar   Errtxt;
      char    errtxt[STRLEN];
      fint    l;
      fint    dfault = -1;

      Errtxt.a = errtxt; Errtxt.l = STRLEN-1; clearstr(Errtxt);
      cancel( "UCAMESSAGE=" );
      l = usertext_c( Errtxt, &dfault, tofchar("UCAMESSAGE="), tofchar("") );
      Errtxt.a[l] = '\0';

      error_status( Errtxt.a, "Enter another box.", "" );

      g->imagesize = 0;
      g->boxexist = NO;
      return;
   }

   /*--------------------------------------------------*/
   /* With the BOX= keyword, we set the number of      */
   /* profiles to be examined.                         */
   /*--------------------------------------------------*/
   {
      int  i;

      g->numprofiles = 1;
      for (i = 0; i < g->subdim; i++)
      {
         g->numprofiles *= g->bhi[i] - g->blo[i] + 1;
      }
   }
   g->boxexist = YES;
   clearstatus();
   return;
}





static void getsetname( ident  id,
                        char  *key,
                        int    code,
                        void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: SETNAME= keyword handler.                         */
/*------------------------------------------------------------*/
{
   fint          r;
   fint          dfault;
   fint          nitems;
   char          dummytxt[STRLEN+1];
   fchar         Dummytxt;
   globaltype   *g=(globaltype*)arg;


   clearstatus();

   g->setnameexist = NO;
   dfault     = HIDDEN;
   Dummytxt.l = STRLEN;
   Dummytxt.a = dummytxt;
   nitems = 1;
   r = userchar_c( Dummytxt, &nitems, &dfault, tofchar("SETNAME="), tofchar(" ") );
   Dummytxt.a[nelc_c(Dummytxt)] = '\0';
   if (r == 0)
      return;
   else
   {
      r = 0;
      if (gds_exist_c(Dummytxt, &r))
      {
         strcpy( g->setname, Dummytxt.a );
         g->setnameexist = YES;
         /* Clear box if a NEW set (i.e. there is a box) is entered */
         if (g->boxexist)
         {
            g->inactive = YES;
            wkeyf( "BOX=" );
         }
         availableaxes( Dummytxt );
         /*--------------------------------------------------*/
         /* If a new set is entered while a profile axis     */
         /* exist then trigger the profile function. Perhaps */
         /* the user wanted to examine the same profile axis.*/
         /* Clear the plot contents.                         */
         /*--------------------------------------------------*/
         if (g->profileexist)
         {
            TriggerKey("PROFILE=");
         }
      }
      else
      {
         rejectf( "SETNAME=", "Cannot find set!" );
         error_status( "Cannot find set!",
                       "Check spelling and directory", "" );
         g->setnameexist = NO;
         return;
      }
   }
}




static void getprofile( ident  id,
                        char  *key,
                        int    code,
                        void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: PROFILE= keyword handler                          */
/*                                                            */
/* If a profile axis name is also entered then the            */
/* set will be read by 'gdsinp'. If then a range is not yet   */
/* specified, a default is generated and the GRIDRANGE=       */
/* keyword is triggered again.                                */
/*------------------------------------------------------------*/
{
   fint          axcount[MAXAXES];
   fint          maxsubs = MAXSUBSETS;
   fint          maxaxes = MAXAXES;          /* Max num. of axes the program can deal with.*/
   fint          class   = 2;                /* Class 2 is for profile applications */
   fint          showdev = DEBUG;
   fint          r, r1, r2;
   fint          dfault;
   char          dummytxt[STRLEN+1];
   fchar         Dummytxt;
   fint          setdim;
   globaltype   *g=(globaltype*)arg;
   int           i;
   fint          grid = 0;


   clearstatus();

   g->profileexist = NO;
   if (!g->setnameexist)
   {
      error_status( "There is no set name",
                    "Enter name of set first", "" );
      return;
   }

   dfault     = HIDDEN;
   Dummytxt.l = STRLEN;
   Dummytxt.a = dummytxt;
   r = usertext_c( Dummytxt, &dfault, tofchar("PROFILE="), tofchar(" ") );
   Dummytxt.a[r] = '\0';
   if (r == 0)
      return;
   else
   {
      strcpy( g->profaxname, Dummytxt.a );
   }
   wkeyf( "DUMMYSET=%s %s", g->setname, g->profaxname );

   dfault = HIDDEN+NOERRORLOOP;          /* Return immediately when error occurs */
   g->axesoutsidesub = 1;                /* Allow ONE axis, i.e. a profile */
   clearstr( g->Setin );

   g->profilelen = gdsinp_c( g->Setin,
                        g->subin,
                        &maxsubs,
                        &dfault,
                        tofchar( "DUMMYSET=" ),
                        tofchar(" "),
                        &showdev,        /* Device number (as in ANYOUT). */
                        g->axnum,
                        axcount,         /* Number of grids on axes in 'axnum' */
                        &maxaxes,
                        &class,
                        &(g->axesoutsidesub) );

   if (g->profilelen <= 0)
   {
      fchar   Errtxt;
      char    errtxt[STRLEN];
      fint    dfault = -1;              /* Special default to get UCA message */
      fint    l;

      Errtxt.a = errtxt; Errtxt.l = STRLEN-1; clearstr(Errtxt);
      g->setexist = NO;
      g->imagesize = 0;
      l = usertext_c( Errtxt, &dfault, tofchar("UCAMESSAGE="), tofchar("") );
      Errtxt.a[l] = '\0';
      error_status( Errtxt.a,
                    "Check name of profile axis, range",
                    "and minimum axis length of profile (length > 1)." );
      g->profileexist = NO;
      return;
   }

   setdim  = gdsc_ndims_c( g->Setin, &setlevel );
   g->subdim = setdim - g->axesoutsidesub;
   g->profaxnum = g->axnum[g->subdim];
   g->profileexist = YES;
   g->profcw = 0;
   for (i = 0; i < g->subdim; i++)
   {
      r = 0;
      g->profcw = gdsc_word_c( g->Setin, &(g->axnum[i]), &grid, &(g->profcw), &r );
   }

   /*--------------------------------------------------*/
   /* Now we have all the set properties needed for    */
   /* conversions from grids to physical coordinates.  */
   /* The actual reading of the set with user given or */
   /* default range on the profiles axis is done later.*/
   /* If user did not enter a grid or physical range   */
   /* on the command line then take the default range  */
   /* in grids and update the physical range.          */
   /*--------------------------------------------------*/


   /* Read minimum and maximum grid values on profile axis */
   {
      fint  cwlo, cwhi;
      gdsc_range_c( g->Setin, &setlevel, &cwlo, &cwhi, &r1 );
      r1 = 0;
      g->profgridlo = gdsc_grid_c( g->Setin, &(g->profaxnum), &cwlo, &r1 );
      r2 = 0;
      g->profgridhi = gdsc_grid_c( g->Setin, &(g->profaxnum), &cwhi, &r2 );
   }


   /* Is a range specified? */

   r1 = usertext_c( Dummytxt, &dfault, tofchar("GRIDRANGE="), tofchar(" ") );
   r2 = usertext_c( Dummytxt, &dfault, tofchar("PHYSRANGE="), tofchar(" ") );

   /* Grid has higher priority than physical coordinates when */
   /* entered on the command line. */

   if (!g->boxexist)
   {
      if (r1)
         TriggerKey("GRIDRANGE=");
      else if (r2)
         TriggerKey("PHYSRANGE=");
   }

   /* Clear box field if a NEW profile is entered */
   if (g->boxexist)
   {
      g->inactive = YES;
      wkeyf( "BOX=" );
   }


   if (g->boxexist || (!r1 && !r2) )
   {
      /* We must set the default */
      g->plohi[0] = g->profgridlo;
      g->plohi[1] = g->profgridhi;
      wkeyf( "GRIDRANGE=%d %d", g->plohi[0], g->plohi[1] );
   }
}




static void getset( ident  id,
                    char  *key,
                    int    code,
                    void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: INSET= keyword handler                            */
/*------------------------------------------------------------*/
{
   fint          axcount[MAXAXES];
   fint          maxsubs = MAXSUBSETS;
   fint          maxaxes = MAXAXES;          /* Max num. of axes the program can deal with.*/
   fint          class   = 2;                /* Class 2 is for profile applications */
   fint          showdev = DEBUG;
   fint          r;
   fint          dfault;
   fint          setdim;
   globaltype   *g=(globaltype*)arg;


   clearstatus();

   if (!g->setnameexist)
   {
      error_status( "There is no set name",
                    "Enter name of set first", "" );
      return;
   }
   if (!g->profileexist)
   {
      error_status( "There is no profile axis name entered",
                    "Enter name of profile axis", "" );
      return;
   }


   r = 0;
   gds_close_c( g->Setin, &r );

   wkeyf( "INSET=%s %s %d:%d", g->setname, g->profaxname, g->plohi[0], g->plohi[1] );
   anyoutf( DEBUG, "Final opening of set! [%s]", g->setname );

   dfault = HIDDEN+NOERRORLOOP;          /* Return immediately when error occurs */
   g->axesoutsidesub = 1;                /* Allow ONE axis, i.e. a profile */
   clearstr( g->Setin );

   g->profilelen = gdsinp_c(
                        g->Setin,
                        g->subin,
                        &maxsubs,
                        &dfault,
                        tofchar( "INSET=" ),
                        tofchar(" "),
                        &showdev,        /* Device number (as in ANYOUT). */
                        g->axnum,
                        axcount,         /* Number of grids on axes in 'axnum' */
                        &maxaxes,
                        &class,
                        &(g->axesoutsidesub) );

   if (g->profilelen <= 0)
   {
      fchar   Errtxt;
      char    errtxt[STRLEN];
      fint    dfault = -1;              /* Special default to get UCA message */
      fint    l;

      Errtxt.a = errtxt; Errtxt.l = STRLEN-1; clearstr(Errtxt);
      g->setexist = NO;
      g->imagesize = 0;
      l = usertext_c( Errtxt, &dfault, tofchar("UCAMESSAGE="), tofchar("") );
      Errtxt.a[l] = '\0';
      error_status( Errtxt.a,
                    "Check name of profile axis, range",
                    "and minimum axis length of profile (length > 1)." );
      g->profileexist = NO;
      return;
   }
   g->setexist = YES;
   setdim  = gdsc_ndims_c( g->Setin, &setlevel );

   /*-------------------------------*/
   /* Determine edges of input      */
   /*-------------------------------*/
   {
      fint cwlo, cwhi;                          /* Local coordinate words */
      int  m;
      fint r2, r1 = 0;
      gdsc_range_c( g->Setin, &setlevel, &cwlo, &cwhi, &r1 );
      r1 = r2 = 0;
      for (m = 0; m < setdim; m++)
      {
         g->flo[m] = gdsc_grid_c( g->Setin, &(g->axnum[m]), &cwlo, &r1 );
         g->fhi[m] = gdsc_grid_c( g->Setin, &(g->axnum[m]), &cwhi, &r2 );
      }
   }

   /*-------------------------------*/
   /* Allocate memory for profile   */
   /* arrays.                       */
   /*-------------------------------*/
   if ( !allocs(g, g->profilelen) )
   {
      error_status( "Cannot allocate enough memory!",
                    "No memory left for profile work arrays.",
                    "Perhaps you should restart with a smaller box." );
      g->setexist = NO;
      g->profileexist = g->boxexist = NO;
      return;
   }

   /*-------------------------------*/
   /* Get the grids on the profiles */
   /*-------------------------------*/
   {
      int i;
      for (i = 0; i < g->profilelen; i++)
         g->xgrids[i] = (double) (i + g->plohi[0]);
   }
   dminmax( g->xgrids, &(g->profilelen), &(g->xgridminmax[0]), &(g->xgridminmax[1]) );


   /*-------------------------------*/
   /* Transform also to physical    */
   /* coordinates.                  */
   /*-------------------------------*/
   {
      int i;

      g->tophys = YES;
      for (i = 0; i < g->profilelen; i++)
      {
         int r;

         r = grid2phys2( g->Setin, g->xgrids[i], g->profcw, g->profaxnum, &(g->xphys[i]) );
         /* If conversion not possible, stop loop immediately */
         if (r == 0)
         {
            g->tophys = NO;
            break;
         }
      }
   }

   /*-------------------------------*/
   /* Get the grid spacing.         */
   /*-------------------------------*/
   g->gridspac = (g->xphys[g->profilelen-1]-g->xphys[0] ) /
                                                (double) (g->profilelen - 1);


   /*-------------------------------*/
   /* Transform physical coordinates*/
   /* to interval [-1,1].           */
   /*-------------------------------*/
   if (g->tophys)
   {
      scalearray( g->xphys, g->xphysscale, g->profilelen );
   }
   if (!TriggerKey("VSYS="))
      wkeyf( "VSYS=%g", g->xphys[g->profilelen/2] );

   /*-------------------------------*/
   /* Get a string with axis names  */
   /* for annotation in gui and plot*/
   /*-------------------------------*/
   {
      int    i;

      strcpy( g->axisnames, "" );
      for (i = 0; i < g->subdim; i++)
      {
         fint   r = 0;
         char   name[FITSLEN+1];
         char   shortname[FITSLEN];
         fchar  Axisname;

         Axisname.a = name;
         Axisname.l = FITSLEN;
         gdsc_name_c( Axisname, g->Setin, &(g->axnum[i]), &r );
         Axisname.a[nelc_c(Axisname)] = '\0';
         sprintf( shortname, " %s", strtok( name, " -" ) ); /* Space and hyphen! */
         strcat( g->axisnames, shortname );
      }
      wkeyf( "BOXAXES=%s", g->axisnames );
   }

   /*-------------------------------*/
   /* Get the units of the ampli-   */
   /* tudes for annotations.        */
   /*-------------------------------*/
   {
      fint    r = 0;
      fchar   Dataunits;

      Dataunits.a = g->dataunits;
      Dataunits.l = FITSLEN;
      gdsd_rchar_c( g->Setin, tofchar("BUNIT"), &setlevel, Dataunits, &r );
      Dataunits.a[nelc_c(Dataunits)]  = '\0';
      if (r < 0)
         strcpy( g->dataunits, "??" );
      wkeyf("DATAUNITS=%s", g->dataunits );
   }

   /*-------------------------------*/
   /* Get the units of the profile  */
   /* axis for annotations.         */
   /*-------------------------------*/
   if (g->tophys)
   {
      fint    r;
      fchar   Axunits;

      Axunits.a = g->axunits;
      Axunits.l = FITSLEN;
      r = axunit_c( g->Setin, &(g->profaxnum), Axunits );
      Axunits.a[nelc_c(Axunits)]  = '\0';
      if (r != 0)
         strcpy( g->axunits, "??" );
   }
   else
   {
      strcpy( g->axunits, "grids" );
   }
   wkeyf("AXUNITS=%s", g->axunits );

   /*--------------------------------------------------*/
   /* 1) If a box does not yet exist and if the user   */
   /* did not enter a box on the command line, then    */
   /* wait for the user to enter a box.                */
   /* 2) If a box does not yet exist and the user      */
   /* entered a box on the command line then force a   */
   /* call to 'getbox'.                                */
   /* 3) If a new set or profile is entered (but not a */
   /* new range.                                       */
   /*--------------------------------------------------*/

   TriggerKey("BOX=");
   TriggerKey("OUTSET=");

#ifdef ERUIT
   if (first)
   {
      first = NO;
      TriggerKey("BOX=");
   }

getbox( NULL, "BOX=", KEYCHANGE, g );     /* pseudo event */
#endif
}



static void physrange( ident id,
                       char  *key,
                       int   code,
                       void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: PHYSRANGE= keyword handler.                       */
/*------------------------------------------------------------*/
{
   globaltype   *g=(globaltype*)arg;
   char          dummytxt[STRLEN+1];
   fchar         Dummytxt;
   fint          r;
   fint          dfault;
   char          *s1, *s2;



   if (g->inactive)
   {
      g->inactive = NO;
      return;
   }

   if (!g->setnameexist)
   {
      error_status( "There is no set name",
                    "Enter name of set first", "" );
      return;
   }
   if (!g->profileexist)
   {
      error_status( "There is no profile axis name entered",
                    "Enter name of profile axis", "" );
      return;
   }

   dfault     = HIDDEN;
   Dummytxt.l = STRLEN;
   Dummytxt.a = dummytxt;
   r = usertext_c( Dummytxt, &dfault, tofchar("PHYSRANGE="), tofchar(" ") );
   Dummytxt.a[r] = '\0';
   s1 = strtok( dummytxt, " ");
   s2 = strtok( NULL, " ");

   wkeyf( "DUMMYKEY=%s %s %s %s", s1, g->axunits, s2, g->axunits );

   {
      fint          boxopt = 16;
      fint          loval, hival;
      fint          showdev = DEBUG;

      dfault     = HIDDEN;
      gdsbox_c( &loval, &hival,
                g->Setin, &(g->profcw),
                &dfault,
                tofchar("DUMMYKEY="),
                tofchar(" "),
                &showdev,
                &boxopt );
      if (boxopt == -1)
      {
          cancel( "UCAMESSAGE=" );
          g->rangeexist = NO;
          return;
      }
      else
      {
         g->plohi[0] = loval;
         g->plohi[1] = hival;
      }
   }

   g->rangeexist = YES;
   g->inactive = YES;
   wkeyf( "GRIDRANGE=%d %d", g->plohi[0], g->plohi[1] );
   getset( NULL, "INSET=", KEYCHANGE, g );                    /* pseudo event */
}



static void gridrange( ident id,
                       char  *key,
                       int   code,
                       void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: GRIDRANGE= keyword handler.                       */
/*------------------------------------------------------------*/
{
   globaltype   *g=(globaltype*)arg;
   char          dummytxt[STRLEN+1];
   fchar         Dummytxt;
   fint          r;
   fint          dfault;


   if (g->inactive)
   {
      g->inactive = NO;
      return;
   }

   if (!g->setnameexist)
   {
      error_status( "There is no set name",
                    "Enter name of set first", "" );
      return;
   }
   if (!g->profileexist)
   {
      error_status( "There is no profile axis name entered",
                    "Enter name of profile axis", "" );
      return;
   }


   dfault     = HIDDEN;
   Dummytxt.l = STRLEN;
   Dummytxt.a = dummytxt;
   r = usertext_c( Dummytxt, &dfault, tofchar("GRIDRANGE="), tofchar(" ") );
   Dummytxt.a[r] = '\0';
   if (r == 0)
   {
      wkeyf( "GRIDRANGE=%d %d", g->profgridlo, g->profgridhi );
      return;
   }
   else
   {
      fint      boxopt = 16;
      fint      loval, hival;
      fint      showdev = DEBUG;             /* Show info only in debug mode */

      dfault = HIDDEN;
      gdsbox_c( &loval, &hival,
             g->Setin, &(g->profcw),
             &dfault,
             tofchar( "GRIDRANGE=" ),
             tofchar(" "),
             &showdev,
             &boxopt );

      if (boxopt == -1)
      {
          cancel( "UCAMESSAGE=" );
          g->rangeexist = NO;
          return;
      }
      else
      {
         g->plohi[0] = loval;
         g->plohi[1] = hival;
      }
   }

   /* Update the physical coordinates input field. */
   {
      double physlo, physhi;
      int    r;
      r = grid2phys2( g->Setin, g->plohi[0], g->profcw, g->profaxnum, &physlo);
      r = grid2phys2( g->Setin, g->plohi[1], g->profcw, g->profaxnum, &physhi);
      g->inactive = YES;
      wkeyf( "PHYSRANGE=%g %g", physlo, physhi );
   }

   getset( NULL, "INSET=", KEYCHANGE, g );     /* pseudo event */
}




static void getsetest(  ident  id,
                        char  *key,
                        int    code,
                        void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: ESTSET= keyword handler.                          */
/* Ask name of a set and check for existence. This set is an  */
/* alternative set with which estimates are made.             */
/*------------------------------------------------------------*/
{
   fint          r;
   fint          dfault;
   fint          nitems;
   char          dummytxt[STRLEN+1];
   fchar         Dummytxt;
   globaltype   *g=(globaltype*)arg;


   clearstatus();
   g->altset = NO;
   if (!g->setexist)
   {
      g->altset = NO;
      error_status( "No input set entered!", "Enter main set before special set", "" );
      return;
   }
   dfault     = HIDDEN;
   Dummytxt.l = STRLEN;
   Dummytxt.a = dummytxt;
   nitems = 1;
   r = usertext_c( Dummytxt, &dfault, tofchar("ESTSET="), tofchar(" ") );

   /*--------------------------------------------------*/
   /* If GgiPrompter is cancelled by the user, this    */
   /* routine is called again with an empty keyword.   */
   /*--------------------------------------------------*/
   if (r == 0)
      return;

   Dummytxt.a[nelc_c(Dummytxt)] = '\0';
   strcpy( g->setestbuf, Dummytxt.a );

   /* To check existence, only the name is needed */
   r = 0;
   if (!gds_exist_c(tofchar(strtok(Dummytxt.a, " ")), &r))
   {
      g->altset = NO;
      error_status( "This set does not exist!", "", "" );
      GgiPrompter( "AGAIN", "Set does not exist!" );
      return;
   }
   else
   {
      fint   maxsubs = MAXSUBSETS;
      fint   showdev = DEBUG;             /* Show info only in debug mode */
      fint   estaxnum[MAXAXES];
      fint   axcount;
      fint   maxaxes = MAXAXES;
      fint   class = 2;
      fint   axesoutsidesub = 1;
      fint   cwlo, cwhi;

      /* Test this set */

      dfault = HIDDEN+NOERRORLOOP;        /* Return immediately when error occurs */
      g->naltsubs = gdsinp_c(
                        g->Setest,
                        g->subest,
                        &maxsubs,
                        &dfault,
                        tofchar( "ESTSET=" ),
                        tofchar(" "),
                        &showdev,        /* Device number (as in ANYOUT). */
                        estaxnum,
                        &axcount,        /* Number of grids on axes in 'axnum' */
                        &maxaxes,
                        &class,
                        &axesoutsidesub );

      if (g->naltsubs <= 0)              /* A reject generated by gdsinp */
      {
         fchar   Errtxt;
         char    errtxt[STRLEN];
         fint    dfault = -1;            /* Special default to get UCA message */
         fint    l;

         Errtxt.a = errtxt; Errtxt.l = STRLEN-1; clearstr(Errtxt);
         g->altset = NO;
         l = usertext_c( Errtxt, &dfault, tofchar("UCAMESSAGE="), tofchar("") );
         Errtxt.a[l] = '\0';
         error_status( Errtxt.a,
                      "Check name of profile axis, range",
                      "and minimum axis length of profile (length > 1)." );
         return;
      }

      g->Setest.a[nelc_c(g->Setest)] = '\0';

      cwlo  = gdsc_fill_c( g->Setest, &(g->subest[0]), g->blo );
      cwhi  = gdsc_fill_c( g->Setest, &(g->subest[0]), g->bhi );

      if (cwlo < 0 || cwhi < 0)
      {
         g->altset = NO;
         error_status( "Box of alternative set is not correct",
                       "Check size of both input sets", " " );
         GgiPrompter( "AGAIN", "Estimate box not ok" );
         return;
      }
      GgiPrompter( "ACCEPT", "" );
      g->altset = YES;
   }
}



static void getsetout(  ident  id,
                        char  *key,
                        int    code,
                        void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: OUTSET= keyword handler.                          */
/* Ask name of a set and check for existence.                 */
/*------------------------------------------------------------*/
{
   fint          r;
   fint          dfault;
   fint          nitems;
   char          dummytxt[STRLEN+1];
   fchar         Dummytxt;
   globaltype   *g=(globaltype*)arg;
   bool          overwrite = NO;



   clearstatus();
   g->outsetexist = NO;
   dfault     = HIDDEN;
   Dummytxt.l = STRLEN;
   Dummytxt.a = dummytxt;
   nitems = 1;
   r = userchar_c( Dummytxt, &nitems, &dfault,  KEY_OUTSET, tofchar(" ") );
   Dummytxt.a[nelc_c(Dummytxt)] = '\0';
   strcpy( g->setoutbuf, Dummytxt.a );
   if (r == 0)
   {
      g->outsetexist = NO;
      return;
   }
   else
   {
      r = 0;
      if (gds_exist_c(Dummytxt, &r))
      {
         error_status( "Output set already exists!", "", "" );
         overwrite = GgiVerify("Set exists, ok to overwrite?", NULL, NULL );
         if (overwrite)
         {
            r = 0;
            gds_delete_c( Dummytxt, &r );     /* Delete set, do not overwrite */
            g->outsetexist = YES;
            error_status( "", "", "" );
         }
         else
         {
            cancel_c( KEY_OUTSET );
            g->outsetexist = NO;
         }
      }
      else
      {
         /* Set does not yet exists, but is created later */
         g->outsetexist = YES;
      }      
   }
}




static void getsetgau(  ident  id,
                        char  *key,
                        int    code,
                        void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: GAUSET= keyword handler.                          */
/* Ask name of a set and check for existence.                 */
/* This is the set to which the fitted profiles are written,  */
/* i.e. not the parameters, but the fit data.                 */
/*------------------------------------------------------------*/
{
   fint          r;
   fint          dfault;
   fint          nitems;
   char          dummytxt[STRLEN+1];
   fchar         Dummytxt;
   globaltype   *g=(globaltype*)arg;
   bool          overwrite = NO;



   clearstatus();
   g->gausetexist = NO;
   dfault     = HIDDEN;
   Dummytxt.l = STRLEN;
   Dummytxt.a = dummytxt;
   nitems = 1;
   r = userchar_c( Dummytxt, &nitems, &dfault,  KEY_GAUSET, tofchar(" ") );
   Dummytxt.a[nelc_c(Dummytxt)] = '\0';
   strcpy( g->setgaubuf, Dummytxt.a );
   if (r == 0)
   {
      g->gausetexist = NO;
      return;
   }
   else
   {
      r = 0;
      if (gds_exist_c(Dummytxt, &r))
      {
         error_status( "Output set for fitted profiles already exists!", "", "" );
         overwrite = GgiVerify("Set exists, ok to overwrite?", NULL, NULL );
         if (overwrite)
         {
            r = 0;
            gds_delete_c( Dummytxt, &r );     /* Delete set, do not overwrite */
            g->gausetexist = YES;
            error_status( "", "", "" );
         }
         else
         {
            cancel_c( KEY_OUTSET );
            g->gausetexist = NO;
         }
      }
      else
      {
         /* Set does not yet exists, but will be created later */
         g->gausetexist = YES;
      }      
   }
}




static void getqs( ident  id,
                   char  *key,
                   int    code,
                   void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: Keyword handler for Q=                            */
/*------------------------------------------------------------*/
{
   globaltype *g = (globaltype*) arg;
   fint       nitems  = MAXQ;
   fint       dfault  = HIDDEN;
   int        i;


   g->numq = userint_c( g->smoothingpar,
                        &nitems,
                        &dfault,
                        tofchar(key),
                        tofchar(" ") );


   /* Smoothing parameters must be > 1 */
   for (i = 0; i < g->numq; i++)
   {
      if (g->smoothingpar[i] < 1)
         g->smoothingpar[i] = 1;
   }
}




static void getfunction( ident  id,
                         char  *key,
                         int    code,
                         void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: Keyword handler for FUNCTION=                     */
/*------------------------------------------------------------*/
{
   globaltype *g=(globaltype*)arg;
   fint       userval;
   fint       nitems = 1;
   fint       dfault = HIDDEN;


   (void) userint_c( &userval,
                     &nitems,
                     &dfault,
                     tofchar(key),
                     tofchar(" ") );
   if (userval < 0 || userval >= MAXFUNCTIONS)
   {
      rejectf( key, "Invalid function selection!" );
      g->function = GAUSS;
      g->ncomp = 1;
      return;
   }
   g->function = (int) userval;
   if (g->function == GAUSS)
      g->nvar = 3;              /* Amp, disp, center */
   else if (g->function == GAUSSHERMITE1)
      g->nvar = 4;
   else if (g->function == GAUSSHERMITE2)
      g->nvar = 5;
   else
      g->nvar = 4;
}




static void getncomp( ident  id,
                      char  *key,
                      int    code,
                      void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: NCOMP= keyword handler for (f)ints.               */
/*------------------------------------------------------------*/
{
   globaltype *g = (globaltype*) arg;
   fint       userval = 1;
   fint       nitems  = 1;
   fint       dfault  = HIDDEN;


   (void) userint_c( &userval,
                     &nitems,
                     &dfault,
                     tofchar(key),
                     tofchar(" ") );

   g->ncomp = userval;
}



static void showhelp( ident  id,
                      char  *key,
                      int    code,
                      void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: HELP= keyword handler.                            */
/* Display the dc1 document or other help documents with the  */
/* document viewer.                                           */
/*------------------------------------------------------------*/
{
   fint         indx   = -1;
   fint         nitems =  1;
   fint         dfault = HIDDEN;


   (void) userint_c( &indx,
                     &nitems,
                     &dfault,
                     tofchar(key),
                     tofchar(" ") );

   if (indx < 0)                                                   /* No item */
      return;

   switch( indx )
   {
      case 0:
      {
         GgiShowDoc(NULL);
         break;
      }
      case 1:
      {
         /* The text line was selected */
         break;
      }
      case 2:
      {
         char       path[FILENAME_MAX];
         char       filename[FILENAME_MAX];
         FILE       *fp;

         strcpy( path, getenv("gip_tsk") );

         (void) sprintf( filename, "%s/%s.dc1", path, taskname );
         fp = fopen( filename, "r" );
         if (fp == NULL)
         {
            (void) sprintf( filename, "./%s.dc1", taskname );
            fp = fopen( filename, "r" );
            if (fp == NULL)
            {
               error_status( "Cannot open dc1 document for help", filename, "" );
               return;
            }
         }
         GgiShowDoc( filename );
         fclose( fp );
         break;         
      }
      case 3:
      {
         char       path[FILENAME_MAX];
         char       filename[FILENAME_MAX];
         FILE       *fp;

         strcpy( path, getenv("gip_root") );

         (void) sprintf( filename, "%s/doc/ggi.doc", path, taskname );
         fp = fopen( filename, "r" );
         if (fp == NULL)
         {
            error_status( "Cannot open document about Gui:", filename, "" );
            return;
         }
         GgiShowDoc( filename );
         fclose( fp );
         break;         
      }      
      case 4:
      {
         char       path[FILENAME_MAX];
         char       filename[FILENAME_MAX];
         FILE       *fp;

         strcpy( path, getenv("gip_root") );

         (void) sprintf( filename, "%s/doc/input.doc", path, taskname );
         fp = fopen( filename, "r" );
         if (fp == NULL)
         {
            error_status( "Cannot open input document:", filename, "" );
            return;
         }
         GgiShowDoc( filename );
         fclose( fp );
         break;         
      }      
   }
}




static void quit( ident  id,
                  char  *key,
                  int    code,
                  void  *arg)
/*------------------------------------------------------------*/
/* PURPOSE: QUIT= keyword handler.                            */
/*------------------------------------------------------------*/
{
   static bool quitting=FALSE;
   globaltype  *g = (globaltype*) arg;


/*   bool quit = UserLog(key);*/
   if (!quitting && UserLog("GO=")) {
      quitting = TRUE;
      ScheduleKeyevent(quit, "GO=", KEYCHANGE, arg);
      gui_error("Wait until processing finished . . .");
      return;
   }
   if (1)
   {
      /* Release memory */
      if (g->xphysscale != NULL)  free( g->xphysscale );
      if (g->wdat != NULL)  free( g->wdat );
      if (g->filtered != NULL)  free( g->filtered );
      if (g->work != NULL)  free( g->work );
      if (g->xphys != NULL)  free( g->xphys );
      if (g->xgrids != NULL) free( g->xgrids );
      cancel( key );                       /* remove QUIT= keyword */
      finis_c();                           /* terminate task */
   }
   wkeyf( key );
}




static void initializestruct( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Initialize values for current status.             */
/*------------------------------------------------------------*/
{
   g->Setin.a = g->setbuf;
   g->Setin.l = STRLEN - 1;
   g->Setout.a = g->setoutbuf;
   g->Setout.l = STRLEN - 1;
   g->Setgau.a = g->setgaubuf;
   g->Setgau.l = STRLEN - 1;   
   g->Setest.a = g->setestbuf;
   g->Setest.l = STRLEN - 1;
   g->setexist = NO;
   g->setnameexist = NO;
   g->outsetexist = NO;
   g->gausetexist = NO;   
   g->profileexist = NO;
   g->boxexist = NO;
   g->inactive = NO;
   g->rangeexist = NO;
   g->blo[0] = g->blo[1] = 0;    /* Low  edge of box in grids */
   g->bhi[0] = g->bhi[1] = 0;    /* High edge of box in grids */
   g->profilelen = 0;                 /* Number of input subsets */
   g->imagesize = 0;
   g->profiles = NULL;
   g->results = NULL;
   g->estfromset = NULL;
   g->altset = NO;
   g->wdat = NULL;
   g->work = NULL;
   g->filtered = NULL;
   g->currentprofilenumber = 0;
   g->xgrids = NULL;
   g->xphys  = NULL;
   g->xphysscale  = NULL;
   g->spline = NULL;
   g->median = NULL;
   g->tophys = NO;
   g->Yres = NULL;
   g->Yfit = NULL;   
   g->sortopt = 0;
   g->gridspac = 0.0;
   strcpy( g->setname, "" );
   strcpy( g->profaxname, "" );
   strcpy( g->dataunits, "" );
   strcpy( g->axunits, "" );
   strcpy( g->axisnames, "" );
   g->ampperc = 50.0;
   g->tolerance = TOLERANCE;
   g->lab = LAMBDA;
   g->maxits = MAXITERS;
   g->smoothingpar[0] = 2;
   g->numq = 1;
   g->rms = 1.0;
   g->critamp = 0.0;
   g->critdisp = 0.0;
   g->lorentz = 0.0;
   g->h3 = g->h4 = 0.0;
   g->ncomp = 1;
   g->gaussiansfound = 0;
   g->function = GAUSS;                 /* Default fit a 1 comp standard Gauss */
   g->nvar = 3;                         /* 3 variable parameters */
   g->estonly = tobool( NO );
   g->fitpars = 4;                      /* Total number of vars per comp in output */
   g->errpars = 4;                      /* And their errors */
   g->filter.amp[0] = g->filter.amp[1] = dblank;
   g->filter.vel[0] = g->filter.vel[1] = dblank;
   g->filter.disp[0] = g->filter.disp[1] = dblank;
}



static void giveinstructions( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Display basic input in gui log.                   */
/*------------------------------------------------------------*/
{
   char   mes[STRLEN];
   int    line = 1;

   sprintf( mes, "            XGAUFIT version %s", RELEASE );
   gui_log( mes, line++ );
   gui_log( "               BASIC RECIPE",  line++ );
   gui_log( " 1. Enter name of set with profiles in 'Setname'",  line++ );
   gui_log( " 2. Select profile axis in 'Profile axis' menu",  line++ );
   gui_log( " 3. Adjust default profile length in 'Grid/Phys range'",  line++ );
   gui_log( " 4. Press Enter in 'Box' for default box",  line++ );
   gui_log( " 5. Enter name of an output set", line++ );
   gui_log( " 6. Select function type with 'FUNCTION' menu",  line++ );
   gui_log( " 7. Press 'GO/STOP' to start fitting",  line++ );
   gui_log( "    (Use same button to abort and wait for event)",  line++ );
   gui_log( " 8. Use output in Hermes log to identify subsets",  line++ );
   gui_log( " 9. View results with task VIEW",  line++ );
   gui_log( "10. Exit program with FILE/Exit",  line++ );   
}



MAIN_PROGRAM_ENTRY
/*-------------------------------------------------------------------------*/
/* The macro MAIN_PROGRAM_ENTRY replaces the C-call main() to start the    */
/* main body of your GIPSY application. Variables defined as 'fchar' start */
/* with a capital.                                                         */
/*-------------------------------------------------------------------------*/
{
   static globaltype   glob;                /* Big, so do not use the stack */


   init_c();                                /* Contact Hermes */

   strcpy( taskname, "xgaufit" );           /* Initialze global variables */
   setfblank_c( &fblank );
   setdblank_c( &dblank );
   setlevel = 0;
   maxarg = log( sqrt(DBL_MAX) );           /* Limit range in exp function */

   initializestruct( &glob );               /* Set default properties */

   gui();                                   /* Start graphical user interface */

   giveinstructions( &glob );

   /* Start scheduling the functions */
   (void) ScheduleKeyevent( quit,        "QUIT=",     KEYCHANGE, &glob );
   (void) ScheduleKeyevent( quit,        "FILE=",     KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getsetname,  "SETNAME=",  KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getprofile,  "PROFILE=",  KEYCHANGE, &glob );
/*   (void) ScheduleKeyevent( getset,      "INSET=",    KEYCHANGE, &glob );*/
   (void) ScheduleKeyevent( physrange,   "PHYSRANGE=",KEYCHANGE, &glob );
   (void) ScheduleKeyevent( gridrange,   "GRIDRANGE=",KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getbox,      "BOX=",      KEYCHANGE, &glob );
   (void) ScheduleKeyevent( processdata, "GO=",       KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getncomp,    "NCOMP=",    KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getfunction, "FUNCTION=", KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getsetout,   "OUTSET=",   KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getsetgau,   "GAUSET=",   KEYCHANGE, &glob );   
   (void) ScheduleKeyevent( getsetest,   "ESTSET=",   KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getdble,     "CRITAMP=",  KEYCHANGE, &(glob.critamp) );
   (void) ScheduleKeyevent( getdble,     "CRITDISP=", KEYCHANGE, &(glob.critdisp) );
   (void) ScheduleKeyevent( showhelp,    "HELP=",     KEYCHANGE, NULL );
   (void) ScheduleKeyevent( getdble,     "TOLERANCE=",KEYCHANGE, &(glob.tolerance) );
   (void) ScheduleKeyevent( getdble,     "LAB=",      KEYCHANGE, &(glob.lab) );
   (void) ScheduleKeyevent( getint,      "MAXITS=",   KEYCHANGE, &(glob.maxits) );
   (void) ScheduleKeyevent( getqs,       "Q=",        KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getdble,     "RMS=",      KEYCHANGE, &(glob.rms) );
   (void) ScheduleKeyevent( getampfilter,"AMPRANGE=", KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getvelfilter,"VELRANGE=", KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getdispfilter,"DISPRANGE=", KEYCHANGE, &glob );
   (void) ScheduleKeyevent( getdble,     "VSYS=",     KEYCHANGE, &(glob.vsys) );
/*   (void) ScheduleKeyevent( setsort,     "SORT=",     KEYCHANGE, &glob );*/


   /* Initialize some fields */
   if (!TriggerKey("TOLERANCE="))  wkeyf( "TOLERANCE=%g", glob.tolerance );
   if (!TriggerKey("LAB="))        wkeyf( "LAB=%g", glob.lab );
   if (!TriggerKey("NCOMP="))      wkeyf( "NCOMP=%d", glob.ncomp );
   if (!TriggerKey("FUNCTION="))   wkeyf( "FUNCTION=%d", glob.function );
   if (!TriggerKey("Q="))          wkeyf( "Q=1 2 3 4 5" );
   if (!TriggerKey("MAXITS="))     wkeyf( "MAXITS=%d", glob.maxits );
   if (!TriggerKey("CRITDISP="))   wkeyf( "CRITDISP=%g", glob.critdisp );
   if (!TriggerKey("CRITAMP="))    wkeyf( "CRITAMP=%g", glob.critamp );
   if (!TriggerKey("RMS="))        wkeyf( "RMS=%g", glob.rms );
   if (!TriggerKey("AMPPERC="))    wkeyf( "AMPPERC=%g", glob.ampperc );

   /* A function is known, fill the parameter fields with background pars. */

   if (!TriggerKey("Z2_EST="))               /* Quadratic term */
      wkeyf( "Z2_EST=%g", 0.0 );
   if (!TriggerKey("Z1_EST="))               /* Linear term */
      wkeyf( "Z1_EST=%g", 0.0 );
   if (!TriggerKey("Z0_EST="))               /* Zero level */
      wkeyf( "Z0_EST=%g", 0.0 );

   if (!TriggerKey("Z2_FIX="))               /* Linear term initially fixed */
      wkeyf( "Z2_FIX=YES" );
   if (!TriggerKey("Z1_FIX="))               /* Quadratic term initially fixed */
      wkeyf( "Z1_FIX=YES" );
   if (!TriggerKey("Z0_FIX="))               /* Constant term initially fixed */
      wkeyf( "Z0_FIX=YES" );

   TriggerKey("AMPRANGE=");
   TriggerKey("VELRANGE=");
   TriggerKey("DISPRANGE=");


   /* Some events are forced with the trigger_key routine */
   /* Order is important! */


   TriggerKey("SETNAME=");
   TriggerKey("PROFILE=");


   /* Create a help menu */
   {
      static char *names[] = { "Close doc viewer",
                               "_LINE",
                               "Application",
                               "Gui",
                               "General input",
                                NULL }; 
      helpmenu( names );
   }

   eventmonitor("MONITOR=", TRUE);
   MainLoop();

   return( EXIT_SUCCESS );
}
#<

#>            gui.c
/* gui.c -XT
                              COPYRIGHT (c) 1999
                        Kapteyn Astronomical Institute
         University of Groningen  -  9700 AV Groningen, The Netherlands

Author: J.P. Terlouw
*/

#define STATUSLINES    2            /* number of status lines */
#define LOGLINES      22            /* number of log lines */
#define INPWIDTH     400            /* input box width */

#include "stdio.h"
#include "stddef.h"
#include "gipsyc.h"
#include "string.h"
#include "userfio.h"
#include "ggi.h"
#include "param.h"
#include "fitopt.h"
#include "filtopt.h"
#include "gui.h"


#define Vsys_H "Reference value (phys. coords.) for sorting.\n\
Default value is the velocity at the center of the profile."

#define OUTSET_H "Name (only) of set to store fit parameters"
#define GAUSET_H "Name (only) of set to store fitted data"
#define ESTSET_H "Optional. Name and subset specification for subsets from\n\
which you want to read estimates e.g.: myestset param 1 2 3"

static char *file_entries[]={"Estimate Set",
                             "_LINE", 
                             "Exit", NULL};
                             
typedef enum                {ESTSET,
                             LINE,
                             EXIT} file_choices;

static char *func_entries[]={"Gauss",
                             "Gauss Hermite (h3)",
                             "Gauss Hermite (h3, h4)",
                             "Voigt", NULL};

static char *sort_entries[]={"None",
                             "Peak",
                             "Disp",
                             "Vsys", NULL};

static char *ncomp_entries[]={"_LINE","1","2","3","4","5","6","7","8",NULL};
static char *defaxes[]={" ", NULL};
static char **axnames=defaxes;

static ident l_status[STATUSLINES], l_error;
static ident l_log[LOGLINES];
static ident msort, axis, help;

/* -------------------------------------------------------------------------- */
/*                                 setnames                                   */
/* -------------------------------------------------------------------------- */
/*  set name selector.
 */
static char *setnames(char *name)
{
   int  namlen=strlen(name);
   static char result[FILENAME_MAX];

   if (namlen<7) return NULL;
   if (strcmp(name+namlen-6, ".descr")) return NULL;
   strcpy(result, name);
   result[namlen-6] = '\0';
   return result;
}  

/* -------------------------------------------------------------------------- */
/*                                 set_browse                                 */
/* -------------------------------------------------------------------------- */
static void set_browse(ident id, char *key, int code, void *arg)
{
   char *fkey=(char*)arg;
   bool button=toflog(FALSE);
   
   (void)userflog(&button, 1, 2, key, " ");
   if (tobool(button)) {
      ident browser=GgiFileBrowser(ggiIdent, fkey);
      wkeyf(key);
      if (browser) {
         (void)GgiFileBrowser(ggiDelete, browser);
      } else {
         (void)GgiFileBrowser(ggiCreate, fkey, "Name of set", setnames);
      }
   }
}  

/* -------------------------------------------------------------------------- */
/*                                 statusbox                                  */
/* -------------------------------------------------------------------------- */
/*  statusbox() creates form containing error and status lines
 */
static ident statusbox(int width)
{
   int   i;
   ident curform, prvform, topref=NULL;
   
   curform = GgiForm("status", 1);
   prvform = GgiUseShell(curform);
   l_error  = GgiSetLabel(GgiLabel("error"), " ", 0);
   GgiSetPosition(l_error, 0, NULL,  0, topref);
   topref = l_error;
   for (i=0; i<STATUSLINES; i++) {
      l_status[i] = GgiSetLabel(GgiLabel("status"), " ", 0);
      GgiSetPosition(l_status[i], 0, NULL,  -1, topref);
      topref = l_status[i];
   }
   GgiSetPosition(GgiCanvas(" ",1,1), width, NULL, 0, NULL); /* reserve space */
   (void)GgiUseShell(prvform);
   return curform;
}

/* -------------------------------------------------------------------------- */
/*                                 logbox                                     */
/* -------------------------------------------------------------------------- */
/*  statusbox() creates form containing error and status lines
 */
static ident logbox(int width)
{
   int   i;
   ident curform, prvform, topref=NULL;
   
   curform = GgiForm("log", 1);
   prvform = GgiUseShell(curform);
   for (i=0; i<LOGLINES; i++) {
      l_log[i] = GgiSetLabel(GgiLabel("log"), " ", 0);
      GgiSetPosition(l_log[i], 0, NULL,  -1, topref);
      topref = l_log[i];
   }
   GgiSetPosition(GgiCanvas(" ",1,1), width, NULL, 0, NULL); /* reserve space */
   (void)GgiUseShell(prvform);
   return curform;
}

/* -------------------------------------------------------------------------- */
/*                                 handle_ncomp                               */
/* -------------------------------------------------------------------------- */
static void handle_ncomp(ident id, char *key, int code, void *arg)
{
   fint ncomp=0;
    
   (void)userfint(&ncomp, 1, 2, key, " ");
   if (ncomp>1) GgiActivate(msort); else GgiDeactivate(msort);
}

/* -------------------------------------------------------------------------- */
/*                                 handle_menu                                */
/* -------------------------------------------------------------------------- */
static void handle_menu(ident id, char *key, int code, void *arg)
{
   fint choice, nread;
    
   nread = userfint(&choice, 1, 2, key, " ");
   if (nread>0) {
      switch (choice) {
         case ESTSET: {
#if 0
            GgiPrompter("ESTSET=", "Name and subset specification for \n\
estimate input e.g.: myestset param 1 2 3");
#else
            ident composer=GgiInset(ggiIdent, "ESTSET=");
            if (composer) {
               (void)GgiInset(ggiDelete, composer);
            } else {
               composer = GgiInset(ggiCreate, "EST_NAME=", "ESTSET=", NULL,
                          "Set/subset for estimate input");
            }
#endif
            break;
         }
         case EXIT: {
            wkeyf("QUIT=YES");
            break;
         }
         default: {
            break;
         }
      }
      wkeyf(key);
   }
}

/* -------------------------------------------------------------------------- */
/*                                 handle_axis                                */
/* -------------------------------------------------------------------------- */
static void handle_axis(ident id, char *key, int code, void *arg)
{
   fint choice, nread;
    
   nread = userfint(&choice, 1, 2, key, " ");
   if (nread>0) {
      wkeyf("PROFILE=%s", axnames[choice]);
      GgiSetLabel(axis, axnames[choice], 0);
   }
}

/* ========================================================================== */
/*                                 gui_set_units                              */
/* -------------------------------------------------------------------------- */
/*  Keyword handler for setting units-like strings.
 */
extern void gui_set_units(ident id, char *key, int code, void *arg)
{
   fchar ftext;
   char  ctext[41];
   char  ltext[43];
   int   nread;
   
   ftext.a = ctext; ftext.l = 40;
   nread = userftext(ftext, 2, key, " ");
   ctext[nread] = '\0';
   sprintf(ltext,"(%s)", ctext);
   GgiSetLabel((ident)arg, ltext, 0);
}

/* ========================================================================== */
/*                                 gui                                        */
/* -------------------------------------------------------------------------- */
/*  Interface routine to initialize the graphical user interface.
 */
extern void gui(void)
{
static char *my_resources[] = {
#if 1
   "*Background: #D8D8D8",
   "*logo.foreground: #B8B8B8",
   "*TextField.background: #FAF8CC",
#endif
   "*error.foreground: red",
   "*log.font: *fixed-bold-r-normal*13*",
   NULL
};

   int   barwidth=2*INPWIDTH+9;
   
   ident file, filtopt, fitopt, param, func, ncomp;
   ident inset, selset, outset, gauset, laxis, gridrange, physrange, axunits,
         box, boxaxes, go,
         vsys;
   ident menubar, status, log, inputs, prbar, logo;
   char  *labels[]={"", NULL};
   
   wkeyf("PROGRESS=0.0");

   GgiAutoLayout(FALSE);
   GgiPostponeRealize(TRUE);
   GgiSetResources(my_resources);

/* ----------------------- top level elements ------------------------------ */
   menubar = GgiForm("top", 1);
   status  = statusbox(INPWIDTH-1);
   log     = logbox(INPWIDTH-1);
   inputs  = GgiForm("inputs", 1);
   prbar   = GgiProgress("PROGRESS=", "Progress", 300);
   logo    = GgiLogo(NULL, 0, 0);    
   
   GgiSetPosition(menubar, 0, NULL, 0, NULL);
   GgiSetPosition(log,     0, NULL, 0, menubar);
   GgiSetPosition(inputs,  0, log,  0, menubar);
   GgiSetPosition(status,  0, log,  0, inputs);
   GgiSetPosition(prbar,   4, log, -GgiHeight(prbar)-2, log);
   GgiSetPosition(logo, -GgiWidth(logo), menubar, -GgiHeight(logo), log);

/* ---------------------------- menu bar ----------------------------------- */
   (void)GgiUseShell(menubar);
   file   = GgiSetLabel(GgiMenu("GUI_FILE=", NULL, file_entries), "FILE", 0);
   (void)ScheduleKeyevent(handle_menu, "GUI_FILE=", KEYCHANGE, NULL);
   help   = GgiMenu("HELP=", NULL, labels);
   filtopt = GgiSetLabel(GgiButton("FILTOPT=", "Filter options"),
                        "FILTERS", 0);                        
   fitopt = GgiSetLabel(GgiButton("OPTIONS=", "Miscellaneous fit options"),
                        "FIT OPTIONS", 0);
   param  = GgiButton("PARAM=", "Parameter window");
   func   = GgiMenu("FUNCTION=", "Function menu", func_entries);
   ncomp  = GgiMenu("NCOMP=", "Number of components", ncomp_entries);
   msort  = GgiMenu("SORT=", "Sort options", sort_entries);
   
   GgiSetPosition(file, 0, NULL, 0, NULL);
   GgiSetPosition(help, barwidth-GgiWidth(help), NULL, 0, NULL);
   GgiSetPosition(msort, -GgiWidth(help)-GgiWidth(msort)-4, help, 0, NULL);
   GgiSetPosition(ncomp, -GgiWidth(msort)-GgiWidth(ncomp)-4, msort, 0, NULL);
   GgiDeactivate(msort);
   (void)ScheduleKeyevent(handle_ncomp, "NCOMP=", KEYCHANGE, NULL);
   GgiSetPosition(func, -GgiWidth(ncomp)-GgiWidth(func)-4, ncomp, 0, NULL);
   GgiSetPosition(param, -GgiWidth(func)-GgiWidth(param)-4, func, 0, NULL);
   GgiSetPosition(fitopt, -GgiWidth(param)-GgiWidth(fitopt)-4, param, 0, NULL);
   GgiSetPosition(filtopt, -GgiWidth(fitopt)-GgiWidth(filtopt)-4, fitopt, 0,NULL);   
   (void)GgiUseShell(NULL);

/* ------------------------- inputs box ------------------------------------ */
   (void)GgiUseShell(inputs);
   inset = GgiSetLabel(GgiTextField("SETNAME=", "Name of set", 10),
                       "Setname", 80);
   selset = GgiSetLabel(GgiButton("SELSET=", NULL), "BROWSE", 0);
   (void)ScheduleKeyevent(set_browse, "SELSET=", KEYCHANGE, "SETNAME=");

   laxis = GgiLabel("Profile axis");
   axis  = GgiSetLabel(GgiMenu("PROFINDX=", "Profile axis name", axnames),
                       "(none)", 0);
   (void)ScheduleKeyevent(handle_axis, "PROFINDX=", KEYCHANGE, NULL);
   gridrange = GgiSetLabel(GgiTextField("GRIDRANGE=", "Profile delimeters in x in grids\n Enter as: -30  32", 10),
                           "Grid range", 80);
   physrange = GgiSetLabel(GgiTextField("PHYSRANGE=", "Profile delimeters in x in physical coordinates.\n Enter as: 200.1  660.4" , 10),
                           "Phys. range", 80);
   axunits = GgiLabel("( )");
   (void)ScheduleKeyevent(gui_set_units, "AXUNITS=", KEYCHANGE, axunits);
   box      = GgiSetLabel(GgiTextField("BOX=", "Each position in this box contains a profile.\n Enter as: position low  position high", 20), "Box", 80);
   boxaxes  = GgiLabel("( )");
   (void)ScheduleKeyevent(gui_set_units, "BOXAXES=", KEYCHANGE, boxaxes);
   
   go       = GgiSetLabel(GgiButton("GO=", NULL), "GO / STOP", 0);
   vsys     = GgiSetLabel(GgiTextField("VSYS=", Vsys_H, 10), "Vsys", 80);

   outset   = GgiSetLabel(GgiTextField("OUTSET=", OUTSET_H, 20),
                          "Parameter output", 140);   
   gauset   = GgiSetLabel(GgiTextField("GAUSET=", GAUSET_H, 20),
                          "Data output", 140);   
    
   GgiSetPosition(inset, 0, NULL, 0, NULL);
   GgiSetPosition(selset, 0, inset, 0, NULL);
   GgiSetPosition(laxis, 0, selset, 0, NULL);
   GgiSetPosition(axis, 0, laxis, 0, NULL);
   
   GgiSetPosition(gridrange, 0, NULL, 0, inset);
   GgiSetPosition(physrange, 0, gridrange, 0, inset);
   GgiSetPosition(axunits, 0, physrange, 0, inset);
   
   GgiSetPosition(box, 0, NULL, 0, gridrange);
   GgiSetPosition(boxaxes, 0, box,  0, gridrange);
   
   GgiSetPosition(vsys,   0, NULL, 0, box);
   GgiSetPosition(outset, 0, NULL, 0, vsys);
   GgiSetPosition(gauset, 0, NULL, 0, outset);
   
   GgiSetPosition(go, INPWIDTH-GgiWidth(go)-1, NULL, 0, box);
   
   GgiUseShell(NULL);

   Param("PARAM=");
   Fitopt("OPTIONS=");
   Filtopt("FILTOPT=");

   GgiRealize();
}

/* ========================================================================== */
/*                                 gui_status                                 */
/* -------------------------------------------------------------------------- */
/*  Interface routine to display a status message.   
 */
void gui_status(char *message, int lineno)
{
   GgiSetLabel(l_status[lineno%STATUSLINES], message, 0); 
}

/* ========================================================================== */
/*                                 gui_log                                    */
/* -------------------------------------------------------------------------- */
/*  Interface routine to display a status message.   
 */
void gui_log(char *message, int lineno)
{
   char buffer[512];
   
   sprintf(buffer, "%s", message);
   GgiSetLabel(l_log[lineno%LOGLINES], buffer, 0); 
}

/* ========================================================================== */
/*                                 gui_error                                  */
/* -------------------------------------------------------------------------- */
/*  Interface routine to display an error message.
 */
void gui_error(char *message)
{
   char buffer[512];
   
   if (*message) {
      sprintf(buffer,"Error: %s",message);
      GgiSetLabel(l_error, buffer, 0);
   } else {
      GgiSetLabel(l_error, " ", 0);
   }
}

/* ========================================================================== */
/*                                 gui_axnames                                */
/* -------------------------------------------------------------------------- */
/*  Interface routine to change the axis names menu.
 */
extern void gui_axnames(char **names)
{
   fchar curname;
   char  curnamec[41];
   fint  nread;
   
   curname.a = curnamec;
   curname.l = 40;
   axnames = names;
   (void)GgiMenuLabels(axis, axnames);
   nread = userftext(curname, 2, "PROFILE=", " ");
   if (nread>0) {
      curnamec[nread] = '\0';
      GgiSetLabel(axis, curnamec, 0);
   } else {
      GgiSetLabel(axis, "(none)", 0);
   }
}

/* ========================================================================== */
/*                                 helpmenu                                   */
/* -------------------------------------------------------------------------- */
/*  Set up a help menu.
 */
void helpmenu( char** names )
{  
   GgiMenuLabels(help, names);
}
   

#if defined(TESTBED)  /* ++++++++++++++++++++++++++++++++++++++++++++++++++ */
#include "init.h"
#include "finis.h"
#include "cmain.h"  
#include "eventmonitor.h"

static void setmessage(ident id, void *arg)
{
   static int mesno=0;
   char message[80];

   sprintf(message,"Status message %d (%o Oct, %X Hex)", mesno, mesno, mesno);
   gui_status(message,mesno);
   mesno++;
}

static void seterror(ident id, void *arg)
{
   static int mesno=1;
   char message[80];
   
   sprintf(message,"Error message %d (%o Oct, %X Hex)", mesno, mesno, mesno);
   mesno++;
   if (mesno%2) gui_error(message); else gui_error("");
}

static void setlog(ident id, void *arg)
{
   static int mesno=0;
   char message[80];
   
   sprintf(message,"This is log message number %d ", mesno*mesno);
   gui_log(message, mesno);
   mesno++;
}

static void quit(ident id, char *key, int code, void *arg)
{
   bool leave=toflog(FALSE);
   (void)userflog(&leave, 1, 2, key, " ");
   if (tobool(leave)) {
      wkeyf(key);
      finis_c();
   }
}
  
MAIN_PROGRAM_ENTRY
{
   static char *axes[]={"RA", "DEC", "FREQ", NULL};

   init_c();
   gui();
   (void)ScheduleTimer(setmessage,  4 SECONDS, NULL);
   (void)ScheduleTimer(seterror,   15 SECONDS, NULL);
   (void)ScheduleTimer(setlog,    10 SECONDS,  NULL);
   (void)ScheduleKeyevent(quit, "QUIT=", KEYCHANGE, NULL);
   eventmonitor("MONITOR=", TRUE);
   gui_axnames(axes);
/*   ggiImmediateKeyevent(TRUE); */
   {
      char *labels[]={"DC1 document", "_LINE", "Interaction", "Theory", NULL};
      helpmenu(labels);
   }
   MainLoop();
}
#endif
#<

#>            dofit.c
/* dofit.c
                           COPYRIGHT (c) 1999
                      Kapteyn Astronomical Institute
                University of Groningen, The Netherlands
                           All Rights Reserved.
*/

#include    "xgaufit.h"
#include    "lsqfitd.h"
#include    "gauestd.h"
#include    "time.h"
#include    "utils.h"
#include    "outputresults.h"


extern  float    fblank;
extern  double   dblank;
extern  double   maxarg;      /* Limit range in exp function */



static char *prefix[5][5] = { {"AMP", "CENTER", "DISP", "", ""},   
                              {"AMP", "CENTER", "DISP", "H3", ""},  
                              {"AMP", "CENTER", "DISP", "H3", "H4"},
                              {"AMP", "CENTER", "DISP", "LORENTZ", ""},
                              {"Z0", "Z1", "Z2", "", ""}  };

typedef struct                /* 'qsort' struct for comparisons */  
{
   double par[MAXPAR];
   double err[MAXPAR];
   double vsys;
} 
qsortstruct;




extern void disp2fwhm( double  *disp )
/*------------------------------------------------------------*/
/* PURPOSE: Convert from dispersion to Full Width Half Max.   */
/*------------------------------------------------------------*/
{
   *disp *= 2.0*sqrt(2.0*log(2.0));
}




static double gausshermiteh3( double  X,
                              double *fpar,
                              int     ncomp )
/*------------------------------------------------------------*/
/* PURPOSE: Calculate Hermite polynomial for X.               */
/*                                                            */
/* The function is a superposition of functions of the form   */
/* f =  A.Exp[-(x-xc)^2/(2s^2] *                              */
/*      (1+h3[c1((x-xc)/s)^3-c2(x-xc)/s]                      */
/*                                                            */
/* ==> F(x) = SUM[ f ] + z0 + B*x + A*x*x                     */
/*                                                            */
/* Amplitude        Center      dispersion    skewness        */ 
/* A = fpar[0],  X0 = fpar[1],  s = fpar[2], h3 = fpar[3]     */
/*                                                            */
/*------------------------------------------------------------*/
{
   int     i;
   int     maxpar = 4;               /* 4 parameters per component */
   double  result;
   double  X_X0 = 0.0;
   double  c1 = 2.0*sqrt(3.0)/3.0, c2 = sqrt(3.0);   

   result = 0.0;      
   for (i = 0; i < ncomp; i++) 
   {
      double A, X0, s, h3;
      double arg;
      double E, F;

      A   = fpar[0+i*maxpar];
      X0  = fpar[1+i*maxpar];
      s   = fpar[2+i*maxpar];
      h3  = fpar[3+i*maxpar];

      X_X0 = X - X0;
      F = X_X0/s;
       
      if (s == 0.0)
         E = 0.0;
      else
      {
         arg = 0.5 * F * F;      
         if (arg < maxarg)         /* 'maxarg' is a global variable */
            E = exp( -arg );
         else
            E = 0.0;     
      }
      E *= (1.0+h3*F*(c1*F*F-c2));
      result += A * E;
   }
   /* Add background components */
   result += fpar[ncomp*maxpar+0] + fpar[ncomp*maxpar+1]*X_X0 + 
             fpar[ncomp*maxpar+2]*X_X0*X_X0;   
   return( result );
}




static void gausshermiteh3derv( double  X,
                                double *fpar,
                                double *epar,
                                int     ncomp ) 
/*------------------------------------------------------------*/
/* PURPOSE: Calculate the derivatives for a skewed gauss at X */
/*------------------------------------------------------------*/
{
   int     i;
   int     maxpar = 4;               /* 4 parameters per component */
   double  c1 = 2.0*sqrt(3.0)/3.0, c2 = sqrt(3.0);
   double  X_X0 = 0.0;

   
   for (i = 0; i < ncomp; i++) 
   {
      double A, X0, s, h3;
      double arg;
      double E, F, Q;

      A   = fpar[0+i*maxpar];
      X0  = fpar[1+i*maxpar];
      s   = fpar[2+i*maxpar];
      h3  = fpar[3+i*maxpar];

      X_X0 = X - X0;
       
      if (s == 0.0)
         E = 0.0;
      else      
      {
         arg = 0.5 * (X_X0/s) * (X_X0/s);      
         if (arg < maxarg)         /* 'maxarg' is a global variable */
            E = exp( -arg );
         else
            E = 0.0;
      }
      if (E == 0.0)
         F = Q = 0.0;
      else
      {
         F = X_X0 / s;
         Q = F * (c1*F*F - c2);
      }
   
      /* Diff amplitude A */
      epar[0+i*maxpar] = E * (1.0 + h3*Q);
     
      /* Diff center X0 */
      epar[1+i*maxpar] = (A*E/s) *( h3*(c2-3.0*c1*F*F) + F*(1.0+h3*Q) );
         
      /* Diff sigma_0 s0 */
      epar[2+i*maxpar] = F * epar[1+i*maxpar];

      /*  Diff h3 */
      epar[3+i*maxpar] = A*E*Q;
   }      
   /* Derivatives for the background */
   epar[ncomp*maxpar+0] = 1.0;
   epar[ncomp*maxpar+1] = X_X0;
   epar[ncomp*maxpar+2] = X_X0*X_X0;   
}




double gausshermiteh3h4( double  X,
                         double *fpar,
                         int     ncomp )
/*------------------------------------------------------------*/
/* PURPOSE: Calculate Hermite polynomial for X.               */
/*                                                            */
/* The function is a superposition of functions of the form   */
/*                                                            */
/* g = (x-xc)/s                                               */
/* f =  A.Exp[-1/2.g^2]*                                      */
/*      { 1 + h3[c1.g+c3.g^3] + h4[c5+c2.g^2+c4.g^4] }        */
/*                                                            */
/* ==> F(x) = SUM[ f ] + z0 + B*x + A*x*x                     */
/*                                                            */
/* Amplitude        Center      dispersion                    */ 
/* A = fpar[0],  X0 = fpar[1],  s = fpar[2], h3 = fpar[3]     */
/* h4 = fpar[4]                                               */
/*                                                            */
/*------------------------------------------------------------*/
{
   int     i;
   int     maxpar = 5;               /* 5 parameters per component */
   double  result;
   double  X_X0 = 0.0;
   double  c0 = sqrt(6.0)/4.0;
   double  c1 = -sqrt(3.0);      
   double  c2 = -sqrt(6.0);
   double  c3 = 2.0*sqrt(3.0)/3.0;
   double  c4 = sqrt(6.0)/3.0;


   result = 0.0;      
   for (i = 0; i < ncomp; i++) 
   {
      double A, X0, s, h3, h4;
      double arg;
      double E, F;

      A   = fpar[0+i*maxpar];
      X0  = fpar[1+i*maxpar];
      s   = fpar[2+i*maxpar];
      h3  = fpar[3+i*maxpar];
      h4  = fpar[4+i*maxpar];

      X_X0 = X - X0;
      F = X_X0/s;
       
      if (s == 0.0)
         E = 0.0;
      else
      {
         arg = 0.5 * F * F;      
         if (arg < maxarg)         /* 'maxarg' is a global variable */
            E = exp( -arg );
         else
            E = 0.0;     
      }
      E *= ( 1.0 + h3*F*(c3*F*F+c1) + h4*(c0+F*F*(c2+c4*F*F)) );
      result += A * E;
   }
   /* Add background components */
   result += fpar[ncomp*maxpar+0] + fpar[ncomp*maxpar+1]*X_X0 + 
             fpar[ncomp*maxpar+2]*X_X0*X_X0;   
   return( result );
}




static void gausshermiteh3h4derv( double  X,
                                double *fpar,
                                double *epar,
                                int     ncomp ) 
/*------------------------------------------------------------*/
/* PURPOSE: Calculate the derivatives for a skewed gauss at X */
/*------------------------------------------------------------*/
{
   int     i;
   int     maxpar = 5;               /* 5 parameters per component */
   double  X_X0 = 0.0;
   double  c0 = sqrt(6.0)/4.0;
   double  c1 = -sqrt(3.0);      
   double  c2 = -sqrt(6.0);
   double  c3 = 2.0*sqrt(3.0)/3.0;
   double  c4 = sqrt(6.0)/3.0;

   
   for (i = 0; i < ncomp; i++) 
   {
      double A, X0, s, h3, h4;
      double arg;
      double E, F, Q3, Q4;

      A   = fpar[0+i*maxpar];
      X0  = fpar[1+i*maxpar];
      s   = fpar[2+i*maxpar];
      h3  = fpar[3+i*maxpar];
      h4  = fpar[4+i*maxpar];


      X_X0 = X - X0;
       
      if (s == 0.0)
         E = 0.0;
      else      
      {
         arg = 0.5 * (X_X0/s) * (X_X0/s);      
         if (arg < maxarg)         /* 'maxarg' is a global variable */
            E = exp( -arg );
         else
            E = 0.0;
      }
      if (E == 0.0)
         F = Q3 = Q4 = 0.0;
      else
      {
         F = X_X0 / s;
         Q3 = F * (c3*F*F+c1);
         Q4 = c0 + F*F*(c2+c4*F*F);
      }
   
      /* Diff amplitude A */
      epar[0+i*maxpar] = E * (1.0 + h3*Q3 + h4*Q4);
     
      /* Diff center X0 */
      epar[1+i*maxpar] = (A*E/s) *( 
                         -h3*(c1+3.0*c3*F*F) - 2.0*h4*F*(c2+2.0*c4*F*F) 
                         + F*(1.0+h3*Q3+h4*Q4) );
         
      /* Diff sigma */
      epar[2+i*maxpar] = F * epar[1+i*maxpar];

      /*  Diff h3 */
      epar[3+i*maxpar] = A*E*Q3;
      
      /*  Diff h4 */
      epar[4+i*maxpar] = A*E*Q4;     
   }      
   /* Derivatives for the background */
   epar[ncomp*maxpar+0] = 1.0;
   epar[ncomp*maxpar+1] = X_X0;
   epar[ncomp*maxpar+2] = X_X0*X_X0;   
}




static double gauss( double  X,
                     double *fpar,
                     int     ngauss )
/*------------------------------------------------------------*/
/* PURPOSE: Fit routine needs evaluation of 1-dim gauss func- */
/*          tion on arbitrary position.                       */
/*                                                            */
/* The function is a superposition of gaussians of the form   */
/* f = A.Exp[-(x-xc)^2/(2s^2] with a zero level z0 added.     */
/* ==> F(x) = SUM[ f ] + z0 + B*x + A*x*x                     */
/*                                                            */
/* Amplitude        Center      dispersion                    */ 
/* A = fpar[0],  x0 = fpar[1],  s = fpar[2];                  */
/*                                                            */
/* etc.                                                       */
/* and the zero level is stored in fpar[ngauss*MAXPAR], the   */
/* last element of the input array.                           */
/*------------------------------------------------------------*/
{
   int      i;
   double   A, X0, s;   
   double   result = 0.0;
   double   X_X0= 0.0;
   int      maxpar = 3;                 /* 3 components for 1 Gauss */


   for (i = 0; i < ngauss; i++)
   {
      double  E, arg;                  /* Exponent and its argument */
      int     offset = i * maxpar;
      
      A  = fpar[offset];  X0 = fpar[1+offset];  s = fpar[2+offset];
      X_X0 = X - X0;
 
      if (s == 0.0)
         E = 0.0;
      else
      {                              
         arg = 0.5*(X_X0/s)*(X_X0/s);
         if (arg < maxarg)         /* 'maxarg' is a global variable */
            E = exp( -arg );
         else
            E = 0.0;
      }
      result += A * E;
   }
   /* Add background components */
   result += fpar[ngauss*maxpar+0] + fpar[ngauss*maxpar+1]*X_X0 + 
             fpar[ngauss*maxpar+2]*X_X0*X_X0;
   return( result );       
}



static void gaussderv( double  X,
                       double *fpar,
                       double *epar,
                       int     ngauss ) 
/*------------------------------------------------------------*/
/* PURPOSE: Fit routine needs evaluation of 1-dim gauss func- */
/*          tion on arbitrary position.                       */
/*                                                            */
/* The function is a superposition of gaussians of the form   */
/* f = A.Exp[-(x-xc)^2/2s^2] with a zero level z0 added.      */
/* ==> F(x) = SUM[ f ] + z0                                   */
/*                                                            */
/* Amplitude        Center      dispersion                    */ 
/* A = fpar[0],  x0 = fpar[1],  s = fpar[2];                  */
/*                                                            */
/* etc.                                                       */
/* and the zero level is stored in fpar[ngauss*MAXPAR], the   */
/* last element of the input array.                           */
/*                                                            */
/* Let E = Exp[-(x-x0)^2(2s^2)], then                         */
/* D[f,A]  =  E                                               */
/* D[f,x0] =  A.(x-x0).E / s^2                                */
/* D[f,a]  =  A.(x-x0)^2.E / s^2                              */
/* D[F,z0] =  1                                               */
/*------------------------------------------------------------*/
{
   int      i;
   double   A, X0, X_X0=0.0, s;   
   int      maxpar = 3;                 /* 3 components for 1 Gauss */


   for (i = 0; i < ngauss; i++)
   {
      int     offset = i * maxpar;
      double  E, arg;                  /* Exponent and its argument */
      double  AEX_X0;                  /* A*(X-Xc)*E  */
      
      A  = fpar[offset];  X0 = fpar[1+offset];  s = fpar[2+offset];

      X_X0 = X - X0; 

      if (s == 0.0)
         E = 0.0;
      else
      {
         arg = 0.5*((X_X0)/s)*((X_X0)/s);     
         if (arg < maxarg)
            E = exp( -arg );
         else
            E = 0.0;
      }

      AEX_X0 = A * E * X_X0;
      
      epar[0+offset] = E;            /* Derivative wrt. Amplitude A */
      epar[1+offset] = AEX_X0/s/s;   /* Derv. wrt Center Xc */
      epar[2+offset] = AEX_X0*X_X0/s/s/s;  /* Derv. wrt. dispersion */
   }
   /* Derivatives for the background */
   epar[ngauss*maxpar+0] = 1.0;
   epar[ngauss*maxpar+1] = X_X0;
   epar[ngauss*maxpar+2] = X_X0*X_X0;
}




extern void w( double Rz, 
               double Iz,
               double *Rr,
               double *Ir )
/*------------------------------------------------------------*/
/* PURPOSE: Error function for complex arguments.             */
/*          W(Z) = EXP(-Z*Z)ERFC(-iZ)                         */
/*          see Abromowitz and Stegun (chapter 7)             */
/*          |error|  <  2 * 10e-6                             */          
/*------------------------------------------------------------*/
{
#define EPS     ( 0.0000001 )
   double	wdx = 0.0;		/* real part result */
   double	wdy = 0.0;		/* imag part result */
   double	x, y;			/* real and imaginary arg */

   x = fabs( Rz );                      /* real part */
   y = fabs( Iz );                      /* imaginairy part */
   if ( x > 6.0 || y > 6.0 ) {		/* approximation for large arguments */
      static double	at[] = { 0.5124242, 0.05176536 };
      static double	bt[] = { 0.2752551, 2.72474500 };
      int		i;

      for ( i = 0; i < 2; i++ ) {
         double	det, dtx, dty, sav;

         sav = ( x * x - y * y - bt[i] );
         det = ( sav * sav + 4.0 * x * x * y * y );
         dtx = ( 2.0 * x * x * y - y * sav ) * at[i];
         dty = ( 2.0 * x * y * y + x * sav ) * at[i];
         wdx += dtx / det;
         wdy += dty / det;
      }
   } else if ( x > 3.9 || y > 3.0 ) {	/* approximation for intermediate arguments */
      static double	at[] = { 0.4613135, 0.09999216, 0.002883894 };
      static double	bt[] = { 0.1901635, 1.78449270, 5.525343700 };
      int		i;

      for ( i = 0; i < 3; i++ ) {
         double	det, dtx, dty, sav;

         sav = ( x * x - y * y - bt[i] );
         det = ( sav * sav + 4.0 * x * x * y * y );
         dtx = ( 2.0 * x * x * y - y * sav ) * at[i];
         dty = ( 2.0 * x * y * y + x * sav ) * at[i];
         wdx += dtx / det;
         wdy += dty / det;
      }
   } else {				/* no approximation */
      double	r;

      wdx = 1.0;
      r = sqrt( x * x + y * y );
      if ( r > 0.0 ) {
         static double	tt[] = { 1.0000000000, 0.5641895835 };
         double		tn[2];
         double		del, csp, snp, tcn, tsn;
         int		n = 0;

         csp = -y / r;
         snp = x / r;
         tcn = 1.0; tsn = 0.0;
         tn[0] = tt[0]; tn[1] = tt[1];
         do {
            double	tc, ts;
            int		i;

            n += 1;				/* increment interation number */
            tc = tcn;				/* save */
            ts = tsn;				/* save */
            tcn = ( tc * csp - ts * snp );	/* next cosine term */
            tsn = ( ts * csp + tc * snp );	/* next sine term */
            i = n%2;				/* argument */
            tn[i] *= ( 2.0 / (double) n );
            tn[0] *= r;				/* multiply with radius */
            tn[1] *= r;				/* multiply with radius */
            del = tn[i];			/* increment */
            wdx += ( tcn * del );		/* add increment */
            wdy += ( tsn * del );		/* add increment */
         } while ( del > EPS );			/* precision is reached */
      }
   }
   if ( Rz >= 0.0 && Iz >= 0.0 ) {
      *Rr = wdx; *Ir = wdy;
   } else if ( Rz >= 0.0 && Iz < 0.0 ) {
      double	csp, snp, sav;

      csp = cos( 2.0 * x * y );
      snp = sin( 2.0 * x * y );
      sav = exp( y * y - x * x );
      *Rr = sav * csp - wdx;
      *Ir = sav * snp + wdy;
   } else if ( Rz < 0.0 && Iz >= 0.0 ) {
      *Rr = wdx;
      *Ir = -wdy;
   } else if ( Rz < 0.0 && Iz < 0.0 ) {
      double	csp, snp, sav;

      csp = cos( 2.0 * x * y );
      snp = sin( 2.0 * x * y );
      sav = exp( y * y - x * x );
      *Rr = sav * csp - wdx;
      *Ir = -sav * snp - wdy;
   }
#undef EPS   
}





static double Voigt( double  X,
                     double *dpar,
                     int     ncomp )
/*------------------------------------------------------------*/
/* PURPOSE: Calculate Voigt function at X.                    */
/*                                                            */
/* We want to use the voigt function from the GIPSY library.  */
/* The parameters of that function are floats and the order   */
/* and meaning of the parameters is different from what we    */
/* are used to. The first parameter is the zero level and     */
/* then for each component the Intensity (?), center, FWHM    */
/* Lorentz factor.                                            */
/*------------------------------------------------------------*/
{
   int      maxpar = 4;                            /* Per component */
   int      i;
   double   x, y;  
   double   result, V;
   double   sqln2 = sqrt( log(2.0) );
   double   Rres, Ires;
   double   X_X0 = 0.0;


   result = 0.0;
   for (i = 0; i < ncomp; i++)
   {
      int      offset = i * maxpar;
      double   Int, X0;
      double   aD, aL;      
      Int = dpar[0+offset];    X0 = dpar[1+offset];  
      aD  = dpar[2+offset];    aL = dpar[3+offset];
      X_X0 = X - X0;
       
      if (aD == 0.0)
         V = 0.0;
      else
      { 
         double ampfct;
         
         x = X_X0 * sqln2 / aD;
         y = aL * sqln2 / aD; 
         w( x, y, &Rres, &Ires );
         ampfct = Int * sqln2/(aD*sqrt(PI));
         V = ampfct * Rres;
      }
      result += V;     
   }
   /* Add background components */
   result += dpar[ncomp*maxpar+0] + dpar[ncomp*maxpar+1]*X_X0 + 
             dpar[ncomp*maxpar+2]*X_X0*X_X0;
   return( result );
}




static void   Voigtderv( double  X,
                         double *dpar,
                         double *epar,
                         int     ncomp )
/*------------------------------------------------------------*/
/* PURPOSE: Calculate Voigt derivatives at X.                 */
/*                                                            */
/*------------------------------------------------------------*/
{
   int      maxpar = 4;                            /* Per component */
   int      i;
   double   x, y;  
   double   sqln2 = sqrt( log(2.0) );
   double   X_X0 = 0.0;
   double   sqpi = sqrt(PI);


   for (i = 0; i < ncomp; i++)
   {
      int      offset = i * maxpar;
      double   Int, X0;
      double   aD, aL;
            
      Int = dpar[0+offset];    X0 = dpar[1+offset];  
      aD  = dpar[2+offset];    aL = dpar[3+offset];
      X_X0 = X - X0;
       
      if (aD == 0.0)
         epar[0+offset] = epar[1+offset] = epar[2+offset] = epar[3+offset] = 0.0; 
      else
      {         
         double Rezwz,Imzwz;
         double ampfct;
         double Rew, Imw;         
         double dxVoigt, dyVoigt;
         
         x = X_X0 * sqln2 / aD;              /* Definitions for x and y */
         y = aL * sqln2 / aD; 
         w( x, y, &Rew, &Imw );
         Rezwz  = ( x*Rew - y*Imw );         /* Real part of x.w(z) */
         Imzwz  = ( x*Imw + y*Rew );         /* Imaginary part of x.w(z) */
         dxVoigt = -2.0 * Rezwz;             /* dVoigt/dx */
         dyVoigt = -2.0/sqpi + 2.0 * Imzwz;  /* dVoigt/dy */
         ampfct = Int * sqln2/(aD*sqpi);     /* Help variable */
         
                 
         /* The partial derivatives: */
         
         epar[0+offset] = sqln2/(aD*sqpi) * Rew;             /* Amplitude */
         epar[1+offset] = ampfct * (-sqln2/aD) * dxVoigt;    /* Center    */
         epar[2+offset] = (-ampfct/aD) * 
                          (Rew + x*dxVoigt+ y*dyVoigt);      /* Doppler factor    */
         epar[3+offset] = ampfct * (sqln2/aD) * dyVoigt;     /* Lorentz factor    */      
      }
   }
   /* Derivatives for the background */
   epar[ncomp*maxpar+0] = 1.0;
   epar[ncomp*maxpar+1] = X_X0;
   epar[ncomp*maxpar+2] = X_X0*X_X0;   
}            





extern double funcd_c( double *xdat,
                       double *fpar,
                       fint   *npar,
                       fint   *fopt )
/*------------------------------------------------------------*/
/* PURPOSE: It depends on the value of the function option    */
/*          which function and how many components will be    */
/*          fitted. This function is called by 'lsqfitd' and  */
/*          returns the function value at certain X.          */
/*------------------------------------------------------------*/
{
   int      ncomp, fieid;
   double   result = 0.0;
   double   X = xdat[0];


   fieid = (*fopt) / FIEID;
   ncomp = *fopt - fieid*FIEID;
   
   if (fieid == GAUSS)
      result = gauss( X, fpar, ncomp );
   else if (fieid == GAUSSHERMITE1)
      result = gausshermiteh3( X, fpar, ncomp );
   else if (fieid == GAUSSHERMITE2)
      result = gausshermiteh3h4( X, fpar, ncomp );      
   else if (fieid == VOIGT)
      result = Voigt( X, fpar, ncomp );


   return( result );       
}




extern void dervd_c( double *xdat,
                     double *fpar,
                     double *epar,
                     fint   *npar,
                     fint   *fopt )
/*------------------------------------------------------------*/
/* PURPOSE: It depends on the value of the function option    */
/*          which function and how many components will be    */
/*          fitted. This function is called by 'lsqfitd' and  */
/*          returns the function derivatives at certain X.    */
/*------------------------------------------------------------*/
{
   int      ncomp, fieid;
   double   X = xdat[0];


   fieid = (*fopt) / FIEID;
   ncomp = *fopt - fieid*FIEID;
   
   if (fieid == GAUSS)
      gaussderv( X, fpar, epar, ncomp );
   else if (fieid == GAUSSHERMITE1)
      gausshermiteh3derv( X, fpar, epar, ncomp );
   else if (fieid == GAUSSHERMITE2)
      gausshermiteh3h4derv( X, fpar, epar, ncomp );      
   else if (fieid == VOIGT)
      Voigtderv( X, fpar, epar, ncomp );
      
}




int comp1( qsortstruct *s1, qsortstruct *s2 )
/*------------------------------------------------------------*/
/* PURPOSE: Compare function for 'qsort' only! Sort wrt.      */
/*          distance to peak.                                 */
/*------------------------------------------------------------*/
{
   if (s1->par[0] == dblank)     /* Sort dblanks to end of array */ 
      return( 1 ); 
   if (s2->par[0] == dblank)
      return( 1 ); 
   if (s1->par[0] == s2->par[0])
      return( 0 );
   if (fabs(s1->par[0]) > fabs(s2->par[0]))
      return( 1 );
   return( -1 );
}




int comp2( qsortstruct *s1, qsortstruct *s2 )
/*------------------------------------------------------------*/
/* PURPOSE: Compare function for 'qsort' only! Sort wrt.      */
/*          distance to the dispersion.                       */
/*------------------------------------------------------------*/
{
   if (s1->par[2] == dblank)     /* Sort dblanks to end of array */ 
      return( 1 ); 
   if (s2->par[2] == dblank)
      return( 1 ); 
   if (s1->par[2] == s2->par[2])
      return( 0 );
   if (fabs(s1->par[2]) > fabs(s2->par[2]))
      return( 1 );
   return( -1 );
}




int comp3( qsortstruct *s1, qsortstruct *s2 )
/*------------------------------------------------------------*/
/* PURPOSE: Compare function for 'qsort' only! Sort wrt.      */
/*          distance to Vsys.                                 */
/*------------------------------------------------------------*/
{
   if (s1->par[1] == dblank)     /* Sort dblanks to end of array */ 
      return( 1 ); 
   if (s2->par[1] == dblank)
      return( 1 ); 
   if (s1->par[1] == s2->par[1])
      return( 0 );
   if (fabs(s1->par[1]-s1->vsys) > fabs(s2->par[1]-s2->vsys))
      return( 1 );
   return( -1 );
}


 

static void mysort( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Sort the array's with (fitted) parameters.        */
/*                                                            */
/* Usually the sorting is about a few (two) components and    */
/* a fancy sorting algorithm is not necessary.                */
/*------------------------------------------------------------*/
{
   int          i, j;
   qsortstruct   comps[MAXGAUSS];

   
   for (i = 0; i < g->ncomp; i++)    /* Copy arrays to qsort struct */
   {
      for (j = 0; j < g->nvar; j++)
      {
         comps[i].par[j] = g->estimates[i*g->nvar+j];
         comps[i].err[j] = g->errorpar[i*g->nvar+j];         
      }   
   }

   comps[i].vsys = g->vsys;
   
   if      (g->sortopt == 1)                           /* peak */
      qsort( comps, g->ncomp, sizeof(qsortstruct), (int(*)())comp1 );
   else if (g->sortopt == 2)                           /* dispersion */
      qsort( comps, g->ncomp, sizeof(qsortstruct), (int(*)())comp2 );
   else if (g->sortopt == 3)                           /* Distance to vsys */
      qsort( comps, g->ncomp, sizeof(qsortstruct), (int(*)())comp3 ); 

   for (i = 0; i < g->ncomp; i++)       /* Copy sorted structs back to arrays */
   {
      for (j = 0; j < g->nvar; j++)
      {
         g->estimates[i*g->nvar+j] = comps[i].par[j];
         g->errorpar[i*g->nvar+j] = comps[i].err[j]; 
      }
   }
}




static int dolsqfit( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Do the real lsqfit.                               */
/* Return the number of iteration or a number < 0 for an      */
/* error.                                                     */
/*------------------------------------------------------------*/
{
   double       tolerance;
   double       mixingpar;
   fint         maxits;
   int          pnum = g->currentprofilenumber;
   fint         fiecode;
   int          i;
   fint         proflen = g->profilelen;
   fint         tpar;
   fint         xdim = 1;                       /* Dimension of fit */
   double       P0, Pmax;


   tolerance = MYMAX( g->tolerance, 0.0 );
   mixingpar = ABS( g->lab );
   maxits    = MYMAX( 1, g->maxits );


   P0   = (double) g->xphys[0];
   Pmax = (double) g->xphys[g->profilelen-1];

   for (i = 0; i < g->ncomp; i++) 
   {
      /* Transform center and factor a to interval [0,1] */
      g->estimates[i*g->nvar+1]  = (g->estimates[i*g->nvar+1] - P0)/(Pmax - P0);   
      g->estimates[i*g->nvar+2] /= (Pmax - P0); 
      if (g->function == VOIGT)
      {
         g->estimates[i*g->nvar+3] /= (Pmax - P0);
      }
   }
   g->estimates[g->ncomp*g->nvar+1] *= (Pmax - P0); 
   g->estimates[g->ncomp*g->nvar+2] *= (Pmax - P0)*(Pmax - P0);
  
   tpar = g->ncomp * g->nvar + 3;                /* Background pars added. */

   /*--------------------------------------------------*/
   /* Copy amplitudes in double array and weight data. */
   /* Blanks get weight == 0.0.                        */
   /*--------------------------------------------------*/   
   for (i = 0; i < proflen; i++) 
   { 
      g->wdat[i] = 1.0; 
      if (g->profiles[pnum][i] == dblank)
         g->wdat[i] = 0.0;
   } 
            
   /*anyoutf( DEBUG, "xdim=%d, proflen=%d, tpar=%d, tolerance=%g, maxits=%d, mixing=%g",
            xdim,proflen,tpar, tolerance, maxits, mixingpar );*/
            
   fiecode = FIEID*g->function + g->ncomp;

   g->iters = lsqfitd_c( 
                     g->xphysscale,          /* The scaled [0,1] X-coordinates */
                     &xdim,                  /* Dimension of fit */
                     g->profiles[pnum],      /* Y data */
                     g->wdat,                /* Weights */
                     &proflen,               /* Number of data points */
                     g->estimates,           /* Initial values */
                     g->errorpar,            /* Error return values */
                     g->fixedfree,           /* Fixed or free */
                     &tpar,                  /* Total number of parameters */
                     &tolerance,             /* Tolerance */ 
                     &maxits,                /* Max. num. of iterations */
                     &mixingpar,             /* Mixing parameter */
                     &fiecode );             /* Number of gaussians & fie. id */

   if (g->iters < 0)
   {
      int   offset = g->ncomp * 3;           /* Reserved for the estimates */
      for (i = offset; i < g->noutsubs; i++)
      {
         g->results[i][pnum] = dblank;
      }
      return( (int) g->iters );           
   }

   /*----------------------------------------------------------*/
   /* THE ARRAY G->ESTIMATES[] CONTAINS NOW FITTED PARAMETERS! */
   /*----------------------------------------------------------*/   
      
   /* Transform center and dispersion back to interval [physmin,physmax] */
   for (i = 0; i < g->ncomp; i++) 
   {
      g->estimates[i*g->nvar+1]  = g->estimates[i*g->nvar+1]*(Pmax - P0) + P0;
      g->estimates[i*g->nvar+2] *= (Pmax - P0); 
      g->errorpar[i*g->nvar+1]  *= (Pmax - P0);
      g->errorpar[i*g->nvar+2]  *= (Pmax - P0);
      if (g->function == VOIGT)
      {
         g->estimates[i*g->nvar+0] *= (Pmax - P0); 
         g->errorpar[i*g->nvar+0]  *= (Pmax - P0);          
         g->estimates[i*g->nvar+3] *= (Pmax - P0);
         g->errorpar[i*g->nvar+3]  *= (Pmax - P0);
      }         
   }
   /* Convert the background parameters B and A */
   g->estimates[g->ncomp*g->nvar+1] /= (Pmax - P0); 
   g->estimates[g->ncomp*g->nvar+2] /= (Pmax - P0)*(Pmax - P0);
   g->errorpar[g->ncomp*g->nvar+1]  /= (Pmax - P0); 
   g->errorpar[g->ncomp*g->nvar+2]  /= (Pmax - P0)*(Pmax - P0);

   /*--------------------------------------------------*/
   /* Here we have to sort the fitted parameters if    */
   /* there is more than one fitted component and the  */
   /* user wants sorting.                              */
   /*--------------------------------------------------*/
   if (g->sortopt)
   {
      mysort( g );
   }

   return( (int) g->iters ); 
}




void setuserguess( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Set user's estimates and set fit parameters to    */
/*          fixed or free.                                    */
/*------------------------------------------------------------*/
{
   int        i;
   fint       nitems = 1;
   fint       dfault = HIDDEN;
   fint       r;
   double     userval;
   char       key[STRLEN];
   bool       fixed;
   int        fie = g->function;    
   

   for (i = 0; i < g->ncomp; i++)    
   {
      int   indx = i*g->nvar;
      int   j;

      for (j = 0; j < g->nvar; j++)
      {                  
         sprintf( key, "%s_FIX%d=", prefix[fie][j], i );
         fixed = UserLog( key );
         if (fixed) 
            g->fixedfree[indx+j] = 0;
         else
            g->fixedfree[indx+j] = 1;

         sprintf( key, "%s_EST%d=", prefix[fie][j], i );         
                           
         r = userdble_c( &userval, &nitems, &dfault, tofchar(key), tofchar(" ") );
        
         /*-----------------------------------------------------*/       
         /* A field must be filled if it is empty and fixed     */
         /*-----------------------------------------------------*/
         if (!r)
         {
            g->fixedvalues[indx+j] = 0.0;                   
            wkeyf( "%s%g", key, g->fixedvalues[indx+j] );            
         }
         else
         {
            g->fixedvalues[indx+j] = userval;
         }
      }         
   }

   /* Background components */
   for (i = 0; i < 3; i++) 
   {
      int  indx = g->ncomp * g->nvar + i;
      
      sprintf( key, "%s_FIX=", prefix[MAXFUNCTIONS][i] );
      fixed = UserLog( key );
      if (fixed) 
         g->fixedfree[indx] = 0;
      else
         g->fixedfree[indx] = 1;
         
      sprintf( key, "%s_EST=", prefix[MAXFUNCTIONS][i] );
      r = userdble_c( &userval, &nitems, &dfault, tofchar(key), tofchar(" ") );
      /*-----------------------------------------------------*/       
      /* A field must be filled if it is empty and fixed     */
      /*-----------------------------------------------------*/
      if (!r)
      {
         g->fixedvalues[indx] = 0.0;                
         wkeyf( "%s%g", key, g->fixedvalues[indx] );
      }
      else
      {
         /* Always use this value, fixed or not fixed. */
         g->fixedvalues[indx] = userval;
      }
   }      
}



static double getresidual( globaltype *g, 
                           double     *Destimates,
                           double      zerolev )
/*------------------------------------------------------------*/
/* PURPOSE: Routine to determine a chi2 for the estimates     */
/*          function. This function is called in a loop over  */
/*          Q, the smoothing parameter.                       */
/*------------------------------------------------------------*/
{
   fint    nvar = 3;            /* It is always a standard gauss -> nvar=3*/
   fint    npar, opt;
   int     i;
   int     profnum = g->currentprofilenumber;
   double  chi2 = 0.0;

   
   /* Note that the number of gaussians found in the estimate */
   /* routine never exceeds the maximum number of components  */
   /* entered by the user */
   
   npar = g->gaussiansfound * nvar + 3;
   
   /* The estimate function is always a standard gauss */
   opt  = FIEID*GAUSS + g->gaussiansfound; 
   
   /* Add the parameters for the background. */
   Destimates[g->gaussiansfound * nvar+0] = zerolev;
   Destimates[g->gaussiansfound * nvar+1] = 0.0;
   Destimates[g->gaussiansfound * nvar+2] = 0.0;   

   
   for (i = 0; i < g->profilelen; i++)
   {
      double   delta;
      double   Yi = g->profiles[profnum][i];
      double   Xi = (double) i;    /*g->xphys[i];*/
    
      if (Yi != dblank)
      {
         delta = Yi - funcd_c( &Xi, Destimates, &npar, &opt );
         chi2 += delta * delta;
      }
   }
   return( chi2 );
}




static void fillfixedvalues( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Fill the estimates that are fixed in the fit.     */
/*------------------------------------------------------------*/
{
   int     i;
   
   for (i = 0; i < g->ncomp*g->nvar+3; i++)
   {
      if (g->fixedfree[i] == 0)               /* i.e. fixed ! */
         g->estimates[i] = g->fixedvalues[i];
   }
}



static int getestimatesfromset( globaltype *g ) 
/*------------------------------------------------------------*/ 
/* PURPOSE: Get the estmates from a set.                      */
/* If the first value is a blank, return immediately.         */
/*------------------------------------------------------------*/
{
   int   i;
   int   pnum = g->currentprofilenumber; 
   
   if (g->estfromset[pnum][0] == dblank)
      return( NO );
      
   for (i = 0; i < 3*g->ncomp; i++)
   {
      g->estimates[i] = g->estfromset[pnum][i];
   }
   
   /* Store these estimates in the results array */
   for (i = 0; i < g->ncomp; i++)  
   {
      g->results[i*3+0][pnum] = g->estimates[i*3+0];  /* Amplitude */
      g->results[i*3+1][pnum] = g->estimates[i*3+1];  /* Phys center */
      g->results[i*3+2][pnum] = g->estimates[i*3+2];  /* Phys dispersion */
   }  
 
   /* Background parameters */
   for (i = 0; i < 3; i++)
   {
      int   indx = g->ncomp*g->nvar+i;
      g->estimates[indx] = g->fixedvalues[indx];
   }
  
   return( YES );
}




static int getestimates( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Prepare for lsqfit.                               */
/*------------------------------------------------------------*/
{
   fint         maxgaussians = g->ncomp;
   double       rms = g->rms;
   double       critamp = g->critamp;
   double       critdisp = g->critdisp;      /* In physical coordinates */
   double       zerolev;
   int          pnum = g->currentprofilenumber;
   fint         proflen = g->profilelen;
   double       Destimates[MAXPARALL];

   

   /*--------------------------------------------------*/ 
   /* We need initial values for the background,       */
   /* with this value blanks are replaced for the      */
   /* estimate routine.                                */
   /*--------------------------------------------------*/ 
   {
      int     i;      
      for (i = 0; i < 3; i++)
      {
         int   indx = g->ncomp*g->nvar+i;
         g->estimates[indx] = g->fixedvalues[indx];
      }
   }
      
   zerolev = g->estimates[g->ncomp*g->nvar+0];
   g->rms = ABS( g->rms );   
   rms = g->rms;      
   g->usedQ = fblank;
   /*--------------------------------------------------*/
   /* Next, check the other parameters. The critical   */
   /* dispersion is entered in physical coordinates    */
   /* and must be transformed to grids.                */      
   /*--------------------------------------------------*/
   critamp      = ABS( critamp );
   maxgaussians = BETWEEN( maxgaussians, 1, MAXGAUSS );

   critdisp    /= ABS(g->gridspac);                    /* Now it is in grids */
                      

   /* Copy profiles[pnum] into filtered */
   {
      int     i;
         
      for (i = 0; i < proflen; i++)
      {         
         if (g->profiles[pnum][i] == dblank)
         {
            g->filtered[i] = dblank;
         }
         else
         {
            g->filtered[i] = g->profiles[pnum][i];
         }
      }
   }

   g->gaussiansfound = 0;


   /*--------------------------------------------------*/ 
   /* The Ulrich Schwarz method is used to search for  */
   /* gaussian components in a profile:                */
   /* The routine uses first an automatic window method*/
   /* to define the signal region. Then the second     */
   /* derivative of the profile in the signal region   */
   /* is calculated by fitting a second degree         */
   /* polynomal. The smoothing parameter Q determines  */
   /* the number of points used for this (=2*Q+1).     */
   /*--------------------------------------------------*/ 
   {
      int     i;
      fint    maxparameters = 3 * maxgaussians;  /* Ampl, center,disp only */
      double  residual;
      int     k, minindx = 0;
      int     specialindx = -1;
      int     special = NO;
      int     signal = 0; 
      
     
      /*------------------------------*/
      /* Subtract zero level first    */
      /*------------------------------*/      
      for (i = 0; i < proflen; i++)              
      {
         if (g->filtered[i] != dblank)
            g->filtered[i] -= zerolev;
         else
            g->filtered[i] = zerolev;
         if (g->filtered[i] != 0.0)
         {
            specialindx = i;
            signal++;
         }
      }
      if (signal == 0)
      {
         /* There is no hope for this profile */
         for (i = 0; i < 3*g->ncomp; i++)
         {
            g->results[i][pnum] = dblank;
         }
         return( 0 );
      } 
      if (signal == 1)
      {
         int     foundQis1 = NO;
         /*--------------------------------------------------*/
         /* This is a special case where all points are zero */
         /* except one. If this is a positive amplitude then */
         /* it is the maximum of a very small peak for which */
         /* gauest will not find initial estimates if Q == 1.*/
         /* Therefore we have to search for Q=1 in the array */
         /* of smoothing factors and if there is one, then   */
         /* it is automatically the best initial estimate Q  */
         /* and we can construct special estimates using     */
         /* FWHM == 1 pixel which is equivalent to a         */
         /* dispersion of 1/2.35 pixels.                     */
         /*--------------------------------------------------*/
         for (k = 0; k < g->numq; k++)
         {
            if (g->smoothingpar[k] == 1)
            {
               foundQis1 = YES;
               minindx = k;
               break;
            }
         }
         if (foundQis1)
         {
            Destimates[0] = g->filtered[specialindx];
            Destimates[1] = specialindx;
            Destimates[2] = 1.0 / 2.35482;   /* FWHM --> dispersion */
            if (Destimates[0] < critamp || Destimates[2] < critdisp)
            {
               for (i = 0; i < 3*g->ncomp; i++)
               {
                  g->results[i][pnum] = dblank;
               }
               return( 0 );
            }
            else 
            {
               g->gaussiansfound = 1;
               g->usedQ = 1;
               special = YES;
            }
         }
      }      
      
      if (g->numq > 1 && !special) 
      {
         double  minres = dblank;        /* Initialize minimum residual */
         minindx = -1;                   /* Flag for not finding a Q that works */
         /*----------------------------------------*/
         /* We must try all Q factors.             */
         /*----------------------------------------*/
         for (k = 0; k < g->numq; k++)
         {
            fint Q = g->smoothingpar[k];
            if (2*Q+1 <= proflen)
            {
               double  zeroamp = 0.0, zerodisp = 0.0;
               g->gaussiansfound = gauestd_c( 
                   g->filtered,          /* Instead of  g->profiles[pnum] */
                   g->work,              /* Used here as a real work array */
                   &proflen,
                   Destimates,           /* Output parameter estimates */
                   &maxparameters,       /* 3 * ncomps */
                   &rms,
                   &zeroamp,
                   &zerodisp,
                   &Q );
               g->gaussiansfound = MYMIN( g->gaussiansfound, g->ncomp );
               if (g->gaussiansfound > 0)
               {
                  residual = getresidual( g, Destimates, zerolev );
                  /* anyoutf( DEBUG, "Residual for q[%d] = %d: %f", k, g->smoothingpar[k], residual );*/
                  if (minres == dblank)
                  {
                     minres = residual;
                     minindx = k;
                  }
                  else
                  {
                     if (residual < minres)
                     {
                        minres = residual;
                        minindx = k;
                     }
                  }
               }
            } /* End if proflen */
         } /* end loop over all Q's */
      } /* end if num Q's > 1 */
      
      /*--------------------------------------------------*/
      /* Now we have either the only Q in an array of 1   */
      /* smoothingfactor, or we have the 'best' Q in an   */
      /* array of 'numq' smoothingfactors. Then it could  */
      /* occur that the program did not find a valid Q    */
      /* which is set by minindx=-1 and there is no need  */
      /* to continue. For other indices we have to repeat */
      /* a call to 'gauest' to find the initial estimates */
      /* corresponding to the current value of Q.         */
      /*--------------------------------------------------*/
      if (minindx == -1)
      {
         /* There is no hope for this profile */
         for (i = 0; i < 3*g->ncomp; i++)
         {
            g->results[i][pnum] = dblank;
         }
         return( 0 );
      }  


      if (!special)   /* For 'special' we already have estimates */
      {
         /* Now 'minindx' is in range 0..numq-1 */   

          g->gaussiansfound = gauestd_c( 
                   g->filtered,          /* Instead of  g->profiles[pnum] */
                   g->work,              /* Used here as a real work array */
                   &proflen,
                   Destimates,           /* Output parameter estimates */
                   &maxparameters,       /* 3 * ncomps */
                   &rms,
                   &critamp,
                   &critdisp,
                   &(g->smoothingpar[minindx]) );
         g->gaussiansfound = MYMIN( g->gaussiansfound, g->ncomp );

         /* Correct number if more were found than wanted */
         if (g->gaussiansfound > maxgaussians)
            g->gaussiansfound = maxgaussians;

         if (g->gaussiansfound == 0)
         {
            /* If there is only one Q and no estimates were found: */ 
            for (i = 0; i < 3*g->ncomp; i++)
            {
               g->results[i][pnum] = dblank;
            }
            return( 0 );
         }
         g->usedQ = (float) g->smoothingpar[minindx];
      }

      /* Copy the  estimates */         
      for (i = 0; i < g->gaussiansfound; i++)
      {
         g->estimates[i*g->nvar+0] = Destimates[i*3+0];
         g->estimates[i*g->nvar+1] = Destimates[i*3+1];
         g->estimates[i*g->nvar+2] = Destimates[i*3+2];
      }
      /* Fill with zero's if less gaussians were found */
      for (i = g->gaussiansfound; i < g->ncomp; i++)
      {
         g->estimates[i*g->nvar+0] = 0.0;
         g->estimates[i*g->nvar+1] = 0.0;
         g->estimates[i*g->nvar+2] = 0.0;            
      }
   }
   
   
   /*--------------------------------------------------*/
   /* Transform x-axis related values from interval    */
   /* [0,profilelen] to interval                       */
   /* [physical_min, physical_max]                     */
   /* Don't do this for the estimates that were set to */
   /* 0.                                               */
   /*--------------------------------------------------*/
   {
      int k;
      for (k = 0; k < g->gaussiansfound; k++)
      {
         g->estimates[k*g->nvar+1] = g->xphys[0] + g->estimates[k*g->nvar+1] * g->gridspac;
         g->estimates[k*g->nvar+2] = g->estimates[k*g->nvar+2] * ABS(g->gridspac);
      }
   }

   /* Store these estimates in the results array */
   {
      int i;
      for (i = 0; i < g->ncomp; i++)  
      {
         g->results[i*3+0][pnum] = g->estimates[i*3+0];  /* Amplitude */
         g->results[i*3+1][pnum] = g->estimates[i*3+1];  /* Phys center */
         g->results[i*3+2][pnum] = g->estimates[i*3+2];  /* Phys dispersion */
      }
   }
   return( 1 );
}



static void adjustestimates( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Adjust the guesses into the estimates array.      */
/*------------------------------------------------------------*/
{
   if (g->function == GAUSS)
   {
      return;
   }
   else
   {
      int i;
      for (i = 0; i < g->gaussiansfound; i++)
      {
         int   offset = i*g->nvar;
         if (g->function == GAUSSHERMITE1)
            g->estimates[offset+3] = g->h3;
         if (g->function == GAUSSHERMITE2)
         {
            g->estimates[offset+3] = g->h3;
            g->estimates[offset+4] = g->h4;
         } 
         if (g->function == VOIGT)
         {
            /* First estimate is area instead of amplitude */
            double amp   = g->estimates[offset+0];
            double sigma = g->estimates[offset+2];
            g->estimates[offset+0] = amp * sigma * sqrt( 2.0*PI ); /* A->Area */
            disp2fwhm( &sigma );                                      /* FWHM */
            g->estimates[offset+2] = 0.5 * sigma;                     /* HWHM */
            g->estimates[offset+3] = g->lorentz;
         }
      }        
   }
}
  
 


void dofit( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Fit current profile and put results in results    */
/*          matrix.                                           */
/*------------------------------------------------------------*/ 
{
   int   r;
      
   /*----------------------------------------*/
   /* Estimates came from a set.             */
   /*----------------------------------------*/   
   if (g->altset)
   {
      if (!getestimatesfromset(g))              /* No fit. Blank output */
      {
         int   offset = g->ncomp * 3;           /* Reserved for the estimates */
         int   i;
         for (i = offset; i < g->noutsubs; i++)
         {
            g->results[i][g->currentprofilenumber] = dblank;
         }
         g->log.noestimates++;
         r = -100;
      }
      else if (!g->estonly)                     /* Try a fit */
      {
         adjustestimates( g );         
         r = dolsqfit( g ); 
      }
   }
   /*----------------------------------------*/
   /* Program must calculates estimates.     */
   /*----------------------------------------*/
   else
   {
      if (!getestimates(g))
      {
         g->log.noestimates++;
      }
      if (!g->estonly)
      {
         adjustestimates( g );      
         fillfixedvalues( g );      /* This is what the user wants to fix */
         r = dolsqfit( g );         /* No estimates? Try a fit anyhow    */
                                    /* (with the previous estimates) */   
      }
   }      
   
   /*-----------------------------------------*/
   /* Besides estimates, user wants fits also */
   /*-----------------------------------------*/
   if (!g->estonly)
   {
      if (r < 0)
      {
         if (r == -4)
            g->log.iterstoosmall++;
         g->log.nofit++;    
      }
      else
      /*----------------------------------------*/
      /* Write the results to the subsets.      */
      /*----------------------------------------*/   
      {
          outputresults( g );
      }
   }
}
#<

#>            utils.c
/* utils.c
                           COPYRIGHT (c) 1999
                      Kapteyn Astronomical Institute
                University of Groningen, The Netherlands
                           All Rights Reserved.
*/

#include    "stdio.h"        /* Defines ANSI C input and output utilities */
#include    "stdlib.h"       /* Defines the ANSI C functions for number */
                             /* conversion, storage allocation, and similar tasks.*/
#include    "string.h"       /* Declares the ANSI C string functions*/
                             /* like:strcpy, strcat etc.*/
#include    "math.h"         /* Declares the mathematical functions and macros.*/
#include    "gipsyc.h"       /* Defines the ANSI-F77 types for Fortran to C intface */
                             /* including def. of char2str,str2char,tofchar,zadd */
                             /* and macros tobool and toflog */
#include    "float.h"        /* Definition of FLT_MAX etc.*/
#include    "ctype.h"        /* Declares ANSI C functions for testing characters */
                             /* like: isalpha, isdigit etc. also tolower, toupper.*/


/* Common includes */

#include    "anyout.h"       /* General character output routine for GIPSY programs.*/
#include    "error.h"        /* User error handling routine. */
#include    "nelc.h"         /* Characters in F-string discarding trailing blanks.*/
#include    "minmax1.h"
#include    "axunit.h"
#include    "cotrans.h"
#include    "status.h"
#include    "gui.h"
#include    "userlog.h"
#include    "events.h"
#include    "userdble.h"
#include    "setdblank.h"

#define FIRSTLINE      0
#define SECONDLINE     1    
#define MAXAXES        10
#define HIDDEN         2


extern  double   dblank;
extern  float    fblank;



void clearstr( fchar Fstr )
/*------------------------------------------------------------*/
/* Purpose: Blank a Fortran string up to 'len' characters     */
/*------------------------------------------------------------*/
{
   int    i;
   fint   len = Fstr.l;


   for (i = 0; i < (int) len; i++)
   {
      Fstr.a[i] = '\0';
   }
} 




void hermesstatus( char *txt )
/*------------------------------------------------------------*/
/* PURPOSE: Update status in hermes status                    */
/*------------------------------------------------------------*/
{
   status_c( tofchar(txt) );
}



void error_status( char *errstr,
                   char *statstr1,
                   char *statstr2 )
/*------------------------------------------------------------*/
/* PURPOSE: Update status and error messages in gui.          */
/*------------------------------------------------------------*/
{
   gui_error( errstr );
   gui_status( statstr1, FIRSTLINE );
   hermesstatus( statstr1 );
   gui_status( statstr2, SECONDLINE );
}



void clearstatus( void )
/*------------------------------------------------------------*/
/* PURPOSE: Clear status and error messages in gui.           */
/*------------------------------------------------------------*/
{
   gui_error( "" );
   gui_status( "", FIRSTLINE );
   gui_status( "", SECONDLINE );
   hermesstatus("");
}  




void dminmax( double *data, fint *n, double *amin, double *amax )
/*------------------------------------------------------------*/
/* PURPOSE: Min and max in double prec. of array data.        */
/*------------------------------------------------------------*/
{
   fint   def = 0;                      /* define mode or not */
   fint   l;                            /* loop counter */
   double blank;                        /* BLANK value */

   setdblank_c( amin );                 /* set to blank, as initial value */
   setdblank_c( amax );                 /* set to blank, as initial value */
   setdblank_c( &blank );               /* set to blank for comparison */
   for (l = 0; l < *n; l++, data++) {   /* loop */
      if (*data != blank) {             /* no BLANK value */
         if (def) {                     /* min/max already defined */
            if (*amin > *data) {        /* new minimum ? */
               *amin = *data;           /* save new minimum */
            } else if (*amax < *data) { /* new maximum ? */
               *amax = *data;           /* save new maximum */
            }
         } else {                       /* define min and max */
            *amin = *amax = *data;      /* save new minimum and maximum */
            def = 1;                    /* turn on define mode */
         }
      }
   }
}    




int grid2phys2( fchar  Setin, 
                fint   grid,
                fint   cowo, 
                int    axnum,
                double *physval )
/*------------------------------------------------------------*/
/* PURPOSE: Return physical value of a grid along profile.    */
/*          Return 0 if no conversion was possible.           */ 
/*------------------------------------------------------------*/
{
   fint    r;
   double  coordout[MAXAXES];        /* Phys. coords after cotrans */ 
   int     cindx;
   fint    grid2phys = 1; 
   double  gridin = (double) grid;
  
 
   /* Axis numbers are 1, 2, ... */
   cindx = axnum - 1;

   r = cotrans_c( Setin, &cowo, &gridin, coordout, &grid2phys );
   if (r != 0)   
      return( 0 );
   else
      *physval = coordout[cindx];
   return( 1 );
}



bool UserLog(char *key)
/*------------------------------------------------------------*/
/* Purpose: Special version of userlog with default.          */
/*------------------------------------------------------------*/
{
   bool result = toflog(FALSE);
   fint fint1 = 1, fint2 = 2;

   (void) userlog_c( &result, &fint1, &fint2, tofchar(key), tofchar(" ") );
   return( tobool(result) );
}




void scalearray( double  *in,
                 double  *out,
                 int     len )
/*------------------------------------------------------------*/
/* PURPOSE: Convert array to double and convert values        */
/*          between 0 and 1.                                  */
/*------------------------------------------------------------*/
{
   int      i;
   double   P0, Pmax, Pi;
            
   P0   = in[0];
   Pmax = in[len-1];

   /* Scale values so that they fit between 0 and 1 */
   for (i = 0; i < len; i++) 
   {
      Pi = (double) in[i];
      out[i] = (Pi - P0)/(Pmax - P0);
   }
}



double **dmatrix( int xlo,
                  int ylo,
                  int xhi,
                  int yhi )
/*------------------------------------------------------------------*/
/* PURPOSE: Allocate space for a matrix with subscript range        */
/*          M[ylo..yhi][xlo..xhi]                                   */
/*------------------------------------------------------------------*/
{
   int      rows = yhi - ylo + 1;
   int      cols = xhi - xlo + 1;
   int      i;
   double   **m;


   /* Allocate memory for pointers to rows */
   m = (double **) malloc( rows * sizeof(double*) );
   if (!m)
      return( NULL );
   m -= ylo;                                /* Adjust subscript */

   /* Pointer to first row allocates memory for box */
   m[ylo] = (double *) malloc( rows * cols * sizeof(double) );
   if (!m[ylo])
       return( NULL );
   m[ylo] -= xlo;                           /* Adjust subscript */

   /* Set pointers to rows */
   for (i = ylo+1; i <= yhi; i++)
      m[i] = m[i-1] + cols;

   /* Return pointer to array of pointers to rows */
   return( m );
}



void freedmatrix( double **m,
                  int   xlo,
                  int   ylo )
/*------------------------------------------------------------------*/
/* PURPOSE: Free space allocated by a matrix made with function     */
/*          'fmatrix'.                                              */
/*------------------------------------------------------------------*/
{
   free( m[ylo] + xlo );
   free( m + ylo );
   m = NULL;
}


float  **fmatrix( int xlo,
                  int ylo,
                  int xhi,
                  int yhi )
/*------------------------------------------------------------------*/
/* PURPOSE: Allocate space for a matrix with subscript range        */
/*          M[ylo..yhi][xlo..xhi]                                   */
/*------------------------------------------------------------------*/
{
   int      rows = yhi - ylo + 1;
   int      cols = xhi - xlo + 1;
   int      i;
   float    **m;


   /* Allocate memory for pointers to rows */
   m = (float **) malloc( rows * sizeof(float*) );
   if (!m)
      return( NULL );
   m -= ylo;                                /* Adjust subscript */

   /* Pointer to first row allocates memory for box */
   m[ylo] = (float *) malloc( rows * cols * sizeof(float) );
   if (!m[ylo])
       return( NULL );
   m[ylo] -= xlo;                           /* Adjust subscript */

   /* Set pointers to rows */
   for (i = ylo+1; i <= yhi; i++)
      m[i] = m[i-1] + cols;

   /* Return pointer to array of pointers to rows */
   return( m );
}



void freefmatrix( float **m,
                  int   xlo,
                  int   ylo )
/*------------------------------------------------------------------*/
/* PURPOSE: Free space allocated by a matrix made with function     */
/*          'fmatrix'.                                              */
/*------------------------------------------------------------------*/
{
   free( m[ylo] + xlo );
   free( m + ylo );
   m = NULL;
}




void getdble( ident id,
                     char  *key,
                     int   code,
                     void  *arg )
/*------------------------------------------------------------*/
/* PURPOSE: Keyword handler for doubles.                      */
/*------------------------------------------------------------*/
{
    double     *userval = (double *) arg;
    fint       nitems = 1;
    fint       dfault = HIDDEN;

    (void) userdble_c( userval,
                       &nitems,
                       &dfault,
                       tofchar(key),
                       tofchar(" ") );
}



float tofloat( double in )
/*------------------------------------------------------------*/
/* PURPOSE: Convert double to float and return blank when     */
/*          outside scope.                                    */
/*------------------------------------------------------------*/ 
{
    if (in == dblank)
       return( fblank );
    if (in  > FLT_MAX || in < -FLT_MAX)
       return( fblank );
    return( (float) in );
}

#<

#>            allocs.c
/* allocs.c
                            COPYRIGHT (c) 1999
                      Kapteyn Astronomical Institute
                University of Groningen, The Netherlands
                           All Rights Reserved.
*/


#include    "stdio.h"
#include    "stdlib.h"       /* Defines the ANSI C functions for number */
                             /* conversion, storage allocation, and similar tasks.*/
#include    "string.h"       /* Declares the ANSI C string functions*/
                             /* like:strcpy, strcat etc.*/

#include    "xgaufit.h"   



int allocs( globaltype *g,
            int         len )
/*------------------------------------------------------------*/ 
/* PURPOSE: Allocate memory for profile arrays.               */ 
/*------------------------------------------------------------*/
{
   /*----------------------------------------*/
   /* Profile coordinates in subset indices. */
   /*----------------------------------------*/
   if (g->xgrids != NULL)    
      free( g->xgrids );   

   g->xgrids = (double *) calloc(len, sizeof(double) );
   if (g->xgrids == NULL)
   {
      return( 0 );
   }


   /*----------------------------------------*/
   /* Profile coordinates in physical        */
   /* coordinates.                           */
   /*----------------------------------------*/
   if (g->xphys != NULL)    
      free( g->xphys );

   g->xphys = (double *) calloc(len, sizeof(double) );
   if (g->xphys == NULL)
   {
      return( 0 );
   }

   
   /*----------------------------------------*/
   /* Work array of profile length.          */
   /*----------------------------------------*/   
   if (g->work != NULL)
      free( g->work ); 
   
   g->work = (double *) calloc(len, sizeof(double) );
   if (g->work == NULL)
   {
     return( 0 );
   }   


   /*----------------------------------------*/
   /* Median filtered version of prof. amps. */
   /*----------------------------------------*/
   if (g->filtered != NULL)
      free( g->filtered ); 
   
   g->filtered = (double *) calloc(len, sizeof(double) );
   if (g->filtered == NULL)
   {
     return( 0 );
   }   


   /*----------------------------------------*/
   /* Resampled profile (Spline).            */
   /*----------------------------------------*/
   if (g->spline != NULL)
      free( g->spline ); 
   
   g->spline = (double *) calloc(len, sizeof(double) );
   if (g->spline == NULL)
   {
     return( 0 );
   }   


   /*----------------------------------------*/ 
   /* Double precision weights array for     */
   /* lsqfit function.                       */
   /*----------------------------------------*/ 
   if (g->wdat != NULL)
      free( g->wdat ); 
   
   g->wdat = (double *) calloc(len, sizeof(double) );
   if (g->wdat == NULL)
   {
     return( 0 );
   }   


   /*----------------------------------------*/
   /* Array to store scaled physical X-coord-*/
   /* inates. The scaled interval is -1<x<1  */
   /*----------------------------------------*/      
   if (g->xphysscale != NULL)    
      free( g->xphysscale );

   g->xphysscale = (double *) calloc(len, sizeof(double) );
   if (g->xphysscale == NULL)
   {
      return( 0 );
   }


   /*----------------------------------------*/
   /* Array to store median filtered profile.*/
   /*----------------------------------------*/      
   if (g->median != NULL)    
      free( g->median );

   g->median = (double *) calloc(len, sizeof(double) );
   if (g->median == NULL)
   {
      return( 0 );
   }


   /*----------------------------------------*/
   /* Array to store residuals after         */
   /* the fit.                               */
   /*----------------------------------------*/      
   if (g->Yres != NULL)    
      free( g->Yres );

   g->Yres = (double *) calloc(len, sizeof(double) );
   if (g->Yres == NULL)
   {
      return( 0 );
   }


   /*----------------------------------------*/
   /* Array to store fitted cube data after  */
   /* the fit.                               */
   /*----------------------------------------*/      
   if (g->Yfit != NULL)    
      free( g->Yfit );

   g->Yfit = (float *) calloc(len, sizeof(float) );
   if (g->Yfit == NULL)
   {
      return( 0 );
   }


   return( 1 );
} 
#<

#>            fitopt.c
/* fitopt.c
                              COPYRIGHT (c) 1999
                        Kapteyn Astronomical Institute
         University of Groningen  -  9700 AV Groningen, The Netherlands


Author: J.P. Terlouw
*/

#define NAME "Miscellaneous fit options"

#define H_CRTA "Critical amplitude of gaussian. Estimates\n\
below this amplitude will be discarded."

#define H_CRTD "Critical dispersion (in physical coordinates) of gaussian.\n\
Smaller estimates will be discarded."

#define H_Q "Number of data points used for second order polynomial in \n\
'gauest' routine is 2Q + 1. The smoothing factor Q with the lowest chi^2 \n\
produces the final estimates."

#define H_INCL "Only data within include range will\n\
be used for fit (default: use all data)"

#define H_RMS "R.m.s. of noise. Used only to\nscale reduced chi-squared."

#include "stddef.h"
#include "stdio.h"
#include "gipsyc.h"
#include "userfio.h"
#include "wkey.h"
#include "ggi.h"
#include "fitopt.h"

#define New(type) ((type *)calloc(1,sizeof(type)))
#define NNew(n,type) ((type *)calloc(n,sizeof(type)))
#define Delete(x)  {free(x); x=NULL;}

static bool  initialized=FALSE;         /* initialized flag */
static ident s_fitopt;                  /* fitopt keywords shell     */
static ident i_box;

/* -------------------------------------------------------------------------- */
/*                                 initialize                                 */
/* -------------------------------------------------------------------------- */
/*  Create pop-up shell and its contents; schedule the appropriate handlers.
 */
static void initialize(char *key)
{
   ident closebut;
   ident l_lsq, l_est;
   ident tol, lab, mxit;
   ident rms;
   ident crta, crtd, smf;
   ident l_ef, estf;

   s_fitopt = GgiShell(NAME);
   (void)GgiUseShell(s_fitopt);
   closebut = GgiSetLabel(GgiButton(key, "Close window"), "CLOSE", 0);
   i_box = GgiForm(NAME, 2);
   (void)GgiUseShell(i_box);
   l_lsq = GgiLabel("Least-squares fit options:");
   tol  = GgiSetLabel(GgiTextField("TOLERANCE=", "Convergence criterion", 9),
                      "Tolerance",85);
   lab  = GgiSetLabel(GgiTextField("LAB=", "Mixing parameter", 9),
                      "Lambda", 0);
   mxit = GgiSetLabel(GgiTextField("MAXITS=",
                                   "Max. number of iterations in lsq. fit", 5),
                      "Max. Its.",0);
   rms  = GgiSetLabel(GgiTextField("RMS=", H_RMS, 9), "RMS", 85);
      
   l_est = GgiLabel("Initial estimates options:");
   crta =  GgiSetLabel(GgiTextField("CRITAMP=", H_CRTA, 9), "Min. ampl.", 85);
   crtd =  GgiSetLabel(GgiTextField("CRITDISP=", H_CRTD, 9), "Min. disp.", 70);
   smf  =  GgiSetLabel(GgiTextField("Q=", H_Q, 9), "Q", 0);
   l_ef =  GgiSetLabel(GgiLabel(" "), "Write estimates only:", 0);
   estf =  GgiSetLabel(GgiButton("ESTONLY=", "Write only estimates to output"), "  ",0);

   GgiSetPosition(l_lsq,  0, NULL,  0, NULL);
   GgiSetPosition(tol,    0, NULL,  0, l_lsq);
   GgiSetPosition(lab,    0, tol,   0, l_lsq);
   GgiSetPosition(mxit,   0, lab,   0, l_lsq);
   GgiSetPosition(rms,    0, NULL,  0, tol);
   GgiSetPosition(l_est,  0, NULL, 10, rms);
   GgiSetPosition(crta,   0, NULL,  0, l_est);
   GgiSetPosition(crtd,   0, crta,  0, l_est);
   GgiSetPosition(smf ,   0, crtd,  0, l_est);
   GgiSetPosition(l_ef,   0, NULL,  0, crta);
   GgiSetPosition(estf,   0, l_ef,  0, crta);

   (void)GgiUseShell(s_fitopt);
   GgiSetPosition(i_box,    0, NULL, 0, NULL);
   GgiSetPosition(closebut, 0, NULL, 0, i_box);
   (void)GgiUseShell(NULL);
   initialized = TRUE;
}

/* -------------------------------------------------------------------------- */
/*                                 fitopt                                     */
/* -------------------------------------------------------------------------- */
/*  Keyword handler causing the shell to be popped up or down.
 */
static void fitopt(ident id, char *key, int code, void *arg)
{
   static bool state=FALSE;
   bool button=FALSE;

   (void)userflog(&button, 1, 2, key, " ");
   if (tobool(button)) {
      GgiShowShell(s_fitopt, !state);
      state = !state;
      wkey_c(tofchar(key));
   }
}

/* ========================================================================== */
/*                                 Fitopt                                     */
/* -------------------------------------------------------------------------- */
/*  Fitopt() creates a pop-up shell for fitopt keywords. The argument 'key'
 *  is the GIPSY keyword (type LOGICAL) which, when set to "YES", will
 *  cause the shell to be popped up or down.
 *  The same keyword is used for the <CLOSE> button in the shell.
 */
extern ident Fitopt(char *key)
{
   initialize(key);
   ScheduleKeyevent(fitopt, key, KEYCHANGE, NULL); 
   return s_fitopt;
}

#<

#>            import.c
/* import.c
                              COPYRIGHT (c) 1999
                        Kapteyn Astronomical Institute
         University of Groningen  -  9700 AV Groningen, The Netherlands

Author: J.P. Terlouw
*/

#define NAME "Import"

#define H_FILE "Text file name"
#define H_ROWS "Row numbers:\n\n first : last\n first :\n : last\n   :  (=all rows)"
#define H_FXC  "X column number"
#define H_FYC  "Y column number"
#define H_FEC  "Errors column number"
#define H_SET  "Set or subset containing table"
#define H_TAB  "Table name"
#define H_SXC  "X column name"
#define H_SYC  "Y column name"
#define H_SEC  "Errors column name"

#include "stddef.h"
#include "stdio.h"
#include "gipsyc.h"
#include "string.h"
#include "userfio.h"
#include "ggi.h"
#include "import.h"

static char *datanames[]={"PROFX", "PROFY", "PROFYERR", NULL};

/* -------------------------------------------------------------------------- */
/*                                 filedata                                   */
/* -------------------------------------------------------------------------- */
/*  Keyword handler for IMP_FILE=, IMP_FROWS= and xxxCOL=.
 *  Convert file, row and column information to 'standard' keywords.
 */
static void filedata(ident id, char *key, int code, void *arg)
{
   static bool initialized=FALSE;
   static char file[64]="";
   static char rows[64]="";
   fchar  data;
   int    nc;
   char   **name;
   char   *dataname=(char*)arg;
   char keybuf[KEYLEN+1];
   char col[16+1];

   data.l = 63;
   if (!initialized) {
      data.a = file;
      nc = userftext(data, 2, "IMP_FILE=", " ");
      data.a[nc] = '\0';
      data.a = rows;
      nc = userftext(data, 2, "IMP_FROWS=", " ");
      data.a[nc] = '\0';
      initialized = TRUE;
   }
   data.a = NULL;
   if (!strcmp(key,"IMP_FILE=")) {
      data.a = file;
   } else if (!strcmp(key,"IMP_FROWS=")) {
      data.a = rows;
   }
   if (data.a) {
                                        /* IMP_FILE= or IMP_FROWS= keyword seen */
      nc = userftext(data, 2, key, " ");
      data.a[nc] = '\0';
      if (!nc && data.a==rows) strcpy(rows,":");
      if (*file && *rows) {
         data.l=16;
         data.a = col;
         for (name=datanames; *name; name++) {
            sprintf(keybuf,"%s_COL=", *name);
            nc = userftext(data, 2, keybuf, " ");
            if (nc) {
               col[nc] = '\0';
               wkeyf("%s=file(%s,%s,%s)", *name, file, col, rows);
            }
         }
      }
   } else {
                                          /* xxx_COL= keyword seen */
      if (*file && *rows) {
         data.l=16; 
         data.a = col;
         nc = userftext(data, 2, key, " ");
         if (nc) {
            col[nc] = '\0';
            wkeyf("%s=file(%s,%s,%s)", dataname, file, col, rows);
         } else {
            wkeyf("%s=", dataname);
         }
      }
   }
}

/* -------------------------------------------------------------------------- */
/*                                 tabledata                                   */
/* -------------------------------------------------------------------------- */
/*  Keyword handler for IMP_SET=, IMP_TAB=, IMP_TROWS= and xxx_TCOL=.
 *  Convert set, table, row and column information to 'standard' keywords.
 */
static void tabledata(ident id, char *key, int code, void *arg)
{
   static bool initialized=FALSE;
   static char set[64]="";
   static char table[64]="";
   static char rows[64]="";
   fchar  data;
   int    nc;
   char   **name;
   char   *dataname=(char*)arg;
   char keybuf[KEYLEN+1];
   char col[16+1];

   data.l = 63;
   if (!initialized) {
      data.a = set;
      nc = userftext(data, 2, "IMP_SET=", " ");
      data.a[nc] = '\0';
      data.a = table;
      nc = userftext(data, 2, "IMP_TAB=", " ");
      data.a[nc] = '\0';
      data.a = rows;
      nc = userftext(data, 2, "IMP_TROWS=", " ");
      data.a[nc] = '\0';
      initialized = TRUE;
   }
   data.a = NULL;
   if (!strcmp(key,"IMP_SET=")) {
      data.a = set;
   } else if (!strcmp(key,"IMP_TAB=")) {
      data.a = table;
   } else if (!strcmp(key,"IMP_TROWS=")) {
      data.a = rows;
   }
   if (data.a) {
                                  /* IMP_SET, IMP_TAB or IMP_TROWS= keyword seen */
      nc = userftext(data, 2, key, " ");
      data.a[nc] = '\0';
      if (!nc && data.a==rows) strcpy(rows,":");
      if (*set && *table && *rows) {
         data.l=16;
         data.a = col;   
         for (name=datanames; *name; name++) {
            sprintf(keybuf,"%s_TCOL=", *name);
            nc = userftext(data, 2, keybuf, " ");
            if (nc) {
               col[nc] = '\0';
               wkeyf("%s=table(%s,%s,%s,%s)",  *name, set, table, col, rows);
            }
         }
      }
   } else {
                                          /* xxx_TCOL= keyword seen */
      if (*set && *table && *rows) {
         data.l=16;
         data.a = col;
         nc = userftext(data, 2, key, " ");
         if (nc) {
            col[nc] = '\0';
            wkeyf("%s=table(%s,%s,%s,%s)", dataname, set, table, col, rows);
         } else {
            wkeyf("%s=", dataname);
         }
      }
   }
}

/* -------------------------------------------------------------------------- */
/*                                 popshell                                   */
/* -------------------------------------------------------------------------- */
/*  Keyword handler causing the shell to be popped up or down.
 */
static void popshell(ident id, char *key, int code, void *arg)
{  
   static bool state=FALSE;
   bool button=toflog(FALSE);
 
   (void)userflog(&button, 1, 2, key, " ");
   if (tobool(button)) {
      GgiShowShell((ident)arg, !state);
      state = !state;
      wkeyf(key);
   }
}

/* ========================================================================== */
/*                                 Import                                     */
/* -------------------------------------------------------------------------- */
/*  Import() creates a pop-up shell with menus for obtaining imported data.
 */
void Import(char *key)
{
   static ident s_import;
   ident s_prev, frame_box, file_box, table_box, generic_box;
   ident closebut;
   char **name;
   char keybuf[KEYLEN+1];
    
   s_import = GgiShell(NAME);
   s_prev = GgiUseShell(s_import);
   frame_box = GgiForm(NAME, 2);
   closebut = GgiSetLabel(GgiButton(key, "Close window"), "CLOSE", 0);
   GgiSetPosition(frame_box, 0, NULL, 0, NULL);
   GgiSetPosition(closebut, 0, NULL, 0, frame_box);
   
   GgiUseShell(frame_box);
   file_box    = GgiForm(NAME, 1);
   table_box   = GgiForm(NAME, 1);
   generic_box = GgiForm(NAME, 1);
   GgiSetPosition(file_box, 0, NULL, 0, NULL);
   GgiSetPosition(table_box,   0, NULL, 0, file_box);
   GgiSetPosition(generic_box, 0, NULL, 0, table_box);

   {
      ident title, file, rows, xcol, ycol, ecol;
      (void)GgiUseShell(file_box);
      title = GgiLabel("Input from text file\n ");
      file = GgiSetLabel(GgiTextField("IMP_FILE=", H_FILE, 40),
                             "Filename", 80);
      xcol = GgiSetLabel(GgiTextField("PROFX_COL=", H_FXC, 2),
                          "Column X", 80);
      ycol = GgiSetLabel(GgiTextField("PROFY_COL=", H_FYC, 2),
                          "Y", 50);
      ecol = GgiSetLabel(GgiTextField("PROFYERR_COL=", H_FEC, 2),
                          "Errors", 50);
      rows = GgiSetLabel(GgiTextField("IMP_FROWS=", H_ROWS, 8),
                          "Rows", 80);

      GgiSetPosition(title, 0, NULL, 0, NULL);
      GgiSetPosition(file, 0, NULL, 0, title);
      GgiSetPosition(xcol, 0, NULL, 0, file);
      GgiSetPosition(ycol, 0, xcol, 0, file);
      GgiSetPosition(ecol, 0, ycol, 0, file);
      GgiSetPosition(rows, 0, NULL, 0, xcol);

      (void)ScheduleKeyevent(filedata, "IMP_FILE=", KEYCHANGE, NULL); 
      (void)ScheduleKeyevent(filedata, "IMP_FROWS=", KEYCHANGE, NULL);
      for (name=datanames; *name; name++) {
         sprintf(keybuf, "%s_COL=", *name);
         (void)ScheduleKeyevent(filedata, keybuf, KEYCHANGE, *name);
      }
   }

   {
      ident title, set, table,  rows, xcol, ycol, ecol;
      (void)GgiUseShell(table_box);
      title = GgiLabel("Input from table\n ");
      set  = GgiSetLabel(GgiTextField("IMP_SET=", H_SET, 40),
                             "Set", 80);
      table= GgiSetLabel(GgiTextField("IMP_TAB=", H_TAB, 40),
                             "Table", 80);
      xcol = GgiSetLabel(GgiTextField("PROFX_TCOL=", H_SXC, 8),
                          "Column X", 80);
      ycol = GgiSetLabel(GgiTextField("PROFY_TCOL=", H_SYC, 8),
                          "Y", 0);
      ecol = GgiSetLabel(GgiTextField("PROFYERR_TCOL=", H_SEC, 8),
                          "Errors", 0);
      rows = GgiSetLabel(GgiTextField("IMP_TROWS=", H_ROWS, 8),
                          "Rows", 80);
       
      GgiSetPosition(title, 0, NULL, 0, NULL);
      GgiSetPosition(set,   0, NULL, 0, title);
      GgiSetPosition(table, 0, NULL, 0, set);
      GgiSetPosition(xcol,  0, NULL, 0, table);
      GgiSetPosition(ycol,  0, xcol, 0, table);
      GgiSetPosition(ecol,  0, ycol, 0, table);
      GgiSetPosition(rows,  0, NULL, 0, xcol);

      (void)ScheduleKeyevent(tabledata, "IMP_SET=",   KEYCHANGE, NULL);
      (void)ScheduleKeyevent(tabledata, "IMP_TAB=",   KEYCHANGE, NULL);
      (void)ScheduleKeyevent(tabledata, "IMP_TROWS=", KEYCHANGE, NULL);
      for (name=datanames; *name; name++) {
         sprintf(keybuf, "%s_TCOL=", *name);
         (void)ScheduleKeyevent(tabledata, keybuf, KEYCHANGE, *name);
      }
   }
   
   {
      ident title, xvalues, yvalues, errors;
      (void)GgiUseShell(generic_box);
      title = GgiLabel("Generic input format\n ");
      xvalues = GgiSetLabel(GgiTextField("PROFX=", NULL, 40),
                            "X values", 80);
      yvalues = GgiSetLabel(GgiTextField("PROFY=", NULL, 40),
                            "Y values", 80);
      errors  = GgiSetLabel(GgiTextField("PROFYERR=",  NULL, 40),
                            "Errors", 80);
      GgiSetPosition(title,   0, NULL, 0, NULL);
      GgiSetPosition(xvalues, 0, NULL, 0, title);
      GgiSetPosition(yvalues, 0, NULL, 0, xvalues);
      GgiSetPosition(errors,  0, NULL, 0, yvalues);
   }
   (void)GgiUseShell(s_import);
   (void)GgiUseShell(s_prev);
   ScheduleKeyevent(popshell, key, KEYCHANGE, s_import);
}
#<

#>            filtopt.c
/* filtopt.c
                              COPYRIGHT (c) 1999
                        Kapteyn Astronomical Institute
         University of Groningen  -  9700 AV Groningen, The Netherlands


Author: J.P. Terlouw
*/

#define NAME "Filter options"

#define H_AMPRANGE "Range in amplitudes. Filter is on if one (clip) or two \n\
(range) values are entered."
#define H_VELRANGE "Range in velocity. Filter is on if one (clip) or two \n\
(range) values are entered (physical coordinates)."
#define H_DISPRANGE "Range in dispersions. Filter is on if one (clip) or two \n\
(range) values are entered (physical coordinates)."

#include "stddef.h"
#include "stdio.h"
#include "gipsyc.h"
#include "userfio.h"
#include "wkey.h"
#include "ggi.h"
#include "filtopt.h"

#define New(type) ((type *)calloc(1,sizeof(type)))
#define NNew(n,type) ((type *)calloc(n,sizeof(type)))
#define Delete(x)  {free(x); x=NULL;}

static bool  initialized=FALSE;         /* initialized flag */
static ident s_filtopt;                  /* hidden keywords shell     */
static ident i_box;

/* -------------------------------------------------------------------------- */
/*                                 initialize                                 */
/* -------------------------------------------------------------------------- */
/*  Create pop-up shell and its contents; schedule the appropriate handlers.
 */
static void initialize(char *key)
{
   ident closebut;
   ident amprange, velrange, disprange;

   s_filtopt = GgiShell(NAME);
   (void)GgiUseShell(s_filtopt);
   closebut = GgiSetLabel(GgiButton(key, "Close window"), "CLOSE", 0);
   i_box = GgiForm(NAME, 2);
   (void)GgiUseShell(i_box);

   amprange = GgiSetLabel(GgiTextField("AMPRANGE=", H_AMPRANGE, 10), "Amp. range", 80);
   velrange = GgiSetLabel(GgiTextField("VELRANGE=", H_VELRANGE, 10), "Vel. range", 80);
   disprange = GgiSetLabel(GgiTextField("DISPRANGE=", H_DISPRANGE, 10), "Disp. range", 80);

   GgiSetPosition(amprange,  0, NULL,    0, NULL);
   GgiSetPosition(velrange,  0, amprange, 0, NULL);
   GgiSetPosition(disprange,  0, velrange, 0, NULL);

   (void)GgiUseShell(s_filtopt);
   GgiSetPosition(i_box,    0, NULL, 0, NULL);
   GgiSetPosition(closebut, 0, NULL, 0, i_box);
   (void)GgiUseShell(NULL);
   initialized = TRUE;
}

/* -------------------------------------------------------------------------- */
/*                                 filtopt                                     */
/* -------------------------------------------------------------------------- */
/*  Keyword handler causing the shell to be popped up or down.
 */
static void filtopt(ident id, char *key, int code, void *arg)
{
   static bool state=FALSE;
   bool button=FALSE;

   (void)userflog(&button, 1, 2, key, " ");
   if (tobool(button)) {
      GgiShowShell(s_filtopt, !state);
      state = !state;
      wkey_c(tofchar(key));
   }
}

/* ========================================================================== */
/*                                 Filtopt                                     */
/* -------------------------------------------------------------------------- */
/*  Filtopt() creates a pop-up shell for miscellaneous options. The argument
 *  is the GIPSY keyword (type LOGICAL) which, when set to "YES", will
 *  cause the shell to be popped up or down.
 *  The same keyword is used for the <CLOSE> button in the shell.
 */
extern ident Filtopt(char *key)
{
   initialize(key);
   ScheduleKeyevent(filtopt, key, KEYCHANGE, NULL); 
   return s_filtopt;
}

#<

#>            param.c
/* param.c
                              COPYRIGHT (c) 1999
                        Kapteyn Astronomical Institute
         University of Groningen  -  9700 AV Groningen, The Netherlands

Author: J.P. Terlouw
*/

#define NAME "Parameters"

#include "stddef.h"
#include "stdio.h"
#include "gipsyc.h"
#include "string.h"
#include "userfio.h"
#include "ggi.h"
#include "radiogroup.h"
#include "param.h"

#define H_EST "Estimate"
#define H_FIX "Fix parameter?"

static char popkey[KEYLEN];

static char *gauss_title[]={"Amplitude", "Center", "Dispersion", NULL};
static char *gh1_title[]={"a", "b", "c", "h3", NULL};
static char *gh2_title[]={"a", "b", "c", "h3", "h4", NULL};
static char *voigt_title[]={"Area", "Center",
                            "Doppler half width", "Lorentz half width", NULL};
static char *back_title[]={"Constant term",
                           "Linear term",
                           "Quadratic term", NULL};
static char **titles[]={gauss_title, gh1_title, gh2_title, voigt_title};

static char *gauss_name[]={"AMP", "CENTER", "DISP"};
static char *gh1_name[]={"AMP", "CENTER", "DISP", "H3"};
static char *gh2_name[]={"AMP", "CENTER", "DISP", "H3", "H4"};
static char *voigt_name[]={"AMP", "CENTER", "DISP", "LORENTZ"};
static char *back_name[]={"Z0", "Z1", "Z2"};
static char **names[]={gauss_name, gh1_name, gh2_name, voigt_name};

static ident shell=NULL;
static ident funid=NULL, ncid=NULL;


/* -------------------------------------------------------------------------- */
/*                                 obsolete                                   */
/* -------------------------------------------------------------------------- */
/*  Force shell to be popped down due to function or # of parameters change.
 */
static void obsolete(ident id, char *key, int code, void *arg)
{
   wkeyf("%sYES", popkey);
}

/* -------------------------------------------------------------------------- */
/*                                 createshell                                */
/* -------------------------------------------------------------------------- */
/*  Create the parameter pop-up shell.
 */
static void createshell(char *key)
{
   int i, j;
   ident s_prev, f_box, colbox = 0, currel, currbox, topref;
   ident closebut;
   fint nread, fun, rows=0;
   char *parname, **title, **prefix;
   char  estkey[KEYLEN], fixkey[KEYLEN];
   ident est, fix;
    
   nread = userfint(&fun, 1, 2, "FUNCTION=", " ");
   if (!nread) return;
   (void)userfint(&rows, 1, 2, "NCOMP=", " ");
   if (rows<=0) return;
      
   title = titles[fun];
   prefix = names[fun];
   shell = GgiShell(NAME);
   s_prev = GgiUseShell(shell);
   closebut = GgiSetLabel(GgiButton(key, NULL), "CLOSE", 0);
   f_box = GgiForm(NAME, 2);

/*
 *   Parameter fields
 */
   currbox = NULL;
   for (i=0; title[i]; i++) {
      (void)GgiUseShell(f_box);
      colbox = GgiForm(" ", 0);
      GgiSetPosition(colbox, 0, currbox, 0, NULL);
      currbox = colbox;
      (void)GgiUseShell(colbox);
      parname = prefix[i];
      currel = GgiLabel(title[i]);
      GgiSetPosition(currel, 0, NULL, 0, NULL);
      for (j=0; j<rows; j++) {
         sprintf(estkey, "%s_EST%d=", parname, j);
         sprintf(fixkey, "%s_FIX%d=", parname, j);
         est = GgiSetLabel(GgiTextField(estkey, H_EST, 10), " ", 1);
         fix = GgiSetLabel(GgiButton(fixkey, H_FIX), "FIX", 0);
         GgiSetPosition(est, 0, NULL, 5, currel);
         GgiSetPosition(fix, 0, est,  5, currel);
         currel = est;
      }
   }

/*
 *   Background fields
 */   
   (void)GgiUseShell(f_box);
   topref = colbox;
   title   = back_title;
   currbox = NULL;
   currel = GgiLabel("Background:");
   GgiSetPosition(currel, 0, NULL, 20, topref);
   topref = currel;
   for (i=0; title[i]; i++) {
      (void)GgiUseShell(f_box);
      colbox = GgiForm(" ", 0);
      GgiSetPosition(colbox, 0, currbox, 0, topref);
      currbox = colbox;
      (void)GgiUseShell(colbox);
      parname = back_name[i];
      currel = GgiLabel(title[i]);
      GgiSetPosition(currel, 0, NULL, 0, NULL);
      sprintf(estkey, "%s_EST=", parname);
      sprintf(fixkey, "%s_FIX=", parname);
      est = GgiSetLabel(GgiTextField(estkey, H_EST, 10), " ", 1);
      fix = GgiSetLabel(GgiButton(fixkey, H_FIX), "FIX", 0);
      GgiSetPosition(est, 0, NULL, 5, currel);
      GgiSetPosition(fix, 0, est,  5, currel);
   }
   
   topref = colbox;
   (void)GgiUseShell(f_box);
   ncid   = ScheduleKeyevent(obsolete, "NCOMP=", KEYCHANGE, NULL);
   funid  = ScheduleKeyevent(obsolete, "FUNCTION=", KEYCHANGE, NULL);
   
   (void)GgiUseShell(shell);
   GgiSetPosition(f_box, 0, NULL, 0, NULL);
   GgiSetPosition(closebut, 0, NULL, 0, f_box);
   (void)GgiUseShell(s_prev);
   GgiShowShell(shell, TRUE);
}

/* -------------------------------------------------------------------------- */
/*                                 popshell                                   */
/* -------------------------------------------------------------------------- */
/*  Keyword handler causing the shell to be popped up or down.
 */
static void popshell(ident id, char *key, int code, void *arg)
{
   static bool state=FALSE;
   bool button=toflog(FALSE);

   (void)userflog(&button, 1, 2, key, " ");
   if (tobool(button)) {
      if (state) {
         DescheduleKeyevent(&ncid);
         DescheduleKeyevent(&funid);
         GgiDelete(&shell);
         shell = NULL;
         state = FALSE;
      } else {
         createshell(key);
         state = (bool)shell;
      }
      wkeyf(key);
   }
}

/* ========================================================================== */
/*                                 Param                                      */
/* -------------------------------------------------------------------------- */
/*  Initialize the parameter pop-up window.
 */
void Param(char *key)
{
   ScheduleKeyevent(popshell, key, KEYCHANGE, NULL);
   wkeyf(key);
   strcpy(popkey,key);
}
    
   
#<

#>            parcol.c
/* parcol.c
                              COPYRIGHT (c) 1999
                        Kapteyn Astronomical Institute
         University of Groningen  -  9700 AV Groningen, The Netherlands

Author: J.P. Terlouw
*/

#include "stddef.h"
#include "stdio.h"
#include "gipsyc.h"
#include "userfio.h"
#include "ggi.h"
#include "parcol.h"

#define H_EST "Estimate"
#define H_FIX "Fix parameter?"
#define H_FIT "Fit result"

ident BuildColumn(int nrows, char *title, char *parname)
{
   int i;
   ident column, current, prvform;
   char  estkey[KEYLEN], fixkey[KEYLEN], fitkey[KEYLEN];
   ident est, fix, fit;
   
#if 1
   column = GgiForm("column", 0);
   prvform = GgiUseShell(column);
#endif
   current = GgiLabel(title);
#if 1
   GgiSetPosition(current, 0, NULL, 0, NULL);
#else
   column = current;
#endif
   
   for (i=0; i<nrows; i++) {
      sprintf(estkey, "%s_EST%d=", parname, i);
      sprintf(fixkey, "%s_FIX%d=", parname, i);
      sprintf(fitkey, "%s_FIT%d=", parname, i);
      est = GgiSetLabel(GgiTextField(estkey, H_EST, 10), " ", 1);
      fix = GgiSetLabel(GgiButton(fixkey, H_FIX), "FIX", 0);
      fit = GgiSetLabel(GgiTextField(fitkey, H_FIT, 10), " ", 1);
      GgiSetPosition(est, 0, NULL, 0, current);
      GgiSetPosition(fix, 0, est,  0, current);
      GgiSetPosition(fit, 0, NULL, -1, est);
      current = fit;
   }
#if 1
   (void)GgiUseShell(prvform);
#endif
   return column;
}
#<

#>            outputresults.c
/* outputresults.c
                            COPYRIGHT (c) 1999
                      Kapteyn Astronomical Institute
                University of Groningen, The Netherlands
                           All Rights Reserved.
*/

#include    "stdio.h"
#include    "stdlib.h"       /* Defines the ANSI C functions for number */
                             /* conversion, storage allocation, and similar tasks.*/
#include    "string.h"       /* Declares the ANSI C string functions*/
                             /* like:strcpy, strcat etc.*/
#include    "gipsyc.h"
#include    "math.h"         /* Declares the mathematical functions and macros.*/
#include    "float.h"
#include    "xgaufit.h"
#include    "keyevents.h"
#include    "time.h"
#include    "minmax1.h"
#include    "userfio.h"

#include    "gui.h"
#include    "ggi.h"
#include    "outputresults.h"


double          gausshermiteh3h4( double, double *, int );
extern  double  dblank; 





static void getposition( int  count, 
                         fint subdim, 
                         fint *blo, 
                         fint *bhi, 
                         fint *gridpos )
/*-----------------------------------------------------------*/
/* PURPOSE: Convert pixel counter to 'sub'-dim grid position.*/
/*-----------------------------------------------------------*/
{
   int  len[MAXAXES];
   int  i;

   len[0] = 1;
     
   for (i = 1; i < subdim; i++)
   {
      len[i] = len[i-1] * (bhi[i-1] - blo[i-1] + 1);
   }
   for (i = subdim - 1; i >= 0; i--) 
   {
      gridpos[i] = count/len[i];      /* Divide by size of (part of) subset */
      count -= gridpos[i] * len[i];   /* Decrease dim of subset */
      gridpos[i] += blo[i];           /* Adjust for starting point */
   }   
}


static void storefitteddata( float  *data,
                             int     count,
                             fint   *blo,
                             fint   *bhi,
                             fint    axlo,
                             fint    axhi,
                             fint    subdim,
                             fchar   Setout )
/*-----------------------------------------------------------*/
/* PURPOSE: Store data of a profile in a set on disk.        */
/*                                                           */
/* The profile has an index number which starts with 0 ends  */
/* with the number of profiles - 1. This number is converted */
/* to a vector in the structure that is left after the       */
/* profile axis is subtracted. The last dimension of this    */
/* vector is the start and end grid of the profile. Note     */
/* that the output set must be a copy of the inset and its   */
/* axis lengths are given by the input box and the profile   */
/* start- and end grids. Not also that this method does not  */
/* contain any buffering and is therefore not the fastest    */
/* solution to the problem.                                  */
/*-----------------------------------------------------------*/
{
   fint   gridpos[MAXAXES];
   int    i;
   fint   setlevel = 0;
   fint   glo[MAXAXES];
   fint   ghi[MAXAXES];
   fint   cwlo, cwhi;
   fint   IOlen;
   fint   writepixels;
   fint   tid = 0;

   
   getposition( count, subdim, blo, bhi, gridpos );
   for (i = 0; i < subdim; i++)
   {
      glo[i] = ghi[i] = gridpos[i];
   }
   glo[subdim] = axlo;                  /* Add last dimension */
   ghi[subdim] = axhi;
    
   cwlo = gdsc_fill_c( Setout, &setlevel, glo );
   cwhi = gdsc_fill_c( Setout, &setlevel, ghi );
   IOlen = axhi - axlo +1;
   gdsi_write_c( Setout,
                 &cwlo, &cwhi,
                 data,                  /* The profile */
                 &IOlen,                /* Lenght of profile axis */
                 &writepixels,          /* Output variable */
                 &tid );
}




static void getredchi2( globaltype *g )

/*------------------------------------------------------------*/
/* PURPOSE: Calculate Chi-squared for the gauss  parameters.  */
/*------------------------------------------------------------*/
{
   int     i;
   int     ndat = 0;
   fint    npar;
   int     freepars;
   fint    opt;
   double  chi, chi2 = 0;
   double  sum = 0.0;
   double  Yfit_dble;
   

   npar = g->gaussiansfound * g->nvar + 3;

   freepars = 0;
   for (i = 0; i < npar; i++)
   {
      if (g->fixedfree[i] != 0)
         freepars++;
   }

   opt = FIEID*g->function  + g->gaussiansfound;

   for (i = 0; i < g->profilelen; i++)
   {
      double   delta;
      double   Yi = g->profiles[g->currentprofilenumber][i];
      double   Xi = g->xphys[i];

      Yfit_dble = funcd_c( &Xi, g->estimates, &npar, &opt );
      g->Yfit[i] = (float) Yfit_dble; 
      delta = Yi - Yfit_dble;
      if (g->wdat[i] != 0.0)
      {
         chi = delta * g->wdat[i];
         sum += delta * delta;
         ndat++;
      }
      else
      {
         delta = 0.0;
         chi = 0.0;
      }
      chi2 += chi * chi;
      g->Yres[i] = delta;
   }
   
   g->chi2  = chi2;
   /*----------------------------------------*/
   /* The weights mode is uniform then       */
   /* the weights are equal to 1. For the    */
   /* chi2 reduced we have to modify the     */
   /* chi2 for sigma == rms.                 */
   /*----------------------------------------*/
   chi2 /= (g->rms)*(g->rms);
   
   g->rchi2 = chi2 / (g->profilelen-freepars-1);

   /*----------------------------------------*/
   /* Redo the loop over the residuals to    */
   /* calculate the standard deviation.      */
   /*----------------------------------------*/
   if (ndat > 1)
   {
      double averageresidual = sum / (double) ndat;
   
      sum = 0.0;                              /* Reset */
      for (i = 0; i < g->profilelen; i++)
      {
         if (g->wdat[i] != 0.0)
            sum += (g->Yres[i] - averageresidual) * (g->Yres[i] - averageresidual);
      }
      g->sigmaresidual = sqrt( sum / (double) (ndat-1) );
   }  
   else
   {
      g->sigmaresidual = 0.0; 
   }  
   
   if (g->gausetexist)
      storefitteddata( g->Yfit,
                       g->count,
                       g->blo,
                       g->bhi,
                       g->plohi[0], g->plohi[1],
                       g->subdim,
                       g->Setgau );
}




static double amp2area( double Amp,
                        double aD,
                        double aL )
/*------------------------------------------------------------*/
/* PURPOSE: Convert an amplitude to the area for a Voigt      */
/*          function.                                         */
/*------------------------------------------------------------*/
{
   double x, y, Rw, Iw;

   x = 0.0;
   y = aL * sqrt(log(2.0)) / aD;
   w( x, y, &Rw, &Iw );

   return(  aD * sqrt(PI/log(2.0)) / Rw );
}






static double fhergau( double  x,
                       double *A )
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate the polynomial for which to find a zero. */
/*------------------------------------------------------------*/
{
    return( x*(x*(x*(x*(x*A[5]+A[4])+A[3])+A[2])+A[1])+A[0] );
}




static double bisection( double *A,
                         double x1,
                         double x2,
                         double xacc,
                         int    maxiters )
/*------------------------------------------------------------*/
/* PURPOSE: Use simple bisection to get zero value.           */
/*------------------------------------------------------------*/
{
   int    j;
   double dx, f, fmid, xmid, rtb;

   f = fhergau( x1, A);
   fmid = fhergau( x2, A );
   if (f*fmid >= 0.0)
   {
      return( dblank );
   }

   rtb = f < 0.0 ? (dx=x2-x1,x1) : (dx=x1-x2,x2);
   for (j = 1; j <= maxiters; j++)
   {
      fmid = fhergau( xmid=rtb+(dx *= 0.5), A );
      if (fmid <= 0.0)
         rtb = xmid;
      if (fabs(dx) < xacc || fmid == 0.0)
         return rtb;
   }
   return( dblank );
}




static void findhermitemax( double a,
                            double b,
                            double c,
                            double h3,
                            double h4,
                            double *gamp,
                            double *gcent )
/*------------------------------------------------------------*/
/* PURPOSE: Find the real max of the GH series.               */
/*                                                            */
/* We use the first derivative of the Gauss-Hermite function  */
/* to find the maximum around 'cent'.                         */
/* If F = (x-b)/c then the function for which we want the     */
/* the zero's is A0+A1*F+A2*F^2+A3*F^3+A4*F^4+A5*F^5 = 0      */
/* where A0 = -c1.h3  A1 = h4*(c0-2.0*c2) + 1.0               */
/* A2 = h3*(c1-3.0*c3) A3 = h4*(c2-4.0*c4) A4=c3.h3 A5=c4.h4  */
/*                                                            */
/* c0 = 1/4sqrt(6) c1 = -sqrt(3) c2 = -sqrt(6)                */
/* c3 = 2/3sqrt(3) c4 = 1/3sqrt(6)                            */
/*                                                            */
/*------------------------------------------------------------*/
{
   double  c0 = sqrt(6.0)/4.0, c1 = -sqrt(3.0), c2 = -sqrt(6.0);
   double  c3 = 2.0*sqrt(3.0)/3.0, c4 = sqrt(6.0)/3.0;
   double  A[6];
   double  xbix1 = b - c/2.0;
   double  xbis2 = b + c/2.0;
   double  acc = 1.0e-8;
   double  fpar[7];
   double  F = 0.0;
   int     maxiters = 100;


   A[0] = -c1*h3;
   A[1] = h4*(c0-2.0*c2) + 1.0;
   A[2] = h3*(c1-3.0*c3);
   A[3] = h4*(c2 - 4.0*c4);
   A[4] = c3*h3;
   A[5] = c4*h4;
   xbix1 = -0.5; xbis2 = 0.5;
   F = bisection( A, xbix1, xbis2, acc, maxiters );
   if (F == dblank)
   {
      *gcent = *gamp = dblank;
   }
   else
   {
      *gcent = c * F + b;
      fpar[0] = a; fpar[1] = b; fpar[2] = c; fpar[3] = h3; fpar[4] = h4;
      fpar[5] = fpar[6] = fpar[7] = 0.0;
      *gamp = gausshermiteh3h4( *gcent, fpar, 1 );        /* One component */
   }
}



static double getflux( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Get the area under the curves.                    */
/*------------------------------------------------------------*/
{
   double     flux= 0.0, sumflux = 0.0;
   int        i;


   for (i = 0; i < g->ncomp; i++)
   {
      int    offset = i*g->nvar;
      if (g->function == GAUSS)
      {
         /*--------------------------------------------------*/
         /* flux = area  = |disp| * Amp * sqrt(2PI)  or      */
         /* flux = 0.5 * A * FWHM * sqrt(PI/ln(2))           */
         /*--------------------------------------------------*/
         flux = g->estimates[offset+0] * ABS(g->estimates[offset+2])
                * sqrt(2.0*PI);
      }
      else if (g->function == GAUSSHERMITE1)
      {
         /*--------------------------------------------------*/
         /* Area = a.c.sqrt(2PI)                             */
         /*--------------------------------------------------*/
         flux = g->estimates[offset+0] * ABS(g->estimates[offset+2])
                * sqrt(2.0*PI);
      }
      else if (g->function == GAUSSHERMITE2)
      {
         /*--------------------------------------------------*/
         /* Area = a.c.sqrt(2PI) * (1+h4*sqrt(6)/4)          */
         /*--------------------------------------------------*/
         double a = g->estimates[offset+0];
         double c = g->estimates[offset+2];
         double h4 = g->estimates[offset+4];

         flux = a*c*sqrt(2.0*PI) * (1.0+h4*sqrt(6.0)/4.0);
      }
      if (g->function == VOIGT)
      {
         /*--------------------------------------------------*/
         /* First parameter is area                          */
         /*--------------------------------------------------*/
         flux = g->estimates[offset+0];
      }
      sumflux += flux;
   }
   return( sumflux );
}



static int filter( double amp,
                   double vel,
                   double disp,
                   double *amprange,
                   double *velrange,
                   double *disprange )
/*------------------------------------------------------------*/
/* PURPOSE: Must these values be filtered?                    */
/*------------------------------------------------------------*/
{
   if (amprange[0] != dblank)
   {
      if (amprange[1] == dblank)
      {
         if (amp < amprange[0])
         {
            return( 1 );            
         }
      }
      else
      {
         if (amp < amprange[0] || amp > amprange[1])
            return( 1 );
      }
   }

   if (velrange[0] != dblank)
   {
      if (velrange[1] == dblank)
      {
         if (vel < velrange[0])
         {
            return( 2 );            
         }
      }
      else
      {
         if (vel < velrange[0] || vel > velrange[1])
            return( 2 );
      }
   }

   if (disprange[0] != dblank)
   {
      if (disprange[1] == dblank)
      {
         if (disp < disprange[0])
         {
            return( 3 );            
         }
      }
      else
      {
         if (disp < disprange[0] || disp > disprange[1])
            return( 3 );
      }
   }

   /* Need not to be filtered out */
   return( NO );
}



extern void outputresults( globaltype *g )
/*------------------------------------------------------------*/
/* PURPOSE: Display fit results.                              */
/*------------------------------------------------------------*/
{
   int          i;
   double       flux;
   int          pnum = g->currentprofilenumber; 
   int          subnr;
   int          dofilter = NO;


   /* The estimates for all the components are written already. */
   subnr = 3 * g->ncomp; 
  
   for (i = 0; i < g->ncomp; i++)
   {
      int    offset = i*g->nvar;

      dofilter = filter( g->estimates[offset+0], 
                         g->estimates[offset+1], 
                         g->estimates[offset+2], 
                         g->filter.amp,
                         g->filter.vel,
                         g->filter.disp );
      if (dofilter) 
      /*--------------------------------------------------*/
      /* This component is filtered by the user. All      */
      /* variables in this component are set to blank.    */
      /*--------------------------------------------------*/      
      {
         int j;
         for (j = 0; j < (g->fitpars+g->errpars); j++)
         {
            g->results[subnr++][pnum] = dblank;            
         }
         if      (dofilter == 1) g->log.ampfiltered++;
         else if (dofilter == 2) g->log.velfiltered++;
         else if (dofilter == 3) g->log.dispfiltered++;         
      }
      else
      {
        switch (g->function)
        {
         case GAUSS:
         {
            /*--------------------------------------------------*/
            /* Line strength Y:                                 */
            /* Y  = A.s.sqrt(2Pi)                               */
            /* dY = sqrt(2Pi) * sqrt{s^2.dA^2+A^2.ds^2}         */
            /*--------------------------------------------------*/
            {
               double A  = g->estimates[offset+0];
               double dA = ABS(g->errorpar[offset+0]);
               double s  = g->estimates[offset+2];
               double ds = ABS(g->errorpar[offset+2]);
               double Y, dY;

               Y  = A * s * sqrt(2.0*PI);
               dY = sqrt(2.0*PI) * sqrt( s*s*dA*dA + A*A*ds*ds );
               
               g->results[subnr++][pnum] =  Y;
               g->results[subnr++][pnum] =  dY;               
            }            
            /*--------------------------------------------------*/
            /* Mean M = x0                                      */
            /*--------------------------------------------------*/
            {
               double M  = g->estimates[offset+1];
               double dM = ABS(g->errorpar[offset+1]);

               g->results[subnr++][pnum] =  M;
               g->results[subnr++][pnum] =  dM;
            }
            /*--------------------------------------------------*/
            /* Dispersion D = s                                 */
            /*--------------------------------------------------*/
            {
               double D  = ABS(g->estimates[offset+2]);
               double dD = ABS(g->errorpar[offset+2]);

               g->results[subnr++][pnum] =  D;
               g->results[subnr++][pnum] =  dD;
            }
            /*--------------------------------------------------*/
            /* Amplitude = A                                    */
            /*--------------------------------------------------*/
            {
               double A  = g->estimates[offset+0];
               double dA = ABS(g->errorpar[offset+0]);

               g->results[subnr++][pnum] =  A;
               g->results[subnr++][pnum] =  dA;
            }
            break;
         }
         case GAUSSHERMITE1:
         {
            /*--------------------------------------------------*/
            /* First list the parameters a,b,c & h3 of the      */
            /* function. Do NOT list in the gui.                */
            /*--------------------------------------------------*/
            {
               double  a = g->estimates[offset+0];
               double  da = ABS(g->errorpar[offset+0]);
                           
               g->results[subnr++][pnum] =  a;
               g->results[subnr++][pnum] =  da;
            }
            {
               double  b = g->estimates[offset+1];
               double  db = ABS(g->errorpar[offset+1]);
                           
               g->results[subnr++][pnum] =  b;
               g->results[subnr++][pnum] =  db;
            }
            {
               double  c = g->estimates[offset+2];
               double  dc = ABS(g->errorpar[offset+2]);
                           
               g->results[subnr++][pnum] =  c;
               g->results[subnr++][pnum] =  dc;
            }
            {
               double  h3 = g->estimates[offset+3];
               double  dh3 = ABS(g->errorpar[offset+3]);
                           
               g->results[subnr++][pnum] =  h3;
               g->results[subnr++][pnum] =  dh3;
            }
            
            
            /*--------------------------------------------------*/
            /* Line strength Y:                                 */
            /* Y  = a.c.sqrt(2Pi) {1+h4.sqrt(6)/4} & h4 = 0     */
            /* dY = sqrt(2Pi) * sqrt{c^2.da^2+a^2.dc^2}         */
            /*--------------------------------------------------*/
            {
               double a = g->estimates[offset+0], c = g->estimates[offset+2];
               double da = ABS(g->errorpar[offset+0]), dc = ABS(g->errorpar[offset+2]);
               double Y, dY;

               Y  = a * c * sqrt(2.0*PI);
               dY =  sqrt(2.0*PI) * sqrt( c*c*da*da + a*a*dc*dc );
               
               g->results[subnr++][pnum] =  Y;
               g->results[subnr++][pnum] =  dY;
            }
            /*--------------------------------------------------*/
            /* Mean M ~= b + h3*sqrt(3)                         */
            /*     dM  = sqrt{db^2+3h3^2dc^2}                   */
            /*--------------------------------------------------*/
            {
               double b = g->estimates[offset+1], c = g->estimates[offset+2];
               double db = ABS(g->errorpar[offset+1]), dc = ABS(g->errorpar[offset+2]);
               double h3 = g->estimates[offset+3], dh3 = ABS(g->errorpar[offset+3]);
               double M, dM;

               M  = b + sqrt(3.0)*h3*c;
               dM = sqrt( db*db + 3.0*h3*h3*dc*dc + 3.0*c*c*dh3*dh3 );
               
               g->results[subnr++][pnum] =  M;
               g->results[subnr++][pnum] =  dM;               
            }

            /*--------------------------------------------------*/
            /* Dispersion D ~= c.{1-3h3^2},  h4 = 0. If we agree*/
            /* to give the values only to the lowest order of   */
            /* h3 then:                                         */
            /*            D = c                                 */
            /*           dD = dc                                */
            /*--------------------------------------------------*/
            {
               double c = g->estimates[offset+2], dc = ABS(g->errorpar[offset+2]);
               double D, dD;

               D = ABS(c);
               dD = dc;
               
               g->results[subnr++][pnum] =  D;
               g->results[subnr++][pnum] =  dD;
            }

            /*--------------------------------------------------*/
            /* Coefficient of skewness: xi1 ~= 4*sqrt(3).h3     */
            /*--------------------------------------------------*/
            {
               double h3 = g->estimates[offset+3], dh3 = ABS(g->errorpar[offset+3]);
               double xi, dxi;
               double fac = 4.0 * sqrt(3.0);

               xi  = fac * h3;
               dxi = fac * dh3;
               
               g->results[subnr++][pnum] =  xi;
               g->results[subnr++][pnum] =  dxi ; 
            }
            /*--------------------------------------------------*/
            /* Real maximum and position of real maximum.       */
            /* Use bisection of first derivative to x.          */
            /*--------------------------------------------------*/
            {
               double amp, cent;
               double damp, dcent;

               /* Max of function and position of max. */
               findhermitemax(g->estimates[offset+0],
                              g->estimates[offset+1],
                              g->estimates[offset+2],
                              g->estimates[offset+3],
                              0.0,
                              &amp, &cent );
               damp = ABS(g->errorpar[offset+0]);
               dcent = ABS(g->errorpar[offset+1]);
        
               g->results[subnr++][pnum] =  amp;
               g->results[subnr++][pnum] =  damp ;
               
               g->results[subnr++][pnum] =  cent;
               g->results[subnr++][pnum] =  dcent ;
            }            
            break;
         }
         case GAUSSHERMITE2:
         {
            /*--------------------------------------------------*/
            /* First list the parameters a,b,c & h3 of the      */
            /* function. Do NOT list in the gui.                */
            /*--------------------------------------------------*/
            {
               double  a = g->estimates[offset+0];
               double  da = ABS(g->errorpar[offset+0]);
                           
               g->results[subnr++][pnum] =  a;
               g->results[subnr++][pnum] =  da;
            }
            {
               double  b = g->estimates[offset+1];
               double  db = ABS(g->errorpar[offset+1]);
                           
               g->results[subnr++][pnum] =  b;
               g->results[subnr++][pnum] =  db;
            }
            {
               double  c = g->estimates[offset+2];
               double  dc = ABS(g->errorpar[offset+2]);
                           
               g->results[subnr++][pnum] =  c;
               g->results[subnr++][pnum] =  dc;
            }
            {
               double  h3 = g->estimates[offset+3];
               double  dh3 = ABS(g->errorpar[offset+3]);
                           
               g->results[subnr++][pnum] =  h3;
               g->results[subnr++][pnum] =  dh3;
            }
            {
               double  h4 = g->estimates[offset+4];
               double  dh4 = ABS(g->errorpar[offset+4]);
                           
               g->results[subnr++][pnum] =  h4;
               g->results[subnr++][pnum] =  dh4;
            }            
            /*--------------------------------------------------*/
            /* Line strength Y:                                 */
            /* Y  = a.c.sqrt(2Pi) {1+h4.sqrt(6)/4}              */
            /* set f = {1+h4.sqrt(6)/4} then:                   */
            /* dY = sqrt(2Pi) * sqrt{c^2.f^2.da^2+a^2.f^2.dc^2+ */
            /*                       3/8*a^2.c^2.dh4^2 }        */
            /*--------------------------------------------------*/
            {
               double a = g->estimates[offset+0], c = g->estimates[offset+2];
               double da = ABS(g->errorpar[offset+0]), dc = ABS(g->errorpar[offset+2]);
               double h4 = g->estimates[offset+4], dh4 = ABS(g->errorpar[offset+4]);
               double Y, dY;
               double f;

               f = 1.0 + h4 * sqrt(6.0) / 4.0;
               Y  = a * c * f * sqrt(2.0*PI);
               dY = sqrt(2.0*PI) * sqrt( c*c*f*f*da*da + a*a*f*f*dc*dc +
                                         3.0*a*a*c*c*dh4*dh4/8.0 );
             
               g->results[subnr++][pnum] =  Y;
               g->results[subnr++][pnum] =  dY; 
            }
            /*--------------------------------------------------*/
            /* Mean M ~= b + h3*sqrt(3)                         */
            /*     dM  = sqrt{db^2+3h3^2dc^2}                   */
            /*--------------------------------------------------*/
            {
               double b = g->estimates[offset+1], c = g->estimates[offset+2];
               double db = ABS(g->errorpar[offset+1]), dc = ABS(g->errorpar[offset+2]);
               double h3 = g->estimates[offset+3], dh3 = ABS(g->errorpar[offset+3]);
               double M, dM;

               M  = b + sqrt(3.0)*h3*c;
               dM = sqrt( db*db + 3.0*h3*h3*dc*dc + 3.0*c*c*dh3*dh3 );

               g->results[subnr++][pnum] =  M;
               g->results[subnr++][pnum] =  dM;
            }

            /*--------------------------------------------------*/
            /* Disersion D ~= c.{1+h4*sqrt(6)}                  */
            /*          dD  = sqrt{ (1+h4.sqrt(6))^2.dc^2 +     */
            /*                       6c^2.dh4^2 }               */
            /*--------------------------------------------------*/
            {
               double c = g->estimates[offset+2], dc = ABS(g->errorpar[offset+2]);
               double h4 = g->estimates[offset+4], dh4 = ABS(g->errorpar[offset+4]);
               double D, dD;
               double fac;

               fac = 1.0 + h4*sqrt(6.0);
               D = ABS(c * fac);
               dD = sqrt( fac*fac*dc*dc + 6.0*c*c*dh4*dh4 );
               
               g->results[subnr++][pnum] =  D;
               g->results[subnr++][pnum] =  dD;
            }

            /*--------------------------------------------------*/
            /* Coefficient of skewness: xi1 ~= 4*sqrt(3).h3     */
            /*--------------------------------------------------*/
            {
               double h3 = g->estimates[offset+3], dh3 = ABS(g->errorpar[offset+3]);
               double xi, dxi;
               double fac = 4.0 * sqrt(3.0);

               xi  = fac * h3;
               dxi = fac * dh3;
               
               g->results[subnr++][pnum] =  xi;
               g->results[subnr++][pnum] =  dxi;
            }

            /*----------------------------------------------------*/
            /* Coefficient of kurtosis: xi2 ~= 3 + 6.5*sqrt(6).h4 */
            /*----------------------------------------------------*/
            {
               double h4 = g->estimates[offset+4], dh4 = ABS(g->errorpar[offset+4]);
               double xi, dxi;
               double f = 8.0 * sqrt(6);

               xi  = f * h4;
               dxi = f * dh4;
               
               g->results[subnr++][pnum] =  xi;
               g->results[subnr++][pnum] =  dxi;
            }

            /*--------------------------------------------------*/
            /* Real maximum and position of real maximum.       */
            /* Use bisection of first derivative to x.          */
            /*--------------------------------------------------*/
            {
               double amp, cent;
               double damp, dcent;
               
               findhermitemax(g->estimates[offset+0],
                              g->estimates[offset+1],
                              g->estimates[offset+2],
                              g->estimates[offset+3],
                              g->estimates[offset+4],
                              &amp, &cent );
                              
               damp = ABS(g->errorpar[offset+0]);
               dcent = ABS(g->errorpar[offset+1]); 
               
               g->results[subnr++][pnum] =  amp;
               g->results[subnr++][pnum] =  damp ;
               
               g->results[subnr++][pnum] =  cent;
               g->results[subnr++][pnum] =  dcent ;
            }
            break;
         }
         case VOIGT:
         {
            /*--------------------------------------------------*/
            /* Line strength Y                                  */
            /*--------------------------------------------------*/
            {
               double Y  = g->estimates[offset+0];
               double dY = ABS(g->errorpar[offset+0]);


               g->results[subnr++][pnum] =  Y;
               g->results[subnr++][pnum] =  dY ;
            }
            /*--------------------------------------------------*/
            /* Center                                           */
            /*--------------------------------------------------*/
            {
               double M  = g->estimates[offset+1];
               double dM = ABS(g->errorpar[offset+1]);
               
               g->results[subnr++][pnum] =  M;
               g->results[subnr++][pnum] =  dM;
            }
            /*--------------------------------------------------*/
            /* Doppler FWHM                                     */
            /* Note that fit parameter was a Half Width         */
            /*--------------------------------------------------*/
            {
               double aD =  2.0*ABS(g->estimates[offset+2]);
               double daD = 2.0*ABS(g->errorpar[offset+2]);

               g->results[subnr++][pnum] =  aD;
               g->results[subnr++][pnum] =  daD ; 
            }
            /*--------------------------------------------------*/
            /* Lorentz FWHM                                     */
            /* Note that fit parameter was a Half Width         */
            /*--------------------------------------------------*/
            {
               double aL =  2.0*ABS(g->estimates[offset+3]);
               double daL = 2.0*ABS(g->errorpar[offset+3]);

               g->results[subnr++][pnum] =  aL;
               g->results[subnr++][pnum] =  daL ;  
            }
            /*--------------------------------------------------*/
            /* Amplitude                                        */
            /* Transform Voigt Area to Amplitude                */
            /*--------------------------------------------------*/
            {
               double    aD    = g->estimates[offset+2];    /* The half widths */
               double    daD   = g->errorpar[offset+2];
               double    aL    = g->estimates[offset+3];
               double    Area  = g->estimates[offset+0];
               double    dArea = g->errorpar[offset+0];
               double    A, dA;

               A = Area / amp2area( Area, aD, aL );
               dA = A * sqrt( (daD*daD/aD/aD) + (dArea*dArea/Area/Area) );

               g->results[subnr++][pnum] =  A;
               g->results[subnr++][pnum] =  dA;
            }            
            break;
         }
         default:
         {
            static int first = YES;
            if (first)
            {
               gui_error( "Unknown function!" );
               first = NO;
            }
         }
      }
     }
   }

   if (dofilter)
      return;
   
   
   /*--------------------------------------------------*/
   /* Finally, list  the background parameters.        */
   /*--------------------------------------------------*/
   {
      double   Z0 = g->estimates[g->gaussiansfound*g->nvar];
      double   dZ0 = ABS(g->errorpar[g->gaussiansfound*g->nvar]);
      
      g->results[subnr++][pnum] =  Z0;
      g->results[subnr++][pnum] =  dZ0 ; 
   }
   {
      double   Z1 = g->estimates[g->gaussiansfound*g->nvar+1];
      double   dZ1 = ABS(g->errorpar[g->gaussiansfound*g->nvar+1]);
      
      g->results[subnr++][pnum] =  Z1;
      g->results[subnr++][pnum] =  dZ1 ; 
   }
   {
      double   Z2 = g->estimates[g->gaussiansfound*g->nvar+2];
      double   dZ2 = ABS(g->errorpar[g->gaussiansfound*g->nvar+2]);
      
      g->results[subnr++][pnum] =  Z2;
      g->results[subnr++][pnum] =  dZ2 ; 
   }

   flux = getflux( g );   
   g->results[subnr++][pnum] =  flux;

   g->results[subnr++][pnum] = (double) g->iters;
   
   getredchi2( g );    
   g->results[subnr++][pnum] =  g->rchi2;

   if (g->rchi2 > g->log.rchimax)
      g->log.rchimax = g->rchi2; 

   g->results[subnr++][pnum] =  g->sigmaresidual;
   g->results[subnr++][pnum] =  g->usedQ;
}
#<
