# define the materials
si = "Si with T dependence";
sio2 = "SiO2 with T dependence";

dispersive_si = "Si (Silicon) - Palik";
dispersive_sio2 = "SiO2 (Glass) - Palik";

# define dT
dT = [ -1,0,1];

# loop over dT and calculate neff
neff = matrix(length(dT));
for(i=1:length(dT)) {
  switchtolayout;
  setnamed("substrate","material",sio2);
  setnamed("waveguide","material",si);
  setanalysis("wavelength",1550e-9);
  setmaterial(si,"x1",dT(i));
  setmaterial(sio2,"x1",dT(i));
  findmodes;
  if(getnamed("FDE","x min bc")=="Symmetric") {
    neff(i) = getdata("mode2","neff");
  } else {
    neff(i) = getdata("mode3","neff");
  }
}

# set back to dispersive materials
switchtolayout;
setnamed("substrate","material",dispersive_sio2);
setnamed("waveguide","material",dispersive_si);

# calculate the change in neff from where dT = 0
neff = neff - neff(find(dT,0));

# export the TM data to file
outdata = [ pinch(dT), real(neff), imag(neff) ];
 
format long;
filename = "delta_neff_vs_dT.txt";
rm(filename);
write(filename,num2str(outdata));

# plot the results
plot(dT,real(neff),"dT (degrees C)","delta neff");
legend("TM");