function [XM_OUT,XOW_OUT,BAI,BAI3,Ts,T,CumFrac,IncFrac,Fc] = ...
    BG15(Tb,Pb,dFdP,JM,JS,Inc) 
%==========================================================================
%  BG15 Model for MORB melting in the spinel and plagioclase fields
%
%  Inputs:
%     Tb - Real temperature (i.e., adiabatic gradient must be added) (ºC)
%     Pb - Pressure (kbars) 
%           Note: Tb & Pb must be vectors of the same length (both values
%           should be either constant or monotonically decreasing as if 
%           ascending along an adiabat).
%     dFdP - Melt productivity (melt fraction/kbar).  Can either be
%           specified either as a single value or vector of length Tb & Pb.
%     JM, JS - Indices for mantle composition and mode, respectively.
%           1 - HZ prim (KG93)
%           2 - HZ dep 1 (KG93)
%           3 - Workman & Hart DMM (EPSL05)
%           4 - HZ dep 2
%           5,6 - WH DMM minus 1%, 2% melt, respectively
%           7 - MM3 (Baker & Stolper 1994; Baker et al., 1995)
%           8 - Kushiro (1996)
%           9 - HZ in CMAS
%     Inc - melt fraction extracted at each step (value between 0.01-1; 
%               where 1 = complete melt extraction (i.e., fractional 
%               melting); 0.01 = batch melting.  Note that a value of 
%               exactly 0 will result in spurious results in BAI.)  
%               Also, Inc only has meaning for calculations with multiple 
%               melt steps.
%
%  Outputs:
%     XM_OUT - Mantle mode [Aug, Opx, Oliv, Plag, Sp]
%     XOW_OUT - Mantle composition [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,
%           K2O,Na2O] (wt%)
%     BAI - Pooled melt composition at top of 1-D column (wt%)
%     BAI3 - Incremental melts [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,
%           Na2O] (wt%)
%     Ts - Solidus temperature (C)
%     T - Temperature corrected for latent heat of melting based on the 
%           amount of melt extracted at each step (C) 
%     CumFrac, IncFrac - cumulative and incrmental melt fraction
%     Fc - max(CumFrac)
%     
%  Example of batch melting with HZ dep 1 mantle composition:
%     [XM,XOW,BAI,BAI3,Ts,T,CumFrac,IncFrac,Fc] = ...
%            BG15([1350,1340],[12,11],dFdP,2,2,0.01);
%  Example of incremental batch melting with 90% melt extraction at each 
%    step using the HZ dep 1 mantle composition:
%     [XM,XOW,BAI,BAI3,Ts,T,CumFrac,IncFrac,Fc] = ...
%            BG15([1350,1340],[12,11],dFdP,2,2,0.9);
%
%
%  BG15 Version 1.0
%  Mark Behn & T.L. Grove [JGR, 2015]
%==========================================================================

%% Initialization Sequence

% Define Oxide Weights
%[SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,MnO,CaO,K2O,Na2O,Fe2O3,P2O5];
FK = [60.0840,79.8988,101.9612,151.9902,71.8464,40.3044,70.9374,56.0794,...
    94.2036,61.9790,159.6920,141.945];  

%Mantle Mode [Aug, Opx, Oliv, Plag, Sp]
ZMS(1,:) = [0.19,0.28,0.50,0,0.03];     %HZ prim
ZMS(2,:) = [0.182,0.28,0.508,0,0.0295]; %HZ dep 1
ZMS(3,:) = [.13,.28,.57,0,.02];         %WH DMM
ZMS(4,:) = [0.182,0.28,0.508,0,0.0295]; %HZ dep 2
ZMS(5,:) = [.13,.28,.57,.00,.02];       %WH DMM
ZMS(6,:) = [.13,.28,.57,.00,.02];       %WH DMM
ZMS(7,:) = [.13,.28,.57,0,.02];         %MM3
ZMS(8,:) = [.13,.28,.57,0,.02];         %Kushiro
ZMS(9,:) = [0.19,0.28,0.50,0,0.03];     %HZ in CMAS

