%垂直磁化情况下
clear;
x=1:1:20;
y=1;
z=1:1:10;
num1=0;    %num是累加和
num2=0;
num3=0;
num4=0;
u0=4*pi*10^(-7);          %u0是绝对磁导率
m=0.9;  %质量m赋初值
z0=5;    %地质体埋深
yq=1;
for  zq=1:1:10;
   for  xq=1:1:20;
      for    i=1:1:20;
            xi=i;
            yi=1;
           zi=1;
                  
      r=[(xq-xi)^2+(yq-yi)^2+(zq-zi)^2]^(1/2);           %r是第q个偶极子到测点的距离
     T(xi,yi,zi)=u0/(4*pi)*m/{[xi^2+yi^2+(zi-z0)^2]^(5/2)}*[2*(zi-z0)^2-xi^2-yi^2];
  %T是测点磁异常
     Bq(xi,yi,zi)=(zq-zi)/[(xq-xi)^2+(yq-yi)^2+(zq-zi)^(3/2)];
  %Bq是基函数
    num1=num1+T(xi,yi,zi)*Bq(xi,yi,zi);
    num2=num2+T(xi,yi,zi)^2;
    num3=num3+Bq(xi,yi,zi)^2;
    num4=(num4*num3)^(1/2);
    cq=num1/num4;
           end
        end 
    end 
    Plot(xq,yq,cq);