########################################
#  This script will calculate the 
#  impedance of a coaxial cable
#
########################################



########################################
# Theory

# From Data Sheet
vr=0.85;          # fraction of c
r=1.024e-3/2;      # inner diameter (18AWG)
Z=75;             # impedance 

# Calculated
epsr=1/vr^2;
n=sqrt(epsr);
R=r*exp( 2*pi*Z*sqrt(eps0/mu0)*sqrt(epsr));
fc=c/( pi*(R+r)*sqrt(epsr) ); 

# Summary
?"  Theory" + endl +
 "  epsr:"+num2str(epsr) + "  n:"+num2str(n) + "  phase velocity:"+num2str(vr)+"c" + endl +
 "  R (mm):"+num2str(R*1e3) +  "  r (mm):"+num2str(r*1e3) +  endl +
 "  Z:"+num2str(Z) +  "  fc (GHz):"+num2str(fc*1e-9) + endl;



#########################################
# Simulated

# get data directly from simulation
r=getnamed("inner","radius");
R=getnamed("ring","inner radius");

m="mode1";
n=real( getdata(m,"neff") );
x=getdata(m,"x");
y=getdata(m,"y");
Ex=pinch(getdata(m,"Ex"));
Ey=pinch(getdata(m,"Ey"));
Ez=pinch(getdata(m,"Ez"));
Hx=pinch(getdata(m,"Hx"));
Hy=pinch(getdata(m,"Hy"));
Hz=pinch(getdata(m,"Hz"));

# calculations
epsr=n^2;
vc=1/n;
fc=c/( pi*(R+r)*sqrt(epsr) ); 

# calculate Er and Hphi
X=meshgridx(x,y);
Y=meshgridy(x,y);
phi=atan2(Y,X);
Er = Ex*cos(phi)+Ey*sin(phi);
Hphi=Hy*cos(phi)-Hx*sin(phi);


# Calculate V=int(E.dr)
# The easiest way to do this integral 
# is integrating Ex along x axis.  

y0=find(y,0);
x0=find(x,0);

E_line=pinch(Ex,2,y0);
E_line=E_line(x0:length(x));
line=x(x0:length(x));

V = integrate(E_line,1,line);


# calculate I=int(H.ds)
# Integrate Hphi around wire
# This requires Hphi be interpolated onto 
# a circular path around the inner conductor

res=1000;                            # number of points in integral
r_int=(R+r)/2;                       # radius of integral
phi_int=linspace(0,360,res)*pi/180;  # vector of angles
x_int=r_int*cos(phi_int);            # x positions coresponding to integral radius and phi
y_int=r_int*sin(phi_int);            # y positions coresponding to integral radius and phi

line=linspace(0,2*pi*r_int,res);     # ds vector (distance around integral line)
H_line=matrix(res);                  # initialize matrix

for (i=1:res) {                      # interpolate H to points on integral line
  H_line(i) = interp(Hphi,x,y,x_int(i),y_int(i));
}

I=integrate(H_line,1,line);


# calculate Z
Z = V/I;


# Summary
?"  Simulation" + endl +
 "  epsr:"+num2str(epsr) + "  n:"+num2str(n) + "  phase velocity:"+num2str(vc) + "c" + endl +
 "  R (mm):"+num2str(R*1e3) + "  r (mm):"+num2str(r*1e3) + endl +
 "  Z:"+num2str(Z) +  "  fc (GHz):"+num2str(fc*1e-9) + endl;

image(x*1e3,y*1e3,Er,"x (mm)","y (mm)","Er");
image(x*1e3,y*1e3,Hphi,"x (mm)","y (mm)","Hphi");


analysis;  # switch back to the analysis window