;+
; NAME:
;   ARG_PRESENT
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   Determine whether output parameter has been passed (IDL4 compatibility)
;
; CALLING SEQUENCE:
;   PRESENT = ARG_PRESENT(ARG)
;
; DESCRIPTION: 
;
;   ARG_PRESENT tests whether an argument to a function or procedure
;   can be used as an output parameter.  The behavior of this function
;   is identical to that of the built-in ARG_PRESENT function in IDL
;   version 5 or greater, and is meant to give the same functionality
;   to programs in IDL 4.
;
;   An IDL procedure or function can use ARG_PRESENT to decide whether
;   the value of a positional or keyword parameter will be returned to
;   the calling procedure.  Generally, if the caller did not pass the
;   parameter then there is no need to compute the value to be
;   returned.
;
;   To be a valid output parameter, the caller must have passed a
;   named variable into which the result is stored.  If the caller
;   passed the parameter by value (e.g., an expression or a
;   subscripted array) the value cannot be returned and ARG_PRESENT
;   returns 0.
;
; INPUTS:
;
;   ARG - the parameter to be tested.  It can be either a positional
;         or a keyword parameter.  Passing a normal local variable
;         (i.e., not a passed parameter) will cause ARG_PRESENT to
;         return zero.
;
; RETURNS:
;
;   Returns a value of 1 if ARG is a valid output parameter, and a
;   value of 0 otherwise.
;
;
; EXAMPLE:
;
;   Consider the following procedure:
;      PRO TESTARG, ARG1
;        print, ARG_PRESENT(ARG1)
;      END
;
;   This procedure will print 1 when an ARG1 can be used as an output
;   parameter.  Here are some examples of the results of TESTARG.
;
;      IDL> testarg
;             0
;      IDL> testarg, x      
;             1
;      IDL> testarg, findgen(10)
;             0
;   
;   In the first case, no argument is passed, so ARG1 cannot be a
;   return variable.  In the second case, X is undefined, but it is
;   still a legal named variable capable of receiving an output
;   parameter.  In the third case, FINDGEN(10) is an expression which
;   cannot receive an output parameter.
;
; SEE ALSO:
;
;   ARG_PRESENT in IDL version 5
;
; MODIFICATION HISTORY:
;   Written, CM, 13 May 2000
;   Small documentation and bug fixes, CM, 04 Jul 2000
;
;-
; Copyright (C) 2000, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
;+
; NAME:
;   ARRDELETE
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   Remove a portion of an existing array.
;
; CALLING SEQUENCE:
;   NEWARR = ARRDELETE(INIT, [AT=POSITION,] [LENGTH=NELEM])
;
; DESCRIPTION: 
;
;   ARRDELETE will remove or excise a portion of an existing array,
;   INIT, and return it as NEWARR.  The returned array will never be
;   larger than the initial array.
;
;   By using the keywords AT and LENGTH, which describe the position
;   and number of elements to be excised respectively, any segment of
;   interest can be removed.  By default the first element is removed.
;
; INPUTS:
;
;   INIT - the initial array, which will have a portion deleted.  Any
;          data type, including structures, is allowed.  Regardless of
;          the dimensions of INIT, it is treated as a one-dimensional
;          array.  If OVERWRITE is not set, then INIT itself is
;          unmodified.
;
; KEYWORDS:
;
;   AT - a long integer indicating the position of the sub-array to be
;        deleted.  If AT is non-negative, then the deleted portion
;        will be NEWARR[AT:AT+LENGTH-1].  If AT is negative, then it
;        represents an index counting from then *end* of INIT,
;        starting at -1L.
;        Default: 0L (deletion begins with first element).
;
;  LENGTH - a long integer indicating the number of elements to be
;           removed.  
;
;  OVERWRITE - if set, then INIT will be overwritten in the process of
;              generating the new array.  Upon return, INIT will be
;              undefined.
;
;  COUNT - upon return, the number of elements in the resulting array.
;          If all of INIT would have been deleted, then -1L is
;          returned and COUNT is set to zero.
;
;  EMPTY1 - if set, then INIT is assumed to be empty (i.e., to have
;           zero elements).  The actual value passed as INIT is
;           ignored.
;
; RETURNS:
;
;   The new array, which is always one-dimensional.  If COUNT is zero,
;   then the scalar -1L is returned.
;
; SEE ALSO:
;
;   STORE_ARRAY in IDL Astronomy Library
;
; MODIFICATION HISTORY:
;   Written, CM, 02 Mar 2000
;   Added OVERWRITE and EMPTY1 keyword, CM 04 Mar 2000
;
;  $Id: arrdelete.pro,v 1.2 2001/03/25 18:10:41 craigm Exp $
;
;-
; Copyright (C) 2000, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
;+
; NAME:
;   ARRINSERT
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   Insert one array into another
;
; CALLING SEQUENCE:
;   NEWARR = ARRINSERT(INIT, INSERT, [AT=POSITION] )
;
; DESCRIPTION: 
;
;   ARRINSERT inserts the contents of one array (INSERT) into
;   another (INIT), and returns the new array (NEWARR).
;
;   ARRINSERT will handle empty lists, which are represented as
;   undefined variables.  If both input arrays are empty, then the
;   scalar -1L is returned, and the keyword COUNT is set to 0L.
;
; INPUTS:
;
;   INIT - the initial array, into which INSERT will be inserted.  Any
;          data type, including structures, is allowed.  Regardless of
;          the dimensions of INIT, it is treated as a one-dimensional
;          array.  If OVERWRITE is not set, then INIT itself is
;          unmodified.
;
;   INSERT - the array to be inserted into INIT, which must be of the
;            same or similar type to INIT.  If INSERT is empty, then
;            INIT is returned unchanged.  Regardless of the dimensions
;            of INSERT, it is treated as a one-dimensional array.
;
; KEYWORDS:
;
;    AT - a long integer indicating the position of the newly inserted
;         sub-array.  If AT is non-negative, then INSERT will appear
;         at NEWARR[AT].  If AT is negative, then INSERT will appear
;         at NEWARR[AT + (N+1)] where N is the number of elements in
;         INIT, which is to say if AT is negative, it indexes from the
;         end side of the array rather than the beginning.
;
;         Default: 0L (INSERT appears at beginning of INIT)
;
;   OVERWRITE - if set, then the initial array INIT will be
;               overwritten by the new array.  Upon exit INIT becomes
;               undefined.
;
;   COUNT - upon return, the number of elements in the resulting
;           array.
;
;  EMPTY1, EMPTY2 - if set, then INIT (for EMPTY1) or INSERT (for
;                   EMPTY2) are assumed to be empty (i.e., to have
;                   zero elements).  The actual values passed as INIT
;                   or INSERT are then ignored.
;
; RETURNS:
;
;   The new array, which is always one-dimensional.  If COUNT is zero,
;   then the scalar -1L is returned.
;
; SEE ALSO:
;
;   ARRDELETE, STORE_ARRAY in IDL Astronomy Library
;
; MODIFICATION HISTORY:
;   Written, CM, 02 Mar 2000
;   Added OVERWRITE and EMPTY keywords, CM, 04 Mar 2000
;   Improved internal docs, and AT keyword docs, CM, 28 Sep 2000
;   Doc clarifications, CM, 29 Sep 2001
;
;  $Id: arrinsert.pro,v 1.3 2001/09/30 03:17:06 craigm Exp $
;
;-
; Copyright (C) 2000,2001, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
;+
; NAME:
;   CHEBCOEF
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Estimate Chebyshev polynomial coefficients of a function on an interval
;
; MAJOR TOPICS:
;   Curve and Surface Fitting
;
; CALLING SEQUENCE:
;   p = CHEBCOEF(FUNC, PRIVATE, FUNCTARGS=functargs, /DOUBLE, /EXPRESSION, $
;                PRECISION=prec, ERROR=err, NMAX=nmax, INTERVAL=interval, $
;                REDUCE_ALGORITHM=, STATUS=)
;
; DESCRIPTION:
;
;   CHEBCOEF estimates the coefficients for a finite sum of Chebyshev
;   polynomials approximating the function FUNC(x) over an interval.
;   The user can choose the desired precision and maximum number of
;   chebyshev coefficients.
;
;   This routine is intended for functions which can be evaluated to
;   full machine precision at arbitrary abcissae, and which are smooth
;   enough to ensure that the coefficients are a decreasing sequence.
;   For already-tabulated or potentially noisy data, the routines
;   CHEBGRID or CHEBFIT should be used instead.
;
;   The function to be approximated may either be the name of an IDL
;   function (the default behavior), or an IDL expression (using the
;   /EXPRESSION keyword).
;
;   The procedure uses a modified form of the classic algorithm for
;   determining the coefficients, which relies the orthogonality
;   relation for Chebyshev polynomials.  The interval [a,b] is
;   subdivided successively into sets of subintervals of length
;   2^(-k)*(b-a),(k = 0,1,2...). After each subdivision the
;   orthogonality properties of the Chebyshev polynomials with respect
;   to summation over equally-spaced points are used to compute two
;   sets of approximate values of the coefficients cj, one set
;   computed using the end-points of the subintervals, and one set
;   using the mid-points.  Certain convergence requirements must be
;   met before terminating.  If the routine fails to converge with 64
;   coefficents, then the current best-fitting coefficients are
;   returned, along with an error estimate in the ERROR keyword.
;   CHEBCOEF never returns more than 64 coefficients.
;
;   The coefficients may be further refined.  If the keyword
;   REDUCE_ALGORITHM is set to a value of 1, then any high order
;   coefficients below a certain threshold are discarded.  If
;   REDUCE_ALGORITHM is set to 2 (the default), then all coefficients
;   below the threshold are discarded rather than just the high order
;   ones.  The threshold is determined by the PRECISION keyword.
;
; INPUTS:
;
;   FUNC - a scalar string, the name of the function to be
;          approximated, or an IDL string containing an expression to
;          be approximated (if /EXPRESSION is set).
;
;  PRIVATE - any optional variable to be passed on to the function to
;            be integrated.  For functions, PRIVATE is passed as the
;            second positional parameter; for expressions, PRIVATE can
;            be referenced by the variable 'P'.  CHEBCOEF does not
;            examine or alter PRIVATE.
;
; RETURNS:
;
;   An array of Chebyshev coefficients which can be passed to
;   CHEBEVAL.  NOTE: the convention employed here is such that the
;   constant term in the expansion is P(0)*T0(x) (i.e., the convention
;   of Luke), and not P(0)/2 * T0(x).
;
; KEYWORD PARAMETERS:
;
;   DOUBLE - if set, then computations are done in double precision
;            rather than single precision.
;
;   ERROR - upon return, this keyword contains an estimate of the
;           maximum absolute error in the approximation.
;
;   EXPRESSION - if set, then FUNC is an IDL expression to be
;                approximated, rather than the name of a function.
;
;   FUNCTARGS - A structure which contains the parameters to be passed
;               to the user-supplied function specified by FUNCT via
;               the _EXTRA mechanism.  This is the way you can pass
;               additional data to your user-supplied function without
;               using common blocks.  By default, no extra parameters
;               are passed to the user-supplied function.
;
;   INTERVAL - a 2-element vector describing the interval over which
;              the polynomial is to be evaluated.
;              Default: [-1, 1]
;
;   NMAX - a scalar, the maximum number of coefficients to be
;          estimated.   This number may not exceed 64. 
;          Default: 64
;
;   PRECISION - a scalar, the requested precision in the
;               approximation.  Any terms which do not contribute
;               significantly, as defined by this threshold, are
;               discarded.  If the function to be estimated is not
;               well-behaved, then the precision is not guaranteed to
;               reach the desired level.  Default: 1E-7
;
;   REDUCE_ALGORITHM - a scalar integer, describes how insignificant
;               terms are removed from the fit.  If 0, then all terms
;               are kept, and none are dicarded.  If 1, then only
;               trailing terms less than PRECISION are discarded.  If
;               2, then both trailing and intermediate terms less than
;               PRECISION are discarded.
;               Default: 2
;
;   STATUS - upon return, this keyword contains information about the
;            status of the approximation.  A value of -1 indicates bad
;            input values; a value of 0 indicates the required
;            accuracy was not obtained; a value of 1 indicates
;            success.
;
; EXAMPLE:
;
;   x = dindgen(1000)/100     ; Range of 0 to 10
;   p = chebcoef('COS(x)', /expr, interval=[0d, 10d])  ;; Compute coefs
;   y = chebeval(x, p, interval=[0d,10d])              ;; Eval Cheby poly
;   plot, x, y - cos(x)       ; Plot residuals
;
; REFERENCES:
;
;   Abramowitz, M. & Stegun, I., 1965, *Handbook of Mathematical
;     Functions*, 1965, U.S. Government Printing Office, Washington,
;     D.C. (Applied Mathematical Series 55)
;   CERN, 1995, CERN Program Library, Function E406
;   Luke, Y. L., *The Special Functions and Their Approximations*,
;     1969, Academic Press, New York
;
; MODIFICATION HISTORY:
;   Written and documented, CM, June 2001
;   Copyright license terms changed, CM, 30 Dec 2001
;   Added usage message, CM, 20 Mar 2002
;   Changed docs slightly, CM, 25 Mar 2002
;
;  $Id: chebcoef.pro,v 1.6 2002/05/03 18:40:27 craigm Exp $
;
;-
; Copyright (C) 2001, 2002, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
;; Evaluate a user-supplied expression
function chebcoef_eval, x, p, expression=expr, _EXTRA=extra
  y = 0
  cmd = 'Y = '+expr
  dummy = execute(cmd)
  return, y
end

