%depress


% Porosity
porosity=0.389;
% Pore volume (cm3)
V_pore=92.92201192;
% Initial brine salinity (prior to methane)
S=0.00035;


% Initial hydrate saturation
Sh=0.13;
% Initial volume of hydrate (cm3)
Vi_hyd=Sh*V_pore;
% Initial mass of hydrate (g)
gi_hyd=Vi_hyd*0.91;
% Initial moles of hydrate (mol)
moli_hyd=0.09101;
% Initial water saturation
Sw = 0.646125795;
% Initial water volume
Vi_H2O=Sw*V_pore;

% Initial salinity 
%Si=S*(Vi_H2O+(108*moli_hyd))/Vi_H2O;
Si=0.000402;


% Load salinity
load p_salinity_Kehua_Liu_NaCl.mat
% T_salinity=(0:0.2:20);
% salinity=(0:0.5:40)/100;
% T=(0:0.2:20);
% N=length(T);

T_salinity=(0:0.2:40);
salinity=(0:0.5:40)/100;
T=(0:0.2:40);
N=length(T)

% Initial equilibrium pressure
%P_eq_init=interpolation_pressure(p_salinity,T_salinity,salinity,Si,1);
P_eq_init=interpolation_pressure(p_salinity,T_salinity,salinity,Si,1)

% Stepwise dissociation
% dissociate over 50 steps
step=50;
CH4_step=moli_hyd/50;

for i=1:step
    mol_hyd(i)=moli_hyd-(i*CH4_step);   % Remaining mol of hydrate
    cum_mol_release(i)=i*CH4_step;      % Cumulative mol of hydrate CH4 release
    gH2O_release(i)=cum_mol_release(i)*108; % Water release g
    V_H2O(i)=Vi_H2O+gH2O_release(i);
    S(i)=Si/(V_H2O(i)/Vi_H2O);
end


for x=1:length(S)
   Peq(x)=interpolation_pressure(p_salinity,T_salinity,salinity,S(x),1);
end

subplot(2,1,1)
plot(cum_mol_release,S)
xlabel('Cumulative hydrate methane(mol)')
ylabel('Salinity (wt.%)')
subplot(2,1,2)
plot(cum_mol_release,Peq)
xlabel('Cumulative hydrate methane(mol)')
ylabel('Pressure (MPa)')






%plot(salinity,PV,'k-')

% figure(1)
% plot(T,P12,'b-',T,P7,'r-',T,P6,'g-',T,P0,'k-')
% ylabel('Pressure (MPa)','FontSize',15)
% xlabel('Temperature (\circC)','FontSize',15)
% ylim([1 5])
% xlim([0 2])
% set(gca,'XTick',(0:0.5:4),'YTick',(0:1:10))
% 
% figure(2)
% plot(T,P12,'b-',T,P7,'r-',T,P6,'g-',T,P0,'k-')
% ylabel('Pressure (MPa)','FontSize',15)
% xlabel('Temperature (\circC)','FontSize',15)
% ylim([0 20])
% xlim([0 20])
% set(gca,'XTick',(0:2:20),'YTick',(0:2:20))