double bessj(int n, double x)
{
	int iacc,j,m;
	double temp,bigno,bigni,ax,tox,bj,bjp,bjm,besj,jsum,sum;
    iacc = 40;
    bigno = 10000000000.0;
    bigni = 0.0000000001;
    if (n < 2 )
	{
        cout<<"bad argument n in bassj";
        _c_exit();
    }
    ax = fabs(x);
    if( ax == 0) 
	{
        temp = 0.0;
	}
    else if (ax > float(n))
	{
        tox = 2.0 / ax;
        bjm = bessj0(ax);
        bj = bessj1(ax);
        for( j = 1;j<=n-1;j++)
		{
            bjp = j * tox * bj - bjm;
            bjm = bj;
            bj = bjp;
        }
        temp = bj;
	}
    else
	{
        tox = 2.0 / ax;
        m = 2 * int(((n + int(sqrt(iacc * n)))) / 2);
        besj = 0.0;
        jsum = 0;
        sum = 0.0;
        bjp = 0.0;
        bj = 1.0;
        for (j = m;j>=1;j--)
		{
            bjm = j * tox * bj - bjp;
            bjp = bj;
            bj = bjm;
            if (fabs(bj) > bigno)
			{
                bj = bj * bigni;
                bjp = bjp * bigni;
                besj = besj * bigni;
                sum = sum * bigni;
            }
            if (jsum != 0) sum = sum + bj;
            jsum = 1 - jsum;
            if (j == n) besj = bjp;
        }
        sum = 2.0 * sum - bj;
        temp = besj / sum;		
    }
	return temp;
}