function C=stripline(h1,eps1,h2,eps2,d,t,l,V0,V1,V2,N,M)
%
%   Calculates the Capacitance of a microstrip
%
%   C=stripline(h1,eps1,h2,eps2,d,t,l,V0,V1,V2,N,M)
%
if nargin<10
  V0=50;    V1=0;    V2=0;
  h1=2e-3; h2=2e-3;   
  d=1e-3;  l=5e-3;  t=10e-6;
  eps1=1;  eps2=10;
  N=60;    M=60;
end
% Define some handy constants
e0=8.8542e-12;
h=h1+h2;
dx=l/N;
dy=h/M;
dx2_inv=1/(2*dx);
dy2_inv=1/(2*dy);
dxsqr_inv=1/(dx^2);
dysqr_inv=1/(dy^2);
dx_inv=1/dx;
dy_inv=1/dy;
two_inv_dxsqr_dysqr=-2*(1/(dx^2)+1/(dy^2));
y1=floor((h1-t)/dy);
y2=ceil(h1/dy);
x1=floor((l-d)/2/dx);
x2=ceil((l+d)/2/dx);
% Display the simulated structure values
disp(['d=',num2str((x2-x1)*dx*1e3),'mm, t=',num2str((y2-y1)*dy*1e6),...
      'um, dx=',num2str(dx*1e6,3),'um, dy=',num2str(dy*1e6,3),...
      'um, h1=',num2str(h1*1e3),'mm, h2=',num2str(h2*1e3),'mm']);
disp(['y1=',num2str(y1),', y2=',num2str(y2),', x1=',num2str(x1),', x2=',num2str(x2)]);
disp('Setting up the equations...');

A=sparse([],[],[],N*M,N*M,0);
HL=zeros(N*M,1);

y=1
% Upper
for x=1:N
  k=(y-1)*N+x;
  A(k,k)=1;
  HL(k)=V1;
end
% Lower
y=
for 
  
end
% Left
x=
for   
   
end
% Right
x=
for 
   
end
for y=2:M-1
  for x=2:N-1
    k=(y-1)*N+x;
    % Contact
    if
       
    % Rand
 	 elseif     
    
    % Inner part
    else
       
    end
  end
end
disp('Solving...');
V=
disp('Plotting...');
V=reshape(V,N,M);
V=V';
% Calculate the electrical field
[Ex,Ey]=grad(-V,dx,dy);
% The electrical field at the contact must be treated separately
Ey(y2,x1:x2)=
Ey(y1,x1:x2)=
Ex(y1:y2,x1)=
Ex(y1:y2,x2)=

eps(1:y2-1,1:N)=eps1;
eps(y2:M,1:N)=eps2;

Dx=
Dy=

Q1=trapz(Dy(y1,x1:x2))*(-dx);
Q2=
Q3=
Q4=
Q=
C=
disp(['C=',num2str(C),' F/m']);

figure(11)
surf(V);
title('V');
xlabel x;ylabel y;
figure(12)
surf(Dy);
title('Dy');
xlabel x;ylabel y;
figure(13)
surf(Ex);
title('Ex');
xlabel x;ylabel y;
figure(14)
surf(Dx);
title('Dx');
xlabel x;ylabel y;
figure(15)
surf(Ey);
title('Ey');
xlabel x;ylabel y
figure(16)
y=y2-5:y2+5;
subplot(121);
plot(y,Ex(y,3:5),'x-')
grid;
title('Ex near the junction');
subplot(122);
plot(y,Dy(y,3:5),'+-');
grid;
title('Dy near the junction');
figure(17)
y=1:M;
subplot(121);
plot(y,Ex(y,3:5),'x',y,Ey(y,3:5),'+')
legend('Ex','Ey',2);
grid;
title('E-field near the junction');
subplot(122);
plot(y,Dx(y,3:5),'x',y,Dy(y,3:5),'+');
legend('Dx','Dy',2);
grid;
title('D-field near the junction');

function [Ex,Ey]=grad(V,dx,dy)
  dx2_inv=1/(2*dx);
  dy2_inv=1/(2*dy);
  dx_inv=1/dx;
  dy_inv=1/dy;
  Xres=size(V,2);
  Yres=size(V,1);
  for x=2:Xres-1
    for y=1:Yres
      Ex(y,x)=(V(y,x+1)-V(y,x-1)).*dx2_inv;
    end
  end
  for y=1:Yres
    Ex(y,1)=(V(y,2)-V(y,1)).*dx_inv;
    Ex(y,Xres)=(V(y,Xres)-V(y,Xres-1)).*dx_inv;
  end
  for y=2:Yres-1
    for x=1:Xres
      Ey(y,x)=(V(y+1,x)-V(y-1,x)).*dy2_inv;
    end
  end
  for x=1:Xres
    Ey(1,x)=(V(2,x)-V(1,x)).*dy_inv;
    Ey(Yres,x)=(V(Yres,x)-V(Yres-1,x)).*dy_inv;
  end
return






