%program to analyze radiation pattern of a cruciform antenna on an aircraft
%lower fuselage & wing by reading in the CAD data
%Adapted from work by Balanis [1] and Mailloux [2]

freq=1*10^9;            %corresponds ot X band
lambda=3*10^8/freq;     %from v=f*lambda
dx=0.5*lambda;          %half wavelngth spacing
dy=0.5*lambda;

uChord=1*lambda;   %useable amound of wing chord in m
uSpan=5*lambda;   %usable amount of span on each wing from root outwards in m
totWidth=0.3048;  %total fuselage width in m
uLength=2*uSpan+totWidth;  %usable amount of fuselage length in m
uWidth=uChord;  %usable amount of fuselage width in m
%totSpan=1.828;    %total span of each wing from root outwards in m

%NOTE: spanwise out the wing is Y-direction, chordwise out the nose is X-direction

nxw=floor(uChord/dx);               %number of elements chordwise on each wing
nyw=floor(uSpan/dx);                %number of elements spanwise on each wing
nxr=floor(uChord/dx);               %number of elements lenghtwise in each root
nyr=floor((totWidth-uWidth)/(2*dx));%number of elements spanwise in each root
nyc=floor(uWidth/dx);               %number of elements chordwise on center section
nxc=floor(uLength/dx);              %number of elements spanwise on center section

Nw=nxw*nyw;
Nr=nxr*nyr;
Nc=nxc*nyc;
N=2*Nw+2*Nr+Nc;  %total number of elements

xmax=uLength/2;
ymax=(2*uSpan+totWidth)/2;          %max distances from center
k=2*pi/lambda;

xw=5*1i*xmax/(400*pi*1i-243)^0.5;       %element weighting variables
yw=5*1i*ymax/(400*pi*1i-243)^0.5;

curve=csvread('S600.csv',8,0);
phasing=csvread('R700.csv',8,0);

CADX=curve(:,3);
CADY=curve(:,2);
CADZ=curve(:,5);
%coordinate systems don't alight, re-assign in surface interpolation
F=scatteredInterpolant(CADZ, CADX, CADY);

%create interpolation ofr erroneous phasing
phaseX=phasing(:,3);
phaseY=phasing(:,2);
phaseZ=phasing(:,5);
Ph=scatteredInterpolant(phaseZ, phaseX, phaseY);

ePhase=zeros(N,1);
coord=zeros(N, 3);
n=1;

%nested loops to position elements in center of cross along fuselage
for row=1:nyc
    for col=1:nxc
        coord(n,1)=dx*(col-0.5*(nxc+1));
        coord(n,2)=dx*(row-0.5*(nyc+1));
        coord(n,3)=0;
        ePhase(n,1)=0;
        n=n+1;
    end
end

%nested loops to position elements in left root area
for row=1:nyr
    for col=1:nxr
        coord(n,1)=dx*(col-0.5*(nxr+1));
        coord(n,2)=dx*(row-0.5*(nyc+1)-nyr);
        coord(n,3)=0;
        ePhase(n,1)=0;
        n=n+1;
    end
end

%nested loops to position elements in right root area
for row=1:nyr
    for col=1:nxr
        coord(n,1)=dx*(col-0.5*(nxr+1));
        coord(n,2)=dx*(-row+0.5*(nyc+1)+nyr);
        coord(n,3)=0;
        ePhase(n,1)=0;
        n=n+1;
    end
end

%nested loops to position elements in left wing
for row=1:nyw
    for col=1:nxw
        coord(n,1)=dx*(col-0.5*(nxw+1));
        coord(n,2)=dx*(row-0.5*(nyc+1)-nyr-nyw);
        %coord(n,3)=F(coord(n, 1), -(coord(n, 2)+totWidth/2));
        coord(n,3)=0;
        %ePhase(n,1)=Ph(coord(n, 1), -(coord(n, 2)+totWidth/2));
        ePhase(n,1)=0;
        n=n+1;
    end
end

%nested loops to position elements in right wing
for row=1:nyw
    for col=1:nxw
        coord(n,1)=dx*(col-0.5*(nxw+1));
        coord(n,2)=dx*(+0.5*(nyc+1)+nyr+nyw-row);
        %coord(n,3)=F(coord(n, 1), (coord(n, 2)-totWidth/2));
        coord(n,3)=0;
        %ePhase(n,1)=Ph(coord(n, 1), (coord(n, 2)-totWidth/2));
        ePhase(n,1)=0;
        n=n+1;
    end
end

%element weighting with 24.3 dB gauss taper
a=zeros(N,1);
for b = 1:N
    a(b,1)=exp(-(coord(b,1)^2)/(2*xw^2))*exp(-(coord(b,2)^2)/(2*yw^2));
end

%calculate total power radiated [1] [2]
cc=zeros(N,1);
u=zeros(900,360);
v=zeros(900,360);
w=zeros(900,360);

