в программе ошибка,найти не могу

#include<stdio.h>
#include<conio.h>
#include<complex.h>
#include<stdlib.h>
#include<graphics.h>
 
#define K 7   
#define M 50  
#define N 16 
#define MM 8  
 
complex iii=complex(0.0,1.0);
double w;
complex const_w;
double r1,r2,h=0.5,L1,L2,L3,M1,M2,M3,R1,R2,R3;
int zad;  
double eps=0.000001; 
double C=0.577215664901532; 
 
double module(complex x)
  { return(sqrt(real(x)*real(x)+imag(x)*imag(x))); }
 
double fact(double n)
  {
   double i,k;
 
   k=1.0;
   for(i=1.0;i<(n+1.0);i++)
     k=k*i;
   return(k);
  }
 
complex J(double x,double n)
  {
   double sum,s;
   double k;
 
   if(n<0.0) return(pow(-1.0,-n)*J(x,-n));
   else
    {
     if(n>1.0) return(2.0*(n-1.0)/x*J(x,n-1.0)-J(x,n-2.0));
     if(n==0.0)
       {
    n=0.0;
    k=-1.0;
    sum=0.0;
    s=0.0;
    do{
       k=k+1.0;
       sum=sum+s;
       s=pow(-1.0,k)/fact(k)/fact(n+k)*pow(x/2.0,2*k+n);
      }while(module(s)>=eps);
    return(sum);
       }
     if(n==1.0)
       {
    n=1.0;
    k=-1.0;
    sum=0.0;
    s=0.0;
    do{
       k=k+1.0;
       sum=sum+s;
       s=pow(-1.0,k)/fact(k)/fact(n+k)*pow(x/2.0,2*k+n);
      }while(module(s)>=eps);
    return(sum);
       }
    }
  }
 
complex J_der(double x,double n)
  { return((J(x,n-1.0)-J(x,n+1.0))/2.0); }
 
complex J_der_der(double x,double n)
  { return((J_der(x,n-1.0)-J_der(x,n+1.0))/2.0); }
 
complex Ne(double x,double n)
  {
   complex sum1,sum2,sum3,s,ss,sss;
   double k,nn,i;
 
   if(n<0.0) return(pow(-1.0,-n)*Ne(x,-n));
   else
    {
     if(n>1.0) return(2.0*(n-1.0)/x*Ne(x,n-1.0)-Ne(x,n-2.0));
     if(n==0.0)
       {
    n=0.0;
    sum1=2.0/M_PI*(C+log(x/2.0))*J(x,n);
    sum2=0.0;
    for(k=0.0;k<n;k=k+1.0)
      sum2=sum2+fact(n-k-1.0)/fact(k)*pow(x/2.0,2.0*k-n);
    sum2=-sum2/M_PI;
    k=-1.0;
    sum3=0.0;
    s=0.0;
    do{
       k=k+1.0;
       sum3=sum3+s;
       s=pow(-1.0,k)/fact(k)/fact(n+k)*pow(x/2.0,2*k+n);
       ss=0.0;
       for(i=0.0;i<(k+1.0);i=i+1.0)
         {
          sss=0.0;
          for(nn=1.0;nn<(n+i+1.0);nn=nn+1.0)
        sss=sss+1.0/nn;
          ss=ss+sss;
         }
       s=s*ss;
      }while(module(s)>=eps);
    sum3=-sum3/M_PI;
    return(sum1+sum2+sum3);
       }
     if(n==1.0)
       {
    n=1.0;
    sum1=2.0/M_PI*(C+log(x/2.0))*J(x,n);
    sum2=0.0;
    for(k=0.0;k<n;k=k+1.0)
      sum2=sum2+fact(n-k-1.0)/fact(k)*pow(x/2.0,2.0*k-n);
    sum2=-sum2/M_PI;
    k=-1.0;
    sum3=0.0;
    s=0.0;
    do{
       k=k+1.0;
       sum3=sum3+s;
       s=pow(-1.0,k)/fact(k)/fact(n+k)*pow(x/2.0,2*k+n);
       ss=0.0;
       for(i=0.0;i<(k+1.0);i=i+1.0)
         {
          sss=0.0;
          for(nn=1.0;nn<(n+i+1.0);nn=nn+1.0)
        sss=sss+1.0/nn;
          ss=ss+sss;
         }
       s=s*ss;
      }while(module(s)>=eps);
    sum3=-sum3/M_PI;
    return(sum1+sum2+sum3);
       }
    }
  }
 
