S=input('Enter incoming energy constant S (default is 100x10^{15}): ');
emis=input('Enter emissivity: ');

if isempty(S), S=1e17; end

Esc=1e15;
T=220:310;
Etop=S;
c=3e7;
A=zeros(size(T));
%A(T<280)=0.05*(280-T(T<280));
%A(T<=270)=0.5;
Ein=Etop*(1-A);

if isempty(emis)
  emisf=input('Enter scaling for dependence of emissivity on T: ');
  emissc=min(ones(size(T)),exp(-(T-220).^2/(1e4*emisf)));
  emis=1;
else
  emissc=1;
end
emis=emis*ones(size(T));

EMIS=emissc.*emis;
Eout=c*EMIS.*T.^4;
Eoutn=c*T.^4;

figure(1), clf
subplot(2,1,1), box on, grid on, hold on
pl1=plot(T,Ein/Esc,'b'); pl2=plot(T,Eout/Esc,'r'); pl3=plot(T,Eoutn/Esc,'r--');
set(pl1,'linewidth',2), set(pl2,'linewidth',2)
legend('E_{in}','E_{out}','E_{out} if emissivity were 1','Location','Northwest')

xlabel('T (K)')
ylabel('Energy flux (PW; 1 PW = 10^{15} W)')

axis([min(T) max(T) 0 250])

subplot(2,1,2), box on, grid on, hold on
plot(T,EMIS,'k',T,ones(size(T)),'k--');
xlabel('T (K)')
ylabel('Emissivity')
axis([min(T) max(T) 0 1.1])