#INPUT

filename = "Piprek2000OQE";

#END INPUT

function differentiate(y,x){

    result = zeros(length(y),1);

    i=1;
    result(i) = (y(i+1) - y(i))/(x(i+1) - x(i));    
    
    for(i=2:length(y)-1){
        result(i) = (y(i) - y(i-1))/(x(i)-x(i-1))/2
                  + (y(i+1) - y(i))/(x(i+1) - x(i))/2;                
    }

    i=length(y);
    result(i) = (y(i) - y(i-1))/(x(i)-x(i-1));
    
    return result;
}

function calculateGroupIndexDeriv(freq,neff){
 
    ng=zeros(length(neff),1);
    dneff = differentiate(neff,freq);
    for(i=1:length(neff)){
        ng(i) = neff(i)+freq(i)*dneff(i);   
    }
    return ng;
}

function cot(x){
    y = 1/tan(x);
    return y;
}

function n_InGaAsP__InP_HHI(comp,freq){
    
    E =h*freq/e;
    if(isstruct(comp)){
        x=comp.x;
        y=comp.y;
        Eg  = 1.35 +  0.668*x - 1.068*y +0.758*x^2 +  0.078*y^2 - 0.069*x*y - 0.322*x^2*y +0.03*x*y^2;
    }
    else{
        x=y=0;
        Eg = h*c/comp/e;      
    }
    a_xy = 5.8688 -0.417*x + 0.1896*y + 0.0125*x*y; # Angstroms
    a_InP = 5.868800000000001e-010/1e-10; # Angstroms, value from Device DB
    f = (a_xy  - a_InP)/a_InP;
    
    R = -0.00115 +0.00191*Eg;
    Gamma = -0.000691+0.00433*Eg;
    A = -0.0453+2.1103*Eg;
    a = 72.32+12.78*Eg;
    b = 4.84+4.66*Eg;
    cc = -0.015+0.02*Eg;
    d = -0.178+1.042*Eg;
        
    epsilon = 1 + a/(b-(E+1i*Gamma)^2)
           + A*sqrt(R)/(E+1i*Gamma)^2*(log(Eg^2/(Eg^2 - (E+1i*Gamma)^2)) 
           + pi*(2*cot(pi*sqrt(R/Eg))-cot(pi*sqrt(R/(Eg-(E+1i*Gamma))))
                                     -cot(pi*sqrt(R/(Eg+(E+1i*Gamma)))) ) );
    
    #n= sqrt(epsilon);
    n= sqrt((abs(epsilon) + real(epsilon))/2.0);
          
    res =struct;
    
    res.Eg = Eg;
    res.R=R;
    res.Gamma = Gamma;
    res.A = A;
    res.a = a;
    res.b = b;
    res.cc = cc;
    res.d = d;
    res.f = f;
    res.n = n;
    res.epsilon = epsilon;
    
    return n;
    
## Revised refractive index and absorption of
## In(1-x)Ga(x)As(y)P(1-y) lattice-matched to InP in
## transparent and absorption IR-region,
## Sten Seifert and Patrick Runge,   
## 1 Feb 2016 | Vol. 6, No. 2 | DOI:10.1364/OME.6.000629 | 
## OPTICAL MATERIALS EXPRESS p629
    
    
}

