/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:04 */
/*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 "ddasin.h"
#include <stdlib.h>
void /*FUNCTION*/ ddasin(
double x,
double xout,
double yout[],
double ypout[],
long neq,
long kold,
double *phi,
double psi[])
{
#define PHI(I_,J_)	(*(phi+(I_)*(neq)+(J_)))
	long int i, j, koldp1;
	double c, d, gamma, temp1;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Psi = &psi[0] - 1;
	double *const Yout = &yout[0] - 1;
	double *const Ypout = &ypout[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 2006, Math a la Carte, Inc.
	 *>> 2001-11-23 ddasin Krogh  Changed many names per library conventions.
	 *>> 2001-11-04 ddasin Krogh  Fixes for F77 and conversion to single
	 *>> 2001-11-01 ddasin Hanson Provide code to Math a la Carte.
	 *--D replaces "?": ?dasin, ?daslx, ?dastp
	 *      IMPLICIT NONE
	 ****BEGIN PROLOGUE  DDASIN
	 ****SUBSIDIARY
	 ****PURPOSE  Interpolation routine for DDASLX.
	 ****LIBRARY   SLATEC (DDASLX)
	 ****TYPE      DOUBLE PRECISION (SDASIN-S, DDASIN-D)
	 ****AUTHOR  Petzold, Linda R., (LLNL)
	 ****DESCRIPTION
	 * ----------------------------------------------------------------------
	 *     THE METHODS IN SUBROUTINE DDASTP USE POLYNOMIALS
	 *     TO APPROXIMATE THE SOLUTION. DDASIN APPROXIMATES THE
	 *     SOLUTION AND ITS DERIVATIVE AT TIME XOUT BY EVALUATING
	 *     ONE OF THESE POLYNOMIALS, AND ITS DERIVATIVE,THERE.
	 *     INFORMATION DEFINING THIS POLYNOMIAL IS PASSED FROM
	 *     DDASTP, SO DDASIN CANNOT BE USED ALONE.
	 *
	 *     THE PARAMETERS ARE:
	 *     X     THE CURRENT TIME IN THE INTEGRATION.
	 *     XOUT  THE TIME AT WHICH THE SOLUTION IS DESIRED
	 *     YOUT  THE INTERPOLATED APPROXIMATION TO Y AT XOUT
	 *           (THIS IS OUTPUT)
	 *     YPOUT THE INTERPOLATED APPROXIMATION TO YPRIME AT XOUT
	 *           (THIS IS OUTPUT)
	 *     NEQ   NUMBER OF EQUATIONS
	 *     KOLD  ORDER USED ON LAST SUCCESSFUL STEP
	 *     PHI   ARRAY OF SCALED DIVIDED DIFFERENCES OF Y
	 *     PSI   ARRAY OF PAST STEPSIZE HISTORY
	 * ----------------------------------------------------------------------
	 ****ROUTINES CALLED  (NONE)
	 ****REVISION HISTORY  (YYMMDD)
	 *   830315  DATE WRITTEN
	 *   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
	 *   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
	 *   901026  Added explicit declarations for all variables and minor
	 *           cosmetic changes to prologue.  (FNF)
	 ****END PROLOGUE  DDASIN
	 * */
 
 
	/****FIRST EXECUTABLE STATEMENT  DDASIN */
	koldp1 = kold + 1;
	temp1 = xout - x;
	for (i = 1; i <= neq; i++)
	{
		Yout[i] = PHI(0,i - 1);
		Ypout[i] = 0.0e0;
	}
	c = 1.0e0;
	d = 0.0e0;
	gamma = temp1/Psi[1];
	for (j = 2; j <= koldp1; j++)
	{
		d = d*gamma + c/Psi[j - 1];
		c *= gamma;
		gamma = (temp1 + Psi[j - 1])/Psi[j];
		for (i = 1; i <= neq; i++)
		{
			Yout[i] += c*PHI(j - 1,i - 1);
			Ypout[i] += d*PHI(j - 1,i - 1);
		}
	}
	return;
 
	/* -----END OF SUBROUTINE DDASIN------ */
#undef	PHI
} /* end of function */