地球物理 课程设计1.docx

【课程设计】联剖装置球体s异常计算姓名： 学号： 专业： 地球物理学指导老师： 一 要求：计算点场源中有高阻球体存在时，主剖面上联合剖面装置sA、sB的分布。如下图所示：二 数学模型先求电位：UMA、UMB、UNA、UNB.计算公式如下：其中，各参数之间的关系如下：其次，勒让德多项式采用递推公式：Pn(y)=2n-1nPn-1(y)- n-1nPn-2(y);且P0(y)=1，P1(y)=y， y=cos；装置系数 K=2AM*ANMN。初始值AM=1.5cm, AN=2.5, MN=1.0。半径r0=1.0m, 埋深h=1.5m, 供电电流密度I=1.0A。 球体电阻率2=0.005m,围岩电阻率1=1.0m。三 C语言程序代码：#include#include#define pi 3.1415926double fun(double n,double y)double Pn1000;Pn0=1;Pn1=y;if(int)n=2)Pn(int)n=(2*n-1)/n*Pn(int)n-1-(n-1)/n*Pn(int)n-2; /Pn(y)计算的递推公式,33公式/ return Pn(int)n;void main()FILE *fp;int i;double cos1,cos2,cos3,cos4,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,pa51,pb51,UMa51,UMb51,UNa51,UNb51,UMNa51,UMNb51; 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(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); cos1=(pow(da,2)+pow(rm,2)-pow(AM,2)/2*da*rm;cos2=(pow(da,2)+pow(rn,2)-pow(AN,2)/2*da*rn;cos3=(pow(db,2)+pow(rm,2)-pow(AN,2)/2*da*rm;cos4=(pow(db,2)+pow(rn,2)-pow(AM,2)/2*db*rn;dotemp1=(p2-p1)*n1/(p1*n1+p2*(n1+1)*pow(r0,2*n1+1)/(pow(da,n1+1)*pow(rm,n1+1)*fun(n1,cos1);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,cos1); t1+=temp1; n1+; while(temp1-temp2)0.01);/由子程序求电位U，第二项开始进行累加，n值以第n项和第n+1项之差小于1%/ dotemp3=(p2-p1)*n2/(p1*n2+p2*(n2+1)*pow(r0,2*n2+1)/(pow(db,n2+1)*pow(rm,n2+1)*fun(n2,cos2); 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,cos2); t2+=temp3;n2+; while(temp3-temp4)0.01); dotemp5=(p2-p1)*n3/(p1*n3+p2*(n3+1)*pow(r0,2*n3+1)/(pow(da,n3+1)*pow(rn,n3+1)*fun(n3,cos3); 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,cos3); 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)*fun(n4,cos4); temp8=(p2-p1)*(n4+1)/(p1*(n4+1)+p2*(n4+1)+1)*pow(r0,2*(n4+1)+1)/(pow(db,(n4+1)+1)*pow(rn,(n4+1)+1)*fun(n4+1,cos4); t4+=temp7;n4+; while(temp7-temp8)0.01);UMai=I*p1/(2*pi)*(1.0/AM+2*t1); UMbi=-I*p1/(2*pi)*(1.0/AN+2*t2); UNai=I*p1/(2*pi)*(1.0/AN+2*t3); UNbi=-I*p1/(2*pi)*(1.0/AM+2*t4); /用x循环分别求四个电位分量/ UMNai=UMai-UNai; UMNbi=UMbi-UNbi;/求两