ZMP(1,:) = [.13,.26,.54,.07,0];         %HZ prim
ZMP(2,:) = [.128,.261,.546,.065,0.0];   %HZ dep 1
ZMP(3,:) = [.128 .261 .546 .065 0.0];   %WH DMM
ZMP(4,:) = [.128,.261,.546,.065,0.0];   %HZ dep 2
ZMP(5,:) = [.128 .261 .546 .065 0.0];   %WH DMM
ZMP(6,:) = [.128 .261 .546 .065 0.0];   %WH DMM
ZMP(7,:) = [.128,.261,.546,.065,0.0];   %MM3
ZMP(8,:) = [.128,.261,.546,.065,0.0];   %Kushiro
ZMP(9,:) = [.13,.26,.54,.07,0];         %HZ in CMAS

%Melting reaction [Aug, Opx, Oliv, Plag, Sp]
YMP = [.31,0.16,-0.06,0.59,0]; %Plag
YMS = [0.97,0.25,-0.30,0.00,0.08]; % Spin (init);  P-dependent spinel 
                                   %   reaction is defined below
   
%Mantle Composition [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] wt%;

% 1 - HZ prim (KG93)
ZO(1,:) = [46.161,0.182,4.08,0.47,7.573,37.946,3.224,0.032,0.333];
% 2 - HZ dep 1 (KG93)
ZO(2,:) = [46.331,0.174,3.935,0.47,7.659,38.452,3.193,0.0074,0.281];
% 3 - WH DMM (EPSL05)
ZO(3,:) = [44.71,0.13,3.98,0.57,8.18,38.73,3.17,0.006,0.28];
% 4 - HZ dep 2
ZO(4,:) = [46.313,0.161,3.639,0.47,7.670,38.890,3.101,0.003,0.223];
% 5 - WH DMM (EPSL05) - 1% melt (calculated from BG15)
ZO(5,:) = [44.812,0.124,3.864,0.471,8.229,39.124,3.136,1.6849e-06,0.24];
% 6 - WH DMM (EPSL05) - 2% melt (calculated from BG15)
ZO(6,:) = [44.758,0.118,3.722,0.474,8.252,39.399,3.082,1.5208e-09,0.195];
% 7 - MM3 (Baker & Stolper 1994; Baker et al., 1995) 
ZO(7,:) = [45.5,0.11,3.98,0.68,7.18,38.3,3.57,0,0.31];
% 8 - Kushiro 1996 
ZO(8,:) = [43.7,0.25,2.75,.28,1.38*.9+8.81,37.22,3.26,.14,.33];
% 9 - HZ in CMAS
ZO(9,:) = [46.381,0.0,4.099,0,0,42.395,3.559,0,0];

% Choose Mantle Mode and Composition to be used.
XO = ZO(JM,:); XOW = XO;
XMS = ZMS(JS,:); XMP = ZMP(JS,:);

%Define Partition Coefficients (Note: some will be updated based on 
%   pressure below).

DMS = [0.22,4.3,0.001,0.25;0.14,5.00,0.001,0.06;0.02,1.45,0.001,0.001;...
    0.14,75.00,0.001,0.001];
DMP = [.001,.15,1.315,.625,2.70,.42;.001,.04,.20,.75...
     3.7,.19;.001,.001,.03,1.35,5.20,.03;.271,1.0,...
     1.32,.04,.05,.001;.001,.001,.018,1.230,2.70,.135]';

% Additional parameters
LH=4e5;             %Latent Heat
Cp=1e3;             %Specific heat


% Initialize dFdP array
if length(dFdP)==1
    dFdP = ones(1,length(Pb))*dFdP;
end

% Initialize melt fraction
MeltFrac = zeros(1,length(Pb));  %Array for melt fractions

for i = 1:length(Pb)
    if length(Pb) == 1
        MeltFrac(i) = dFdP;
    else
        if i == 1
            MeltFrac(i) = -(Pb(i+1)-Pb(i)).*dFdP(i);
        else
            if Pb(i)-Pb(i-1)==0
                MeltFrac(i) = .002;
            else
                MeltFrac(i) = -(Pb(i)-Pb(i-1)).*dFdP(i);
            end
        end
    end
end

