/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:41 */
/*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 "delefi.h"
#include <float.h>
#include <stdlib.h>
void /*FUNCTION*/ delefi(
double phi,
double k,
double *f,
double *e,
long *ierr)
{
	long int _l0;
	double an, cn, dn, hn, k2, l2, p, ph, pn, px, qn, qx,
	 qxdr, r, r0, r2, ri, rj, rk, rlpx, rm, rn, s0, s1, s2, s3, s4,
	 sgn, si, sj, sk, sn, ss, td, tr, ts;
	static double ln2 = 0.693147180559945309417232121458176568075e0;
	static double ln4 = 1.386294361119890618834464242916353136151e0;
	static double pi2 = 1.570796326794896619231321691639751442099e0;
	static double th1 = .7615941559557648881194582826047935904128e0;
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1998-10-29 DELEFI Krogh  Moved external statement up for mangle.
	 *>> 1995-11-17 DELEFI Krogh  Converted SFTRAN to Fortran 77.
	 *>> 1995-10-24 DELEFI Krogh  Removed blanks in numbers for C conversion.
	 *>> 1994-10-19 DELEFI Krogh  Changes to use M77CON
	 *>> 1990-11-27 DELEFI WV Snyder Convert code from NSWC to MATH77
	 *
	 *--D replaces "?": ?ELEFI, ?LNREL, ?ERM1
	 *
	 *-- Begin mask code changes
	 *          REAL ELLIPTIC INTEGRALS OF THE FIRST AND SECOND KINDS
	 *
	 *     PHI [IN]   = ARGUMENT                    (ABS(PHI) .LE. PI/2)
	 *     K [IN]     = MODULUS                       (ABS(K) .LE. 1.0)
	 *     F [OUT[    = ELLIPTIC INTEGRAL OF FIRST KIND = F(PHI, K)
	 *     E [OUT]    = ELLIPTIC INTEGRAL OF SECOND KIND = E(PHI, K)
	 *     IERR [OUT] = ERROR INDICATOR
	 *            0 = NO ERRORS
	 *            1 = ARGUMENT, PHI, > PI/2; SEE AMS 55 17.4.3 OR BYRD AND
	 *                 FRIEDMAN 113.02
	 *            2 = MODULUS, ABS(K), > 1; SEE AMS 55 17.4.15 OR BYRD AND
	 *                 FRIEDMAN 114.01;
	 *            3 = ARGUMENT, PHI, = PI/2 AND MODULUS, K = 1 (F INFINITE).
	 *-- End mask code changes
	 * */
 
	/*     *****     External References     ********************************
	 *
	 * D1MACH  Provides machine constants.
	 * DERM1   Displays an error message.
	 * DLNREL  Computes LOG(1+X) given X.
	 * */
 
	/*     *****     Local Variables     ************************************
	 * */
 
	/*     *****     Data Statements     ************************************
	 *
	 *     LN2 = LN(2)
	 *     LN4 = LN(4)
	 *     TH1 = TANH(1)
	 *     PI2 = PI/2
	 * */
 
	/*     *****     Executable Statements     ******************************
	 * */
	if (phi < 0.0e0)
	{
		sgn = -1.0e0;
		ph = -phi;
	}
	else
	{
		sgn = 1.0e0;
		ph = phi;
	}
	if (ph > pi2)
	{
		derm1( "DELEFI", 1, 0, "ABS(PHI) .GT. PI/2", "PHI", phi, '.' );
		*ierr = 1;
		return;
	}
	if (fabs( k ) > 1.0e0)
	{
		derm1( "DELEFI", 1, 0, "ABS(K) .GT. 1.0", "K", k, '.' );
		*ierr = 2;
		return;
	}
	*ierr = 0;
	if (ph == 0.0e0)
	{
		*f = 0.0e0;
		*e = 0.0e0;
		return;
	}
 
	sn = sin( ph );
	cn = sqrt( 0.5e0 + (0.5e0 - sn*sn) );
 
	k2 = k*k;
	px = fabs( k*sn );
	if (px < th1)
	{
 
		/*     SERIES EXPANSION FOR ABS(K*SIN(PHI)) .LE. TANH(1)
		 * */
		ss = sn*sn;
		pn = 1.0e0;
		qn = 2.0e0;
		an = ph;
		hn = 1.0e0;
		s1 = 0.0e0;
		s2 = 0.0e0;
		tr = ph*ss;
		ts = sn*cn;
L_100:
		;
		an = (pn*an - ts)/qn;
		r = k2*hn/qn;
		s2 += r*an;
		hn = pn*r;
		s0 = s1;
		s1 += hn*an;
		if (fabs( tr ) < fabs( an ))
			goto L_120;
		if (fabs( s1 ) <= fabs( s0 ))
			goto L_120;
		pn = qn + 1.0e0;
		qn = pn + 1.0e0;
		tr *= ss;
		ts *= ss;
		goto L_100;
L_120:
		;
		*f = ph + s1;
		*e = ph - s2;
	}
	else
	{
 
		/*     SERIES EXPANSION FOR ABS(K*SIN(PHI)) .GT. TANH(1)
		 * */
		r2 = 0.5e0 + (0.5e0 - px*px);
		if (r2 == 0.0e0)
		{
			*ierr = 3;
			*f = DBL_MAX;
			*e = 1.0e0;
			return;
		}
		l2 = 0.5e0 + (0.5e0 - k2);
		r = sqrt( r2 );
		si = 1.0e0;
		sj = l2;
		sk = 0.0e0;
		rm = 0.0e0;
		rn = 0.0e0;
		s1 = 0.0e0;
		s2 = 0.0e0;
		s3 = 0.0e0;
		s4 = 0.0e0;
		qx = fabs( k*cn );
		td = qx*r;
		dn = 2.0e0;
L_140:
		;
		pn = (dn - 1.0e0)/dn;
		qn = (dn + 1.0e0)/(dn + 2.0e0);
		ri = pn*si;
		rj = pn*pn*sj;
		rk = sk + 2.0e0/(dn*(dn - 1.0e0));
		r0 = td/dn;
		rm = qn*qn*l2*(rm - r0*ri);
		rn = pn*qn*l2*(rn - r0*si);
		r0 = s3;
		s1 += rj;
		s2 += qn*rj;
		s3 += rm - rj*rk;
		s4 += rn - sj*(pn*rk - 1.0e0/(dn*dn));
		if (s3 >= r0)
			goto L_160;
		si = ri;
		sj = rj*l2;
		sk = rk;
		dn += 2.0e0;
		td *= r2;
		goto L_140;
L_160:
		;
		qxdr = qx/r;
		rlpx = dlnrel( px );
		p = ln4 - 0.5e0*(rlpx + dlnrel( -px )) - dlnrel( qxdr );
		*f = (1.0e0 + s1)*p + qxdr*(rlpx - ln2) + s3;
		*e = (0.5e0 + s2)*l2*p + (1.0e0 - qxdr*(1.0e0 - px)) + s4;
	}
	*f *= sgn;
	*e *= sgn;
	return;
} /* end of function */