complex Ne_der(double x,double n)
  { return((Ne(x,n-1.0)-Ne(x,n+1.0))/2.0); }
 
complex Ne_der_der(double x,double n)
  { return((Ne_der(x,n-1.0)-Ne_der(x,n+1.0))/2.0); }
 
complex H1(double x,double n)
  { return(J(x,n)+iii*Ne(x,n)); }
 
complex H1_der(double x,double n)
  { return(J_der(x,n)+iii*Ne_der(x,n)); }
 
complex H1_der_der(double x,double n)
  { return(J_der_der(x,n)+iii*Ne_der_der(x,n)); }
 
void mod_upr(void)
  {
   if(zad==1)
     {
      L1=11.2*pow(10.0,11.0);
      M1=8.1*pow(10.0,11.0);
      R1=7.7;
      L3=5.3*pow(10.0,11.0);
      M3=2.6*pow(10.0,11.0);
      R3=2.7;
      L2=11.2*pow(10.0,11.0);
      M2=8.1*pow(10.0,11.0);
      R2=7.7;
     }
   if(zad==2)
     {
      L1=5.3*pow(10.0,11.0);
      M1=2.6*pow(10.0,11.0);
      R1=2.7;
      L3=11.2*pow(10.0,11.0);
      M3=8.1*pow(10.0,11.0);
      R3=7.7;
      L2=5.3*pow(10.0,11.0);
      M2=2.6*pow(10.0,11.0);
      R2=2.7;
     }
  }
 
double k1,xi1,k2,xi2,k3,xi3;
complex A1_n[K+K+2],A2_n[K+K+2],A3_n[K+K+2],A4_n[K+K+2];
complex B1_n[K+K+2],B2_n[K+K+2],B3_n[K+K+2],B4_n[K+K+2];
 
complex A[MM][MM]; 
complex F[MM]; 
complex X[MM]; 
float a[N][N]; 
float f[N]; 
float x[N]; 
 
