%%%%%%%%%%%%%%% calculating saturations and gas released and salinity during
%%%%%%%%%%%%%%% depressurization from known initial saturations and
%%%%%%%%%%%%%%% salinity
clc
clear
%% =====================Setting up P,T,S diagram===========================
%%%%%%%%%%%%% sodium chloride %%%%%%%%%%%%%%%%%%%%
% load p_salinity_agu.mat
% T_salinity=(0:0.25:20);
% salinity=(1:0.5:18.5)/100;

%%%%%%%%%%%%%%% sodium chloride %%%%%%%%%%%%%%%%%%%%
load p_salinity_Kehua_Liu_NaCl.mat
T_salinity=(0:0.2:20);
salinity=(0:0.5:40)/100;
%% ==============calculate initial or constant parameters=================================
T_room=273.15+24.73;    % room temperature
P_room=1.01325e5;  % room pressure
R=8.314;           % ideal gas constant
% T=2.15;
Pi=12e6;%8.6e6;        % the initial pressure in Pa
Ti=3; %0;            % the initial temperature in degree centigrade

% P=(3.21e6:0.01e6:6.94e6);         % the final pressure in Pa
P=(12e6:-0.01e6:0e6);%(8.6e6:-0.01e6:0e6);
T=1;%0;              % the final temperature in degree centigrade

N=5.75;          % the average number of methane in one hydrate, dimensionless
Mw=0.018;        % the molar weight of water, the salt or the methane in water only slightly change this number, so we neglect the influence of salt and methane, kg/mol
Mm=0.016;        % the molar weight of methane, kg/mol
Mh=Mm+N*Mw;      % the molar weight of hydrate, kg/mol

cli=0.07;%3188;%0.105;%0.1633;       % the initial mass fraction of salt in water, kg/kg
% cmi=methane_solubility(Pi/10^6, Ti,2*cli/0.11099/(1-cli))*0.016;      % the initial  methane molality in water, mol/kg
% denwi=brine_density(Pi/10^6,Ti,2*cli/0.11099/(1-cli),cmi/0.016);      % initial water density NaCl
cmi=methane_solubility(Pi/10^6, Ti,cli/0.1028938/(1-cli))*0.016; 
denwi=brine_density(Pi/10^6,Ti,cli/0.1028938/(1-cli),cmi/0.016);      % initial water density NaBr
denh=912;        % hydrate density, kg/m^3
dengi=gas_density(Pi/10^6,Ti);  % the gas density under initial temperature, pressure condition, kg/m^3

Shi=0.80;    % initial hydrate saturation
Sgi=0.10;    % initial gas saturation
Swi=1-Shi-Sgi;  % initial water saturation
V_pore=100e-6; % the total volume of the pore in the sample

%load experiment data
expData=load('experiment.csv');
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Num=length(P);
cle=zeros(Num,1);
Swf=zeros(Num,1);
Shf=zeros(Num,1);
Sgf=zeros(Num,1);
dn=zeros(Num,1);
V=zeros(Num,1);
for i=1:Num
    cle(i)=interpolation2(p_salinity,T_salinity,salinity,P(i)./1e6,T);    % the three phase equilibrium salt mass fraction under pressure P, temperature T, kg/kg
    %%%%%%%%%%%%%%%%%% sodium chloride %%%%%%%%%%%%%%%%%%%%%
