clf
clc
clear all
a=0.1;
A=0.2:0.1:1;
s=size(A);
N1=50;
e=8.854e-12;
mu=pi*4e-7;
for K=1:s(2)
  AA=A(K);
    dal_n=2*pi*a/N1;
    N2=2*pi*AA/dal_n;
    div1=2*pi/N1;  
    div2=2*pi/N2;
    n=0;
    for B=0:div1:2*pi
        n=n+1;
        X(n)=a*cos(n*div1);
        Y(n)=a*sin(n*div1);
        dx(n)=dal_n*cos(abs(pi-B-(pi-div1)/2));  
        dy(n)=dal_n*sin(abs(pi-B-(pi-div1)/2));
    end
    for n=1:N1         
        x(n)=(X(n)+X(n+1))/2;
        y(n)=(Y(n)+Y(n+1))/2;
    end
    n=N1;
    for B=0:div2:2*pi
        n=n+1;
        X(n)=AA*cos(n*div2);
        Y(n)=AA*sin(n*div2);
        dx(n)=dal_n*cos(abs(pi-B-(pi-div2)/2));   
        dy(n)=dal_n*sin(abs(pi-B-(pi-div2)/2));
    end
    for n=N1+1:N1+N2         
        x(n)=(X(n)+X(n+1))/2;
        y(n)=(Y(n)+Y(n+1))/2;
    end
    t1=-0.5;
    t2=0.5;
    for m=1:N1+N2
        for n=1:N1+N2
            if m==1
                s(m,n)=dal_n;
            else
                s(m,n)=0;
            end
        end
    end
    for m=2:N1+N2
        for n=1:N1+N2
            a1=(x(1)-x(n))^2+(y(1)-y(n))^2;
            b1=-2*(x(1)-x(n))*dx(n)+2*(y(1)-y(n))*dy(n);
            c1=dx(n)^2+dy(n)^2;
            a2=(x(m)-x(n))^2+(y(m)-y(n))^2;
            b2=-2*(x(m)-x(n))*dx(n)-2*(y(m)-y(n))*dy(n);
            c2=c1;
            d1=b1^2-4*a1*c1;
            d2=b2^2-4*a2*c2;
            if d1>0
                s11=(t1+b1/2/c1)*log(a1+b1*t1+c1*t1^2)-2*t1+(sqrt(d1)/c1)*atanh((2*c1*t1+b1)/sqrt(d1));
                s12=(t2+b1/2/c1)*log(a1+b1*t2+c1*t2^2)-2*t2+(sqrt(d1)/c1)*atanh((2*c1*t2+b1)/sqrt(d1));
            else if d1==0
                    s11=(t1+b1/2/c1)*log(a1+b1*t1+c1*t1^2)-2*t1;
                    s12=(t2+b1/2/c1)*log(a1+b1*t2+c1*t2^2)-2*t2;
                else
                    s11=(t1+b1/2/c1)*log(a1+b1*t1+c1*t1^2)-2*t1+(sqrt(-d1)/c1)*atan((2*c1*t1+b1)/sqrt(-d1));
                    s12=(t2+b1/2/c1)*log(a1+b1*t2+c1*t2^2)-2*t2+(sqrt(-d1)/c1)*atan((2*c1*t2+b1)/sqrt(-d1));
                end
            end
            s1=s12-s11;
            if d2>0
                s21=(t1+b2/2/c2)*log(a2+b2*t1+c2*t1^2)-2*t1+(sqrt(d2)/c2)*atanh((2*c2*t1+b2)/sqrt(d2));
                s22=(t2+b2/2/c2)*log(a2+b2*t2+c2*t2^2)-2*t2+(sqrt(d2)/c2)*atanh((2*c2*t2+b2)/sqrt(d2));
            else if d2==0
                    s21=(t1+b2/2/c2)*log(a2+b2*t1+c2*t1^2)-2*t1;
                    s22=(t2+b2/2/c2)*log(a2+b2*t2+c2*t2^2)-2*t2;
                else
                    s21=(t1+b2/2/c2)*log(a2+b2*t1+c2*t1^2)-2*t1+(sqrt(-d2)/c2)*atan((2*c2*t1+b2)/sqrt(-d2));
                    s22=(t2+b2/2/c2)*log(a2+b2*t2+c2*t2^2)-2*t2+(sqrt(-d2)/c2)*atan((2*c2*t2+b2)/sqrt(-d2));
                end
            end
            s2=s22-s21;
            s(m,n)=dal_n*s1-dal_n*s2;
        end
    end
    v(1:N1)=0;
    v(N1+1:N1+N2)=-1;
    Q=4*pi*e*inv(s)*v';
    C1(K)=sum(Q(1:N1))*dal_n;
    L1(K)=mu*e/C1(K);
    Z1(K)=sqrt(L1(K)/C1(K));
    C2(K)=2*pi*e/log(AA/a);
    L2(K)=mu*log(AA/a)/(2*pi);
    Z2(K)=sqrt(mu/e)*log(AA/a)/(2*pi);
end
figure(1)
plot(A,C1,'b',A,C2,'r*');
xlabel('D/d(外直徑與內(nèi)直徑之比)');
ylabel('單位長度電容(F/m)');
legend('矩量解','解析解');
title('單位長度電容隨D/d的變化曲線');
grid on
figure(2)
plot(A,L1,'b',A,L2,'r*');
xlabel('D/d(外直徑與內(nèi)直徑之比)');
ylabel('單位長度電感(H/m)');
legend('矩量解','解析解');
title('單位長度電感隨D/d的變化曲線');
grid on
figure(3)
plot(A,Z1,'b',A,Z2,'r*');
xlabel('D/d(外直徑與內(nèi)直徑之比)');
ylabel('特性阻抗(Ω)');
legend('矩量解','解析解');
title('特性阻抗隨D/d的變化曲線');
grid on