长方体重力异常正演模拟

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

function gg=cftmx（）
%     均匀长方体重力异常     %
%质心坐标（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） 万有引力常数
%数据保存在‘均匀长方体重力异常.txt’

clear
clc

%长方体模型参数%
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））;
gg=-（g1-g2-g3+g4-g5+g6+g7-g8）*10^5;    %单位mGal

%成图%
figure（1）
mesh（x1y1gg）%三维
xlabel（‘‘）
ylabel（‘‘）
title（‘均匀长方体重力异常‘）
figure（2）
contourf（x1y1gg）%二维
title（‘均匀长方体重力异常‘）


%数据生成文本%
%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