function?gg=cftmx（）
%?????均勻長方體重力異常?????%
%質(zhì)心坐標（x0y0z0）
%埋深H（m）
%長為a（m），寬為b（m），高為c（m），剩余密度p（kg/m^3）
%xyz為采樣點
%重力異常gg?（m/s^2）
%G=6.67e-11（m^3/kg.s^2）?萬有引力常數(shù)
%數(shù)據(jù)保存在‘均勻長方體重力異常.txt’

clear
clc

%長方體模型參數(shù)%
a=2000;???????%長
b=200;???????%寬
c=100;???????%高
x0=0;????????%質(zhì)心坐標9x0y0z0）
y0=0;
H=1000;???????%質(zhì)心埋深z0
ph=2*10^3;???%剩余密度

%采樣區(qū)間%
x=（-2000:50:2000）;
y=（-2000:50:2000）;
z=0;

%常數(shù)%
G=6.67e-11;

%計算異常%
[x1y1]=meshgrid（xy）;???%生成網(wǎng)線節(jié)點矩陣
r1=sqrt（（x0+a/2-x1）.^2+（y0+b/2-y1）.^2+（H+c/2-z）.^2）;
r2=sqrt（（x0+a/2-x1）.^2+（y0+b/2-y1）.^2+（H-c/2-z）.^2）;
r3=sqrt（（x0+a/2-x1）.^2+（y0-b/2-y1）.^2+（H+c/2-z）.^2）;
r4=sqrt（（x0+a/2-x1）.^2+（y0-b/2-y1）.^2+（H-c/2-z）.^2）;
r5=sqrt（（x0-a/2-x1）.^2+（y0+b/2-y1）.^2+（H+c/2-z）.^2）;
r6=sqrt（（x0-a/2-x1）.^2+（y0+b/2-y1）.^2+（H-c/2-z）.^2）;
r7=sqrt（（x0-a/2-x1）.^2+（y0-b/2-y1）.^2+（H+c/2-z）.^2）;
r8=sqrt（（x0-a/2-x1）.^2+（y0-b/2-y1）.^2+（H-c/2-z）.^2）;
g1=G*ph*（（x0+a/2-x1）.*log（（y0+b/2-y1）+r1）+（y0+b/2-y1）.*log（（x0+a/2-x1）+r1）+（H+c/2-z）.*atan（（x0+a/2-x1）.*（y0+b/2-y1））./（（H+c/2-z）.*r1））;
g2=G*ph*（（x0+a/2-x1）.*log（（y0+b/2-y1）+r2）+（y0+b/2-y1）.*log（（x0+a/2-x1）+r2）+（H-c/2-z）.*atan（（x0+a/2-x1）.*（y0+b/2-y1））./（（H-c/2-z）.*r2））;
g3=G*ph*（（x0+a/2-x1）.*log（（y0-b/2-y1）+r3）+（y0-b/2-y1）.*log（（x0+a/2-x1）+r3）+（H+c/2-z）.*atan（（x0+a/2-x1）.*（y0-b/2-y1））./（（H+c/2-z）.*r3））;
g4=G*ph*（（x0+a/2-x1）.*log（（y0-b/2-y1）+r4）+（y0-b/2-y1）.*log（（x0+a/2-x1）+r4）+（H-c/2-z）.*atan（（x0+a/2-x1）.*（y0-b/2-y1））./（（H-c/2-z）.*r4））;
g5=G*ph*（（x0-a/2-x1）.*log（（y0+b/2-y1）+r5）+（y0+b/2-y1）.*log（（x0-a/2-x1）+r5）+（H+c/2-z）.*atan（（x0-a/2-x1）.*（y0+b/2-y1））./（（H+c/2-z）.*r5））;
g6=G*ph*（（x0-a/2-x1）.*log（（y0+b/2-y1）+r6）+（y0+b/2-y1）.*log（（x0-a/2-x1）+r6）+（H-c/2-z）.*atan（（x0-a/2-x1）.*（y0+b/2-y1））./（（H-c/2-z）.*r6））;
g7=G*ph*（（x0-a/2-x1）.*log（（y0-b/2-y1）+r7）+（y0-b/2-y1）.*log（（x0-a/2-x1）+r7）+（H+c/2-z）.*atan（（x0-a/2-x1）.*（y0-b/2-y1））./（（H+c/2-z）.*r7））;
g8=G*ph*（（x0-a/2-x1）.*log（（y0-b/2-y1）+r8）+（y0-b/2-y1）.*log（（x0-a/2-x1）+r8）+（H-c/2-z）.*atan（（x0-a/2-x1）.*（y0-b/2-y1））./（（H-c/2-z）.*r8））;
gg=-（g1-g2-g3+g4-g5+g6+g7-g8）*10^5;????%單位mGal

%成圖%
figure（1）
mesh（x1y1gg）%三維
xlabel（‘‘）
ylabel（‘‘）
title（‘均勻長方體重力異?！?br/>figure（2）
contourf（x1y1gg）%二維
title（‘均勻長方體重力異?！?br/>

%數(shù)據(jù)生成文本%
%t=[x1（:）‘
%???y1（:）‘
%???gg（:）‘];
%fid=fopen（‘均勻長方體重力異常.txt‘‘wt‘）;?%wt以文本格式寫入
%fprintf（fid‘%4.2f?%4.2f?%.2e\n‘t）;
%fclose（fid）;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件????????2487??2015-12-01?21:13??cftmx.m