c This file contains the FORTRAN 77 source code for HAMSCATT
c-----------------------------------------------------------------------
c
c
c             HAMILTON BACKGROUND ILLUMINATION CORRECTION 
c                             (HAMSCATT)
c
c                         Copyright (c) 1994 
c                       Christopher W. Churchill
c                 Board of Astronomy and Astrophysics
c                 University of California, Santa Cruz
c
c **********************************************************************
c *  TERMS OF USAGE:  See the COPYRIGHT file in this directory         *
c *                   do not modify text contained in this box         *
c *                                                                    *
c *  The usage of this software is free, but any published statements  *
c *  based upon usage of this software must reference the paper        *
c *  Churchill, C.W., and Allen, S.L. 1995, PASP, 107, 193             *
c *                                                                    *
c *  This code was originally packaged for distribution from:          *
c *  internet URL=http://ucowww.ucsc.edu/~cwc)                         *
c *  anon ftp URL=ftp://lick.ucsc.edu/pub/outgoing/cwc                 *
c *                                                                    *
c **********************************************************************
c A brief accompanying usage manual is available, see the INSTALL file
c and the README and README2 files
c
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c  
c  
c  
      program             HAMSCATT
c  
c  
c.......................................................................
c  
c
c  author: Christopher W. Churchill
c
c  history: 20-Jan94   created, CWC
c           02-Oct94   moved error/warning messages to subroutines, CWC
c           see individual histories of subroutines
c
c  purpose: background illumination correction of echelle data
c           (1) input model parameters, aperture locations,
c               and input image
c           (2) grid the data for the least squares fitter
c           (3) mask bad pixels (sigma rejection)
c           (4) call to the least squares fitter
c           (5) inform user of convergence criterion satsified
c           (6) compute Chi-Sq, fit variance, data variance, etc.
c           (7) construct the background surface
c           (8) output a background image
c
c  input:   (1) input run file with fitting parameters and flags
c           (2) apmap file with order maps
c           (3) input image to be background corrected
c
c  output:  (1) fitting parameters 
c           (2) background illumination image (IRAF)          
c           (3) data quality file '.dqf' (rejected pixels)
c           (4) fitting file '.fit' (fitting parameters and stats)
c
c  rquired files: 
c
c        hamscatt.par - the parameter file, defines working space memory
c        hamscatt.com - the common blocks, defines common variables
c
c        fcn.f      routines initsol, fcn, and poly
c        gridim.f   routine  gridim
c        io.f       routines inputdeck, ordermap, inform, and output
c        iraf.f     routines inputimag, outputimag, and alarm
c        linlsf.f   routines linlsf, lsfsolve, and function lsffunc
c        mask.f     routines mask and checks
c        results.f  routines scatsurf and fitstats
c
c  linked calls: dnl1se:
c
c  required libraries: IMFORT, SLATEC (partial)
c
c  description: 
c
c  2D two-component scattered light model for echelle spectrographs
c  based upon studies of the UCO/Lick Hamilton spectrograph conducted 
c  between Sep93 - Jan94
c
c  this work is based upon the formalism introduced by Gehren and Ponz
c  Astronomy and Astrophysics _168_ 386 (1986).  
c
c  the model attempts to recover the inter-order minima by modeling 
c  these data as arising from both 
c           (1) a global scattered light surface 
c           (2) a local lambda dependent component from 
c               the nearest neighbor orders
c
c  the resulting solution is a 2D surface with local "undulations"
c  from the contamination of the adjacent orders.  the local component
c  assumes symmetric tails for the cross dispersion illumination
c  functions that follow a power law.  this power law is an adjustable 
c  parameter.  To treat asymmetric profiles, the inter-order minima of 
c  each order can be found using an "optimization" scheme, which, for 
c  each order, finds a weighted mean fractional position between the 
c  adjacent order tracing solutions of IRAF.  the tracing solutions are
c  output from task _apmap_, which is a hacked version of _ecmap_, a 
c  Lick local task written by Frank Valdes.
c
c  the polynomial orders are adjustable input parameters, denoted by
c  xorder and yorder.  the code determines the n coefficients; n-1 for 
c  the 2D polynomial and 1 for the local term.  using the notation 
c  x=row# and y=column#, where y2=y**2 etc, the coefficients of the
c  terms for a xorder=3 yorder=4 polynomial would be 
c 
c             a1 = const      a2 = x          a3 = x2      
c             a4 = y          a5 = y2         a6 = y3      
c             a7 = x*y        a8 = x*y2       a9 = x*y3
c            a10 = x2*y      a11 = x2*y2     a12 = x2*y3
c           a13 = local
c
c  the work-horse code is a multi-dimensional non-linear least squares
c  modified Levenberg-Marquardt algorithm called dnls1.f.  here, we use
c  the upper level driver called dnls1e.f. with the merit function
c  taking the form of chi-square, thus, the solution is a maximum
c  likelihood result.
c
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c  
      implicit undefined (a-z)
      logical            error
      external           fcn
      include            'hamscatt.par'
      integer            m,n,lwa,info,iw(nn),nrej
      parameter         (lwa = nn*(nvec+5)+nvec)
      double precision   a(nn),fvec(nvec),wa(lwa),tol,chisq,chisqnu,
     @                   varfit,vardata,nu
      character*80       inimage,outimage,xmaxfile
      include            'hamscatt.com'
