/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:46 */
/*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 "selefi.h"
#include <float.h>
#include <stdlib.h>
void /*FUNCTION*/ selefi(
float phi,
float k,
float *f,
float *e,
long *ierr)
{
	long int _l0;
	float 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 float ln2 = 0.693147180559945309417232121458176568075e0;
	static float ln4 = 1.386294361119890618834464242916353136151e0;
	static float pi2 = 1.570796326794896619231321691639751442099e0;
	static float th1 = .7615941559557648881194582826047935904128e0;
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1998-10-29 SELEFI Krogh  Moved external statement up for mangle.
	 *>> 1995-11-17 SELEFI Krogh  Converted SFTRAN to Fortran 77.
	 *>> 1995-10-24 SELEFI Krogh  Removed blanks in numbers for C conversion.
	 *>> 1994-10-19 SELEFI Krogh  Changes to use M77CON
	 *>> 1990-11-27 SELEFI WV Snyder Convert code from NSWC to MATH77
	 *
	 *--S 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     ********************************
	 *
	 * R1MACH  Provides machine constants.
	 * SERM1   Displays an error message.
	 * SLNREL  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)
	{
		serm1( "SELEFI", 1, 0, "ABS(PHI) .GT. PI/2", "PHI", phi, '.' );
		*ierr = 1;
		return;
	}
	if (fabsf( k ) > 1.0e0)
	{
		serm1( "SELEFI", 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 = sinf( ph );
	cn = sqrtf( 0.5e0 + (0.5e0 - sn*sn) );
 
	k2 = k*k;
	px = fabsf( 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 (fabsf( tr ) < fabsf( an ))
			goto L_120;
		if (fabsf( s1 ) <= fabsf( 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 = FLT_MAX;
			*e = 1.0e0;
			return;
		}
		l2 = 0.5e0 + (0.5e0 - k2);
		r = sqrtf( 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 = fabsf( 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 = slnrel( px );
		p = ln4 - 0.5e0*(rlpx + slnrel( -px )) - slnrel( 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 */