f = load('oled-flux.dat');

% total power emitted from sources
total_power = f(:,2)+f(:,3)+f(:,4)+f(:,5)+f(:,6)+f(:,7);

% total power in glass
glass_power = f(:,8);

% total power in waveguide mode
waveguide_power = f(:,9)+f(:,10)+f(:,11)+f(:,12);

if (length(find(f(:,2:12) < 0)) > 0)
  disp(sprintf("warning: flux is negative"));
endif

% fraction of total power in each of three regions
glass = glass_power./total_power;
waveguide = waveguide_power./total_power;
aluminum = 1-glass-waveguide;

% wavelengths (um)
lambdas = 1./f(:,1);

if (sum(aluminum < 0) > 0)
  disp(sprintf("warning: aluminum absorption is negative"));
endif

if (sum(aluminum > 1) > 0)
  disp(sprintf("warning: aluminum absorption is larger than 1"));
endif

lambdas_linear = linspace(0.4,0.8,100).';
glass_linear = interp1(lambdas,glass,lambdas_linear,"spline","extrap");
waveguide_linear = interp1(lambdas,waveguide,lambdas_linear,"spline","extrap");
aluminum_linear = interp1(lambdas,aluminum,lambdas_linear,"spline","extrap");
plot(lambdas_linear,glass_linear,'b-',lambdas_linear,aluminum_linear,'r-',lambdas_linear,waveguide_linear,'g-');

xlabel("wavelength (um)"); ylabel("fraction of total power");
legend("glass","aluminum","organic + ITO");
axis([0.4 0.8 0 1]);

disp("power in each region averaged over all wavelengths");
disp(sprintf("glass: %0.6f, %0.6f",mean(glass_linear),std(glass_linear,1)));
disp(sprintf("aluminum: %0.6f, %0.6f",mean(aluminum_linear),std(aluminum_linear,1)));
disp(sprintf("organic + ITO: %0.6f, %0.6f",mean(waveguide_linear),std(waveguide_linear,1)));			    