c  
c  
c  
c
c  set the input error flag low
      error = .false.
c  
c  I/O routine:
c  load input deck ; error trap ; return number of fit coefficients 
c  and the convergence tolerance
c  if an error has been flagged then ABORT
      call inputdeck(inimage,outimage,xmaxfile,error)
      if (error) then
       write(stdout,*) ' ERROR: (routine inputdeck)'
       stop ' abort routine HAMSCATT'
      end if
c
c  I/O routine:
c  load image ; returns pix values, ncols, and nrows ; error trap 
c  if an error has been flagged then ABORT
      call inputimag(inimage,error)
      if (error) then
       write(stdout,*) ' ERROR: (routine inputimag)'
       stop ' abort routine HAMSCATT'
      end if
c
c  I/O routine:
c  load the order map ; error trap; return ndata and the number of 
c  functions to minimize
c  if an error has been flagged then ABORT 
      call ordermap(xmaxfile,error)
      if (error) then
       write(stdout,*) ' ERROR: (routine ordermap)'
       stop ' abort routine HAMSCATT'
      end if
c
c  set the physical spaces based upon inputs and perform some
c  consitency checks ; the info=1 error should have been flagged
c  when the image was opened, but just in case we check it again
c  if an error has been flagged then ABORT
      call checks(n,m,tol,error,info)
      if (error) then
       write(stdout,*) ' ERROR: (routine checks)'
       if (info.eq.1) then
        write(stdout,*) ' physical dimensions of working space arrays'
        write(stdout,*) ' too small for image.  increase physical'
        write(stdout,*) ' dimension parameters MAXCOLS and/or MAXROWS'
        write(stdout,*) ' and re-compile.'
       end if
       if (info.eq.2) then
        write(stdout,*) ' polynomial orders too high for physical'
        write(stdout,*) ' dimensions of working space arrays. '
        write(stdout,*) '(1) lower orders and re-execute or '
        write(stdout,*) '(2) increase physical parameter NN '
        write(stdout,*) '    and re-compile.'
       end if
       write(stdout,*) ' modifications are made in file: hamscatt.par'
       stop ' abort routine HAMSCATT'
      end if
c
c  inputs check OK
c  stuff fitting arrays ; logical input parameter minflag dictates 
c  scheme used to calculate xmin, Imin, and Imax ; stuff initial 
c  "solution" a(i=1,...,n) ; if reject then mask bad pixels
      call gridim(n,a)
      if (reject.ne.0.0d0) then
       call mask(a,wa,m,error,info)
       if (error) then
        write(stdout,*) ' ERROR:  (routine mask)'
        if (info.eq.1) write(stdout,*) ' N > MAX in routine lsfsolve' 
        if (info.eq.2) write(stdout,*) ' pivoting   - singular matrix'
        if (info.eq.3) write(stdout,*) ' diagonal=0 - singular matrix'
        stop ' abort routine HAMSCATT'
       end if
      end if
c  
c  obtain an initial estimate for the solution, the hope is to
c  significantly reduce run-time by reducing time spent in 
c  routine fcn
      call initsol(a,wa,error,info)
      if (error) then
       write(stdout,*) ' ERROR: (routine initsol)'
       if (info.eq.1) write(stdout,*) ' N > MAX in routine lsfsolve' 
       if (info.eq.2) write(stdout,*) ' pivoting   - singular matrix'
       if (info.eq.3) write(stdout,*) ' diagonal=0 - singular matrix'
       stop ' abort routine HAMSCATT'
      end if
c
c  call the least-squares solver ; communicate convergence criterion
c  if an error has occured then ABORT
      write(stdout,*) ' [proceeding with the fit... ]'
      call dnls1e(fcn,2,m,n,a,fvec,tol,0,info,iw,wa,lwa)
      call inform(info,error)
      if (error) then
       write(stdout,*) ' ERROR: (routine dnls1e)'
       stop ' abort routine HAMSCATT'
      end if
c
c  fit is successful!   
c
c  compute fit statistics
      call fitstats(m,n,a,chisq,nu,chisqnu,varfit,vardata,nrej)
      call output(chisq,nu,chisqnu,varfit,vardata,m,a,n,nrej,outimage)
c  
c  compute scattered light surface
      call scatsurf(n,a,wa)
c  
c  output scattered light surface image ; the noclobber irafism is 
c  checked in the subroutine
      call outputimag(inimage,outimage,error)
      if (error) then
       write(stdout,*) ' ERROR: (routine outputimag)'
       stop ' abort routine HAMSCATT'
      end if
c  
c  we are done
      stop ' RUN STATUS = normal termination of routine HAMSCATT'
      end
c
c.......................................................................
c
      include 'fcn.f'
      include 'gridim.f'
      include 'io.f'
      include 'iraf.f'
      include 'linlsf.f'
      include 'mask.f'
      include 'results.f'
c
c.......................................................................
c..eof