c  this file contains subroutine linlsf, lsfsolve, and function lsffunc
c
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
c
c
      subroutine linlsf(line,a,wa,chisq,fitsig,error,dir)   
c
c
c
c.......................................................................
c  
c  author: Christopher W. Churchill
c
c  history: 20-Feb94   created, CWC
c
c  purpose:  to compute a least squares fit for to 1D polynomial that is
c            evaluated in subroutine lsffunc(x,fx,ncoeffs) ; on output 
c            wa contains the fit residuals
c
c  input:  line   = the column of order to be fit
c          dir    = 1=XD column 2=dispersion interorder
c
c  output: a      = the polynomial coefficients
c          wa     = working array, contains fit residuals
c          chisq  = chi-squared of the fit
c          fitsig = standard deviation of the fit
c          error  = (logical) error flag
c
c  I/O: none
c
c  options: dir flags column or interorder fitting
c
c  routines called: lsfsolve
c
c  description: see subroutine comments
c
c.......................................................................
c
      implicit undefined (a-z)
      logical            error
      external           lsffunc
      integer            line,npts,na,max,i,j,k,dir,nrej
      parameter          (max=15)   
      double precision   x,y,sigy,a(1),wa(1),normeq(max),fx(max),
     @                   wrkmtrx,chisq,weight,sigbar,fitsig,ymod
      dimension          wrkmtrx(max,max)
      include            'hamscatt.par'
      include            'hamscatt.com'
c  
c  
c
c
c
      na   = 0
      npts = 0
      nrej = 0
c
c  cross dispersion initialization
      if (dir.eq.1) then
       na   = xorder
       npts = ndata
      end if
c
c  dispersion initialization
      if (dir.eq.2) then
       na    = yorder
       npts  = ncols
      end if
c
c  just be sure we input dir = 1 or 2
      if (na.eq.0) then
       error = .true.
       return
      end if
c
c  initialize the work matrix and zero the normal equation coeffs
      do 1 j=1,na
       do 2 k=1,na  
         wrkmtrx(j,k) = 0.0d0 
 2     continue
       normeq(j) = 0.0d0  
 1    continue  
c  
c  stuff the work matrix and the coeffs of the normal equations
c  x, y, and sig are dir dependent
      do 5 i=1,npts   
       if (dir.eq.1) then
        x    = xmin(i,line)
        y    = Imin(i,line)
        sigy = sig(i,line)
       else
        x    = float(i)
        y    = Imin(line,i)
        sigy = sig(line,i)
       end if
       call lsffunc(x,fx,na)   
c
c  check for rejected pixels, vanish their weight 
       do 4 j=1,na  
        if (sigy.eq.0) then
         weight = 0.0d0
         nrej   = nrej + 1
        else
         weight = fx(j)/(sigy**2)
        end if
        do 3 k=1,j   
         wrkmtrx(j,k) = wrkmtrx(j,k) + weight*fx(k)
 3      continue  
        normeq(j) = normeq(j) + weight*y
 4     continue
 5    continue  
c  
c  mirror the work matrix about diagonal
       do 7 j=2,na  
        do 6 k=1,j-1 
         wrkmtrx(k,j) = wrkmtrx(j,k)   
 6      continue   
 7     continue
c  
c  solve for the coefficients, upon return the normal equation coeffs 
c  are the fit coeffs, copy them over 
c  if an error has occured, replace line with the error type and return
      call lsfsolve(wrkmtrx,normeq,na,error,line)
      if (error) return
      do 8 j=1,na
       a(j) = normeq(j) 
 8    continue  
c  
c  compute reduced chi square, stuff the working array with
c  the fit residuals, compute the fit standard deviation 
      chisq  = 0.0d0  
      sigbar = 0.0d0
      fitsig = 0.0d0
      do 10 i=1,npts   
       if (dir.eq.1) then
        x    = xmin(i,line)
        y    = Imin(i,line)
        sigy = sig(i,line)
       else
        x    = float(i)
        y    = Imin(line,i)
        sigy = sig(line,i)
       end if
       ymod = 0.0d0  
       call lsffunc(x,fx,na)   
       do 9 j=1,na
        ymod = ymod + a(j)*fx(j) 
 9     continue
       wa(i)  = y - ymod
       if (sigy.ne.0) then
        chisq  = chisq + (wa(i)/sigy)**2 
        sigbar = sigbar + 1.0d0/(sigy**2)
       end if
 10   continue  
      sigbar = sigbar/(npts-nrej)
      fitsig = sqrt(chisq/(sigbar*(npts-na-nrej)))
