d10r1.cpp

Graph of the Bessel function J0 Graph of the Bessel function J0

#include <process.h>
#include <math.h>
#include <iomanip.h>
#include <iostream.h>

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  func(double x)
{
 double t;
 t=bessj0( x);
 return t;
}

void scrsho()
{
	int iscr,jscr,j,i,jz;
	char scr[61][22],blank,zero,yy,xx,ff;
	double x,dx,x1,x2,ybig,ysml,dyj, y[60];
    iscr = 60;
    jscr = 21;
    
    blank = ' ';
    zero = '-';
    yy = '1';
    xx = '-';
    ff = 'x';
    //enter x1,x2 [x1=x2 to stop]
    x1 = -5.0; 
    x2 = 5.0 ;
    if (x1 == x2) _c_exit();
    for (j = 1 ; j<=jscr; j++)
	{
        scr[1][ j] = yy;
        scr[iscr][ j] = yy;
    }
    for (i = 2; i<=iscr - 1; i++)
	{
        scr[i][ 1] = xx;
        scr[i][ jscr] = xx;
        for (j = 2; j<=(jscr - 1); j++)
		{
            scr[i][ j] = blank;
        }
    }
    dx = double(x2 - x1) / double(iscr - 1);
    x = x1;
    ybig = 0.0 ;
    ysml = ybig;
    for (i = 1; i<=iscr; i++)
	{
        y[i] = func(x);
        if (y[i] < ysml) ysml = y[i];
        if (y[i] > ybig) ybig = y[i];
        x = x + dx;
    }
    if (ybig == ysml) ybig = ysml + 1.0; 
    dyj = (jscr- 1) /(ybig - ysml);
    jz = int(1 - ysml * dyj);
    for (i = 1; i<= iscr; i++)
	{
        scr[i][jz] = zero;
        j = int(1 + (y[i] - ysml) * dyj);
        scr[i][ j] = ff;
    }
    cout<<ybig;
    cout<< " ";
    for (i = 1; i<=iscr; i++)
	{
        cout<< scr[i][jscr];
    }
    cout<<endl;
    for (j = jscr - 1; j>=2; j--)
	{
        cout<< "         ";
         for (i = 1; i<=iscr; i++)
		 {
            cout<< scr[i][j];
		 }
        cout<<endl;
    }
    cout<<ysml, "  .      ";
    cout<< " ";
   for (i = 1; i<=iscr; i++)
   {
        cout<< scr[i][1];
   }
    cout<<endl;
    cout<< "         ";
    cout<< x1<< "                                                          ";
    cout<< x2<< ".       "<<endl;
}

void main()
	{
    //program d10r1
    //driver for routine scrsho

    cout<<  "Graph of the Bessel function J0:";
    cout<< endl;
    scrsho();
	}