%SIMPLE MATLAB CODE FOR ILLUSRATIONS INCLUDED IN THE PAPER:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
"Electrostatics of crossed arrays of strips" by E.J.Danicki
% LaTeX notations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The applied notations generally agree with the text of paper. 
The differences concern the strip parameters $K$ and $\Delta$ 
which can be different on both sides of the layer: 
{\tt Ku, Kb} are $K_u, K_b$ and {\tt cu, cb} are $\Delta_u, 
\Delta_b$, respectively; the function $P_v(x)$ is {\tt pvx(v,x)}. 
The field spectra names start with the layer side initials, 
correspondingly {\tt u} and {\tt b} for upper and bottom side, 
followed by the names of unknowns: {\tt a} for $\alpha$ 
and {\tt b} for $\beta$.
{\small\begin{verbatim}
%%%% CROSSED STRIPS %%%%
d=.25;dx=2*pi*d;Ku=1;Kb=1;cu=-.9;cb=.9;s=.0001;L=256;
N=4;M=4;NN=2*N+1;MM=2*M+1;qr=zeros(L,NN);o=1e-4;
for i=1:L;r=Ku*(i-1)/L+o;
ua=zeros((MM-1)*NN);ub=zeros((MM-1)*NN,(NN-1)*MM);
ba=zeros((NN-1)*MM,(MM-1)*NN);bb=zeros((NN-1)*MM);
for n=-N:N-1;rn=r+n*Ku;nN=n+N;for m=-M:M-1;mM=m+M;sm=s+m*Kb;
k=sqrt(rn^2+sm^2);ct=coth(k*dx);sh=sinh(k*dx);
for n1=-N:N;nn=n1+N;
ua(mM*NN+NN,mM*NN+nn+1)=(-1)^n1*pvx(-n1-r/Ku,-cu);end;
for n1=-N:N;nn=n1+N;
ua(mM*NN+nN+1,mM*NN+nn+1)=(sign(n-n1+.5)-ct*rn/k)*pvx(n-n1,cu);end;
for m1=-M:M;mm=m1+M;
ub(mM*NN+nN+1,nN*MM+mm+1)=-rn/(k*sh)*pvx(m-m1,cb);end;
for m1=-M:M;mm=m1+M;
bb(nN*MM+MM,nN*MM+mm+1)=(-1)^m1*pvx(-m1-s/Kb,-cb);end;
for m1=-M:M;mm=m1+M;
bb(nN*MM+mM+1,nN*MM+mm+1)=(sign(m-m1+.5)-ct*sm/k)*pvx(m-m1,cb);end;
for n1=-N:N;nn=n1+N;
ba(nN*MM+mM+1,mM*NN+nn+1)=-sm/(k*sh)*pvx(n-n1,cu);end; end;end;
X=[ua ub;ba bb];ia=size(ua);ia=ia(1);ib=size(bb);ib=ib(1);
x=X\[zeros(M*NN+NN-1,1);sin(pi*r/Ku)*Ku/pi;zeros(ia+ib-(M+1)*NN,1)];
for n=-N:N-1;q=0;
for m1=-M:M;mm=m1+M;q=q+x(ia+(n+N)*MM+mm+1)*pvx(-m1-s/Kb,cb);end;
qr(i,n+N+1)=q;end;o=0;end;q1=real(reshape((qr),L*NN,1));
q1f=real(fftshift(ifft(fftshift(q1(1:8*L)))));
subplot(2,1,1);plot([-4*L:4*L-1]'/L,q1(1:8*L),'k');
subplot(2,1,2);plot([-4*L:4*L-1]'/8,8*q1f,'k');
%%%% LEGENDRE FUNCTION %%%%
function f=pvx(v,x);
z=(1-x)/2;a=1;b=1;nn=100;while nn<1/z;nn=100+nn;end;nn=2*nn;f=1;
if abs(z)<1;for n=1:nn;a=z*a*(1-(1+v)/n);b=b*(1+v/n);f=f+a*b;end;end;
%%%% STRIPS ON METALIZED LAYER %%%%
d=.3;dx=2*pi*d;K=1;N=1+round(1/(K*d));NN=2*N+2;u=zeros(NN);M=2^8;s=.0001;
x=0;qr=zeros(M,1);for l=1:M;r=1e-6+(l-1)*K/M;for n=-N:N;nN=n+N+1;rn=r+n*K;
t=tanh(abs(rn)*dx); for n1=-N:N+2;nn=n1+N+1;u(nN,nn)=(sign(n-n1+.5)-t*...
rn/sqrt(rn^2+s^2))*pvx(n-n1,x);end;end;for n1=-N:N+2;nn=n1+N+1;
u(NN,nn)=(-1)^n1*pvx(-n1-r/K,-x);end;b=[zeros(NN-1,1);2*sin(pi*r/K)];
a=u\b;for n1=-N:N+2;nn=n1+N+1;b(nn)=pvx(-n1-r/K,x);end;qr(l)=b.'*a;end;
Q=ifft(qr);plot([-M/2:M/2-1]',real(fftshift(Q)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% IMPROVED CODE FOR s==0 %%%%%%%%%%%%%%%%%%% plots charge on bottom strip
d=.5;dx=2*pi*d;Ku=1;Kb=1;cu=0;cb=0;L=256;o=1e-4;                 N=3;M=3;NN=2*N+1;MM=M+1;
qr=zeros(L+1,NN-1);qP=zeros(1,2*M+2);for m1=-M:M+1;mm=m1+M;qP(m1+M+1)=pvx(-m1,cb);end;
ia=MM*NN;ib=(NN-1)*MM;ua=zeros(ia);ub=zeros(ia,ib);ba=zeros(ib,ia);bb=zeros(ib);
for i=1:L/2+1;r=Ku*(i-1)/L+o;o=0;for n=-N:N-1;rn=r+n*Ku;nN=n+N; for m=-M:0;mM=m+M;sm=m*Kb;
k=sqrt(rn^2+sm^2);th=tanh(k*dx);ch=cosh(k*dx);                %m'=(-M,0);m=(-M,0);
for n1=-N:N;nn=n1+N;ua(mM*NN+NN,mM*NN+nn+1)=(-1)^n1*pvx(-n1-r/Ku,-cu);end;
for n1=-N:N;nn=n1+N;ua(mM*NN+nN+1,mM*NN+nn+1)=(th*sign(n-n1+.5)-rn/k)*pvx(n-n1,cu);end;
for m1=-M:0;mm=m1+M;ub(mM*NN+nN+1,nN*MM+mm+1)=-rn/(k*ch)*(pvx(m-m1,cb)+pvx(-m-m1,cb));end;
for m1=-M:0;mm=m1+M;if m~=0;    bb(nN*MM+mM+1,nN*MM+mm+1)=th*(sign(m-m1+.5)* ...
         pvx(m-m1,cb)+sign(m+m1-.5)*pvx(-m-m1,cb))-(pvx(m-m1,cb)+pvx(-m-m1,cb))*sm/k;
         else;bb(nN*MM+mM+1,nN*MM+mm+1)=(-1)^m1*th* ...
         2e6*(pvx(-m1+.5e-6,-cb)-pvx(-m1-.5e-6,-cb))/Kb-2*pvx(-m1,cb)/k;end;end;
for n1=-N:N;nn=n1+N;ba(nN*MM+mM+1,mM*NN+nn+1)=-pvx(n-n1,cu)/(k*ch);     if m~=0; ...
ba(nN*MM+mM+1,mM*NN+nn+1)=sm*ba(nN*MM+mM+1,mM*NN+nn+1);end;end;end;end;X=[ua ub;ba bb];
x=X\[zeros((M+1)*NN-1,1);(1+exp(8i*pi*r/Ku))*sin(pi*r/Ku)*Ku/pi;zeros(ib,1)];%2 strips!
b=x(ia+1:ia+ib);b=x(ia+1:ia+ib);b=reshape(b,M+1,NN-1);b=[b;flipud(b)];qr(i,:)=qP*b;
qr(L+2-i,:)=qP*flipud(fliplr(conj(b)));end;   q=(reshape(qr(1:L,:),L*(NN-1),1));
plot([-4*L:4*L-1]'/8,real(fftshift(ifft(fftshift([zeros(L,1);q;zeros(L,1)])))));
\end{verbatim}}