/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:56 */
/*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 "ssi.h"
#include <stdlib.h>
float /*FUNCTION*/ ssi(
float x)
{
	long int n;
	float ft, fz, gt, gz, ssi_v, z, zw;
	static float pi2 = 1.57079632679489662e0;
	static float s[23]={0.5e0,0.5e0,4.052926477680623e0,-4.063980844911986e0,
	 2.778756381742663e0,-1.926565091150656e0,1.389308771171888e0,
	 -0.968322236987086e0,0.530148847916522e0,-0.211263780976555e0,
	 0.062033679432003e0,-0.013867445589417e0,0.002436221404749e0,
	 -0.000345469155569e0,0.000040420271419e0,-0.000003972908746e0,
	 0.000000332988589e0,-0.000000024100076e0,0.000000001522370e0,
	 -0.000000000084710e0,0.000000000004185e0,-0.000000000000185e0,
	 0.000000000000007e0};
	static float 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 float 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 */
	float *const F = &f[0] - 1;
	float *const G = &g[0] - 1;
	float *const S = &s[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1998-10-29 SSI Krogh  Moved external statement up for mangle.
	 *>> 1995-11-03 SSI Krogh  Removed blanks in numbers for C conversion.
	 *>> 1994-11-11 SSI Krogh  Declared all vars.
	 *>> 1994-10-20 SSI Krogh  Changes to use M77CON
	 *>> 1989-03-14 SSI Original W. V. Snyder at JPL
	 *
	 *     COMPUTE THE SINE INTEGRAL OF X =
	 *     INTEGRAL FROM 0 TO X OF (SIN(T)/T DT).
	 *
	 *     FOR ABS(X)<16, USE A CHEBYSHEV SERIES WITH ARGUMENT 2*Z*Z-1 WHERE
	 *     Z=X/16 TO EVALUATE SI(X)/Z, THEN MULTIPLY THE RESULT BY Z.  THIS
	 *     AVOIDS STORING ZERO COEFFICIENTS FOR EVEN ORDERS, AND PRESERVES
	 *     ACCURACY FOR SMALL Z.
	 *
	 *     FOR 16.LE.ABS(X).LE.100, USE CHEBYCHEV SERIES WITH ARGUMENT
	 *     2*Z*Z-1, WHERE Z=16/X ARE USED TO COMPUTE F(X)/X AND G(X)/(X*X).
	 *     THEN SI(X)=0.5*PI*SIGN(X)-F(X)/X*COS(X)-G(X)/(X*X)*SIN(X).
	 *
	 *     WHEN X.GT.100, USE ASYMPTOTIC APPROXIMATIONS FOR F(X)/X AND
	 *     G(X)/(X*X) AND COMPUTE SI(X) AS ABOVE.
	 *
	 *     THIS ALGORITHM YIELDS AT MOST 15 DIGITS OF PRECISION.
	 *
	 *--S replaces "?": ?SI, ?CPVAL
	 * */
 
	if (fabsf( x ) < 16.0)
	{
		z = x/16.0;
		zw = z*z;
		z *= scpval( s, 20, zw );
	}
	else
	{
		if (fabsf( x ) <= 100.0)
		{
			/*           16.LE.ABS(X).LE.100 */
			z = 16.0/x;
			zw = z*z;
			fz = z*scpval( f, 10, zw );
			gz = zw*scpval( g, 10, zw );
		}
		else
		{
			/*           ABS(X).GT.100 */
			fz = 1.0/x;
			ft = fz;
			gz = fz/x;
			gt = gz;
			z = gz;
			for (n = 2; n <= 16; n += 2)
			{
				ft *= -(float)( n*(n - 1) )*z;
				gt *= -(float)( n*(n + 1) )*z;
				fz += ft;
				gz += gt;
			}
		}
		z = signf( pi2, x ) - fz*cosf( x ) - gz*sinf( x );
	}
 
	ssi_v = z;
	return( ssi_v );
} /* end of function */