//Program /home/stebla/programs/scilab/method-of-moments-v5.sci
//to solve the scalar Helmholtz equation for the E-field on a pcb
//in the presence of decoupling capacitors. Lines ending in ? 
//can be editted to change the program parameters.
//
//References
//[1] "The impedance of a pcb with discrete decoupling capacitors
//calculated using the method of moments of R.F. Harrington", 4th 
//March 2007, Works Volume XXII.
//History
//6th Mar 2007: Created.
//7th Mar 2007: Copied the original method-of-moments.sci to 
//method-of-moments-v2.sci because the original had a fundamental
//error which calculated the low frequency impedance as inductive, 
//which is unphysical. This version solves the new integral equation
//in section 11 of reference [1].  
//8th Mar 2007: Copied the method-of-moments-v2.sci to 
//method-of-moments-v3.sci because the v2 variant does not have the 
//correct behaviour at low k. The v3 variant uses Scilab's canned
//integration routine intg to calculate the matrix elements move exactly,
//in the hope that this will fix the problem at low k - it does!
//8th Mar 2007: Copied the method-of-moments-v3.sci to 
//method-of-moments-v4.sci . The v3 variant has the correct behaviour
//at low k, but it achieves this by doing the "boundary field due to
//boundary field" integral using Scilab's canned integration routine
//intg; this is slow, and we don't need this precision at higher k
//when the Hankel function is better-behaved. Consequently, the v4
//variant has a k cut-off; below the cut-off we use intg, above the
//cut-off we use the crude approximation in the v2 variant.
//8th Mar 2007: method-of-moments-v4-inst1.sci adds the ability
//to set up port coords to put several species of capacitors on a
//rectangular pcb.
//4th Feb 2008: Tidied up the code to make one canonical version which
//can be put on cvs and so all the previous versions can be forgotten.
// The new version is method-of-moments-v5.sci .
//Usage
//-->exec('method-of-moments-v5.sci');
//-->plot2d(rvec,image(Ezvec));//Plot calculated E-field versus radius
//clear();
//Program controls
circular_board=1;//Circular board shape - do not edit
rectangular_board=2;//Rectangular board shape - do not edit
arbitrary_board=3;//Arbitrary board shape - do not edit
board_shape=rectangular_board;//Generate vertices automatically for a circular board ?
if board_shape==circular_board then
  r1=100.0e-3;//radius for a circular board ?
  N=16;//No of boundary edges ?
elseif board_shape==rectangular_board then  
  a=100.0e-3;//length for a rectangular board ?
  b=100.0e-3;//width for a rectangular board ?
  Nx=5;//Segments along x edge for rect brd ?
  Ny=5;//Segments along y edge for rect brd ?
else //
  N=20;//No of boundary vertices for an arbitrary board ? 
end;
er=4.2;//Relative permittivity of board dielectric ?
r0=125.0e-6;//Radius of the via-posts ?
h=100.0e-6;//Separation of the planes ?
plot_boundary_and_ports=%T;// ?
board_display=[-0.1,-0.1,0.1,0.1];//board display rect [xmin,ymin,xmax,ymax] ?
calc_field=%F;// ?
f=1000.0e6;// Fixed frequency for E-field calc ?
if board_shape==circular_board then
  rvec=r0:0.01*r1:r1;//Radial points to calculate E-field solution ?
else
  rvec=r0:0.001:0.1;//Radial points to calculate E-field solution ?
end;
calc_impedance=%T;// ?
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=100;//No of data points for impedance calc ?
plot_bare_board_impedance=%T;// ?
plot_power_plane_impedance=%T;// ?
fcutoff=50.0e9;//Use intg below this frequency ?
// ************** Decoupling caps ? **************
//These arrays specify the decoupling caps by species. The array sn is the 
//number of caps in each species on the board. C is the cap value for each species,
//Lesl is the effective series inductance for each species and Resr is the effective
//series resistance for each species. These values are just used in the argument 
//list of the function standard_cap to calculate the discrete impedance Zd. 
//If the model for the TDK MLCC caps is required, the if statement in the code
//has to be modified to call the mlcc function instead of standard_cap, because
//the MLCC model uses extra parameters.
sn=  [1       ,3      ,3      ,0       ,37      ,0       ,0       ,0       ,0       ];
C=   [470.0e-6,10.0e-6,1.0e-6 ,330.0e-9,100.0e-9,33.0e-9 ,10.0e-9 ,3.3e-9  ,1.0e-9  ];
Lesl=[4.0e-9  ,1.5e-9 ,1.5e-9 ,1.5e-9  ,1.5e-9  ,1.5e-9  ,1.5e-9  ,1.5e-9  ,1.5e-9  ];
Resr=[19.0e-3 ,20.0e-3,20.0e-3,40.0e-3 ,60.0e-3 ,130.0e-3,200.0e-3,300.0e-3,400.0e-3];
//***********************************************
function v=vertices_for_circular_board(N,r1)
  //Returns the N vertices for a circular board of radius r1.
  v=zeros(N,2);//Boundary vertices
  for i=1:N,
    phi=(i-1)*2*%pi/N;
    v(i,:)=[r1*cos(phi),r1*sin(phi)];
  end;
