/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:17 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_sivx s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_sivx.h"
#include <string.h>
		/* PARAMETER translations */
#define	IFDIM	(16*INEQ + 1)
#define	INEQ	2
#define	IYDIM	(4*INEQ)
#define	MEDDIG	12
#define	MERET	51
#define	NDIG	4
#define	TOL	(powif(10.e0,-NDIG))
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int id[5], _i, _r;
	static long int iopt[12], kord[6];
	float rd[3];
	static float f[IFDIM], y[IYDIM];
	static long mact[3]={MEDDIG,7,MERET};
	static char mdummy[1][3]={"  "};
	static long neq = 2;
	static int _aini = 1;
	/* EQUIVALENCE translations */
	static float  _es0[4];
	float *const delt = (float*)((float*)_es0 + 2);
	float *const h = (float*)((float*)_es0 + 1);
	float *const t = (float*)_es0;
	float *const tfinal = (float*)((float*)_es0 + 3);
	float *const tspecs = (float*)_es0;
	/* end of EQUIVALENCE translations */
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const F = &f[0] - 1;
	long *const Id = &id[0] - 1;
	long *const Iopt = &iopt[0] - 1;
	long *const Kord = &kord[0] - 1;
	long *const Mact = &mact[0] - 1;
	float *const Rd = &rd[0] - 1;
	float *const Tspecs = &tspecs[0] - 1;
	float *const Y = &y[0] - 1;
		/* end of OFFSET VECTORS */
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		*t = 0.e0;
		*h = 1.e0;
		*delt = 100.e0;
		*tfinal = 4.e0;
		Y[1] = 1.e0;
		Y[2] = 0.e0;
		Y[3] = 0.e0;
		Y[4] = 1.e0;
		Iopt[1] = 16;
		Iopt[2] = 6;
		Iopt[3] = 3;
		Kord[6] = 2;
		F[3] = TOL;
		Iopt[4] = 17;
		Iopt[5] = 2;
		Iopt[6] = 6;
		Iopt[7] = 1;
		Iopt[8] = 10;
		Iopt[9] = 5;
		Iopt[10] = 0;
		Iopt[11] = 11;
		Iopt[12] = 0;
		_aini = 0;
	}
 
	/*>> 2009-11-03 DRSIVX  Krogh -- Used option 11
	 *>> 2006-04-10 DRSIVX  Krogh -- * dim for Y and F in sivao.
	 *>> 2001-07-16 DRSIVX  Krogh -- Declared type for parameter TOL.
	 *>> 2001-05-25 DRSIVX  Krogh Minor change for making .f90 version.
	 *>> 1996-06-14 DRSIVX  Krogh   Small change in output format
	 *>> 1994-09-12 DRSIVX  Krogh   Fixed for CHGTYP.
	 *>> 1994-08-19 DRSIVX  Krogh   Specified no. of digits in the output.
	 *>> 1993-05-05 DRSIVX  Krogh   Adjusted to simplify conversion to C.
	 *>> 1992-03-10 DRSIVX Fred T. Krogh
	 *
	 *--S replaces "?": DR?IVX, ?IVA, ?IVACO, ?IVADB, ?IVAF, ?IVAG, ?IVAO,
	 *--&              ?MESS
	 *
	 * Sample driver for SIVA --  Set up to solve two second order equations.
	 * Test SIVAG, SIVADB, and SIVACO
	 * */
	/*++S Default NDIG = 4
	 *++  Default NDIG = 10
	 *++ Substitute for NDIG below */
 
	/* Parameters for setting up message processor to specify no. of digits. */
 
 
	/* Set option for error control, local absolute error < 1.E-10. */
	/* Group the system to be treated as a single unit, set tolerance value
	 * Set option for second order equations */
	/* Set option for a G-Stop */
	/* Get diagnostic output for 5 steps to test MESS */
	/* Set option to initialize some space to 0, and end of options. */
 
 
	/* Specify number of digits in the output. */
	smess( mact, (char*)mdummy,3, mact, f );
 
	/* Do the integration
	 * */
	Kord[1] = 0;
L_100:
	;
	siva( tspecs, y, f, kord, neq, sivaf, sivao, 4, IYDIM, IFDIM,
	 6, iopt );
	if (Kord[1] != 1)
		goto L_100;
	sivadb( 45, tspecs, y, f, kord, "0Sample of SIVA Debug" );
	sivaco( id, rd );
	printf(" KEMAX=%1ld  KSTEP=%4ld  NUMDT=%2ld  EMAX=%14.7e\n", Id[1], Id[2], Id[3],
	   Rd[1]);
	exit(0);
} /* end of function */
 
void /*FUNCTION*/ sivaf(
float t[],
float y[],
float f[],
long kord[])
{
	float tp;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const F = &f[0] - 1;
	long *const Kord = &kord[0] - 1;
	float *const T = &t[0] - 1;
	float *const Y = &y[0] - 1;
		/* end of OFFSET VECTORS */
 
 
	/* Sample derivative subroutine for use with SIVA
	 * This evaluates derivatives for a simple two body problem.
	 * */
 
	/* Evaluate the derivatives
	 * */
	tp = Y[1]*Y[1] + Y[3]*Y[3];
	tp = 1.e0/(tp*sqrtf( tp ));
	F[1] = -Y[1]*tp;
	F[2] = -Y[3]*tp;
	return;
} /* end of function */
 
void /*FUNCTION*/ sivao(
float tspecs[],
float y[],
float f[],
long kord[])
{
	long int iflag, nstop;
	static float g6[1], gt6[1];
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const F = &f[0] - 1;
	float *const G6 = &g6[0] - 1;
	float *const Gt6 = &gt6[0] - 1;
	long *const Kord = &kord[0] - 1;
	float *const Tspecs = &tspecs[0] - 1;
	float *const Y = &y[0] - 1;
		/* end of OFFSET VECTORS */
 
 
	/* Sample output subroutine for use with SIVA.
	 * This subroutine gives output for a simple two body problem.
	 * */
 
 
	/* Do the output
	 * */
	iflag = 0;
	if (Kord[1] == 1)
	{
		printf("            RESULTS FOR A SIMPLE 2-BODY PROBLEM\n\n     T/IFLAG          U/V           UP/VP         UPP/VPP\n");
	}
	if (Kord[1] == 6)
	{
L_100:
		G6[1] = Y[3];
		sivag( tspecs, y, f, kord, &iflag, &nstop, g6, gt6 );
		if ((iflag == 1) || (iflag == 3))
			return;
		if (iflag == 4)
			goto L_100;
	}
	printf("%15.6e%15.6e%15.6e%15.6e\n        %3ld    %15.6e%15.6e%15.6e\n \n", Tspecs[1],
	   Y[1], Y[2], F[1], iflag, Y[3], Y[4], F[2]);
	return;
} /* end of function */