Курсовая работа: Моделирование рассеяния плоской упругой продольной волны на упругом однородном изотропном цилиндрическом слое

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]);

Страницы: 1, 2, 3, 4, 5