%均勻球體與長方體異常疊加%

clear
clc

%球體參數%
r=100;????????%球的半徑
pd=1*10^3;????%剩余密度
x0=0;?????????%球心坐標（x0y0z0）
y0=0;
D=200;????????%球心埋深z0

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

%常數%
G=6.67e-11;

%計算異常%
[x1y1]=meshgrid（xy）;???%生成網線節點矩陣
gg1=G*（（4/3）*pi*r^3*pd）*（D-z）./（（（x1-x0）.^2+（y1-y0）.^2+（D-z）^2）.^（3/2））*10^5;%單位mGal

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

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

%常數%
G=6.67e-11;

%計算異常%
[x1y1]=meshgrid（xy）;???%生成網線節點矩陣
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））;
gg2=-（g1-g2-g3+g4-g5+g6+g7-g8）*10^5;????%單位mGal

%計算總異常%
gg=gg1+gg2;

%成圖%
x=（-2000:50:2000）;
y=（-2000:50:2000）;
[x1y1]=meshgrid（xy）;%生成網線節點矩陣
figure（1）%圖1
mesh（x1y1gg）%三維
xlabel（‘‘）
ylabel（‘‘）
title（‘均勻球體長方體疊加異常‘）
figure（2）
contourf（x1y1gg‘levelstep‘0.1）%二維
title（‘均勻球體長方體疊加異常‘）

%數據生成文本%
%t=[x1（:）‘
%???y1（:）‘
%???gg（:）‘];
%fid=fopen（‘球體_長方體.txt‘‘wt‘）;?%wt以文本格式寫入
%fprintf（fid‘%4.2f?%4.2f?%.2e\n‘t）;
%fclose（fid）;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件????????2762??2015-12-01?21:51??qt_cft.m