#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define PI 3.1415926
#define XN 101
#define ZN 101
#define DH 5
#define KN 200
#define DT 0.001
#define FM 30
#define R 3
int main(void)
{
int m,n,k,i,j;
float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],u5[XN][ZN],u6[XN][ZN],uu0,uu1,uu2,f[XN][ZN];
float s[KN],v[XN][ZN];
FILE *fp;
for(k=0;k<KN;k++)
s[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);
for(i=0;i<XN;i++)
{
for(j=0;j<ZN;j++)
{
if(i==50&&j==50)
f[i][j]=1;
else
f[i][j]=0;
}
}
for(i=0;i<XN;i++)
{
for(j=0;j<ZN;j++)
{
u1[i][j]=0;
u2[i][j]=0;
u3[i][j]=0;
u4[i][j]=0;
u6[i][j]=0;
}
}
for(i=0;i<XN;i++)
{
for(j=0;j<ZN/4;j++)
v[i][j]=2400;
for(j=ZN/4;j<ZN;j++)
v[i][j]=2400;
}
for(k=0;k<KN;k++)
{
for(i=2;i<XN-2;i++)
{
for(j=2;j<ZN-2;j++)
{
uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);
uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];
uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];
u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+s[k]*f[i][j];
u5[i][k]=u3[i][10];
}
}
for(m=0;m<XN;m++)
{
for(n=0;n<ZN;n++)
{
u1[m][n]=u2[m][n];
u2[m][n]=u3[m][n];
}
}
if(k==100)
for(m=0;m<XN;m++)
for(n=0;n<ZN;n++)
u4[m][n]=u3[m][n];//记录波前快照,中间点
}
if(i=40)
{
for(n=0;n<ZN;n++)
u6[i][n]=u3[i][n];
}
if((fp=fopen("wavefront.dat","w"))!=NULL)
{
fprintf(fp,"%d\n",XN);
fprintf(fp,"%d\n",ZN);
for(i=0;i<ZN;i++)
{
for(j=0;j<ZN;j++)
fprintf(fp,"%f",u4[i][j]);
fprintf(fp,"\n");
}
fclose(fp);
}
if((fp=fopen("record.dat","w"))!=NULL)
{fprintf(fp,"%d\n",XN/2);
fprintf(fp,"%d\n",KN);
for(i=0;i<XN;i++)
for(j=0;j<KN;j++)
fprintf(fp,"%f\n",u5[i][j]);
fclose(fp);
}
return 0;
}
评论0