clear
l=1;
for p=0.3:0.05:0.7
N=80;
b=10;
h=b/N;
a=p*b;
s=(1-p)*N;
s0=fix(s);
x=(1-p)*b-s0*(b/N);
if x>=0.5*h
    s=s0+1;
else
    s=s0;
end
w=2/(1+sin(pi/N));
We=2*N-1;
H(1:We)=0;
I=0;
U(1:N+1,1:N+1)=0;
U(1,1:s+1)=1;


% 求電位分布
for A=1:60
    for c=1:N+1
        for d=1:N+1
            if c==1 & d>=1 & d<=s+1
                U(c,d)=1;
            elseif c==N+1
                U(c,:)=0;
            elseif d==N+1
                U(:,d)=0;
            elseif c==1 & d>s+1 & d<N+1
                U(c,d)=U(c,d)+0.25*w*(2*U(c+1,d)+U(c,d+1)+U(c,d-1)-4*U(c,d));
            elseif d==1 & c>1 & c<N+1
                U(c,d)=U(c,d)+0.25*w*(U(c+1,d)+2*U(c,d+1)+U(c-1,d)-4*U(c,d));
            elseif c==N+1 | d==N+1
                U(c,d)=0;
            else
                U(c,d)=U(c,d)+0.25*w*(U(c+1,d)+U(c,d+1)+U(c-1,d)+U(c,d-1)-4*U(c,d));
            end
        end
    end
end

% 求電場、磁場

for e=1:(We+1)/2
    E(e)=0.5*((U(e,N)-U(e,N+1))+(U(e+1,N)-U(e+1,N+1)))/(b/N);
end
k=1;
for e=(We+1)/2+1:We
    E(e)=0.5*((U(N,k)-U(N+1,k))+(U(N,k+1)-U(N+1,k+1)))/(b/N);
    k=k+1;
end
H=E/(120*pi);
    
% 求電流、特性阻抗
for f=1:We
I=I+H(f)*(b/N);
end

Zc(l)=1/(4*I);
l=l+1;
end

p=0.3:0.05:0.7;
plot(p,Zc)
hold on