线阵分析-matlab

%%
%linear array
% number of elements to mainlobe and sidelobe  beamwidth, directivity
% different amplitude taper
% 1. W(n) = 1 ;
% 2. W(n) = abs((2*n-1)/N-1);
% 3. W(n) = abs(((2*n-1)/N)^2-1);
% 4. W(n) = power(cos(((2*n-1)/N-1)*pi/2),1)
% 5. W(n) = power(cos(((2*n-1)/N-1)*pi/2),2)
% 6. W(n) = power(cos(((2*n-1)/N-1)*pi/2),3)
% 7. W(n) = power(cos(((2*n-1)/N-1)*pi/2),4)
clc;
clear all;
close all;
N = 9; %
N = power(2,2:1:N);%elements number array
Ka =7 ; %antenna element amplitude taper conf
M = 2000; % theta sample points
AF = zeros(length(N),M);% array factor buffer
E0 = zeros(length(N),M);%
H0 = zeros(length(N),M);%
U0= zeros(length(N),M);%
f0 = 1e9;
lambda = 3e8/f0;
d = lambda/2;% elements spacing
k =2*pi/lambda;% wave constant
imp = 120*pi;%wave impdance
I0= 1;
r = 100;
l = 0.01;
BeamWidth = zeros(Ka,length(N));
D = zeros(Ka,length(N));
tau = zeros(Ka,length(N));
for m = 1:1:Ka
    figure;
    for i = 1:1:length(N)
        theta = linspace(0,pi,M)+pi/10e10;
        phi = linspace(0,2*pi,M);
        af(i,:) = zeros(1,M);
        a = 0;
        b = 0;
        for j = 1:1:N(i)
            switch(m)
                case 1
                    a= a+ 1;
                    b = b+1;
                    af(i,:) = af(i,:)+ exp(1i*(j-1)*k*d*cos(theta));                    
                  
                case 2
                    a= a+power(cos(((2*j-1)/N(i)-1)*pi/2),1);
                    b = b + power(cos(((2*j-1)/N(i)-1)*pi/2),2);
                    af(i,:) = af(i,:)+ exp(1i*(j-1)*k*d*cos(theta))*power(cos(((2*j-1)/N(i)-1)*pi/2),1);
                    
                case 3
                    a= a+power(cos(((2*j-1)/N(i)-1)*pi/2),2);
                    b = b + power(cos(((2*j-1)/N(i)-1)*pi/2),4);
                    af(i,:) = af(i,:)+ exp(1i*(j-1)*k*d*cos(theta))*power(cos(((2*j-1)/N(i)-1)*pi/2),2);
                    
                case 4
                    a= a+power(cos(((2*j-1)/N(i)-1)*pi/2),3);
                    b = b + power(cos(((2*j-1)/N(i)-1)*pi/2),6);
                    af(i,:) = af(i,:)+ exp(1i*(j-1)*k*d*cos(theta))*power(cos(((2*j-1)/N(i)-1)*pi/2),3);
                    
                case 5
                    a= a+power(cos(((2*j-1)/N(i)-1)*pi/2),4);
                    b = b + power(cos(((2*j-1)/N(i)-1)*pi/2),8);
                    af(i,:) = af(i,:)+ exp(1i*(j-1)*k*d*cos(theta))*power(cos(((2*j-1)/N(i)-1)*pi/2),4);
                case 6
                    a= a+power(0.5+0.5*cos(((2*j-1)/N(i)-1)*pi/2),1);
                    b = b + power(0.5+0.5*cos(((2*j-1)/N(i)-1)*pi/2),2);
                    af(i,:) = af(i,:)+ exp(1i*(j-1)*k*d*cos(theta))*power(0.5+0.5*cos(((2*j-1)/N(i)-1)*pi/2),1);
                    
                case 7
                    a= a+power(0.33+0.67*cos(((2*j-1)/N(i)-1)*pi/2),1);
                    b = b + power(0.33+0.67*cos(((2*j-1)/N(i)-1)*pi/2),2);
                    af(i,:) = af(i,:)+ exp(1i*(j-1)*k*d*cos(theta))*power(0.33+0.67*cos(((2*j-1)/N(i)-1)*pi/2),1);
                    
                    
            end
        end
        af(i,:) = abs(af(i,:)/N(i));
        E0(i,:) = af(i,:) *imp;
        Wav = real(E0(i,:) .* E0(i,:) )/2/imp;
        Prad = imp*pi*3*(I0*l/lambda)^2;
        theta = theta/pi*180;
        af(i,:)= 10*log10(af(i,:)/max(af(i,:)));
        U(i,:)= real(E0(i,:) .*E0(i,:))*r^2/2/imp;
        U(i,:)= 10*log10(U(i,:)/max(U(i,:)));
        plot(theta,U(i,:));
        hold on;
        axis([0 180 -120 0]);
        error = 0.01;
        af(i,:)=U(i,:);
        pos1_3dB = [];
        pos_max = find(max(af(i,:))==af(i,:));
        while(isempty(pos1_3dB))
            pos1_3dB = find(abs(((af(i,1:pos_max)-af(i,pos_max)))+3) < error);
            error = error + 0.005;
        end
        error = 0.01;
        pos2_3dB = [];
        while(isempty(pos2_3dB))
            pos2_3dB = find(abs(((af(i,pos_max:end)-af(i,pos_max)))+3) < error);
            error = error + 0.005;
        end
        BeamWidth(m,i)= (theta(pos2_3dB(1)+pos_max)-theta(pos1_3dB(end)));
        D(m,i) = 10*log10(a^2/b);
        tau(m,i) = a^2/b/N(i);
    end
    switch(m)
        case 1
%           title('\[{W_n} = 1\]1');
             text(80,5,'$$W_{n}= 1$$','interpreter','latex','fontsize',10);
% text(0.5,0.5,'$$f(x)=\frac{sin(x)}{x}$$','interpreter','latex','fontsize',20)
        case 2
           text(50,5,'$$W_{n}=cos\left [ \left ( \frac{2n-1}{N} -1\right )\frac{\pi}{2}  \right ]$$','interpreter','latex','fontsize',10)
         case 3
           text(50,5,'$$W_{n}=cos^{2}\left [ \left ( \frac{2n-1}{N} -1\right )\frac{\pi}{2}  \right ]$$','interpreter','latex','fontsize',10)
case 4
           text(50,5,'$$W_{n}=cos^{3}\left [ \left ( \frac{2n-1}{N} -1\right )\frac{\pi}{2}  \right ]$$','interpreter','latex','fontsize',10)
case 5
           text(50,5,'$$W_{n}=cos^{4}\left [ \left ( \frac{2n-1}{N} -1\right )\frac{\pi}{2}  \right ]$$','interpreter','latex','fontsize',10)
case 6
           text(50,5,'$$W_{n}=0.5+0.5cos\left [ \left ( \frac{2n-1}{N} -1\right )\frac{\pi}{2}  \right ]$$','interpreter','latex','fontsize',10)
case 7
           text(50,5,'$$W_{n}=0.33+0.67cos\left [ \left ( \frac{2n-1}{N} -1\right )\frac{\pi}{2}  \right ]$$','interpreter','latex','fontsize',10)
    end
end

说好的可读性呢?