/*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_sivdb s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_sivdb.h"
/*      program DRSIVDB
 *>> 2009-11-03 DRSIVDB  Krogh -- Used option 11
 *>> 2001-07-16 DRSIVDB  Krogh -- Declared type for paramter TOL.
 *>> 2001-05-25 DRSIVDB  Krogh -- Added comma to format.
 *>> 1996-06-14 DRSIVDB  Krogh  Small change in output format
 *>> 1994-09-12 DRSIVDB  Krogh  Fixed for CHGTYP.
 *>> 1993-07-18 DRSIVDB  Krogh  Fixed to get same output in S.P. and D.P.
 *>> 1993-05-05 DRSIVDB  Krogh  Adjusted to simplify conversion to C.
 *>> 1989-05-04 Fred T. Krogh
 *
 *--S replaces "?": DR?IVDB, ?IVA, ?IVADB, ?IVAF, ?IVAO, ?MESS
 *
 * Sample driver for SIVA --  Set up to solve two second order equations.
 * Tests debug output
 * */
		/* PARAMETER translations */
#define	IFDIM	(16*INEQ + 1)
#define	INEQ	2
#define	IYDIM	(4*INEQ)
#define	NDIG	4
#define	TOL	(powif(10.e0,-NDIG))
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int _i, _r;
	static long int iopt[7], kord[6];
	static float f[IFDIM], y[IYDIM];
	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 Iopt = &iopt[0] - 1;
	long *const Kord = &kord[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 = 6.283185307179586477e0;
		*tfinal = 2.e1;
		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] = 11;
		Iopt[7] = 0;
		_aini = 0;
	}
 
	/*++ Code for SHORT_LINE is inactive
	 *      integer          MACT(5)
	 *++ END */
	/*++S Default NDIG = 4
	 *++  Default NDIG = 10
	 *++ Substitute for NDIG below */
 
	/* Parameters for setting up message processor to specify the number of
	 * digits, and the line length. */
 
 
	/* 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 to initialize some space to 0, and end of options. */
 
 
	/* 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;
	/*++ Code for SHORT_LINE is inactive
	 *      MACT(1) = MEDDIG
	 *      MACT(2) = 7
	 *      MACT(3) = MEMLIN
	 *      MACT(4) = 79
	 *      MACT(5) = MERET
	 *      call SMESS(MACT, ' ', MACT, F)
	 *++ END */
	sivadb( 44, tspecs, y, f, kord, " Test of SIVADB" );
	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;
	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)
{
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const F = &f[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
	 * */
	if (*kord == 1)
	{
		printf("            RESULTS FOR A SIMPLE 2-BODY PROBLEM\n\n        T             U/V           UP/VP         UPP/VPP\n");
	}
	printf("%15.6e%15.6e%15.6e%15.6e\n               %15.6e%15.6e%15.6e\n \n", Tspecs[1],
	   Y[1], Y[2], F[1], Y[3], Y[4], F[2]);
 
	return;
} /* end of function */