######################################################################################
# leaky_wave.lsf
#
# Calculates the directivity pattern, max directivity, gain, and radiation efficiency
# of a leaky wave antenna. The script plots the near field and directivity pattern at
# the chosen frequencies. The script also shows that all the power going through the
# feed can be accounted for on the directivity monitor box with sufficient accuracy to 
# perform a meaningful far field projection. Lastly, a summary of the radiation properties 
# of the antenna is provided on the script prompt.
#
# Copyright 2017 Lumerical Solutions Inc
######################################################################################

clear;
closeall;

fPlot = [11.5;12.5;14.5]*1.0e9;

######################
# Port 1 Power Results
######################

# power on the selected port mode
port1Exp = getresult("::model::FDTD::ports::port 1","expansion for port monitor");
port1TIn = abs(port1Exp.T_in); # input power
port1TOut = abs(port1Exp.T_out); # reflected power
port1TTot = abs(port1Exp.T_total); # total power (input minus reflected)
# total power without mode expansion
port1T = getresult("::model::FDTD::ports::port 1","T");
port1T = abs(port1T.T);
# plot
f = port1Exp.f;
fGHz = f*1.0e-9;
plot(fGHz,port1TIn,port1TOut,port1TTot,port1T,"frequency (GHz)","power","Port 1 Power");
legend("T_in","T_out","T_total","T");
setplot("x min",11.0);
setplot("x max",15.0);
setplot("y max",1.05);
setplot("y min",-0.05);

######################
# Port 2 Power Results
######################

# power on the selected port mode
port2Exp = getresult("::model::FDTD::ports::port 2","expansion for port monitor");
port2TIn = abs(port2Exp.T_int); # input power
port2TOut = abs(port2Exp.T_out); # reflected power
port2TTot = abs(port2Exp.T_total); # total power (input minus reflected)
# total power without mode expansion
port2T = getresult("::model::FDTD::ports::port 2","T");
port2T = abs(port2T.T);
# plot
plot(fGHz,port2TIn,port2TOut,port2TTot,port2T,"frequency (GHz)","power","Port 2 Power");
legend("T_in","T_out","T_total","T");
setplot("x min",11.0);
setplot("x max",15.0);
setplot("y max",1.05);
setplot("y min",-0.05);

#######################
# Scattering Parameters
#######################

S11 = port1Exp.S;
port2Exp = getresult("::model::FDTD::ports::port 2","expansion for port monitor");
S21 = port2Exp.S;
plot(fGHz,20.0*log10(abs(S11)),20*log10(abs(S21)),"frequency (GHz)","amplitude (dB)","Scattering Parameters");
legend("S11","S21");
setplot("x min",11.0);
setplot("x max",15.0);
setplot("y max",2.0);
setplot("y min",-40.0);

##################
# Near Field Plots
##################

# frequencies to plot
f = getresult("::model::XY","f");
fGHz = f*1.0e-9;
fIndex = find(f,fPlot);
# XY plot
E = getresult("::model::XY","E");
H = getresult("::model::XY","H");
XYNearField = rectilineardataset("Near Field",E.x,E.y,E.z);
XYNearField.addparameter("lambda",c/f(fIndex),"f",f(fIndex));
XYNearField.addattribute("E", E.E(1:length(E.x),1:length(E.y),1:length(E.z),fIndex,1:3) );
XYNearField.addattribute("H", H.H(1:length(H.x),1:length(H.y),1:length(H.z),fIndex,1:3) );
visualize(XYNearField);

##################
# Near Field Power
##################

# power going through the box
monX1 = getresult("::model::directivity::x1","T");
monS1 = getresult("::model::directivity::s1","T");
TXMin = abs(monS1.T-monX1.T);
monX2 = getresult("::model::directivity::x2","T");
TXMax = abs(monX2.T);
monY1 = getresult("::model::directivity::y1","T");
TYMin = abs(monY1.T);
monY2 = getresult("::model::directivity::y2","T");
TYMax = abs(monY2.T);
monZ2 = getresult("::model::directivity::z2","T");
TZMax = abs(monZ2.T);
TBox = TXMin+TXMax+TYMin+TYMax+2.0*TZMax;
# power going through the feed
TFeed = abs(monX1.T);
# plot
plot(fGHz,TBox,TFeed,"frequency (GHz)","power","Near Field Power");
legend("directivity box","waveguide feed");
setplot("x min",11.0);
setplot("x max",15.0);
setplot("y max",1.1);
setplot("y min",0.75);

#############
# Directivity
#############

leg = cell(length(fIndex));
# sweep radiation patterns, directivity and radiated power
for(II = 1:length(fPlot)){
    select("directivity");
    set("frequency",f(fIndex(II)));
    runanalysis;
    dirMon = getresult("directivity","Directivity");
    theta = dirMon.theta;
    phi = dirMon.phi;
    nt = length(theta);
    np = length(phi);
    p1 = find(phi,0.0);
    p2 = find(phi,pi);
    phi = phi(p1:p2);
    t1 = find(theta,pi/2.0);#
    t1 = find(theta,80*pi/180.0);#
    theta = theta(t1);
    Dtheta = 10.0*log10(dirMon.Dtheta);
    Dphi = 10.0*log10(dirMon.Dphi);
    Dir = 10.0*log10(dirMon.Dphi+dirMon.Dtheta);
    leg{II} = num2str(round(fGHz(fIndex(II))*100.0)/100.0)+"GHz";
    PradII = getresult("directivity","Prad")*TFeed(fIndex(II))/TBox(fIndex(II));
    DmaxII = max(Dir);
    if(II == 1){
        D = Dir(p1:p2,t1);
        Prad = PradII;
        Dmax = DmaxII;
    }else{
        D = [D,Dir(p1:p2,t1)];
        Prad = [Prad;PradII];
        Dmax = [Dmax;DmaxII];
    }
}
# plot
plot(phi*180.0/pi,D,"phi (deg)","amplitude (dB)","Directivity on Azimuth Plane");
legend(leg);
setplot("x min",0.0);
setplot("x max",180.0);
setplot("y max",15.0);
setplot("y min",-5.0);

#####################################################
# Power, Efficiency and Directivity for Script Prompt
#####################################################

Pin = port1Exp.T_in(fIndex)*sourcepower(f(fIndex)); # input power
Pacc = port1Exp.T_total(fIndex)*sourcepower(f(fIndex)); # accepted power
radEffIn = Prad/Pin; # radiated efficiency from input power
radEffAcc = Prad/Pacc; # radiated efficiency from accepted power

?msg = "=====Radiation Performance=====";
?msg = "Frequency (GHz):\n "+num2str(round(fGHz(fIndex)*100.0)/100.0);
?msg = "Directivity (dB):\n "+num2str(round(Dmax*100.0)/100.0)+" dB ";
?msg = "Radiation Efficiency (%):\n "+num2str(round(radEffIn*1.0e3)/10.0);
?msg = "Total Realized Gain (dB):\n "+num2str(round((Dmax+10.0*log10(radEffIn))*1.0e2)/1.0e2);