%==========================================================================
% Example script to generate the isobaric melt calculations shown in 
%   Figure 4 of Behn & Grove [JGR, 2015].
% Note different values of dF/dP are used to produce different melt steps 
%   going from 1% to 20% melting for the different pressures of melting.
%
%  Mark Behn, January 2015
%==========================================================================


%% Calculate Isobaric Melting at Pressures of 5, 10, 15, and 20 kbar

Pinit = 5; %5 kbar 
Pfin = 4.9;
Pinc = .1;
Pb = Pinit:-Pinc:Pfin;
Tb = 1650*ones(1,length(Pb)); %Note temperature is set high to ensure 
                              %  melting; melt step is set by dFdP*Pinc

dFdP = [.1,.5,1,1.5,2];

for i = 1:length(dFdP)
    [XM,XOW,BAI,BAI3,Ts,T,CumFrac,IncFrac,Fc] = ...
        BG15(Tb,Pb,dFdP(i),2,2,0.0); 
        BG15_batch_5kb(i,:) = BAI3(1,:);
end    


Pinit = 10; %10 kbar 
Pfin = 9.9;
Pinc = .1;
Pb = Pinit:-Pinc:Pfin;
Tb = 1650*ones(1,length(Pb)); 

dFdP = [.1,.5,1,1.5,2];

for i = 1:length(dFdP)
    [XM,XOW,BAI,BAI3,Ts,T,CumFrac,IncFrac,Fc] = ...
        BG15(Tb,Pb,dFdP(i),2,2,0.0); 
        BG15_batch_10kb(i,:) = BAI3(1,:);
end    

Pinit = 15; %15 kbar 
Pfin = 14.9;
Pinc = .1;
Pb = Pinit:-Pinc:Pfin;

dFdP = [.1,.5,1,1.5,2];

for i = 1:length(dFdP)
    [XM,XOW,BAI,BAI3,Ts,T,CumFrac,IncFrac,Fc] = ...
        BG15(Tb,Pb,dFdP(i),2,2,0.0); 
        BG15_batch_15kb(i,:) = BAI3(1,:);
end    


Pinit = 20; %20 kbar 
Pfin = 19.9;
Pinc = .1;
Pb = Pinit:-Pinc:Pfin;

dFdP = [.1,.5,1,1.5,2];

for i = 1:length(dFdP)
    [XM,XOW,BAI,BAI3,Ts,T,CumFrac,IncFrac,Fc] = ...
        BG15(Tb,Pb,dFdP(i),2,2,0.0); 
        BG15_batch_20kb(i,:) = BAI3(1,:);
end    

%% Plot isobaric results
figure(1)
clf

subplot(221)
plot(BG15_batch_5kb(:,5),BG15_batch_5kb(:,9),'b-o','markerfacecolor','w')
hold on
plot(BG15_batch_10kb(:,5),BG15_batch_10kb(:,9),'b-d','markerfacecolor','b')
plot(BG15_batch_15kb(:,5),BG15_batch_15kb(:,9),'b-^','markerfacecolor','b')
plot(BG15_batch_20kb(:,5),BG15_batch_20kb(:,9),'b-s','markerfacecolor','b')

xlabel('FeO (wt%)')
ylabel('Na_2O (wt%)')
legend('5 kbar','10 kbar','15 kbar','20 kbar')
axis([3 12 1 7.5])

subplot(222)
plot(BG15_batch_5kb(:,1),BG15_batch_5kb(:,9),'b-o','markerfacecolor','w')
hold on
plot(BG15_batch_10kb(:,1),BG15_batch_10kb(:,9),'b-d','markerfacecolor','b')
plot(BG15_batch_15kb(:,1),BG15_batch_15kb(:,9),'b-^','markerfacecolor','b')
plot(BG15_batch_20kb(:,1),BG15_batch_20kb(:,9),'b-s','markerfacecolor','b')

xlabel('SiO_2 (wt%)')
ylabel('Na_2O (wt%)')
axis([40 56 1 7.5])

subplot(223)
plot(BG15_batch_5kb(:,5),BG15_batch_5kb(:,7)./BG15_batch_5kb(:,3),...
    'b-o','markerfacecolor','w')
hold on
plot(BG15_batch_10kb(:,5),BG15_batch_10kb(:,7)./BG15_batch_10kb(:,3),...
    'b-d','markerfacecolor','b')
plot(BG15_batch_15kb(:,5),BG15_batch_15kb(:,7)./BG15_batch_15kb(:,3),...
    'b-^','markerfacecolor','b')
plot(BG15_batch_20kb(:,5),BG15_batch_20kb(:,7)./BG15_batch_20kb(:,3),...
    'b-s','markerfacecolor','b')

xlabel('FeO (wt%)')
ylabel('CaO/Al_2O_3')
axis([3 12 0.2 1])

