%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^   
%***********************************************************************
%***          HALLEN'S INTEGRAL EQUATION                             ***
%***********************************************************************

case 2
%******************************************************
%     SOME CONSTANTS AND INPUT DATA
%******************************************************
nm=input('NUMBER OF SUBDIVISIONS <ODD OR EVEN NUMBER> =');
tl=input('TOTAL DIPOLE LENGTH <WAVELENGTHS> =');
ra=input('RADIUS OF DIPOLE <WAVELENGTHS> =');
n=nm;
l=tl;
rho=ra;
rtod=180/pi;
eta=120*pi;
bk=2*pi;
cj=0.0+j*1.0;
dz=l/(2*(n-1));
dz1=l/(2*n);
%C***********************************************************************
%C    PERFORM COMPLEX INTEGRATION USING SIXTEEN POINT 
%C    GAUSSIAN QUADRATURE WITH INCREASING ACCURACY SET 
%C    BY INTEGER NO
%C***********************************************************************
absica=[-0.095012509837637,-0.281603550779259;
        -0.458016777657227,-0.617876244402644;
        -0.755404408355003,-0.865631202387832;
        -0.944575023073233,-0.989400934991650;
         0.095012509837637,0.281603550779259;
         0.458016777657227,0.617876244402644;
         0.755404408355003,0.865631202387832;
         0.944575023073233,0.989400934991650];
wght=[0.189450610455068,0.182603415044924;
         0.169156519395002,0.149595988816577;
         0.124628971255534,0.095158511682493;
         0.062253523938648,0.027152459411754;
         0.189450610455068,0.182603415044924;
         0.169156519395002,0.149595988816577;
         0.124628971255534,0.095158511682493;
         0.062253523938648,0.027152459411754];
      %************************************************************
      %     FILL THE MATRIX AND EXCITATION VECTOR OF THE SYSTEM: 
      %      [ZMATRX] and [ELECUR]
      %************************************************************
for i=1:n
    z=(2*i-1)*dz1/2.0;
    zmatrx(i,n)=-cos(bk*z);
    elecur(i)=-cj*sin(bk*z)/(2.0*eta);
    for j=1:n-1
        lower=(j-1)*dz1;
        upper=(j)*dz1;
        %*****************************************************************
        %     PERFORM NUMERICAL INTEGRATION OF THE KERNEL FOR HALLEN'S
        %     INTEGRAL EQUATION
        %*****************************************************************
        no=10;
        del=(upper-lower)/(2.0*no);
        sum=0.0+j*0.0;
        for t=1:no
            s=lower+(2*t-1)*del;
            for q=1:16
               x=s+absica(q)*del;
               %***********************************************************************
               %C    KERNEL PROVIDES THE KERNEL OF HALLEN'S EQN FOR INTEGRATION 
               %C    SYMMETRY IS USED TO REDUCE THE SYSTEM OF EQUATIONS AND 
               %C    HENCE ALTERS THE KERNEL
               %***********************************************************************
                r1=sqrt(rho*rho+(z-x)*(z-x));
                r2=sqrt(rho*rho+(z+x)*(z+x));
                kernel=exp(-cj*bk*r1)/(4.0*pi*r1)+exp(-cj*bk*r2)/(4.0*pi*r2);
                sum=sum+wght(q)*kernel;
            end
        end
        res=sum*del;
    zmatrx(i,j)=res;
  end
end

%*********************************************************************************
%   DECOMPOSE AND SOLVE THE SYSTEM FOR THE CURRENT DISTRIBUTION
%*********************************************************************************

%****************************
%    GET SCALING INFO.
%****************************
for i=1:n
    zmax=0.0;
    for j=1:n
        caz=abs(zmatrx(i,j));
        if caz >= zmax
            zmax=caz;
        end
    end
    scal(i)=1.0/zmax;
 end