function n_InGaAsP__InP_HHI_Eg(Eg,freq){
    
    E =h*freq/e; 
    
    #Eg  = 1.35 +  0.668*x - 1.068*y +0.758*x^2 +  0.078*y^2 - 0.069*x*y - 0.322*x^2*y +0.03*x*y^2;
    x=y=0.0;    
    a_xy = 5.8688 -0.417*x + 0.1896*y + 0.0125*x*y; # Angstroms
    a_InP = 5.868800000000001e-010/1e-10; # Angstroms, value from Device DB
    f = (a_xy  - a_InP)/a_InP;

    
    R = -0.00115 +0.00191*Eg;
    Gamma = -0.000691+0.00433*Eg;
    A = -0.0453+2.1103*Eg;
    a = 72.32+12.78*Eg;
    b = 4.84+4.66*Eg;
    cc = -0.015+0.02*Eg;
    d = -0.178+1.042*Eg;
        
    epsilon = 1 + a/(b-(E+1i*Gamma)^2)
           + A*sqrt(R)/(E+1i*Gamma)^2*(log(Eg^2/(Eg^2 - (E+1i*Gamma)^2)) 
           + pi*(2*cot(pi*sqrt(R/Eg))-cot(pi*sqrt(R/(Eg-(E+1i*Gamma))))
                                     -cot(pi*sqrt(R/(Eg+(E+1i*Gamma)))) ) );
    
    #n= sqrt(epsilon);
    n= sqrt((abs(epsilon) + real(epsilon))/2.0);
          
    res =struct;
    
    res.Eg = Eg;
    res.R=R;
    res.Gamma = Gamma;
    res.A = A;
    res.a = a;
    res.b = b;
    res.cc = cc;
    res.d = d;
    res.f = f;
    res.n = n;
    res.epsilon = epsilon;
    
    return n;
    
## Revised refractive index and absorption of
## In(1-x)Ga(x)As(y)P(1-y) lattice-matched to InP in
## transparent and absorption IR-region,
## Sten Seifert and Patrick Runge,   
## 1 Feb 2016 | Vol. 6, No. 2 | DOI:10.1364/OME.6.000629 | 
## OPTICAL MATERIALS EXPRESS p629
    
    
}

function createOpticalDataPiprek1550(freq){

result= struct;

result.freq=freq;

#InP;
comp=struct;
comp.x=0.0;
comp.y=0.0;
index=zeros(length(freq),1);
for (i=1:length(freq)){
    index(i) = n_InGaAsP__InP_HHI(comp,freq(i));
}
result.InP = index;
clear(comp);

#InGaAsP_1_15um_SCL
comp=1.15e-6;
index=zeros(length(freq),1);
for (i=1:length(freq)){
    index(i) = n_InGaAsP__InP_HHI(comp,freq(i));
}
result.InGaAsP_1_15um_SCL = index;
clear(comp);

#InGaAsP_ActiveMQW;
compBarr=struct;
compBarr.x=0.29;
compBarr.y=0.55;
depthEdgeBarr = 17e-9;
depthBarr=5.5e-9;
depthWell=6.4e-9;
compWell=struct;
compWell.x=.24;
compWell.y=.79;
index=zeros(length(freq),1);
for (i=1:length(freq)){
    index(i) = (  5*depthBarr*n_InGaAsP__InP_HHI(compBarr,freq(i))
                + 6*depthWell*n_InGaAsP__InP_HHI(compWell,freq(i))
                + 2*depthEdgeBarr*n_InGaAsP__InP_HHI(compBarr,freq(i)) )
                / (5*depthBarr + 6*depthWell +2*depthEdgeBarr);
}
result.InGaAsP_ActiveMQW = index;

return result;

}

function createOpticalData(freq,mat){

result = cell(length(mat));

j=0;
for(i=1:(length(mat))){

    sumindex = zeros(length(freq),1);
	sumthick = 0.0;
	
	if( iscell(mat{i}) ) {
		sublay = mat{i};
		numsublay = length(mat{i});
	}
	else{
		sublay = cell(1);
		sublay{1} = mat{i};
		numsublay=1;
	}
	
    for(k=1:numsublay){		
		if(sublay{k}.mat=="InGaAsP"){
		    comp = sublay{k}.comp;
		    thick = sublay{k}.thick;
		    for(m=1:length(freq)){
				sumindex(m) = sumindex(m)+thick*n_InGaAsP__InP_HHI(comp,freq(i));
		    }
	    }
		else if( sublay{k}.mat=="Air"){
		    comp = sublay{k}.comp;
		    thick = sublay{k}.thick;
		    for(m=1:length(freq)){
				sumindex(m) = sumindex(m)+thick*1.0;
			}
		}
		sumthick = sumthick+thick;
	}
	result{i}.thick = sumthick;
	result{i}.index = sumindex/sumthick;
	result{i}.name =  mat{i}.name;
	result{i}.numsublay =  numsublay;

}
}

