#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);
}

function confinementFactor(kvec,domain){

### DOCS
#
# kvec = [0;1;0] = propagation direction
# domain = 3 = domain ID with respect to which the confinement 
#               factor is calculated
#
### END DOCS

myfield =getresult("FEEM","fields");
numFieldVertices = size(myfield.x,1);
rField = [myfield.x,myfield.y, myfield.z];

# make 2D vector in cross sectional plane
rFieldTrans = matrix(numFieldVertices,2);
j=0;
for(i=1:3){
    if(kvec(i) < 0.5){
        j=j+1;
        rFieldTrans(:,j) = rField(:,i);
    }
}
    
SField = myfield.S(:,1,1,:);
SField = pinch(SField);
powerField = real(mult(SField,kvec));#parallel

mygrid =getresult("FEEM","grid");
numGridVertices = size(mygrid.x,1);

id = mygrid.ID;    
idGain = find(id==domain);

powerGrid = matrix(numGridVertices,1);
rGrid = [mygrid.x,mygrid.y, mygrid.z];

# make 2D vector in cross sectional plane
rGridTrans = matrix(numGridVertices,2);
j=0;
for(i=1:3){
    if(kvec(i) < 0.5){
        j=j+1;
        rGridTrans(:,j) = rGrid(:,i);
    }
}

domainElements = mygrid.elements(idGain,:);

for(i=1:numGridVertices){
    xequal = find(almostequal(rFieldTrans(:,1),rGridTrans(i,1)));
    xyequal = find(almostequal(rFieldTrans(xequal,2),rGridTrans(i,2)));
    powerGrid(i) = mean(powerField(xequal(xyequal)));
}  
    

totalPowerGrid = quadtri(mygrid.elements,rGridTrans,powerGrid);
domainPowerGrid = quadtri(domainElements,rGridTrans,powerGrid);

result = domainPowerGrid/totalPowerGrid;

return result;

}

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

wgSweep = struct;
wgSweep.input_file = common.name;
wgSweep.feem_template_filename ="feem.ldev";
wgSweep.feem_run_filename = common.name+".ldev";
wgSweep.feem_result_filename=common.name+"_opt";

if(fileexists(wgSweep.feem_run_filename)) {rm(wgSweep.feem_run_filename);}
cp(wgSweep.feem_template_filename,wgSweep.feem_run_filename);

load(wgSweep.feem_run_filename);

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

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

wgSweep.kvec = opt.propagation_direction; # propagation direction
wgSweep.domain =3; # gain material domain ID

wgSweep.f_basename = filebasename(currentfilename);

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);

for(i=1:length(wgSweep.freq)){
    
    switchtolayout;
    setnamed("FEEM","wavelength",c/wgSweep.freq(i));
    
    setnamed("materials::InP - Seifert::InP - Seifert",
    "refractive index",opticalData.InP(i));
    
    setnamed("materials::InGaAsP - 1.15um::InGaAsP - 1.15um",
    "refractive index",opticalData.InGaAsP_1_15um_SCL(i));
    
    setnamed("materials::InGaAsP - 1.40um::InGaAsP - 1.40um",
    "refractive index",opticalData.InGaAsP_ActiveMQW(i));
    
    #output which simulation is running
    ?"saving simulation " + num2str(i) + " of " + num2str(length(wgSweep.freq));
    
    wgSweep.f_name=wgSweep.f_basename+"_"+num2str(i);
    save(wgSweep.f_name);
    addjob(wgSweep.f_name);
}

runjobs;


wgSweep.neff = wgSweep.loss = 
wgSweep.confinement_factor = zeros(length(wgSweep.freq),1);
for(i=1:length(wgSweep.freq)){
    
    wgSweep.f_name=wgSweep.f_basename+"_"+num2str(i);
    load(wgSweep.f_name);

    
    myr=getresult("FEEM", "modeproperties");
    y=myr.getattribute("neff");
    wgSweep.neff(i) = y(1);
    y=myr.getattribute("loss");
    wgSweep.loss(i) = y(1);
    
    wgSweep.confinement_factor(i) = confinementFactor(wgSweep.kvec,
                                                   wgSweep.domain);
    
}

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

clear(i,myr,y);

load(wgSweep.feem_run_filename);
jsonsave(wgSweep.feem_result_filename,wgSweep);