function chebcoef, f0, priv, functargs=fa, double=double, error=err, $
                   nmax=nmax, interval=interval, precision=prec0, $
                   expression=expr, reduce_algorithm=redalg0, $
                   status=status, indices=igood

  if n_params() EQ 0 then begin
      message, 'USAGE:', /info
      message, 'P = CHEBCOEF(FUNCT, [PRIV,] INTERVAL=[a,b], NMAX=...)', /info
      return, !values.d_nan
  endif

  sz = size(f0)
  err = -1
  if sz(sz(0)+1) NE 7 OR n_elements(f0) NE 1 then begin
      NO_FUNCT:
      message, 'ERROR: FUNCT must be a scalar string', /info
      return, 0
  endif

  ;; Check for empty string
  f = strtrim(f0(0),2)
  if f EQ '' then goto, NO_FUNCT

  ;; Prepare for EXPRESSION if requested
  if keyword_set(expr) then begin
      f = 'CHEBCOEF_EVAL'
      fa = {expression: strtrim(f0(0),2)}
  endif else begin
      f = strtrim(f0(0),2)
  endelse

  ;; Handle error conditions gracefully
  if NOT keyword_set(nocatch) then begin
      catch, catcherror
      if catcherror NE 0 then begin
          catch, /cancel
          message, 'Error detected while approximating '+f, /info
          message, !err_string, /info

          errmsg = 0
          if NOT keyword_set(expr) then begin
              f1 = byte(strupcase(strtrim(f0(0),2)))
              ca = (byte('A'))(0) 
              cz = (byte('Z'))(0) 
              c0 = (byte('0'))(0) 
              c9 = (byte('9'))(0) 
              c_ = (byte('_'))(0) 
              wh = where((f1 GE ca AND f1 LE cz) EQ 0 AND f1 NE c_ $
                         AND (f1 GE c0 AND f1 LE c9) EQ 0, ct)
              if ct GT 0 OR (f1(0) GE c0 AND f1(0) LE c9) then begin
                  message, ('FUNCT appears to be an expression.  Did you '+$
                            'intend to pass the /EXPRESSION keyword?'), /info
                  errmsg = 1
              endif
          endif
          if errmsg EQ 0 then $
            message, ('Please verify that function works and conforms to '+$
                      'the documentation'), /info
          ier = -1L
          return, 0L
      endif
  endif

  if n_elements(prec0) EQ 0 then prec = 1e-7 else prec = prec0(0)
  zero = prec*0.
  if keyword_set(double) then zero = 0D
  if n_elements(interval) LT 2 then interval = zero + [-1., 1.]
  if n_elements(redalg0) EQ 0 then redalg = 2 else redalg = floor(redalg0(0))
  status = -1
  
  a = interval(0)
  b = interval(1)
  hf = zero + 0.5
  eps = prec
  z1 = zero + 1
  z2 = zero + 2
  sz = size(zero)
  if sz(sz(0)+1) EQ 5 then pi = !dpi else pi = !pi

  x0 = [a, b]
  if n_elements(priv) GT 0 then begin
      if n_elements(fa) GT 0 then fv = call_function(f, x0, priv, _EXTRA=fa) $
      else                        fv = call_function(f, x0, priv)
  endif else begin
      if n_elements(fa) GT 0 then fv = call_function(f, x0, _EXTRA=fa) $
      else                        fv = call_function(f, x0)
  endelse
  

  ALFA=HF*(B-A)
  BETA=HF*(B+A)
  C1=fv(0)
  C2=fv(1)
  AC = [C2+C1, C2-C1]
  BC = AC*0

  for i = 1, 7 do begin
      I1=2^(I-1)
      I2=I1-1
      I3=2*I1
      C1=Z2/I1
      C2=PI/I1

      jj = dindgen(i2+1)
      x = alfa*cos((jj+hf)*c2)+beta
      if n_elements(priv) GT 0 then begin
          if n_elements(fa) GT 0 then fv = call_function(f, x, priv, $
                                                         _EXTRA=fa) $
          else                        fv = call_function(f, x, priv)
      endif else begin
          if n_elements(fa) GT 0 then fv = call_function(f, x, _EXTRA=fa) $
          else                        fv = call_function(f, x)
      endelse
      c = fv

      ;; Compute B-coefficients
      for j = 0L, i2 do begin
          F1=J*C2
          F2=-HF*F1
      
          C3=2*COS(F1)
          A2=zero
          A1=zero
          A0=C(I2)

          for K = I2-1,0L,-1 do begin
              A2=A1
              A1=A0
              A0=C(K)+C3*A1-A2
          endfor
          BC(J)=C1*(A0*COS(F1+F2)-A1*COS(F2))
          BC(I1)=zero
      endfor

      c = hf*[ac(0:i1-1)+bc(0:i1-1), rotate(ac(0:i1)-bc(0:i1),2)]
      cc = abs(c)
      cmx = max(cc)
      
      if (CMX GT 0) THEN begin
          CMX=1/CMX
          CC(I3)=HF*CC(I3)
          A0=CC(I2)*CMX
          A1=CC(I1)*CMX
          for J = I1+2,I3 do begin
              A2=CC(J)*CMX
              IF(A0 LE EPS AND A1 LE EPS AND A2 LE EPS) THEN $
                goto, CHEB9
              A0=A1
              A1=A2
          endfor
      ENDIF

      ;; DOUBLE THE NUMBER OF COEFFICIENTS.      
      if i LT 7 then begin
          ac = c(0:i3)
          bc = ac*0
      endif
  endfor

  ;; REQUIRED ACCURACY NOT OBTAINED

  NC=64
  DELTA=total(abs(c(60:nc)))
  message, 'WARNING: Required accuracy not obtained', /info
  status = 0
  goto, CLEANUP

  CHEB9:
  ;; REQUIRED ACCURACY OBTAINED
  ;; SUM NEGLECTED TERMS IN EXPANSION

  status = 1
  DELTA=total(cc(j:i3))

  ;;  CHECK IF FURTHER REDUCTION OF COEFFICIENTS IS POSSIBLE.

  NC=J-1
  REST=EPS-DELTA
  IF (REST GT 0) AND redalg GT 0 THEN begin
      while (CC(NC) LT REST) do begin
          DELTA=DELTA+CC(NC)
          REST=REST-CC(NC)
          NC=NC-1
      endwhile
  ENDIF

  CLEANUP:
  C(0)=HF*C(0)
  p = c(0:nc)
  rest = eps - delta

  if redalg EQ 2 then begin
      wh = where(cc(0:nc) LT prec, ct)
      i = ct-1
      while (i GE 0) AND ((rest GT 0) OR (status EQ 1)) do begin
          delta = delta + cc(wh(i))
          rest  = rest - cc(wh(i))
          p(wh(i)) = 0
          i = i - 1
      endwhile
  endif

  DONE:
  igood = where(p NE 0)
  err = delta
  RETURN, p
end

;+
; NAME:
;   CHEBEVAL
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Evaluate a Chebyshev polynomial on an interval, given the coefficients
;
; MAJOR TOPICS:
;   Curve and Surface Fitting
;
; CALLING SEQUENCE:
;   y = CHEBEVAL(X, P, INTERVAL=interval, DERIVATIVE=deriv)
;
; DESCRIPTION:
;
;   CHEBEVAL computes the values of a Chebyshev polynomial function at
;   specified abcissae, over the interval [a,b].  The user must supply
;   the abcissae and the polynomial coefficients.  The function is of
;   the form:
;
;               N
;       y(x) = Sum p_n T_n(x*)     x in [a,b]
;              i=0
;
;   Where T_n(x*) are the orthogonal Chebyshev polynomials of the
;   first kind, defined on the interval [-1,1] and p_n are the
;   coefficients.  The scaled variable x* is defined on the [-1,1]
;   interval such that (x*) = (2*x - a - b)/(b - a), and x is defined
;   on the [a,b] interval.
;
;   The derivative of the function may be computed simultaneously
;   using the DERIVATIVE keyword.  
;
;   The is some ambiguity about the definition of the first
;   coefficient, p_0, namely, the use of p_0 vs. the use of p_0/2.
;   The p_0 definition of Luke is used in this function.
;
; INPUTS:
;
;   X - a numerical scalar or vector, the abcissae at which to
;       evaluate the polynomial.  If INTERVAL is specified, then all
;       values of X must lie within the interval.
;
;   P - a vector, the Chebyshev polynomial coefficients, as returned
;       by CHEBFIT or CHEBCOEF.
;
; RETURNS:
;
;   An array of function values, evaluated at the abcissae.  The
;   numeric precision is the greater of X or P.
;
; KEYWORD PARAMETERS:
;
;   DERIVATIVE - upon return, a vector containing the derivative of
;                the function at each abcissa is returned in this
;                keyword.
;
;   INTERVAL - a 2-element vector describing the interval over which
;              the polynomial is to be evaluated.
;              Default: [-1, 1]
;
; EXAMPLE:
;
;   x = dindgen(1000)/100     ; Range of 0 to 10
;   p = chebcoef('COS(x)', /expr, interval=[0d, 10d])  ;; Compute coefs
;   y = chebeval(x, p, interval=[0d,10d])              ;; Eval Cheby poly
;   plot, x, y - cos(x)       ; Plot residuals
;
; REFERENCES:
;
;   Abramowitz, M. & Stegun, I., 1965, *Handbook of Mathematical
;     Functions*, 1965, U.S. Government Printing Office, Washington,
;     D.C. (Applied Mathematical Series 55)
;   CERN, 1995, CERN Program Library, Function E407
;   Luke, Y. L., *The Special Functions and Their Approximations*,
;     1969, Academic Press, New York
;
; MODIFICATION HISTORY:
;   Written and documented, CM, June 2001
;   Copyright license terms changed, CM, 30 Dec 2001
;   Added usage message, CM, 20 Mar 2002
;   Return a vector even when P has one element, CM, 22 Nov 2004
;   Fix bug in evaluation of derivatives, CM, 22 Nov 2004
;
;  $Id: chebeval.pro,v 1.6 2004/11/22 07:08:00 craigm Exp $
;
;-
; Copyright (C) 2001, 2002, 2004, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
;+
; NAME:
;   CHEBFIT
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Fit Chebyshev polynomial coefficients to a tabulated function
;
; MAJOR TOPICS:
;   Curve and Surface Fitting
;
; CALLING SEQUENCE:
;   p = CHEBFIT(X, Y, ERR, INTERVAL=interval, NMAX=nmax,
;               PRECISION=prec, /EVEN, /ODD, REDUCE_ALGORITHM=)
;
; DESCRIPTION:
;
;   CHEBFIT fits a series of Chebyshev polynomials to a set of
;   tabulated and possibly noisy data points.  The functions MPFIT and
;   CHEBEVAL, available from the above web page, must also be in your
;   IDL path for this function to work properly.  The user can choose
;   the desired precision and maximum number of chebyshev
;   coefficients.
;
;   This function is intended for use on already-tabulated data which
;   are potentially noisy.  The user should never expect more than
;   NPOINTS terms, where NPOINTS is the number of (x,y) pairs.  For
;   functions which can be evaluated to full machine precision at
;   arbitrary abcissae, the routine CHEBCOEF should be used instead.
;   For exact data tabulated on a regular grid, the routine CHEBGRID
;   should be tried.
;
;   The user can also specify that the function is even or odd, using
;   the keywords EVEN or ODD.  This saves computation time because
;   certain terms in the expansion can be ignored.  For the purposes
;   of this function even and odd refer to the symmetry about the
;   center of the interval.
;
;   The algorithm is employed in three steps.  In the first step, the
;   coefficients are estimated at a crude level.  In the second step,
;   it is determined whether certain coefficients are deemed
;   "ignoreable", i.e., they do not contribute significantly to the
;   function and are discarded.  The operation of this step is
;   determined by the REDUCE_ALGORITHM keyword.  Finally, the
;   remaining "good" coefficients are re-fitted to achieve the best
;   fit.
;
; INPUTS:
;
;   X, Y - the x- and y- tabulated values to be fitted.
;
;   ERR - (optional) the y-error bar associated with each (x,y) pair.
;         Default: 1
;
; RETURNS:
;
;   An array of Chebyshev coefficients which can be passed to
;   CHEBEVAL.  NOTE: the convention employed here is such that the
;   constant term in the expansion is P(0)*T0(x) (i.e., the convention
;   of Luke), and not P(0)/2 * T0(x).
;
; KEYWORD PARAMETERS:
;
;   EVEN, ODD - if set, then the fitting routine assumes the function
;               is even or odd, about the center of the interval.
;
;   INTERVAL - a 2-element vector describing the interval over which
;              the polynomial is to be evaluated.
;              Default: [-1, 1]
;
;   NMAX - a scalar, the maximum number of polynomial terms to be
;          fitted at one time.
;          Default: 16
;
;   PRECISION - a scalar, the requested precision in the fit.  Any
;               terms which do not contribute significantly, as
;               defined by this threshold, are discarded.  If the
;               function to be fitted is not well-behaved, then the
;               precision is not guaranteed to reach the desired
;               level.
;               Default: 1E-7
;
;   REDUCE_ALGORITHM - a scalar integer, describes how insignificant
;               terms are removed from the fit.  If 0, then all terms
;               are kept, and none are dicarded.  If 1, then only
;               trailing terms less than PRECISION are discarded.  If
;               2, then both trailing and intermediate terms less than
;               PRECISION are discarded.
;               Default: 2
;
; EXAMPLE:
;
;   x = dindgen(1000)/100     ; Range of 0 to 10
;   y = cos(x) + randomn(seed,1000)*0.01  ; Function with some noise
;   p = chebfit(x, y, interval=[0d,10d])
;   plot, x, y - chebeval(x,p, interval=[0d,10d])
;
; REFERENCES:
;
;   Abramowitz, M. & Stegun, I., 1965, *Handbook of Mathematical
;     Functions*, 1965, U.S. Government Printing Office, Washington,
;     D.C. (Applied Mathematical Series 55)
;   CERN, 1995, CERN Program Library, Function E407
;   Luke, Y. L., *The Special Functions and Their Approximations*,
;     1969, Academic Press, New York
;
; MODIFICATION HISTORY:
;   Written and documented, CM, June 2001
;   Copyright license terms changed, CM, 30 Dec 2001
;   Added usage message, CM, 20 Mar 2002
;   Slight docs change, CM, 25 Mar 2002
;
;  $Id: chebfit.pro,v 1.7 2003/07/20 05:53:44 craigm Exp $
;
;-
; Copyright (C) 2001, 2002, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
;; Compute residuals for MPFIT
function chebfit_eval, p, interval=interval, nterms=nterms, igood=igood, $
                       _EXTRA=extra

  common chebfit_common, x, y, err

  if n_elements(igood) EQ 0 then begin
      p1 = p
  endif else begin
      p1 = replicate(p(0)*0, nterms)
      p1(igood) = p
  endelse

  ;; Compute the Chebyshev polynomial
  f = chebeval(x, p1, interval=interval)

  ;; Compute the deviates, applying either errors or weights
  if n_elements(err) GT 0 then begin
      result = (y-f)/err
  endif else if n_elements(wts) GT 0 then begin
      result = (y-f)*wts
  endif else begin
      result = (y-f)
  endelse
  
  ;; Make sure the returned result is one-dimensional.
  result = reform(result, n_elements(result), /overwrite)
  return, result
end

function chebfit, x, y, err, nmax=nterms0, interval=interval, $
                  precision=prec, even=even, odd=odd, quiet=quiet, $
                  initialize=init, reduce_algorithm=redalg0, $
                  indices=igood, nocatch=nocatch, $
                  yfit=yfit, perror=perror, bestnorm=bestnorm, dof=dof

  if n_params() EQ 0 then begin
      message, 'USAGE:', /info
      message, 'P = CHEBFIT(X, Y, ERR, INTERVAL=[a,b], NMAX=, ...)', /info
      return, !values.d_nan
  endif

  if n_elements(nterms0) EQ 0 then nterms = 16L $
  else                             nterms = floor(nterms0(0)) > 2L
  nterms = nterms < n_elements(x)
  if n_elements(interval) LT 2 then interval = [-1., 1.]
  if n_elements(prec) EQ 0 then prec = 1.e-7
  if n_elements(redalg0) EQ 0 then redalg = 2 else redalg = floor(redalg0(0))
  if n_elements(quiet) EQ 0 then quiet = 1

  ;; Handle error conditions gracefully
  if NOT keyword_set(nocatch) then begin
      catch, catcherror
      if catcherror NE 0 then begin
          catch, /cancel
          message, 'Error detected while fitting', /info
          message, !err_string, /info
          ier = -1L
          return, 0L
      endif
  endif

  if n_elements(p) LT nterms OR keyword_set(init) then begin
      p = replicate(x(0)*0 + 1, nterms) / (findgen(nterms)+1)^2
      p(0) = total(y)/n_elements(y)
      ;; If mean is *exactly* zero, then shift it off slightly
      if p(0) EQ 0 then p(0) = sqrt(total(y^2))/n_elements(y)/10
  endif
  p0 = p
  igood = lindgen(nterms)

  if keyword_set(even) OR keyword_set(odd) then $
    igood = lindgen(n_elements(p)/2)*2 + keyword_set(odd)
  nt = min([nterms, max(igood)+1])

  ;; Cancel out old common entries
  common chebfit_common, xc, yc, errc
  xc   = 0 & dummy = temporary(xc)
  yc   = 0 & dummy = temporary(yc)
  errc = 0 & dummy = temporary(errc)
  
  xc = x
  yc = y
  if n_elements(err) GT 0 then begin
      errc = err
  endif

  fa = {interval: interval, igood: igood, nterms: nt}
  p1 = mpfit('CHEBFIT_EVAL', p0(igood), functargs=fa, maxiter=5, quiet=quiet)
  p0(igood) = p1

  ;; Look for and remove the insignificant terms from the fit
  if redalg GT 0 then begin
      wh = where(abs(p1) GT prec(0), ct)
      if ct EQ 0 then begin
          ALL_ZERO:
          message, 'WARNING: no significant Chebyshev terms were detected', $
            /info
          p = p0*0
          return, 0L
      endif
      if max(wh) LT n_elements(igood)-1 then begin
          imax = max(wh)
          igood = igood(0:imax)
          p1 = p1(0:imax)
      endif

      if redalg EQ 2 then begin
          wh = where(abs(p1) GT 0.1*prec, ct)
          if ct EQ 0 then goto, ALL_ZERO
          igood = igood(wh)
          p1 = p1(wh)
      endif
  endif

  nt = min([nterms, max(igood)+1])
  fa = {interval: interval, igood: igood, nterms: nt}
  p2 = mpfit('CHEBFIT_EVAL', p1, functargs=fa, maxiter=10, quiet=quiet, $
             perror=dp2, bestnorm=bestnorm, dof=dof)

  xc = 0 & yc = 0 & errc = 0
  p = p0*0 
  perror = p
  p(igood) = p2
  perror(igood) = dp2

  if arg_present(yfit) then $
    yfit = chebeval(x, p, interval=interval)

  return, p
