#########################################################
# PatchAntenna.lsf
#
# Uses FDTD to calculates the return loss, directivity 
# pattern, max directivity, gain, and radiation efficiency
# of a coaxial-fed monopole antenna ontop of an infinite 
# ground plane. The script generates the plots and radiation 
# properties in the script prompt.
#
# Copyright 2016 Lumerical Solutions Inc
##########################################################

clear;
closeall;

Le = 10.6822e-3; # effective length of patch
sub_h = getnamed("::model::antenna::Substrate","z span"); # substrate height
Wpatch = getnamed("::model::antenna::Patch","y span"); # width of patch

####################
# 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);
# plot
f = port1Exp.f;
fGHz = f*1.0e-9;
plot(fGHz,port1TIn,port1TOut,port1TTot,port1T,"frequency (GHz)","normalized power","Port Power");
legend("T_in","T_out","T_total","T");
setplot("x min",8.0);
setplot("x max",12.0);
setplot("y max",1.05);
setplot("y min",-0.05);

#############
# Return Loss
#############

retLoss = 10.0*log10(port1TOut);
plot(fGHz,retLoss,"frequency (GHz)","return loss (dB)","Return Loss from Reflected Power");
setplot("show legend",false);
setplot("x min",8.0);
setplot("x max",12.0);
setplot("y max",2.0);
setplot("y min",-20.0);

#################
# Input Impedance
#################

S11 = port1Exp.S;
rin = getnamed("::model::antenna::Pin","radius");
rout = getnamed("::model::antenna::Coax Outer Conductor","inner radius");
port1Neff = getresult("::model::FDTD::ports::port 1","neff");
neff = port1Neff.neff;
neff = interp(neff,port1Neff.f,f);
Z0 = sqrt(mu0/eps0)/neff/(2.0*pi)*log(rout/rin); # theoretical coax feed Z0
Zin = Z0*(1+S11)/(1-S11);
plot(fGHz,Zin,"frequency (GHz)","Zin","Input Impedance");
setplot("x min",8.0);
setplot("x max",12.0);
setplot("y1 max",45.0);
setplot("y1 min",-20.0);
setplot("y2 max",45.0);
setplot("y2 min",-20.0);
setplot("show legend",false);

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

# power going through the directivity box
monZ2 = getresult("::model::directivity::z2","T");
TZMax = abs(monZ2.T);
monX1 = getresult("::model::directivity::x1","T");
TXMin = abs(monX1.T);
monX2 = getresult("::model::directivity::x2","T");
TXMax = abs(monX2.T);
monY2 = getresult("::model::directivity::y2","T");
TYMax = abs(monY2.T);
TBox = TZMax+TXMax+TXMin+2.0*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",8.0);
setplot("x max",12.0);
setplot("y max",1.1);
setplot("y min",-0.1);

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

fMinIndex = find(port1TOut,min(port1TOut));
f0 = f(fMinIndex);
k0 = 2.0*pi*f0/c;
setnamed("::model::directivity","frequency",f0);
runanalysis;
dirMon = getresult("directivity","Directivity");
Dtheta = 10.0*log10(dirMon.Dtheta);
Dphi = 10.0*log10(dirMon.Dphi);
Dir = 10.0*log10(dirMon.Dphi+dirMon.Dtheta);
Dmax = max(Dir);
Phi1 = find(dirMon.phi,0.0);
Phi2 = find(dirMon.phi,pi);
Phi3 = find(dirMon.phi,pi/2.0);
Phi4 = find(dirMon.phi,3.0*pi/2.0);
nt = length(dirMon.theta);
np = length(dirMon.phi);
# E-plane plot
theta = [-flip(dirMon.theta(2:nt),1);dirMon.theta];
Eplane = sin(k0*sub_h*cos(theta)/2.0)*cos(k0*Le*sin(theta)/2.0)/(k0*sub_h*cos(theta)/2.0);
Hplane = sin(theta+pi/2.0)*sin(k0*sub_h*sin(theta+pi/2.0)/2.0)*sin(k0*Wpatch*cos(theta+pi/2.0)/2.0)/(k0*sub_h*sin(theta+pi/2.0)/2.0*k0*Wpatch*cos(theta+pi/2.0)/2.0);
plot(theta*180.0/pi,
     [flip(Dtheta(Phi2,2:nt),2),Dtheta(Phi1,1:nt)],
     [flip(Dphi(Phi2,2:nt),2),Dphi(Phi1,1:nt)],
     20.0*log10(abs(Eplane))+6.78118,
     "angle (deg)","D (dB)","Directivity: E-Plane (X-Z Cut)");
legend("D(theta) - simulation", "D(phi)- simulation","D(theta) - theory");
setplot("x max",90.0);
setplot("x min",-90.0);
setplot("y min",-20.0);
setplot("y max",10.0);
# H-plane plot
plot(theta*180.0/pi,
     [flip(Dtheta(Phi4,2:nt),2),Dtheta(Phi3,1:nt)],
     [flip(Dphi(Phi4,2:nt),2),Dphi(Phi3,1:nt)],
     20.0*log10(abs(Hplane))+6.7811,
     "angle (deg)","D (dB)","Directivity: H-Plane (Y-Z Cut)");
legend("D(theta) - simulation","D(phi) - simulation", "D(phi) - theory");
setplot("x max",90.0);
setplot("x min",-90.0);
setplot("y min",-20.0);
setplot("y max",10.0);
# polar plot
Dthetap = Dtheta+(20.0-max(Dtheta));  # total directivitiy normalized such that 20dB is max direct
Dphip = Dphi+(20.0-max(Dphi));
# filter out results below 0dB.
index = find(real(Dthetap) > 0.0);
filter = matrix(np,nt);
filter(index) = 1;
Dthetap = Dthetap*filter; # normalized and filtered total directivitity
index = find(real(Dphip) > 0.0);
filter = matrix(np,nt);
filter(index) = 1;
Dphip = Dphip*filter; #normalized and filtered total directivitity
D_el = [transpose(flip(Dthetap(Phi2,2:nt),2));transpose(Dthetap(Phi1,1:nt))];
D_el2 = [transpose(flip(Dphip(Phi4,2:nt),2));transpose(Dphip(Phi3,1:nt))];
polar2(theta,D_el,theta,D_el2,"","","Directivity");
legend("D(theta) - E-plane","D(phi) - H-plane");
setplot("legend position",1);

#####################################################
# 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")*TFeed(fMinIndex)/TBox(fMinIndex); # 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 = "=======================================";