% greenhouse factor gamma for varying tau_s when T = T_s (tau/tau_s)^(R/M_a/c_p)
% IJH 22 Sept 2018
clear;

%% Parameters
M_a = 29e-3; %kg/mol
R = 8.3; %J/K/mol
g = 9.8; %m/s^2
c_p = 10^3; %J/kg/K

%% Calculate gamma(tau_s)
z = linspace(0,1,101); % pts to discretise integral
tau_ss = linspace(0,10,101); % values of tau_s
gammas = NaN*tau_ss;
for j = 1:length(tau_ss)
    tau_s = tau_ss(j);
    gammas(j) = exp(-2*tau_s) + tau_s*trapz(z,2*z.^(4*R/M_a/c_p).*exp(-2*tau_s*z)); % calculate integral expression for gamma using trapezium rule
end

%% Plot gamma(tau_s)
figure(1); clf; set(gcf,'Paperpositionmode','auto','units','centimeters','position',[2 4 12 8]);
set(gcf,'DefaultAxesFontSize',12,'DefaultTextFontSize',12);
    plot(tau_ss,gammas,'b-','linewidth',2);
    ylim([0 1]);
    hold on;
    plot(tau_ss,1-2*tau_ss*4*R/(4*R+M_a*c_p),'r--'); % small tau_s approx
    plot(tau_ss,(2*tau_ss).^(-4*R/M_a/c_p).*gamma(1+4*R/M_a/c_p),'r--'); % large tau_s approx
    
    xlabel('Optical thickness $\tau_s$','interpreter','latex');
    ylabel('Greenhouse factor $\gamma$','interpreter','latex');
   
%% Save figure    
% print(gcf,'-depsc2',mfilename,'-loose');