% Initialize BAI, BIW, and BAI3 arrays that will hold melt compositions
BAI = zeros(1,9);   %Accumulated Melt (cation units)
BIW = zeros(1,9);   %Incremental Melt (wt%)
BAI3 = zeros(length(Pb),9);   %Incremental Melt (cation units)

% Calculate initial pyroxene modes in spinel field based on P-dependent 
%   cpx-opx join
Pyx = XMS(1) + XMS(2);
XMS(1) = (XO(7)-2.*Pyx)./(15-0.17*Pb(1));
XMS(2) = Pyx-XMS(1);

% Define mantle mode arrays
XOW_OUT = zeros(length(Pb),length(XOW));
XM_OUT = zeros(length(Pb),length(XMS));
XOW_OUT(1,:) = XOW;
XM_OUT(1,:) = XMS;

T = zeros(1,length(Pb));         %Array for corrected temperature
Ts = zeros(1,length(Pb));        %Array for solidi

%% Iteration Sequence

Frac = 0;           % Melt fraction is set to 0 at first step
melton = 0;         % Flag to determine whether melting is turned on.
if length(Pb) <= 2  % Note that usually melting is turned off to start, 
    melton = 1;     %   but for special cases with a single melt step
end                 %   it is turned on for the initial step.


CumFrac = zeros(1,length(Pb));
cumFrac = 0;
IncFrac = zeros(1,length(Pb));

