lmin = 2.0;          # source min wavelength
lmax = 5.0;          # source max wavelength
fmin = 1/lmax;       # source min frequency
fmax = 1/lmin ;      # source max frequency
fcen = 0.5*(fmin+fmax);

thetas = [0:5:30];
kx = fcen*sind(thetas);

Refl = zeros(50,length(thetas));
Abs = zeros(50,length(thetas));
theta_out = zeros(50,length(thetas));

for k = 1:length(thetas)
  f0 = dlmread(sprintf("flux0_a4.3_theta%d.dat",thetas(k)),",");
  f = dlmread(sprintf("flux_a4.3_r1.72_theta%d.dat",thetas(k)),",");
  Refl(:,k) = -f(:,2)./f0(:,2);
  theta_out(:,k) = asind(kx(k)./f0(:,1));
  Abs(:,k) = (1-Refl(:,k)).*cosd(theta_out(:,k));
endfor

Abs(Abs < 0) = 0;
wvl = 1./f0(:,1);
wvls = repmat(wvl,1,length(thetas));

pcolor(theta_out, wvls, Abs);
shading flat; c = colormap("hot"); colormap(flipud(c)); colorbar;
eval(sprintf("axis([%0.2g %0.2g %0.2g %0.2g])",min(min(theta_out)),max(max(theta_out)),min(wvl),max(wvl)));
set(gca, 'xtick', [0:10:50]);
set(gca, 'ytick', [2.0:0.5:5.0]);
xlabel("emission angle (degrees)")
ylabel("wavelength (um)")
title("emissivity: a=4.3 um, r=1.72 um")