三度体球体磁场正演

一 、 试验要求：

根据正演公式，计算给定参数球体磁场水平X分量、水平Y分量，垂直z分量，总磁异常，并总结异常特征。

二 、实验课时：2课时

三 、实验成果提交：

各磁场分量及总场平面等值线图、主剖面图。

四 、球体正演公式

五 、球体正演参数

六 、球体正演图(平面等值线图和主剖面图)

1. 水平X分量图
   1. 平面等值线图

1. 主剖面

1. 水平Y分量
   1. 平面等值线图

   

   1. 主剖面

1. 垂直Z分量图
   1. 平面等值线图
   2. 主剖面

1. 总场异常图
   1. 平面等值线图

1. 主剖面图

异常特征分析：Hax,Hay,Za,主剖面过正负异常中心。 Za等值线呈现等轴状，负异常包着正异常：极大值和极小值的连线对应磁化强度M在平面上的投影方向：极小值位于正异常的北侧，极大值位于坐标原点的南侧，用Za平面等值线图中的极大值和极小值的连线确定主剖面。由图可知，受磁化倾角的影响比Za大，在相同的磁化倾角其负值较大。

七、源代码

#include<stdio.h>

#include<math.h>

#include<stdlib.h>

#define dz1 "file1.dat"

main()

{

double m,M=1;int x,y,z=0,R=200,R0=100;float I=0.785,A=0.785,Ha,Hb,Za,T;int i,j;double s1,s2,s3;

FILE *fp;

fp=fopen(dz1,"w+");

fclose(fp);

fp=fopen(dz1,"a+");//文件操作

m=3.14*R0*R0*R0*4/3;

for(x=-500;x<=500;x=x+10)

{

for(y=-500;y<=500;y=y+10)//计算开始

{

Ha=100*m/(pow((x*x+y*y+R*R),(2.5)))*((2*x*x-y*y-R*R)*cos(I)*cos(A)-3*R*x*sin(I)+3*x*y*cos(I)*sin(A));

Hb=100*m/(pow((x*x+y*y+R*R),(2.5)))*((2*y*y-x*x-R*R)*cos(I)*sin(A)-3*R*y*sin(I)+3*x*y*cos(I)*cos(A));

Za=100*m/(pow((x*x+y*y+R*R),(2.5)))*((2*R*R-y*y-x*x)*sin(I)-3*R*x*cos(I)*cos(A)-3*R*y*cos(I)*sin(A));

T =100*m/(pow((x*x+y*y+R*R),(2.5)))*((2*R*R-y*y-x*x)*sin(I)*sin(I)+(2*x*x-y*y-R*R)*cos(I)*cos(I)*cos(A)*

cos(A)+(2*y*y-x*x-R*R)*cos(I)*cos(I)*sin(A)*sin(A)-3*R*x*sin(2*I)*cos(A)+3*x*y*cos(I)*cos(I)*cos(I)*sin(2*A)

-3*R*y*sin(2*I)*sin(A) );//计算

printf("%f\t%f\t%f\t%f\n",Ha,Hb,Za,T);//输出到屏幕

fprintf(fp,"%d\t%d\t%f\t%f\t%f\t%f\n",x,y,Ha,Hb,Za,T);//输出到文件

}

}

fclose(fp);

}