###############################################################################################
# quarter_wave_monopole.lsf
# 
# Calculates the return loss, directivity pattern, maximum directivity, gain and radiation 
# efficiency of a coaxial fed quarter wave monopole antenna on top of an infinite ground plane.
# The script also plots the near field at the resonant frequency of the antenna and shows that
# all the power going through the directivity monitor box can be accounted for with sufficient 
# accuracy to perform a useful far field projection. Finally, a summary of the radiation properties
# of the antenna is provided on the script prompt.
#
# Copyright 2017 Lumerical Solutions Inc.
###############################################################################################

clear;
closeall;

####################
# Port 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); # injected power (input minus reflected)

# total power without mode expansion
port1T = getresult("::model::FDTD::ports::port 1","T");
port1T = abs(port1T.T); # recorded by the port DFT

# plot
f = port1Exp.f;
fGHz = f*1.0e-9;
plot(fGHz,port1TIn,port1TOut,port1TTot,port1T,"frequency (GHz)","power","Port Power");
legend("T_in","T_out","T_total","T");
setplot("x min",1.0);
setplot("x max",3.5);
setplot("y max",1.05);
setplot("y min",-0.05);




#############
# Return Loss
#############
#S11=port1Exp.S;
#plot(fGHz,20*log10(abs(S11)),"frequency (GHz)"," (dB)","Quater-Wave Monopole -S11");

retLoss = 10.0*log10(port1TOut);
plot(fGHz,retLoss,"frequency (GHz)","return loss (dB)","Return Loss");
setplot("show legend",false);
setplot("x min",1.0);
setplot("x max",3.5);
setplot("y max",2.0);
setplot("y min",-40.0);

#################
# Near Field Plot
#################

fMinIndex = find(port1TOut,min(port1TOut));
# XZ plot
E2 = getelectric("::model::XZ",2);
x = getdata("::model::XZ","x")*1.0e3;
y = getdata("::model::XZ","y")*1.0e3;
z = getdata("::model::XZ","z")*1.0e3;
E2 = pinch(E2(1:length(x),1:length(y),1:length(z),fMinIndex));
E2 = 10.0*log10(E2);
image(x,z,E2,"x (mm)","z (mm)","|E|^2 at f = "+num2str(round(fGHz(fMinIndex)*1e2)/1e2)+" GHz");
setplot("x min",-75.0);
setplot("x max",75.0);
setplot("y min",-65.0);
setplot("y max",85.0);
setplot("colorbar min",-40.0);
setplot("colorbar max",0.0);

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

# power going through the directivity box
monZ2 = getresult("::model::directivity::z2","T");
TZMax = abs(monZ2.T);
monX2 = getresult("::model::directivity::x2","T");
TXMax = abs(monX2.T);
monY2 = getresult("::model::directivity::y2","T");
TYMax = abs(monY2.T);
TBox = abs(TZMax)+2.0*abs(TXMax)+2.0*abs(TYMax);
# power going through the coax feed
TFeed = port1T;
# power plot
plot(fGHz,TBox,TFeed,"frequency (GHz)","power","Near-Field Power");
legend("directivity box","coaxial feed");
setplot("x min",1.0);
setplot("x max",3.5);
setplot("y max",1.1);
setplot("y min",-0.1);

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

# calculate directivity from the far field
select("directivity");
set("frequency",f(fMinIndex));
runanalysis;
dirMon = getresult("directivity","Directivity");
theta = dirMon.theta;
phi = dirMon.phi; 
nt = length(theta);
np = length(phi);
p1 = find(phi,0.0);
t1 = find(theta,pi/2.0);
Dtheta = 10.0*log10(dirMon.Dtheta);
Dphi = 10.0*log10(dirMon.Dphi);
Dmax = max(Dtheta); 
Dtheory = 10.0*log10((cos(pi/2.0*cos(theta))/sin(theta))^2.0)+5.167; # theoretical directivity
# plot directivity on elevation plane
plot([-flip(theta(2:nt),1);theta]*180.0/pi,
     [flip(Dtheory(2:nt),1);Dtheory(1:nt)],
     [flip(Dtheta(p1,2:nt),2),Dtheta(p1,1:nt)],
     "theta (deg)","D(theta) (dB)","Directivity on Elevation Plane");
legend("theory","simulation");
setplot("x min",-90.0);
setplot("x max",90.0);
setplot("y max",7.0);
setplot("y min",-25.0);

# plot directivity on azimuth plane
plot(phi*180.0/pi,
     Dtheory(t1)+matrix(np,1),
     Dtheta(1:np,t1),
     "phi (deg)","D(phi) (dB)","Directivity on Azimuth Plane");
legend("theory","simulation");
setplot("x min",0.0);
setplot("x max",360.0);
setplot("y max",7.0);
setplot("y min",-25.0);
# combined polar plot
Dtheta = Dtheta+(20.0-max(Dtheta));  # directivitiy normalized such that 20dB is max direct
# filter out results below 0dB.
index = find(Dtheta>0.0);
filter = matrix(np,nt);
filter(index) = 1.0;
Dtheta = Dtheta*filter; # normalized and filtered total directivitiy
thetaplot = [-flip(theta(2:nt),1);theta];
Delev=[transpose(flip(Dtheta(1,2:nt),2));transpose(Dtheta(1,1:nt))];
polar2(thetaplot,Delev,phi,Dtheta(1:np,t1),"","","Directivity");
legend("elevation plane","azimuth plane");
setplot("legend position",5);

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

Pin = port1TIn(fMinIndex)*sourcepower(f(fMinIndex)); # input power
Pacc = port1TTot(fMinIndex)*sourcepower(f(fMinIndex)); # accepted power
Prad = getresult("directivity","Prad"); # radiated power
#Prad = getresult("directivity","Prad")*TFeed(fMinIndex)/TBox(fMinIndex); # correction for radiated power

radEffIn = Prad/Pin; # radiated efficiency from input power
radEffAcc = Prad/Pacc; # radiated efficiency from accepted power
?msg = "============Radiation Performance==============";
?msg = "Resonant Frequency: "+num2str(round(fGHz(fMinIndex)*100.0)/100.0)+" GHz";
?msg = "Input Power: "+num2str(round(Pin*1.0e11)/1.0e2)+" nW";
?msg = "Accepted Power: "+num2str(round(Pacc*1.0e11)/1.0e2)+" nW";
?msg = "Radiated Power: "+num2str(round(Prad*1.0e11)/1.0e2)+" nW";
?msg = "Radiation Efficiency from Input Power: "+num2str(round(radEffIn*1.0e3)/10.0) + " %"; 
?msg = "Radiation Efficiency from Accepted Power: "+num2str(round(radEffAcc*1.0e3)/10.0) + " %";
?msg = "Maximum Directivity: "+num2str(round(Dmax*100.0)/100.0)+" dB ";
?msg = "Total Realized Gain: "+num2str(round((Dmax+10.0*log10(radEffIn))*100.0)/100.0)+" dB ";
?msg = "=======================================";