/*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 = >[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 */