/* A Program to compute Derivatives using Bessel's Formula*/
/*Its a executeable file just paste this code in turbo c and run*/

#include <iostream.h>
#include <conio.h>
#include <dos.h>

void main ()
{

int n;
float x[10],h,xp,x0,y[10],p,fp1D,fp2D,sf1,sf2,sf3,sf4,sf5,sf6,sf7,
      sf8,sf1e2,sf2e2,sf1e3,sf1e4,sf3e2,sf2e3,sf4e2,sf4e4,
      sf3e3,sf2e4,sf5e2,sf4e3,sf3e4,sf6e2,sf5e3;

/* Here fp1D & fp2D denotes the 1st & 2nd derivatives of fp at a given point*/
/* s represents here delta*/
/* sf 2,3,4,5,6,7,8 represents the differences of y i.e., y1-y0,y2-y1 etc. */
/* sf1e2 means 1st value of delta square in the table, similarly sf3e2 means
   3rd value of delta square in the table. */
/* similarly sf1e3,sf2e3 means 1st & 2nd value of delta cube in the table. */

clrscr();
cout<<"***** A Program to Compute Derivatives using Bessel's Formula *****\n";
delay (5000);
clrscr();
cout<<"            Muhammad Ali Javed \n\n";
delay (5000);
cout<<"                   040267      \n";
delay (5000);
clrscr();
cout<<"			NOTE \n\n";
cout<<"Here fp1D & fp2D denotes the 1st & 2nd derivatives of fp at a given poiny\n\n";
cout<<"s represents here delta\n\n";
cout<<"sf 2,3,4,5,6,7,8 represents the differences of y i.e., y1-y0,y2-y1 etc.\n\n";
cout<<"sf1e2 means 1st value of delta square in the table, similarly sf3e2 means\
       3rd value of delta square in the table. \n\n";
cout<<"similarly sf1e3,sf2e3 means 1st & 2nd value of delta cube in the table.\n\n";
cout<<"Press Enter to continue\n";
getch();
start:
clrscr();
cout<<"Enter the value of n between 4 and 8.";
cout<<"\n n = ";
cin>>n;
if (n<4||n>8)
{
	cout<<"Invalid entry, try again.";
	getch();
	goto start;
}
else
{

for (int i=0;i!=n;i++)
{
cout<<"\n Enter x"<<i+1<<": \t";
cin>>x[i];
cout<<"\n Enter y("<<i+1<<"): \t";
cin>>y[i];
h=x[1]-x[0];
}
if (n==4)
{
sf1=y[1]-y[0];
cout<<"sf1=\t"<<sf1<<"\n";
sf2=y[2]-y[1];
cout<<"sf2=\t"<<sf2<<"\n";
sf3=y[3]-y[2];
cout<<"sf3=\t"<<sf3<<"\n";
sf1e2=sf2-sf1;
cout<<"\t sf1e2="<<sf1e2<<"\n";
sf2e2=sf3-sf2;
cout<<"\t sf2e2="<<sf2e2<<"\n";
sf1e3=sf2e2-sf1e2;
cout<<"\t\t sf1e3="<<sf1e3<<"\n";
sf1e4=00;
cout<<"Enter the value of xp=\t";
cin>>xp;
if (xp>x[1]&&xp<x[2])
{
 x0=x[1];
 cout<<"x0=\t"<<x0<<endl;

p=(xp-x0)/h;
cout<<"Value of p is= "<< p<<endl;
fp1D=(1/h)*(sf2+(.25*(2*p-1)*(sf1e2+sf2e2))+((1/6)*(3*p*p-3*p+.5)*sf1e3)+((1/24)*(2*p*p*p-3*p*p-p+2)*0));
cout<<"1st Derivative is:\t"<<fp1D<<endl;
fp2D=(1/(h*h))*((.5*(sf1e2+sf2e2))+(.5*(2*p-1)*sf1e3)+((1/24)*(6*p*p-6*p-1)*0));
cout<<"2nd Derivative is:\t"<<fp2D<<endl;
}
if (xp<x[1]||xp>x[2])
{
goto start;
}
}

if (n==5)
{
sf1=y[1]-y[0];
cout<<"sf1=\t"<<sf1<<"\n";
sf2=y[2]-y[1];
cout<<"sf2=\t"<<sf2<<"\n";
sf3=y[3]-y[2];
cout<<"sf3=\t"<<sf3<<"\n";
sf4=y[4]-y[3];
cout<<"sf4=\t"<<sf4<<"\n";
sf1e2=sf2-sf1;
cout<<"\t sf1e2="<<sf1e2<<"\n";
sf2e2=sf3-sf2;
cout<<"\t sf2e2="<<sf2e2<<"\n";
sf3e2=sf4-sf3;
cout<<"\t sf3e2="<<sf3e2<<"\n";
sf1e3=sf2e2-sf1e2;
cout<<"\t\t sf1e3="<<sf1e3<<"\n";
sf2e3=sf3e2-sf2e2;
cout<<"\t\t sf2e3="<<sf2e3<<"\n";
sf1e4=sf2e3-sf1e3;
cout<<"\t\t\t sf1e4="<<sf1e4<<"\n";
cout<<"Enter the value of xp=\t";
cin>>xp;
if (xp>x[1]&&xp<x[2])
{
x0=x[1];
cout<<"x0=\t"<<x0<<endl;
p=(xp-x0)/h;
cout<<"Value of p is= "<<p<<endl;
fp1D=(1/h)*(sf2+(.25*(2*p-1)*(sf1e2+sf2e2))+((1/6)*(3*p*p-3*p+.5)*sf1e3)+((1/24)*(2*p*p*p-3*p*p-p+2)*sf1e4));
cout<<"1st Derivative is:\t"<<fp1D<<endl;
fp2D=(1/(h*h))*((.5*(sf1e2+sf2e2))+(.5*(2*p-1)*sf1e3)+((1/24)*(6*p*p-6*p-1)*sf1e4));
cout<<"2nd Derivative is:\t"<<fp2D<<endl;
}
if (xp>x[2]&&xp<x[3])
{
x0=x[2];
cout<<"x0=\t"<<x0<<endl;
p=(xp-x0)/h;
cout<<"Value of p is="<<p<<endl;
fp1D=(1/h)*(sf3+(.25*(2*p-1)*(sf2e2+sf3e2))+((1/6)*(3*p*p-3*p+.5)*sf2e3)+((1/24)*(2*p*p*p-3*p*p-p+2)*(sf1e4)));
cout<<"1st Derivative is:\t"<<fp1D<<endl;
fp2D=(1/(h*h))*((.5*(sf2e2+sf3e2))+(.5*(2*p-1)*sf2e3)+((1/24)*(6*p*p-6*p-1)*(sf1e4)));
cout<<"2nd Derivative is:\t"<<fp2D<<endl;
}
if (xp<x[1]||xp>x[3])
{
goto start;
}
}

if (n==6)
{
sf1=y[1]-y[0];
cout<<"sf1=\t"<<sf1<<"\n";
sf2=y[2]-y[1];
cout<<"sf2=\t"<<sf2<<"\n";
sf3=y[3]-y[2];
cout<<"sf3=\t"<<sf3<<"\n";
sf4=y[4]-y[3];
cout<<"sf4=\t"<<sf4<<"\n";
sf5=y[5]-y[4];
cout<<"sf5=\t"<<sf5<<"\n";
sf1e2=sf2-sf1;
cout<<"\t sf1e2="<<sf1e2<<"\n";
sf2e2=sf3-sf2;
cout<<"\t sf2e2="<<sf2e2<<"\n";
sf3e2=sf4-sf3;
cout<<"\t sf3e2="<<sf3e2<<"\n";
sf4e2=sf5-sf4;
cout<<"\t sf4e2="<<sf4e2<<"\n";
sf1e3=sf2e2-sf1e2;
cout<<"\t\t sf1e3="<<sf1e3<<"\n";
sf2e3=sf3e2-sf2e2;
cout<<"\t\t sf2e3="<<sf2e3<<"\n";
sf3e3=sf4e2-sf3e2;
cout<<"\t\t sf3e3="<<sf3e3<<"\n";
sf1e4=sf2e3-sf1e3;
cout<<"\t\t\t sf1e4="<<sf1e4<<"\n";
sf2e4=sf3e3-sf2e3;
cout<<"\t\t\t sf2e4="<<sf2e4<<"\n";
cout<<"Enter the value of xp=\t";
cin>>xp;
if (xp>x[2]&&xp<x[3])
{
x0=x[2];
cout<<"x0=\t"<<x0<<endl;
p=(xp-x0)/h;
cout<<"Value of p is="<<p<<endl;
fp1D=(1/h)*(sf3+(.25*(2*p-1)*(sf2e2+sf3e2))+((1/6)*(3*p*p-3*p+.5)*sf2e3)+((1/24)*(2*p*p*p-3*p*p-p+2)*(sf1e4+sf2e4)));
cout<<"1st Derivative is:\t"<<fp1D<<endl;
fp2D=(1/(h*h))*((.5*(sf2e2+sf3e2))+(.5*(2*p-1)*sf2e3)+((1/24)*(6*p*p-6*p-1)*(sf1e4+sf2e4)));
cout<<"2nd Derivative is:\t"<<fp2D<<endl;
}
if (xp>x[3]&&xp<x[4])
{
x0=x[3];
cout<<"x0=\t"<<x0<<endl;
p=(xp-x0)/h;
cout<<"Value of p is="<<p<<endl;
fp1D=(1/h)*(sf4+(.25*(2*p-1)*(sf3e2+sf4e2))+((1/6)*(3*p*p-3*p+.5)*sf3e3)+((1/24)*(2*p*p*p-3*p*p-p+2)*(sf2e4)));
cout<<"1st Derivative is:\t"<<fp1D<<endl;
fp2D=(1/(h*h))*((.5*(sf3e2+sf4e2))+(.5*(2*p-1)*sf3e3)+((1/24)*(6*p*p-6*p-1)*(sf2e4)));
cout<<"2nd Derivative is:\t"<<fp2D<<endl;
}
if (xp<x[2]||xp>x[4])
{
goto start;
}
}

if (n==7)
{
sf1=y[1]-y[0];
cout<<"sf1=\t"<<sf1<<"\n";
sf2=y[2]-y[1];
cout<<"sf2=\t"<<sf2<<"\n";
sf3=y[3]-y[2];
cout<<"sf3=\t"<<sf3<<"\n";
sf4=y[4]-y[3];
cout<<"sf4=\t"<<sf4<<"\n";
sf5=y[5]-y[4];
cout<<"sf5=\t"<<sf5<<"\n";
sf6=y[6]-y[5];
cout<<"sf6=\t"<<sf6<<"\n";
sf1e2=sf2-sf1;
cout<<"\t sf1e2="<<sf1e2<<"\n";
sf2e2=sf3-sf2;
cout<<"\t sf2e2="<<sf2e2<<"\n";
sf3e2=sf4-sf3;
cout<<"\t sf3e2="<<sf3e2<<"\n";
sf4e2=sf5-sf4;
cout<<"\t sf4e2="<<sf4e2<<"\n";
sf5e2=sf6-sf5;
cout<<"\t sf5e2="<<sf5e2<<"\n";
sf1e3=sf2e2-sf1e2;
cout<<"\t\t sf1e3="<<sf1e3<<"\n";
sf2e3=sf3e2-sf2e2;
cout<<"\t\t sf2e3="<<sf2e3<<"\n";
sf3e3=sf4e2-sf3e2;
cout<<"\t\t sf3e3="<<sf3e3<<"\n";
sf4e3=sf5e2-sf4e2;
cout<<"\t\t sf4e3="<<sf4e3<<"\n";
sf1e4=sf2e3-sf1e3;
cout<<"\t\t\t sf1e4="<<sf1e4<<"\n";
sf2e4=sf3e3-sf2e3;
cout<<"\t\t\t sf2e4="<<sf2e4<<"\n";
sf3e4=sf4e3-sf3e3;
cout<<"\t\t\t sf3e4="<<sf3e4<<"\n";
cout<<"Enter the value of xp=\t";
cin>>xp;
if (xp>x[2]&&xp<x[3])
{
x0=x[2];
cout<<"x0=\t"<<x0<<endl;
p=(xp-x0)/h;
cout<<"Value of p is="<<p<<endl;
fp1D=(1/h)*(sf3+(.25*(2*p-1)*(sf2e2+sf3e2))+((1/6)*(3*p*p-3*p+.5)*sf2e3)+((1/24)*(2*p*p*p-3*p*p-p+2)*(sf1e4+sf2e4)));
cout<<"1st Derivative is:\t"<<fp1D<<endl;
fp2D=(1/(h*h))*((.5*(sf2e2+sf3e2))+(.5*(2*p-1)*sf2e3)+((1/24)*(6*p*p-6*p-1)*(sf1e4+sf2e4)));
cout<<"2nd Derivative is:\t"<<fp2D<<endl;
}
if (xp>x[3]&&xp<x[4])
{
x0=x[3];
cout<<"x0=\t"<<x0<<endl;
p=(xp-x0)/h;
cout<<"Value of p is="<<p<<endl;
fp1D=(1/h)*(sf4+(.25*(2*p-1)*(sf3e2+sf4e2))+((1/6)*(3*p*p-3*p+.5)*sf3e3)+((1/24)*(2*p*p*p-3*p*p-p+2)*(sf2e4+sf3e4)));
cout<<"1st Derivative is:\t"<<fp1D<<endl;
fp2D=(1/(h*h))*((.5*(sf3e2+sf4e2))+(.5*(2*p-1)*sf3e3)+((1/24)*(6*p*p-6*p-1)*(sf2e4+sf3e4)));
cout<<"2nd Derivative is:\t"<<fp2D<<endl;
}
if (xp<x[2]||xp>x[4])
{
goto start;
}
}

if (n==8)
{
sf1=y[1]-y[0];
cout<<"sf1=\t"<<sf1<<"\n";
sf2=y[2]-y[1];
cout<<"sf2=\t"<<sf2<<"\n";
sf3=y[3]-y[2];
cout<<"sf3=\t"<<sf3<<"\n";
sf4=y[4]-y[3];
cout<<"sf4=\t"<<sf4<<"\n";
sf5=y[5]-y[4];
cout<<"sf5=\t"<<sf5<<"\n";
sf6=y[6]-y[5];
cout<<"sf6=\t"<<sf6<<"\n";
sf7=y[7]-y[6];
cout<<"sf7=\t"<<sf7<<"\n";
sf1e2=sf2-sf1;
cout<<"\t sf1e2="<<sf1e2<<"\n";
sf2e2=sf3-sf2;
cout<<"\t sf2e2="<<sf2e2<<"\n";
sf3e2=sf4-sf3;
cout<<"\t sf3e2="<<sf3e2<<"\n";
sf4e2=sf5-sf4;
cout<<"\t sf4e2="<<sf4e2<<"\n";
sf5e2=sf6-sf5;
cout<<"\t sf5e2="<<sf5e2<<"\n";
sf6e2=sf7-sf6;
cout<<"\t sf6e2="<<sf6e2<<"\n";
sf1e3=sf2e2-sf1e2;
cout<<"\t\t sf1e3="<<sf1e3<<"\n";
sf2e3=sf3e2-sf2e2;
cout<<"\t\t sf2e3="<<sf2e3<<"\n";
sf3e3=sf4e2-sf3e2;
cout<<"\t\t sf3e3="<<sf3e3<<"\n";
sf4e3=sf5e2-sf4e2;
cout<<"\t\t sf4e3="<<sf4e3<<"\n";
sf5e3=sf6e2-sf5e2;
cout<<"\t\t sf5e3="<<sf5e3<<"\n";
sf1e4=sf2e3-sf1e3;
cout<<"\t\t\t sf1e4="<<sf1e4<<"\n";
sf2e4=sf3e3-sf2e3;
cout<<"\t\t\t sf2e4="<<sf2e4<<"\n";
sf3e4=sf4e3-sf3e3;
cout<<"\t\t\t sf3e4="<<sf3e4<<"\n";
sf4e4=sf5e3-sf4e3;
cout<<"\t\t\t sf4e4="<<sf4e4<<"\n";
cout<<"Enter the value of xp=\t";
cin>>xp;
if (xp>x[3]&&xp<x[4])
{
x0=x[3];
cout<<"x0=\t"<<x0<<endl;
p=(xp-x0)/h;
cout<<"Value of p is="<<p<<endl;
fp1D=(1/h)*(sf4+(.25*(2*p-1)*(sf3e2+sf4e2))+((1/6)*(3*p*p-3*p+.5)*sf3e3)+((1/24)*(2*p*p*p-3*p*p-p+2)*(sf2e4+sf3e4)));
cout<<"1st Derivative is:\t"<<fp1D<<endl;
fp2D=(1/(h*h))*((.5*(sf3e2+sf4e2))+(.5*(2*p-1)*sf3e3)+((1/24)*(6*p*p-6*p-1)*(sf2e4+sf3e4)));
cout<<"2nd Derivative is:\t"<<fp2D<<endl;
}
if (xp>x[4]&&xp<x[5])
{
x0=x[4];
cout<<"x0=\t"<<x0<<endl;
p=(xp-x0)/h;
cout<<"Value of p is="<<p<<endl;
fp1D=(1/h)*(sf5+(.25*(2*p-1)*(sf4e2+sf5e2))+((1/6)*(3*p*p-3*p+.5)*sf4e3)+((1/24)*(2*p*p*p-3*p*p-p+2)*(sf3e4+sf4e4)));
cout<<"1st Derivative is:\t"<<fp1D<<endl;
fp2D=(1/(h*h))*((.5*(sf4e2+sf5e2))+(.5*(2*p-1)*sf4e3)+((1/24)*(6*p*p-6*p-1)*(sf3e4+sf4e4)));
cout<<"2nd Derivative is:\t"<<fp2D<<endl;
}
if (xp<x[3]||xp>x[5])
{
goto start;
}
}

getch();
}
}