地球物理 课件设计2.docxVIP

《课程设计》程序及成果图件报告 姓名： 学号： 专业：地球物理学 指导老师： 2014.1.7 前言： 重磁异常正演是在给定地质体的形状大小、空间位置及物性参数的条件下，求它在它外部空间任意点上的物理场值。 电子计算机的应用使正演方法进入了一个新的阶段，它可以对任意形状物体进行计算；对物性“非均匀”的地质体，也可以通过将其分解成若干个“均匀”物性的物体的方法和分区块均匀的数值方法来计算，使正演方法可以达到相当高的精度。 以下是我对截面为多边形的水平柱体模型的重磁异常的正演程序。 C语言代码： #include stdio.h #include math.h #define pi 3.14159 #define n 101 /**********测点数********/ #define N 360 /*********拟合多边形边数********/ double fun(double a,double b) /********关于atan(b/a)的讨论********/ { doublethet; if(fabs(a)1.0e-15) { if((b/a)0.0) { if(a0.0b0.0) thet=atan(b/a); else if(a0.0b0.0) thet=atan(b/a)-pi; } else if((b/a)0.0) thet=atan(b/a)+pi; } else if(fabs(a)=1.0e-15) { if(b0.0) thet=-pi/2; else if(b=0.0) thet=pi/2; } return(thet); } void main() { FILE *fp; fp=fopen(输出.txt,w+); double x[N]={0.0},z[N]={0.0},P[n]={0.0},Q[n]={0.0}; deltag[n]={0.0},deltax[n]={0.0},deltaz[n]={0.0}; double M=2000.0,G=6.67e-5,sigma=2.67,r=200.0,x0=1000.0,z0=1000.0; inti,k; doublexk[101],x1[N],z1[N]; double g[n]={0.0},alpha=(2*pi)/N; doubleMx=M*sin(pi/4),Mz=M*cos(pi/4); for(i=0;iN;i++) { x1[i]=x0+r*sin(alpha*i); z1[i]=z0-r*cos(alpha*i); } for(k=0;kn;k++) { xk[k]=20.0*k; for(i=0;iN-1;i++) { x[i]=x1[i]-xk[k]; x[i+1]=x1[i+1]-xk[k]; z[i]=z1[i]; z[i+1]=z1[i+1]; Q[k]+=((z[i+1]-z[i])*(x[i]-x[i+1])/(pow((z[i+1]-z[i]),2)+pow((x[i]-x[i+1]),2))) *(fun(x[i],z[i])-fun(x[i+1],z[i+1]))-(1/2)*pow((z[i+1]-z[i]),2)*log((pow(x[i+1],2)+pow(z[i+1],2))/(pow(x[i],2)+pow(z[i],2)))/(pow((z[i+1]-z[i]),2)+pow((x[i]-x[i+1]),2)); P[k]+=(pow((z[i+1]-z[i]),2)/(pow((z[i+1]-z[i]),2)+pow((x[i]-x[i+1]),2))) *(fun(x[i],z[i])-fun(x[i+1],z[i+1]))-(1/2)*(z[i+1]-z[i])*(x[i]-x[i+1])*log((pow(x[i+1],2)+pow(z[i+1],2))/(pow(x[i],2)+pow(z[i],2)))/(pow((z[i+1]-z[i]),2)+pow((x[i]-x[i+1]),2)); deltax[k]=(Mx*P[k]+Mz*Q[k])/(2*pi); deltaz[k]=(Mx*Q[k]-Mz*P[k])/(2*