#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>

int kmax, kount;
double dxsav;
int sgn(double x)
{
	int temp;
	if(x>0)temp=1;
	if(x=0)temp=0;
	if(x<0)temp=-1;
	return temp;
}

double bessj1(double x)
{
	double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6;
	double bbb,ccc,aaa,temp,ax,xx,z,y;
    r1 = 72362614232.0;           r2 = -7895059235.0;
    r3 = 242396853.1;            r4 = -2972611.439;
    r5 = 15704.4826;             r6 = -30.16036606;
    s1 = 144725228442.0;          s2 = 2300535178.0;
    s3 = 18583304.74;            s4 = 99447.43394;
    s5 = 376.9991397;            s6 = 1.0;
    p1 = 1.0;                     p2 = 0.00183105;
    p3 = -0.00003516396496;      p4 = 0.000002457520174;
    p5 = -0.000000240337019;
    q1 = 0.04687499995;          q2 = -0.0002002690873;
    q3 = 0.000008449199096;      q4 = -0.00000088228987;
    q5 = 0.000000105787412;
    if (fabs(x) < 8.0 )
	{
        y = x*x;
        aaa = r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))));
        bbb = s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))));
        temp = x * aaa / bbb;
	}
    else
	{
        ax = fabs(x);
        z = 8.0 / ax;
        y = z*z;;
        xx = ax - 2.356194491;
        aaa = p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)));
        bbb = q1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
        ccc = sqrt(0.636619772 / ax);
        temp = ccc * (cos(xx) * aaa - z * sin(xx) * bbb * sgn(x));
	}  
	return temp;
}

double bessj0(double x)
{
	double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6;
	double y,bbb,ccc,aaa,temp,eee,ddd,ax,xx,z;
    p1 = 1.0;                  p2 = -0.001098628627;
    p3 = 0.00002734510407;    p4 = -0.000002073370639;
    p5 = 2.093887211e-07;
    q1 = -0.01562499995;      q2 = 0.0001430488765;
    q3 = -0.000006911147651;  q4 = 7.621095161e-07;
    q5 = -9.34945152e-08;
    r1 = 57568490574.0;        r2 = -13362590354.0;
    r3 = 651619640.7;         r4 = -11214424.18;
    r5 = 77392.33017;         r6 = -184.9052456;
    s1 = 57568490411.0;        s2 = 1029532985.0;
    s3 = 9494680.718;         s4 = 59272.64853;
    s5 = 267.8532712;        s6 = 1.0;
    if (fabs(x) < 8.0)
	{
       y = x * x;
       bbb = y * (r4 + y * (r5 + y * r6));
       aaa = r1 + y * (r2 + y * (r3 + bbb));
       ccc = y * (s3 + y * (s4 + y * (s5 + y * s6)));
       temp = aaa / (s1 + y * (s2 + ccc));
	}
    else
	{
       ax = fabs(x);
       z = 8.0 / ax;
       y = z * z;
       xx = ax - 0.785398164;
       ccc = y * (p3 + y * (p4 + y * p5));
       aaa = p1 + y * (p2 + ccc);
       ddd = y * (q3 + y * (q4 + y * q5));
       eee = z * sin(xx) * (q1 + y * (q2 + ddd));
       temp = sqrt(0.636619772 / ax) * (cos(xx) * aaa - eee);
	}
	return temp;
}

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;
}

void derivs(double x, double y[],double dydx[])
{
    dydx[1] = -y[2];
    dydx[2] = y[1] - (1.0 / x) * y[2];
    dydx[3] = y[2] - (2.0 / x) * y[3];
    dydx[4] = y[3] - (3.0 / x) * y[4];
}

void rk4(double y[], double dydx[], int n, double x, double &h, double yout[])
{
	int i;
    double yt[11], dyt[11], dym[11],hh,h6,xh;
    hh = h * 0.5;
    h6 = h / 6.0;
    xh = x + hh;
    for( i = 1;i<=n;i++)
	{
        yt[i] = y[i] + hh * dydx[i];
    }
    derivs(xh, yt, dyt);
    for( i = 1;i<=n;i++)
	{
        yt[i] = y[i] + hh * dyt[i];
    }
    derivs(xh, yt, dym);
    for( i = 1;i<=n;i++)
	{
        yt[i] = y[i] + h * dym[i];
        dym[i] = dyt[i] + dym[i];
    }
    derivs(x + h, yt, dyt);
    for (i = 1;i<=n;i++)
	{
        yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + 2.0 * dym[i]);
    }
	for(i=1;i<=10;i++)
	{
		dym[i]=0.0;
		dyt[i]=0.0;
		yt[i]=0.0;
	}
}

