/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:45 */
/*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 "scpltk.h"
#include <float.h>
#include <stdlib.h>
float /*FUNCTION*/ scpltk(
float em)
{
	long int _l0, j;
	static long int maxt;
	float emp, s1, s2, scpltk_v, u;
	static float fac;
	static float vl = 1.3862943611198906188344642429163531361510002687205e0;
	static float zero = .0e0;
	static float one = 1.e0;
	static float half = .5e0;
	static long m = 0;
	static float a[3][10]={6.63980111468e-03,2.59628884526e-02,2.37612248576e-02,
	 3.15594316275e-02,9.65786196226e-02,0.0e0,0.0e0,0.0e0,0.0e0,0.0e0,
	 3.00725199036864838e-04,3.96847090209897819e-03,1.07959904905916349e-02,
	 1.05899536209893585e-02,7.51938672180838102e-03,8.92664629455646620e-03,
	 1.49420291422820783e-02,3.08851730018997099e-02,9.65735903017425285e-02,
	 0.0e0,1.3930878570066467279e-04,2.2966348983969586869e-03,8.0030039806499853708e-03,
	 9.8489293221768937682e-03,6.8479092826245051197e-03,6.1796274460533176084e-03,
	 8.7898018745550646778e-03,1.4938013532687165242e-02,3.0885146271305189866e-02,
	 9.6573590280856255384e-02};
	static float b[3][10]={1.84723416323e-03,1.87516602769e-02,4.49838755399e-02,
	 7.01487577829e-02,1.24999295975e-01,0.0e0,0.0e0,0.0e0,0.0e0,0.0e0,
	 6.66317524646073151e-05,1.72161470979865212e-03,9.28116038296860419e-03,
	 2.06902400051008404e-02,2.95037293486887130e-02,3.73355466822860296e-02,
	 4.88271550481180099e-02,7.03124954595466082e-02,1.24999999997640658e-01,
	 0.0e0,2.9700280966555612066e-05,9.2155463496324984638e-04,5.9739042991554291551e-03,
	 1.5530941631977203877e-02,2.3931913323110790077e-02,3.0124849012898930266e-02,
	 3.7377739758623604144e-02,4.8828041906862397978e-02,7.0312499739038352054e-02,
	 1.2499999999990808051e-01};
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 2001-06-18 SCPLTK Krogh  Remove a ' ' inside a floating point num.
	 *>> 1998-10-29 SCPLTK Krogh  Moved external statement up for mangle.
	 *>> 1996-03-30 SCPLTK Krogh   Changed MAX to MAXT
	 *>> 1994-10-27 SCPLTK Snyder  Correct missing type declarations
	 *>> 1994-10-20 SCPLTK Krogh  Changes to use M77CON
	 *>> 1994-05-16 SCPLTK Snyder  Make SP and DP codes alike using CHGTYP
	 *>> 1989-06-16 SCPLTK Snyder  Remove implied DO from DATA stmts for CFT
	 *>> 1985-08-02 SCPLTK Lawson  Initial code.
	 *--S replaces "?": ?CPLTK, ?ERM1
	 *
	 *     Complete elliptic integral K(EM). Method given by
	 *     W.J.CODY, Math. of Comp.,Vol.19, (1965), pp. 105-111.
	 *     Data are taken from Cody's table in reverse order.
	 *     We are using N = 5 or 9 for single precision and
	 *     N = 10 for double precision.
	 *
	 *     C.L.LAWSON & STELLA CHAN,JPL,1983 JUNE.
	 *
	 *     ----------------------------------------------------
	 *
	 *     POLY DEGREE     NEGATIVE LOG BASE 10 OF MAX ABS
	 *                     ERROR OF APPROXIMATION
	 *
	 *         N           FOR K           FOR E
	 *         -           -----           -----
	 *         5            9.50            9.44
	 *         9           15.87           15.84
	 *        10           17.45           17.42
	 *
	 *     ----------------------------------------------------
	 * */
 
 
	/*     VL = LOG(4) */
 
 
 
 
 
	if (m == 0)
	{
		if (FLT_EPSILON/FLT_RADIX < 6.31e-17)
		{
			m = 3;
			maxt = 10;
			/*          -LOG10(6.31E-17) = 16.2 */
		}
		else if (FLT_EPSILON/FLT_RADIX < 6.31e-9)
		{
			/*          -LOG10(6.31E-9) = 8.2 */
			m = 2;
			maxt = 9;
		}
		else
		{
			m = 1;
			maxt = 5;
		}
		fac = one + 2*FLT_EPSILON/FLT_RADIX;
	}
	if (em >= one)
	{
		serm1( "SCPLTK", 1, 0, "INFINITE VALUE WHEN EM .EQ. 1 AND UNDEFINED FOR EM .GT. 1"
		 , "EM", em, '.' );
		u = zero;
	}
	else if (em < zero)
	{
		serm1( "SCPLTK", 2, 0, "UNDEFINED FOR EM .LT. ZERO", "EM",
		 em, '.' );
		u = zero;
	}
	else
	{
		emp = half + (half - em);
		s1 = emp*a[m - 1][0];
		s2 = emp*b[m - 1][0];
		for (j = 2; j <= maxt; j++)
		{
			s1 = (s1 + a[m - 1][j - 1])*emp;
			s2 = (s2 + b[m - 1][j - 1])*emp;
		}
		u = vl + s1 + logf( one/emp )*(half + s2);
	}
 
	scpltk_v = u*fac;
	return( scpltk_v );
 
} /* end of function */