clear;
basename = "twlm_selfheating";

# Sweep name and number of current points in the sweep.
sweep_name = "sweep_370K";
numpoints = 9;

# Go into the sweep folder, open sweep files for each current and extract recombination
# rates and thermionic leakage current. Transform the recombination rates into currents.
cd(basename + "_" + sweep_name);
T = zeros(numpoints);
temperature = zeros(numpoints);
current = zeros(numpoints);
counter = 0;
plotIndex = zeros(numpoints);
for(i=1:numpoints){
    load(sweep_name + "_" + num2str(i));
    power = getresult("OPWM_1","sum/power");
    if(power > 2e-3){
        temp = getresult("TWLM_1","active layer temperature profile");
        temp = temp.getattribute("temperature (K)");
        temperature(i) = mean(temp(end,:)); #average over position for the last temperature iteration
        current(i) = getnamed("DC_1","amplitude");
        counter = counter + 1;
        plotIndex(counter) = i;
    }
}
plotIndex = plotIndex(1:counter);
temperature = temperature(plotIndex);
current = current(plotIndex);

# Plot different components of the total injected current.
plot(current,temperature,"","","","linewidth=2");

setplot("x label","current [mA]");
setplot("y label","average temperature [K]");

# Exit the sweep folder and reopen the main project file.
cd("../");
load('twlm_selfheating.icp');