/* @(#)zima.c	17.1.1.1 (ES0-DMD) 01/25/02 17:40:37 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.COPYRIGHT:  Copyright (c) 1994 European Southern Observatory,
                                         all rights reserved
.IDENTIFIER  Czima
.LANGUAGE    C
.AUTHOR      K. Banse			IPG-ESO Garching
.KEYWORDS    bulk data frame, pixel values
.PURPOSE     get data values at given real pixel no's. of an image
.ALGORITHM   a list of real pixel no's. is scanned and corresponding data
             values of input array are stored in output array
             bilinear interpolation is used to obtain values for real pixels 
             not falling exactly in the grid of coordinates for pixel no's. 
             falling on the upward border lines in x or y, interpolation
             is done only in x or y, i.e. the border grid points are duplicated
             to the outside... 
             calculate also dynamic range on the fly
.INPUT/OUTPUT
  call as    Czima( p_in, npix, xpixl, ypixl, ndim, p_out, fmin, fmax )
  
  input:
             float *p_in	input array of data
             int   *npix	dimension of input data
             float *xpixl	buffer with real pixel no's. in x
             float *ypixl	buffer with real pixel no's. in y
             int   ndim         dimension of XPIXL + YPIXL
  
  output:
             float *p_out	output array of data
             float *fmin	minimum value of output array
             float *fmax	max value ...

.ENVIRONment MIDAS

.VERSIONS    1.00       940606  from ZIMA.FOR    R.M. van Hees
.VERSIONS    1.10       941122  fix it , KB
------------------------------------------------------------*/

/* Define _POSIX_SOURCE to indicate that this is a POSIX program */

#define  _POSIX_SOURCE 1

#include <midas_def.h>


/*

*/

#ifdef __STDC__
void Czima(float *p_in,int *npix,float *xpixl,float *ypixl,int ndim, 
           float *p_out,float *fmin,float *fmax)

#else
void Czima(p_in,npix,xpixl,ypixl,ndim,p_out,fmin,fmax)
int   *npix, ndim;
float *p_in, *p_out, *xpixl, *ypixl, *fmin, *fmax;
#endif

{
register int nn;
register float cindx;

int   indxx, indxy, indx1, indx2, indx3, indx4;
float val, xfrac, yfrac;



if (*npix <= 1)				/* only 1-D in Y */
   {
   nn = ndim/2;
   indx1 = ypixl[nn] - 1;
   if (indx1 < 0) 
      indx1 = 0;
   else if (indx1 > (npix[1]-1))
      indx1 = npix[1] - 1;
   *fmin = *fmax = p_in[indx1];

   for (nn=0; nn<ndim; nn++)
      {
      cindx = ypixl[nn] - 1.0;          /* move to [0,npix-1] */
      if (cindx <= 0.0)
         {
         indx1 = 0;
         yfrac = 0.0;
         }
      else if (cindx > (float) (npix[1]-1))
         {
         indx1 = npix[1] - 1;
         yfrac = 0.0;
         }
      else
         {
         indx1 = (int) cindx;
         yfrac = cindx - indx1;
         }
      indx2 = indx1 + 1;

      if (indx2 > (npix[1]-1))
         val = p_in[indx1];
      else
         val = p_in[indx1] + (yfrac * (p_in[indx2] - p_in[indx1]));

      if ( *fmin > val )
         *fmin = val;
      else if ( *fmax < val )
         *fmax = val;

      *p_out++ = val;
      }
   }

else if (npix[1] <= 1)			/* only 1-D in X */
   {
   nn = ndim/2;
   indx1 = xpixl[nn] - 1;
   if (indx1 < 0) 
      indx1 = 0;
   else if (indx1 > (npix[0]-1))
      indx1 = npix[0] - 1;
   *fmin = *fmax = p_in[indx1];

   for (nn=0; nn<ndim; nn++)
      {
      cindx = xpixl[nn] - 1.0;		/* move to [0,npix-1] */
      if (cindx <= 0.0)
         {
         indx1 = 0;
         xfrac = 0.0;
         }
      else if (cindx > (float) (npix[0]-1))
         {
         indx1 = *npix - 1;
         xfrac = 0.0;
         }
      else
         {
         indx1 = (int) cindx;
         xfrac = cindx - indx1;
         }
      indx2 = indx1 + 1;

      if (indx2 > (*npix-1))
         val = p_in[indx1];
      else
         val = p_in[indx1] + (xfrac * (p_in[indx2] - p_in[indx1]));

      if ( *fmin > val )
         *fmin = val;
      else if ( *fmax < val )
         *fmax = val;

      *p_out++ = val;
      }
   }

else						/* 2-D image */
   {
   int nsize;

   nsize = npix[0] * npix[1];

   nn = ndim/2;
   indx1 = xpixl[nn] - 1;
   if (indx1 < 0) 
      indx1 = 0;
   else if (indx1 > (npix[0]-1))
      indx1 = npix[0] - 1;
   indx2 = ypixl[nn] - 1;
   if (indx2 < 0) 
      indx2 = 0;
   else if (indx2 > (npix[1]-1))
      indx2 = npix[1] - 1;
   indx1 += (indx2*npix[0]);
   *fmin = *fmax = p_in[indx1];

   for (nn=0; nn<ndim; nn++)
      {
      cindx = xpixl[nn] - 1.0;          /* move to [0,npix-1] */
      if (cindx <= 0.0)
         {
         indxx = 0;
         xfrac = 0.0;
         }
      else if (cindx >= (float) (npix[0]-1))
         {
         indxx = npix[0] - 1;
         xfrac = 0.0;
         }
      else
         {
         indxx = (int) cindx;
         xfrac = cindx - indxx;
         }
      
      cindx = ypixl[nn] - 1.0;          /* move to [0,npix-1] */
      if (cindx <= 0.0)
         {
         indxy = 0;
         yfrac = 0.0;
         }
      else if (cindx >= (float) (npix[1]-1))
         {
         indxy = npix[1] - 1;
         yfrac = 0.0;
         }
      else
         {
         indxy = (int) cindx;
         yfrac = cindx - indxy;
         }

      indx1 = indxx + ((*npix) * indxy);
      indx2 = indx1 + 1;
      indx3 = indx1 + (*npix);
      if ((indxx + 1) >= *npix )             /* are we out of grid in x? */
         {
         if (indx2 >= nsize)			/* also out in y? */
            val = p_in[indx1];
         else
            val = p_in[indx1] + ( yfrac * (p_in[indx3] - p_in[indx1]) );
         }
      else
         {
         if (indx3 >= nsize)			/* are we out of grid in y? */
            val = p_in[indx1] + ( xfrac * (p_in[indx2] - p_in[indx1]) );
         else
            {
            indx4 = indx3 + 1;
            val = p_in[indx1] 
                  + xfrac * (p_in[indx2] - p_in[indx1])
                  + yfrac * (p_in[indx3] - p_in[indx1])
                  + xfrac * yfrac * (p_in[indx1] - p_in[indx2]
                  - p_in[indx3] + p_in[indx4]);
            }
         }

      if ( *fmin > val )
         *fmin = val;
      else if ( *fmax < val )
         *fmax = val;

      *p_out++ = val;
      }
   }
}