长安大学地磁学磁化曲线报告.docx

一、 实验要求1、 编程计算并绘制倾斜磁化的球体Za与AT等值线平面图2、 编程计算并绘制垂直磁化的球体Za与AT等值线平面图3、 编程计算并绘制不同磁化强度倾角的Za、Ha剖面图4、 编程计算并绘制倾斜磁化的球体Za与AT空间等值线图5、 分析图件，总结收获二、 实验原理在本次试验中，我们用磁化强度倾角或有效磁化强度倾角、剖面与磁化强度水平投影夹角表达磁场Za、Hax和Hay对AT, 一般讨论磁化强度方向与地磁场方向一直情况下的表 达，即磁化强度倾角与地磁场倾角一致，均以I,剖面磁方位角以A表示当剩磁与感磁不 一致时，则地磁倾角用I0表示，剖面磁方位角 用A0表示；磁化强度倾角和剖面与磁化强度 水平投影夹角分别用I，A表示1）球体的磁场表达式在自然界中，当地质体埋深远大于其直径 时，可以近似为球体进行研究如右图，设球体埋深为R，磁化强度为M，体积为v，磁矩m=Mv；球心坐标为Q （0，0，R），设空间任意一点的坐标为P （x，y，z），通过泊松公式或此偶极子磁场公式可以得出其 磁场表达式，由此我们可以将球体的引力位带入泊松公式得到均匀磁化球体磁场G8v球体的引力位为:v=fc^一d对上式求二次导数后，令Z=0，，=R，磁化强度倾角为I，剖面与磁化强度水平投影夹 角为A。