void Matrix_A_F(float n)
  {
   A[0][0]=k1*H1_der(k1*r1,n);
   A[0][1]=iii*n/r1*H1(xi1*r1,n);
   A[0][2]=0.0;
   A[0][3]=0.0;
   A[0][4]=-k3*J_der(k3*r1,n);
   A[0][5]=-k3*Ne_der(k3*r1,n);
   A[0][6]=-iii*n/r1*J(xi3*r1,n);
   A[0][7]=-iii*n/r1*Ne(xi3*r1,n);
   F[0]=-pow(iii,n)*k1*J_der(k1*r1,n);
 
   A[1][0]=0.0;
   A[1][1]=0.0;
   A[1][2]=k2*J_der(k2*r2,n);
   A[1][3]=iii*n/r2*J(xi2*r2,n);
   A[1][4]=-k3*J_der(k3*r2,n);
   A[1][5]=-k3*Ne_der(k3*r2,n);
   A[1][6]=-iii*n/r2*J(xi3*r2,n);
   A[1][7]=-iii*n/r2*Ne(xi3*r2,n);
   F[1]=0.0;
 
   A[2][0]=iii*n/r1*H1(k1*r1,n);
   A[2][1]=-xi1*H1_der(xi1*r1,n);
   A[2][2]=0.0;
   A[2][3]=0.0;
   A[2][4]=-iii*n/r1*J(k3*r1,n);
   A[2][5]=-iii*n/r1*Ne(k3*r1,n);
   A[2][6]=xi1*J_der(xi3*r1,n);
   A[2][7]=xi3*Ne_der(xi3*r1,n);
   F[2]=-pow(iii,n+1.0)*n/r1*J(k1*r1,n);
 
 
   A[3][0]=0.0;
   A[3][1]=0.0;
   A[3][2]=iii*n/r2*J(k2*r2,n);
   A[3][3]=-xi2*J_der(xi2*r2,n);
   A[3][4]=-iii*n/r2*J(k3*r2,n);
   A[3][5]=-iii*n/r2*Ne(k3*r2,n);
   A[3][6]=xi3*J_der(xi3*r2,n);
   A[3][7]=xi3*Ne_der(xi3*r2,n);
   F[3]=0.0;
 
   A[4][0]=2.0*M1*k1*k1*H1_der_der(k1*r1,n)-L1*k1*k1*H1(k1*r1,n);
   A[4][1]=2.0*M1*iii*n/r1*(xi1*H1_der(xi1*r1,n)-H1(xi1*r1,n)/r1);
   A[4][2]=0.0;
   A[4][3]=0.0;
   A[4][4]=-2.0*M3*k3*k3*J_der_der(k3*r1,n)-L3*k3*k3*J(k3*r1,n);
   A[4][5]=-2.0*M3*k3*k3*Ne_der_der(k3*r1,n)-L3*k3*k3*Ne(k3*r1,n);
   A[4][6]=-2.0*M3*iii*n/r1*(xi3*J_der(xi3*r1,n)-J(xi3*r1,n)/r1);
   A[4][7]=-2.0*M3*iii*n/r1*(xi3*Ne_der(xi3*r1,n)-Ne(xi3*r1,n)/r1);
   F[4]=-pow(iii,n)*(2.0*M1*k1*k1*J_der_der(k1*r1,n)-L1*k1*k1*J(k1*r1,n));
 
   A[5][0]=2.0*M1*iii*n/r1*(k1*H1_der(k1*r1,n)-H1(k1*r1,n)/r1);
   A[5][1]=M1*(-xi1*xi1*H1_der_der(xi1*r1,n)-n*n/r1/r1*H1(xi1*r1,n)+xi1/r1*H1_der(xi1*r1,n));
   A[5][2]=0.0;
   A[5][3]=0.0;
   A[5][4]=-2.0*M3*iii*n/r1*(k3*J_der(k3*r1,n)-J(k3*r1,n)/r1);
   A[5][5]=-2.0*M3*iii*n/r1*(k3*Ne_der(k3*r1,n)-Ne(k3*r1,n)/r1);
   A[5][6]=-M3*(-xi3*xi3*J_der_der(xi3*r1,n)-n*n/r1/r1*J(xi3*r1,n)+xi3/r1*J_der(xi3*r1,n));
   A[5][7]=-M3*(-xi3*xi3*Ne_der_der(xi3*r1,n)-n*n/r1/r1*Ne(xi3*r1,n)+xi3/r1*Ne_der(xi3*r1,n));
   F[5]=-2.0*M1/r1*pow(iii,n+1.0)*n*(k1*J_der(k1*r1,n)-J(k1*r1,n)/r1);
 
   A[6][0]=0.0;
   A[6][1]=0.0;
   A[6][2]=2.0*M2*k2*k2*J_der_der(k2*r2,n)-L2*k2*k2*H1(k2*r2,n);
   A[6][3]=2.0*M2*iii*n/r2*(xi2*H1_der(xi2*r2,n)-H1(xi2*r2,n)/r2);
   A[6][4]=-2.0*M3*k3*k3*J_der_der(k3*r2,n)-L3*k3*k3*J(k3*r2,n);
   A[6][5]=-2.0*M3*k3*k3*Ne_der_der(k3*r2,n)-L3*k3*k3*Ne(k3*r2,n);
   A[6][6]=-2.0*M3*iii*n/r2*(xi3*J_der(xi3*r2,n)-J(xi3*r2,n)/r2);
   A[6][7]=-2.0*M3*iii*n/r2*(xi3*Ne_der(xi3*r2,n)-Ne(xi3*r2,n)/r2);
   F[6]=0.0;
 
   A[7][0]=0.0;
   A[7][1]=0.0;
   A[7][2]=2.0*M2*iii*n/r2*(k2*H1_der(k2*r2,n)-H1(k2*r2,n)/r2);
   A[7][3]=M2*(-xi2*xi2*H1_der_der(xi2*r2,n)-n*n/r2/r2*H1(xi2*r2,n)+xi2/r2*H1_der(xi2*r2,n));
   A[7][4]=-2.0*M3*iii*n/r2*(k3*J_der(k3*r2,n)-J(k3*r2,n)/r2);
   A[7][5]=-2.0*M3*iii*n/r2*(k3*Ne_der(k3*r2,n)-Ne(k3*r2,n)/r2);
   A[7][6]=-M3*(-xi3*xi3*J_der_der(xi3*r2,n)-n*n/r2/r2*J(xi3*r2,n)+xi3/r2*J_der(xi3*r2,n));
   A[7][7]=-M3*(-xi3*xi3*Ne_der_der(xi3*r2,n)-n*n/r2/r2*Ne(xi3*r2,n)+xi3/r2*Ne_der(xi3*r2,n));
   F[7]=0.0;
  }
 
