###########################################################################################
# rect_oewg.lsf
#
# Calculates the return loss, directivity pattern, maximum directivity, gain and radiation 
# efficiency of a a rectangular open-ended waveguide (OEWG) probe antenna mounted on both a 
# finite and an infinite ground plane. The script also shows that all the power in the near 
# field going through the directivity monitor box can be accounted for with sufficient accuracy 
# to perform a meaningful 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;

simType = cell(2);
    simType{1} = "Infinite";
    simType{2} = "Finite";

for(II = 1:length(simType)){

    ############
    # Simulation
    ############

    switchtolayout;
    f0 = getglobalsource("center frequency");
    lambda0 = getglobalsource("center wavelength");
    setnamed("::model","lambda0",lambda0);
    a = getnamed("::model","a");
    b = getnamed("::model","b");
    if(simType{II} == "Infinite"){
        ground_x = 200.0e-3;
        ground_y = 2000.0e-3;
        setnamed("::model","ground_x",ground_x);
        setnamed("::model","ground_y",ground_y);
        setnamed("::model","FDTD_x",a+lambda0/2.0);
        setnamed("::model","FDTD_y",b+lambda0/2.0);
        setnamed("::model","FDTD_z",lambda0/128.0);
        setnamed("::model::directivity","inf gp",true);
        setnamed("::model::directivity","source window",false);
    }else{
        if(simType{II} == "Finite"){
            ground_x = 40.0e-3;
            ground_y = 20.0e-3;
            setnamed("::model","ground_x",ground_x);
            setnamed("::model","ground_y",ground_y);
            setnamed("::model","FDTD_x",ground_x+lambda0/2.0);
            setnamed("::model","FDTD_y",ground_y+lambda0/2.0);
            setnamed("::model","FDTD_z",lambda0/64.0);
            setnamed("::model::directivity","inf gp",false);
            setnamed("::model::directivity","source window",true);
        }else{
            ?msg = "Invalid simulation type!";        
        }
    }
    run;
    
    ####################
    # 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); # injected power
    port1TOut = abs(port1Exp.T_out); # reflected power
    port1TTot = abs(port1Exp.T_total); # total power (injected 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 Power: "+simType{II}+" Ground Plane");
    legend("T_in","T_out","T_total","T");
    setplot("x min",9.5);
    setplot("x max",10.5);
    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: "+simType{II}+" Ground Plane");
    setplot("show legend",false);
    setplot("x min",9.5);
    setplot("x max",10.5);
    setplot("y max",2.0);
    setplot("y min",-40.0);
    
    #################
    # Near Field Plot
    #################
    
    f0Index = find(f,f0);
    E2 = getelectric("::model::aperture",2);
    x = getdata("::model::aperture","x")*1.0e3;
    y = getdata("::model::aperture","y")*1.0e3;
    z = getdata("::model::aperture","z")*1.0e3;
    E2 = pinch(E2(1:length(x),1:length(y),1:length(z),f0Index));
    #E2 = 10.0*log10(E2);
    image(x,y,E2,"x (mm)","z (mm)","|E|^2 at f = "+num2str(round(fGHz(f0Index)*1e2)/1e2)+" GHz");
    
    ##################
    # 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 waveguide feed
    TFeed = port1T;
    # power plot
    plot(fGHz,TBox,TFeed,"frequency (GHz)","power","Near-Field Power: "+simType{II}+" Ground Plane");
    legend("directivity box","waveguide feed");
    setplot("x min",9.5);
    setplot("x max",10.5);
    setplot("y max",1.1);
    setplot("y min",-0.1);
    
    #############
    # Directivity
    #############
    
    k0 = 2.0*pi*f0/c;
    lambda0 = c/f0;
    select("directivity");
    set("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);
    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);
    theta = [-flip(dirMon.theta(2:nt),1);dirMon.theta];
    # theory
    Dtheory_hplane = 20.0*log10(cos(theta)*sin((k0*a/2.0)*sin(theta))/((k0*a/2.0)*sin(theta)));
    Dtheory_eplane = 20.0*log10(sin(pi*b/lambda0*sin(theta))/(pi*b/lambda0*sin(theta)));
    Dtheory = 10.0*log10(0.81*4.0*pi*(a*b)/lambda0^2);
    # polar plot
    Dthetap = Dtheta+(40.0-max(Dtheta)); # total directivitiy normalized such that 40dB is max direct.
    Dphip = Dphi+(40.0-max(Dphi));
    Dtheory_hplane = Dtheory_hplane+40.0;
    Dtheory_eplane = Dtheory_eplane+40.0;
    # 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
    index = find(real(Dtheory_hplane) > 0.0);
    filter = matrix(1,nt*2.0-1);
    filter(index) = 1;
    Dtheory_hplane = Dtheory_hplane*filter;  # normalized and filtered total directivitity
    index = find(real(Dtheory_eplane) > 0.0);
    filter = matrix(1,nt*2.0-1.0);
    filter(index) = 1;
    Dtheory_eplane = Dtheory_eplane*filter;  # normalized and filtered total directivitity
    Dtheta_eplane = [transpose(flip(Dthetap(phi4,2:nt),2));transpose(Dthetap(phi3,1:nt))];
    Dphi_hplane = [transpose(flip(Dphip(phi2,2:nt),2));transpose(Dphip(phi1,1:nt))];
    if(simType{II} == "Infinite"){
        # H-Plane
        polar2(theta,Dtheory_hplane,theta,Dphi_hplane,"","","Directivity H-Plane: "+simType{II}+" Ground plane");
        legend("theory","simulation");
        # E-Plane
        polar2(theta,Dtheory_eplane,theta,Dtheta_eplane,"","","Directivity E-Plane: "+simType{II}+" Ground plane");
        legend("theory","simulation");
    }else{
        polar2(theta,Dtheta_eplane,theta,Dphi_hplane,"","","Directivity: "+simType{II}+" Ground plane");    
        legend("E plane","H plane");
    }
    
    #####################################################
    # Power, Efficiency and Directivity for Script Prompt
    #####################################################

    Pin = port1TIn(f0Index)*sourcepower(f(f0Index)); # input power
    Pacc = port1TTot(f0Index)*sourcepower(f(f0Index)); # accepted power
    Prad = getresult("directivity","Prad"); # radiated power
    radEffIn = Prad/Pin; # radiation efficiency from input power
    radEffAcc = Prad/Pacc; # radiation efficiency from accepted power
    ?msg = "=====Radiation Performance: "+simType{II}+" Ground Plane=====";
    ?msg = "Design Frequency: "+num2str(round(fGHz(f0Index)*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*log10(radEffIn))*100.0)/100.0)+" dB ";
    ?msg = "====================================";

}