/*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 "dbespq.h"
#include <float.h>
#include <stdlib.h>
		/* PARAMETER translations */
#define	ZERO	0.e0
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ dbespq(
double x,
double v,
double *p,
double *q)
{
	long int _l0;
	double a0, a1, a2, b, emu, psum, q1, qsum, sig, term,
	 x8;
	static double small;
	static double one = 1.e0;
	static double two = 2.e0;
	static double four = 4.e0;
	static double eight = 8.e0;
	static double c9 = 9.e0;
	static double xpq = ZERO;
	static double c11293 = 1.1293e0;
	static double cp3 = 0.30103e0;
	static double c59 = 0.59e0;
	static double half = 0.5e0;
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 2001-05-25 DBESPQ Krogh Minor change for making .f90 version.
	 *>> 1998-10-29 DBESPQ Krogh  Moved external statement up for mangle.
	 *>> 1995-11-13 DBESPQ Krogh  Converted SFTRAN to Fortran
	 *>> 1994-10-19 DBESPQ Krogh  Changes to use M77CON
	 *>> 1994-04-19 DBESPQ CLL  Edited to make DP & SP files similar.
	 *>> 1992-03-13 DBESPQ FTK  Removed implicit statements.
	 *>> 1986-03-18 DBESPQ Lawson  Initial code.
	 *--D replaces "?": ?BESPQ, ?ERM1, ?ERV1
	 *
	 *     This subr evaluates asymptotic series for P and Q.
	 *     These can be used to compute Bessel functions by the formulas
	 *
	 *          J = sqrt(2/(pi*X)) * (P * cos(chi) - Q * sin(chi))
	 *
	 *          Y = sqrt(2/(pi*X)) * (P * sin(chi) + Q * cos(chi))
	 *     where
	 *          chi = X - (0.5 * V + 0.25) * pi
	 *
	 *          Reference: NBS AMS55 Eqs 9,2.9 and 9,2.10
	 *
	 *          We assume V is limited to the range [0,2].
	 *          To compute P with a relative accuracy of at least
	 *     10**(-s), X must be restricted to be not less than XPQ,
	 *     where XPQ = 1.1293 * s - 0.59
	 *     (This formula for XPQ was determined for s in the range from
	 *     5 to 25, and limiting V to [0,2].)
	 *          Let s0 = -log10( machine_eps )
	 *     and let s1 = s0 + .3,  s2 = s0 + .6
	 *     We will set XPQ using s2, and then sum the series till a
	 *     term less than 10**(-s1) is reached.
	 *     By setting XPQ using s2 we provide some tolerance to assure
	 *     that the series will contain a term less than 10**(-s1).
	 *        We add the constant term of each series in last to
	 *     reduce the amount of accumulated rounding error.
	 *
	 *     1984 Apr 2, JPL, C. L. Lawson and S. Chan.
	 *     ------------------------------------------------------------------
	 *
	 *     >     Let the terms of these two series be numbered
	 *     1, 3, 5,... in the P series, and 2, 4, 6,... in the Q series.
	 *     For X >> V these terms decrease in magnitude in the order 1, 2,
	 *     3, 4,... to some smallest term, say number N, and then
	 *     following terms increase in magnitude.
	 *     >     For V = 2 and for given X, let N be the number of the
	 *     smallest term, and let SIZE be the magitude of this
	 *     smallest term.  Here are some values of X, N, and -LOG10(SIZE):
	 *
	 *               X = 5    10     15     20     25     30     35
	 *               N =12    22     32     42     52     62     72
	 *     -LOG10(SIZE)= 4.79  9.35  13.81  18.23  22.63  27.02  31.40
	 *     ------------------------------------------------------------------ */
	/*     ----------
	 * */
	/*     ------------------------------------------------------------------
	 * */
	if (xpq == ZERO)
	{
		small = half*DBL_EPSILON/FLT_RADIX;
		xpq = c11293*(cp3 - log10( small )) - c59;
	}
	/*     ------------------------------------------------------------------ */
	if ((x < xpq || v < ZERO) || v > two)
	{
		derm1( "DBESPQ", 1, 0, "Require X .ge. XPQ and V in [0.,2.]."
		 , "X", x, ',' );
		derv1( "V", v, ',' );
		derv1( "XPQ", xpq, '.' );
		*p = ZERO;
		*q = ZERO;
		return;
	}
	/*     ------------------------------------------------------------------ */
	emu = four*(v*v);
	x8 = eight*x;
	a0 = c9;
	a1 = eight;
	a2 = eight;
	b = two;
	sig = -one;
	psum = ZERO;
	term = (emu - one)/x8;
	q1 = term;
	qsum = ZERO;
 
L_20:
	;
	term = term*(emu - a0)/(b*x8);
	if (fabs( term ) <= small)
		goto L_40;
	psum += sig*term;
	a1 += a2;
	a0 += a1;
	b += 1;
 
	term = term*(emu - a0)/(b*x8);
	if (fabs( term ) <= small)
		goto L_40;
	qsum += sig*term;
	a1 += a2;
	a0 += a1;
	b += 1;
 
	sig = -sig;
	goto L_20;
 
L_40:
	;
	*p = half + (half + psum);
	*q = q1 + qsum;
	return;
} /* end of function */