/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:07 */
/*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 "dtgc0.h"
#include <stdlib.h>
void /*FUNCTION*/ dtgc0(
double s[][3],
double *zout,
LOGICAL32 wantdz,
double dzout[])
{
	double delz[2], denom;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Delz = &delz[0] - 1;
	double *const Dzout = &dzout[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.
	 *>> 1996-02-02 DTGC0 CLL
	 *>> 1996-01-11 DTGC0 CLL
	 *>> 1995-11-01 DTGC0 CLL
	 *>> 1995-09-26 DTGC0 CLL Editing for inclusion into MATH77.
	 *     C.L.LAWSON, JPL, 1976 DEC 7
	 *  This subr interpolates over a triangle using linear interpolation.
	 *  This method gives C0 continuity with neighboring triangles, i.e.,
	 *  continuity of the value but generally not continuity of first and
	 *  higher order partial derivatives.
	 *  Optionally the subr also computes first partial derivatives.
	 *
	 *   (S(1:3, 1:4) [inout]  Columns 1-4 contain data set by the user to
	 *          specify the interpolation problem.  Col 1 contains
	 *          unnormalized barycentric coordinates of the interpolation
	 *          point.  The other cols contain data depending only on the
	 *          triangle and its vertex data and not on the interpolation
	 *          point.
	 *
	 *         S( ,1) = UNNORMALIZED BARYCENTRIC COORDS OF INTERP POINT.
	 *         S( ,2) = U = X COORD OF EDGE VECTOR
	 *         S( ,3) = V = Y COORD OF EDGE VECTOR
	 *         S( ,4) = Z = FCN VALUE AT VERTEX
	 *
	 *   ZOUT [out]  INTERPOLATED VALUE COMPUTED BY SUBR.
	 *
	 *   WANTDZ [in]  =.TRUE. MEANS COMPUTE DZOUT() AS WELL AS ZOUT.
	 *              =.FALSE. MEANS COMPUTE ONLY ZOUT AND NOT DZOUT().
	 *
	 *   DZOUT(1:2) [out]  First partial derivs w.r.t. X and Y of the
	 *      interpolated surface at the interpolation point.
	 *      Note that since this subroutine computes only a linear
	 *      interpolant over a triangle, the partial derivs w.r.t. X and Y
	 *      will be the same at every point of the triangle.  In general
	 *      these derivatives will jump to different values when moving to
	 *      an adjacent triangle.
	 *     ------------------------------------------------------------------
	 *         Details of the contents of TRI(1:7) and S(1:3, 1:4).
	 *
	 *  We assume the point Q is in 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(i) to P(i+1) is indexed by TRI(i) for i = 1, 2,
	 *  and 3.  If there is no adjacent triangle across this edge then
	 *  TRI(i) = 0.
	 *
	 *  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).
	 *
	 *  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)).
	 *  The function value Z at P(i+2) is stored in S(i,4) for i = 1, 2, & 3.
	 *     ------------------------------------------------------------------
	 *--D replaces "?": ?TGC0
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------
	 *        To compute zout:
	 *        Combine the vertex function values using the unnormalized
	 *        barycentric coordinates as weights, and normalize by the
	 *        sum of the unnormalized barycentric coordinates.
	 * */
	*zout = (s[0][0]*s[3][0] + s[0][1]*s[3][1] + s[0][2]*s[3][2])/
	 (s[0][0] + s[0][1] + s[0][2]);
	if (wantdz)
	{
		Delz[1] = s[3][2] - s[3][1];
		Delz[2] = s[3][0] - s[3][2];
		denom = s[1][0]*s[2][1] - s[1][1]*s[2][0];
		if (denom != 0.0e0)
		{
			Dzout[1] = (Delz[1]*s[2][1] - Delz[2]*s[2][0])/denom;
			Dzout[2] = (s[1][0]*Delz[2] - s[1][1]*Delz[1])/denom;
		}
		else
		{
			Dzout[1] = 0.0e0;
			Dzout[2] = 0.0e0;
		}
	}
	return;
} /* end of function */