/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:40 */
/*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 "dci.h"
#include <stdlib.h>
double /*FUNCTION*/ dci(
double x)
{
	double dci_v;
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1998-10-29 DCI Krogh  Moved external statement up for mangle.
	 *>> 1996-03-30 DCI Krogh  Added external statements.
	 *>> 1995-11-22 DCI Krogh  Removed multiple entry for C conversion.
	 *>> 1995-11-03 DCI Krogh  Removed blanks in numbers for C conversion.
	 *>> 1994-10-20 DCI Krogh  Changes to use M77CON
	 *>> 1990-01-23 DCI CLL Using name DCIN for result in entry DCIN.
	 *>> 1989-03-14 Original W. V. Snyder at JPL
	 *
	 *     COMPUTE THE COSINE INTEGRAL OF X =
	 *     INTEGRAL FROM X TO INFINITY OF -(COS(T)/T) DT =
	 *     GAMMA + LOG(ABS(X)) + INTEGRAL FROM 0 TO X OF ((COS(T)-1)/T DT),
	 *     WHERE GAMMA IS EULER'S CONSTANT.
	 *
	 *     FOR ABS(X)<16, USE A CHEBYSHEV SERIES WITH ARGUMENT 2*Z*Z-1 WHERE
	 *     Z=X/16 TO EVALUATE (CI(X)-LOG(ABS(X))-GAMMA)/(Z*Z), THEN MULTIPLY
	 *     THE RESULT BY Z*Z AND ADD LOG(ABS(X)) AND GAMMA.  LOSS OF ABOUT
	 *     TWO DIGITS OCCURS NEAR ABS(X)=16 DUE TO CANCELLATION.
	 *
	 *     FOR ABS(X).GE.16, USE CHEBYSHEV SERIES WITH ARGUMENT 2*Z*Z-1 WHERE
	 *     WHERE Z=16/X ARE USED TO COMPUTE F(X)/Z AND G(X)/(Z*Z).  THEN
	 *     CI(X)=F(X)*SIN(X)-G(X)*COS(X).
	 *
	 *     THIS ALGORITHM YIELDS AT MOST 15 DIGITS OF PRECISION.
	 *
	 *--D replaces "?": ?CI, ?CII, ?CIN, ?CPVAL, ?ERM1
	 * */
	dci_v = dcii( TRUE, x );
	return( dci_v );
} /* end of function */
 
double /*FUNCTION*/ dcin(
double x)
{
	double dcin_v;
 
	dcin_v = dcii( FALSE, x );
	return( dcin_v );
} /* end of function */
 
		/* PARAMETER translations */
#define	ERRLEV	0
#define	GAMMA	0.577215664901533e0
		/* end of PARAMETER translations */
 
double /*FUNCTION*/ dcii(
LOGICAL32 ldci,
double x)
{
	double dcii_v, z, zw;
	static double c[23]={0.5e0,0.5e0,14.992589367813409e0,-19.386124096607770e0,
	 12.741870869758071e0,-8.107903970562531e0,4.862022348500627e0,
	 -2.497505088539025e0,1.008660787358110e0,-0.312080924825428e0,
	 0.074678255294576e0,-0.014110865253535e0,0.002152046752074e0,
	 -0.000270212331184e0,0.000028416945498e0,-0.000002540125611e0,
	 0.000000195437144e0,-0.000000013084020e0,0.000000000769379e0,
	 -0.000000000040066e0,0.000000000001861e0,-0.000000000000078e0,
	 0.000000000000003e0};
	static double f[13]={0.5e0,0.5e0,0.062263729028927e0,-0.000233756041393e0,
	 0.000002453755677e0,-0.000000058670317e0,0.000000002356196e0,
	 -0.000000000136096e0,0.000000000010308e0,-0.000000000000964e0,
	 0.000000000000107e0,-0.000000000000014e0,0.000000000000002e0};
	static double g[13]={0.5e0,0.5e0,0.003862856096703e0,-0.000042644182622e0,
	 0.000000724995950e0,-0.000000023468225e0,0.000000001169202e0,
	 -0.000000000079604e0,0.000000000006875e0,-0.000000000000717e0,
	 0.000000000000087e0,-0.000000000000012e0,0.000000000000002e0};
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const C = &c[0] - 1;
	double *const F = &f[0] - 1;
	double *const G = &g[0] - 1;
		/* end of OFFSET VECTORS */
 
 
	if (ldci)
	{
		if (x <= 0.0)
		{
			derm1( "DCI", 1, ERRLEV, "Argument not positive", "X",
			 x, '.' );
			/*        Provide a value in case the error message processor returns. */
			dcii_v = 0.0;
		}
		else if (x < 16.0)
		{
			z = x/16.0;
			zw = z*z;
			dcii_v = log( x ) - zw*dcpval( c, 20, zw ) + GAMMA;
		}
		else
		{
			z = 16.0/x;
			zw = z*z;
			dcii_v = z*(dcpval( f, 10, zw )*sin( x ) - z*dcpval( g,
			 10, zw )*cos( x ));
		}
	}
	else
	{
 
		/*     Evaluate the entire function
		 *     Cin(x) = Integral from 0 to x of ((cos(t) - 1) / t) dt.
		 * */
		if (fabs( x ) < 16.0)
		{
			z = x/16.0;
			zw = z*z;
			dcii_v = zw*dcpval( c, 20, zw );
		}
		else
		{
			z = 16.0/x;
			zw = z*z;
			dcii_v = log( fabs( x ) ) + GAMMA - z*(dcpval( f, 10,
			 zw )*sin( x ) - z*dcpval( g, 10, zw )*cos( x ));
		}
	}
	return( dcii_v );
} /* end of function */