c  this file contains subroutines scatsurf and fitstats
c  
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c  
c  
c  
      subroutine         scatsurf(n,a,wa)
c  
c  
c  
c.......................................................................
c  
c  author: Christopher W. Churchill
c
c  history: 20-Jan94   created, CWC
c           20-Feb94   global only option added, CWC
c
c  purpose: create the background illumination surface ; on output the
c           array pix is replaced with the background data
c
c  input:   n        = number of fitting parameters
c           a        = array of fitting parameters
c           wa       = working array for polynomial storage
c           pix      = pixel values [DN] of input image
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:  pix      = pixel values [DN] of background surface
c
c  I/O: none
c
c  options: corrflag = (logical) flag for local term correction
c           .true.   = background surface is global + local 
c           .false.  = background surface is global component only
c
c  routines called: poly
c
c  description: see subroutine comments
c
c.......................................................................
c  
      implicit undefined (a-z)
      integer            n,j,icol,irow,z1,z2
      double precision   a(1),x,xjm1,xjp1,fjm1,fjp1,y,poly,
     @                   phijm1,phijp1,f1,b2,wa(1)
      include            'hamscatt.par'
      include            'hamscatt.com'
c  
c  
c  
c  
c  
c  --- begin global + local term block ---
c
      if (corrflag) then
c  
c  loop over the pixels and compute the background, sample the inter-
c  order minima as a combined fraction by saving the upper limit 
c  position as b2; note that primary loop is over columns and not the 
c  usual orders, this simplifies book keeping for the pix(icol,irow) 
c  values at the actual inter-order minima
c
      do 11 icol=1,ncols
       y = real(icol)
c  
c  to minimize the calls to the polynomial, stuff the values of poly
c  evaluated at the maxmima for the current column into array p(j)
       do 09 j=1,norders
        x     = xmax(j,icol)
        wa(j) = poly(n,a,x,y)
 09    continue
c  
c  treat the first aperture, j=1
       j    = 1
       z2   = int(xmin(j,icol))
       xjp1 = xmax(j+1,icol)
       fjp1 = Imax(j+1,icol) - wa(j+1)
       do 15 irow=1,z2
         x      = real(irow)
         phijp1 = abs(x-xjp1)**gamma
         pix(icol,irow) = (poly(n,a,x,y) + a(n)*fjp1/phijp1)/gain
         if (irow.eq.z2) b2 = pix(icol,irow) 
 15    continue 
c  
c  treat the 2nd to the norders-1 apertures; 
       do 13 j=2,ndata
        z1   = int(xmin(j-1,icol))
        z2   = int(xmin(j,icol))
        f1   = xmin(j-1,icol) - real(z1)
        xjm1 = xmax(j-1,icol)
        xjp1 = xmax(j+1,icol)
        fjm1 = Imax(j-1,icol) - wa(j-1)
        fjp1 = Imax(j+1,icol) - wa(j+1)
        do 17 irow=z1,z2
         x            = real(irow) + f1
         phijm1       = abs(x-xjm1)**gamma
         phijp1       = abs(x-xjp1)**gamma
         pix(icol,irow) = (poly(n,a,x,y) + 
     @                  a(n)*(fjm1/phijm1+fjp1/phijp1))/gain
         if (irow.eq.z1) pix(icol,irow) = 
     @                   (1.0d0-f1)*pix(icol,irow) + f1*b2
         if (irow.eq.z2) b2 = pix(icol,irow) 
 17     continue
 13    continue
c  
c  treat the last aperture, j=norders
       j    = norders
       z1   = int(xmin(j-1,icol))
       f1   = xmin(j-1,icol) - real(z1)
       xjm1 = xmax(j-1,icol)
       fjm1 = Imax(j-1,icol) - wa(j-1)
       do 19 irow=z1,nrows
        x            = real(irow)
        phijm1       = abs(x-xjm1)**gamma
        pix(icol,irow) = (poly(n,a,x,y) + a(n)*fjm1/phijm1)/gain
        if (irow.eq.z1) pix(icol,irow) = 
     @                  (1.0d0-f1)*pix(icol,irow) + f1*b2
 19   continue 
 11   continue
c  
      end if
c
c  --- end global + local term block ---
c
c
c
c  --- begin global term only block ---
c
      if (.not.corrflag) then
c  
c  output the global component only as the background
c
      do 21 irow=1,nrows
       x = real(irow)
       do 23 icol=1,ncols
        y = real(icol)
        pix(icol,irow) = poly(n,a,x,y)/gain
 23    continue
 21   continue
c  
      end if
c  
c  --- end global term only block ---
c
c  
      return
      end
c  
c.......................................................................
c  
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c  
c  
c  
      subroutine         fitstats(m,n,a,chisq,nu,chisqnu,varfit,
     @                            vardata,nrej)
c  
c  
c  
c.......................................................................
c  
c  author: Christopher W. Churchill
c
c  history: 30-Aug94   created, CWC
c
c  input:  m       = number of fitted points
c          n       = number of fittd parameters (coefficients) 
c          a       = solution vector (coefficients)
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: chisq   = the fit Chi-Squared
c          nu      = number of degrees of freedom
c          chisqnu = reduced Chi-Squared
c          varfit  = variance in the fit
c          vardata = variance in the data (population variance)
c          nrej    = number of rejected pixels
c
c  I/O: none
c
c  options: termflag = (logical) output surface flag
c           .true.   = global + local term on output
c           .false.  = global term only on output
c
c  routines called: poly
c
c  description: see subroutine comments, this routine was included for
c               computations to help quantify the fits to various 
c               polynomial order and profile power-laws ; it may serve
c               the same purpose for a general user
c
c  NB. chisq and chisqnu are gain independent, varfit and vardata are
c      computed for data in electrons (statistics have little meaning 
c      in units of DN)
c
c.......................................................................
c  
      implicit undefined (a-z)
      integer            j,k,m,n,iord,nrej
      double precision   a(1),y,xa,xb,xc,pa,pb,pc,ymod,phi,poly
      double precision   nu,chisq,vardata,chisqnu,varfit
      include            'hamscatt.par'
      include            'hamscatt.com'
c  
c  
c  
c  
c
c  compute degrees of freedom and number of rejected pixels  
      nu   = m - n
      nrej = ndata*ncols - m
c
c  initialize
      j       = 0
      chisq   = 0.0d0
      vardata = 0.0d0
c  
c  compute chi-squared, this is painfully ugly, but easy to follow
      do 02 iord=1,ndata
       do 03 k=1,ncols
        if (sig(iord,k).ne.0.) then
         y  = real(k)
         xa = xmin(iord,k)
         pa = poly(n,a,xa,y)
         if (termflag) then
          xb   = xmax(iord,k)
          xc   = xmax(iord+1,k)
          phi  = (0.5*(xc-xb))**gamma
          pb   = poly(n,a,xb,y)
          pc   = poly(n,a,xc,y)
          ymod = pa + a(n)*((Imax(iord,k)-pb)+(Imax(iord+1,k)-pc))/phi
         else
          ymod = pa
         end if
         chisq   = chisq + ((ymod - Imin(iord,k))/sig(iord,k))**2
         vardata = vardata + 1.0d0/(sig(iord,k)**2)
        end if
 03    continue
 02   continue
      vardata = m/vardata
      chisqnu = chisq/nu    
      varfit  = chisqnu*vardata
c  
c  
      return
      end
c  
c.......................................................................
c..eof