| 查看: 604 | 回复: 1 | ||
[求助]
请教兰姆波频散曲线问题
|
|
学过固体的超声波的高手，帮帮忙，你们谁有兰姆波频散曲线计算的程序，能够传给我看看吗？很着急。dzxukai@sina.cn 谢谢！ 非高级版内容，请勿放入高级版中 [ Last edited by 华丽的飘过 on 2012-9-22 at 21:55 ] |
回复此楼» 猜你喜欢
回收溶剂求助
已经有7人回复
-
职称评审没过，求安慰
已经有40人回复
-
硝基苯如何除去
已经有3人回复
-
A期刊撤稿
已经有4人回复
-
垃圾破二本职称评审标准
已经有17人回复
-
投稿Elsevier的Neoplasia杂志，到最后选publishing options时页面空白，不能完成投稿
已经有22人回复
-
申请26博士
已经有5人回复
-
EST投稿状态问题
已经有7人回复
-
毕业后当辅导员了，天天各种学生超烦
已经有4人回复
-
求助文献
已经有3人回复
wilde2540
木虫 (正式写手)
- 应助: 3 (幼儿园)
- 金币: 3736.8
- 散金: 7
- 红花: 4
- 帖子: 590
- 在线: 188.9小时
- 虫号: 835104
- 注册: 2009-08-27
- 专业: 声学
【答案】应助回帖
|
#include #include #include #include #include #include #include #define PI 3.1415926 // double cd=5.8e3,ct=3.1e3; //钢 // double cd=6.42e3,ct=3.04e3; //铝 double cd=2.68e3,ct=1.1e3; //有机玻璃 int BisectRoot(double a,double b,double h,double eps,double fb,double* x,int n,int *m); double Func(double cp,double fb); void main() { int i,j,n=20,m,N=300; double a,b,h,eps,*x,*cp,*dcp,*cg,fb,D=1e-3,step=0.02e3; FILE *fp; char str[10],str1[10],cgfile[50][10]; printf("Input the cg file name,please:\n"  ;scanf("%s",str1); for(j=0;j sprintf(str,"%d.dat",j); strcpy(cgfile[j],str1); strcat(cgfile[j],str); } a=1.0e0; b=15.0e3; h=0.7; eps=1e-6; for (j=0;j<10;j++) { if((fp=fopen(cgfile[j],"w" )==NULL){ printf("cg file %s can't open.\n",cgfile[j]); exit(0);} else printf("file %s is written.\n",cgfile[j]); cp=(double*)calloc(N,sizeof(double)); if(cp==NULL) exit(1); dcp=(double*)calloc(N,sizeof(double)); if(dcp==NULL) exit(1); cg=(double*)calloc(N,sizeof(double)); if(cg==NULL) exit(1); m=0; i=0; fb=0.0; do { x=(double*)calloc(n,sizeof(double)); if(x==NULL) exit(1); BisectRoot(a,b,h,eps,fb,x,n,&m); cp=x[j]; if(cp>0 && i>2) { dcp=(cp-cp[i-1])/step; cg=cp/(1.0-fb*dcp/cp); fprintf(fp,"%6.3f %6.3f\n ",2*(fb-step)*D, cg*D); //fb*D is frequency with unit kHz // c[j] is velocity with unit m/s // and c[j]*D is velocity with unit km/s } i++; fb+=step; // fb is frequency with unit Hz, step=50 free(x); }while(i printf("m=%d\n\n",m); fclose(fp); free(cp);free(dcp);free(cg); } } double Func(double cp,double fb) { double ab,bt,kk,at,kd,kt; double fs,fa; ab=sqrt(fabs(cp*cp/(cd*cd)-1.0)); bt=sqrt(fabs(cp*cp/(ct*ct)-1.0)); kk=(cp*cp/(ct*ct)-2.0)*(cp*cp/(ct*ct)-2.0); at=4*ab*bt; kd=ab*PI*2*fb/cp; kt=bt*PI*2*fb/cp; if(cp>=cd) { fs=kk*sin(kt)*cos(kd)+at*sin(kd)*cos(kt); fa=kk*sin(kd)*cos(kt)+at*sin(kt)*cos(kd); } if(cp>ct && cp fs=kk*sin(kt)*cosh(kd)-at*sinh(kd)*cos(kt); fa=kk*sinh(kd)*cos(kt)+at*sin(kt)*cosh(kd); } if(cp fs=kk*tanh(kt)-at*tanh(kd); fa=kk*tanh(kd)-at*tanh(kt); } // return (fs); // 对称模式 return (fa); // 反对称模式 } int BisectRoot(double a,double b,double h,double eps,double fb,double* x,int n,int *m) { double z,z0,z1,y,y0,y1; *m=0; z=a; y=Func(z,fb); while(1) { if((z>b+h/2)||(*m==n)) return(1); if(fabs(y) *m+=1; x[*m-1]=z; z+=h/2; y=Func(z,fb); continue; } z1=z+h; y1=Func(z1,fb); if(fabs(y1) *m+=1; x[*m-1]=z1; z=z1+h/2; y=Func(z,fb); continue; } if(y*y1>0) { y=y1; z=z1; continue; } while(1) { if(fabs(z1-z) *m+=1; x[*m-1]=(z1+z)/2; z=z1+h/2; y=Func(z,fb); break; } z0=(z1+z)/2; y0=Func(z0,fb); if(fabs(y0) *m+=1; x[*m-1]=z0; z=z0+h/2; y=Func(z,fb); break; } if(y*y0<0) { z1=z0; y1=y0; } else { z=z0; y=y0; } } } } N年之前的代码，不负责对程序进行解释，仅供参考喔。图片也是N之前计算的一个例子  |
2楼2011-10-31 15:21:11