地球物理课程设计.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*pi);g[k]+=(1/2)*(z[i+1]-z[i])*(x[i]*z[i+1]-x[i+1]*z[i])*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))-((x[i+1]-x[i])*(x[i]*z[i+1]-x[i+1]*z[i])/(pow(