##############################################################
# dipole_direct.lsf
#
# Uses FDTD to calculates the directivity pattern, max
# directivity, and radiated power of an ideal dipole 
# radiating into unbounded free space and an ideal dipole
# radiating into a bounded free space where z<0 is an infinite 
# PEC ground plane. The script runs both simulation setups 
# and generates the radiation plots and outputs the input 
# power, radiated power, radiation efficiency, and max 
# directivity in the script prompt.
#
# Copyright 2016 Lumerical Solutions Inc
##########################################################

closeall;clear;
freq=c/0.45e-6;    #Define frequency of interest

##############################################################################################
#Unbounded antenna simulation setup
switchtolayout;
select("Ground");set("enabled",0);
select("PEC");set("enabled",1);
select("Directivity");set("inf gp", 0);set("z span",0.2e-6);
run;

##############################################################################################
#Retrieve simulation results from analysis group
select("Directivity");set("frequency",freq);
runanalysis;
DirMon=getresult("Directivity","Directivity");
Prad=getresult("Directivity","Prad");
theta=DirMon.theta;phi=DirMon.phi; nt=length(theta);np=length(phi);p1=find(phi,0);t1=find(theta,pi/2);
Dtheta=10*log10(DirMon.Dtheta); Dphi=10*log10(DirMon.Dphi);Dmax=max(Dtheta);
#Define theoretical directivity
Dtheory=10*log10(abs(sin(theta))^2)+1.76;
Dtheory(1)=-316; ##Replace NaN with s

##############################################################################################
#Plot directivity patterns
holdoff;plot([-flip(theta(2:nt),1);theta]*180/pi,[flip(Dtheory(2:nt),1);Dtheory(1:nt)]);    
holdon;plot([-flip(theta(2:nt),1);theta]*180/pi,[flip(Dtheta(p1,2:nt),2),Dtheta(p1,1:nt)]);
plot(phi*180/pi-180,Dtheory(t1)+matrix(np,1));plot(phi*180/pi-180,Dtheta(1:np,t1));
legend("Theory - Elevation","FDTD - Elevation","Theory - Azimuth","FDTD - Azimuth");setplot("title", "Unbounded Ideal Dipole - Dtheta");
setplot("x min",-180);setplot("x max",180);setplot("x label", "(deg)");setplot("y label", "(dB)");
setplot("y max",5);setplot("y min",-25);

##############################################################################################
#Make polar plot
Dtheta=Dtheta+(20-max(Dtheta));  #Total directivitiy normalized such that 20dB is max direct.
#Filter out results below 0dB.
index = find(Dtheta>0);filter=matrix(np,nt);filter(index)=1;Dtheta=Dtheta*filter;  #Normalized and filtered total directivitiy
thetaplot=[-flip(theta(2:nt),1); theta];D_el=[transpose(flip(Dtheta(1,2:nt),2));transpose(Dtheta(1,1:nt))];
polar2(thetaplot,D_el,phi,Dtheta(1:np,t1));
legend("Elevation Plane","Azimuth Plane");setplot("title","Unbounded Ideal Dipole - Dtheta");

##############################################################################################
#Generate power quantities and directivity in script prompt
Pin=dipolepower(freq);
?msg = "==========Unbounded Dipole==============";
?msg = "Incident Power: " + num2str(round(Pin/1e-19)*1e-3) + " fW";
?msg = "Radiated Power: " + num2str(round(Prad/1e-19)*1e-3) + " fW";
?msg = "Radiation Efficiency: " + num2str(round(Prad*1e3/Pin)/1e1) + " %" ;
?msg = "Directivity: " + num2str(round(max(Dmax)*1e3)/1e3) + " dB" +endl ;

##############################################################################################
##############################################################################################
#Bounded Antenna simulation setup
#clear;
switchtolayout;
select("Ground");set("enabled",1);
select("PEC");set("enabled",1);
select("Directivity");set("inf gp",1);set("z span",0.1e-6);
run;

##############################################################################################
#Retrieve simulation results from analysis group
select("Directivity");runanalysis;
DirMon=getresult("Directivity","Directivity");
Prad=getresult("Directivity","Prad");
theta=DirMon.theta;phi=DirMon.phi; nt=length(theta);np=length(phi);p1=find(phi,0);t1=find(theta,pi/2);
Dtheta=10*log10(DirMon.Dtheta); Dphi=10*log10(DirMon.Dphi);Dmax=max(Dtheta);
Dtheory=10*log10(abs(sin(theta))^2)+4.76;
Dtheory(1)=-316; ##Replace NaN with small value

##############################################################################################
#Plot directivity patterns
holdoff;plot([-flip(theta(2:nt),1);theta]*180/pi,[flip(Dtheory(2:nt),1);Dtheory(1:nt)]);    
holdon;plot([-flip(theta(2:nt),1);theta]*180/pi,[flip(Dtheta(p1,2:nt),2),Dtheta(p1,1:nt)]);
plot(phi*180/pi-180,Dtheory(t1)+matrix(np,1));plot(phi*180/pi-180,Dtheta(1:np,t1));
legend("Theory - Elevation","FDTD - Elevation","Theory - Azimuth","FDTD - Azimuth");setplot("title", "Bounded Ideal Dipole - Dtheta");
setplot("x min",-180);setplot("x max",180);setplot("x label", "(deg)");setplot("y label", "(dB)");
setplot("y max",5);setplot("y min",-25);

##############################################################################################
#Make polar plot
Dtheta=Dtheta+(20-max(Dtheta));  #Total directivitiy normalized such that 20dB is max direct.
#Filter out results below 0dB.
index = find(Dtheta>0);filter=matrix(np,nt);filter(index)=1;Dtheta=Dtheta*filter;  #Normalized and filtered total directivitiy
thetaplot=[-flip(theta(2:nt),1); theta];D_el=[transpose(flip(Dtheta(1,2:nt),2));transpose(Dtheta(1,1:nt))];
polar2(thetaplot,D_el,phi,Dtheta(1:np,t1));
legend("Elevation Plane","Azimuth Plane");setplot("title","Bounded Ideal Dipole - Dtheta");

##############################################################################################
#Generate power quantities and directivity in script prompt
Pin=dipolepower(freq);
?msg = "==========Bounded Dipole==============";
?msg = "Incident Power: " + num2str(round(Pin/1e-19)*1e-3) + " fW";
?msg = "Radiated Power: " + num2str(round(Prad/1e-19)*1e-3) + " fW";
?msg = "Radiation Efficiency: " + num2str(round(Prad*1e3/Pin)/1e1) + " %" ;
?msg = "Directivity: " + num2str(round(max(Dmax)*1e3)/1e3) + " dB"+endl;
clear;