c  this file contains subroutines initsol, fcn, and poly
c
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
c
c
      subroutine         initsol(a,wa,error,etype)
c
c
c
c.......................................................................
c
c  author: Christopher W. Churchill
c
c  history: created 30-Aug94
c 
c  purpose: obtain an initial estimate of the dispersion and cross
c           dispersion polynomial coefficients, not including cross 
c           terms ; hopefully this reduces the number of calls to fcn
c
c  input:  a     =  fitting coefficients 
c          wa    =  working array
c
c  output: a     =  fitting coefficients (inital estimate)
c          wa    =  fitting residuals
c          error = (logical) error flag
c          etype =  error type
c
c  I/O : none
c
c  options: none
c
c  routines called: linlsf
c
c  description: see subroutine comments
c
c.......................................................................
c
      implicit undefined (a-z)
      logical            error
      integer            i,line,etype
      double precision   a(1),wa(1),chisq,fitsig
      include            'hamscatt.par'
      include            'hamscatt.com'
c
c
c
c
c
c  obtain dispersion direction cofficients along the center order
c  on return the first yorder a contain the initial guess
c  these will need to be "moved" to their proper index 
c  do not include the constant term, which we will obtain
c  from the XD solution
c  on error, line contains the error type
      line = ndata/2
      call linlsf(line,a,wa,chisq,fitsig,error,2)   
      if (error) then
       etype = line
       return
      end if
      do 01 i=1,yorder-1
       a(xorder+i) = a(i+1)
 01   continue
c
c  obtain XD coefficients along the center column 
c  on return the first xorder a contain the initial guess
c  on error, line contains the error type
      line = ncols/2
      call linlsf(line,a,wa,chisq,fitsig,error,1)   
      if (error) then
       etype = line
       return
      end if
c
c  the cross terms are still zero
      return
      end
c  
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c  
c  
c  
      subroutine         fcn(iflag,m,n,a,fvec,fjac,ldfjac)
c  
c  
c  
c.......................................................................
c  
c  author: Christopher W. Churchill
c
c  history: 20-Jan94   created, CWC
c           20-Feb94   local term only option added, CWC
c
c  purpose: update the merit function or the jacobian for the next
c           iteration to the least squares fitter DNLS1 (slatec)
c
c  input:   m        = number of fitting functions
c           n        = number of fitting parameters
c           a        = array of fitting parameters
c           ldfjac   = (hidden slatec parameter) = m
c           xmax     = order map [pixel locations] (APEXTRACT tracings)
c           xmin     = interorder map [pixel locations]
c           Imin     = 2D grid to fit [e-] at xmin locations
c           Imax     = order Maxima [e-] at xmax locations
c           sig      = Imin uncertainty ; fitting weights
c           gamma    = power law of cross dispersion profile 
c
c  output:  fvec     = m dimensional merit function
c           fjac     = m by n dimensional jacobian
c
c  I/O: none
c
c  options: iflag    = (integer) flag determines action taken
c             1      = update fvec = (ymodel - ydata)/sig
c             2      = update fjac = dfvec/da
c           termflag = (logical) flag for local term inclusion to fit
c           .true.   = include local term in fit
c           .false.  = exclude local term from fit
c
c
c  routines called: poly
c 
c  description: fcn is the user-supplied subroutine which calculates 
c               the functions for the least squares fitter DNLS1.  
c               DNLS1 is a modification of the  Levenberg-Marquardt 
c               algorithm.  Two of its main characteristics involve the
c               proper use of implicitly scaled variables and an optimal
c               choice for the correction.  The use of implicitly scaled
c               variables achieves scale invariance of DNLS1 and limits
c               the size of the correction in any direction where the 
c               functions are changing rapidly.  The optimal choice of 
c               the correction guarantees (under reasonable conditions)
c               global convergence from starting points far from the 
c               solution and a fast rate of convergence for problems 
c               with small residuals.
c               
c  
c.......................................................................
c  
      implicit undefined (a-z)
      integer            i,j,k,ix,iy,icol
      integer            m,n,ldfjac,iflag
      double precision   a(1),fvec(1),fjac(ldfjac,n)
      double precision   xa,xb,xc,y,pa,pb,pc,phi,poly,ymod
      include            'hamscatt.par'
      include            'hamscatt.com'
c  
c  
c  
c  
c  
c
c  zero the vector index and loop over the number of interorder
c  minima (j) and the number of columns (icol)
      i  = 0
c
c
      if (termflag) then
c
c  --- begin local + global component ---
c
c  local + global component function computation ; check for rejected 
c  pixels ; increment the index and compute the polynomial at the 
c  interorder and adjacent orders; stuff the vector to be minimized
c
       if (iflag.eq.1) then
        do 07 j=1,ndata
         do 06 icol=1,ncols
          if (sig(j,icol).ne.0.) then 
           i    = i + 1  
           y    = real(icol)
           xa   = xmin(j,icol)
           xb   = xmax(j,icol)
           xc   = xmax(j+1,icol)
           pa   = poly(n,a,xa,y)
           pb   = poly(n,a,xb,y)
           pc   = poly(n,a,xc,y) 
           phi  = (0.5*(xc-xb))**gamma
           ymod = pa + a(n)*((Imax(j,icol)-pb)+(Imax(j+1,icol)-pc))/phi
           fvec(i) = (Imin(j,icol) - ymod)/sig(j,icol)
          end if
 06      continue
 07     continue
        return
       end if
