/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:13 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_dxrk8 s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_dxrk8.h"
/*      program DRDXRK8
 *>> 2008-02-24 DRDXRK8  Krogh -- Change dimenstions for new usage.
 *>> 2001-05-25 DRDXRK8  Krogh -- Added comma to format.
 *>> 1997-12-18 DRDXRK8  Krogh   Initial code
 *--D replaces "?": DR?XRK8, ?XRK8, ?XRK8F, ?XRK8O
 * Sample driver for DXRK8 --  Set up to solve simple two body problem
 * with circular motion.
 * */
		/* PARAMETER translations */
#define	LDAT	(32 + 2*NEQ + 2 + 2)
#define	LIDAT	(49 + 2)
#define	LWORK	(9*4 + 6 + 8*4)
#define	NEQ	4
#define	TOL	1.e-10
#define	XRKAE	2.e0
#define	XRKEND	0.e0
#define	XRKTDT	24.e0
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int _i, _r;
	static long int idat[LIDAT];
	double dat[LDAT], work[LWORK];
	static double opt[6]={XRKAE,TOL,XRKTDT,0.e0,6.283185307179586477e0,
	 XRKEND};
	static double ts[3]={0.e0,0.e0,20.e0};
	static double y[NEQ]={1.e0,0.e0,0.e0,1.e0};
	static int _aini = 1;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Dat = &dat[0] - 1;
	long *const Idat = &idat[0] - 1;
	double *const Opt = &opt[0] - 1;
	double *const Ts = &ts[0] - 1;
	double *const Work = &work[0] - 1;
	double *const Y = &y[0] - 1;
		/* end of OFFSET VECTORS */
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		Idat[4] = NEQ;
		Idat[5] = LIDAT;
		Idat[6] = LDAT;
		Idat[7] = LWORK;
		_aini = 0;
	}
 
	/*++S Default SETTOL = 'TOL = 1.E-4'
	 *++  Default SETTOL = 'TOL = 1.D-10'
	 *++  Replace 'TOL = 1.D-10' = SETTOL */
	/*            Abs. Err. tol.   Output at intervals of 2 Pi       End opt. */
	/*      Starting:  T     H;  Final T. */
	/*      Initial Y */
 
	/* Do the integration */
	Idat[1] = 0;
L_100:
	;
	dxrk8( ts, y, opt, idat, dat, work );
	if (Idat[1] != 3)
		goto L_100;
	exit(0);
} /* end of function */
 
void /*FUNCTION*/ dxrk8f(
double *t,
double y[],
double f[],
long idat[])
{
	double rqbi;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const F = &f[0] - 1;
	long *const Idat = &idat[0] - 1;
	double *const Y = &y[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Sample derivative subroutine for use with DXRK8
	 * This evaluates derivatives for a simple two body problem.
	 * */
 
	/* Evaluate the derivatives */
	rqbi = 1.e0/(powi(sqrt( SQ(Y[1]) + SQ(Y[3]) ),3));
	F[1] = Y[2];
	F[2] = -Y[1]*rqbi;
	F[3] = Y[4];
	F[4] = -Y[3]*rqbi;
	return;
} /* end of function */
 
void /*FUNCTION*/ dxrk8o(
double ts[],
double y[],
long idat[],
double dat[])
{
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Dat = &dat[0] - 1;
	long *const Idat = &idat[0] - 1;
	double *const Ts = &ts[0] - 1;
	double *const Y = &y[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Sample output subroutine for use with DXRK8.
	 * This subroutine gives output for a simple two body problem.
	 * */
 
	/* Do the output */
	if (Ts[1] == 0.e0)
	{
		printf("      RESULTS FOR A SIMPLE 2-BODY PROBLEM\n\n        T             U/V           UP/VP\n \n");
	}
	printf("%15.6e%15.6e%15.6e\n               %15.6e%15.6e\n \n", Ts[1], Y[1], Y[2],
	   Y[3], Y[4]);
	return;
} /* end of function */