/* @(#)pixlin.c	17.1.1.1 (ES0-DMD) 01/25/02 17:40:35 */
/*===========================================================================
  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  Cpixlin
.LANGUAGE    C
.AUTHOR      K. Banse		ESO - Garching
.KEYWORDS    bulk data frame, pixel values
.PURPOSE     calculate real pixel no's of a line across a 2-dim frame
.ALGORITHM   analytical geometry
.INPUT/OUTPUT
  call as    nlin = Cpixlin(xa,ya,xb,yb,step,ndim,xindx,yindx)

  input:
             float xa		x-start (in pixel no.)
             float ya		y-start (in pixel no.)
             float xb		x-end (in pixel no.)
             float yb		y-end (in pixel no.)
             double step	step size along the line
             int   ndim		dimension of XINDX (= dimension YINDX)
  output:
             float *xindx	buffer with real pixel no's. in x-direction
             float *yindx	buffer with real pixel no's. in y-direction

.RETURNS     actual number of points on line

.ENVIRONment MIDAS

.VERSIONS    1.00       941021  converted from ZIMA.FOR    
------------------------------------------------------------*/

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

#define  _POSIX_SOURCE 1


#include <midas_def.h>			/* for proto types */
#include <math.h>

 
/*

*/

#ifdef __STDC__
int Cpixlin(float xa,float ya,float xb,float yb,
            double step,int ndim,float *xindx,float *yindx)
#else
int Cpixlin(xa,ya,xb,yb,step,ndim,xindx,yindx)
int    ndim;
float  xa, ya, xb, yb, *xindx, *yindx;
double step;
#endif

{
int    nindx = 0;

register double pix, rm, turbo;


/* in order to make step not too small we move along the x-axis if the angle
   of the line is less than 45 degrees, else we move along the y-axis */

pix = xb - xa;
if (fabs(pix) < 10.e-10)			/* is it parallel to y-dir? */
   {
   rm = 0.0;
   goto move_y;				/* Yes, we move in y-direction */
   }

rm = (yb-ya) / pix;
if (fabs(rm) > 1.0001) 
   {
   rm = 1.0 / rm;
   step *= cos( atan( rm ));
   goto move_y;				/* we move in y-direction */
   }


/* here we move along in x-direction */

step *= cos( atan( rm ));
turbo = ya - (rm * xa);
pix = xa;

if (xa > xb)
   {
   while ( (nindx < ndim) && (pix >= xb) )
      {
      *xindx++ = (float) pix;
      *yindx++ = (float) ((rm * pix) + turbo);
      pix -= step;
      nindx ++;
      }
   }

else
   {
   while ( (nindx < ndim) && (pix <= xb) )
      {
      *xindx++ = (float) pix;
      *yindx++ = (float) ((rm * pix) + turbo);
      pix += step;
      nindx ++;
      }
   }

return (nindx);



/* here we move along in y-direction */

move_y:
turbo = xa - (rm * ya);
pix = ya;

if (ya > yb)
   {
   while ( (nindx < ndim) && (pix >= yb) )
      {
      *yindx++ = (float) pix;
      *xindx++ = (float) ((rm * pix) + turbo);
      pix -= step;
      nindx ++;
      }
   }

else
   {
   while ( (nindx < ndim) && (pix <= yb) )
      {
      *yindx++ = (float) pix;
      *xindx++ = (float) ((rm * pix) + turbo);
      pix += step;
      nindx ++;
      }
   }

return (nindx);
}