/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:46 */
/*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 "sivag.h"
#include <stdlib.h>
#include <string.h>
		/* PARAMETER translations */
#define	KDIM	16
#define	MAXORD	2
#define	MAXSTF	1
#define	MEEMES	52
#define	MEMDA1	27
#define	MERET	51
		/* end of PARAMETER translations */
 
	/* COMMON translations */
struct t_sivaev {
	float eeps2, eept75, eovep2, ovtm75, ovd10, eeps10, eeps16, erov10;
	}	sivaev;
struct t_sivasc {
	float tn, xi[KDIM];
	long int iopst, kordi, kqmaxd, kqmaxi, ldt, maxdif, maxint, nkdko,
	 nte, nyny, ndtf, numdt;
	}	sivasc;
struct t_sivamc {
	float tg[2], tgstop[2], tmark, tmarkx, tout, tolg, hc, hdec, hinc,
	 hincc, hmax, hmaxp9, hmin, alpha[KDIM], beta[KDIM + 1], d[MAXORD][MAXSTF + MAXORD],
	 g[MAXORD][KDIM], v[KDIM + MAXORD], ds[MAXORD][MAXSTF + MAXORD],
	 gs[KDIM], sigma[KDIM], rbq[KDIM], dnoise, eave, eimax, eimin,
	 emax, erep, robnd, snoise, fdat[11];
	long int icf, ics, igflg, igtype[2], igstop[2], ilgrep, ings,
	 iop3, iop4, iop5, iop6, iop7, iop8, iop9, iop10, iop11, iop12,
	 iop13, iop14, iop15, iop16, iop17, iop18, iop19, iop20, iop21,
	 iop22, iop21s, itolep, iy, kemax, kis, kmark, kord1i, kord2i,
	 kpred, kqdcon, kqicon, kqmaxs, kqmxds, kqmxil, kqmxip, kqmxis,
	 ksc, ksout, ksstrt, kstep, lex, linc, lincd, lincq, lsc, maxkqd,
	 maxkqi, method, ne, neptol, ng, ngtot, noiseq, noutko, ntolf,
	 ny, idat[6];
	}	sivamc;
	/* end of COMMON translations */