void Real_Gauss(void)
  {
   int i,j,k,l,maxk;
   float max,w[N],v[N][N],sum,e,c;
 
   for(i=0;i<N;i++)
     {
      for(j=0;j<N;j++)
    v[i][j]=a[i][j];
      w[i]=f[i];
     }
   for(k=0;k<N;k++)
     {
      maxk=k;
      max=fabs(a[k][k]);
      for(i=k;i<N;i++)
    if(fabs(a[i][k])>max)
       {
        maxk=i;
        max=fabs(a[i][k]);
       }
      for(i=0;i<N;i++)
    {
     e=a[k][i];
     a[k][i]=a[maxk][i];
     a[maxk][i]=e;
    }
      e=f[k];
      f[k]=f[maxk];
      f[maxk]=e;
      for(i=k+1;i<N;i++)
    {
     c=a[i][k]/a[k][k];
     f[i]=f[i]-f[k]*c;
     for(j=k;j<N;j++)
       a[i][j]=a[i][j]-a[k][j]*c;
    }
     }
   for(i=0;i<N;i++)
     x[i]=0.0;
   for(i=N-1;i>=0;i--)
     {
      c=0.0;
      for(j=i+1;j<N;j++)
    c=c+a[i][j]*x[j];
      x[i]=(f[i]-c)/a[i][i];
     }
  }
 
 
void Complex_Gauss(void)
 {
  int i,j;
  complex sum;
 
  for(i=0;i<MM;i++)
   {
     a[2*i][0]=real(A[i][0]);  a[2*i][1]=-imag(A[i][0]);
     a[2*i][2]=real(A[i][1]);  a[2*i][3]=-imag(A[i][1]);
     a[2*i][4]=real(A[i][2]);  a[2*i][5]=-imag(A[i][2]);
     a[2*i][6]=real(A[i][3]);  a[2*i][7]=-imag(A[i][3]);
     a[2*i][8]=real(A[i][4]);  a[2*i][9]=-imag(A[i][4]);
    a[2*i][10]=real(A[i][5]); a[2*i][11]=-imag(A[i][5]);
    a[2*i][12]=real(A[i][6]); a[2*i][13]=-imag(A[i][6]);
    a[2*i][14]=real(A[i][7]); a[2*i][15]=-imag(A[i][7]);
    f[2*i]=real(F[i]);
     a[2*i+1][0]=imag(A[i][0]);  a[2*i+1][1]=real(A[i][0]);
     a[2*i+1][2]=imag(A[i][1]);  a[2*i+1][3]=real(A[i][1]);
     a[2*i+1][4]=imag(A[i][2]);  a[2*i+1][5]=real(A[i][2]);
     a[2*i+1][6]=imag(A[i][3]);  a[2*i+1][7]=real(A[i][3]);
     a[2*i+1][8]=imag(A[i][4]);  a[2*i+1][9]=real(A[i][4]);
    a[2*i+1][10]=imag(A[i][5]); a[2*i+1][11]=real(A[i][5]);
    a[2*i+1][12]=imag(A[i][6]); a[2*i+1][13]=real(A[i][6]);
    a[2*i+1][14]=imag(A[i][7]); a[2*i+1][15]=real(A[i][7]);
    f[2*i+1]=imag(F[i]);
   }
   Real_Gauss();
   X[0]=complex(x[0],x[1]);
   X[1]=complex(x[2],x[3]);
   X[2]=complex(x[4],x[5]);
   X[3]=complex(x[6],x[7]);
   X[4]=complex(x[8],x[9]);
   X[5]=complex(x[10],x[11]);
   X[6]=complex(x[12],x[13]);
   X[7]=complex(x[14],x[15]);
  }
 
