% Part of Supplementary Information for "New Constraints on the Thermal Conductivity of the Upper Mantle from Numerical Models of Radiation Transport"
% C. J. Grose and J. C. Afonso
% Published in Geochemistry, Geophysics, Geosystems
% doi: xxx
% 
% This code calculates the radiative thermal conductivity (krad, W/m/K) using our effective medium theory (Eq. 3,4,6,9,15,16,19)
% The main input parameters are the phase fractions, inter-phase contact statistics, grain size (cm), temperature range (K), and pressure (GPa).

%% VARIABLE INPUTS
phaseFracs           =  [0 0.56   0.15  0.11  0.18];                           % phase fractions [melt, olivine, opx, cpx, grt]
fkl                  =  [0 0      0     0     0;                               % inter-phase contact statistics (contact area between phase k and l)
                         0 0.0660 0     0     0;
                         0 0.227  0.002 0     0;
                         0 0.188  0.076 0.001 0;
                         0 0.251  0.094 0.093 0.002];

rcm                  =   0.0001*exp(0.0138293399*(0:999)');                    % grain size array, cm
T                    =   (500:10:2000)';                                       % temperature array in K
P                    =   2;                                                    % pressure, GPa

spectralRes          =   500;                                                  % resolution of spectral integration
dnn                  =   0.0002;                                               % resolution of integration over n

%% CONSTANT INPUTS
% Absorption spectra model coefficients
olmodel              =   [500,700,800,1000,1200,1400,1600,2000,2500,3000;100,100,100,100,100,100,100,100,100,100;19.5000000000000,24,25,27.5000000000000,28,30,31,33,36,36;30,28,29,27.2000000000000,28,27.2000000000000,26,22,22,29;15,20,20,20,20,20,20,15,14,10;0.400000000000000,0.500000000000000,0.700000000000000,1.10000000000000,1.60000000000000,2.10000000000000,2.70000000000000,3.80000000000000,4.90000000000000,5.70000000000000;0,0,0,0,0,0,0,0,0,0;7800,7600,7450,7350,7100,6920,6800,6600,6200,5200;9475,9350,9350,9300,9150,9050,8950,8750,8200,7900;11050,10800,10800,10800,10800,10800,10800,10800,10800,10800;26800,26800,26800,26800,26800,26800,26800,26800,26800,26800;1200,1200,1200,1200,1200,1200,1200,1200,1200,1200;760,890,950,1150,1250,1400,1600,2100,2600,3000;565,600,720,770,900,950,1000,1300,1400,1700;1250,1330,1330,1330,1330,1330,1330,1550,1700,1700;20,20,20,20,20,20,20,20,20,20];
opxcpxgrt            =   [27,150,100,15,70,400,0,0,0;34000,30000,11000,5200,0,0,0,0,0;13000,2000,650,800,1050,500,0,0,0;0,0,0,0,0,0,0,0,0;60,250,19,18,6,5,400,0,0;33000,30000,13000,9500,4000,1500,300,0,0;10000,2300,3500,1100,500,1300,600,0,0;0.00025,0.001833333,-0.000333333,-0.000333333,0,0,0,0,0;16,150,6,28,17,14,25,1500,225;33000,30000,18700,7800,5800,4100,1700,500,1825;8000,1000,2000,900,700,700,650,400,75;0.000666667,0.001333333,0.000666667,0.000666667,0.001333333,0.001,0.001,0,0;0.11,2,0,0,0,0,0,0,0;19000,36000,0,0,0,0,0,0,0;10000,6200,0,0,0,0,0,0,0;0.6,4.3,0,0,0,0,0,0,0;19000,37250,0,0,0,0,0,0,0;8500,8100,0,0,0,0,0,0,0;0.95,5.1,0,0,0,0,0,0,0;19000,38500,0,0,0,0,0,0,0;7000,10000,0,0,0,0,0,0,0];

c                    =   299792458;                                            % speed of light in vacuum
h                    =   6.6260695729*10^-34;                                  % Planck constant
kB                   =   1.380648813*10^-23;                                   % Boltzmann constant

% Array properties/formatting/corrections
numr                 =   numel(rcm);
numT                 =   numel(T);
numphase             =   numel(phaseFracs);
phaseFracs           =   bsxfun(@rdivide,phaseFracs,sum(phaseFracs,2));        % phase fractions must sum to 1

%% REFRACTIVE INDEX AND SCATTERING INTERACTION COEFFICIENTS
%                            refractive indices assumed for each phase:
%                            melt  ol     opx    cpx    grt
n0                   =      [1.500 1.640  1.680  1.685  1.731];                % min n
n1                   = n0 + [0.001 0.045  0.045  0.045  0.030];                % max n
n2                   = n1(2:end)'.^2;

% calculate n_kl matrix of mean mismatches among all phases
nkl = zeros(numphase,numphase);
for k=1:numphase
    for l=1:numphase
        refk         =   n0(k):dnn:n1(k);
        refl         =   n0(l):dnn:n1(l);
        
        refractmatk  =   repmat(refk,numel(refl),1);
        refractmatl  =   repmat(refl',1,numel(refk));
        nkl(k,l)     =   mean(mean(abs(refractmatk-refractmatl)));             % Eq. 22
    end
end
dnev                 =   sum(bsxfun(@times,nkl(:)',fkl(:)'),2);                % Eq. 19

% effective n
avgnp                =   (n0+n1)./2;                                           % average n for each phase
avgn                 =   sum(bsxfun(@times,phaseFracs,avgnp),2);               % Mean n (Eq. 17)

% scattering interaction term
eev                  =   dnev./avgn;                                           % Eq. 16

% effective n^2
avgn2p               =   ((n1.^3-n0.^3)./3)./(n1-n0);                          % average n^2 for each phase
avgn2                =   sum(bsxfun(@times,phaseFracs,avgn2p),2);              % arithmetic average for all models (eq. 12)

%% ABSORPTION SPECTRA
% some olivine model bits (originates from fits to Hofmeister, 2005)
Tind                 =   [500,700,800,1000,1200,1400,1600,2000,2500,3000];     % temperature index for olivine (for interpolating olivine Abs)
Tindr                =   repmat(Tind,numT,1);
Ti                   =   repmat(T,1,10);
Tmin                 =   sum((Tindr-Ti)<0,2);
Tmin(Tmin<1)         =   1;
Tdif                 =   (T-Tind(Tmin)')./(Tind(Tmin+1)'-Tind(Tmin)');
Tdif(Tdif<0)         =   0;

% Calculate harmonic and arithmetic means for each temperature and spectral band
absMean              =   zeros(numT,spectralRes);                              % Result{temperature}(grainsize,waveband)dkdp
absHarm              =   zeros(numT,spectralRes);                              % Result{temperature}(grainsize,waveband)
for t=1:numT
    [absorbance,absT,freq_min,freq_max,freq,dfreq(t)] = radspectra(spectralRes,olmodel,opxcpxgrt,T(t),Tmin(t),Tdif(t),c,3000);
    for a=1:spectralRes
        absMean(t,a) =    sum(phaseFracs.*absorbance(a,:));                    % arithmetic model
        absHarm(t,a) = 1./sum(phaseFracs./absorbance(a,:));                    % harmonic model
    end
end

%% RADIATION SOURCE FUNCTION
dT                   =   0.0001;                                               % finite difference in T
I_vb                 =   zeros(numT,spectralRes);                              % black body radiation at T
I_vb_dT              =   zeros(numT,spectralRes);                              % black body radiation at T+dT
for t=1:numT
    I_vb(t,:)        =   ((2*h*avgn2).*freq.^3./c.^2./(exp(h*freq./(T(t)*kB))-1)).*dfreq(t);
    I_vb_dT(t,:)     =   ((2*h*avgn2).*freq.^3./c.^2./(exp(h*freq./((T(t)+dT)*kB))-1)).*dfreq(t);
end
I_vb_k               =   (I_vb_dT-I_vb)./dT;                                   % d(n^2I_vb)/dT, W/m^3/K

%% CALCULATE KRAD USING ROSSELAND DIFFUSION EQUATION
krad_mat             =   zeros(numr,numT,spectralRes);   % matrix for krad for each (grain size, temperature,spectral band)
for i=1:numr
    r                =   rcm(i)/100;                                           % grain radius, cm-->m
    sigma            =   eev./r;                                               % Eq. 15
    
    b                =   2./absMean;                                           % Eq. 11
    u                =   b./rcm(i);
    effAbs           =   u.*(1-exp(-1./u)).*(absMean-absHarm)+absHarm;         % Eq. 9
    
    extinct          =   effAbs.*100+repmat(sigma,1,spectralRes);              % Extinction coefficient
    meanfree         =   1./extinct;                                           % mean free path
    
    krad_mat(i,:,:)  =   (4*pi/3).*(meanfree.*I_vb_k);                         % actual krad
end
krad                 =   sum(krad_mat,3)';

%% ADD LATTICE CONDUCTIVITY TO GET TOTAL
% pressure dependence of lattice conductivity
dKdP                 =   135+5.1*P;                                            % bulk modulus
KK                   =   5.1/dKdP;
dkdp                 =   exp(KK.*P);

% klat models for each phase
klatmat              =   ones(numT,5);
klatmat(:,1)         =   1.5;                                                  % melt
klatmat(:,2)         =  -0.11 + 1930./T + 0.00125.*T;                          % olivine
klatmat(:,3)         =   0.40 + 1520./T + 0.00074.*T;                          % opx
klatmat(:,4)         =   2.35 +  747./T - 0.00030.*T;                          % cpx
klatmat(:,5)         =   2.44 +  420./T + 0.00004.*T;                          % grt

klat                 =   prod(klatmat.^phaseFracs,2).*dkdp;                    % Eq. 23
klat                 =   repmat(klat,[1,numr]);
rmmat                =   repmat(rcm',[numT,1])./100;
klat                 =   1./(1./klat+(1./rmmat).*10^-8);                       % Eq. 26

%% TOTAL THERMAL CONDUCTIVITY
k                    =   klat+krad;                                            

%% VISUALIZATION

figure(1);
subplot(1,3,1); imagesc([min(log10(rcm)) max(log10(rcm))],[min(T) max(T)],klat);   colorbar;   title('Lattice thermal conductivity, W/m/K');     xlabel('Grain size, log_{10}(cm)');   ylabel('Temperature, K'); caxis([0 7]);
subplot(1,3,2); imagesc([min(log10(rcm)) max(log10(rcm))],[min(T) max(T)],krad);   colorbar;   title('Radiative thermal conductivity, W/m/K');   xlabel('Grain size, log_{10}(cm)');   ylabel('Temperature, K'); caxis([0 7]);
subplot(1,3,3); imagesc([min(log10(rcm)) max(log10(rcm))],[min(T) max(T)],k);      colorbar;   title('Total thermal conductivity, W/m/K');       xlabel('Grain size, log_{10}(cm)');   ylabel('Temperature, K'); caxis([3 10]);
colormap(jet(128));

    