for theta = 1:900
    for phi = 1:360
        u(theta,phi)=sind(theta/10)*cosd(phi);
        v(theta,phi)=sind(theta/10)*sind(phi);
        w(theta,phi)=cosd(theta/10);
    end
end

f=zeros(900,360);  %value of ration pattern electric field in theta and phi [1]

for theta = 1:900
    for phi = 1:360       
        for b = 1:N
              cc(b,1)=(k*u(theta,phi)*coord(b,1))+(k*v(theta,phi)*coord(b,2))+(k*w(theta,phi)*coord(b,3));
              f(theta,phi)=(f(theta,phi)+abs(a(b,1))*exp(-1i*k*ePhase(b,1))*exp(1i*cc(b,1)));        %phased array on
              %f(theta,phi)=(f(theta,phi)+(a(b,1))*exp(1i*cc(b,1)));                                 %phasing off
        end
    end
end

Ftot=zeros(900,360);
F1=zeros(900,360);

for theta=1:900
    for phi=1:360
        Ftot(theta,phi)=abs(f(theta,phi))^2;
        F1(theta,phi)=abs(f(theta,phi))^2*sind(theta/10);
    end
end

P0=(pi/900)*(pi/360)*sum(F1(:))/(2*pi);
P=10*log10(Ftot./P0);

%calculate power in 10 degree cone [2]
C=zeros(N,1);
U1=zeros(1000,360);
V1=zeros(1000,360);
W1=zeros(1000,360);

for theta = 1:1000
    for phi = 1:360
        U1(theta,phi)=sind(theta/100)*cosd(phi);
        V1(theta,phi)=sind(theta/100)*sind(phi);
        W1(theta,phi)=cosd(theta/100);
    end
end

f=zeros(1000,360);  %value of ration pattern electric field in theta and phi [2]

for theta = 1:1000
    for phi = 1:360       
        for b = 1:N
              C(b,1)=(k*U1(theta,phi)*coord(b,1))+(k*V1(theta,phi)*coord(b,2))+(k*W1(theta,phi)*coord(b,3));
              f(theta,phi)=(f(theta,phi)+abs(a(b,1))*exp(-1i*k*ePhase(b,1))*exp(1i*C(b,1)));        %phased array on
              %f(theta,phi)=(f(theta,phi)+(a(b,1))*exp(1i*cc(b,1)));                                 %phasing off
        end
    end
end

Fd=zeros(1000,360);

for theta=1:1000
    for phi=1:360
        Fd(theta,phi)=abs(f(theta,phi))^2;
        %F1(theta,phi)=abs(f(theta,phi))^2*sind(theta/10);
    end
end

Pb=10*log10(Fd./P0);       %total gain direcivity
P1=Pb-max(Pb(:));             %net directivity for finding angle

plotX=coord(:,1);
plotY=coord(:,2);
plotZ=coord(:,3);

%plots of beam from 0 to 90
axis=linspace(0,90,900);

figure(1)
plot(axis,P(:,45))
title('principal pattern cut at 45 deg')
xlabel('Angle \theta from perpendicular (degrees)')
ylabel('Gain (dBi)')
ylim([-50 40])

figure(2)
plot(axis,P(:,90))
title('principal pattern cut at 90 deg')
xlabel('Angle \theta from perpendicular (degrees)')
ylabel('Gain (dBi)')
ylim([-50 40])

figure(3)
plot(axis,P(:,180))
title('principal pattern cut at 0 deg')
xlabel('Angle \theta from perpendicular (degrees)')
ylabel('Gain (dBi)')
ylim([-50 40])

saveas(1,'45.jpg');
saveas(2,'90.jpg');
saveas(3,'0.jpg');

%plots of beam from 0 to 10
axis2=linspace(0,10,1000);

figure(4)
plot(axis2,Pb(:,45))
title('principal pattern cut at 45 deg')
xlabel('Angle \theta from perpendicular (degrees)')
ylabel('Gain (dBi)')
ylim([0 40])

figure(5)
plot(axis2,Pb(:,90))
title('principal pattern cut at 90 deg')
xlabel('Angle \theta from perpendicular (degrees)')
ylabel('Gain (dBi)')
ylim([0 40])

figure(6)
plot(axis2,Pb(:,180))
title('principal pattern cut at 0 deg')
xlabel('Angle \theta from perpendicular (degrees)')
ylabel('Gain (dBi)')
ylim([0 40])

saveas(4,'45 10.jpg');
saveas(5,'90 10.jpg');
saveas(6,'0 10.jpg');

%figure()
%surf(P)
%title('principal pattern cut')

%figure()
%scatter3(PX1, PY1, PZ1)

%figure()
%scatter3(plotX, plotY, plotZ)

Gain=max(Pb(:))      %max gain

deflectError=coord(:,3)-ePhase(:,1);

Width90=(sum(P1(:,90)>-3)+1)/50
Width0=(sum(P1(:,180)>-3)+1)/50
Width45=(sum(P1(:,45)>-3)+1)/50
maxError=max(max(deflectError), abs(min(deflectError)))/lambda*100