function Delta_T = qiuti_Delta_T(x,y)
% 程序功能:正演Delta_T
% 日期:2018-12-10
%------------------------------------------------%
D = x^2+y^2;
h = 15; %球体中心埋深
B = 50000; %本地磁场(nT)
k = 0.015; %磁化率(SI)
r = 10; %球体半径
vol = 4 * pi * r^3 / 3; %球体体积
mu_0 = 4 * pi * 1e-7; %导磁率
M = k * B / mu_0; %球体磁化强度
m = M * vol; %球体磁矩
i = pi / 4; %磁倾角
I0 = pi / 4; %磁偏角
% i = 0;
% I0 = pi/4;
Za = mu_0 * m * ((2 * h ^ 2 - D) * sin(I0) - 3 * x * h * cos(I0) * cos(i) - 3 * y * h * cos(I0) * sin(i)) / (4 * pi * (D + h ^ 2) ^ 2.5);
Hax = mu_0 * m * ((2 * x ^ 2 - y ^ 2 - h ^ 2) * cos(I0) * cos(i) - 3 * x * h * sin(I0) + 3 * x * y * cos(I0) * sin(i))/(4 * pi * (D + h ^ 2) ^ 2.5);
Hay = mu_0 * m * ((2 * y ^ 2 - x ^ 2 - h ^ 2) * cos(I0) * cos(i) - 3 * y * h * sin(I0) + 3 * x * y * cos(I0) * cos(i))/(4 * pi * (D + h ^ 2) ^ 2.5);
Delta_T = Hax * cos(I0) * cos(i) + Hay * cos(I0) * sin(i) + Za * sin(I0);
end