//Program /home/stebla/programs/vector-helmholtz.sci to study
//transients on a pcb. The solution is obtained in reference [1].
//References
//[1] "An attempt to estimate the effect of a switch-mode power
//supply on a low-level circuit", 13th February 2007, Works Volume XXI.
//[2] "Impedance of a circular printed circuit board", 26th February 2007,
//Works Volume XXI.
//[3] "Solution for a circular pcb with distributed inductance and
//capacitance", 13th March 2007, Works Volume XXII.
//History
//22nd Feb 2007: Created
//26th Feb 2007: Solved for a finite pcb by putting the bc that the 
//surface current density is zero at r=r1.
//26th Feb 2007: The solution was tidied-up by following the note [2]
//because the solution in [1] was a rather messy.
//13th Mar 2007: Copied original to vector-helmholtz-v2.sci and modified
//it to use a frequency-dependent permittivity in order to test the
//theory in reference [3].  
//8th Feb 2008: The solution was reworked from the Helmholtz equation
//for the magnetic field following appendix A of reference [2] because
//we get an expression for the electric field in terms of the port
//current (the previous version was in terms of the port voltage)
//which can then be directly compared with the electric field in the 
//numerical meth in method-of-moments-v5.sci and so is an additional
//way of validating the numerical code.
//11th Feb 2008: Incorporated the solution for a large number of 
//decoupling caps using the method of the effective permittivity.
//Usage
//-->exec('vector-helmholtz.sci');
//-->bode(fvec,Z0);//Bode plot of impedance
//-->plot2d(r,imag(E));//Plot of E-field versus radius
//-->plot2d(r,real(H));//Plot of H-field versus radius
//Program controls
//Board parameters
er=4.2;//Relative permittivity of board dielectric ?
r0=125.0e-6;//Radius of via-post for the impressed voltage ?
r1=100.0e-3;//Radius of pcb ?
h=150.0e-6;//Separation of the planes ?
bare_board_impedance=%F;//%T bare board, %F to add decoupling caps ?
M=100;//Number of decoupling capacitors ?
Cd=100.0e-9;//Nominal decoupling capacitance per device ?
Lesl=1.5e-9;//Effective series inductance per device ?
Resr=60.0e-3;//Effective series resistance per device ?
//Impedance Bode plot parameters
f0=100.0e3;// Start frequency for impedance calc?
no_of_decades=5;//No of decades f0 to f0*10^no_of_decades ?
no_of_data_points=1000;//No of data points for impedance calc ?
//Calculation of E-field versus radius parameters
f=1000.0e6;// Fixed frequency for E-field calc ?
Ie=1;//Impressed current up via-post (port current) ?
r=r0:0.01*r1:r1;//Radial points to calculate E-field solution ?
//SI units
//Physical constants
c=3.0e8;//Speed of light
e0=1.0e7/(4*%pi*c^2);//Permittivity of vacuum
u0=4*%pi*1.0e-7;//Permeability of vacuum
//PCB parameters
//Effective permittivity from decoupling capacitors
s=poly(0,"s");//Variable for the effective permittivity polynomial
Pden=s^2*Lesl*Cd+s*Resr*Cd+1;//Freq-dependence of effective permittivity
ed0=M*Cd*h/(%pi*r1^2);//Zero-frequency addition to permittivity
//Frequency-dependent permittivity function
function eeff=eff_perm(w,ed0,Pden)
  //Computes the effective permittivity for a single species
  //of decoupling capacitor. ed0 is the zero-frequency addition
  //to the usual permittivity and Pden is the polynomial in the
  //denominator.
  eeff=e0*er+ed0/horner(Pden,%i*w);
endfunction;
// ************ Start of the code *******************
w=2*%pi*f;
if bare_board_impedance then
  zeta0=sqrt(u0/(e0*er));//Free space impedance
  k=sqrt(u0*e0*er)*w;//Wavenumber without decoupling
else
  e_eff=eff_perm(w,ed0,Pden);
  zeta0=sqrt(u0/e_eff);//Effective free space impedance
  k=sqrt(u0*e_eff)*w;//Wavenumber with freq-dep permittivity
end;
H10kr0=besselh(0,1,k*r0);//H1nkr=H^1_n(kr) Hankel function of 1st kind
H20kr0=besselh(0,2,k*r0);//H2nkr=H^2_n(kr) Hankel function of 2nd kind
H11kr1=besselh(1,1,k*r1);
H21kr1=besselh(1,2,k*r1);
H11kr0=besselh(1,1,k*r0);
H21kr0=besselh(1,2,k*r0);
for j=1:length(r);
  H10kr=besselh(0,1,k*r(j));
  H20kr=besselh(0,2,k*r(j));
  H11kr=besselh(1,1,k*r(j));
  H21kr=besselh(1,2,k*r(j));
  H(j)=(Ie/(2*%pi*r0))*(H11kr*H21kr1-H21kr*H11kr1)/(H11kr0*H21kr1-H21kr0*H11kr1);
  E(j)=-%i*zeta0*(Ie/(2*%pi*r0))*..
               (H10kr*H21kr1-H20kr*H11kr1)/(H11kr0*H21kr1-H21kr0*H11kr1);
end;
n=no_of_data_points;//No of data points for impedance calc ?
m=no_of_decades;//No of decades f0 to f0*10^m
for j=1:n;
  fvec(j)=f0*10^(m*(j-1)/(n-1));
  w=2*%pi*fvec(j);
  if bare_board_impedance then
    zeta0=sqrt(u0/(e0*er));//Free space impedance
    k=sqrt(u0*e0*er)*w;//Wavenumber without decoupling
  else
    e_eff=eff_perm(w,ed0,Pden);
    zeta0=sqrt(u0/e_eff);//Effective free space impedance
    k=sqrt(u0*e_eff)*w;//Wavenumber with freq-dep permittivity
  end;
  H10kr0=besselh(0,1,k*r0);//H1nkr=H^1_n(kr) Hankel function of 1st kind
  H20kr0=besselh(0,2,k*r0);//H2nkr=H^2_n(kr) Hankel function of 2nd kind
  H11kr1=besselh(1,1,k*r1);
  H21kr1=besselh(1,2,k*r1);
  H11kr0=besselh(1,1,k*r0);
  H21kr0=besselh(1,2,k*r0);
  Z0(j)=(1/(2*%pi*%i))*zeta0*(h/r0)*..
        (H20kr0*H11kr1-H10kr0*H21kr1)/(H11kr0*H21kr1-H21kr0*H11kr1);
end;