end
;+
; NAME:
;   CHEBGRID
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Estimate Chebyshev polynomial coefficients of a function on a grid
;
; MAJOR TOPICS:
;   Curve and Surface Fitting
;
; CALLING SEQUENCE:
;   p = CHEBGRID(T, X, [ DXDT, NPOINTS=, NPOLY=, NGRANULE= , $
;                RMS=, DRMS=, RESIDUALS=, DRESIDUALS= , $
;                XMATRIX=, DXMATRIX=, RESET=, 
;                DERIV_WEIGHT= ] )
;
; DESCRIPTION:
;
;   CHEBGRID estimates the coefficients for a finite sum of Chebyshev
;   polynomials approximating a continuous tabulated function over an
;   interval.  The function (and optionally its derivative) must be
;   tabulated on a regularly sampled grid.  The implementation of this
;   function is taken from a method described by X. X. Newhall, used
;   in estimating coefficients for ephemerides in the solar system.
;
;   The tabulated function is assumed to be continuous over the entire
;   interval.  A Chebyshev series is fitted to the function over small
;   segments, called granules.  The size of each granule, the number
;   of points in each granule, and the number of Chebyshev polynomials
;   are all configurable.
;
;   Users may specify either the function alone, or the function and
;   its first derivative.  By also giving the tabulated derivative, a
;   more accurate Chebyshev polynomial can be developed.  Aside from
;   the constraints mentioned in the next paragraph, the polynomial
;   that is returned is the best-fit polynomial in a least-squares
;   sense.
;
;   Here is a definition of terms:
;
;      GRANULE - a single continuous fitted segment.  The length of the
;         granule, NGRANULE, is specified in units of the tabulated
;         grid size.  Because of the continuity requirements developed
;         below, granules will always overlap at their endpoints.
;         Thus, then length of a granule should be a factor of
;         N_ELEMENTS(X)-1.  For simple functions over short intervals,
;         the granule size can be equal to N_ELEMENTS(X)-1
;
;      NUMBER OF POINTS the number of points, NPOINTS, within a
;         granule to be fitted to the polynomial, not necessarily
;         equal to the granule size.  The greater the number of
;         points, the more computation time and storage is required.
;         This number *must* be a factor of NGRANULE.  Typically
;         NPOINTS is a number between 8 and 12.  Because of the
;         single-point overlap between granules (see below), the
;         actual number of points per fit is NPOINTS+1.
;
;      NUMBER OF POLYNOMIALS the number of Chebyshev polynomial terms,
;         NPOLYNOMIAL, to be fitted per granule.  The greater the
;         number of polynomial terms, the more computation time and
;         storage is required, but also the greater the approximating
;         precision of the fit.
;
;   The particular set of Chebyshev polynomial coefficients developed
;   by this function have some special properties.  If both the
;   function and its derivative are specified, then the value and
;   derivative of the interpolating polynomial at the granule
;   endpoints will be exactly equal to the tabulated endpoint values.
;   This feature allows many approximations to be strung together
;   piecewise, and the function value and first derivative will be
;   continuous across granule boundaries.
;
;   If only the function value is specified, then only the function
;   value will be continuous at the granule endpoints, and not the
;   derivative.
;
;   An extensive set of statistics are computed to assess the quality
;   of the Chebyshev polynomial fit.  The keywords RESIDUALS and
;   DRESIDUALS return the residuals of the fit after subtracting the
;   interpolation.  The RMS and DRMS keywords return the root mean
;   squared deviations between data and model.
;
;   If the user does not know how many granules, points, or polynomial
;   coefficients to use, then he or she should try several
;   combinations and see which minimizes the r.m.s. value with the
;   fewest number of coefficients.
;
;   If the XMATRIX and DXMATRIX keywords are passed, then CHEBGRID
;   attempts to avoid recomputing several of the matrices it uses in
;   estimating the coefficients.  If multiple calls to CHEBGRID are to
;   be made, some compution time savings can be made.  In the first
;   call CHEBGRID the required matrices are computed and returned.  In
;   subsequent calls, CHEBGRID detects the XMATRIX and DXMATRIX
;   keyword values and uses those values if it can.
;
;   The user can also estimate their own coefficients.  The matrices
;   returned are (NPOINTS+1)x(NPOLYNOMIAL).  The coefficients from a
;   NPOINTS+1 tabulation, X, are found by:
;
;      PCHEB = XMATRIX ## X  +  DXMATRIX ## DXDT
;
;   if derivative information is known, or
;
;      PCHEB = XMATRIX ## X
;
;   if no derivative information is known.  [ Note: the matrices are
;   different, depending on whether derivative information is known or
;   not. ]
;
;   
; INPUTS:
;
;   T - array of regularly sampled *independent* variables.  The number
;       of elements in T should be a multiple of NGRANULE, plus one.
;
;   X - array of regularly sampled *dependent* variables.  The number
;       of elements in X should be equal to the number of elements in
;       T.
;
;   DXDT - optionally, a tabulated array of first derivatives of X
;          with respect to T, at the same grid points.
;
; KEYWORD PARAMETERS:
;
;   NGRANULE - size of a "granule", in grid intervals.  NGRANULE must
;              be at least 2, and a factor of N_ELEMENTS(T)-1.
;              Default: 8
;
;   NPOINTS - number of points per granule that are fitted.  NPOINTS
;             must be at least 2, and a factor of NGRANULE. 
;             Default: NGRANULE
;
;   NPOLYNOMIAL - number of Chebyshev polynomial terms per fit.
;                 NPOLYNOMIAL must be at least 2 and less than
;                 2*(NPOINTS+1), when derivative information is
;                 specified; or less than NPOINTS+1, when no
;                 derivative information is specified.
;                 Default: 7
;
;  RESIDUALS - upon return, an array of size N_ELEMENTS(T), with
;              residuals of the tabulated function minus the
;              interpolated function.
;
;  DRESIDUALS - same as RESIDUALS, but for the function's first
;               derivative.
;
;  RMS - upon return, the root mean square of the function value
;        residuals.
;
;  DRMS - same as RMS, but for the function's first derivative.
;
;  XMATRIX - upon return, the matrix used to compute Chebyshev
;            polynomial coefficients from the function value.
;
;            Upon input, CHEBGRID determines if XMATRIX will apply to
;            the data, and if so, XMATRIX is reused rather than
;            computed.  If XMATRIX cannot be reused, then it is
;            computed afresh, and the new value is returned in the
;            XMATRIX keyword.
;
;            The user should not modify the contents of this array.
;
;  DXMATRIX - same as XMATRIX, but for the function's first
;             derivative.
;
;  RESET - if set, force a recomputation of XMATRIX and/or DXMATRIX.
;
;  DERIV_WEIGHT - amount of weight to give to function derivative,
;                 relative to the function value.
;                 Default: 0.16d
;
;
; RETURNS:
;
;   An array of coefficient values.  The dimensions of the array are
;   NPOLYNOMIALxNSEGS, where NSEGS is the number of granules in the
;   entire interval.
;
;
; EXAMPLE:
;
;   ;; Estimate Chebyshev coefficients for the function SIN(X), on the
;   ;; interval [-1,+1].  
;   xx = dindgen(9)/4d - 1d   ;; Regular grid from -1 to 1 (9 points)
;   yy = sin(xx)              ;; Function values, sin(x), ...
;   dy = cos(xx)              ;; ... and derivatives
;
;   ;; Estimate coefficients using CHEBGRID (single granule of 8 intervals)
;   p = chebgrid(xx, yy, dy, npoints=8, ngranule=8, npoly=10)
;
;   xxx = dindgen(1001)/500 - 1d   ;; New grid for testing
;   res = sin(xxx) - chebeval(xxx, p)
;   plot, xxx, res
;
;   ;; Same as example above, except extended range to [-1, +15],
;   using eight granules.
;   xx2 = dindgen(65)/4d - 1
;   yy2 = sin(xx2)
;   dy2 = cos(xx2)
;   p = chebgrid(xx2, yy2, dy2, ngranule=8, npoint=8, npoly=10)
;   help, p
;       P               DOUBLE    = Array[10, 8]
;   ;; (i.e., 10 polynomial coefficients over 8 granules)
;
;
; REFERENCES:
;
;   Abramowitz, M. & Stegun, I., 1965, *Handbook of Mathematical
;     Functions*, 1965, U.S. Government Printing Office, Washington,
;     D.C. (Applied Mathematical Series 55)
;   Newhall, X. X. 1989, Celestial Mechanics, 45, p. 305-310
;
; MODIFICATION HISTORY:
;   Written, CM, Feb 2002
;   Documented, CM, 24 Mar 2002
;   Corrected documentation, CM, 28 Apr 2002
;   Typo correction, CM, 10 Oct 2002
;
;  $Id: chebgrid.pro,v 1.5 2002/11/07 00:15:23 craigm Exp $
;
;-
; Copyright (C) 2002, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
;; Utility function: compute XMATRIX and DXMATRIX using Newhall approach
pro chebpcmat, npts, npoly, xmat, vmat, dweight=weight0

  ;; n0 is the number of intervals in Cheb approx.
  n0 = npts - 1
  if n_elements(weight0) EQ 0 then $
    weight = 0.16d $
  else $
    weight = weight0(0)

  tmat = dblarr(npoly, npts)
  tdot = tmat

  cj = dblarr(npoly)
  xj = 1d - 2d*dindgen(npts)/n0
  for i = 0, npoly-1 do begin
      cj(*) = 0 & cj(i) = 1
      tmat(i,*) = chebeval(xj, cj, deriv=v)
      tdot(i,*) = v
  endfor

  ;; Form matrix T*W
  tw = dblarr(2,npts,npoly)
  tw(0,*,*) = transpose(tmat)
  tw(1,*,*) = transpose(tdot) * weight

  ;; Form matrix T*WT
  twt = reform(tw(0,*,*),npts,npoly) ## tmat + $
    reform(tw(1,*,*),npts,npoly) ## tdot

  tw  = reform(tw, 2*npts, npoly, /overwrite)
  twt = reform(twt, npoly, npoly, /overwrite)

  ;; Augment matrix T*W to get matrix C2
  c2 = dblarr(2*npts,npoly+4)
  c2(*,0:npoly-1) = tw
  c2(0,npoly)          = 1 & c2(1,npoly+1)        = 1
  c2(2*npts-2,npoly+2) = 1 & c2(2*npts-1,npoly+3) = 1

  ;; Augment matrix T*WT to get the matrix C1
  c1 = dblarr(npoly+4,npoly+4)
  c1(0:npoly-1,0:npoly-1) = twt

  c1(0:npoly-1,npoly+0) = tmat(*,0)
  c1(0:npoly-1,npoly+1) = tdot(*,0)
  c1(0:npoly-1,npoly+2) = tmat(*,npts-1)
  c1(0:npoly-1,npoly+3) = tdot(*,npts-1)

  c1(npoly:*,0:npoly-1) = transpose(c1(0:npoly-1,npoly:*))

  ;; Compute matrix C1^(-1)
  c1inv = invert(c1)
  ;; Compute matrix C1^(-1) C2
  c1c2 = c1inv ## c2

  c1c2 = reform(c1c2, 2,npts,npoly+4)
  c1c2 = reverse(c1c2,2)
  c1c2 = reform(c1c2, 2*npts,npoly+4)

  ii = lindgen(npts)*2
  xmat = c1c2(ii,0:npoly-1)   ;; Split into terms multiplying Y and VY
  vmat = c1c2(ii+1,0:npoly-1)

  return
end

;; Utility function: compute XMATRIX only, using only the constraint
;; on the function values at the endpoints.
pro chebpcmat_xonly, npts, npoly, xmat

  ;; n0 is the number of points in Cheb approx.
  n0 = npts - 1
  
  tmat = dblarr(npoly, npts)

  cj = dblarr(npoly)
  xj = 1d - 2d*dindgen(npts)/n0
  for i = 0, npoly-1 do begin
      cj(*) = 0 & cj(i) = 1
      tmat(i,*) = chebeval(xj, cj, deriv=v)
  endfor

  ;; Augment matrix T to get matrix C2
  c2 = dblarr(npts,npoly+2)
  c2(*,0:npoly-1) = transpose(tmat)
  c2(0,npoly)     = 1
  c2(npts-1,npoly+1) = 1

  ;; Augment matrix T*WT to get the matrix C1
  c1 = dblarr(npoly+2,npoly+2)
  c1(0:npoly-1,0:npoly-1) = transpose(tmat) ## tmat

  c1(0:npoly-1,npoly+0) = tmat(*,0)
  c1(0:npoly-1,npoly+1) = tmat(*,npts-1)

  c1(npoly:*,0:npoly-1) = transpose(c1(0:npoly-1,npoly:*))

  ;; Compute matrix C1^(-1)
  c1inv = invert(c1)
  ;; Compute matrix C1^(-1) C2
  c1c2 = c1inv ## c2

  c1c2 = reform(c1c2, npts,npoly+2)
  c1c2 = reverse(c1c2,1)
  xmat = c1c2(*,0:npoly-1)

  return
end

function chebgrid, t, x, dxdt, ngranule=ngran0, npoints=npts0, $
                   npolynomial=npoly0, deriv_weight=dweight0, $
                   rms=rms, drms=drms, residuals=resid, dresiduals=dresid, $
                   xmatrix=xmatrix, dxmatrix=dxmatrix, reset=reset

  ;; Default processing
  if n_elements(ngran0) EQ 0 then ngran = 8 $
  else ngran = round(ngran0(0)) > 2

  if n_elements(npts0) EQ 0 then npts = ngran $
  else npts = round(npts0(0)) > 2

  if n_elements(npoly0) EQ 0 then npoly = 7 $
  else npoly = round(npoly0(0)) > 2

  ;; Error checking
  if ngran LT npts then begin
      message, 'ERROR: Granule size ('+strtrim(ngran,2)+') is too '+ $
        'small for number of samples ('+strtrim(npts,2)+')'
      return, !values.d_nan
  endif
  
  ;; Be sure NGRAN is a multiple of NPTS - or not.  Instead, a warning
  ;; message is printed in the loop.
;  if abs(double(ngran)/npts - round(ngran/npts)) GT 1d-5 then begin
;      message, 'ERROR: NPOINTS must be a multiple of NGRANULE'
;      return, !values.d_nan
;  endif
  
  ;; Be sure we are solving a least-squares problem.  If the number of
  ;; polynomials is too great then it becomes underconstrained, not
  ;; overconstrained.
  if n_elements(dxdt) GT 0 then begin
      if npoly GE 2*(npts+1) then $
        message, 'ERROR: NPOLYNOMIAL must be less than 2*(NPOINTS+1)'
  endif else begin
      if npoly GE npts+1 then $
        message, 'ERROR: NPOLYNOMIAL must be less than NPOINTS+1'
  endelse
  
  ;; Begin size checking of input matrices - we may be able to use the
  ;; previously computed version.
  szx = size(xmatrix)
  szv = size(dxmatrix)

  ;; Cases: recompute because existing X matrix is wrong size;
  ;;        recompute because existing V matrix is wrong size;
  ;;        recompute because a V matrix was passed, but no DXDT was
  redo_x = (szx(0) NE 2 OR szx(1) NE npts+1 OR szx(2) NE npoly)
  redo_v = (n_elements(dxdt) GT 0 AND $
            (szv(0) NE 2 OR szv(1) NE npts+1 OR szv(2) NE npoly))
  no_v = (n_elements(dxdt) EQ 0 AND n_elements(dxmatrix) GT 0)

  ;; Actual recomputation of matrices
  if redo_x OR redo_v OR no_v OR keyword_set(reset) then begin
      COMPUTE_CHEBMAT:
      xmatrix = 0 & dummy = temporary(xmatrix)
      dxmatrix = 0 & dummy = temporary(dxmatrix)

      if n_elements(dxdt) GT 0 then $
        chebpcmat, npts+1, npoly, xmatrix, dxmatrix, dweight=dweight0 $
      else $
        chebpcmat_xonly, npts+1, npoly, xmatrix
  endif

  rms = 0.*x(0)
  drms = rms
  chebm = dblarr(npoly, (n_elements(x)-1)/ngran)
  resid = x*0.
  dresid = resid

  ispan = lindgen(npts+1)*(ngran/npts)
  imax = max(ispan)
  ng = 0L
  for ibase = 0, n_elements(x)-1, ngran do begin
      if ibase EQ n_elements(x)-1 then goto, DONE
      if n_elements(x)-ibase LT ngran+1 then begin
          nlost = n_elements(x)-ibase
          message, 'WARNING: last '+strtrim(nlost,2)+' elements of X '+$
            'were discarded because they formed only a fractional granule.', $
            /info
          goto, DONE
      endif

      tspan  = [t(ibase), t(ibase+imax)]
      tgran = t(ibase:ibase+imax-1)-t(ibase)
      dt     = tspan(1) - tspan(0)
      tspan  = tspan - tspan(0)

      ;; Compute the X portion of the coefficients
      xgran  = x(ibase+ispan)
      chebi  = xmatrix ## xgran

      ;; Compute the DXDT portion if it is available
      if n_elements(dxdt) GT 0 then begin
          dxgran = dxdt(ibase+ispan) * dt/2.
          chebi = chebi + dxmatrix ## dxgran

          ;; Statistics - V first, then X comes later
          xmod = chebeval(tgran, chebi, interval=tspan, derivative=dxmod)

          ;; DXDT portion of statistics
          dresid(ibase:ibase+imax-1) = dxdt(ibase:ibase+imax-1) -  dxmod
          diff_dx = (dxdt(ibase:ibase+imax-1) - dxmod)^2
          drms = drms + total(diff_dx)

      endif else begin
          ;; Statistics - X only

          xmod = chebeval(tgran, chebi, interval=tspan)
      endelse

      ;; Finish statistics with X portion
      resid(ibase:ibase+imax-1) = x(ibase:ibase+imax-1) -  xmod
      diff_x  = (   x(ibase:ibase+imax-1) -  xmod)^2
      rms  =  rms + total(diff_x)

      ;; Append to existing coefficient list
      chebm(*,ng) = chebi(*)
      ng = ng + 1L
  endfor
  
  DONE:
  ;; Final adjustments to statistics
  rms  = sqrt( rms / ngran)
  if n_elements(dxdt) GT 0 then drms = sqrt(drms / ngran)

  return, chebm
