function isotest
clear
clc
close all
p=1e5 ;
Mg=1;
Mf=1;
Tg=273.15+300;
Tf=273.15+20;
A=0.1;
lambda=0.58;

Sf=Mf*XSteam('s_pT',p/1e5,Tf-273.15)*1000;
Sg=Mg*XSteam('s_pT',p/1e5,Tg-273.15)*1000;
tend=3600*24*10;
[T,Y]=ode45(@derivation,[0 tend],[Mf Mg Sf Sg]);
T=T/3600/24;

function dydt=derivation(t,y)
dydt=zeros(4,1);
Mf=y(1);
Mg=y(2);
Sf=y(3); sf=Sf/Mf;
Sg=y(4); sg=Sg/Mg;
rhof=XSteam('rho_ps',p/1e5,sf/1000);
rhog=XSteam('rho_ps',p/1e5,sg/1000);
delta=(Mf/rhof+Mg/rhog)/2/A;
k=lambda/delta;

%% Wie bekomme ich Tf und Tg als Vektoren in den Workspace?
Tf=XSteam('T_ps',p/1e5,sf/1000)+273.15;
Tg=XSteam('T_ps',p/1e5,sg/1000)+273.15;
assignin('caller','Temp',Tg)

plot(t/3600/24,Tf,'r*')
xlabel('Zeit [Tage]')
ylabel('Temperatur [K]')
hold on
% Temp(1)=Tf
% assignin('base','Temp',Tf)
dotQ=k*A*(Tg-Tf);Tf-273.15
hf=1000*XSteam('h_pT',p/1e5,20);
hg=1000*XSteam('h_pT',p/1e5,200);
dotm=dotQ/(hg-hf);
dydt(1)=dotm;
dydt(2)=-dotm;
if dotm>0
dotS=dotm*sg; % Kondensationsfall
else
dotS=dotm*sf; % Verdampfungsfall
end
dydt(3)=dotS+dotQ/Tf ;
dydt(4)=-dotS-dotQ/Tg ;

end

end