clear
 format compact
% program pulbf.m   (rev 2/99)
% scattering from a wire using the method of moments with pulse
% basis functions.  wire is along the z axis. TM pol assumed at an 
% angle theta.  Nseg is the number of segments -- assumed to be equal
% in length and the width is w (<< wavelength).
% Perfect electric conductor assumumed.   
% uses delta function approximation for Zmn
% 
      start=0; stop=90; dt=1;
      rad=pi/180;
      bk=2*pi;
      conk=4*pi;
      eta=377;
      sigc=bk^2/conk*eta^2;
      lwz=input('enter length: ');
      nseg=input('enter nseg: ');
	   Rs=input('resistivity: ');
% calculations use an equivalent strip width (radius*4) thus use 4*aw
      ae=lwz/2/74.2/2;      
      delw=lwz/nseg;
% wire patch center points, zn
     for k=1:nseg
      zn(k)=(-lwz/2+(2*k-1)*delw/2);
     end
%111111111111111111111 COMPUTE IMPEDANCE ELEMENTS 1111111111111111111
% compute impedance elements
    nint=5;
    dint=delw/nint;
    for m=1:nseg
     for n=1:nseg
      sum=0;
      zlo=zn(n)-delw/2;
        for i=1:nint
         z=zlo+dint/2+(i-1)*dint;
         r=sqrt((zn(m)-z)^2+ae^2);
         sum=sum+exp(-j*bk*r)/r;
        end
      gmn=delw*sum*dint/conk;
      r1=sqrt((zn(m)-zn(n))^2+ae^2);
      r2=sqrt((zn(n)-zn(m)-delw)^2+ae^2);
      r3=sqrt((zn(n)-zn(m)+delw)^2+ae^2);
      s1=exp(-j*bk*r1)/r1/conk;
      s2=exp(-j*bk*r2)/r2/conk;
      s3=exp(-j*bk*r3)/r3/conk;
      sum=j*(gmn-(2*s1-s2-s3)/bk^2);
      Z(m,n)=eta*bk*(gmn-(2*s1-s2-s3)/bk^2);
     end
    end
	ZR=eye(size(Z))*delw*Rs/2/pi/ae;
	Z=Z+ZR;
    disp('impedance elements computed')
    Zinv=inv(Z);
    disp('impedance matrix inverted')
%2222222222222222222222222 PATTERN LOOP 2222222222222222222222222222
      it=floor((stop-start)/dt)+1;
    for i=1:it           
      ang(i)=(i-1)*dt+start;
      thr=ang(i)*rad;
      ct=cos(thr);
      st=sin(thr);
      et=0;           
% plane wave excitation elements
     for n=1:nseg
      arg=bk*delw*ct/2;
      fn=1;
      if abs(arg) > 1e-5, fn=sin(arg)/arg; end
      rw(n)=st*delw*exp(j*bk*zn(n)*ct)*fn;
     end
% receive vector same as excitation vector
      B=rw;
% solve matrix equation for expansion coefficients
      C=Zinv*B.';
% scattered field
      et=rw*C;
      etm=abs(et)^2*sigc;
      sigdb(i)=10*log10(etm+1e-5);
%     disp([num2str(i),' of ',num2str(it),' angles completed'])
    end
   figure(1),clf
   plot(ang,sigdb),grid
   xlabel('Angle, Degrees')
   ylabel('Sigma/wavl^2, dB')
   title(['length in wavlengths= ',num2str(lwz),', number of segments= ',num2str(nseg)])