% model spectra for melt/ol/opx/cpx/grt
function [absorbance,absT,freq_min,freq_max,freq,dfreq] = radspectra(numphotons,olmodel,opxcpxgrt,absT,Tmini,Tdifi,c,T4maxfreq) % initial photon positions

% ABSORBANCE = [MELT OL OPX CPX GRT]
% [500,700,800,1000,1200,1400,1600,2000,2500,3000]

freq_min    = 0.5*10^13;
freq_max    = 3.034*10^11*(T4maxfreq)-1.75*10^10; % BOUNDS BASED ON BLACK BODY SPECTRAL DISTRIBUTION

dfreq       = (freq_max-freq_min)/(numphotons-1);
freq        = (freq_min:dfreq:freq_max)';

wavenums_cm = 1./((c./freq).*100);

%% ABSORBANCE SPECTRA 

% PHASE 0: UPPER MANTLE 0 GPa MELT (M-GLASS, Murakami 2014)
a1 = opxcpxgrt(13,1)*1000;  b1 = opxcpxgrt(13,2)*1000;
a2 = opxcpxgrt(14,1);       b2 = opxcpxgrt(14,2);
a3 = opxcpxgrt(15,1);       b3 = opxcpxgrt(15,2);

% PHASE 0: MID-MANTLE, 40 GPa MELT (M-GLASS, Murakami 2014)
% a1 = auxp(16,1);  b1 = auxp(16,2);
% a2 = auxp(17,1);  b2 = auxp(17,2);
% a3 = auxp(18,1);  b3 = auxp(18,2);

% PHASE 0: DEEP MANTLE 80 GPa MELT (M-GLASS, Murakami 2014)
% a1 = auxp(16,1);  b1 = auxp(16,2);
% a2 = auxp(17,1);  b2 = auxp(17,2);
% a3 = auxp(18,1);  b3 = auxp(18,2);

absorbance0  = a1*exp(-((wavenums_cm-a2).^2)/2/a3^2)...
             + b1*exp(-((wavenums_cm-b2).^2)/2/b3^2);

% PHASE 1: OL
    for h=0:1
        ab_a = olmodel(2,Tmini+h);
        ab_d = olmodel(3,Tmini+h);
        ab_g = olmodel(4,Tmini+h);
        ab_j = olmodel(5,Tmini+h);
        ab_m = olmodel(6,Tmini+h);
        ab_b = olmodel(7,Tmini+h);
        ab_e = olmodel(8,Tmini+h);
        ab_h = olmodel(9,Tmini+h);
        ab_k = olmodel(10,Tmini+h);
        ab_n = olmodel(11,Tmini+h);
        ab_c = olmodel(12,Tmini+h);
        ab_f = olmodel(13,Tmini+h);
        ab_i = olmodel(14,Tmini+h);
        ab_l = olmodel(15,Tmini+h);
        ab_o = olmodel(16,Tmini+h);

        absorbance10{h+1} = ab_a*exp(-((wavenums_cm-ab_b).^2)/2/ab_c^2)+ab_d*exp(-((wavenums_cm-ab_e).^2)/2/ab_f^2)+ab_g*exp(-((wavenums_cm-ab_h).^2)/2/ab_i^2)+ab_j*exp(-((wavenums_cm-ab_k).^2)/2/ab_l^2);
        absorbance20{h+1} = ab_m*exp((wavenums_cm/ab_n).^ab_o);
        absorbance20{h+1}(wavenums_cm<10000) = 0;
        absorbance10{h+1}(absorbance20{h+1}>absorbance10{h+1}) = absorbance20{h+1}(absorbance20{h+1}>absorbance10{h+1});
    end

    absorbance1 = absorbance10{1}.*(1-Tdifi)+absorbance10{2}.*Tdifi;

% PHASE 2: OPX
a1 = opxcpxgrt(1,1);  b1 = opxcpxgrt(1,2);  c1 = opxcpxgrt(1,3);  d1 = opxcpxgrt(1,4);  e1 = opxcpxgrt(1,5);  f1 = opxcpxgrt(1,6);
a2 = opxcpxgrt(2,1);  b2 = opxcpxgrt(2,2);  c2 = opxcpxgrt(2,3);  d2 = opxcpxgrt(2,4);  e2 = opxcpxgrt(2,5);  f2 = opxcpxgrt(2,6);
a3 = opxcpxgrt(3,1);  b3 = opxcpxgrt(3,2);  c3 = opxcpxgrt(3,3);  d3 = opxcpxgrt(3,4);  e3 = opxcpxgrt(3,5);  f3 = opxcpxgrt(3,6);
a4 = opxcpxgrt(4,1);  b4 = opxcpxgrt(4,2);  c4 = opxcpxgrt(4,3);  d4 = opxcpxgrt(4,4);  e4 = opxcpxgrt(4,5);  f4 = opxcpxgrt(4,6);

absorbance2  = a1*exp(-((wavenums_cm-a2).^2)/2/a3^2).*(1+(absT-500)*a4)...
             + b1*exp(-((wavenums_cm-b2).^2)/2/b3^2).*(1+(absT-500)*b4)...
             + c1*exp(-((wavenums_cm-c2).^2)/2/c3^2).*(1+(absT-500)*c4)...
             + d1*exp(-((wavenums_cm-d2).^2)/2/d3^2).*(1+(absT-500)*d4)...
         	 + e1*exp(-((wavenums_cm-e2).^2)/2/e3^2).*(1+(absT-500)*e4)...
             + f1*exp(-((wavenums_cm-f2).^2)/2/f3^2).*(1+(absT-500)*f4);

