资源描述：
《球体重力异常.doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Gravityball.m%%功能：球体模型的重力异常%%入口参数：x,y----平面坐标r--球半径rd--球体剩余密度%h--球体中心距测点距离%出口参数：g--重力异常Vxz,Vyz,Vzz--二阶导数Vzzz--三阶导数(剖面）%gg--重力异常%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;x=-30:0.5:30;
2、y=-30:0.5:30;[X,Y]=meshgrid(x,y);网格采样点（x,y）r=4;rd=0.8;h=3;rm=(4/3)*pi*r^3*rd;G=6.67;%%-----------------计算重力异常,用矩阵的方法进行计算(方法一的语句)---------g1=G*rm*h./(X.^2+Y.^2+h^2).^(3/2);gg=g1*10^-2;%%%--------------循环的方法进行计算(方法二的语句)------------%form=1:1:121%forn=1:1:121%gg(m,n)=G*rm*h/(x(n)^2+y(m)^2+
3、h^2)^(3/2);%end%end%gg=gg*10^-2;%(gg:g.u.)%---------------------------------------------------------%--------------以下均用方法一来计算其他分量,也可改写成循环(方法二)来计算--------------V1xz=-3*G*rm*h*X./(X.^2+Y.^2+h^2).^(5/2);Vxz=V1xz*10;%(Vxz:E)1E=10^-91/s^2V1yz=-3*G*rm*h*Y./(X.^2+Y.^2+h^2).^(5/2);Vyz=V1yz*10;
4、V1zz=G*rm*(2*h^2-X.^2-Y.^2)./(X.^2+Y.^2+h^2).^(5/2);Vzz=V1zz*10;V1zzz=3*G*rm*(2*h^2-3*X.^2-3*Y.^2)./(X.^2+Y.^2+h^2).^(7/2);Vzzz=V1zzz*10^4;%(Vzzz:pMKS)1pMKS=10^-121/m*s^2%-----------------画等值线图-----------figure(1);%[C,H]=contour(gg,20);等值线colormap(hsv);控制曲面图的颜色clabel(C,H);可以设置字体颜色字号tit
5、le('Gravitycontour');[C,H]=contour(gg,100);title('Gravitycontour');xlabel('x/m');ylabel('△g/g.u.');%--------------------画立体图-----------figure(2);subplot(2,2,1);将多个图画在一个平面上surf(Vxz);画立体图title('Vxz');....subplot(2,2,2);surf(Vyz);title('Vyz');....subplot(2,2,3);surf(Vzz);title('Vzz');....
6、subplot(2,2,4);surf(Vzzz);title('Vxz');....%--------------------画剖面图-----------figure(3);subplot(2,1,1);plot(x,gg(61,:),'r-');画函数图legend('gg(11)');加图例title('Gravitymainprofile')set(gca,'XLim',[-3030],'YLim',[02]);xlabel('x/m');ylabel('△g/g.u.');subplot(2,1,2);Vxz_max=max(Vxz(61,:));Vzz
7、_max=max(Vzz(61,:));Vzzz_max=max(Vzzz(61,:));plot(x,Vxz(61,:)/Vxz_max,'b',x,Vzz(61,:)/Vzz_max,'r',x,Vzzz(61,:)/Vzzz_max,'g');set(gca,'XLim',[-3030],'YLim',[-1.51.5]);xlabel('x/m');title('Derivatesofmaingravityprofile')legend('Vxz(11)','Vzz(11)','Vzzz(11)');ans='OK!'67
