資源描述：
《球體重力異常.doc》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在教育資源-天天文庫(kù)。
1、%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Gravityball.m%%功能：球體模型的重力異常%%入口參數(shù)：x,y----平面坐標(biāo)r--球半徑rd--球體剩余密度%h--球體中心距測(cè)點(diǎn)距離%出口參數(shù)：g--重力異常Vxz,Vyz,Vzz--二階導(dǎo)數(shù)Vzzz--三階導(dǎo)數(shù)(剖面）%gg--重力異常%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;x=-30:0.5:30;
2、y=-30:0.5:30;[X,Y]=meshgrid(x,y);網(wǎng)格采樣點(diǎn)（x,y）r=4;rd=0.8;h=3;rm=(4/3)*pi*r^3*rd;G=6.67;%%-----------------計(jì)算重力異常,用矩陣的方法進(jìn)行計(jì)算(方法一的語(yǔ)句)---------g1=G*rm*h./(X.^2+Y.^2+h^2).^(3/2);gg=g1*10^-2;%%%--------------循環(huán)的方法進(jìn)行計(jì)算(方法二的語(yǔ)句)------------%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.)%---------------------------------------------------------%--------------以下均用方法一來(lái)計(jì)算其他分量,也可改寫成循環(huán)(方法二)來(lái)計(jì)算--------------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);可以設(shè)置字體顏色字號(hào)tit
5、le('Gravitycontour');[C,H]=contour(gg,100);title('Gravitycontour');xlabel('x/m');ylabel('△g/g.u.');%--------------------畫立體圖-----------figure(2);subplot(2,2,1);將多個(gè)圖畫在一個(gè)平面上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-');畫函數(shù)圖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