% Calculate melting looping through all pressures (P)
for ii = 1:length(Pb)
    if melton == 0
        Frac = 0;
    elseif melton == 1
        Frac = MeltFrac(ii);
    end
    
    % Calculate P-dependent melting reaction [Aug, Opx, Oliv, Plag, Sp]
    YMS(1) = 0.0016*Pb(ii)^2 + 0.0006*Pb(ii) + 0.5;
    YMS(2) = -0.0022*Pb(ii)^2 - 0.0223*Pb(ii) + 0.9365;
    YMS(3) = 0.0015*Pb(ii)^2 - 0.008*Pb(ii) - 0.3414;
    YMS(4) = 0;
    YMS(5) = 0.0024*Pb(ii) + 0.0738;
    YMS = YMS./sum(YMS);     

    % Calculate bulk composition of mantle 
    for i = 1:9
        XOW(i) = (XOW(i) - Frac.*BIW(i))./(1-Frac);
    end
    %XOW(XOW<=0) = 0;
    
    %Check to ensure no oxide is exhausted in residue and reset residue
    %  and melt
    ineg = find(XOW<0);
    BIW(ineg) = BIW(ineg) + (XOW(ineg)./Frac).*(1-Frac);
    XOW(ineg) = 0;
    
    for i = 1:9
        XOW(i) = (1-Inc).*Frac.*BIW(i)+(1-(1-Inc)*Frac)*XOW(i);
    end
    XOW = 100.*XOW./sum(XOW);
       
    % Estimage Mg# of melt from Sp Lherz as function of FeO, MgO, F, Phase
    % Proportions and mineral-melt Fe-Mg Kds

    %Re-order elements in original KG92 formulation for BG15_Mgnum routine:
    %  [K,Na,Ca,Fe,Mg,Ti,Si,Al]
    iXOW = [XOW(8),XOW(9),XOW(7),XOW(5),XOW(6),XOW(2),XOW(1),XOW(3)];
       
    %For reference the original Fe-Mg KD values from KG92 are:
    % XOlK = 0.3;           %Olivine Fe-Mg Kd     (0.30)
    % XCpxK = 0.28;         %Cpx Fe-Mg Kd   (0.98*XOlK = 0.294)
    % XOpxK = 0.27;         %Opx Fe-Mg Kd   (0.96*XOlK = 0.288)
    % XSpK = 0.42;          %Spinel Fe-Mg Kd (1.52*XOlK = 0.456) 

    if Pb(ii) >= 9  %Choose between spinel and plag field
        XOlK = 0.31509+0.0062362*(Pb(ii)/10);   %Olivine Fe-Mg Kd
        XCpxK = 0.27652+0.021516*(Pb(ii)/10);   %Cpx Fe-Mg Kd 
        XOpxK = 0.27262+0.042031*(Pb(ii)/10);   %Opx Fe-Mg Kd
        %XSpK = 1.05210-0.311430*(Pb(ii)/10);    %Spinel Fe-Mg Kd 
        %  Note that given the uncertainty we chose a single value for
        %    the Spinel Fe-Mg Kd rather than a P-dependent value.
        XSpK = 0.6;    %Spinel Fe-Mg Kd 
        XMG = BG15_Mgnum(Frac,XMS,YMS,XOlK,XCpxK,XOpxK,XSpK,iXOW);
    else
        XOlK = 0.23701+0.0527690*(Pb(ii)/10);   %Olivine Fe-Mg Kd
        XCpxK = 0.24015+0.030003*(Pb(ii)/10);   %Cpx Fe-Mg Kd 
        XOpxK = 0.21347+0.039536*(Pb(ii)/10);   %Opx Fe-Mg Kd
        XSpK = 0.456; %Spinel Fe-Mg Kd (This value is not used but
                      % required as a place-holder.
        XMG = BG15_Mgnum(Frac,XMP,YMP,XOlK,XCpxK,XOpxK,XSpK,iXOW);
    end
    
    % Calculate Mg# of olivine in residual
    XMGOL = (XMG/(XOlK*(1-XMG)))./(1+(XMG/(XOlK*(1-XMG))));
    
    % Define P-dependent coefficients in D matrix
    DMS(1,1) = 22.570*(Pb(ii)^(-1.53));    % Dticpx
    DMS(1,4) = -0.440+0.557*log10(Pb(ii)); % Dnacpx
    DMS(2,1) =  0.926*(Pb(ii)^(-0.565));   % Dtiopx
    DMS(2,4) = -0.151+0.166*log10(Pb(ii)); % Dnaopx
    
    %Spinel field
    XCS = (DMS'*[XMS(1:3),XMS(5)]')'; %bulk partition coefficient in solid 
    XPS = (DMS'*[YMS(1:3),YMS(5)]')'; %bulk partition coefficient in melt
    
    %Plag field
    XCP = (DMP*XMP')';  %bulk partition coefficient in solid
    XPP = (DMP*YMP')';  %bulk partition coefficient in melt
    
    % Calculate concentration of oxide in melt (wt%)
    %Spinel
    XLFS(1) = XOW(2)/(XCS(1)+Frac*(1-XPS(1))); %Ti
    XLFS(2) = XOW(4)/(XCS(2)+Frac*(1-XPS(2))); %Cr
    XLFS(3) = XOW(8)/(XCS(3)+Frac*(1-XPS(3))); %K
    XLFS(4) = XOW(9)/(XCS(4)+Frac*(1-XPS(4))); %Na
    %Plagioclase    
    XLFP(1) = XOW(8)/(XCP(1)+Frac*(1-XPP(1))); %K
    XLFP(2) = XOW(9)/(XCP(2)+Frac*(1-XPP(2))); %Na
    XLFP(3) = XOW(7)/(XCP(3)+Frac*(1-XPP(3))); %Ca
    XLFP(4) = XOW(5)/(XCP(4)+Frac*(1-XPP(4))); %Fe
    XLFP(5) = XOW(6)/(XCP(5)+Frac*(1-XPP(5))); %Mg
    XLFP(6) = XOW(2)/(XCP(6)+Frac*(1-XPP(6))); %Ti
   
    % Caculate Ca content in melt
    %  Plag (based on Appendix A.4 KG92)
    DMP(3,1) = 1.07+0.09*XLFP(2);   %Dcacpx
    DMP(3,2) = .16+0.021*XLFP(2);   %Dcaopx
    XCP(3) = DMP(3,:)*XMP';
    XPP(3) = DMP(3,:)*YMP';
    XLFP(3) = XOW(7)/(XCP(3)+Frac*(1-XPP(3))); %Ca     
    
    %Spinel Field
    if Pb(ii)<=23
        DCA = 1.6 - 0.0196*Pb(ii)+0.071*(1-XMG);
        PCA = 24.49-0.469*Pb(ii)-12.92*(1-XMG);
        QCA = PCA/DCA;  %Ca content (XLFS)
    else % for later development with garnet melting
        DCA = 1.6 - 0.0196*23+0.071*(1-XMG);
        PCA = 24.49-0.469*23-12.92*(1-XMG);
        QCA = PCA/DCA;
    end
        
    % Calculate NaK# wt, Mg# mol, assign wt% TiO2, K2O, Na/Na+Ca wt
    if Pb(ii) >= 9
        %Spinel Field
        WNA=(XLFS(3)+XLFS(4))/(XLFS(3)+XLFS(4)+QCA); %NaK# = (Na + K) / (Na + K + Ca) 
        WNACA=XLFS(4)/(XLFS(4)+QCA); %Na# = Na / (Na + Ca)
        % assign wt.% TiO2 and K2O  
        WT=XLFS(1);
        WK=XLFS(3);
    else
        %Plag Field
        WNA=(XLFP(1)+XLFP(2))/(XLFP(1)+XLFP(2)+QCA);
        WNACA=XLFP(2)/(XLFP(2)+QCA);
        % assign wt.% TiO2 and K2O  
        WT=XLFP(6);
        WK=XLFP(1);
    end

    % Calculate mineral components and T (V. 8 MIT+CMASN)
    %   Coefficients from Till et al. [2012]
    if Pb(ii) >= 9
        %Spinel Field
        Ts(ii)=1212.+11.99*Pb(ii)-97.33*(1-XMG)-87.76*WNA+3.44*WT-4.58*WK;
        S(1)=.123+.0085*Pb(ii)+.203*(1-XMG)-.268*WNA-.019*WT+.012*WK;  %Olivine
        S(2)=.180+.0005*Pb(ii)+.075*(1-XMG)-.264*WNA+.007*WT+.012*WK;  %CPX
        S(3)=.554-.0043*Pb(ii)-.192*(1-XMG)+.949*WNA-.007*WT-.041*WK;  %Plag
        S(4)=.143-.0046*Pb(ii)-.085*(1-XMG)-.416*WNA +.020*WT+.018*WK; %Qtz
    else
        %Plag Field
        Ts(ii)=1216.+10.44*Pb(ii)-72.83*(1-XMG)-194.9*WNA+24.08*WT-1.55*WK;
        S(1)=.132+.0069*Pb(ii)+.173*(1-XMG)-.261*WNA-.009*WT+.010*WK;  %Olivine
        S(2)=.238-.0051*Pb(ii)-.065*(1-XMG)-.097*WNA+.004*WT-.013*WK;  %CPX
        S(3)=.408+.0146*Pb(ii)-.136*(1-XMG)+.630*WNA+.014*WT-.028*WK;  %Plag
        S(4)=.222-.0164*Pb(ii)-.029*(1-XMG)-.272*WNA -.009*WT+.031*WK; %Qtz
    end       
    
    
    % Save mantle mode 
    if Pb(ii) >= 9;
        XM_OUT(ii,:) = XMS;
    else
        XM_OUT(ii,:) = XMP;
    end
    XOW_OUT(ii,:) = XOW;

    
    % CHECK FOR MELTING
    T(ii) = Tb(ii) - ((cumFrac).*LH./Cp);
    if isnan(Ts(ii)) == 1;
        Ts(ii) = Ts(ii-1);
    end
    
    if T(ii)<Ts(ii)
        melton = 0;
        Frac = 0;
        BAI3(ii,:) = zeros(1,9);
        cumFrac = cumFrac+Frac;
        CumFrac(ii) = cumFrac;

    else
        melton = 1;
        Frac = MeltFrac(ii);
           
        S = S./sum(S);   % Renormalize Ol, Cpx, Plag, Qtz to 1.0
        TA = zeros(1,8);
        TA(1) = S(4);    %Qtz
        TA(3) = S(3);    %Plag
        TA(4) = S(1);    %Olivine
        TA(5) = S(2);    %Cpx

        % Recalculate to mole components
        TA(1)=TA(1);     %Qtz
        TA(2)=TA(2)/4.;  %Albite
        TA(3)=TA(3)/4.;  %Anorthite
        TA(4)=TA(4)/2.;  %Olivine
        TA(5)=TA(5)/3.;  %Cpx
        TA(6)=TA(6)/1.5; %TiO2-bearing mineral?
        TA(7)=TA(7)/4.;  %K-bearing mineral?
        TA(8)=TA(8)/6.;  %??
        TA = TA./sum(TA);

        %Recalculate albite & anorthite contents in plag component
        TPL=TA(2)+TA(3);  %plag
        XBRT=2.*WNACA/61.979;    %ratio1 of KG92
        XNRT=(1.-WNACA)/56.0794; %ratio2 of KG92
        XNA=XBRT/(XBRT+XNRT);
        TA(2)=XNA*(TPL+TA(5)); %albite 
        TA(3)=TPL-TA(2); %anorthite
        GP(9)=TA(7);     %K2O (cation frac)
        GP(10)=TA(2);    %Na2O (cation frac)

        % Transform from molar to cation fraction
        GP(3)=2.*TA(3)+TA(2)+TA(7);  %Al2O3 
        GP(8)=TA(5)+TA(3);   %CaO
        GP(2)=TA(6);  %TiO2
        GP(5)=(1.-XMG)*(2.*TA(4)+TA(6)+TA(5)); %FeO
        GP(6)=XMG*(2.*TA(4)+TA(6)+TA(5));  %MgO
        GP(1)=TA(1)+TA(4)+2.*TA(5)+2.*TA(3)+3.*TA(7)+3.*TA(2); %SiO2
        GP(4)=0.0; %Cr2O3
        GP(7)=0.0; %MnO
        GP(11)=0.0;  %Fe2O3
        GP(12)=0.0;  %P2O5
        GP = GP./sum(GP);

        % Transform from cation fraction to mole %
        GN(1)=GP(1);        %SiO2
        GN(2)=GP(2);        %TiO2
        GN(3)=0.5*GP(3);    %Al2O3
        GN(4)=GP(4);        %Cr2O3 
        GN(5)=GP(5);        %FeO
        GN(6)=GP(6);        %MgO
        GN(7)=GP(7);        %MnO
        GN(8)=GP(8);        %CaO
        GN(9)=0.5*GP(9);    %K2O
        GN(10)=0.5*GP(10);  %Na2O
        GN(11)=0.5*GP(11);  %Fe2O3
        GN(12)=0.5*GP(12);  %P2O5
        GN = 100.*GN./sum(GN);

        % Calculate composition in wt%
        CN = 100.*GN.*FK./sum(GN.*FK);

        % Add Ti, Cr, & K
        CN(2) = WT;             %Ti
        if Pb(ii)>=9
            CN(4) = XLFS(2);    %Cr
        else
            CN(4) = XLFP(2);    %Cr
        end
        CN(9) = WK;             %K
        CN = 100.*CN./sum(CN);

        % Transform to mole fraction
        GN = 100.*(CN./FK)./sum(CN./FK);

        % Recaculate molar oxides with addition of ilmenite and orthoclase
        % components
        GN(1)=GN(1)+6.*GN(9);  %SiO2 = SiO2 + 6 Or
        GN(3)=GN(3)+GN(9);     %Al2O = Al2O3 + Or
        GN(5)=GN(5)+(1.-XMG)*GN(2);  %FeO + (1-Mg#)*Ilm
        GN(6)=GN(6)+XMG*GN(2); %MgO + Mg# * Ilm
        GN = 100.*GN./sum(GN);

        %Recalculate composition in wt%
        CN = GN.*FK;
        CN = 100.*CN./sum(CN);

        % Rearrange to remove MnO, Fe2O3, P2O5 from saved melt composition
        BIW = CN;
        BIW(7) = CN(8); %cao
        BIW(8) = CN(9); %k2o
        BIW(9) = CN(10); %na2o

        % Add incremental melts
        BAI=BAI+(1-max(CumFrac))*Frac*Inc.*BIW(1:9);
        
        % Track composition at each depth
        %[SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O]
        BAI3(ii,:)=100*BIW(1:9)./sum(BIW(1:9));
        
        % Calculate cumulative melt fraction
        cumFrac = cumFrac+Frac;
        CumFrac(ii) = cumFrac;
        IncFrac(ii) = Frac;

    end

    % Calculate new phase proportions in residue based on melting reaction
    for i=1:5
        XMS(i)=(XMS(i)-Frac*YMS(i));
        XMP(i)=(XMP(i)-Frac*YMP(i));
        if XMS(i) <= 0.; XMS(i) = 0.0; end
        if XMP(i) <= 0.; XMP(i) = 0.0; end
    end
    XMS = XMS./sum(XMS);
    XMP = XMP./sum(XMP);

    % Convert remaining Sp to Plag at P = 9 kbar
    if Pb(ii) == 9
        XMP(4) = 1.5*XMS(5);
        XMP(5) = 0;
        XMP(1) = XMS(1) - XMS(5);
        XMP(2) = XMS(2) - XMS(5);
        XMP(3) = XMS(3) + 1.5*XMS(5);
    end
    % If plag exhausted correct melting reaction
    if XMP(4) == 0
        YMP(1) = 0.66;
        YMP(2) = -0.02;
        YMP(3) = 0.36;
        YMP(4) = 0;
        YMP(5) = 0;
    end 
end

BAI = 100.*BAI./sum(BAI);
Fc = max(CumFrac);
iFc = find(CumFrac == Fc);
if length(iFc) > 1
    CumFrac(iFc(2:end)) = 0;
end

%% Additional Subroutines
function [RT4BIS] = BG15_Mgnum(FL,XM,YM,XOlK,XCpxK,XOpxK,XSpK,XO)

%==========================================================================
% This is a modified version of the rtbis routine found in
% Numerical Recipes.  Its purpose is to be an elegant root finder
% This version is modified to estimate the mg# of melt from mantle
% assumed to consist of pl or sp lherz as function of F (extent of melt)
% and depends on ol/melt KD and cpx/melt, opx/melt, and sp/melt KDs
% and phase proportions in mantle and on bulk mgo and feo
% 
% FL is melt fraction (WT), X1 and X2 are guesses of mgnum of melt, XACC is
% tolerance, XOlK, XCpxK, XOpxK, XSpK are the ol/melt cpx/melt, opx/melt,
% and sp/melt KDs, repsectively, XM(5) is vector of phase proportions (WT)
% MMGNUM is function that solves for melt mgnum
% 
% This version of the RT4BIS subroutine from Kinzler & Grove 92 takes XO,
%   XM, YM, and FL in weight fractions.  Note that the order of elements in 
%   XO is [K,Na,Ca,Fe,Mg,Ti,Si,Al] as in the original KG92 routines.
% An important component are the assumptions that go into the ol,px,&sp
%   composition in FUNC1-FUNC4 in the MMGNUM subroutine.
% This routine is valid in both the plag & spinel fields because plag will
%   not affect the Mg-Fe partitioning (i.e., XMM(4) does not come into the
%   final calculation of Mg#)
% M. Behn & T. Grove 11/13/12
% M. Behn adjusted Ca/Mg+Fe for cpx from 25 (KG92) to 40 (K97) 11/14
%==========================================================================

XACC = 1e-5;
X1 = 0.80;
X2 = 0.68;
JMAX=40;

% Use melt fraction (FL), initial phase proportions (XM) and melt 
% reaction (YM) to estimate phase proportions at end of melting (all 
% in wt. units)

SUM=0.0;
for I = 1:5
    XMN(I)=XM(I)-FL*YM(I);
    if XMN(I) <= 0.0; XMN(I)=0.0; end
    SUM=SUM+XMN(I);
end
for I = 1:5
    XMN(I)=XMN(I)/SUM;
end

% convert wt. fraction FL to cation fraction FLM
FLM=0.0002+0.926*FL;

% convert wt.% FeO (XO(4)) and MgO (XO(5)) in mantle to CATION %

XFEOM=100*(XO(4)/71.846)/((2*XO(1)/94.201)+(2*XO(2)/61.979)+...
    (XO(3)/56.794)+(XO(4)/71.846)+(XO(5)/40.311)+(XO(6)/79.89)...
    +(XO(7)/60.0848)+(2*XO(8)/101.961));
XMGOM=100*(XO(5)/40.311)/((2*XO(1)/94.201)+(2*XO(2)/61.979)+...
    (XO(3)/56.794)+(XO(4)/71.846)+(XO(5)/40.311)+(XO(6)/79.89)...
    +(XO(7)/60.0848)+(2*XO(8)/101.961));

F=MMGNUM(X1,FLM,XMGOM,XFEOM,XMN,XOlK,XCpxK,XOpxK,XSpK);

if F < 0
    RT4BIS=X1;
    DX=X2-X1;
else
    RT4BIS=X2;
    DX=X1-X2;
end
for J = 1:JMAX
    DX=DX*.5;
    XMID=RT4BIS+DX;
    FMID=MMGNUM(XMID,FLM,XMGOM,XFEOM,XMN,XOlK,XCpxK,XOpxK,XSpK);
    if FMID < 0; RT4BIS=XMID; end
    if abs(DX) < XACC | FMID == 0; return; end
end

function FUNC = MMGNUM(MGNUM,FLM,XMGOM,XFEOM,XMN,XOlK,XCpxK,XOpxK,XSpK);
% Mg number of melt from LHERZOLITE mantle function
% MGNUM is the Mg# of MELT
% FLM is melt fraction (CATION FRACTION)
% XOlK,XCpxK,XOpxK,XSpK are KDs for oliv/melt, cpx/melt, opx/melt, &
%    sp/melt
% XM is vector containing phase proportions  (WT. FRACTION)
%  at end of FL melt production
% XMGOM is MGO (CATION MOLE FRACTION) IN bulk mantle prior to melting
% XFEOM IS FEO (CATION MOLE FRACTION) IN bulk mantle prior to melting

MGNM(1)=MGNUM/(MGNUM+XCpxK*(1-MGNUM));
MGNM(2)=MGNUM/(MGNUM+XOpxK*(1-MGNUM));
MGNM(3)=MGNUM/(MGNUM+XOlK*(1-MGNUM));
MGNM(5)=MGNUM/(MGNUM+XSpK*(1-MGNUM));
% Divison by 1.5 below is to correct everything to 4 oxygens (cpx & opx
% both have 6 oxygens).
SUM=XMN(1)*100/(MGNM(1)*216.56/1.5+(1-MGNM(1))*248.095/1.5)+...
    XMN(2)*100/(MGNM(2)*200.792/1.5+(1-MGNM(2))*263.862/1.5)+...
    XMN(3)*100/(MGNM(3)*140.707+(1-MGNM(3))*203.777)+...
    XMN(4)*100/275.012+XMN(5)*100/(MGNM(5)*142.272+(1-MGNM(5))*173.807);
XMM(1)=(XMN(1)*100/(MGNM(1)*216.56/1.5+(1-MGNM(1))*248.095/1.5)/SUM);
XMM(2)=(XMN(2)*100/(MGNM(2)*200.792/1.5+(1-MGNM(2))*263.862/1.5)/SUM);
XMM(3)=(XMN(3)*100/(MGNM(3)*140.707+(1-MGNM(3))*203.777)/SUM);
XMM(4)=(XMN(4)*100/275.012/SUM);
XMM(5)=(XMN(5)*100/(MGNM(5)*142.272+(1-MGNM(5))*173.807)/SUM);
% Note that plagioclase (here the 4th component XMM(4)) does not have Fe or
% Mg its structure -- so it does not factor into the calculation of Mg#
% below.
FUNC1=XMM(3)*(66.7*MGNUM/(MGNUM+XOlK*(1-MGNUM)));
FUNC2=XMM(2)*(50.*MGNUM/(MGNUM+XOpxK*(1-MGNUM)));
FUNC3=XMM(1)*(40.*MGNUM/(MGNUM+XCpxK*(1-MGNUM)));
FUNC4=XMM(5)*(33.33*MGNUM/(MGNUM+XSpK*(1-MGNUM)));
FUNC=XMGOM-(1-FLM)*(FUNC1+FUNC2+FUNC3+FUNC4);
FUNC=FUNC/(XMGOM+XFEOM-(1-FLM)*(66.67*XMM(3)+50.0*XMM(2)...
    +40.0*XMM(1)+33.33*XMM(5)));
FUNC=FUNC-MGNUM;