%     cmf=hydrate_solubility(P(i)/10^6, T,2*cle(i)/0.11099/(1-cle(i)))*0.016;        % final  methane molality in water, mol/kg
%     denwf=brine_density(P(i)/10^6,T,2*cle(i)/0.11099/(1-cle(i)),cmf/0.016); % NaCl
   %%%%%%%%%%%%%%% sodium bromide %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    cmf=hydrate_solubility(P(i)/10^6, T,cle(i)/0.1028938/(1-cle(i)))*0.016;  
    denwf=brine_density(P(i)/10^6,T,cle(i)/0.1028938/(1-cle(i)),cmf/0.016); % NaBr
    
    dengf=gas_density(P(i)/10^6,T);  % gas density at P, T
    
    
    Swf(i)=Swi*denwi*cli/(denwf*cle(i));           % water saturation at P, T , from mass balance for salt 
    Shf(i)=(Swi*denwi*(1-cli)*(1-cmi)+Shi*denh*(N*Mw)/Mh-Swf(i)*denwf*(1-cle(i))*(1-cmf))*Mh/(N*Mw*denh);  % hydrate saturation at P, T, from mass  balance for H2O
    Sgf(i)=1-Swf(i)-Shf(i);   % gas saturation at P, T
    
    
    %%%%%%%%%%%%% if initial salinity is less than the 3-phase salinity, no
    %%%%%%%%%%%%% hydrate melting, liquid, hydrate, gas saturation does not
    %%%%%%%%%%%%% change, but liquid, gas densities change, it is gas
    %%%%%%%%%%%%% expansion
   if cle(i)>=cli
       Swf(i)=Swi;
        Shf(i)=Shi;
        Sgf(i)=Sgi;     
        cle(i)=cli;
   end
    
   
   %%%%%%%%%%%% if hydrate saturation at P, T is less than zero, this means
   %%%%%%%%%%%% all hydrate has been melted, hydrate saturation should be
   %%%%%%%%%%%% zero, liquid and gas saturation does not change, but their
   %%%%%%%%%%%% density could change, it is gas expansion
    if Shf(i)<0
        Swf(i)=Swf(i-1);
        Shf(i)=0;
        Sgf(i)=Sgf(i-1);     
        cle(i)=cle(i-1);
    end
    
    %%%%%%%%%%%%%%%%%%% molar or volume of collected gas at room pressure
    %%%%%%%%%%%%%%%%%%% and temperature, from mass balance  for CH4
    dn(i)=Swi*denwi*cmi/Mm+Shi*denh/Mh+Sgi*dengi/Mm-Shf(i)*denh/Mh-Sgf(i)*dengf/Mm-Swf(i)*denwf*cmf/Mm;
    V(i)=dn(i)*R*T_room/P_room;
    
end

%%%%%%%%%%%%%%%%% plot data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% hold on
% subplot(2,2,1)
% hold on
% plot(Shf,P./1e6,'r','LineWidth',4)
% xlabel('Hydrate saturation','FontSize',30)
% ylabel('Pressure, MPa','FontSize',30)
% ylim([0 8])
% set(gca,'FontSize',30,'LineWidth',2)
% box on
% 
% hold on
% subplot(2,2,2)
% hold on
% plot(Swf,P./1e6,'b','LineWidth',4)
% xlabel('Water saturation','FontSize',30)
% ylabel('Pressure, MPa','FontSize',30)
% set(gca,'FontSize',30,'LineWidth',2)
% ylim([0 8])
% box on
% 
% 
% hold on
% subplot(2,2,3)
% hold on

% 
% 
% 
% hold on
% subplot(2,2,4)
% hold on
% plot(V*V_pore,P./1e6,'r','LineWidth',4)
% xlabel('Cummulative volume of gas released, mL','FontSize',30)
% ylabel('Pressure, MPa','FontSize',30)
% ylim([0 8])
% set(gca,'FontSize',30,'LineWidth',2)
% box on

figure(1)
% subplot(2,1,1);
hold on
plot(dn*V_pore*1000,P./1e6,'r-','LineWidth',1.5,'MarkerSize',5)    
xlabel('Cumulative CH_4 (mmol)','FontSize',16)
ylabel('Pressure (MPa)','FontSize',16)
ylim([0 14])
xlim([0 1000])
set(gca,'XTick',(0:50:1000),'YTick',(0:2:14))
set(gca,'FontSize',16,'LineWidth',1.5)
box on
grid off
legend('Equilibrium model for NaCl at 1 \circC')

mmolOut=dn*V_pore*1000;
Pout=P'./1e6;

% subplot(2,1,2);
% plot(dn*V_pore*1000,cle.*100,'b','LineWidth',1.5)
% xlabel('Salinity,wt.%','FontSize',16)
% ylabel('Salinity (wt.%)','FontSize',16)
% set(gca,'FontSize',16,'LineWidth',1.5)
% ylim([6 11])
% xlim([0 300])
% box on