/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:56 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "stgfnd.h"
#include <stdlib.h>
void /*FUNCTION*/ stgfnd(
float x[],
float y[],
long triang[],
long nt,
float q[],
long *indtri,
long tri[],
float s[][3],
long *mode)
{
	long int icount, iskip, j, j1, j2, j2save;
	float p1, p2;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Q = &q[0] - 1;
	long *const Tri = &tri[0] - 1;
	float *const X = &x[0] - 1;
	float *const Y = &y[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1997-06-15 STGFND CLL In args: Removed NP. Changed meaning of MODE.
	 *>> 1996-03-30 STGFND  Krogh  Removed Fortran 90 comments.
	 *>> 1996-02-02 STGFND CLL
	 *>> 1995-09-26 STGFND CLL Editing for inclusion into MATH77.
	 *
	 *     C.L.LAWSON, JPL, 1976 DEC 3
	 *         LOOK UP POINT Q IN TRIANGULAR GRID.
	 *
	 *   X(1:NP), Y(1:NP) [in]  (x,y) coordinates of vertices of the
	 *                          triangular grid.
	 *
	 *         {NP is the number of vertices in the triangular grid, but is
	 *          not explicitly used in this subroutine.}
	 *
	 *   TRIANG(1:6*NT) [in]  Array of integer pointers defining the
	 *            connectivity of the triangular grid.
	 *
	 *   NT [in]  No. of triangles in the triangular grid.
	 *
	 *   Q(1:2) [in]  The (x,y) coordinates of the point for which this
	 *             subr will attempt to find an enclosing triangle.
	 *
	 *   INDTRI [in]  DESIGNATES TRIANGLE AT WHICH SEARCH WILL BEGIN.
	 *              ON RETURN INDTRI IS THE INDEX OF THE LAST TRIANGLE
	 *              TESTED.  THIS IS THE TRIANGLE CONTAINING Q IF MODE=0.
	 *              IF MODE = 1, 2, or 3, THIS IS A BOUNDARY TRIANGLE AND Q
	 *              IS OUTSIDE A BOUNDARY EDGE OF THIS TRIANGLE.
	 *
	 *   TRI(1:7) [out]  INTEGER ARRAY OF POINTERS ASSOCIATED
	 *              WITH TRIANGLE INDTRI.
	 *
	 *   S(1:3, 1:3) [out]  If MODE=0 on return, S(,) CONTAINS
	 *              S(:,1) = UNNORMALIZED BARYCENTRIC COORDINATES OF Q.
	 *              S(:,2) = X COORDINATES OF EDGE VECTORS.
	 *              S(:,3) = Y COORDINATES OF EDGE VECTORS.
	 *
	 *   MODE [out]
	 *      = 0  MEANS OK.  Q IS INTERIOR OR ALMOST SO.
	 *      = 1, 2, or 3.   Means Q is exterior by more than the built-in
	 *           tolerance. The value of MODE identifies the edge of the
	 *           triangle relative to which Q is outside.  TRI(MODE) will be
	 *           zero, indicating a boundary edge.  Q is outside this edge.
	 *           The vertices at the ends of this edge, in counterclockwise
	 *           order are TRI(3+MODE) and TRI(4+MODE).
	 *
	 *      = -1  BAD.  SUBR IS CYCLING IN THE SEARCH.  THIS SHOULD NEVER
	 *              HAPPEN.
	 *
	 *     ------------------------------------------------------------------
	 *         Details of the contents of TRI(1:7) and S(1:3, 1:3).
	 *
	 *  For descriptive convenience we regard the subscript of P() and the
	 *  1st subsubscript of S(,) as always being reduced modulo 3 to 1, 2,
	 *  or 3.  Also for convenience we shall write P(i) to mean the vertex
	 *  indexed by P(i).
	 *
	 *  TRI() contains integer pointers defining the triangle indexed by
	 *  INDTRI.  The indices of the vertices of this triangle, in counter-
	 *  clockwise order are P(1), P(2), and P(3) which may be obtained as
	 *  P(i) = TRI(3+i) for i = 1, 2, and 3.  TRI(7) contains the same
	 *  value as TRI(4).  The triangle adjacent to this triangle across
	 *  the edge from P(3+i) to P(4+i) is indexed by TRI(i) for i = 1, 2,
	 *  and 3.  If there is no adjacent triangle across this edge then
	 *  TRI(i) = 0.
	 *
	 *  The unnormalized barycentric coordinate that is zero along the edge
	 *  from P(i) to P(i+1) and has a positive value at P(i+2) is stored in
	 *  S(i,1).  The (x,y) coordinates of the vector from P(i) to P(i+1) are
	 *  stored in (S(i,2), S(i,3)).
	 *     ------------------------------------------------------------------
	 *--S replaces "?": ?TGFND, ?TGGET
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------ */
	if (*indtri <= 0 || *indtri > nt)
		*indtri = 1;
	*mode = -1;
	iskip = 0;
	j2save = 0;
	for (icount = 1; icount <= nt; icount++)
	{
		stgget( *indtri, tri, triang );
		Tri[7] = Tri[4];
 
		for (j = 1; j <= 3; j++)
		{
			j1 = Tri[j + 3];
			if (j1 == j2save)
			{
				iskip = j;
			}
			else
			{
				j2 = Tri[j + 4];
				s[1][j - 1] = X[j2] - X[j1];
				p2 = s[1][j - 1]*(Q[2] - Y[j1]);
				s[2][j - 1] = Y[j2] - Y[j1];
				p1 = s[2][j - 1]*(Q[1] - X[j1]);
				s[0][j - 1] = -p1 + p2;
 
				if (s[0][j - 1] < -1.0e-4*(fabsf( p1 ) + fabsf( p2 )))
				{
 
					/*                        Q IS OUTSIDE TRIANGLE INDTRI BY MORE THAN
					 *                        TOLERANCE.  MOVE TO NEIGHBORING TRIANGLE IF
					 *                        THERE IS ONE.
					 * */
					if (Tri[j] != 0)
					{
						j2save = j2;
						*indtri = Tri[j];
						goto L_40;
					}
					else
					{
						*mode = j;
						goto L_80;
					}
				}
			}
		}
		/*                                  Q IS INSIDE OR ALMOST INSIDE
		 *                                  THIS TRIANGLE */
		*mode = 0;
		goto L_90;
L_40:
		;
	}
L_80:
	;
L_90:
	;
	/*                             COMPUTE SKIPPED INFO FOR USE IN
	 *                             INTERPOLATION SUBR. */
	if (iskip != 0)
	{
		j1 = Tri[iskip + 3];
		j2 = Tri[iskip + 4];
		s[1][iskip - 1] = X[j2] - X[j1];
		s[2][iskip - 1] = Y[j2] - Y[j1];
		s[0][iskip - 1] = -s[2][iskip - 1]*(Q[1] - X[j1]) + s[1][iskip - 1]*
		 (Q[2] - Y[j1]);
	}
	return;
} /* end of function */