%**************************** 
%    CROUT's algorithm.
%****************************
for j=1:n
    for i=1:j-1
        for k=1:i-1
            zmatrx(i,j)=zmatrx(i,j)-zmatrx(i,k)*zmatrx(k,j);
        end
    end
    zmax=0.0;
    %*****************************************
    %    SEARCH FOR LARGEST PIVOT ELEMENT.
    %*****************************************
    for i=j:n
        for k=1:j-1
            zmatrx(i,j)=zmatrx(i,j)-zmatrx(i,k)*zmatrx(k,j);
        end
        problem=scal(i)*abs(zmatrx(i,j));
        if problem >=zmax
            imax=i;
            zmax=problem;
        end
     end
     %***********************************
     %    INTERCHANGE THE ROWS.
     %***********************************
    if j~=imax
        for k=1:n
            temp=zmatrx(imax,k);
            zmatrx(imax,k)=zmatrx(j,k);
            zmatrx(j,k)=temp;
        end
        scal(imax)=scal(j);
    end
    iperm(j)=imax;
    %***********************************
    %    DIVIDE BY PIVOT ELEMENT.
    %***********************************
    if j~=n
        for i=j+1:n
            zmatrx(i,j)=zmatrx(i,j)/zmatrx(j,j);
        end
    end
 end
%****************************************************************** 
%    SOLVES LINEAR SYSTEM GIVEN THE LU DECOMPOSITION FROM LUDEC
%    FORCING VECTOR IS REPLACED WITH SOLUTION VECTOR UPON EXIT
%*****************************
%    FORWARD SUBSTITUTION.
%******************************************************************
for i=1:n
    temp=elecur(iperm(i));
    elecur(iperm(i))=elecur(i);
    for j=1:i-1
        temp=temp-zmatrx(i,j)*elecur(j);
    end
    elecur(i)=temp;
 end
%******************************** 
%    BACKWARD SUBSTITUTION.
%********************************
for i=1:n
    ii=n-i+1;
    temp=elecur(ii);
    for j=ii+1:n
        temp=temp-zmatrx(ii,j)*elecur(j);
    end
    elecur(ii)=temp/zmatrx(ii,ii);
 end
%*************************************************************************** 
%     COMPUTATION OF THE INPUT IMPEDANCE AND CURRENT DISTRIBUTION
%***************************************************************************
zin=1.0/elecur(1);
%******************************************************************
%     Output files for curr-MoM_m.dat
%******************************************************************
fid=fopen('Curr-MoM_m.dat','wt');
fprintf(fid,'CURRENT DISTRIBUTION ALONG ONE HALF OF THE DIPOLE\n');
fprintf(fid,'POSITION Z:   CURRENT MAGNITUDE:    CURRENT PHASE:    SECTION:\n');
for i=1:n-1
   %******************************************************************
   %    THIS FUNCTION IS COMPUTES THE ARCTANGENT GIVEN X,Y.  IT IS 
   %    SIMILAR TO ATAN2 EXCEPT IT AVOIDS THE RUN TIME ERRORS ON 
   %    SOME MACHINES FOR SMALL ARGUMENTS.
   %******************************************************************
    smlt=0.0000001;
    if (abs(x)<=smlt) & (abs(y)<=smlt)
        btan2=0.0;
    else
        btan2=atan2(real(elecur(i)),imag(elecur(i)));
    end
    cur=abs(elecur(i));
    pha=rtod*btan2;
    table=[i*dz-dz/2.0,cur,pha,i];
    fprintf(fid,'%2.6f       %2.6f              %3.6f          %3.1f\n\n\n\n',table);
 end
%******************************************************************
%     Output figure for curr-MoM_m.dat
%******************************************************************

i=[1:n-1];
smlt=0.0000001;
    if (abs(x)<=smlt) & (abs(y)<=smlt)
        btan2=0.0;
    else
        btan2=atan2(real(elecur(i)),imag(elecur(i)));
    end
     cur=abs(elecur(i));
    pha=rtod*btan2;
stem(i*dz-dz/2.0,cur)
grid on
xlabel('Position (z)')
ylabel('Current magnitude')
title('CURRENT DISTRIBUTION ALONG ONE HALF OF THE DIPOLE')
fclose(fid);
figure;
%*********************************************************************************
%   CALCULATE THE RADIATION PATTERN OF THE ANTENNA
%*********************************************************************************
maxang=181;
pmax=-1.0;
for j=1:maxang
    theta=(j-1)/(maxang-1)*pi;
    %********************************************************************
    %    CALCULATES THE RADIATED POWER LEVEL AT ANGLE THETA RADIANS
    %    TO THE DIPOLE AXIS.  SINCE THE PATTERN IS NORMALIZED TO THE 
    %    MAXIMUM RADIATED POWER, COMMON CONSTANTS ARE REMOVED
    %********************************************************************
    sth=sin(theta);
    cth=cos(theta);
    arg=pi*dz*cth;
    if abs(arg) <= 0.001
        ft=1.0;
    else
        ft=sin(arg)/arg;
    end
    n2=n-1;;
    crt=0.0+j*0.0;
    for k=1:n2
        argp=pi*(-l+(k-1)*2.0*dz+dz)*cth;
        crt=crt+exp(cj*argp)*ft*elecur(k);
        argp=pi*(-l+(n2+k-1)*2.0*dz+dz)*cth;
        crt=crt+exp(cj*argp)*ft*elecur(k);
    end
    power=abs(crt)*sth*sth;
    pwr(j)=power;
    if pwr(j)>=pmax
        pmax=pwr(j);
    end
 end
 
