##  This script file will load temperature data from a user defined .mat file
##  containing a 'temperature' dataset with attribute 'T' in units of K
##  and calculate index perturbation due to the temperature variation
##  for silcion at a wavelength of 1550 nm using a linear model
##  The script will also save the index data in a dataset and
##  (optionally) load it into an n,k material attribute defined by user

clear;
closeall;
format long;

## User inputs
lambda = 1550e-9;			# options: 1550 nm only
data_filename = "T.mat";		# dataset must be named 'charge' with attributes 'n' and 'p'
data_type = "finite element";		# options: "finite element" or "rectilinear"
bias_parameter_name = "P";		# leave blank ("") if no bias sweep is present
nk_attribute_name = "FEEM::nk import";	# leave blank ("") to not load data
data_assignment_behaviour = "accumulate";	# options: "accumulate" or "overwrite"
background_index = 0;			# (keep at zero) can be used to define a constant background index

# This script is specifically for Silicon at 1550 nm

# Load charge data
matlabload(data_filename);
T = pinch(temperature.T);
x = temperature.x;
y = temperature.y;
z = temperature.z;
if(data_type == "finite element"){
    elements = temperature.elements;
}
if(length(bias_parameter_name)>0){
    eval("P = temperature."+bias_parameter_name+";");
}

## Index perturbation model (linear) parameter

dneff_dT = 0.00018;  	# dneff/dT at 1550 nm for Silicon
Tref = 300;		# at 300 K

# Calculate delta_n and delta_k
delta_n = dneff_dT * (T - Tref);

temperature.addattribute("dT",T - Tref);
temperature.addattribute("delta_n",delta_n);

index = background_index + delta_n;

if(data_type == "finite element"){
    perturbation = unstructureddataset("perturbation",x,y,z,elements);
}
else if(data_type == "rectilinear"){
    perturbation = rectilineardataset("perturbation",x,y,z);
}
if(length(bias_parameter_name)>0){
    perturbation.addparameter("P",P);
}
else{
    perturbation.addparameter("a",1);
}
perturbation.addattribute("index",index);

if(length(nk_attribute_name)>0){
    switchtolayout;
    select(nk_attribute_name);
    importdataset(perturbation);
    set("selected attribute","index");
    set("data assignment behaviour",data_assignment_behaviour);
}
format short;