c  this file contains subroutine gridim
c  
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c  
c  
c  
      subroutine         gridim(n,a)
c  
c  
c  
c.......................................................................
c  
c
c  author: Christopher W. Churchill
c
c  history: 20-Jan94   created, CWC
c           30-Jan94   optimization option added, CWC
c
c  purpose: (1) stuff the initial guess for the solution 
c           (2) stuff the 2D grid for fitting
c           (3) stuff the sigma weights 
c
c  input:   n        = number of fitting parameters
c           a        = array of fitting parameters
c           xmax     = order map [pixel locations] (APEXTRACT tracings)
c           gain     = instrumental gain [e-/DN]
c           rdnoise  = instrumental read-noise [e-]
c           pix      = input image pixel values [DN]
c
c  output:  a        = fitting parameters "initialization"
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
c  I/O: none
c
c  options: minflag  = (logical) flag for grid optimization
c           .true.   = "optimimzed" xmin, Imin, and Imax 
c                       calculates weighted mean interorder map  
c                       calculates Imax and Imin at fractional pix
c           .false.  =  xmin are xmax mid-points ; calculates Imax 
c                       and Imin to be extreme pixel value in radius 
c                       +/- 1pix from integer xmax and xmin
c
c  routines called: none
c
c  description: see subroutine comments
c
c.......................................................................
c  
      implicit undefined (a-z)
      integer            n,i,j,icol,irow,x1,x2
      double precision   a(1),fpix,xlo,xhi,y1,y2,fi,fsum,wi,wsum,dx,phi
      include            'hamscatt.par'
      include            'hamscatt.com'
c  
c  
c  
c  
c  
c  initialize the fitting coefficients to zero, except the local term
c  which is set to unity
      do 10 j=1,n-1
        a(j) = 0.0d0
 10   continue
      if (termflag) then
       a(n) = 1.0d0
      else
       a(n) = 0.0d0
      end if
c  
c
c  
      if (minflag) then
c
c  --- begin optimization block ---
c
c  account for fractional pixels for Imax(j,icol) and Imin(j,icol)
c  perform optimal determination of the values of xmin(j,icol), 
c  the weights are the difference 1/{xhi-xlo}, which are designed to 
c  measure the uncertainty in the true minima ; this should account 
c  for non-symmetric profiles and (with a little doctoring) ghost 
c  apertures
c  search across inter-order region for the minimum y1 at xlo and the 
c  second minimum y2 at xhi of column icol of aperture j
c  compute the weighted sum and the sum of weights for aperture j
      do 01 j=1,ndata
       fsum = 0.0d0
       wsum = 0.0d0
       do 03 icol=1,ncols
        x1  = nint(xmax(j,icol))  
        x2  = nint(xmax(j+1,icol))
        y1  = pix(icol,x1)
        xlo = float(x1)
        do 05 i=x1+1,x2-1
         if (pix(icol,i).le.y1) then
          y2  = y1
          xhi = xlo
          y1  = pix(icol,i)
          xlo = float(i)
         end if            
 05     continue
        wi   = 1.0d0/(abs(xlo-xhi))
        fi   = (xlo-xmax(j,icol))/(xmax(j+1,icol)-xmax(j,icol))
        wsum = wsum + wi
        fsum = fsum + wi*fi
 03    continue
c  
c  stuff xmin(j,icol), Imin(j,icol) Imax(j,icol) and sig(j,icol)
        do 02 icol=1,ncols
          xmin(j,icol) = xmax(j,icol) + 
     @                   fsum*(xmax(j+1,icol)-xmax(j,icol))/wsum
          irow         = int(xmin(j,icol))
          fpix         = xmin(j,icol) - float(irow)
          Imin(j,icol) = gain * ((1.0d0-fpix)*pix(icol,irow) + 
     @                   fpix*pix(icol,irow+1))
          irow         = int(xmax(j,icol))
          fpix         = xmax(j,icol) - float(irow)
          Imax(j,icol) = gain * ((1.0d0-fpix)*pix(icol,irow) + 
     @                    fpix*pix(icol,irow+1))
          dx           = abs(xmax(j,icol)-xmax(j+1,icol))
          phi          = (0.5*dx)**gamma
          sig(j,icol)  = sqrt(abs(Imin(j,icol)) + rdnoise**2)
 02     continue
 01   continue
c  
c  hand stuff the last aperture and bail
      do 06 icol=1,ncols
        irow               = int(xmax(norders,icol))
        fpix               = xmax(norders,icol) - float(irow)
        Imax(norders,icol) = gain * ((1.0d0-fpix)*pix(icol,irow) + 
     @                       fpix*pix(icol,irow+1))
 06   continue
c  
      end if
c
c  --- end optimization block ---
c
c
c  
c  --- begin non optimization block ---
c  
      if (.not.minflag) then
c
c  do not account for fractional pixels for the modeling data, assume 
c  1) Imax(j,icol) lies within +/- 1 pixel of xmax(j,icol)
c  2) xmin(j,icol) is the half-way point between apertures
c  3) Imin(j,icol) lies within +/- 1 pixel of xmin(j,icol)
c  note use of intrinsic function int and not nint
      do 11 j=1,ndata
        do 13 icol=1,ncols
          irow         = int(xmax(j,icol))
          Imax(j,icol) = gain * max(pix(icol,irow),pix(icol,irow+1))
          xmin(j,icol) = 0.5*(xmax(j,icol)+xmax(j+1,icol))
          irow         = int(xmin(j,icol))
          Imin(j,icol) = gain * min(pix(icol,irow),pix(icol,irow+1))
          dx           = abs(xmax(j,icol)-xmax(j+1,icol))
          phi          = (0.5*dx)**gamma
          sig(j,icol)  = sqrt(abs(Imin(j,icol)) + rdnoise**2)
 13     continue
 11   continue
c  
c  hand stuff the last aperture and bail
      do 15 icol=1,ncols
        irow               = int(xmax(norders,icol))
        Imax(norders,icol) = gain * max(pix(icol,irow),pix(icol,irow+1))
 15   continue
c  
      end if
c  
c  --- end non optimization block ---
c  
c
      return
      end
c  
c  
c  
c.......................................................................
c..eof