void rkqc(double y[], double dydx[], int n, double &x, double htry, double eps, double yscal[], double &hdid, double &hnext)
{
	int i,flag;
	double one,safety,errcon,fcor,ytemp[11], ysav[11], dysav[11],xsav,h,hh,pgrow,pshrnk,errmax;
    one = 1.0;
    safety = 0.9;
    errcon = 0.0006;
	fcor = 0.066667;
    pgrow = -0.2;
    pshrnk = -0.25;
    xsav = x;
    for( i = 1;i<=n;i++)
	{
        ysav[i] = y[i];
        dysav[i] = dydx[i];
    }
    h = htry;
    do 
	{
        hh = 0.5 * h;
        rk4(ysav, dysav, n, xsav, hh, ytemp);
        x = xsav + hh;
        derivs(x, ytemp, dydx);
        rk4(ytemp, dydx, n, x, hh, y);
        x = xsav + h;
        if( x == xsav)
		{
            cout<<" stepsize not significant in rkqc."<<endl;
            _c_exit();
        }
        rk4(ysav, dysav, n, xsav, h, ytemp);
        errmax = 0.0;
        for( i = 1;i<=n;i++)
		{
            ytemp[i] = y[i] - ytemp[i];
            if (errmax < fabs(ytemp[i] / yscal[i]))
			{
                errmax = fabs(ytemp[i] / yscal[i]);
			}
        }
        errmax = errmax / eps;
        if (errmax > one)
		{
            h = safety * h * pow(errmax , pshrnk);
            flag = 1;
		}
        else
		{
            hdid = h;
            if( errmax > errcon )
			{
                hnext = safety * h * pow(errmax , pgrow);
			}
            else
			{
                hnext = 4.0 * h;
			}
            flag = 0;
		}
	}
	while( flag == 1);
    for (i=1;i<=n;i++)
	{
        y[i] = y[i] + ytemp[i] * fcor;
    }
	for(i=1;i<=10;i++)
	{
		dysav[i]=0.0;
		ysav[i]=0.0;
		ytemp[i]=0.0;
	}
}

void odeint(double ystart[5], int &nvar, double &x1, double &x2, double &eps, double &h1, double &hmin, double &nok, double &nbad, double xp[201], double yp[11][201])
{
	int i,nstp,maxstp;
	double yscal[11], y[11], dydx[11],two,zero,tiny,x,h,xsav,hdid,hnext;
    maxstp = 10000;
    two = 2.0;
    zero = 0.0;
    tiny = 1e-30;   
    x = x1;
    h = fabs(h1) * sgn(x2 - x1);
    nok = 0;
    nbad = 0;
    kount = 0;
	hdid=0.0;
	hnext=0.0;
    for (i = 1; i<=nvar; i++)
	{
        y[i] = ystart[i];
    }
    if (kmax > 0)
		xsav = x - dxsav * two;
    for( nstp = 1; nstp<=maxstp; nstp++)
	{
        derivs(x, y, dydx);
        for( i = 1; i<=nvar; i++)
		{
            yscal[i] = fabs(y[i]) + fabs(h * dydx[i]) + tiny;
        }
        if( kmax > 0 )
		{
            if( fabs(x - xsav) > fabs(dxsav))
			{
                if( kount < kmax - 1)
				{
                    kount = kount + 1;
                    xp[kount] = x;
                    for( i = 1; i<=nvar; i++)
					{
                        yp[i][kount] = y[i];
                    }
                    xsav = x;
                }
            }
        }
        if( (x + h - x2) * (x + h - x1) > zero )
		{
			h = x2 - x;
		}
        rkqc(y, dydx, nvar, x, h, eps, yscal, hdid, hnext);
        if( hdid == h)
		{
            nok = nok + 1;
		}
        else
		{
            nbad = nbad + 1;
		}
        if( (x - x2) * (x2 - x1) >= zero )
		{
            for (i = 1;i<=nvar;i++)
			{
                ystart[i] = y[i];
            }
            if( kmax != 0 )
			{
                kount = kount + 1;
                xp[kount] = x;
                for( i = 1; i<=nvar; i++)
				{
                    yp[i][kount] = y[i];
                }
            }
			for(i=1; i<=10; i++)
			{
				dydx[i]=0.0;
				y[i]=0.0;
				yscal[i]=0.0;
			}
            goto over;
        }
        if (fabs(hnext) < hmin)
		{
            cout<<" stepsize smaller than minimum."<<endl;
            _c_exit();
        }
        h = hnext;
    }
    cout<<" too many steps."<<endl;
over:;
}

void main()
{
    //program d14r4
    //driver for routine odeint
	int i,nvar;
	double vstart[5], xp[201], yp[10][201],x1,x2,eps,h1,hmin,nok,nbad;
    nvar = 4;
    x1 = 1.0;
    x2 = 10.0;
    vstart[1] = bessj0(x1);
    vstart[2] = bessj1(x1);
    vstart[3] = bessj(2, x1);
    vstart[4] = bessj(3, x1);
    eps = 0.0001;
    h1 = 0.1;
    hmin = 0.0;
    kmax = 100;
    dxsav = (x2 - x1) / 20.0;
    odeint(vstart, nvar, x1, x2, eps, h1, hmin, nok, nbad, xp, yp);
    cout<<endl;
    cout<<"Successful step: "<<nok<<endl;
    cout<<"Bad step: "<<nbad<<endl;
    cout<<"Stored intermediate values: "<<kount<<endl;
    cout<<endl;
    cout<<"       x         Integral      Bessj(3,x)"<<endl;
    for( i = 1; i<=kount; i++)
	{
        cout<<setprecision(4)<<setw(10)<<xp[i];
        cout<<setprecision(6)<<setw(15)<<yp[4][i];
        cout<<setprecision(6)<<setw(15)<<bessj(3, xp[i])<<endl;
    }
}