% PHASE 3: CPX
a1 = opxcpxgrt(5,1);  b1 = opxcpxgrt(5,2);  c1 = opxcpxgrt(5,3);  d1 = opxcpxgrt(5,4);  e1 = opxcpxgrt(5,5);  f1 = opxcpxgrt(5,6);  g1 = opxcpxgrt(5,7);
a2 = opxcpxgrt(6,1);  b2 = opxcpxgrt(6,2);  c2 = opxcpxgrt(6,3);  d2 = opxcpxgrt(6,4);  e2 = opxcpxgrt(6,5);  f2 = opxcpxgrt(6,6);  g2 = opxcpxgrt(6,7);
a3 = opxcpxgrt(7,1);  b3 = opxcpxgrt(7,2);  c3 = opxcpxgrt(7,3);  d3 = opxcpxgrt(7,4);  e3 = opxcpxgrt(7,5);  f3 = opxcpxgrt(7,6);  g3 = opxcpxgrt(7,7);
a4 = opxcpxgrt(8,1);  b4 = opxcpxgrt(8,2);  c4 = opxcpxgrt(8,3);  d4 = opxcpxgrt(8,4);  e4 = opxcpxgrt(8,5);  f4 = opxcpxgrt(8,6);  g4 = opxcpxgrt(8,7);

absorbance3  = a1*exp(-((wavenums_cm-a2).^2)/2/a3^2).*(1+(absT-600)*a4)...
             + b1*exp(-((wavenums_cm-b2).^2)/2/b3^2).*(1+(absT-600)*b4)...
             + c1*exp(-((wavenums_cm-c2).^2)/2/c3^2).*(1+(absT-600)*c4)...
             + d1*exp(-((wavenums_cm-d2).^2)/2/d3^2).*(1+(absT-600)*d4)...
             + e1*exp(-((wavenums_cm-e2).^2)/2/e3^2).*(1+(absT-600)*e4)...
             + f1*exp(-((wavenums_cm-f2).^2)/2/f3^2).*(1+(absT-600)*f4)...
             + g1*exp(-((wavenums_cm-g2).^2)/2/g3^2).*(1+(absT-600)*g4);

% PHASE 4: GRT
a1 = opxcpxgrt(9,1);   b1 = opxcpxgrt(9,2);   c1 = opxcpxgrt(9,3);   d1 = opxcpxgrt(9,4);   e1 = opxcpxgrt(9,5);   f1 = opxcpxgrt(9,6);   g1 = opxcpxgrt(9,7);      h1 = opxcpxgrt(9,8);     i1 = opxcpxgrt(9,9);
a2 = opxcpxgrt(10,1);  b2 = opxcpxgrt(10,2);  c2 = opxcpxgrt(10,3);  d2 = opxcpxgrt(10,4);  e2 = opxcpxgrt(10,5);  f2 = opxcpxgrt(10,6);  g2 = opxcpxgrt(10,7);     h2 = opxcpxgrt(10,8);    i2 = opxcpxgrt(10,9);
a3 = opxcpxgrt(11,1);  b3 = opxcpxgrt(11,2);  c3 = opxcpxgrt(11,3);  d3 = opxcpxgrt(11,4);  e3 = opxcpxgrt(11,5);  f3 = opxcpxgrt(11,6);  g3 = opxcpxgrt(11,7);     h3 = opxcpxgrt(11,8);    i3 = opxcpxgrt(11,9);
a4 = opxcpxgrt(12,1);  b4 = opxcpxgrt(12,2);  c4 = opxcpxgrt(12,3);  d4 = opxcpxgrt(12,4);  e4 = opxcpxgrt(12,5);  f4 = opxcpxgrt(12,6);  g4 = opxcpxgrt(12,7);     h4 = opxcpxgrt(12,8);    i4 = opxcpxgrt(12,9);

absorbance4  = a1*exp(-((wavenums_cm-a2).^2)/2/a3^2).*(1+(absT-600)*a4)...
             + b1*exp(-((wavenums_cm-b2).^2)/2/b3^2).*(1+(absT-600)*b4)...
             + c1*exp(-((wavenums_cm-c2).^2)/2/c3^2).*(1+(absT-600)*c4)...
             + d1*exp(-((wavenums_cm-d2).^2)/2/d3^2).*(1+(absT-600)*d4)...
             + e1*exp(-((wavenums_cm-e2).^2)/2/e3^2).*(1+(absT-600)*e4)...
             + f1*exp(-((wavenums_cm-f2).^2)/2/f3^2).*(1+(absT-600)*f4)...
             + g1*exp(-((wavenums_cm-g2).^2)/2/g3^2).*(1+(absT-600)*g4)...
             + h1*exp(-((wavenums_cm-h2).^2)/2/h3^2).*(1+(absT-600)*h4)...
             + i1*exp(-((wavenums_cm-i2).^2)/2/i3^2).*(1+(absT-600)*i4);

            % melt        olivine     opx         cpx         grt
absorbance = [absorbance0 absorbance1 absorbance2 absorbance3 absorbance4]; % absorbances in 1/cm

absmins = (absorbance<0.5);
absorbance(absmins) = 0.5;

absmins = (absorbance>4000);
absorbance(absmins) = 4000;

