联剖装置球体 异常的计算5.docVIP

课程设计报告 ——联剖装置球体异常的计算 专业：地球物理学 姓名： 学号： 指导教师： 二、实验C++程序： #includestdio.h8 #includemath.h #define pi 3.1416 double fun(double n,double x) { double Pn[1000]; Pn[0]=1;Pn[1]=x; if((int)n=2) Pn[(int)n]=(2*n-1)/n*Pn[(int)n-1]-(n-1)/n*Pn[(int)n-2];//Pn(y)计算的递推公式,3—3公式// return Pn[(int)n]; } int main() { FILE *fp; int i; double k,x,n1,n2,n3,n4,t1,t2,t3,t4,temp1,temp2,temp3,temp4,temp5,temp6,temp7,temp8,da,db,rm,rn,Ao,oB,Mo,No, AM=1.5,AN=2.5,MN=1.0,r0=1.0,h=1.5,I=1.0,p1=1.0,p2=0.005,pa[51],pb[51],UMa[51],UMb[51],UNa[51],UNb[51],UMNa[51],UMNb[51]; fp=fopen(data.txt,w+);//定义每个参数的初值// k=2*pi*AM*AN/MN;//计算装置系数// for(i=0,x=-25.0;x26.0;x++,i++)//定义测线长度及步长// { t1=0.0;n1=0.0;t2=0.0;n2=0.0;t3=0.0;n3=0.0;t4=0.0;n4=0.0; No=fabs(x+MN/2);Mo=fabs(x-MN/2); if(x0.0) {Ao=AN+No;oB=AM-No;} if(x=0.0) {Ao=AN-No;oB=AM+No;}//求测线上供电电极、测量电极和球体中心投影的关系// da=pow(pow(Ao,2)+pow(h,2),0.5); db=pow(pow(oB,2)+pow(h,2),0.5); rm=pow(pow(Mo,2)+pow(h,2),0.5);//求低阻体到供电电极和测量电极的距离// rn=pow(pow(No,2)+pow(h,2),0.5);//printf(%lf\n,fabs(-5.0))// do { temp1=(p2-p1)*n1/(p1*n1+p2*(n1+1))*pow(r0,2*n1+1)/(pow(da,n1+1)*pow(rm,n1+1))*fun(n1,x); temp2=(p2-p1)*(n1+1)/(p1*(n1+1)+p2*((n1+1)+1))*pow(r0,2*(n1+1)+1)/(pow(da,(n1+1)+1)*pow(rm,(n1+1)+1))*fun(n1+1,x); t1+=temp1; n1++; } while((temp1-temp2)0.01);//由子程序求电位U，第二项开始进行累加，n值以第n项和第n+1项之差小于1%// do { temp3=(p2-p1)*n2/(p1*n2+p2*(n2+1))*pow(r0,2*n2+1)/(pow(db,n2+1)*pow(rm,n2+1))*fun(n2,x); temp4=(p2-p1)*(n2+1)/(p1*(n2+1)+p2*((n2+1)+1))*pow(r0,2*(n2+1)+1)/(pow(db,(n2+1)+1)*pow(rm,(n2+1)+1))*fun(n2+1,x); t2+=temp3; n2++; } while((temp3-temp4)0.01); do { temp5=(p2-p1)*n3/(p1*n3+p2*(n3+1))*pow(r0,2*n3+1)/(pow(da,n3+1)*pow(rn,n3+1))*fun(n3,x); temp6=(p2-p1)*(n3+1)/(p1*(n3+1)+p2*((n3+1)+1))*pow(r0,2*(n3+1)+1)/(pow(da,(n3+1)+1)*pow(rn,(n3+1)+1))*fun(n3+1,x); t3+=temp5; n3++; } while((temp5-temp6)0.01); do { temp7=(p2-p1)*n4/(p1*n4+p2*(n4+1))*pow(r0,2*n4+1)/(pow(db,n4+1)*pow(rn,n4+1))*