function getIndex(opticalData,name,freq){

	nData = length(opticalData);
	result=-1*ones(nData,1);
	
	if(nData ==0){
        ?"no optical data";
        return result;		
	}
	  
    found = false;
	j=0;	  
	for(0 ; (j < nData) & !found ; 0) {
			  
	    j=j+1;
		found = (opticalData{j}.name == name);
	}
	return interp(opticalData{j}.index,opticalData{j}.freq,freq);
}


if(!exist('fname')){
    fname = filename;
}
jsonload(fname);

wgSweep = struct;
wgSweep.input_file = common.name;
wgSweep.eigenmode_template_filename ="mode.lms";
wgSweep.eigenmode_run_filename = common.name+".lms";
wgSweep.eigenmode_result_filename=common.name+"_opt";

if(fileexists(wgSweep.eigenmode_run_filename)) {rm(wgSweep.eigenmode_run_filename);}
cp(wgSweep.eigenmode_template_filename,wgSweep.eigenmode_run_filename);

load(wgSweep.eigenmode_run_filename);

switchtolayout;
setnamed("MQW Waveguide","n_wells",common.n_wells);
setnamed("MQW Waveguide","t_barrier",common.barrier_thickness);
setnamed("MQW Waveguide","t_barrier_bottom",common.edge1_barrier_thickness);
setnamed("MQW Waveguide","t_barrier_top",common.edge2_barrier_thickness);
setnamed("MQW Waveguide","t_well",common.well_thickness);
setnamed("MQW Waveguide","w_ridge",common.active_region_width);
save;

wgSweep.freqCentre = opt.centre_frequency;
wgSweep.freqBW = 2*opt.frequency_step;
wgSweep.numpoints = 3;

wgSweep.freqStart = wgSweep.freqCentre-wgSweep.freqBW;
wgSweep.freqEnd = wgSweep.freqCentre+wgSweep.freqBW;
wgSweep.freq=linspace(wgSweep.freqStart,wgSweep.freqEnd,wgSweep.numpoints);

opticalData=createOpticalDataPiprek1550(wgSweep.freq);

wgSweep.neff = wgSweep.loss = wgSweep.confinement_factor = zeros(length(wgSweep.freq),1);

for(i=1:length(wgSweep.freq)){
    
    switchtolayout;
    setanalysis('wavelength',c/wgSweep.freq(i));
    
    setmaterial("InP - Seifert","Refractive Index",opticalData.InP(i));
    setmaterial("InGaAsP - 1.15um","Refractive Index",opticalData.InGaAsP_1_15um_SCL(i));
    setmaterial("InGaAsP - 1.40um","Refractive Index",opticalData.InGaAsP_ActiveMQW(i));

    findmodes;

    wgSweep.neff(i) = getresult('FDE::data::mode1','neff');
    wgSweep.loss(i) = getresult('FDE::data::mode1','loss');
    
    setanalysis('integration shape','Rectangular');
    setanalysis('x1',opt.gain_domain_x1);
    setanalysis('x2',opt.gain_domain_x2);
    setanalysis('z1',opt.gain_domain_z1);
    setanalysis('z2',opt.gain_domain_z2); 
    
    wgSweep.confinement_factor(i) = get('fraction integrated');
    
}

wgSweep.ng=calculateGroupIndexDeriv(wgSweep.freq,wgSweep.neff);

jsonsave(wgSweep.eigenmode_result_filename,wgSweep);