c
c  local + global component jacobian computation ; check for rejected 
c  pixels
       if (iflag.eq.2) then
        do 09 j=1,ndata
         do 08 icol=1,ncols
          if (sig(j,icol).ne.0.) then 
           i   = i + 1
           y   = real(icol)
           xa  = xmin(j,icol)
           xb  = xmax(j,icol)
           xc  = xmax(j+1,icol)
           phi = (0.5*(xc-xb))**gamma
           k  = xorder
           do 13 ix=1,xorder
            fjac(i,ix) = ( a(n)*(xb**(ix-1)+xc**(ix-1))/phi - 
     @                   (xa**(ix-1)) ) / sig(j,icol)
            k          = k + 1
            fjac(i,k)  = y*fjac(i,ix)
            do 14 iy=2,yorder-1
             k         = k + 1
             fjac(i,k) = y*fjac(i,k-1)
 14          continue
 13        continue
           fjac(i,n) = ((pb-Imax(j,icol)) + (pc-Imax(j+1,icol)))
     @               / (phi*sig(j,icol))
          end if
 08      continue
 09     continue
        return
       end if
c
c  --- end local + global component ---
c
      else
c
c  --- begin global component only ---
c
c  global component only function computation ; check for rejected pixel
c  increment the index and compute the polynomial at the interorder
c  stuff the vector to be minimized
       if (iflag.eq.1) then
        do 03 j=1,ndata
         do 02 icol=1,ncols
          if (sig(j,icol).ne.0.) then 
           i    = i + 1  
           y    = real(icol)
           xa   = xmin(j,icol)
           ymod = poly(n,a,xa,y)
           fvec(i) = (Imin(j,icol) - ymod)/sig(j,icol)
          end if
 02      continue
 03     continue
       end if
c
c  global component only jacobian computation ; check for rejected pixel
       if (iflag.eq.2) then
        do 05 j=1,ndata
         do 04 icol=1,ncols
          if (sig(j,icol).ne.0.) then 
           i  = i + 1
           y  = real(icol)
           xa = xmin(j,icol)
           k  = xorder
           do 11 ix=1,xorder
            fjac(i,ix) = -(xa**(ix-1))/sig(j,icol)
            k          = k + 1
            fjac(i,k)  = y*fjac(i,ix)
            do 12 iy=2,yorder-1
             k         = k + 1
             fjac(i,k) = y*fjac(i,k-1)
 12         continue
 11        continue
          end if
 04      continue
 05     continue
        return
       end if
c
c  --- end local component only ---
c
      end if
      end
c  
c.......................................................................
c  
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c  
c  
c  
      double precision   function poly(n,a,x,y)
c  
c  
c  
c.......................................................................
c  
c  author: Christopher W. Churchill
c
c  history: 20-Jan94   created, CWC
c
c  purpose: evaluate the polynomial (global component) at the x,y pixel
c           position
c
c  input:   n      = number of fitting parameters
c           a      = array of fitting parameters
c           x      = pixel location in cross dispersion direction
c           y      = pixel location in dispersion direction
c           xorder = order of x dependence
c           yorder = order of y dependence
c
c  output:  poly   = the value of the polynomial at x,y
c
c  I/O: none
c
c  options: none
c
c  routine called: none
c
c  description: the polynomial evaluated at the input x,y position ; 
c               this coding is for a generalized polyomial of xorder by 
c               yorder, the first xorder coefficients are x terms, the 
c               next yorder-1 coefficients are y terms, the following
c               xorder*yorder-(xorder+yorder-1) terms are the cross 
c               terms evaluated such the fitting parameter index 
c               increments with the y degree dependence varying for 
c               each step in x degree dependence
c
c.......................................................................
c  
      implicit undefined (a-z)
      integer            n,i,j,k
      double precision   a(1),x,y,p,poly
      include            'hamscatt.par'
      include            'hamscatt.com'
c  
c  
c  
c  
c  
c  set constant term ; sum x (row) terms ; 
c  sum y (column) terms ; sum cross terms
      p = a(1)
      do 01 i=2,xorder
       p = p + a(i)*(x**(i-1))
 01   continue
      do 03 j=1,yorder-1
       i = xorder + j 
       p = p + a(i)*(y**j)
 03   continue
      i = xorder + yorder  
      do 05 k=1,xorder-1
       do 04 j=1,yorder-1
        p = p + a(i)*(x**k)*(y**j)
        i = i + 1
 04    continue
 05   continue
c  
      poly = p
c  
c  
      return
      end
c
c.......................................................................
c..eof