endfunction;
function [v,N]=vertices_for_rectangular_board(Nx,Ny,a,b)
  //Returns the N vertices for a rectangular board of size a x b
  //with Nx segments along x and Ny segments along y. 
  N=2*(Nx+Ny);
  da=a/Nx;
  db=b/Ny;
  v=zeros(N,2);//Boundary vertices
  for i=1:Nx,
    v(i,:)=[(i-1)*da-0.5*a,-0.5*b];
  end;
  for j=1:Ny,
    v(Nx+j,:)=[0.5*a,(j-1)*db-0.5*b];
  end;
  for i=1:Nx,
    v(Nx+Ny+i,:)=[0.5*a-(i-1)*da,0.5*b];
  end;
  for j=1:Nx,
    v(2*Nx+Ny+j,:)=[-0.5*a,0.5*b-(j-1)*db];
  end;
endfunction;

function pe=random_ports_for_circular_board(M,r1)
  //A function to automatically generate M ports uniformly
  //over a circular board of radius r1.
  pe=zeros(M,2);
  pe(1,:)=1.0e-3*[0,0];//Driver port coords at centre
  for i=2:M,
    u=rand();//Two random number on [0,1]
    v=rand();
    r=r1*sqrt(u);
    theta=2*%pi*v;
    pe(i,:)=[r*cos(theta),r*sin(theta)];
  end;
endfunction;
function pe=random_ports_for_rectangular_board(M,a,b)
  //A function to automatically generate M ports uniformly
  //over a rectangular board of dimensions a x b.
  pe=zeros(M,2);
  pe(1,:)=1.0e-3*[0,0];//Driver port coords at centre
  for i=2:M,
    u=rand();//Two random number on [0,1]
    v=rand();
    x=a*(u-0.5);
    y=b*(v-0.5);
    pe(i,:)=[x,y];
  end;
endfunction;
function pe=hand_crafted_ports(M)
  // A function to encapsulate the edits for hand-crafted 
  //ports. Clearly M must agree with the vertices 
  //actually assigned.
  pe=zeros(M,2);
  pe(1,:)=1.0e-3*[0,0];//Driver port coords 
  pe(2,:)=1.0e-3*[-15,-30];
  pe(3,:)=1.0e-3*[-15,-35];
  pe(4,:)=1.0e-3*[-10,-40];
  pe(5,:)=1.0e-3*[-10,-30];
  pe(6,:)=1.0e-3*[-10,-35];
  pe(7,:)=1.0e-3*[-5,-40];
  pe(8,:)=1.0e-3*[-5,-35];
  pe(9,:)=1.0e-3*[-5,-30];
endfunction;
function v=hand_crafted_vertices(N)
  // A function to encapsulate the edits for hand-crafted 
  //vertices. Clearly N must agree with the vertices 
  //actually assigned.
  v=zeros(N,2);//Boundary vertices
  v(1,:)=1.0e-3*[-50,-50];
  v(2,:)=1.0e-3*[-30,-50];
  v(3,:)=1.0e-3*[-10,-50];
  v(4,:)=1.0e-3*[10,-50];
  v(5,:)=1.0e-3*[30,-50];
  v(6,:)=1.0e-3*[50,-50];
  v(7,:)=1.0e-3*[50,-30];
  v(8,:)=1.0e-3*[50,-10];
  v(9,:)=1.0e-3*[50,10];
  v(10,:)=1.0e-3*[50,30];
  v(11,:)=1.0e-3*[50,50];
  v(12,:)=1.0e-3*[30,50];
  v(13,:)=1.0e-3*[10,50];
  v(14,:)=1.0e-3*[-10,50];
  v(15,:)=1.0e-3*[-30,50];
  v(16,:)=1.0e-3*[-50,50];
  v(17,:)=1.0e-3*[-50,30];
  v(18,:)=1.0e-3*[-50,10];
  v(19,:)=1.0e-3*[-50,-10];
  v(20,:)=1.0e-3*[-50,-30];