可以求得其个方向磁化强度分量H 二("牝 \2+y2+R22)•〔2 工2- y 2 - Rjcos I cos A 一 3Rx sin I + 3 xy cos I sin AHa^y)• [2 y 2 - X2 一 R2) cos I cos A - 3Ry sin I + 3xy cos I sin A八三 12 + y 2+R 22)• (2 R2 一 y 2 - X2) sin I - 3Rx cos I cos A - 3Ry cos I sin A以及AT的表达式2R2 - X2 - y2) simi + (2X2 - y2 - R2)cos2i cos2 a + \ 2 R2 - X2 - y2) sim acos2i一 3Rx sin 21 cos A + 3xy cos 21 sin 2 A - 3Ry sin 21 sin A]（2）作图参数设置绘制倾斜磁化等值线图时，I=A=45°绘制垂直极化等值线图时1=90 °三、计算程序(1) 等值线平面图计算程序 program ballreal(8)I,A,M,k,T,v,x,y,z,R,u0,Za,Hax,Hay,DT,pi,ce,Z e,Hsz(4),Zsz(4),b,c!I磁化强度倾角，A为剖面与磁化强度水平 投影夹角，M为磁化强度，v为球体体积， x, y，z为P点坐标，R为球体中心埋藏深 度，u0为真空磁导率,Ze为充零数值 open(1,file='input.txt')open(2,file='output.dat') open(3,file='output2.dat') read(1,*)k,T,vA=45;I=90;R=15print*,m,vpi=3.14159u0=4.0*pi*10.0**(-7.0)print*,u0M=k*(T/u0)!输入A，I，RZe=0Za=0DT=01 format(f10.3,'\t'c,f10.3,'\t'c,f25.10)2 format(f10.3,'\t'c,f25.10)do x=-80,80,1do y=-80,80,1ce=4*pi*(x**2+y**2+r**2)**2.5Za=(u0*m*v/ce)*((2.0*R**2.0-x**2.0-y**2.0 )*sind(I)-3*R*x*cosd(I)*cosd(A)-3*R*y*cos d(I)*sind(A))Hay=(u0*m*v)/(ce)*((2*y**2-x**2-R**2)*c osd(I)*sind(A)-3*R*y*sind(I)+3*x*y*cosd(I) *cosd(A))DT=(u0*m*v)/(ce)*((2*R**2-x**2-y**2)*(si nd(I)**2)+(2*x**2-y**2-R**2)*cosd(I)**2*c osd(A)**2+(2*y**2-x**2-R**2)*Cosd(I)**2 *sind(A)**2-3*x*R*sind(2*I)*cosd(A)+3*x* y*cosd(I)**2*sind(2*A)-3*y*R*sind(2*I)*sin d(A))write(2,1)x,y,Za write(3,1)x,y,DT Za=0DT=0enddoHax=0Za=0 enddoend(2) 剖面图计算程序program ballreal(8)I,A,M,k,T,v,x,y,z,R,u0,Za,Hax,Hay,DT,pi,ce,Z e,Hsz(4),Zsz(4),b,c!I磁化强度倾角，A为剖面与磁化强度水平 投影夹角，M为磁化强度，v为球体体积， x, y，z为P点坐标，R为球体中心埋藏深 度，u0为真空磁导率,Ze为充零数值 open(1,file='input.txt')open(2,file='output.dat') open(3,file='output2.dat') open(4,file='output3.dat') open(5,file='output4.dat') open(6,file='output5.dat')open(7,file='output6.dat') open(8,file='output7.dat') open(9,file='output8.dat') open(10,file='output9.dat') open(11,file='output10.dat')1 format(f10.3,'\t'c,f10.3,'\t'c,f25.10)2 format(f10.3,'\t'c,f25.10)read(1,*)k,T,vA=45R=15print*,m,vpi=3.14159(3) 空间等值线图计算程序 program ballreal(8)I,A,M,k,T,v,x,y,z,R,u0,Za,Hax,Hay,DT,pi,ce,Z e,Hsz(4),Zsz(4),b,c!I磁化强度倾角，A为剖面与磁化强度水平 投影夹角，M为磁化强度，v为球体体积， x，y，z为P点坐标，R为球体中心埋藏深 度，u0为真空磁导率,Ze为充零数值 open(1,file='input.txt')open(2,file='output.dat') open(3,file='output2.dat')1 format(f10.3,'\t'c,f10.3,'\t'c,f25.10)2 format(f10.3,'\t'c,f25.10)u0=4.0*pi*10.0**(-7.0)print*,u0M=k*(T/u0)hsz=0zsz=0do x=-100,100,1b=1do I=0,90,30Hsz(b)=(u0*m*v)/(4*pi*(x**2+r**2)**2.5)*((2.0*x**2-R**2)*cosd(I)-3*R*x*sind(I))Zsz(b)=(u0*m*v)/(4*pi*(x**2+r**2)**2.5)*((2.0*R**2.0-x**2.0)*sind(I)-3*R*x*cosd(I))b=b+1enddowrite(4,2)x,Hsz(1)write(5,2)x,Hsz(2)write(6,2)x,Hsz(3)write(7,2)x,Hsz(4)write(8,2)x,Zsz(1)write(9,2)x,Zsz(2)write(10,2)x,Zsz(3)write(11,2)x,Zsz(4)Zsz=0Hsz=0enddoend read(1,*)k,T,vA=45I=45R=15print*,m,vpi=3.14159u0=4.0*pi*10.0**(-7.0)print*,u0M=k*(T/u0)Ze=0Za=0DT=0y=0do x=-80,80,1do Z=-80,80,1R=zce=4*pi*(x**2+y**2+Z**2)**2.5Za=(u0*m*v/ce)*((2.0*Z**2.0-x**2.0-y**2.0)*sind(I)-3*Z*x*cosd(I)*cosd(A)-3*Z*y*cos d(I)*sind(A))Hay=(u0*m*v)/(ce)*((2*y**2-x**2-Z**2)*co sd(I)*sind(A)-3*Z*y*sind(I)+3*x*y*cosd(I)* cosd(A))DT=(u0*m*v)/(ce)*((2*Z**2-x**2-y**2)*(si nd(I)**2)+(2*x**2-y**2-Z**2)*cosd(I)**2*c osd(A)**2+(2*y**2-x**2-Z**2)*Cosd(I)**2 *sind(A)**2-3*x*Z*sind(2*I)*cosd(A)+3*x* y*cosd(I)**2*sind(2*A)-3*y*Z*sind(2*I)*sin d(A))write(2,1)x,Z,Zawrite(3,1)x,z,DTZa=0DT=0enddoHax=0Za=0enddoend四、 图件（见附图）五、 图件分析（1） 等值线平面图分析① 倾斜磁化（图一、图二）从图中可以看出，Za与AT均在X，Y零点附近取得极值。

Za的等值线为等轴状，负异常包围着正异常极大值与极小值的连线对应着磁化强度矢 量在平面上的投影方向△T的等值线也为等轴状，负值所占区域稍大于正值极大值点与极小值点的连线即为主剖面② 垂直极化（图三、图四）从图中看出，垂直极化的Za与AT等值线均为以XY坐标原点为圆心的同心圆外侧 为负值，内部为正值，在XY轴交点处取得极大值2） 空间等值线图分析（图五、图六）可以看出，在y=0的铅垂面内Za与AT的空间等值线均呈现双纽线型，其零值线均 为交于原点的直线，极值点位于原点月靠近原点的等值线对应的绝对值越大3） Za，Ha剖面图分析（图七）可以看出，在远离磁化球体时，两者均趋于零在原点处，is相同时，Za值均大于Ha 值，is=0时的Za曲线与is=90°的Ha曲线基本重合两者均为不对称曲线两者的极大 值都随着is的增大而增大六、感想和收获1. 在研究地球物理问题时，实际情况往往较为复杂，应根据实际地质和物理条件，合理的进 行简化，用极限的思想进行类比2. 计算中计算量较大时可以通过使用快速算法，改进循环方式，注意关闭打开程序的方法减 少内存占用提高程序运行速度80X80n0-80」图一球体Za等值线平面图A=I=45°, R=15图二球体AT等值线平面图A=I=45°, R=15X80-|图三球体Za等值线平面图1=90°。