format compact
clear;clc;

atm=1; %1 atm
cp=4.186; %KJ/kg/deg

load('flow2.dat')
Aflow=flow2(1);
Tylo=flow2(2);
Tylowet=flow2(3);
Tyhi=flow2(4);
Txhi=flow2(5);
Txlo=flow2(6);
Wflow=flow2(7);

% Tylo=23.149;
% Tylowet=14.949;
% Tyhi=30.621;
% Txhi=35.339;
% Txlo=25.505;
% Aflow=2.549;
% Wflow=12;

[hr,rh,enthal,dp,vh]=psychro2(Tylo,Tylowet,atm);
Gyprime=Aflow/vh;
Gx=Wflow*0.995;
Hy(1,:)=ones(1,3)*enthal;

Tx(1,:)=ones(1,3)*Txlo;
[hr,enthal,dp,wb,vh]=psychro(Txlo,100,atm);
Hstar(1,:)=ones(1,3)*enthal;

htulo=0.01;
htuhi=5;
htumid=(htulo+htuhi)/2;
htu=[htulo htumid htuhi];
num=10;
dL=1/num;

dT=[100 100 100];
while abs(dT(2)) > 0.01
%for k=1:50
    for i=2:num
        dHdL=htu.*(Hstar(i-1,:)-Hy(i-1,:));
        dTdL=dHdL*Gyprime/Gx/cp;
        Tx(i,:)=Tx(i-1,:)+dTdL*dL;
        for j=1:3
            [hr,enthal,dp,wb,vh]=psychro(Tx(i,j),100,atm);
            Hstar(i,j)=enthal;
            Hy(i,j)=Gx*cp/Gyprime*(Tx(i,j)-Tx(i-1,j))+Hy(i-1,j);
        end
    end
    dT=Tx(num,:)-Txhi
    if dT(1)*dT(2)<0
        htuhi=htumid;
    elseif dT(2)*dT(3)<0
        htulo=htumid;
    end
    htumid=(htuhi+htulo)/2;
    htu=[htulo htumid htuhi];   
    figure(5)
    plot(Tx)
    hold on
    plot(num,Txhi,'ok')
    hold off
    pause;
end
htu