void /*FUNCTION*/ sivag(
float tspecs[],
float y[],
float f[],
long kord[],
long *iflag,
long *nstop,
float gnew[],
float gt[])
{
	long int i, ig;
	static float gold, told, tsave;
	static char mtxtaa[1][96]={"SIVAG$BCall with bad values of KORD.  KORD(1)=$I, KORD(2)=$I, when TSPECS(1)=$F and KSTEP=$M.$E"};
	static long mact[7]={MEMDA1,0,MEEMES,68,24,0,MERET};
	float *const hh = (float*)sivamc.g;
	long int *const igflgs = (long*)&sivamc.itolep;
	long int *const izflag = (long*)&sivamc.iy;
	long int *const kexit = (long*)&sivamc.iop17;
	long int *const ngstop = (long*)&sivamc.iop6;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const F = &f[0] - 1;
	float *const Gnew = &gnew[0] - 1;
	float *const Gt = &gt[0] - 1;
	long *const Igstop = &sivamc.igstop[0] - 1;
	long *const Igtype = &sivamc.igtype[0] - 1;
	long *const Kord = &kord[0] - 1;
	long *const Mact = &mact[0] - 1;
	long *const Ngstop = &ngstop[0] - 1;
	float *const Tg = &sivamc.tg[0] - 1;
	float *const Tgstop = &sivamc.tgstop[0] - 1;
	float *const Tspecs = &tspecs[0] - 1;
	float *const Xi = &sivasc.xi[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.
	 *>> 2001-09-07 SIVAG  Krogh  Changes to allow user tol on G-Stops.
	 *>> 1995-06-20 SIVAG  Krogh  Fixed problem introduced with last change.
	 *>> 1995-05-09 SIVAG  Krogh  Fixed G-Stop/discontinuity code interaction
	 *>> 1994-11-11 SIVAG  Krogh  Declared all vars.
	 *>> 1994-10-20 SIVAG  Krogh  Changes to use M77CON
	 *>> 1994-09-12 SIVAG  Krogh  Added CHGTYP code.
	 *>> 1994-08-17 SIVAG  Krogh  Modified internal comment.
	 *>> 1994-03-07 SIVAG  Krogh  Allow larger order in single precision.
	 *>> 1993-04-27 SIVAG  Krogh  Additions for Conversion to C.
	 *>> 1992-10-12 SIVAG  Krogh  Fixed G-Stop/discontinuity code interaction
	 *>> 1992-09-17 SIVAG  Krogh  Slight change in check for sign change.
	 *>> 1992-04-08 SIVAG  Krogh  Unused labels 140,150,230, and 250 removed.
	 *>> 1992-03-10 SIVAG  Krogh  Fixed value for KDIM in single p. version.
	 *>> 1990-01-29 SIVAG  Krogh  Added arg to call to DERMN.
	 *>> 1988-03-04 SIVAG  Krogh  Initial code.
	 *
	 *--S replaces "?": ?IVAG,?IVABU,?IVAEV,?IVAIN,?IVAMC,?IVASC,?MESS,?ZERO
	 *
	 *     -XXXG(TSPECS, Y, F, KORD, IFLAG, NSTOP, GNEW, GT)
	 *
	 *  SUBROUTINE TO LOCATE OUTPUT POINTS AT ZEROS OF ARBITRARY
	 *  FUNCTIONS  **** GSTOPS **** FOR DIFFERENTIAL EQUATION
	 *  INTEGRATOR -ODE (OR -IVA).
	 * */
	/*--S Next line special: P=>D, X=>Q */
 
	/*++SP Default KDIM = 16
	 *++  Default KDIM = 20
	 *++  Default MAXORD = 2, MAXSTF = 1 */
	/*++ Substitute for KDIM, MAXORD, MAXSTF below */
	/*--S Next line special: P=>D, X=>Q */
 
	/*--S Next line special: P=>D, X=>Q */
 
 
	/*.    SPECIFICATION OF ENVIRONMENTAL CONSTANTS. */
 
 
 
	/*                      Declarations for error message processing.
	 * */
	/* ********* Error message text ***************
	 *[Last 2 letters of Param. name]  [Text generating message.]
	 *AA SIVAG$B
	 *AB Call with bad values of KORD.  KORD(1)=$I, KORD(2)=$I, when $C
	 *   TSPECS(1)=$F and KSTEP=$M.$E */
 
	/* ****************** START OF EXECUTABLE CODE **************
	 * */
	*iflag = 1;
	ig = Kord[2];
	if ((ig != 0) && (ig != 1))
		goto L_500;
	switch (IARITHIF(sivamc.igflg - 3))
	{
		case -1: goto L_10;
		case  0: goto L_80;
		case  1: goto L_70;
	}
L_10:
	switch (IARITHIF(sivamc.igflg - 1))
	{
		case -1: goto L_20;
		case  0: goto L_210;
		case  1: goto L_60;
	}
 
	/* ******************** INITIAL POINT ***********************
	 * */
L_20:
	if (sivamc.igflg == -3)
	{
		sivamc.igflg = 5;
		return;
	}
	Igtype[ig + 1] = 0;
	if ((Ngstop[ig + 1] <= 0) || (ig + sivamc.igflg == -1))
		goto L_30;
	sivamc.igflg = ig - 2;
	goto L_40;
L_30:
	sivamc.igflg = 5;
L_40:
	sivamc.ng = Ngstop[2 - ig];
	for (i = 1; i <= sivamc.ng; i++)
	{
		Gt[i] = Gnew[i];
	}
	goto L_480;
 
	/*     **** USER HAS BEEN TOLD THAT A GSTOP WAS FOUND
	 *          TEST IF CALLED FROM OUTPUT WHEN CALL SHOULD
	 *          BE FROM DERIVS */
L_60:
	if ((Igtype[1] == 0) && (ig != 0))
		goto L_420;
	/*     **** PROTECT AGAINST NOISEY G NEAR THE ZERO */
	if (Gnew[sivamc.ings] == 0.e0)
		Gnew[sivamc.ings] = Gt[sivamc.ings];
 
	/* ****** TEST FOR CHANGE IN THE SIGN OF A G ****************
	 * */
L_70:
	sivamc.ng = Ngstop[2 - ig];
	sivamc.ings = 0;
L_80:
	sivamc.ings += 1;
	if (sivamc.ings > sivamc.ng)
		switch (IARITHIF(sivamc.igflg - 4))
		{
			case -1: goto L_400;
			case  0: goto L_380;
			case  1: goto L_480;
		}
	switch (SARITHIF(Gnew[sivamc.ings]))
	{
		case -1: goto L_90;
		case  0: goto L_100;
		case  1: goto L_110;
	}
L_90:
	switch (SARITHIF(Gt[sivamc.ings]))
	{
		case -1: goto L_120;
		case  0: goto L_120;
		case  1: goto L_130;
	}
L_100:
	switch (SARITHIF(Gt[sivamc.ings]))
	{
		case -1: goto L_130;
		case  0: goto L_80;
		case  1: goto L_130;
	}
L_110:
	switch (SARITHIF(Gt[sivamc.ings]))
	{
		case -1: goto L_130;
		case  0: goto L_120;
		case  1: goto L_120;
	}
L_120:
	Gt[sivamc.ings] = Gnew[sivamc.ings];
	goto L_80;
 
	/* ********* A SIGN CHANGE HAS BEEN FOUND *******************
	 * */
L_130:
	*nstop = sivamc.ings;
	if (ig == 0)
		*nstop = -sivamc.ings;
	if (sivamc.igflg != 5)
		goto L_200;
	/*     **** USUAL CASE -- TEST IF OUTPUT POINT PRECEDES THE
	 *          SIGN CHANGE, OR IF PREDICTING, CORRECTING, OR
	 *          NOT THROUGH THE FIRST STEP.
	 *     **** TEST IF AN INTERPOLATED G WAS WHAT CHANGED SIGN */
	if (ig != 0)
		goto L_180;
	/*     **** BACK UP DIFFERENCES AND OTHER STUFF TO BEGINNING
	 *          OF THE STEP */
	sivabu( f, kord );
	/*     **** TEST IF CORRECTING */
	if (sivamc.kord1i == 2)
		goto L_170;
	/*     **** TEST IF THROUGH THE FIRST STEP */
	if (sivamc.lsc < 4)
		goto L_180;
	/*     **** IF FIRST DERIVATIVE EVALUATION OF THE FIRST
	 *          STEP, FIND THE GSTOP, AND USE IT TO GET A NEW
	 *          INITIAL STEPSIZE */
	if (sivamc.lsc == 7)
		goto L_200;
	/*     **** SET NEW STEPSIZE AFTER SIGN CHANGE WHEN STARTING */
L_160:
	*hh = Tspecs[1] - sivasc.tn;
	/*     **** SET KEXIT TO TRY NEW STEPSIZE */
L_170:
	*kexit = 1;
	Tspecs[1] = sivasc.tn;
	goto L_460;
	/*     **** SET KEXIT FOR USUAL CASE */
L_180:
	*kexit = ig + 2;
	/*     **** TEST IF SIGN CHANGE IN G PRECEDES NEXT OUTPUT PT. */
	switch (SARITHIF(*hh*(Tspecs[1] - sivamc.tmark)))
	{
		case -1: goto L_200;
		case  0: goto L_200;
		case  1: goto L_190;
	}
	/*     **** SET UP TO EVALUATE G AT OUTPUT POINT */
L_190:
	sivamc.igflg = 4;
	Tspecs[1] = sivamc.tmark;
	*nstop = 0;
	goto L_240;
 
	/* ***************** FIND THE ZERO OF G *********************
	 *
	 *     **** INITIALIZE ZERO FINDER */
L_200:
	told = Tg[2 - ig];
	gold = Gt[sivamc.ings];
	tsave = Tspecs[1];
	*igflgs = sivamc.igflg;
	sivamc.igflg = 1;
	*izflag = 0;
	goto L_220;
	/*     **** TEST IF ZERO ALREADY FOUND */
L_210:
	switch (IARITHIF(*izflag - 1))
	{
		case -1: goto L_350;
		case  0: goto L_220;
		case  1: goto L_310;
	}
L_220:
	;
	szero( &Tspecs[1], &Gnew[sivamc.ings], &told, &gold, izflag, sivamc.tolg );
	/*     **** TEST FOR CONVERGENCE */
	if (*izflag != 1)
		goto L_260;
	/*     **** INTERPOLATE NEW Y, AND GO COMPUTE G AGAIN */
L_240:
	sivain( &Tspecs[1], y, f, kord );
	*iflag = 4;
	sivamc.kord2i = ig - 3;
	return;
	/*     **** CONVERGENCE -- CHOOSE TOLD TO GIVE A CHANGE
	 *          IN SIGN */
L_260:
	if (Gnew[sivamc.ings] == 0.e0)
		goto L_290;
	switch (SARITHIF(Tspecs[1] - told))
	{
		case -1: goto L_270;
		case  0: goto L_300;
		case  1: goto L_280;
	}
L_270:
	switch (SARITHIF(*hh))
	{
		case -1: goto L_290;
		case  0: goto L_300;
		case  1: goto L_300;
	}
L_280:
	switch (SARITHIF(*hh))
	{
		case -1: goto L_300;
		case  0: goto L_300;
		case  1: goto L_290;
	}
L_290:
	told = Tspecs[1];
	/*     **** CHECK IF SIGN CHANGE DUE TO NOISE */
L_300:
	Tspecs[1] = told + Xi[1];
	goto L_240;
L_310:
	Tspecs[1] = told;
	switch (SARITHIF(Gnew[sivamc.ings]))
	{
		case -1: goto L_320;
		case  0: goto L_340;
		case  1: goto L_330;
	}
L_320:
	switch (SARITHIF(Gt[sivamc.ings]))
	{
		case -1: goto L_340;
		case  0: goto L_360;
		case  1: goto L_360;
	}
L_330:
	switch (SARITHIF(Gt[sivamc.ings]))
	{
		case -1: goto L_360;
		case  0: goto L_360;
		case  1: goto L_340;
	}
	/*     **** ZERO WAS EVIDENTLY DUE TO NOISE */
L_340:
	Tspecs[1] = tsave;
	*izflag = 0;
	goto L_370;
L_350:
	sivamc.igflg = *igflgs;
	/*     SET KORD2I TO INITIAL VALUE TO AVOID LOOP */
	sivamc.kord2i = ig;
	goto L_80;
	/*     **** SAVE INFORMATION ABOUT THE STOP */
L_360:
	sivamc.igflg = 3;
	Tgstop[2 - ig] = Tspecs[1];
	Igtype[2 - ig] = *izflag + 3;
	Igstop[2 - ig] = *nstop;
L_370:
	*nstop = 0;
	goto L_240;
 
	/* ************** AFTER SEARCH FOR A SIGN CHANGE ************
	 *
	 *     **** NO SIGN CHANGE AT A T OUTPUT POINT
	 *     TEST IF CALLED FROM OUTPUT */
L_380:
	if (ig != 0)
		goto L_390;
	/*     SET UP FOR CALL TO OUTPUT */
	sivamc.kord1i = 7;
	sivamc.kord2i = -2;
	*iflag = 3;
	goto L_480;
	/*     **** ADJUST KEXIT AND SET UP TO GIVE OUTPUT */
L_390:
	*kexit += 2;
	sivamc.kord1i = min( sivamc.kmark, 5 );
	Kord[3] = sivamc.kmark;
	Kord[1] = sivamc.kord1i;
	sivamc.igflg = 5;
	*iflag = 2;
	goto L_470;
	/*     **** TEST IF USER HAS BEEN TOLD OF GSTOP */
L_400:
	if (sivamc.igflg == 2)
		goto L_450;
	/*     **** A GSTOP HAS BEEN FOUND
	 *     TEST IF STARTING */
	if (sivamc.lsc == 7)
		goto L_160;
	*iflag = Igtype[2 - ig];
	*nstop = Igstop[2 - ig];
	sivamc.ings = labs( *nstop );
	if (sivamc.ings == 0)
		goto L_410;
	Gt[sivamc.ings] = -Gt[sivamc.ings];
	if (ig == 0)
		goto L_430;
	sivamc.igflg = 2;
	/*     If interpolated GSTOP was found set to check again in case of
	 *     multiple stops at exactly the same point. */
	if (Igtype[1] != 0)
		goto L_440;
	/*     **** TELL USER OF AN EXTRAPOLATED GSTOP */
L_410:
	*iflag = Igtype[2];
	*nstop = Igstop[2];
	sivamc.ings = labs( *nstop );
	/*     **** SET SO DERIVS IS CALLED WITH KORD(1) = KPRED */
L_420:
	sivamc.kord1i = sivamc.kpred;
	sivamc.kord2i = -3;
	return;
	/*     **** AN EXTRAPOLATED GSTOP WAS FOUND, SET UP TO CHECK
	 *          INTERPOLATED STOPS (IF ANY) */
L_430:
	sivamc.ng = Ngstop[1];
	sivamc.ings = 0;
	*iflag = 3;
	*nstop = Igstop[2];
	/*     **** SET SO OUTPUT IS CALLED WITH KORD(1) = 7 */
L_440:
	sivamc.kord1i = 7;
	sivamc.kord2i = -2;
	goto L_490;
	/*     **** CHECK THAT AN EXTRAPOLATED G-STOP IS NOT MISSED */
L_450:
	if ((ig == 0) || (Igtype[2] == 0))
		goto L_460;
	/*     SET TO CHECK FOR INTERPOLATED G-S. */
	Tg[1] = Tspecs[1];
	Igtype[1] = 0;
	Tspecs[1] = Tgstop[2];
	sivamc.ings = 0;
	sivamc.igflg = 3;
	goto L_240;
	/*     **** SET SO INTEGRATOR GOES TO PLACE DIRECTED BY KEXIT */
L_460:
	*nstop = 0;
	sivamc.igflg = 5;
	*iflag = 3;
L_470:
	sivamc.kord2i = -7;
	/*     **** STORE INFO. ON LAST G COMPUTED */
L_480:
	Igtype[2 - ig] = 0;
L_490:
	Tg[2 - ig] = Tspecs[1];
	return;
 
	/* ********************** ERROR PROCESSING ******************
	 * */
L_500:
	Mact[2] = sivamc.kstep;
	/*--S Next line special: P=>S, X=>D */
	smess( mact, (char*)mtxtaa,96, kord, tspecs );
	*iflag = 8;
	*kexit = 6;
	goto L_470;
} /* end of function */