/*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_siva s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_siva.h"
		/* PARAMETER translations */
#define	IFDIM	(16*INEQ + 1)
#define	IKDIM	6
#define	INEQ	2
#define	ITDIM	4
#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[6], kord[IKDIM];
	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] = 0;
		_aini = 0;
	}
 
	/*>> 2010-06-09 DRSIVA  Krogh Used parameters for all dimenssions.
	 *>> 2001-05-25 DRSIVA  Krogh Minor change for making .f90 version.
	 *>> 1996-06-14 DRSIVA  Krogh  Small change in output format
	 *>> 1994-11-02 DRSIVA  Krogh  Changes to use M77CON
	 *>> 1994-07-18 DRSIVA  Krogh   Last change.
	 *--S replaces "?": DR?IVA, ?IVA, ?IVAF, ?IVAO
	 * Sample driver for SIVA --  Set up to solve two second order equations.
	 * */
	/*--S Next line special: P=>D, X=>Q */
	/*++SP Default NDIG = 4
	 *++  Default NDIG = 10
	 *++ Substitute for NDIG below */
 
 
	/* Set option for error control, local absolute error < TOL. */
	/* Group the system to be treated as a single unit, set tolerance value
	 * Set option for second order equations */
	/*  Flag end of options */
 
	/* Do the integration
	 * */
	Kord[1] = 0;
L_100:
	;
	siva( tspecs, y, f, kord, neq, sivaf, sivao, ITDIM, IYDIM, IFDIM,
	 IKDIM, iopt );
	if (Kord[1] != 1)
		goto L_100;
	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.
	 * */
	/*--S Next line special: P=>D, X=>Q */
 
	/* 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.
	 * */
	/*--S Next line special: P=>D, X=>Q */
 
	/* 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 */