end
;+
; NAME:
;   CMAPPLY
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   Applies a function to specified dimensions of an array
;
; MAJOR TOPICS:
;   Arrays
;
; CALLING SEQUENCE:
;   XX = CMAPPLY(OP, ARRAY, DIMS, [/DOUBLE], [TYPE=TYPE])
;
; DESCRIPTION:
;   CMAPPLY will apply one of a few select functions to specified 
;   dimensions of an array.  Unlike some IDL functions, you *do* have
;   a choice of which dimensions that are to be "collapsed" by this
;   function.  Iterative loops are avoided where possible, for 
;   performance reasons.
;
;   The possible functions are:             (and number of loop iterations:)
;     +     - Performs a sum (as in TOTAL)       number of collapsed dimensions
;     AND   - Finds LOGICAL "AND" (not bitwise)  same
;     OR    - Finds LOGICAL "OR"  (not bitwise)  same
;     *     - Performs a product                 LOG_2[no. of collapsed elts.]
;
;     MIN   - Finds the minimum value            number of collapsed dimensions
;     MAX   - Finds the maximum value            same
;     MEDIAN- Finds the median  value            same
;
;     USER  - Applies user-defined function      no. of output elements
;
;
;   It is possible to perform user-defined operations arrays using
;   CMAPPLY.  The OP parameter is set to 'USER:FUNCTNAME', where
;   FUNCTNAME is the name of a user-defined function.  The user
;   defined function should be defined such that it accepts a single
;   parameter, a vector, and returns a single scalar value.  Here is a
;   prototype for the function definition:
;
;      FUNCTION FUNCTNAME, x, KEYWORD1=key1, ...
;         scalar = ... function of x or keywords ...
;         RETURN, scalar
;      END
;
;   The function may accept keywords.  Keyword values are passed in to
;   CMAPPLY through the FUNCTARGS keywords parameter, and passed to
;   the user function via the _EXTRA mechanism.  Thus, while the
;   definition of the user function is highly constrained in the
;   number of positional parameters, there is absolute freedom in
;   passing keyword parameters.
;
;   It's worth noting however, that the implementation of user-defined
;   functions is not particularly optimized for speed.  Users are
;   encouraged to implement their own array if the number of output
;   elements is large.
;
;
; INPUTS:
;
;   OP - The operation to perform, as a string.  May be upper or lower
;        case.
;
;        If a user-defined operation is to be passed, then OP is of
;        the form, 'USER:FUNCTNAME', where FUNCTNAME is the name of
;        the user-defined function.
;
;   ARRAY - An array of values to be operated on.  Must not be of type
;           STRING (7) or STRUCTURE (8).
;
; OPTIONAL INPUTS:
;
;   DIMS - An array of dimensions that are to be "collapsed", where
;          the the first dimension starts with 1 (ie, same convention
;          as IDL function TOTAL).  Whereas TOTAL only allows one
;          dimension to be added, you can specify multiple dimensions
;          to CMAPPLY.  Order does not matter, since all operations
;          are associative and transitive.  NOTE: the dimensions refer
;          to the *input* array, not the output array.  IDL allows a 
;          maximum of 8 dimensions.
;          DEFAULT: 1 (ie, first dimension)
;
; KEYWORDS:
;
;   DOUBLE - Set this if you wish the internal computations to be done
;            in double precision if necessary.  If ARRAY is double 
;            precision (real or complex) then DOUBLE=1 is implied.
;            DEFAULT: not set
;
;   TYPE - Set this to the IDL code of the desired output type (refer
;          to documentation of SIZE()).  Internal results will be
;          rounded to the nearest integer if the output type is an
;          integer type.
;          DEFAULT: same is input type
;
;   FUNCTARGS - If OP is 'USER:...', then the contents of this keyword
;               are passed to the user function using the _EXTRA
;               mechanism.  This way you can pass additional data to
;               your user-supplied function, via keywords, without
;               using common blocks.
;               DEFAULT: undefined (i.e., no keywords passed by _EXTRA)
;
; RETURN VALUE:
;
;   An array of the required TYPE, whose elements are the result of
;   the requested operation.  Depending on the operation and number of
;   elements in the input array, the result may be vulnerable to
;   overflow or underflow.
;
; EXAMPLES:
;   Shows how CMAPPLY can be used to total the second dimension of the
;   array called IN.  This is equivalent to OUT = TOTAL(IN, 2)
;
;   IDL> IN  = INDGEN(5,5)
;   IDL> OUT = CMAPPLY('+', IN, [2])
;   IDL> HELP, OUT
;   OUT             INT       = Array[5]
;
;   Second example.  Input is assumed to be an 5x100 array of 1's and
;   0's indicating the status of 5 detectors at 100 points in time.
;   The desired output is an array of 100 values, indicating whether
;   all 5 detectors are on (=1) at one time.  Use the logical AND
;   operation.
;
;   IDL> IN = detector_status             ; 5x100 array
;   IDL> OUT = CMAPPLY('AND', IN, [1])    ; collapses 1st dimension
;   IDL> HELP, OUT
;   OUT             BYTE      = Array[100]
;
;   (note that MIN could also have been used in this particular case,
;   although there would have been more loop iterations).
;
;   Third example.  Shows sum over first and third dimensions in an
;   array with dimensions 4x4x4:
;
;   IDL> IN = INDGEN(4,4,4)
;   IDL> OUT = CMAPPLY('+', IN, [1,3])
;   IDL> PRINT, OUT
;        408     472     536     600
;
;   Fourth example. A user-function (MEDIAN) is used:
;
;   IDL> IN = RANDOMN(SEED,10,10,5)
;   IDL> OUT = CMAPPLY('USER:MEDIAN', IN, 3)
;   IDL> HELP, OUT
;   OUT             FLOAT     = Array[10, 10]
;
;   (OUT(i,j) is the median value of IN(i,j,*))
;
; MODIFICATION HISTORY:
;   Mar 1998, Written, CM
;   Changed usage message to not bomb, 24 Mar 2000, CM
;   Signficant rewrite for *, MIN and MAX (inspired by Todd Clements
;     <Todd_Clements@alumni.hmc.edu>); FOR loop indices are now type
;     LONG; copying terms are liberalized, CM, 22, Aug 2000
;   More efficient MAX/MIN (inspired by Alex Schuster), CM, 25 Jan
;     2002
;   Make new MAX/MIN actually work with 3d arrays, CM, 08 Feb 2002
;   Add user-defined functions, ON_ERROR, CM, 09 Feb 2002
;   Correct bug in MAX/MIN initialization of RESULT, CM, 05 Dec 2002
;   Correct bug in CMAPPLY_PRODUCT implementation when there are an
;     odd number of values to combine, CM 26 Jul 2006
;   Add fast IDL versions of '*', 'MEDIAN', 'MIN' and 'MAX', where
;     IDL supports it, CM 26 Jul 2006
;
;  $Id: cmapply.pro,v 1.6 2006/07/26 19:34:24 craigm Exp $
;
;-
; Copyright (C) 1998, 2000, 2002, 2006, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
;; Utility function, adapted from CMPRODUCT
function cmapply_product, x
  sz = size(x)
  n = sz(1)

  while n GT 1 do begin
      if (n mod 2) EQ 1 then x(0,*) = x(0,*) * x(n-1,*)
      n2 = floor(n/2)
      x = x(0:n2-1,*) * x(n2:n2*2-1,*)
      n = n2
  endwhile
  return, reform(x(0,*), /overwrite)
end

;; Utility function, used to collect collaped dimensions
pro cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep
  sz = size(newarr)
  ;; First task: rearrange dimensions so that the dimensions
  ;; that are "kept" (ie, uncollapsed) are at the back
  dimkeep = where(histogram(dimapply,min=1,max=sz(0)) ne 1, nkeep)
  if nkeep EQ 0 then return

  newarr = transpose(temporary(newarr), [dimapply-1, dimkeep])
  ;; totcol is the total number of collapsed elements
  totcol = sz(dimapply(0))
  for i = 1, n_elements(dimapply)-1 do totcol = totcol * sz(dimapply(i))
  totkeep = sz(dimkeep(0)+1)
  for i = 1, n_elements(dimkeep)-1 do totkeep = totkeep * sz(dimkeep(i)+1)

  ;; this new array has two dimensions:
  ;;   * the first, all elements that will be collapsed
  ;;   * the second, all dimensions that will be preserved
  ;; (the ordering is so that all elements to be collapsed are
  ;;  adjacent in memory)
  newarr = reform(newarr, [totcol, totkeep], /overwrite)
end

;; Main function
function cmapply, op, array, dimapply, double=dbl, type=type, $
                  functargs=functargs, nocatch=nocatch

  if n_params() LT 2 then begin
      message, "USAGE: XX = CMAPPLY('OP',ARRAY,2)", /info
      message, '       where OP is +, *, AND, OR, MIN, MAX', /info
      return, -1L
  endif
  if NOT keyword_set(nocatch) then $
    on_error, 2 $
  else $
    on_error, 0

  version = double(!version.release)
;  version = 0
;  print, 'version = ',version

  ;; Parameter checking
  ;; 1) the dimensions of the array
  sz = size(array)
  if sz(0) EQ 0 then $
    message, 'ERROR: ARRAY must be an array!'

  ;; 2) The type of the array
  if sz(sz(0)+1) EQ 0 OR sz(sz(0)+1) EQ 7 OR sz(sz(0)+1) EQ 8 then $
    message, 'ERROR: Cannot apply to UNDEFINED, STRING, or STRUCTURE'
  if n_elements(type) EQ 0 then type = sz(sz(0)+1)

  ;; 3) The type of the operation
  szop = size(op)
  if szop(szop(0)+1) NE 7 then $
    message, 'ERROR: operation OP was not a string'

  ;; 4) The dimensions to apply (default is to apply to first dim)
  if n_params() EQ 2 then dimapply = 1
  dimapply = [ dimapply ]
  dimapply = dimapply(sort(dimapply))   ; Sort in ascending order
  napply = n_elements(dimapply)

  ;; 5) Use double precision if requested or if needed
  if n_elements(dbl) EQ 0 then begin
      dbl=0
      if type EQ 5 OR type EQ 9 then dbl=1
  endif

  newop = strupcase(op)
  newarr = array
  newarr = reform(newarr, sz(1:sz(0)), /overwrite)
  case 1 of

      ;; *** Addition
      (newop EQ '+'): begin
          for i = 0L, napply-1 do begin
              newarr = total(temporary(newarr), dimapply(i)-i, double=dbl)
          endfor
      end

      ;; *** Multiplication
      (newop EQ '*'): begin ;; Multiplication (by summation of logarithms)
          forward_function product

          cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep
          if nkeep EQ 0 then begin
              newarr = reform(newarr, n_elements(newarr), 1, /overwrite)
              if version GT 5.55 then return, product(newarr)
              return, (cmapply_product(newarr))(0)
          endif

          if version GT 5.55 then begin
              result = product(newarr,1)
          endif else begin
              result = cmapply_product(newarr)
          endelse
          result = reform(result, sz(dimkeep+1), /overwrite)
          return, result
      end

      ;; *** LOGICAL AND or OR
      ((newop EQ 'AND') OR (newop EQ 'OR')): begin
          newarr = temporary(newarr) NE 0
          totelt = 1L
          for i = 0L, napply-1 do begin
              newarr = total(temporary(newarr), dimapply(i)-i)
              totelt = totelt * sz(dimapply(i))
          endfor
          if newop EQ 'AND' then return, (round(newarr) EQ totelt)
          if newop EQ 'OR'  then return, (round(newarr) NE 0)
      end

      ;; Operations requiring a little more attention over how to
      ;; iterate
      ((newop EQ 'MAX') OR (newop EQ 'MIN') OR (newop EQ 'MEDIAN')): begin
          cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep
          if nkeep EQ 0 then begin
              if newop EQ 'MAX' then return, max(newarr)
              if newop EQ 'MIN' then return, min(newarr)
              if newop EQ 'MEDIAN' then return, median(newarr)
          endif

          ;; IDL 5.5 introduced the DIMENSION keyword to MAX() and MIN()
          if version GT 5.45 then begin
              extra = {dimension:1}
              if newop EQ 'MAX' then result = max(newarr, _EXTRA=extra)
              if newop EQ 'MIN' then result = min(newarr, _EXTRA=extra)
              if newop EQ 'MEDIAN' then result = median(newarr, _EXTRA=extra)
          endif else begin
              
              ;; Next task: create result array
              result = make_array(totkeep, type=type)
              
              ;; Now either iterate over the number of output elements, or
              ;; the number of collapsed elements, whichever is smaller.
              if (totcol LT totkeep) AND newop NE 'MEDIAN' then begin
                  ;; Iterate over the number of collapsed elements
                  result(0) = reform(newarr(0,*),totkeep,/overwrite)
                  case newop of 
                      'MAX': for i = 1L, totcol-1 do $
                        result(0) = result > newarr(i,*)
                      'MIN': for i = 1L, totcol-1 do $
                        result(0) = result < newarr(i,*)
                  endcase
              endif else begin
                  ;; Iterate over the number of output elements
                  case newop of 
                      'MAX': for i = 0L, totkeep-1 do result(i) = max(newarr(*,i))
                      'MIN': for i = 0L, totkeep-1 do result(i) = min(newarr(*,i))
                      'MEDIAN': for i = 0L, totkeep-1 do result(i) = median(newarr(*,i))
                  endcase
              endelse
          endelse

          result = reform(result, sz(dimkeep+1), /overwrite)
          return, result
      end

      ;; User function
      (strmid(newop,0,4) EQ 'USER'): begin
          functname = strmid(newop,5)
          if functname EQ '' then $
            message, 'ERROR: '+newop+' is not a valid operation'

          cmapply_redim, newarr, dimapply, dimkeep, nkeep, totcol, totkeep
          if nkeep EQ 0 then begin
              if n_elements(functargs) GT 0 then $
                return, call_function(functname, newarr, _EXTRA=functargs)
              return, call_function(functname, newarr)
          endif
          
          ;; Next task: create result array
          result = make_array(totkeep, type=type)
          
          ;; Iterate over the number of output elements
          if n_elements(functargs) GT 0 then begin
              for i = 0L, totkeep-1 do $
                result(i) = call_function(functname, newarr(*,i), _EXTRA=functargs)
          endif else begin
              for i = 0L, totkeep-1 do $
                result(i) = call_function(functname, newarr(*,i))
          endelse

          result = reform(result, sz(dimkeep+1), /overwrite)
          return, result
      end

              
  endcase

  newsz = size(newarr)
  if type EQ newsz(newsz(0)+1) then return, newarr

  ;; Cast the result into the desired type, if necessary
  castfns = ['UNDEF', 'BYTE', 'FIX', 'LONG', 'FLOAT', $
             'DOUBLE', 'COMPLEX', 'UNDEF', 'UNDEF', 'DCOMPLEX' ]
  if type GE 1 AND type LE 3 then $
    return, call_function(castfns(type), round(newarr)) $
  else $
    return, call_function(castfns(type), newarr)
end
  
