Bessels Interpolation Formula

1:  #include<stdio.h>  
2:  #include<math.h>  
3:  float factorial(int);  
4:  void main()  
5:  {  
6:       float ax[10],y[10];  
7:       float dy[10][10],u,h,x,a,sum=0,sum2,sum3,l,prod=1,prod2=1;  
8:       int i,j,k,m,n,fact,mid=0;  
9:       int diff1,diff2,num_minus;  
10:       printf("\n\t\t\tBessel Interpolation");  
11:       printf("\n=========================================================================================================================\n");  
12:       printf("Enter the number of terms \n");  
13:       scanf("%d",&n);  
14:       printf("\n=========================================================================================================================\n");  
15:       printf("\t\t\tEnter the array x\n");  
16:       for(i=0;i<n;i++)  
17:       {  
18:            printf("Enter the value of ax[%d] = ",i);  
19:            scanf("%f",&ax[i]);  
20:       }  
21:       printf("\n=========================================================================================================================\n");  
22:       printf("\t\t\tEnter the array Y\n");  
23:       for(i=0;i<n;i++)  
24:       {  
25:            printf("Enter the value of y[%d] = ",i);  
26:            scanf("%f",&y[i]);  
27:       }  
28:       printf("\n=========================================================================================================================\n");  
29:       printf("\nEnter the interpolation point ");  
30:       scanf("%f",&x);  
31:       h=ax[1]-ax[0];  
32:       for(i=0;i<n;i++)  
33:       {  
34:            for(j=0;j<=n;j++)  
35:                 dy[i][j]=-1;  
36:       }  
37:       for(i=0;i<n;i++)  
38:       {  
39:            dy[i][0]=y[i];  
40:       }  
41:       k=1;  
42:       j=1;  
43:       while(j<n)  
44:       {  
45:            i=1;  
46:            while(i<=n-k)  
47:            {  
48:                 dy[i-1][j]=dy[i][j-1]-dy[i-1][j-1];  
49:                 i++;  
50:            }  
51:            j++;  
52:            k++;  
53:       }  
54:       i=0;  
55:       printf("\n=========================================================================================================================\n");  
56:       printf("\ty");  
57:       for(i=1;i<n;i++)  
58:            printf("   d^%dy",i);  
59:       printf("\n");  
60:       for(i=0;i<n;i++)  
61:       {  
62:            for(j=0;j<n-i;j++)  
63:            {  
64:                 if(dy[i][j]==-1.0)  
65:                      printf(" _ ");  
66:                 else  
67:                      printf("  %-7.4f",dy[i][j]);  
68:            }  
69:            printf("\n");  
70:       }  
71:       i=0;  
72:       do{  
73:            i++;  
74:       }while(ax[i]<x);  
75:       i--;  
76:       u=(x-ax[i])/h;  
77:       sum=(y[i]+y[i+1])/2;  
78:       sum2=sum3=0;  
79:       diff1=i;  
80:       diff2=i-1;  
81:       num_minus=1;  
82:       prod2=1;  
83:       prod=1;  
84:       for(i=1;i<n;i++)  
85:       {  
86:            a=u;  
87:            if(i%2!=0)  
88:            {  
89:                 if(diff1<0)  
90:                      prod=0;  
91:                 else  
92:                      prod=(((a-0.5)*prod2)/factorial(i))*dy[diff1][i];  
93:                 sum=sum+prod;  
94:            }  
95:            else  
96:            {  
97:                 prod2=a;  
98:                 if(i>2)  
99:                      prod2=prod2*(u+1);  
100:                 for(j=1;j<=num_minus;j++)  
101:                 {  
102:                      prod2=prod2*(a-j);  
103:                 }  
104:                 num_minus++;  
105:                 prod=prod2/factorial(i);  
106:                 if(diff2<0)  
107:                      sum3=0;  
108:                 else  
109:                      sum3=(dy[diff1][i]+dy[diff2][i])/2;  
110:                 prod=prod*sum3;  
111:                 diff1=diff2;  
112:                 diff2=diff2-1;  
113:                 sum=sum+prod;  
114:            }  
115:       }  
116:       printf("\nValue at %6.4f is %f ",x,sum);  
117:       getchar();  
118:  }  
119:  float factorial(int x)  
120:  {  
121:       int i,j=1;  
122:       for(i=1;i<=x;i++)  
123:            j=j*i;  
124:       return j;  
125:  }