endfunction;
function u=direction(p1,p2)
  //unit vector from p2 to p1; u=(p1-p2)/|p1-p2|
  u=p1-p2;
  len=sqrt(u*u');
  u=u/len;
endfunction;
function d=distance(p1,p2)
  //distance from p2 to p1; d=|p1-p2|
  u=p1-p2;
  d=sqrt(u*u');
endfunction;
function D1_kernel=D1_int(t,k,p,dl,e,p0)
  //Function for integral over a boundary segment.
  //The midpoint of the segment is p, the segment edge 
  //vector is dl, the outward unit normal to the segment 
  //is e, the field point is p0. This function must not 
  //be used when the field point p0 taken on the segment 
  //itself, that case has to be done analytically because 
  //of the singularity in the Hankel function. The path is 
  //parameterised by t which goes from -1/2 to 1/2. The
  //wavenumber is k. Scilab's intg canned integration
  //routine only seems to do real integrals, so we have to
  //split the integral up into real and imaginary parts.
  //This does the real part.
  pt=p+t*dl;
  d=distance(pt,p0);
  u=direction(pt,p0);
  dircos=e*u';
  H21kd=besselh(1,2,k*d);
  D1_kernel=real(H21kd*dircos);
endfunction;
function D2_kernel=D2_int(t,k,p,dl,e,p0)
  //Function for integral over a boundary segment.
  //The midpoint of the segment is p, the segment edge 
  //vector is dl, the outward unit normal to the segment 
  //is e, the field point is p0. This function must not 
  //be used when the field point p0 taken on the segment 
  //itself, that case has to be done analytically because 
  //of the singularity in the Hankel function. The path is 
  //parameterised by t which goes from -1/2 to 1/2. The
  //wavenumber is k. Scilab's intg canned integration
  //routine only seems to do real integrals, so we have to
  //split the integral up into real and imaginary parts.
  //This does the imaginary part.
  pt=p+t*dl;
  d=distance(pt,p0);
  u=direction(pt,p0);
  dircos=e*u';
  H21kd=besselh(1,2,k*d);
  D2_kernel=imag(H21kd*dircos);
endfunction;
function [CNxM,DNxN,CMxM,DMxN]=port_matrices(k)
  //Calculates the port matrices for a given wavenumber k
  //Assumes that the geometry has been set up in the
  //arrays p,pe,e,dl.
  ea=1.0e-10;//abolute error for intg
  er=1.0e-3;//realtive error for intg 
  //Construct the CNxM matrix (boundary fields due to port currents)
  CNxM=zeros(N,M);
  for i=1:N,
    for l=1:M,
      H20kd=besselh(0,2,k*distance(p(i,:),pe(l,:)));
      CNxM(i,l)=-0.25*k*zeta0*H20kd;
    end;
  end;
  //Construct the DNxN matrix (boundary fields due to boundary fields)
  DNxN=zeros(N,N);
  for i=1:N,
    for j=1:N,
      if i==j then
        DNxN(j,j)=0.5;
      else
        Dl=sqrt(dl(j,:)*dl(j,:)');//Segment length
        if k<kcutoff then
          int1=intg(-0.5,0.5,list(D1_int,k,p(j,:),dl(j,:),e(j,:),p(i,:)),ea,er);
          int2=intg(-0.5,0.5,list(D2_int,k,p(j,:),dl(j,:),e(j,:),p(i,:)),ea,er);
          DNxN(i,j)=-0.25*%i*k*Dl*(int1+%i*int2);
        else
          d=distance(p(i,:),p(j,:));
          u=direction(p(i,:),p(j,:));
          dircos=e(j,:)*u';
          H21kd=besselh(1,2,k*d);
          DNxN(i,j)=0.25*%i*k*Dl*H21kd*dircos;
        end;
      end;
    end;
  end;
  //Construct the CMxM matrix (port fields due to port currents)
  CMxM=zeros(M,M);
  //Start by filling in the diagonal
  for m=1:M,
      CMxM(m,m)=-0.25*k*zeta0*(1-(2*%i/%pi)*(log(k*r0/2)+gEM-0.5));
  end;
  for m=1:M-1,
    for l=m+1:M,
      H20kd=besselh(0,2,k*distance(pe(l,:),pe(m,:)));
      CMxM(l,m)=-0.25*k*zeta0*H20kd;
      CMxM(m,l)=CMxM(l,m);//Fill in the symmetric matrix
    end;
  end;
  //Construct the DMxN matrix (port fields due to boundary fields)
  DMxN=zeros(M,N);
  for m=1:M,
    for j=1:N,
      Dl=sqrt(dl(j,:)*dl(j,:)');//Segment length
      u=direction(pe(m,:),p(j,:));
      dircos=e(j,:)*u';
      H21kd=besselh(1,2,k*distance(pe(m,:),p(j,:)));
      DMxN(m,j)=0.25*%i*k*Dl*H21kd*dircos;
    end;
  end;
endfunction;
function Ez=E_field(k,EzNx1,IeMx1,p0)
  //Calculates the interior E-field at a point p0
  //given the wavenumber k, the boundary fields EzNx1, 
  //and the exciter currents IeMx1. The boundary fields
  //EzNx1 are found by getting the CNxM and DNxN matrices (at the 
  //given k),
  //[CNxM,DNxN,CMxM,DMxN]=port_matrices(k);//Wasteful, don't use last pair
  //and then calculating the boundary fields EzNx1 as,
  //EzNx1=inv(eye(N,N)-DNxN)*CNxM*IeMx1;
  //
  //Solution for the E-field in the interior
  Ez=0;
  //Contribution from the port currents
  for l=1:M,
    d=distance(pe(l,:),p0);
    //If the field point p0 is on a port, we have to use the 
    //formula worked out in connection with the CMxM matrix that
    //integrates through the singularity.
    if d < r0 then
      Ez=Ez-0.25*k*zeta0*(1-(2*%i/%pi)*(log(k*r0/2)+gEM-0.5))*IeMx1(l);
    else 
      H20kd=besselh(0,2,k*d);
      Ez=Ez-0.25*k*zeta0*H20kd*IeMx1(l);
    end;
  end;
  //Contribution from the boundary currents
  //Note, we don't deal with the singularity when we get close to
  //the boundary
  for j=1:N,
    Dl=sqrt(dl(j,:)*dl(j,:)');//Segment length
    d=distance(p(j,:),p0);
    u=direction(p0,p(j,:));
    dircos=e(j,:)*u';
    H21kd=besselh(1,2,k*d);
    Ez=Ez+0.25*%i*k*Dl*H21kd*dircos*EzNx1(j);
  end;
endfunction;
function Z=standard_cap(w,C,Lesl,Resr)
  //Equivalent for a standard cap
  Z=Resr+%i*w*Lesl+1/(%i*w*C);
endfunction;
function Z=mlcc(w,Cp,Lesl,Re,R,Ri)
  //Equivalent circuit for a TDK MLCC (Multilayer
  //Ceramic Capacitor) from http://www.component.tdk.com
  //w=2*%pi*f
  //Cp=capacitance
  //Lesl=equivalent series inductance
  //Re=electrode resistance
  //R=value in Ohm.Hz from TDK tech notes
  //Ri=insulation resistanceRh=2*%pi*R/w;//Hysteris-derived resistance
  Z=Re+%i*w*Lesl+1/(1/Ri+1/Rh+%i*w*Cp);
endfunction;
//************ Useful constants ****************
//SI units
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
zeta0=sqrt(u0/(e0*er));//Free space impedance
gEM=-dlgamma(1);//Euler-Mascheroni constant (0.5772)
//Cut-off frequency. The integrals which compute the boundary
//E-field due to a boundary E-field need to be done by the intg
//canned routine at low k. The criterion for "low k" is when
//kd ~ 1 where d is the board length because then all effects
//are computed in the regime where the Hankel function changes 
//rapidly towards its singularity. Above the cut-off the Hankel
//function is better behaved and we use a crude approximation to 
//the integrals which reduces the time for the calculations.
wcutoff=2*%pi*fcutoff;
kcutoff=sqrt(u0*e0*er)*wcutoff;//Wavenumber
// ************ Start of the code *******************
if board_shape==circular_board then
  v=vertices_for_circular_board(N,r1);
elseif board_shape==rectangular_board then
  [v,N]=vertices_for_rectangular_board(Nx,Ny,a,b);
elseif board_shape==arbitrary_board then
  v=hand_crafted_vertices(N);// ?
else 
  error('No board specified');
end;
M=1;//Plus 1 for driving port
for s=1:length(sn),
  M=M+sn(s);//Get number of ports
end;
if board_shape==circular_board then
  pe=random_ports_for_circular_board(M,r1);
elseif board_shape==rectangular_board then
  pe=random_ports_for_rectangular_board(M,a,b);
elseif board_shape==arbitrary_board then
  pe=hand_crafted_ports(M);// ?
else
  error('No board specified');
end;

species=zeros(M,1);//Label each port with the species.
ic=1;//first port is driver port 
species(ic)=0;//Use 0 for the driver "species"
for s=1:length(sn),
  for sic=1:sn(s),
    ic=ic+1;//Increment port
    species(ic)=s;//This ic is species s
  end;
end;
//
//Boundary points p, edges dl, and outward unit normals e; this 
//section generates automatically from the vertices v
p=zeros(N,2);
dl=zeros(N,2);
e=zeros(N,2);
for i=1:N,
  if i==N then
    p(i,:)=(v(i,:)+v(1,:))/2;
    dl(i,:)=v(1,:)-v(i,:);
  else
    p(i,:)=(v(i,:)+v(i+1,:))/2;
    dl(i,:)=v(i+1,:)-v(i,:);
  end;
  e(i,:)=[dl(i,2),-dl(i,1)]/sqrt(dl(i,:)*dl(i,:)');
end;
if plot_boundary_and_ports then
  //Plot the shape of the pcb
  xbasc();
  printf('Plotting boundary and ports\n');
  printf('Type return to continue\n');
  for i=1:N,
    xv(1)=p(i,1)-dl(i,1)/2;
    xv(2)=p(i,1)+dl(i,1)/2;
    yv(1)=p(i,2)-dl(i,2)/2;
    yv(2)=p(i,2)+dl(i,2)/2;
    plot2d4(xv,yv,rect=board_display,frameflag=3,style=1);
  end;
  plot2d(pe(:,1)',pe(:,2)',style=-3,rect=board_display);//Plot ports
  pause;
  xbasc();//Clear plot of board ready for bode plot
end;
//E-field at a fixed frequency f
if calc_field then
  Ie=zeros(M,1);
  Ie(1)=1;//Port currents
  w=2*%pi*f;
  k=sqrt(u0*e0*er)*w;//Wavenumber
  [CkNxM,DkNxN,CkMxM,DkMxN]=port_matrices(k);//Wasteful, don't use last pair
  EzNx1=inv(eye(N,N)-DkNxN)*CkNxM*Ie;//Solve for the boundary fields
  phi=0;
  for j=1:length(rvec),
    p0=[rvec(j)*cos(phi),rvec(j)*sin(phi)];
    Ezvec(j)=E_field(k,EzNx1,Ie,p0);
  end;
end;
if calc_impedance then
  //Matrix solution for an M-port system of discrete components
  //connected to the board. This section needs to be hand-crafted 
  //to set up the particular mix of capacitors required.
  Zd=zeros(M,M);//The diagonal matrix of discrete impedances
  Zc=zeros(M,M);//The symmetric impedance matrix from the field solution
  Id=zeros(M,1);//Row vector of currents into each port.
  Vd=zeros(M,1);//Port voltages
  ed=zeros(M,1);//Row vector of discrete exciter emfs
  //Exciter emfs
  ed(1)=1;
  n=no_of_data_points;//No of data points for impedance calc 
  fvec=zeros(1,n);//Frequency vector (for Bode plot)
  Zvec=zeros(1,n);//Impedance looking into port 1
  Z11vec=zeros(1,n);//Impedance looking into port 1 for bare board
  m=no_of_decades;//No of decades f0 to f0*10^m
  for j=1:n,//Loop on frequency
    fvec(j)=f0*10^(m*(j-1)/(n-1));
    w=2*%pi*fvec(j);
    k=sqrt(u0*e0*er)*w;//Wavenumber
    [CkNxM,DkNxN,CkMxM,DkMxN]=port_matrices(k);//Using all matrices now
    Zc=-h*(CkMxM+DkMxN*inv(eye(N,N)-DkNxN)*CkNxM);//Impedance matrix
    Z11vec(j)=Zc(1,1);//Bare board impedance into port 1 is Zc(1,1)
    //Calculate the impedance matrix of the discretes
    Zd(1,1)=0;//Port 1 is the exciter port
    for i=2:M,
      is=species(i);//Get the species of cap at the ith port
      //Modify this if statement to explicitly call other cap models ?
      if is<=9 then
        Zd(i,i)=standard_cap(w,C(is),Lesl(is),Resr(is));
      else 
        error('This should not happen');
      end
    end;
    //Solve for the currents of the M-port system 
    Id=inv(Zd+Zc)*ed;
    Vd=Zc*Id;//Get the port voltages
    Zvec(j)=Vd(1)/Id(1);//Impedance at port 1 for each frequency
    printf('%s%f%s%+f%+fi\n','Zvec(',fvec(j),')=',real(Zvec(j)),imag(Zvec(j)));
  end;
  if plot_power_plane_impedance then
    bode(fvec,Zvec);
  end;
  if plot_bare_board_impedance then
    bode(fvec,Z11vec);
  end;
end;