;+
; NAME:
;   ARG_PRESENT
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   Determine whether output parameter has been passed (IDL4 compatibility)
;
; CALLING SEQUENCE:
;   PRESENT = ARG_PRESENT(ARG)
;
; DESCRIPTION: 
;
;   ARG_PRESENT tests whether an argument to a function or procedure
;   can be used as an output parameter.  The behavior of this function
;   is identical to that of the built-in ARG_PRESENT function in IDL
;   version 5 or greater, and is meant to give the same functionality
;   to programs in IDL 4.
;
;   An IDL procedure or function can use ARG_PRESENT to decide whether
;   the value of a positional or keyword parameter will be returned to
;   the calling procedure.  Generally, if the caller did not pass the
;   parameter then there is no need to compute the value to be
;   returned.
;
;   To be a valid output parameter, the caller must have passed a
;   named variable into which the result is stored.  If the caller
;   passed the parameter by value (e.g., an expression or a
;   subscripted array) the value cannot be returned and ARG_PRESENT
;   returns 0.
;
; INPUTS:
;
;   ARG - the parameter to be tested.  It can be either a positional
;         or a keyword parameter.  Passing a normal local variable
;         (i.e., not a passed parameter) will cause ARG_PRESENT to
;         return zero.
;
; RETURNS:
;
;   Returns a value of 1 if ARG is a valid output parameter, and a
;   value of 0 otherwise.
;
;
; EXAMPLE:
;
;   Consider the following procedure:
;      PRO TESTARG, ARG1
;        print, ARG_PRESENT(ARG1)
;      END
;
;   This procedure will print 1 when an ARG1 can be used as an output
;   parameter.  Here are some examples of the results of TESTARG.
;
;      IDL> testarg
;             0
;      IDL> testarg, x      
;             1
;      IDL> testarg, findgen(10)
;             0
;   
;   In the first case, no argument is passed, so ARG1 cannot be a
;   return variable.  In the second case, X is undefined, but it is
;   still a legal named variable capable of receiving an output
;   parameter.  In the third case, FINDGEN(10) is an expression which
;   cannot receive an output parameter.
;
; SEE ALSO:
;
;   ARG_PRESENT in IDL version 5
;
; MODIFICATION HISTORY:
;   Written, CM, 13 May 2000
;   Small documentation and bug fixes, CM, 04 Jul 2000
;
;-
; Copyright (C) 2000, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
; $Id: cmcongrid.pro,v 1.3 2007/03/29 13:52:20 craigm Exp $

;+
; NAME:
;	CMCONGRID
;
; PURPOSE:
;       Shrink or expand the size of an array by an arbitrary amount.
;       This IDL procedure simulates the action of the VAX/VMS
;       CONGRID/CONGRIDI function.
;
;	This function is similar to "REBIN" in that it can resize a
;       one, two, or three dimensional array.   "REBIN", however,
;       requires that the new array size must be an integer multiple
;       of the original size.   CONGRID will resize an array to any
;       arbitrary size (REBIN is somewhat faster, however).
;       REBIN averages multiple points when shrinking an array,
;       while CONGRID just resamples the array.
;
; CATEGORY:
;       Array Manipulation.
;
; CALLING SEQUENCE:
;	array = CONGRID(array, x, y, z)
;
; INPUTS:
;       array:  A 1, 2, or 3 dimensional array to resize.
;               Data Type : Any type except string or structure.
;
;       x:      The new X dimension of the resized array.
;               Data Type : Int or Long (greater than or equal to 2).
;
; OPTIONAL INPUTS:
;       y:      The new Y dimension of the resized array.   If the original
;               array has only 1 dimension then y is ignored.   If the
;               original array has 2 or 3 dimensions then y MUST be present.
;
;       z:      The new Z dimension of the resized array.   If the original
;               array has only 1 or 2 dimensions then z is ignored.   If the
;               original array has 3 dimensions then z MUST be present.
;
; KEYWORD PARAMETERS:
;       INTERP: If set, causes linear interpolation to be used.
;               Otherwise, the nearest-neighbor method is used.
;
;	CUBIC:	If set, uses "Cubic convolution" interpolation.  A more
;		accurate, but more time-consuming, form of interpolation.
;		CUBIC has no effect when used with 3 dimensional arrays.
;
;       MINUS_ONE:
;               If set, will prevent CONGRID from extrapolating one row or
;               column beyond the bounds of the input array.   For example,
;               If the input array has the dimensions (i, j) and the
;               output array has the dimensions (x, y), then by
;               default the array is resampled by a factor of (i/x)
;               in the X direction and (j/y) in the Y direction.
;               If MINUS_ONE is present (AND IS NON-ZERO) then the array
;               will be resampled by the factors (i-1)/(x-1) and
;               (j-1)/(y-1).
;
;       HALF_HALF:
;               If set, will tell CONGRID to extrapolate a *half* row
;               and column on either side, rather than the default of
;               one full row/column at the ends of the array.  If you
;               are interpolating images with few rows, then the
;               output will be more consistent with this technique.
;               This keyword is intended as a replacement for
;               MINUS_ONE, and both keywords probably should not be
;               used in the same call to CONGRID.
;
; OUTPUTS:
;	The returned array has the same number of dimensions as the original
;       array and is of the same data type.   The returned array will have
;       the dimensions (x), (x, y), or (x, y, z) depending on how many
;       dimensions the input array had.
;
; PROCEDURE:
;       IF the input array has three dimensions, or if INTERP is set,
;       then the IDL interpolate function is used to interpolate the
;       data values.
;       If the input array has two dimensions, and INTERP is NOT set,
;       then the IDL POLY_2D function is used for nearest neighbor sampling.
;       If the input array has one dimension, and INTERP is NOT set,
;       then nearest neighbor sampling is used.
;
; EXAMPLE:
;       ; vol is a 3-D array with the dimensions (80, 100, 57)
;       ; Resize vol to be a (90, 90, 80) array
;       vol = CONGRID(vol, 90, 90, 80)
;
; MODIFICATION HISTORY:
;       DMS, Sept. 1988.
;       DMS, Added the MINUS_ONE keyword, Sept. 1992.
; 	Daniel Carr. Re-wrote to handle one and three dimensional arrays
;                    using INTERPOLATE function.
; 	DMS, RSI, Nov, 1993.  Added CUBIC keyword.
;       Craig Markwardt, Dec, 1997.  Added halfhalf keyword to
;                        more evenly distribute "dead" pixel row
;       Use uniformly spaced grid points for half_half W. Landsman   Feb. 2000
;              (and slightly modified by C. Markwardt 14 Feb 2000)
;       Fix in case where INTERP=0 (nearest neighbor interp) and
;         expanding the image (thanks to Larry Bradley) 28 Mar 2007
;
;  $Id: cmcongrid.pro,v 1.3 2007/03/29 13:52:20 craigm Exp $
;-
;  Supply defaults = no interpolate, and no minus_one.
if n_elements(int) le 0 then int = 0 else int = keyword_set(int)
if n_elements(m1) le 0 then m1 = 0 else m1 = keyword_set(m1)

; Compute offsets pixel offsets for half_half
halfx = 0.0 & halfy = 0.0 & halfz = 0.0
if keyword_set(hh) then begin
    if s(0) GE 1 then halfx = -0.5 + (float(s(1))/x)
    if s(0) GE 2 then halfy = -0.5 + (float(s(2))/y)
    if s(0) GE 3 then halfz = -0.5 + (float(s(3))/z)
endif
cub = KEYWORD_SET(cubic)
if cub THEN int = 1	;Cubic implies interpolate
	

CASE s(0) OF
   1: BEGIN				; *** ONE DIMENSIONAL ARRAY
	srx = float(s(1) - m1)/(x-m1) * findgen(x) + halfx
      IF int THEN $
         RETURN, INTERPOLATE(arr, srx, CUBIC = cub) ELSE $
         RETURN, arr(ROUND(srx))
   ENDCASE
   2: BEGIN ; *** TWO DIMENSIONAL ARRAY
	IF int THEN BEGIN
	  srx = float(s(1) - m1) / (x-m1) * findgen(x) + halfx
	  sry = float(s(2) - m1) / (y-m1) * findgen(y) + halfy
          RETURN, INTERPOLATE(arr, srx, sry, /GRID, CUBIC=cub)
	ENDIF ELSE BEGIN
          ;; match IDL's CONGRID function
          expand = (x gt s[1])
          xm1 = (m1 or expand) ? x-1 : x
          RETURN, POLY_2D(arr, $
               [[0,0],[(s(1)-m1)/float(xm1),0]], $ ;Use poly_2d
               [[0,(s(2)-m1)/float(y-m1)],[0,0]],int,x,y)
        ENDELSE

   ENDCASE
   3: BEGIN ; *** THREE DIMENSIONAL ARRAY
	srx = float(s(1) - m1) / (x-m1) * findgen(x) + halfx
	sry = float(s(2) - m1) / (y-m1) * findgen(y) + halfy
	srz = float(s(3) - m1) / (z-m1) * findgen(z) + halfz
	RETURN, interpolate(arr, srx, sry, srz, /grid)
   ENDCASE
ENDCASE