void grafic(double *k_1, double *k_2, double *k_3, double *k_4, double a, double b, double c, double d, double col_x, double col_y)
  {
   double h,hx,hy,dx,dy;
   int i,maxx,maxy;
   int borderx_left=0,borderx_right=0;
   int bordery_up=0,bordery_down=0;
   int gdriver=DETECT, gmode, errorcode;
 
   clrscr();
   initgraph(&gdriver,&gmode," ");
   errorcode=graphresult();
   if(errorcode!=grOk)
     {
      printf("Не могу открыть графический экран!\n");
      printf("Нажмите любую клавишу!\n");
      getch();
      exit(1);
     }
   setfillstyle(SOLID_FILL,WHITE);
   floodfill(0,0,WHITE);
   maxx=getmaxx();
   maxy=getmaxy();
   h=(double)(maxx-(borderx_left+borderx_right));
   hx=(b-a)/h;
   h=(double)(maxy-(bordery_up+bordery_down));
   hy=(d-c)/h;
   setcolor(BLACK);
   line(borderx_left,bordery_up,borderx_left,maxy-bordery_down);
   line(borderx_left,maxy-bordery_down,maxx-borderx_right,maxy-bordery_down);
   line(maxx-borderx_right,maxy-bordery_down,maxx-borderx_right,bordery_up);
   line(maxx-borderx_right,bordery_up,borderx_left,bordery_up);
   line(0,0,0,maxy);
   line(0,maxy,maxx,maxy);
   line(maxx,maxy,maxx,0);
   line(maxx,0,0,0);
   dx=(b-a)/col_x;
   dy=(d-c)/col_y;
   setcolor(BLACK);
   for(i=1;i<col_x;i++)
     line(borderx_left+i*dx/hx,bordery_up,borderx_left+i*dx/hx,maxy-bordery_down);
   for(i=1;i<col_y;i++)
     line(borderx_left,bordery_up+i*dy/hy,maxx-borderx_right,bordery_up+i*dy/hy);
   setcolor(BLACK);
   for(i=0;i<M;i++)
     line(borderx_left+(k_1[i]-a)/hx, maxy-bordery_down-(k_2[i]-c)/hy,
      borderx_left+(k_1[i+1]-a)/hx, maxy-bordery_down-(k_2[i+1]-c)/hy);
   setcolor(BLACK);
   for(i=0;i<M;i++)
     line(borderx_left+(k_3[i]-a)/hx, maxy-bordery_down-(k_4[i]-c)/hy,
      borderx_left+(k_3[i+1]-a)/hx, maxy-bordery_down-(k_4[i+1]-c)/hy);
   getch();
   closegraph();
  }
 
double F_rass(double fi)
  {
   complex sum;
   int i;
 
   sum=0.0;
   for(i=-K;i<=K;i=i+1.0)
      sum=sum+pow(iii,i)*A1_n[K-i]*exp(iii*i*fi);
   sum=2.0/sqrt(M_PI*const_w)*module(sum);
   return(module(sum));
  }
 
void main(void)
  {
   int j;
   double k,n;
   double k_0,k_n,dk;
   double k_1[M+1],k_2[M+1],k_3[M+1],k_4[M+1];
 
   clrscr();
   const_w=2.0;
   r1=3.5;
   r2=1.0;
   for(j=0;j<(M+1);j++)
     {
      k_1[j]=0.0;
      k_2[j]=0.0;
      k_3[j]=0.0;
      k_4[j]=0.0;
     }
   clrscr();
   k_0=M_PI;
   k_n=2.0*M_PI;
   dk=(k_n-k_0)/M;
   j=0;
   zad=1;
   mod_upr();
   w=module(const_w*sqrt((L1+2.0*M1)/R1)/(r1-r2));
   k1=w/sqrt((L1+2.0*M1)/R1);
   k2=w/sqrt((L2+2.0*M2)/R2);
   k3=w/sqrt((L3+2.0*M3)/R3);
   xi1=w/sqrt(M1/R1);
   xi2=w/sqrt(M2/R2);
   xi3=w/sqrt(M3/R3);
   for(n=-K;n<=K;n=n+1)
     {
      Matrix_A_F(n);
      Complex_Gauss();
      A1_n[K-n]=X[0];
      B1_n[K-n]=X[1];
      A2_n[K-n]=X[2];
      B2_n[K-n]=X[3];
      A3_n[K-n]=X[4];
      A4_n[K-n]=X[5];
      B3_n[K-n]=X[6];
      B4_n[K-n]=X[7];
     }
   for(j=0,k=k_0;k<=k_n;k=k+dk,j++)
     {
      k_1[j]=-F_rass(k)*cos(k);
      k_2[j]=-F_rass(k)*sin(k);
      k_3[j]=-F_rass(k)*cos(k);
      k_4[j]=F_rass(k)*sin(k);
     }
   grafic(k_1,k_2,k_3,k_4,-2.0,6.0,-2.0,2.0,8.0,4.0);
  }