ballps.m

应用地球物理

function ballPs()
%计算球体的标准曲线
%
type = input('输入装置类型是（1、二极；2、联合剖面；3、对称四极；4、偶极）：');
ptype = input('输入（1、高阻；2、低阻）：');
x = linspace(-50,50,100);
z0 = 6;
r0 = 2;
p1 = 100;
if ptype == 1
    p2 = 10*p1;
else
    p2 = p1/10;
end
ps = zeros(5,100);
psa = zeros(5,100);
psb = zeros(5,100);
px = zeros(5,100);
py = zeros(5,100);
switch type
    case 1
        figure(1);
        for i = 1:5
            xa = x - r0*i/2;
            xm = x + r0*i/2;
            k = 2*pi*r0*i;
            u = ballu(0,0,z0,xa,0,0,xm,0,0,p1,p2,r0,1);
            ps(i,:) = k*u;
        end
        plot(x,ps','LineWidth',5)
    case 2
        figure(1);
        for i = 1:5
            xa = x - r0*(i+i/2);
            xm = xa + r0*i;
            xn = xm + r0*i;
            xb = xn + r0*i;
            k = 4*pi*r0*i;
            uam = ballu(0,0,z0,xa,0,0,xm,0,0,p1,p2,r0,1);
            uan = ballu(0,0,z0,xa,0,0,xn,0,0,p1,p2,r0,1);
            ubm = ballu(0,0,z0,xb,0,0,xm,0,0,p1,p2,r0,-1);
            ubn = ballu(0,0,z0,xb,0,0,xn,0,0,p1,p2,r0,-1);
            psa(i,:) = k*(uam - uan);
            psb(i,:) = k*(ubm - ubn);
            subplot(5,1,i); plot(x,psa(i,:),'-',x,psb(i,:),'--','LineWidth',5);
        end
    case 3
        figure(1);
        for i = 1:5
            xa = x - r0*(i+i/2);
            xm = xa + r0*i;
            xn = xm + r0*i;
            xb = xn + r0*i;
            k = 2*pi*r0*i;
            uam = ballu(0,0,z0,xa,0,0,xm,0,0,p1,p2,r0,1);
            uan = ballu(0,0,z0,xa,0,0,xn,0,0,p1,p2,r0,1);
            ubm = ballu(0,0,z0,xb,0,0,xm,0,0,p1,p2,r0,-1);
            ubn = ballu(0,0,z0,xb,0,0,xn,0,0,p1,p2,r0,-1);
            ps(i,:) = k*((uam+ubm) - (uan+ubn));
        end
        plot(x,ps,'LineWidth',5)
    case 4
        figure(1);
        for i = 1:5
            xa = x - r0*(1+i*2/2);
            xb = xa + r0;
            xm = x + r0*i*2/2;
            xn = xm + r0;
            k = pi*r0*2*i*(i*2+1)*(i*2+2);
            uam = ballu(0,0,z0,xa,0,0,xm,0,0,p1,p2,r0,1);
            uan = ballu(0,0,z0,xa,0,0,xn,0,0,p1,p2,r0,1);
            ubm = ballu(0,0,z0,xb,0,0,xm,0,0,p1,p2,r0,-1);
            ubn = ballu(0,0,z0,xb,0,0,xn,0,0,p1,p2,r0,-1);
            ps(i,:) = -k*((uam+ubm) - (uan+ubn));
            px(i,:) = x;
            py(i,:) = r0*(1+i*2/2);
        end
        plot(x,ps,'LineWidth',5),
        figure(2)
        contour(px,py,ps);
end



function u = ballu(x0,y0,z0,xa,ya,za,xp,yp,zp,p1,p2,r0,I)
%球体地表电位
%
R = distance(xa,ya,za,xp,yp,zp);
d = distance(xa,ya,za,x0,y0,z0);
r = distance(xp,yp,zp,x0,y0,z0);
sumleg = zeros(1,100);
for i = 1:5
    sumleg = sumleg +2*(p2-p1)*i/((i+1)*p2+i*p1)*r0^(2*i+1)./(d.*r).^(i+1).*mylegendre(i,((d.^2+r.^2-R.^2)./(2*d.*r)));
end
u = I*p1*(1./R+sumleg)/(2*pi);


function d=distance(xa,ya,za,xp,yp,zp)
%
d = sqrt((xa-xp).^2+(ya-yp).^2+(za-zp).^2);

function legd = mylegendre(n,x)
%
if n == 0
    legd  = ones(size(x));
elseif n == 1
    legd = x;
else
    legd = ((2*(n-1)+1).*x.*mylegendre(n-1,x) - (n-1)*mylegendre(n-2,x))/n;
end