%Calcula la tª de burbuja en ºC
%F={Caudal(kmol/h),composición(5x1),Presión(bar),Temperatura(ºC)}

function Tb=Tb(F)
[CF,z,PF,TF]=deal(F{:});
z=z/sum(z);
z=z';
Tc=TF;%ºC
yc=K(PF,Tc,z,z).*z;
h=0.01;
DP=5;

% Para asegurarnos que calculamos el coeficente de fugacidad del vapor
% cuando la PF es alta vamos incrementando poco a poco la presión para 
% obtener buenos valores iniciales de la composición
% Pp es la presión que se prueba Pc  a la que se ha conseguido la convergencia

Pc=min(0.5*PF,1);
Pp=Pc;%.Pc
is=0;
while Pc ~= PF
   % Bucle externo de composición
   
   ype=yc;
   for i = 1:30 %Se fija un número máximo de iteraciones 
      % yp es la composición de prueba, ype la salida del bucle de temperatura
      % e yc la calculada satisfactoriamente
      yp=ype/sum(ype);
      % Bucle interno de temperatura
      Tp=Tc;
      a=0.5;
      for j=1:30
         [k,iz]=K(Pp,Tp,z,yp);
   		[k1,iz1]=K(Pp,Tp+h,z,yp);
         [k2,iz2]=K(Pp,Tp,z,yp);
         Tp=Tp-a*log(sum(k.*z))*2*h/(log(sum(k1.*z))-log(sum(k2.*z))); 
         %Fijamos dos cotas para la temperatura. Si se superan se aumenta la amortiguación
   		if Tp>=290
      		% error('T>=290,Tªlimite')
      		Tp=290;
        	 	a=a/2;
     		else 
         	if Tp<-100;
            	Tp=-100;
            	a=a/2;
         	end     
         end 
         %--------------------------------
         
          % Iteración de presión
      if a < 0.01
         a2=1;
         
         Tp=Tc+1;
         yp=ype/sum(ype);
         for j1=1:30
            [k,iz]=k(Pp,Tp,z,yp);
   			[k1,iz1]=k(Pp+h,Tp,z,yp);
            [k2,iz2]=k(Pp-h,Tp,z,yp);
            Pp=Pp-a2*log(sum(k.*z))*2*h/(log(sum(k1.*z))-log(sum(k2.*z)));
            if Pp>50
               Tp=Tp-10;
               a2=a2/1.1;
            end
            
            ype=K(Pp,Tp,z,yp).*z;
            errP=abs(log(sum(ype)));
         	if errP < 0.001
            	break;
            	else if j1==30
               'Máximo número de iteraciones sin converger P'
            	end
         	end   
         end % for j1
       end

         %--------------------------------
         
         
      	ype=K(Pp,Tp,z,yp).*z;
         errT=abs(log(sum(ype)));
         if errT < 0.001
            break;
            else if j==30
               disp('Máximo número de iteraciones sin converger T');
            end
         end   
      end % for j
      %Comprobamos la convergencia de composiciones
      erry=sum(abs(ype - yp));
      if erry < 0.001 
                  break;
        else if i==30
               disp('Máximo número de iteraciones sin converger y');
        end
      end   
   end % for i
   % Corregimos las presiones si el factor de compresibilidad es de un líquido
   if iz >1
      if is==0
         Pp=Pp-DP;
      else
         DP=DP/2;
         Pp=Pp-DP ;
      end
      
      ype=yc;
   else
      is=1;
      Pc=Pp;
      Tc=Tp;
      yc = ype;
      if abs(PF-Pp)<DP
         Pp=PF;
      else
         Pp=Pc+DP*sign(PF-Pp);
      end    
   end
end


   
      

      
   
 'Comprobación';
Tc;
Pc;
[k,iz]=K(Pc,Tc,z,yc);
yc;
yc=K(Pc,Tc,z,yc).*z;
sum(yc);
Tb=Tc;