RETURN, arr_r
END
;+
; NAME:
;   CMPRODUCT
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   CMPRODUCT() is the multiplicative equivalent of TOTAL().
;
; CALLING SEQUENCE:
;   Result = CMPRODUCT(ARRAY)
;
; DESCRIPTION: 
;
;   Calculates the product of all the elements of an array.  Vector
;   multiplication in groups of powers of two make this operation
;   faster than a simple FOR loop.  The number of actual
;   multiplications is still N_ELEMENTS(ARRAY).  Double precision
;   should be used for the highest accuracy when multiplying many
;   numbers.
;
; INPUTS:
;
;   ARRAY - Array of elements to multiply together.  For instance,
;           ARRAY could contain the dimensions of another array--then
;           CMPRODUCT(ARRAY) would be the total number of elements of
;           that other array.
;
; RETURNS:
;  The result of the function is the total product of all the elements
;  of ARRAY.
;
; EXAMPLE:
;
; SEE ALSO:
;
;   TOTAL, PRODUCT (from Astronomy User's Library)
;
; MODIFICATION HISTORY:
;   Written, CM, 28 Mar 2000
;     (based on outline of PRODUCT by William Thompson)
;
;  $Id: cmproduct.pro,v 1.2 2001/03/25 18:10:42 craigm Exp $
;
;-
; Copyright (C) 2000, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
;
;  Check the number of parameters.
;
  IF N_PARAMS() NE 1 THEN MESSAGE,'Syntax:  Result = PRODUCT(ARRAY)'
;
;  Check the type of ARRAY.
;
  SZ = SIZE(ARRAY)
  TYPE = SZ(SZ(0)+1)
  IF TYPE EQ 0 THEN MESSAGE,'ARRAY not defined'
  IF TYPE EQ 7 THEN MESSAGE,'Operation illegal with string arrays'
  IF TYPE EQ 8 THEN MESSAGE,'Operation illegal with structures'
;
;  Calculate the product.
;
  X = ARRAY
  N = N_ELEMENTS(X)
  WHILE N GT 1 DO BEGIN
      IF (N MOD 2) EQ 1 THEN X(0) = X(0) * X(N-1)
      N2 = FLOOR(N/2)
      X = X(0:N2-1) * X(N2:*)
      N = N2
  ENDWHILE
;
  RETURN,X(0)
END
;+
; NAME:
;   CMPS_FORM
;
; PURPOSE:
;   This function puts up a form the user can configure a PostScript
;   device driver. The function result (if the user selects either the
;   ACCEPT or CREATE FILE buttons) can be sent directly to the DEVICE
;   procedure by means of its _Extra keyword.  User's predefined
;   configurations may be retrieved from a common block.
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   $Id: cmps_form.pro,v 1.5 2004/10/03 09:40:08 craigm Exp $
;
;   Based almost entirely on, but a totally revamped version of, CMPS_FORM by 
;   FANNING SOFTWARE CONSULTING (David Fanning Ph.D.) http://www.dfanning.com
;
; MAJOR TOPICS:
;   Device Drivers, Hardcopy Output, PostScript Output
;
; PROCEDURE:
;   This is a pop-up form widget. It is a modal or blocking widget.
;   Keywords appropriate for the PostScript DEVICE command are returned.
;
;   Use your LEFT mouse button to move the "Plot Window" around the page.
;   Use your RIGHT mouse button to draw your own "Plot Window" on the page.
;
; HELP:
;   formInfo = CMPS_FORM(/Help)
;
; CALLING SEQUENCE:
;    formInfo = CMPS_FORM(xoffset, yoffset, Cancel=cancelButton)
;
; OPTIONAL INPUTS:
;
;    XOFFSET -- Optional xoffset of the top-level base of cmps_form. Default is
;    to try to center the form on the display.
;
;    YOFFSET -- Optional yoffset of the top-level base of cmps_form. Default is
;    to try to center the form on the display.
;
; INPUT KEYWORD PARAMETERS:
;
;    BITS_PER_PIXEL -- The initial configuration of the bits per pixel button.
;
;    BLOCKING -- Set this keyword to make this a blocking widget under IDL 5.0.
;    (All widget programs block under IDL 4.0.)
;
;    COLOR -- The initial configuration of the color switch.
;
;    DEFAULTS -- A stucture variable of the same type and structure as the
;    RETURN VALUE of cmps_form. It will set initial conditions. This makes
;    it possible to start cmps_form up again with the same values it had the
;    last time it was called. For example:
;
;       mysetup = cmps_form()
;       newsetup = cmps_form(Defaults=mysetup)
;
;    ENCAPSULATED -- The initial configuration of the encapsulated switch.
;
;    FILENAME -- The initial filename to be used on the form.
;
;    HELP -- Prints a helpful message in the output log.
;
;    INCHES -- The initial configuration of the inches/cm switch.
;
;    INITIALIZE -- If this keyword is set, the program immediately returns the
;    "localdefaults" structure. This gives you the means to configue the
;    PostScript device without interrupting the user.
;
;    SELECT -- used only when INITIALIZE is set.  Set SELECT to a
;              string which identifies the predefined configuration to
;              be returned by cmps_form when INITIALIZE is set.  This is
;              a convenient way to select a predefined config
;              non-interactively.
;
;    LANDSCAPE -- The initial configuration of the landscape/portrait switch.
;
;    LOCALDEFAULTS -- A structure like the DEFAULTS structure. If specified,
;    then it is added as a predefined configuration entry called "Local".
;    See below for a further discussion of predefined configurations.
;
;    PREDEFINED -- An alternate way to specify predefined
;                  configurations.  Pass an array of structures to
;                  populate the "predefined" dropbox in the
;                  dialog. This array, if specified, overrides the the
;                  common block technique.
;
;    XOFFSET -- The initial XOffSet of the PostScript window.
;
;    YOFFSET -- The initial YOffSet of the PostScript window.
;
;    XSIZE -- The initial XSize of the PostScript window.
;
;    YSIZE -- The initial YSize of the PostScript window.
;
;    ASPECT -- The aspect ratio of the window (Y/X).  This keyword can
;              substitute for one of XSIZE or YSIZE.
; 
;    PRESERVE_ASPECT -- Set this keyword if you want to hold the
;                       aspect ratio constant.
;
;    PAPERSIZE -- If set, allows user to specify the size of the paper
;                 media to be printed on, as a scalar string.  NOTE:
;                 this specification cannot be passed to DEVICE, but
;                 can be selected for completeness's sake.  Default is
;                 'Letter'.
;
;    MARGINSIZE -- Size of the margins on all sides.  Default is 0.25 inches.
;                  When MARGINSIZE is non-zero, a graphic cannot directly
;                  abut the edge of the page.  This is normally a good thing,
;                  since there is often a  non-printable region which borders
;                  the page.
;
;   DEFAULTPAPER -- Default paper size to use, when it is unspecified
;                   in a predefined, "local", or "default"
;                   configuration.  This value also overrides any
;                   configuration from common blocks.  European users
;                   will probably set this to 'A4'.
;
;   PARENT -- if this widget is invoked by another widget program,
;             then this keyword parameter must be set to the top level
;             widget which is to serve as the group leader.  Failure
;             to do so will result in unexpected behavior.  IDL 4
;             programs do not need to pass this parameter.  Default:
;             NONE.
;
; OUTPUT KEYWORD PARAMETERS
;
;    CANCEL -- This is an OUTPUT keyword. It is used to check if the user
;    selected the "Cancel" button on the form. Check this variable rather
;    than the return value of the function, since the return value is designed
;    to be sent directly to the DEVICE procedure. The varible is set to 1 if
;    the user selected the "Cancel" button. Otherwise, it is set to 0.
;
;    CREATE -- This output keyword can be used to determine if the user
;    selected the 'Create File' button rather than the 'Accept' button.
;    The value is 1 if selected, and 0 otherwise.
;
;    PAPERSIZE -- If set to a named variable, any newly selected paper
;    size is returned in that variable.
;
;    XPAGESIZE -- Size of paper in "X" dimension, in units given by
;                 the returned config structure.
;
;    YPAGESIZE -- Size of paper in "Y" dimension, in units given by
;                 the returned config structure.
;
;    PAGEBOX -- specifies the page rectangle relative to the plot
;               window, in normalized units.  A 4-vector of the form
;               [XLL, YLL, XUR, YUR] is returned, giving the positions
;               of the lower left (LL) and upper right (UR) corners of
;               the page with respect to the plot window.  Thus, the
;               following command:
;
;                    PLOT, x, y, position=PAGEBOX
;
;               will construct a graphic whose plot region exactly
;               fills the page (with no margin around the edges).
;
;               Naturally, the page is usually larger than the
;               graphics window, so the normalized coordinates will
;               usually fall outside the range [0,1].
;
;               However, the bounding box constructed by the
;               Postscript driver includes only the graphics window.
;               Anything drawn outside of it may be clipped or
;               discarded.
;
; RETURN VALUE:
;
;     formInfo = { cmps_form_INFO, $
;                  xsize:0.0, $        ; The x size of the plot
;                  xoff:0.0, $         ; The x offset of the plot
;                  ysize:0.0, $        ; The y size of the plot
;                  yoff:0.0 $          ; The y offset of the plot
;                  filename:'', $      ; The name of the output file
;                  inches:0 $          ; Inches or centimeters?
;                  color:0, $          ; Color on or off?
;                  bits_per_pixel:0, $ ; How many bits per image pixel?
;                  encapsulated:0,$    ; Encapsulated or regular PostScript?
;                  isolatin1:0,$       ; Encoded with ISOLATIN1?
;                  landscape:0 }       ; Landscape or portrait mode?
;
; USAGE:
;
;  The calling procedure for this function in a widget program will
;  look something like this:
;
;     info.ps_config = cmps_form(/Initialize)
;
;     formInfo = cmps_form(Cancel=canceled, Create=create, $
;                          Defaults=info.ps_config)
;
;     IF NOT canceled THEN BEGIN
;        IF create THEN BEGIN
;           thisDevice = !D.Name
;           Set_Plot, "PS"
;           Device, _Extra=formInfo
;
;           Enter Your Graphics Commands Here!
;
;           Device, /Close
;           Set_Plot, thisDevice
;           info.ps_config = formInfo
;        ENDIF ELSE
;           info.ps_config = formInfo
;     ENDIF
;
; MAJOR FUNCTIONS and PROCEDURES:
;
;   None. Designed to work originally in conjunction with XWindow, a
;   resizable graphics window.  [ NOTE: this modified version of
;   cmps_form, by Craig Markwardt, is incompatible with the original
;   version of XWINDOW. ]
;
; MODIFICATION HISTORY:
;
;   Based on cmps_form of : David Fanning, RSI, March 1995.
;   Major rewrite by: Craig Markwardt, October 1997.
;     - Drawing and updating of form and sample box are now modular
;     - Option of storing more than one predefined postscript configuration
;     - Selection of paper size by name
;     - Access to predfined configurations through (optional) common
;       block
;   Several additions, CM, April 1998  VERSION CM2.0
;     - better integration of paper sizes throughout program.
;       Predefined configurations now also know about paper.
;     - allow passing predefined configurations instead of using
;       common block
;     - addition of ISOLATIN selection, and streamlining of dialog
;       appearance
;   Fixed bug in INITIALIZE w.r.t. paper sizes, CM, Nov 1998
;   Added SELECT keyword, CM, 09 Dec 1998
;   Added Parent keyword to allow modal widgets in IDL 5, 19 Jan 1999
;   Added "Choose" button for filename selection, 19 Sep 1999
;   Added ability to program different button names, 19 Sep 1999
;   Added ASPECT and PRESERVE_ASPECT, based on work by Aaron Barth, 18
;     Oct 1999
;   Removed NOCOMMON documentation and logic, 19 Oct 1999, CM
;   Added aspect to cmps_form_numevents (per Aaron Barth), 18 Oct 1999
;   Corrected small bug under Initialize keyword (inches), 18 Oct 1999
;   Made call to *_pscoord more consistent, 18 Oct 1999
;   Added XPAGESIZE, YPAGESIZE and PAGEBOX keywords, 19 Oct 1999
;   Small cosmetic cleanup, CM, 01 Feb 2000
;   Fix for IDL 5.5's handling of structures with arrays, CM, 11 Dec 2001
;   Replaced obsolete PICKFILE call with DIALOG_PICKFILE, Jeff Guerber,
;     24 Sep 2004
;   Transfer DEFAULTS and LOCALDEFAULTS values via STRUCT_ASSIGN,/NOZERO
;     instead of EXECUTE, Jeff Guerber, 24 Sep 2004.
;   Set CANCELBUTTON and CREATEBUTTON immediately on entry, so they're
;     defined even if user kills the window, Jeff Guerber, 24 Sep 2004.
;
; COMMON BLOCKS:
; 
;   The user may store frequently used or helpful configurations in a 
;   common block, and cmps_form() will attempt to access them.  This
;   provides a way for the user to have persistent, named,
;   configurations.
;
;   NOTE: this format has changed since the last version.  You may
;   have to quit your IDL session for the changes to take effect
;   properly.  If you have place a predefined configuration in your
;   startup file, you should review the new format.
;   
;     COMMON CMPS_FORM_CONFIGS, cmps_form_DEFAULT_PAPERSIZE, $
;                               cmps_form_STDCONFIGS
;
;        cmps_form_DEFAULT_PAPERSIZE - a string designating the default
;                                    paper size, when none is given.
;                                    The predefined configurations
;                                    offerred by this program will
;                                    respect the default value.  (See
;                                    also the DEFAULTPAPER keyword.)
;
;        cmps_form_STDCONFIGS - An array of cmps_form_CONFIG structures,
;                             each containing information about one
;                             predefined configuration, such as its
;                             name and size of paper.  Each "config"
;                             element is a cmps_form_INFO structure,
;                             which contains the actual postscript
;                             configuration.
;
;   See the IDL source code cmps_form_LOAD_CONFIGS for an example of how
;   to make a list of configurations.  One possibility would be to
;   declare and populate the common block from within the user's
;   start-up script, allowing the same configurations to appear in
;   every session.
;
;   cmps_form() takes its initial list of configurations from this
;   common block if it exists.  A default list is provided ala the
;   procedure cmps_form_LOAD_CONFIGS.  Any modifications that take place
;   during the cmps_form() widget session are not transferred back to
;   the common block upon return.  It might be useful to be able to do
;   this, through some form of 'save' procedure.
;
;   Also, if the PREDEFINED keyword is used, then the common block is
;   not consulted.
;
;  $Id: cmps_form.pro,v 1.5 2004/10/03 09:40:08 craigm Exp $
;
;-
; Copyright (C) 1996-1997, David Fanning
; Copyright (C) 1997-2001, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
; ***************************************************************
; Utility routines

forward_function filepath

; Convert from inches and centimeters to WIDGET_DRAW pixels 
pro cmps_form_Draw_Coords, drawpixperunit, xoff, yoff, xsize, ysize

  if n_elements(xoff) GT 0 then xoff   = xoff  * drawpixperunit + 1
  if n_elements(yoff) GT 0 then yoff   = yoff  * drawpixperunit + 1
  if n_elements(xsize) GT 0 then xsize = xsize * drawpixperunit
  if n_elements(ysize) GT 0 then ysize = ysize * drawpixperunit

  return
end

; Perform the opposite conversion of cmps_form_DRAW_COORDS
pro cmps_form_Real_Coords, drawpixperunit, xoff, yoff, xsize, ysize

  if n_elements(xoff) GT 0 then xoff   = (xoff-1) / drawpixperunit
  if n_elements(yoff) GT 0 then yoff   = (yoff-1) / drawpixperunit
  if n_elements(xsize) GT 0 then xsize = xsize / drawpixperunit
  if n_elements(ysize) GT 0 then ysize = ysize / drawpixperunit

  return
end

Pro cmps_form_Select_File, event

   ; Allows the user to select a filename for writing.

Widget_Control, event.top, Get_UValue=info, /No_Copy

   ; Start with the name in the filename widget.

Widget_Control, info.idfilename, Get_Value=initialFilename
initialFilename = initialFilename(0)
filename = Dialog_Pickfile(/Write, File=initialFilename)
IF filename NE '' THEN $
   Widget_Control, info.idfilename, Set_Value=filename
Widget_Control, event.top, Set_UValue=info, /No_Copy
END ;*******************************************************************

; Calculate a list of vertices to be plotted as a box in the 
; draw widget.
Function cmps_form_PlotBox_Coords, xsize, ysize, xoff, yoff, drawpixperunit

   ; This function converts sizes and offsets to appropriate
   ; Device coordinates for drawing the PLOT BOX on the PostScript
   ; page. The return value is a [2,5] array.

returnValue = IntArr(2,5)
xs = xsize
ys = ysize
xof = xoff
yof = yoff
cmps_form_draw_coords, drawpixperunit, xof, yof, xs, ys

; Add one because we do for the page outline
xcoords = Round([xof, xof+xs, xof+xs, xof, xof]) + 1
ycoords = Round([yof, yof, yof+ys, yof+ys, yof]) + 1

returnValue(0,*) = xcoords
returnValue(1,*) = ycoords

RETURN, returnValue
END ;*******************************************************************

; Convert between the IDL-form of PS coordinates (including the 
; strange definition of YOFFSET and XOFFSET) to a more
; "human-readable" form where the Xoffset and YOFFSET always refer to
; the lower-left hand corner of the output
pro cmps_form_conv_pscoord, info, xpagesize, ypagesize, $
      toidl=toidl, tohuman=tohuman

  if info.landscape EQ 1 then begin
      ixoff=info.xoff 
      iyoff=info.yoff
      if keyword_set(tohuman) then begin
          info.yoff = ixoff
          info.xoff = xpagesize - iyoff
      endif else if keyword_set(toidl) then begin
          info.xoff = iyoff
          info.yoff = xpagesize - ixoff
      endif
  endif
  return
end

; Return names of paper sizes
function cmps_form_papernames

  return, ['Letter','Legal','Tabloid','Ledger','Executive','Monarch', $
            'Statement','Folio','Quarto','C5','B4','B5','Dl','A0','A1', $
            'A2','A3','A4','A5','A6']
end

; Select a paper size based on number or string.  Returns x and 
; y page sizes, accouting for the units of measurement and the
; orientation of the page.
pro cmps_form_select_papersize, papertype, xpagesize, ypagesize, $
      inches=inches, landscape=landscape, index=index

;           Letter Legal Tabloid Ledger Executive Monarch Statement Folio
  xpaper = [612.,  612,  792,    792,   540,      279,    396,      612, $
$;          Quarto C5  B4  B5  Dl  A0   A1   A2   A3  A4  A5  A6
	    610,   459,729,516,312,2380,1684,1190,842,595,420,297]

;           Letter Legal Tabloid Ledger Executive Monarch Statement Folio
  ypaper = [792.,  1008, 1224,   1224,  720,      540,    612,      936, $
$;          Quarto C5  B4   B5  Dl  A0   A1   A2   A3   A4  A5  A6
	    780,   649,1032,729,624,3368,2380,1684,1190,842,595,421]
  names  = cmps_form_papernames()

  sz = size(papertype)
  tp = sz(sz(0) + 1)
  if tp GT 0 AND tp LT 6 then begin
      index = fix(papertype)
  endif else if tp EQ 7 then begin
      index = where(strupcase(papertype) EQ strupcase(names), ict)
      if ict EQ 0 then index = 0
  endif else $
    index = 0

  index = index(0)

  xpagesize = xpaper(index) / 72.  ; Convert to inches
  ypagesize = ypaper(index) / 72.
  xpagesize = xpagesize(0)
  ypagesize = ypagesize(0)
  if NOT keyword_set(inches) then begin
      xpagesize = xpagesize * 2.54
      ypagesize = ypagesize * 2.54
  endif
  if keyword_set(landscape) then begin
      temp = xpagesize
      xpagesize = ypagesize
      ypagesize = temp
  endif
  
  return
end

; cmps_form_LOAD_CONFIGS
;
; Loads a set of default configurations into the output variables,
;
;   CONFIGNAMES - array of names for configurations.
;
;   CONFIGS - array of cmps_form_INFO structures, each with a separate
;             configuration in it, and corresponding to the
;             configuration name.
;
; Intended as an intelligent default when no other is specified.
;
pro cmps_form_load_configs, defaultpaper, configs

  ; This is the default paper size, when none is given
  defaultpaper = 'Letter'

  ; Here is how the cmps_form_INFO structure is defined.  Refer to it
  ; when creating new structures.
  template = { cmps_form_INFO, $
               xsize:0.0, $     ; The x size of the plot
               xoff:0.0, $      ; The x offset of the plot
               ysize:0.0, $     ; The y size of the plot
               yoff:0.0, $       ; The y offset of the plot
               filename:'', $   ; The name of the output file
               inches:0, $       ; Inches or centimeters?
               color:0, $       ; Color on or off?
               bits_per_pixel:0, $ ; How many bits per image pixel?
               encapsulated:0,$ ; Encapsulated or regular PostScript?
	       isolatin1:0,$    ; Encoding is not ISOLATIN1
               landscape:0 }    ; Landscape or portrait mode?

  pctemplate = { cmps_form_CONFIG,  $
                 config:{cmps_form_INFO},  $
                 configname: '', $    ; Name of configuration
                 papersize: '' }      ; Size of paper for configuration
                 
  
  ; Set of default configurations (no ISOLATIN1 encoding)
  ;   1.  7x5    inch color plot region in portrait
  ;   2.  7.5x10 inch centered color plot region, covering almost whole
  ;                   portrait page (0.5 inch margins)
  ;   3.  10x7.5 inch centered color plot region, covering almost whole
  ;                   landscape page (0.5 inch margins)
  ;   4.  7x5    inch gray plot region in portrait (IDL default config)
  configs = [{cmps_form_CONFIG, config:$
              {cmps_form_INFO, 7.0, 0.75, 5.0, 5.0, 'idl.ps', 1, 1, 8, 0, 0, 0},$
              configname:'Half Portrait (color)', papersize:defaultpaper}, $
             {cmps_form_CONFIG, config:$
              {cmps_form_INFO, 7.5, 0.50, 10., 0.5, 'idl.ps', 1, 1, 8, 0, 0, 0},$
              configname:'Full Portrait (color)', papersize:defaultpaper}, $
             {cmps_form_CONFIG, config:$
              {cmps_form_INFO, 10., 0.50, 7.5, 10.5,'idl.ps', 1, 1, 8, 0, 0, 1},$
              configname:'Full Landscape (color)', papersize:defaultpaper}, $
             {cmps_form_CONFIG, config:$
              {cmps_form_INFO, 18., 1.5, 26.7, 1.5, 'idl.ps', 0, 1, 8, 0, 0, 0},$
              configname:'A4 Portrait (color)', papersize:'A4'}, $
             {cmps_form_CONFIG, config:$
              {cmps_form_INFO, 26.7, 1.5, 18.,28.2039,'idl.ps',0,1, 8, 0, 0, 1},$
              configname:'A4 Landscape (color)', papersize:'A4'}, $
             {cmps_form_CONFIG, config:$
              {cmps_form_INFO, 17.78,1.91,12.70,12.70,'idl.ps',0,1, 4, 0, 0, 0},$
              configname:'IDL Standard', papersize:defaultpaper} ]


  return
end

; 
; cmps_form_Update_Info
;
; This procedure modifies an "info" structure, according to new
; specifications about the PS configuration.  This is the central
; clearing house for self-consistent modification of the info structure.
;
; INPUTS
;   info    - info structure to be modified
;   keywords- IDL keywords are contain information is folded
;             into the "info" structure.
;             Valid keywords are:
;                XSIZE, YSIZE, 
;                XOFF, YOFF    - size and offset of plotting region in
;                                "human" coordinates.  This is the
;                                natural size as measured from the
;                                lower-left corner of the page in its
;                                proper orientation (not the IDL
;                                definition!).  These are the same
;                                values that are printed in the form's
;                                Size and Offset fields.
;                INCHES        - whether dimensions are in inches or
;                                centimeters (1=in, 0=cm)
;                COLOR         - whether output is color (1=y, 0=n)
;                BITS_PER_PIXEL- number of bits per pixel (2,4,8)
;                ENCAPSULATED  - whether output is EPS (1=EPS, 0=PS)
;                LANDSCAPE     - whether output is portrait or
;                                landscape (1=land, 0=port)
;                FILENAME      - output file name (with respect to
;                                current directory)
;
Pro cmps_form_Update_Info, info, set=set, _EXTRA=newdata

  if n_elements(newdata) GT 0 then $
    names = Tag_Names(newdata)
  set   = keyword_set(set)
  centerfactor = 1.0

  FOR j=0, N_Elements(names)-1 DO BEGIN
      
      case strupcase(names(j)) of
          'XSIZE':       info.devconfig.xsize    = float(newdata.xsize)
          'YSIZE':       info.devconfig.ysize    = float(newdata.ysize)
          'XOFF':        info.devconfig.xoff     = float(newdata.xoff)
          'YOFF':        info.devconfig.yoff     = float(newdata.yoff)
          'INCHES':      BEGIN
              inches   = fix(newdata.inches)
              if inches NE 0 then inches = 1
              if set NE 1 then begin
                  convfactor = 1.0
                  if info.devconfig.inches EQ 0 AND inches EQ 1 then $
                    convfactor = 1.0/2.54 $ ; centimeters to inches
                  else if info.devconfig.inches EQ 1 AND inches EQ 0 then $
                    convfactor = 2.54 ; inches to centimeters

                  info.devconfig.xsize = info.devconfig.xsize * convfactor
                  info.devconfig.ysize = info.devconfig.ysize * convfactor
                  info.devconfig.xoff  = info.devconfig.xoff  * convfactor
                  info.devconfig.yoff  = info.devconfig.yoff  * convfactor
                  info.xpagesize       = info.xpagesize       * convfactor
                  info.ypagesize       = info.ypagesize       * convfactor
                  info.marginsize      = info.marginsize      * convfactor
                  info.drawpixperunit  = info.drawpixperunit  / convfactor
                    
              endif
                  
              info.devconfig.inches = inches
          end
          
          'LANDSCAPE':   begin
              landscape= fix(newdata.landscape)
              if landscape NE 0 then landscape = 1
              
              if landscape NE info.devconfig.landscape AND $
                set NE 1 then begin
                  temp = info.xpagesize
                  info.xpagesize = info.ypagesize
                  info.ypagesize = temp
                  
                  ; Since the margins are bound to be way out of wack,
                  ; we could recenter here.
                  
                  xsize = info.devconfig.xsize
                  ysize = info.devconfig.ysize

                  centerfactor = 2.0
                  
                  ; We will have to redraw the reserve pixmap
                  info.pixredraw = 1
                  
              endif 
              
              info.devconfig.landscape = landscape
          end
          
          'COLOR':       begin
              info.devconfig.color    = fix(newdata.color)
              if info.devconfig.color        NE 0 then info.devconfig.color = 1
          end
          'ENCAPSULATED':  begin
              info.devconfig.encapsulated = fix(newdata.encapsulated)
              if info.devconfig.encapsulated NE 0 then $
                info.devconfig.encapsulated = 1
          end
          'ISOLATIN1': begin
              info.devconfig.isolatin1 = fix(newdata.isolatin1)
              if info.devconfig.isolatin1 NE 0 then $
                info.devconfig.isolatin1 = 1
          end
          'BITS_PER_PIXEL': begin
              bpp = fix(newdata.bits_per_pixel)
              if bpp LT 1              then bpp = 2
              if bpp GT 2 AND bpp LT 4 then bpp = 4
              if bpp GT 4 AND bpp LT 8 then bpp = 8
              if bpp GT 8              then bpp = 8
              info.devconfig.bits_per_pixel = bpp
          end
          'FILENAME': begin
              if string(newdata.filename) NE info.devconfig.filename then $
                info.filechanged = 1
              info.devconfig.filename = string(newdata.filename)
          end
      endcase
  endfor
  
  ; Now check the sizes and offsets, to be sure they are sane for the 
  ; particular landscape/portrait and inch/cm settings that have been
  ; chosen.
  pgwid = info.xpagesize
  pglen = info.ypagesize
  pgmar = info.marginsize

  if set NE 1 then begin
      info.devconfig.xsize = (pgmar) > info.devconfig.xsize < (pgwid-2.*pgmar)
      info.devconfig.ysize = (pgmar) > info.devconfig.ysize < (pglen-2.*pgmar)
      info.devconfig.xoff  = (pgmar) > info.devconfig.xoff  < (pgwid-info.devconfig.xsize - pgmar)
      info.devconfig.yoff  = (pgmar) > info.devconfig.yoff  < (pglen-info.devconfig.ysize - pgmar)
      if info.devconfig.xsize + info.devconfig.xoff GT (pgwid-pgmar) then $
        info.devconfig.xoff = (pgwid - info.devconfig.xsize) / centerfactor
      if info.devconfig.ysize + info.devconfig.yoff GT (pglen-pgmar) then $
        info.devconfig.yoff = (pglen - info.devconfig.ysize) / centerfactor
  endif

  ; Preserve aspect ratio if necessary
  if (info.preserve_aspect EQ 1) then begin
      sizeratio = info.aspect / (info.ypagesize / info.xpagesize)
      if (sizeratio GE 1) then $
        info.devconfig.xsize = info.devconfig.ysize / info.aspect $
      else $
        info.devconfig.ysize = info.devconfig.xsize * info.aspect
  endif

  return
end

; 
; PRO cmps_form_DRAW_BOX
;
; Draw the "sample" box in the draw widget.  If necessary, also
; redraws the backing reserve pixmap.
;
pro cmps_form_draw_box, xsize, ysize, xoff, yoff, info

  ; First order of business is to make a new reserve pixmap, if
  ; necessary.

  if info.pixredraw EQ 1 then begin
      
      ; Operate on the pixmap first
      wset, info.idpixwid
      erase
      ; Make background ...
      tv, replicate(info.bkgcolor, info.xpixwinsize, info.ypixwinsize)
      ; ... and page outline
      coords = cmps_form_plotbox_coords(info.xpagesize, info.ypagesize, $
                                      0.,0., info.drawpixperunit)
      plots, coords(0,*), coords(1,*), /device, color=info.pagecolor

      info.pixredraw = 0
  endif

  ; Now, we transfer the reserve pixmap to the screen

  wset, info.idwid
  device, copy=[0, 0, info.xpixwinsize, info.ypixwinsize, 0, 0, $
                info.idpixwid]

  ; Finally we overlay the plot region
  coords = cmps_form_plotbox_coords(xsize, ysize, xoff, yoff,info.drawpixperunit)
  plots, coords(0,*), coords(1,*), color=info.boxcolor, /device

  return
end

; 
; cmps_form_DRAW_FORM
;
; Update the widget elements of the cmps_form form, using the INFO structure.
;
; If the NOBOX keyword is set, then the draw widget is not updated.
;
pro cmps_form_draw_form, info, nobox=nobox

  ; Draw the DRAW widget if needed
  if NOT keyword_set(nobox) then $
    cmps_form_draw_box, info.devconfig.xsize, info.devconfig.ysize, $
    info.devconfig.xoff, info.devconfig.yoff, info

  ; Update the numeric text fields
  xsizestr = strtrim(string(info.devconfig.xsize, format='(F6.2)'), 2)
  ysizestr = strtrim(string(info.devconfig.ysize, format='(F6.2)'), 2)
  xoffstr  = strtrim(string(info.devconfig.xoff, format='(F6.2)'), 2)
  yoffstr  = strtrim(string(info.devconfig.yoff, format='(F6.2)'), 2)

  widget_control, info.idxsize, set_value=xsizestr
  widget_control, info.idysize, set_value=ysizestr
  widget_control, info.idxoff, set_value=xoffstr
  widget_control, info.idyoff, set_value=yoffstr

  widget_control, info.idaspect, set_button=(info.preserve_aspect EQ 1)

  ; Set EPS (encapsulated ps) buttons
  Widget_Control, info.idencap, Set_Button=(info.devconfig.encapsulated EQ 1)

   ; Set color buttons.
  Widget_Control, info.idcolor, Set_Button=(info.devconfig.color EQ 1)

  ; Set inch/cm buttons.
  Widget_Control, info.idinch, Set_Button=(info.devconfig.inches EQ 1)
  Widget_Control, info.idcm,   Set_Button=(info.devconfig.inches EQ 0)

  ; Set bits_per_pixel buttons.
  Widget_Control, info.idbit2, Set_Button=(info.devconfig.bits_per_pixel EQ 2)
  Widget_Control, info.idbit4, Set_Button=(info.devconfig.bits_per_pixel EQ 4)
  Widget_Control, info.idbit8, Set_Button=(info.devconfig.bits_per_pixel EQ 8)
  Widget_Control, info.idbitbase, Sensitive=(info.devconfig.color EQ 1)

  ; Set encoding button
  widget_control, info.idisolatin1, Set_Button=(info.devconfig.isolatin1 EQ 1)

  ; Set default filename.
  Widget_Control, info.idfilename, Get_Value=wfilename
  if string(wfilename(0)) NE info.devconfig.filename then begin
      Widget_Control, info.idfilename, Set_Value=info.devconfig.filename
      ; Put caret at end of pathname text so that filename itself is visible
      Widget_Control, info.idfilename, $
        Set_Text_Select=[ strlen(info.devconfig.filename), 0 ]
  endif

  ; Set protrait/landscape button.
  Widget_Control, info.idland, Set_Button=(info.devconfig.landscape EQ 1)
  Widget_Control, info.idport, Set_Button=(info.devconfig.landscape EQ 0)

  ; Set Paper
  pn = cmps_form_papernames()
  xp = strtrim(string(info.xpagesize, format='(F10.2)'),2)
  yp = strtrim(string(info.ypagesize, format='(F10.2)'),2)
  un = 'in'
  if NOT info.devconfig.inches then un = 'cm'

  paperlab = string(pn(info.paperindex), xp, un, yp, un, $
                      format='("   Paper: ",A0," (",A0,A0," x ",A0,A0,")   ")')
  Widget_Control, info.idpaperlabel, set_value=paperlab
  
  return
end

Pro cmps_form_Null_Events, event
END ;*******************************************************************

Function cmps_form_What_Button_Type, event

   ; Checks event.type to find out what kind of button
   ; was clicked in a draw widget. This is NOT an event handler.

type = ['DOWN', 'UP', 'MOTION', 'SCROLL']
Return, type(event.type)
END ;*******************************************************************

Function cmps_form_What_Button_Pressed, event

   ; Checks event.press to find out what kind of button
   ; was pressed in a draw widget.  This is NOT an event handler.

button = ['NONE', 'LEFT', 'MIDDLE', 'NONE', 'RIGHT']
Return, button(event.press)
END ;*******************************************************************

Function cmps_form_What_Button_Released, event

   ; Checks event.release to find out what kind of button
   ; was released in a draw widget.  This is NOT an event handler.

button = ['NONE', 'LEFT', 'MIDDLE', 'NONE', 'RIGHT']
Return, button(event.release)
END ;*******************************************************************



;
; cmps_form_NUMEVENTS
;
; Events sent to the numeric text field widgets are sent here.  We
; harvest the data values from the text field and update the screen.
;
Pro cmps_form_NumEvents, event

   ; If an event comes here, read the offsets and sizes from the
   ; form and draw the appropriately sized box in the draw widget.

Widget_Control, event.top, Get_UValue= info, /No_Copy

   ; Get current values for offset and sizes

Widget_Control, info.idxsize, Get_Value=xsize
Widget_Control, info.idysize, Get_Value=ysize
Widget_Control, info.idxoff, Get_Value=xoff
Widget_Control, info.idyoff, Get_Value=yoff

xsize = xsize(0)
ysize = ysize(0)
xoff = xoff(0)
yoff = yoff(0)

if info.preserve_aspect EQ 1 then begin
    if event.id EQ info.idysize then xsize = ysize / info.aspect $
    else                             ysize = xsize * info.aspect
endif

; Fold this information into the "info" structure
cmps_form_update_info, info, xsize=xsize, ysize=ysize, xoff=xoff, yoff=yoff

; Update form and redraw sample box
cmps_form_draw_form, info

; Put the info structure back into the top-level base

Widget_Control, event.top, Set_UValue=info, /No_Copy

END ;*******************************************************************


Pro cmps_form_Move_Box, event

   ; This is the event handler that allows the user to "move"
   ; the plot box around in the page window. It will set the
   ; event handler back to "cmps_form_Box_Events" when it senses an
   ; "UP" draw button event and it will also turn cmps_form_Draw_Motion_Events
   ; OFF.

   ; Get the info structure out of the top-level base.

  Widget_Control, event.top, Get_UValue=info, /No_Copy

  whatButtonType = cmps_form_What_Button_Type(event)

  dpu = info.drawpixperunit

  ixmin = 0.
  iymin = 0.
  ixsize = info.devconfig.xsize
  iysize = info.devconfig.ysize
  cmps_form_draw_coords, dpu, ixmin, iymin, ixsize, iysize
  ; Now ixmin,iymin have the minimum values of x and y, in pixels
  ; ixsize and iysize are the size of the box, in pixels

  ixmax = info.xpagesize
  iymax = info.ypagesize
  cmps_form_draw_coords, dpu, ixmax, iymax
  ; ixmax and iymax are the max values of x and y, in pixels

  ; info.ideltx/y contains the offset of the lower left corner of the box,
  ; with respect to the mouse's position
  ixoff = event.x + info.ideltx
  iyoff = event.y + info.idelty

  ; Keep box inside the page
  if ixoff LT ixmin then ixoff = ixmin
  if iyoff LT iymin then iyoff = iymin
  if (ixoff+ixsize) GT ixmax then ixoff = ixmax - ixsize
  if (iyoff+iysize) GT iymax then iyoff = iymax - iysize

  IF whatButtonType EQ 'UP' THEN Begin

      ; When the button is "up" the moving event is over.  We reset the
      ; event function and update the information about the box's position
      
      Widget_Control, info.iddraw, Draw_Motion_Events=0, $ ; Motion events off
        Event_Pro='cmps_form_Box_Events' ; Change to normal processing
      
      cmps_form_real_coords, dpu, ixoff, iyoff, ixsize, iysize
      
      ; Update the info structure
      cmps_form_update_info, info, xoff=ixoff, yoff=iyoff
      ; Draw it
      cmps_form_draw_form, info
      
      ; Put the info structure back in the top-level base and RETURN
      Widget_Control, event.top, Set_UValue=info, /No_Copy
      Return

  ENDIF

   ; You come to this section of the code for all events except
   ; an UP button event. Most of the action in this event handler
   ; occurs here.

  cmps_form_real_coords, dpu, ixoff, iyoff, ixsize, iysize

  ; Simply draw the new box
  cmps_form_draw_box, ixsize, iysize, ixoff, iyoff, info

  ; Put the info structure back into the top-level base.
  Widget_Control, event.top, Set_UValue=info, /No_Copy

END ;*******************************************************************


Pro cmps_form_Grow_Box, event

   ; This event handler is summoned when a RIGHT button is clicked
   ; in the draw widget. It allows the user to draw the outline of a
   ; box with the mouse. It will continue drawing the new box shape
   ; until an UP event is detected. Then it will set the event handler
   ; back to cmps_form_Box_Events and turn cmps_form_Draw_Motion_Events to OFF.

   ; Get the info structure out of the top-level base.

  Widget_Control, event.top, Get_UValue=info, /No_Copy

  whatButtonType = cmps_form_What_Button_Type(event)

  dpu = info.drawpixperunit

  ixmin = 0.
  iymin = 0.
  ixsize = info.devconfig.xsize
  iysize = info.devconfig.ysize
  cmps_form_draw_coords, dpu, ixmin, iymin, ixsize, iysize
  ; Now ixmin,iymin have the minimum values of x and y, in pixels
  ; ixsize and iysize are the size of the box, in pixels

  ixmax = info.xpagesize
  iymax = info.ypagesize
  cmps_form_draw_coords, dpu, ixmax, iymax
  ; ixmax and iymax are the max values of x and y, in pixels

  ; Keep box inside the page
  if event.x LT ixmin then event.x = ixmin
  if event.x GT ixmax then event.x = ixmax
  if event.y LT iymin then event.y = iymin
  if event.y GT iymax then event.y = iymax

  ; Decide on which corner is the lower left (it's arbitrary)
  ixoff  = min([info.imousex, event.x])
  iyoff  = min([info.imousey, event.y])
  ixsize = max([info.imousex, event.x]) - ixoff
  iysize = max([info.imousey, event.y]) - iyoff

  ;; Enforce the aspect ratio
  if info.preserve_aspect EQ 1 then begin
      sizeratio = info.aspect / (info.ypagesize / info.xpagesize)
      if (sizeratio GE 1) then ixsize = iysize / info.aspect $
      else iysize = ixsize * info.aspect
      if info.imousex GT event.x then ixoff = info.imousex - ixsize
      if info.imousey GT event.y then iyoff = info.imousey - iysize
  endif

  IF whatButtonType EQ 'UP' THEN Begin

      ; When the button is "up" the moving event is over.  We reset the
      ; event function and update the information about the box's position

      Widget_Control, info.iddraw, Draw_Motion_Events=0, $ ; Motion events off
        Event_Pro='cmps_form_Box_Events' ; Change to normal processing

      cmps_form_real_coords, dpu, ixoff, iyoff, ixsize, iysize

      ; Update the info structure
      cmps_form_update_info, info, xoff=ixoff, yoff=iyoff, $
        xsize=ixsize, ysize=iysize
      ; Draw it
      cmps_form_draw_form, info

      ; Put the info structure back in the top-level base and RETURN
      Widget_Control, event.top, Set_UValue=info, /No_Copy
      Return
      
  ENDIF

   ; This is the portion of the code that handles all events except for
   ; UP button events. The bulk of the work is done here. Basically,
   ; you need to erase the old box and draw a new box at the new
   ; location. Just keep doing this until you get an UP event.

  cmps_form_real_coords, dpu, ixoff, iyoff, ixsize, iysize

  ; Simply draw the new box
  cmps_form_draw_box, ixsize, iysize, ixoff, iyoff, info

  ; Put the info structure back in the top-level base.
  Widget_Control, event.top, Set_UValue=info, /No_Copy

END ;*******************************************************************



; 
; Buttondown events sent to this procedure at first.  This is sets up
; the initial move/drag elements and hands off the events to the more
; specialized procedures cmps_form_grow_box and cmps_form_move_box above.
;
Pro cmps_form_Box_Events, event

  whatButtonType = cmps_form_What_Button_Type(event)
  IF whatButtonType NE 'DOWN' THEN Return

  ; Get info structure out of TLB

  Widget_Control, event.top, Get_UValue=info, /No_Copy

  whatButtonPressed = cmps_form_What_Button_Pressed(event)

  dpu = info.drawpixperunit
  ixmin = 0.
  iymin = 0.
  ixsize = info.devconfig.xsize
  iysize = info.devconfig.ysize
  cmps_form_draw_coords, dpu, ixmin, iymin, ixsize, iysize
  ixmax = info.xpagesize
  iymax = info.ypagesize
  cmps_form_draw_coords, dpu, ixmax, iymax
  ixoff = info.devconfig.xoff
  iyoff = info.devconfig.yoff
  cmps_form_draw_coords, dpu, ixoff, iyoff

  if event.x LT ixmin OR event.x GT ixmax $
    OR event.y LT iymin OR event.y GT iymax then begin
      widget_control, event.top, set_uvalue=info, /no_copy
      return
  endif
  
  CASE whatButtonPressed OF
      
      'RIGHT': Begin
          
         ; Resize the plot box interactively. Change the event handler
         ; to cmps_form_Grow_Box. All subsequent events will be handled by
         ; cmps_form_Grow_Box until an UP event is detected. Then you will
         ; return to this event handler. Also, turn motion events ON.
          
          Widget_Control, event.id, Event_Pro='cmps_form_Grow_Box', $
            Draw_Motion_Events=1
          
          cmps_form_draw_box, 1./dpu, 1./dpu, ixoff, iyoff, info
          
          info.imousex = event.x
          info.imousey = event.y
          
      End
      
      'LEFT': Begin
          
         ; Resize the plot box interactively. Change the event handler
         ; to cmps_form_Move_Box. All subsequent events will be handled by
         ; cmps_form_Move_Box until an UP event is detected. Then you will
         ; return to this  event handler. Also, turn motion events ON.


         ; Only move the box if the cursor is inside the box.
         ;If it is NOT, then RETURN.

          if event.x LT ixoff OR event.x GT (ixoff+ixsize) OR $
            event.y LT iyoff OR event.y GT (iyoff+iysize) then begin 
              
              Widget_Control, event.top, Set_UValue=info, /No_Copy
              Return
          ENDIF
          
       ; Relocate the event handler and turn motion events ON.

          Widget_Control, event.id, Event_Pro='cmps_form_Move_Box', $
            Draw_Motion_Events=1
          
          ; ideltx and idelty contain the offset of the lower left 
          ; corner of the plot region with respect to the mouse.
          info.ideltx  = ixoff - event.x
          info.idelty  = iyoff - event.y
          
      End

      ELSE:                     ; Middle button ignored in this program

  ENDCASE

   ; Put the info structure back into the top-level base

  Widget_Control, event.top, Set_UValue=info, /No_Copy

END ;*******************************************************************

; 
; Handle events to the drop-list widgets, which contain predefined
; configurations.
;
pro cmps_form_predef_events, event

  name = tag_names(event, /structure_name)
  if strupcase(name) NE 'WIDGET_DROPLIST' then return

   ; Get the info structure out of the top-level base
  Widget_Control, event.top, Get_UValue=info, /No_Copy

  Widget_Control, event.id, Get_UValue=thislist

  ; Pre-read the values from the text fields
  Widget_Control, info.idfilename, Get_Value=filename
  cmps_form_update_info, info, filename=filename

  case thislist of
      'PAPER':  info.paperindex = event.index   ; Paper change 
      'PREDEF': begin
          old_filename = info.devconfig.filename         ; Keep old filename
          info.devconfig = info.predefined(event.index)  ; New config
          info.paperindex = info.papersizes(event.index) ; New paper too
          if info.filechanged then $
            info.devconfig.filename = old_filename $
          else begin
              cd, current=thisdir
              l = strlen(thisdir)
              if strmid(info.devconfig.filename, 0, l) NE thisdir then $
                info.devconfig.filename = old_filename $
              else $
                info.devconfig.filename = filepath(info.devconfig.filename, $
                                                   root_dir=thisdir)
          endelse
      end
  endcase

  ; Be sure to select a pristine set of paper
  cmps_form_select_papersize, info.paperindex, xpagesize, ypagesize, $
    landscape=info.devconfig.landscape, inches=info.devconfig.inches
  info.xpagesize = xpagesize
  info.ypagesize = ypagesize

  widget_control, info.idpaperlist, set_droplist_select=info.paperindex

  ; Reset the drawpixperunit value
  convfactor = 1.0
  if info.devconfig.inches EQ 0 then convfactor = convfactor * 2.54
  info.marginsize = 0.25 * convfactor

  ; The conversion between length and pixels cannot always be set precisely,
  ; depending on the size of the paper
  dpp = 10.0 / convfactor  ; Desire 10 pixels per inch
  if dpp * info.xpagesize GT info.xpixwinsize OR $
    dpp * info.ypagesize GT info.ypixwinsize then $
    dpp = min( [ float(info.xpixwinsize-2)/info.xpagesize, $
                 float(info.ypixwinsize-2)/info.ypagesize ])
  info.drawpixperunit = dpp
  info.pixredraw = 1

  ; Update the info structure and draw it
  cmps_form_update_info, info, xoff=info.devconfig.xoff
  cmps_form_draw_form, info

  Widget_Control, event.top, Set_UValue=info, /No_Copy
  return
end

;
; Handle events sent to any of the button elements of the form.
;
Pro cmps_form_Event, event

   ; This is the main event handler for cmps_form. It handles
   ; the exclusive buttons on the form. Other events on the form
   ; will have their own event handlers.

   ; Get the name of the event structure

  name = Tag_Names(event, /Structure_Name)

   ; Get the User Value of the Button
  Widget_Control, event.id, Get_UValue=thisButton

   ; If name is NOT "WIDGET_BUTTON" or this is not a button
   ; selection event, RETURN.
  nonexclusive = ( thisButton EQ 'ISOLATIN1' OR $
                   thisButton EQ 'COLOR' OR $
                   thisButton EQ 'ENCAPSULATED' OR $
                   thisButton EQ 'ASPECT' )

  IF name NE 'WIDGET_BUTTON' OR $
     (NOT nonexclusive AND event.select NE 1)  THEN Return

   ; Get the info structure out of the top-level base
  Widget_Control, event.top, Get_UValue=info, /No_Copy


  redraw_form = 0
  redraw_box = 0

  ; Pre-read the values from the text fields
  Widget_Control, info.idxsize, Get_Value=xsize
  Widget_Control, info.idysize, Get_Value=ysize
  Widget_Control, info.idxoff, Get_Value=xoff
  Widget_Control, info.idyoff, Get_Value=yoff
  Widget_Control, info.idfilename, Get_Value=filename
  cmps_form_update_info, info, filename=filename

   ; Respond appropriately to whatever button was selected
  CASE thisButton OF

      'INCHES': Begin
          cmps_form_update_info, info, xsize=xsize, ysize=ysize, $
            xoff=xoff, yoff=yoff
          cmps_form_update_info, info, inches=1
          redraw_form = 1
      end
          
      'CENTIMETERS': Begin
          cmps_form_update_info, info, xsize=xsize, ysize=ysize, $
            xoff=xoff, yoff=yoff
          cmps_form_update_info, info, inches=0
          redraw_form = 1
      End
      
      'COLOR': Begin
          cmps_form_update_info, info, color=(1-info.devconfig.color)
          redraw_form = 1
      End

      'BITS2': Begin
          cmps_form_update_info, info, bits_per_pixel=2
          redraw_form = 1
      End

      'BITS4': Begin
          cmps_form_update_info, info, bits_per_pixel=4
          redraw_form = 1
      End

      'BITS8': Begin
          cmps_form_update_info, info, bits_per_pixel=8
          redraw_form = 1
      End

      'ISOLATIN1': Begin
          cmps_form_update_info, info, isolatin1=(1-info.devconfig.isolatin1)
       End

      'ASPECT': begin
          if info.preserve_aspect EQ 0 then $
            info.aspect = info.devconfig.ysize / info.devconfig.xsize
          info.preserve_aspect = (1 - info.preserve_aspect)
      end

      'LANDSCAPE': Begin
          cmps_form_update_info, info, xsize=xsize, ysize=ysize, $
            xoff=xoff, yoff=yoff
          cmps_form_update_info, info, landscape=1
          redraw_form = 1
          redraw_box = 1
      End

      'PORTRAIT': Begin
          cmps_form_update_info, info, landscape=0
          cmps_form_update_info, info, xsize=xsize, ysize=ysize, $
            xoff=xoff, yoff=yoff
          redraw_form = 1
          redraw_box = 1
      End

      'ENCAPSULATED': Begin
          cmps_form_update_info, info, encapsulated=(1-info.devconfig.encapsulated)
      End

      'ACCEPT': Begin

         ; The user wants to accept the information in the form.
         ; The procedure is to gather all the information from the
         ; form and then fill out a formInfo structure variable
         ; with the information. The formInfo structure is stored
         ; in a pointer. The reason for this is that we want the
         ; information to exist even after the form is destroyed.

         ; Gather the information from the form

          Widget_Control, info.idfilename, Get_Value=filename
          cmps_form_update_info, info, xsize=xsize, ysize=ysize, $
            xoff=xoff, yoff=yoff
          cmps_form_update_info, info, filename=filename
          widget_control, event.id, get_value=buttonname

          formInfo = { $
                       cancel:0, $            ; CANCEL flag
                       create:0, $            ; CREATE flag
                       buttonname: buttonname, $
                       xpagesize:info.xpagesize, $
                       ypagesize:info.ypagesize, $
                       paperindex:info.paperindex, $
                       result:info.devconfig $; Results are ready-made
                     }

          goto, FINISH_DESTROY
      End

      'CREATE': Begin

          Widget_Control, info.idfilename, Get_Value=filename
          cmps_form_update_info, info, xsize=xsize, ysize=ysize, $
            xoff=xoff, yoff=yoff
          cmps_form_update_info, info, filename=filename

          formInfo = { $
                       cancel:0, $            ; CANCEL flag
                       create:1, $            ; CREATE flag
                       buttonname: 'Create File', $ 
                       xpagesize:info.xpagesize, $
                       ypagesize:info.ypagesize, $
                       paperindex:info.paperindex, $
                       result:info.devconfig $; Results are ready-made
                     }

          goto, FINISH_DESTROY

      End

      'CANCEL': Begin

         ; The user wants to cancel out of this form. We need a way to
         ; do that gracefully. Our method here is to set a "cancel"
         ; field in the formInfo structure.
          formInfo = {cancel:1, create:0}

          goto, FINISH_DESTROY
      End

  ENDCASE

  if redraw_form EQ 1 then $
    cmps_form_draw_form, info, nobox=(1-redraw_box)

   ; Put the info structure back into the top-level base if the
   ; base is still in existence.

  If Widget_Info(event.top, /Valid) THEN $
    Widget_Control, event.top, Set_UValue=info, /No_Copy
  return


  ; We only reach this stage if we are ending the cmps_form widget
  ; These commands store the results, restore colors, and destroy
  ; the form widget.
  FINISH_DESTROY:

         ; Put the formInfo structure into the location pointer
         ; to by the pointer
  Handle_Value, info.ptrresult, formInfo, /Set, /No_Copy

         ; Delete the pixmap window
  WDelete, info.idpixwid

         ; Restore the user's color table
  TVLct, info.red, info.green, info.blue

         ; Destroy the cmps_form widget program
  Widget_Control, event.top, /Destroy

  return

END ;*******************************************************************

Function cmps_form, xoffset, yoffset, Cancel=cancelButton, Help=help, $
     XSize=xsize, YSize=ysize, XOffset=xoff, YOffset=yoff, $
     Inches=inches, Color=color, Bits_Per_Pixel=bits_per_pixel, $
     Encapsulated=encapsulated, Landscape=landscape, Filename=filename, $
     Defaults=defaults, LocalDefaults=localDefaults, Initialize=initialize, $
     select=select, parent=parent, $
     Create=createButton, NoCommon=nocommon, PaperSize=paperSize, $
     button_names=buttons, button_sel=button_sel, $
     PreDefined=predefined, DefaultPaper=defaultpaper, $
     aspect=aspect, preserve_aspect=preserve_aspect, $
     xpagesize=xpagesize, ypagesize=ypagesize, pagebox=pagebox

   ; If the Help keyword is set, print some help information and return

  IF Keyword_Set(help) THEN BEGIN
      Doc_Library, 'cmps_form'
      RETURN, 0
  ENDIF

  ; Set cancelButton and createButton as if canceled, so will be defined
  ; (and with appropriate values) even if user kills the window instead of
  ; using the buttons.  Normal exit will reassign them later on.
  cancelButton = 1
  createButton = 0

  ; Load default setups via a common block, if they are available
  if n_elements(predefined) EQ 0 then begin
      common cmps_form_configs, cmps_form_default_papersize, $
        cmps_form_stdconfigs
      if n_elements(cmps_form_stdconfigs) GT 0 then $
        predefined = cmps_form_stdconfigs
  endif

  ; If the user has not set up a common block, then get some pre
  if n_elements(predefined) EQ 0 then $
    cmps_form_load_configs, cmps_form_default_papersize, predefined

  ; Transfer to local copies so that we don't overwrite
  confignames = predefined(*).configname
  configs     = predefined(*).config
  configs     = configs(*) ;; IDL 5.5 will make a 1xN array -- collapse it now
  papernames  = predefined(*).papersize
  if n_elements(defaultpaper) EQ 0 $
    AND n_elements(cmps_form_default_papersize) GT 0 then $
    defaultpaper = cmps_form_default_papersize
  if n_elements(defaultpaper) EQ 0 then $
    defaultpaper = 'Letter'

  papersizes = intarr(n_elements(papernames))

  ; If localdefaults exist, then enter them into a new first entry of 
  ; the configuration list
  if n_elements(localDefaults) NE 0 then begin
      configs     = [ configs(0), configs ]
      confignames = [ 'Local',    confignames ]
      papernames  = [defaultpaper, papernames ]
      papersizes  = [ 0,          papersizes ]

      tmpc = configs(0)
      struct_assign, localdefaults, tmpc, /nozero
      configs(0) = tmpc
  endif


  ; Generate a new entry at the beginning, which will be the initial, 
  ; default configuration.
  configs     = [ configs(0), configs ]
  confignames = [ 'Default',  confignames ]
  papernames  = [defaultpaper, papernames ]
  papersizes  = [ 0,          papersizes ]

  filechanged = 0
  defaultset = 0
  if n_elements(defaults) NE 0 then begin
      defaultset = 1
      tmpc = configs(0)
      struct_assign, defaults, tmpc, /nozero
      configs(0) = tmpc
      void = where( strupcase(Tag_Names(defaults)) EQ 'FILENAME', count )
      if (count NE 0) then filechanged = 1
  endif

  ; Next, enter in the keyword defaults
  IF NOT defaultset OR N_ELEMENTS(inches) GT 0 then begin
      if n_elements(inches) EQ 0 then inches = 1
      configs(0).inches    = keyword_set(inches)
  endif
  IF NOT defaultset OR n_ele