%=========================================================================
% Example script to generate the fractional melting calculations for
%   different mantle adiabats (1350 & 1450C; blue & red respectively) 
%   shown in Figures 7a,b of Behn & Grove [JGR, 2015].  
%   dF/dP is fixed to 1%/kbar.  Large symbols show pooled melts; small 
%   symbols show incremental melts plotted every 1 kbar.
%
%  Mark Behn, January 2015
%==========================================================================


clear

%% Caculate Fractional Melting Along 1350 and 1450ºC adiabats
Pinit = 30;
Pfin = 3;
Pinc = .1;
Pb = Pinit:-Pinc:Pfin;
Tb = 1350*ones(1,length(Pb))+1.5*Pb; %T + rough adiabatic gradient

dFdP = .01;  %Note you can add variable dFdP if you like -- just make sure
             %  that the dFdP vector is the same length as Pb & Tb.

[XM,XOW,BAI,BAI3,Ts,T,CumFrac,IncFrac,Fc] = BG15(Tb,Pb,dFdP,2,2,0.9); 
BAI3_BG15_1350 = BAI3(2:end,:);  %Be careful of the first melt it can be a bit weird.
BAI_BG15_1350 = mean(BAI3_BG15_1350(BAI3_BG15_1350(:,1)>0.01,:));  %Mean melt


Pinit = 30;
Pfin = 3;
Pinc = .1;
Pb2 = Pinit:-Pinc:Pfin;
Tb = 1450*ones(1,length(Pb2))+1.5*Pb2;

[XM,XOW,BAI,BAI3,Ts,T,CumFrac,IncFrac,Fc] = BG15(Tb,Pb2,dFdP,2,2,0.9); 
BAI3_BG15_1450 = BAI3(2:end,:);
BAI_BG15_1450 = mean(BAI3_BG15_1450(BAI3_BG15_1450(:,1)>0.01,:));


%% Plot fractional melting results

figure(1)
clf

subplot(221)
plot(BAI3_BG15_1350(8:10:end,5),BAI3_BG15_1350(8:10:end,9),'b^','markerfacecolor','w')
hold on
plot(BAI_BG15_1350(1,5),BAI_BG15_1350(1,9),'b^','markerfacecolor','b','markersize',15)

xlabel('FeO (wt%)')
ylabel('Na_2O (wt%)')

text(8.1,5.25,'T_p = 1350^oC','fontsize',12)

grid on
axis([6 9.5 0 6])

subplot(222)
plot(BAI3_BG15_1350(8:10:end,1),BAI3_BG15_1350(8:10:end,9),'b^','markerfacecolor','w')
hold on
plot(BAI_BG15_1350(1,1),BAI_BG15_1350(1,9),'b^','markerfacecolor','b','markersize',15)


xlabel('SiO_2 (wt%)')
ylabel('Na_2O (wt%)')

grid on
axis([45 53 0 6])


subplot(223)
plot(BAI3_BG15_1450(7:10:end,5),BAI3_BG15_1450(7:10:end,9),'bo','markerfacecolor','w')
hold on
plot(BAI_BG15_1450(1,5),BAI_BG15_1450(1,9),'bo','markerfacecolor','b','markersize',15)

xlabel('FeO (wt%)')
ylabel('Na_2O (wt%)')

text(8.1,5.25,'T_p = 1450^oC','fontsize',12)

grid on
axis([6 9.5 0 6])

subplot(224)
plot(BAI3_BG15_1450(7:10:end,1),BAI3_BG15_1450(7:10:end,9),'bo','markerfacecolor','w')
hold on
plot(BAI_BG15_1450(1,1),BAI_BG15_1450(1,9),'bo','markerfacecolor','b','markersize',15)


xlabel('SiO_2 (wt%)')
ylabel('Na_2O (wt%)')

grid on
axis([45 53 0 6])

orient tall