%******************************************************************
%     Output files for Patt-MoM_m.dat
%******************************************************************
hel=fopen('Patt-MoM_m.dat','wt');
fprintf(hel,'RADIATION PATTERN   vs  OBSERVATION ANGLE THETA\n');
fprintf(hel,'THETA (in degrees)      MAGNITUDE (in dB)\n');
%*******************************************************************
%   WRITE THE RADIATION PATTERN IN dB
%*******************************************************************
for i=1:maxang
    theta=(i-1)/(maxang-1)*180;
    pattrn=pwr(i)/pmax;
    if pattrn <=0.00001
        pattrn=-100;
    else
        pattrn=20*log10(pattrn);
    end
    table1=[theta,pattrn];
    fprintf(hel,'%3.1f                         %3.3f         \n',table1);
end
%******************************************************************
%     Output figure for Patt-MoM_m.dat
%******************************************************************
i=[1:181];
theta=(i-1);
pattrn=pwr(i)/pmax;
pattrn=20*log10(pattrn);
% Polar Plot
pattrn=[pattrn,fliplr(pattrn(1:180))];
q=polar_dB([0:360],pattrn,-40,0,4,'-');
set(q,'linewidth',1.5);
% plot(theta,pattrn,'linewidth',2);
% axis ([0 180 -60 0]);
% grid on;
% xlabel('Theta(degrees)');
% ylabel('Magnitude(dB)');
title('RADIATION PATTERN   vs   OBSERVATION ANGLE')
for i=182:361
    theta=(i-1)/(maxang-1)*180;
    pattrn=pwr(362-i)/pmax;
    if pattrn <=0.00001
        pattrn=-100;
    else      
        pattrn=20*log10(pattrn);
    end
    table2=[theta,pattrn];
    fprintf(hel,'%3.1f                         %3.3f         \n',table2);
end
fclose(hel);
end
%*******************************************
%        OUTPUT FILE
%*******************************************
if(option_a==2)
   diary(filename);
end
%*******************************************
%         FORMAT STATEMENTS 
%*******************************************     
disp(sprintf('\n\n\n\n\n\n\n\n'));
disp(strvcat('WIRE ANTENNA PROBLEM','================'));
if (option==1)&(iex==1)
    disp(sprintf('POCKLINGTON''S EQUATION AND MAGNETIC FRILL MODEL'));
elseif (option==1)&(iex==2)
    disp(sprintf('POCKLINGTON''S EQUATION AND DELTA GAP MODEL'));
elseif (option==2)
    disp(sprintf('HALLEN''S EQUATION'));
end
disp(sprintf('\nLENGTH = %4.4f (WLS)',tl));
disp(sprintf('RADIUS OF THE WIRE = %4.4f (WLS)',ra));
disp(sprintf('NUMBER OF SUBSECTIONS = %2.2f\n',nm));
ZR=real(zin);
ZI=imag(zin);
if ZI>0
    output=sprintf('INPUT IMPEDANCE: Z= %5.1f +j %5.1f  (OHMS)\n',ZR,ZI);
else
    output=sprintf('INPUT IMPEDANCE: Z= %5.1f -j %5.1f  (OHMS)\n',ZR,abs(ZI));
end
disp(output);
disp(strvcat('*** NOTE:',...
   '    THE DIPOLE CURRENT DISTRIBUTION IS STORED IN Curr-MoM_m.dat',...
   '    THE AMPLITUDE RADIATION PATTERN IS STORED IN Patt-MoM_m.dat',...
   '    ========================================================='));
diary off;
warning off;