wiregrid.m

用矩量法计算细线天线的散射场

clear   
format compact
% program wiregrid.m   (rev 3/00)
% scattering from wires 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
%
      rad=pi/180;
      conk=4*pi;
      eta=377;
	   C0=3e8;
% data for wire endpoints -- noninteractive case (if interact=1)
      nwires=6; 
      P1(1,1:3)=[-.2 -.2 0]; P2(1,1:3)=[-.1 0 0]; nsegs(1)=8;  % left wing
 	   P1(2,1:3)=[-.2 .2 0];  P2(2,1:3)=[-.1 0 0]; nsegs(2)=8;  % right wing
 	   P1(3,1:3)=[-.55 0 0]; P2(3,1:3)=[0 0 0]; nsegs(3)=20;  % fuselage
 	   P1(4,1:3)=[-.55 -.1 0]; P2(4,1:3)=[-.5 0 0]; nsegs(4)=5; % left tail
 	   P1(5,1:3)=[-.55 .1 0]; P2(5,1:3)=[-.5 0 0]; nsegs(5)=5;  % right tail
      P1(6,1:3)=[-.50 0 0]; P2(6,1:3)=[-.55 0 .1]; nsegs(6)=5;  % vert tail
      Lm=0;
	  interact=1; 
	  if interact==0, nwires=input('enter nwires: '); end
	  kseg=0; 
    for i=1:nwires
	    if interact==0
			R1(1:3)=input(['enter [x,y,z] of end 1 of wire ',num2str(i),': ']);
			R2(1:3)=input(['enter [x,y,z] of end 2 of wire ',num2str(i),': ']);
            nsegs(i)=input(['enter number of segments for wire ',num2str(i),': ']');
		else
			R1(1:3)=P1(i,:);
			R2(1:3)=P2(i,:);
		end
            dR=R2-R1; dx=dR(1); dy=dR(2); dz=dR(3);
            L=norm(dR);
            delw=L/nsegs(i);
            Lm=max(L,Lm);
% wire direction cosines
            Lh=dR/L;
            thw=acos(dz/L); phw=atan2(dy,dx+1e-10);
            uu=sin(thw)*cos(phw);
            vv=sin(thw)*sin(phw);
            ww=cos(thw);
% wire segment center points
         for k=1:nsegs(i)
		  kseg=kseg+1;
          dt=(delw/2+(k-1)*delw);
		  dseg(kseg)=delw;
		  Lhat(kseg,:)=Lh;
		  U(kseg)=uu; V(kseg)=vv; W(kseg)=ww;
          x(kseg)=R1(1)+dt*uu;
          y(kseg)=R1(2)+dt*vv;
          z(kseg)=R1(3)+dt*ww;
	     end
     end
     disp(['maximum segment length is: ',num2str(Lm)])
	 tlsegs=kseg;
% calculations use an equivalent radius (same for all wires)
      aw=0.005;   
%	  aw=L/2/74.2; 
%000000000000000000000 PLOT THE WIRE 0000000000000000000000000000000
  figure(1)
  plot3(x,y,z,'o'),grid,xlabel('x'),ylabel('y'),zlabel('z')
  title('wire segment center points')
  Zmax=max([x,y,z]); Zmin=min([x,y,z]); Lim=max([Zmax,abs(Zmin)]);
  axis([-Lim,Lim,-Lim,Lim,-Lim,Lim]),axis square
%111111111111111111111 COMPUTE IMPEDANCE ELEMENTS 1111111111111111111
  fstart=input('start frequency in MHz: ');
  fstop=input('stop frequency in MHz: ');
  df=0; nf=1;
  if fstart~=fstop, df=input('frequency step in MHz: '); end
  if df~=0, nf=floor((fstop-fstart)/df)+1; end
  if nf>1
	  phd=input('phi (deg) for backscatter calculation: ');
	  thd=input('theta (deg) for backscatter calculation: ');
	  it=1; dang=0; start=thd;
  end
  if nf==1
      start=0; stop=360; dang=2;
      it=floor((stop-start)/dang)+1;
	  phd=0;
  end
  phr=phd*rad;
   for nn=1:nf
	frq=fstart+(nn-1)*df;
	F(nn)=frq;
	freq=frq*1e6;
	wave=C0/freq;
	bk=2*pi/wave;
	sigc=bk^2/conk*eta^2;
	disp(['calculating frequency number ',num2str(nn),'   of ',num2str(nf)])
% compute impedance elements
    nint=5;
    for m=1:tlsegs
     for n=1:tlsegs
	  dint=dseg(n)/nint;
      sum=0;
      hlo=-dseg(n)/2;
        for i=1:nint
	     h=hlo+dint/2+(i-1)*dint;
         xi=x(n)+h*U(n);
         yi=y(n)+h*V(n);
         zi=z(n)+h*W(n);
         r=sqrt((xi-x(m))^2+(yi-y(m))^2+(zi-z(m))^2+aw^2);
         sum=sum+exp(-j*bk*r)/r;
	    end
	  tt=U(m)*U(n)+V(m)*V(n)+W(m)*W(n);
      gmn=dseg(n)*sum*dint/conk;
	  r1=sqrt((x(n)-x(m))^2+(y(n)-y(m))^2+(z(n)-z(m))^2+aw^2);
	  r2=sqrt(sqrt((x(n)-x(m)-dseg(m)*U(m))^2+(y(n)-y(m)-dseg(m)*V(m))^2+...
	      (z(n)-z(m)-dseg(m)*W(m))^2)^2+aw^2);
	  r3=sqrt(sqrt((x(n)-x(m)+dseg(m)*U(m))^2+(y(n)-y(m)+dseg(m)*V(m))^2+...
	      (z(n)-z(m)+dseg(m)*W(m))^2)^2+aw^2);
      s1=exp(-j*bk*r1)/r1/conk;
      s2=exp(-j*bk*r2)/r2/conk;
      s3=exp(-j*bk*r3)/r3/conk;
      Z(m,n)=eta*bk*(gmn*tt-(2*s1-s2-s3)/bk^2);
     end
    end
    disp('impedance elements computed')
    Zinv=inv(Z);
    disp('impedance matrix inverted')
%2222222222222222222222222 PATTERN LOOP 2222222222222222222222222222

	  cp=cos(phr); sp=sin(phr);
    for i=1:it           
      A(i)=(i-1)*dang+start;
      thr=A(i)*rad;
      ct=cos(thr);
      st=sin(thr);
	  That=[ct*cp ct*sp -st];
      Phat=[-sp cp 0];
      et=0; ep=0;       
% plane wave excitation elements
      Ui=st*cp; Vi=st*sp; Wi=ct;
     for n=1:tlsegs   
       sum1=0; sum2=0;
	   Tt=dot(That,Lhat(n,:));
       Tp=dot(Phat,Lhat(n,:));
	   dint=dseg(n)/nint;
       for k=1:nint
            h=-dseg(n)/2+dint/2+(k-1)*dint;
            xi=x(n)+h*U(n);
            yi=y(n)+h*V(n);
            zi=z(n)+h*W(n);
            Exp=exp(j*bk*(xi*Ui+yi*Vi+zi*Wi));
            sum1=sum1+Exp*Tt;
            sum2=sum2+Exp*Tp;
       end
        rwt(n)=sum1*dint;
        rwp(n)=sum2*dint;
     end
% receive vector same as excitation vector
      B=rwt;
% solve matrix equation for expansion coefficients
      C=Zinv*B.';
% scattered field
      et=rwt*C;
      etm=abs(et)^2*sigc;
      if nf==1, sigdb(i)=10*log10(etm+1e-5); end
	  if nf>1, sigfdb(nn)=10*log10(etm+1e-5); F(nn)=frq; end
    end    % end of angle loop
   end     % end of frequency loop
      if nf==1
         figure(2)
         plot(A,sigdb),grid, Ax=axis; axis([start,stop,Ax(3),Ax(4)])
         xlabel('Angle, Degrees')
         ylabel('RCS, dBsm')
         title(['number of wires in model= ',num2str(nwires),...
          '     number of segments= ',num2str(tlsegs),'  phi= ',...
		  num2str(phd)])
	  end
	  if nf>1
         figure(2)
         plot(F,sigfdb),grid, Ax=axis; axis([fstart,fstop,Ax(3),Ax(4)])
         xlabel('Frequency, MHz')
         ylabel('RCS, dBsm')
         title(['number of wires in model= ',num2str(nwires),...
          '     number of segments= ',num2str(tlsegs),'  theta= ',num2str(thd)...
		  '     phi= ',num2str(phd)])
	  end