##################################################################

# Graded index fiber (FEEM):

# This script calculates the error of a set of modes of a graded-index fiber as a function of edges per wavelength


# Copyright 2018 Lumerical Solutions, Inc.

##################################################################


clear;
#closeall;

# Modal properties:
numberOfTrialModes = 12; # number of trial modes used for the test
wavelength = 1e-6; # wavelength for mode calculation

# Fiber index parameters:
peakRefractiveIndex = 2.0; #n0
h_parameter = 40e-6; 

# Mesh properties:
edgesPerWavelength = [0.1,0.2,0.4,0.6,0.8];
polyOrder = [1,2,3];

# ---------------

load("graded_index_fiber.ldev");

switchtolayout;
select("FEEM");
set("wavelength", wavelength);
set("number of trial modes", numberOfTrialModes);

# ----------------
# Analytic Result:
# ----------------

?"Computing analytic effective indices.";

# TE mode indices
M = N = [0:ceil(sqrt(numberOfTrialModes))-1];
# effective index for each mode index pair
peakWaveNumber = (2.0*pi/wavelength)*peakRefractiveIndex;
exactNeffTE = matrix(length(M),length(N));
for(m = M){
    for(n = N){
        alpha = 2.0*(1+m+n);
        exactNeffTE(m+1,n+1) = peakRefractiveIndex*sqrt(1.0-alpha/(h_parameter*peakWaveNumber));
    }
}
clear(m,n,alpha);

# create a flat list with sorted effective indices
exactNeffList = [exactNeffTE(:);exactNeffTE(:)]; #duplicate effective indices since TE and TM are degenerate
numberOfExactModes = length(exactNeffList);

# sort the effective indices by decreasing real part
for( I = 1:(numberOfExactModes-1) ){
    for( J = (I+1):numberOfExactModes ){       
        if( real(exactNeffList(I)) < real(exactNeffList(J)) ){
            tmp = exactNeffList(J);
            exactNeffList(J) = exactNeffList(I);
            exactNeffList(I) = tmp;
            clear(tmp);
        }
    }
}
exactNeffList = exactNeffList(1:numberOfTrialModes); #Trim list of analytical modes to desired number of modes

# -----------------
# Convergence Test:
# -----------------

?"Performing convergence test.";

numericNeffAll = matrix(numberOfTrialModes, length(edgesPerWavelength),length(polyOrder));
maxNeffRelErr = matrix(length(edgesPerWavelength),length(polyOrder));

for(epw_ind = 1:length(edgesPerWavelength))
{    
    ?"Edges per wavelength = " + num2str(edgesPerWavelength(epw_ind));    
    
    switchtolayout;
    setnamed("FEEM","mesh lock",false);
    setnamed("FEEM","edges per wavelength",edgesPerWavelength(epw_ind));
    mesh;
    setnamed("FEEM","mesh lock",true);
    
    for(po_ind = 1:length(polyOrder)) 
    {
        ?"* Polynomial order = " + num2str(polyOrder(po_ind));
        switchtolayout;
        setnamed("FEEM","polynomial order",polyOrder(po_ind));
        run;
        numericNeffList = getattribute(getresult("FEEM","modeproperties"),"neff");
        
        if(length(numericNeffList)!=numberOfTrialModes){
            ?"Number of modes found numerically is less than the number of trial modes";
            break;
        }
          
        for(I = 1:(numberOfTrialModes-1)){
            for(J = (I+1):numberOfTrialModes){
                if( real(numericNeffList(I)) < real(numericNeffList(J)) ){
                    tmp = numericNeffList(J);
                    numericNeffList(J) = numericNeffList(I);
                    numericNeffList(I) = tmp;
                    clear(tmp);
                }
            }
        }
        
        maxNeffRelErr(epw_ind,po_ind) = max(abs(numericNeffList-exactNeffList)/exactNeffList);
        numericNeffAll(:,epw_ind,po_ind) = numericNeffList;
    }
  
}

switchtolayout;
setnamed("FEEM","mesh lock",false);

for(po_ind = 1:length(polyOrder))
{
    title = "Polynomial order = " + num2str(polyOrder(po_ind));
    legStr = cell(length(edgesPerWavelength)+1);
    legStr{1} = "Analytical";
    for(epw_ind = 1:length(edgesPerWavelength)){ legStr{epw_ind+1} = "edges per wavelength = "+num2str(edgesPerWavelength(epw_ind));}    
    
    plot(1:numberOfTrialModes, exactNeffList, real(numericNeffAll(:,:,po_ind)), "mode number", "neff", title);
    legend(legStr);
}

plot(edgesPerWavelength,maxNeffRelErr,"Edges per wavelength","max neff relative error","Graded index fiber: convergence test","plot points,log10y");
legStr = cell(length(polyOrder));
for(I=1:length(polyOrder)){ legStr{I} = "p = "+num2str(polyOrder(I)); }
legend(legStr);