%% Porenkompressibilität aus Volumenänderung
%% ----------------------------------------------------------------------
clc; clear all;
%% Inputdaten bereitstellen (besser automatisiert aus Datei lesen)
%  Druck   
p=[4.999;12.07;13.62;15.13;19.13;23.259];
% Pipette oder Quizixwerte
mvp=[18.3896;18.3323;18.3166;18.2986;18.2433;18.1926];
% Anfangsplugvolumen
vp0=22.961;
% Sample-Name (wegen Titel und Bezeichnung des PS-Files)
sample='NN';
% -----------------------------------------------------------------------                    
% Meßwert bei p=0 (Extrapolation des Splines)
mvp0=spline(p,mvp,0);
% Meßwerte (Druck, Pipette) um diese Angaben (vorne) ergänzen
p=[0;p]; mvp=[mvp0;mvp];
% Porenvolumen aus Pipettenwerten (vp(0)=vp0)
vp=vp0-(mvp0-mvp);
% Porenvolumen interpolieren 
% Die Interpolation würde auch mit yy=spline(p,vp,xx) gehen. die pp-Form
% benötigt man aber zur Differentiation!
pp = csape(p,vp);
x=p(1):0.2:p(end); x=x';
y=ppval(pp,x);
% Porenvolumen differenzieren
dpp=fnder(pp);
% Porenkompressibilität (an den Meßpunkten)
dvp=ppval(dpp,p);
pk=-dvp/vp0;
% Porenkompressibilität (für alle x)
dvp=ppval(dpp,x);
yy=-dvp/vp0;
%% Grafik
%%------------------------------------------------------------------------
clf
FN='Times New Roman'; FS=13;
% Meßdaten und Interpolation
subplot(2,1,1);
plot(p,vp,'bo'); hold on;
plot(x,y,'r-')
xlabel('pressure p [MPa]','FontSize',FS,'LineWidth',0.5,'FontName',FN);
ylabel('pore volume V_p [ml]','FontSize',FS,'LineWidth',0.5,'FontName',FN);
set(gca,'box','on','FontSize',FS,'LineWidth',0.5,'FontName',FN);
title(sample,'FontSize',FS+2,'FontName',FN);
% Porenkompressibilität
subplot(2,1,2);
plot(p,pk,'bo'); hold on;
h=plot(x,yy,'r-');
hh=legend(h,'$$-\frac{1}{V_p(0)}\,\frac{\partial V_p(p)}{\partial p}$$');
set(hh,'interpreter','latex')
xlabel('pressure p [MPa]','FontSize',FS,'LineWidth',0.5,'FontName',FN);
ylabel('pore compressibility [1/MPa]','FontSize',FS,'LineWidth',0.5,'FontName',FN);
set(gca,'box','on','FontSize',FS,'LineWidth',0.5,'FontName',FN);
%% Textdatei mit daten an den Meßpunkten schreiben
sample = strrep(sample, '\_', '_');
fid=fopen([sample '.txt'],'w');
fprintf(fid,'pressure   measured data   pore volume   pore compressibility\r\n');
fprintf(fid,'%6.2f %14.5f %14.5f %20.8f \r\n',[p,mvp,vp,pk]');
fclose(fid);
%% PS-File mit Grafik erzeugen
set(gcf,'PaperOrientation','portrait',...
   'PaperUnits','centimeter',...
   'PaperPosition',[1.5,3,18,23]);
print ('-deps2c',[sample '.ps']);