c  
      return
      end   
c
c.......................................................................
c
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
c
c
      subroutine lsfsolve(wrkmtrx,normeq,n,error,etype)
c
c
c
c.......................................................................
c  
c  author: Christopher W. Churchill
c
c  history: 20-Feb94   created, CWC
c
c  purpose: classical Gauss-Jordan elimination with pivoting, on output
c           wrkmtrx is destroyed and normeq is replaced by the fit 
c           coefficents
c
c  input:  wrkmtrx = coefficients of normal equations
c          normeq  = RHS normal equations
c          n       = number of coefficients, number of normal equations
c
c  output: normeq  = solution vector, the polynomial coefficients
c          error   = (logical) error flag
c          etype   = error type
c
c  I/O: none
c
c  options: none
c
c  routines called: none
c
c  description: see subroutine comments, three possible error are
c               flagged (1) exceed local working space 
c                       (2) pivoting error - singular matrix
c                       (3) zero diagonal  - singular matrix
c
c  NB. wrkmtrx dimension - must equal physical dimension in linlsf
c      ie. max in linlsf = max in lsfsolve
c
c.......................................................................
c  
      implicit undefined (a-z)
      logical            error
      integer            n,max,i,j,k,irow,icol,etype
      parameter          (max=15)
      integer            pivot(max)
      double precision   wrkmtrx,normeq,colmax,tmp,pivinv
      dimension          wrkmtrx(max,max),normeq(1)
      include            'hamscatt.par'
c
c
c
c
c
c  n must be less than max
      if (n.gt.max) then
       error = .true.
       etype = 1
       return
      end if
c  
c  initialize
      do 1 j=1,n
       pivot(j) = 0
 1    continue
c  
c  loop over columns
      do 11 i=1,n
       colmax = 0.0d0
c  
c  search for the pivot element
       do 3 j=1,n
        if (pivot(j).ne.1) then
         do 2 k=1,n
          if (pivot(k).eq.0) then
           if (abs(wrkmtrx(j,k)).ge.colmax) then
            colmax  = abs(wrkmtrx(j,k))
            icol    = k
            irow    = j
           endif
          else 
           if (pivot(k).gt.1) then
            error = .true.
            etype = 2
            return
           end if
          end if
 2       continue
        endif
 3     continue
       pivot(icol) = pivot(icol) + 1
c  
c  index the pivot element to be on the diagonal
        if (irow.ne.icol) then
         do 4 j=1,n
          tmp             = wrkmtrx(irow,j)
          wrkmtrx(irow,j) = wrkmtrx(icol,j)
          wrkmtrx(icol,j) = tmp
 4       continue
         tmp          = normeq(irow)
         normeq(irow) = normeq(icol)
         normeq(icol) = tmp
        endif
c  
c  divide the pivot row by the pivot element 
        if (wrkmtrx(icol,icol).eq.0.0d0) then
         error = .true.
         etype = 3
         return
        end if
        pivinv             = 1.0d0/wrkmtrx(icol,icol)
        wrkmtrx(icol,icol) = 1.0d0
        do 5 j=1,n
          wrkmtrx(icol,j) = wrkmtrx(icol,j)*pivinv
 5      continue
        normeq(icol) = normeq(icol)*pivinv
c  
c  reduce all the rows, barring the pivot row
       do 7 j=1,n
        if (j.ne.icol) then
         tmp           = wrkmtrx(j,icol)
         wrkmtrx(j,icol) = 0.0d0
         do 6 k=1,n
          wrkmtrx(j,k) = wrkmtrx(j,k) - tmp*wrkmtrx(icol,k)
 6       continue
         normeq(j) = normeq(j) - tmp*normeq(icol)
        endif
 7     continue
c  
c  next column
 11   continue
c  
      return
      end
c  
c.......................................................................
c
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
c
c
      subroutine         lsffunc(x,fx,na)
c
c
c
c.......................................................................
c
c this is the subroutuine for linlsf, it returns the functions for each
c of the polynomial coefficients
c
c
c.......................................................................
c
      implicit undefined (a-z)
      integer            na,i
      double precision   x,fx(1)
c
c
c
c
c
c  a polynomial of order na
      fx(1) = 1.0d0
      do 1 i=2,na
       fx(i) = x**(i-1)
  1   continue
c      
      return
      end
c
c.......................................................................
c..eof