%%                                                   CALCULATING C-FACTORS
%       For testing a sticky-air approach
%               as described in 
%          Crameri et al. (GJI,2012)
%                                                Fabio Crameri, 01.02.2013
clear

%% INPUT OF MODEL VALUES
eta_st      = 1e18;    	% sticky air viscosity [Pa·s]
h_st        = 300e3;   	% sticky air thickness [m]

h_model     = 600e3; 	% depth of the model [m]
L_model     = 2800e3; 	% width of the model [m]

eta_ch      = 1e21;     % viscosity controlling the relaxation [Pa·s]

drho        = 100;      % density difference [kg/m3]
rho_ch      = 3300;     % crustal density [kg/m3]


%% CRITERION 1: isostatic
C_isost = 3/(16*pi^3) * ( L_model/h_st )^3 * (eta_st / eta_ch) ;


%% CRITERION 2: Stokes
C_Stokes  = 1/16 * drho/rho_ch * (h_model./h_st)^3 * (eta_st / eta_ch) ;


%% OUTPUT OF C-VALUES
display('**************')
display(['Using eta_st = ',num2str(eta_st,3),' Pa·s and h_st = ',num2str(h_st/1e3,3),' km'])
display(' gives:')
if C_isost<1
    display(['          C_isost  = ',num2str(C_isost,3),'    ok!'])
else
    display(['          C_isost  = ',num2str(C_isost,3),'    bad!'])
end
if C_Stokes<1
    display(['          C_Stokes = ',num2str(C_Stokes,3),'   ok!'])
else
    display(['          C_Stokes = ',num2str(C_Stokes,3),'   bad!'])
end
display('**************')