均匀球体与长方体重力异常正演模拟Matlab代码

均匀球体与长方体重力异常正演模拟Matlab代码 地球物理勘探

%均匀球体与长方体异常叠加%

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