c-----------------------------------------------------------------------
c Flash Distillation
c Given the feed composition and volatility, find the liquid and vapor fraction
c and composition.
c Variables used ...
c   n    ... number of components
c   V    ... vapor phase fraction
c   x(i) ... mole fraction of component i in the liquid phase
c   y(i) ... mole fraction of component i in the vapor phase
c   z(i) ... mole fraction of component i in the feed
c   xtotal... total mole fraction in the liquid phase (should be 1)
c   ytotal... total mole fraction in the vapor phase (should be 1)
c   ztotal... total mole fraction in the feed phase (should be 1)
c   xK(i)... relative volatility of component i
c Use RTBIS of Numerical Recipes to find roots to f(x)=0
c Instructor: Nam Sun Wang
c-----------------------------------------------------------------------
      parameter(nsize=10)
      real x(nsize), y(nsize)
      external f
      common /cblock/n, z(10), xK(10)

c Program header -------------------------------------------------------
      write(*,*)'*********** FLASH DISTILLATION ************'
      write(*,*)
      write(*,*)'                       +----------> V, y(i)'
      write(*,*)'                       |                   '
      write(*,*)'                 +-----+------+            '
      write(*,*)'                 |            |            '
      write(*,*)'                 |   Flash    |            '
      write(*,*)'     z(i) ----->-+            |            '
      write(*,*)'                 |Distillation|            '
      write(*,*)'                 |            |            '
      write(*,*)'                 +-----+------+            '
      write(*,*)'                       |                   '
      write(*,*)'                       +----------> L, x(i)'
      write(*,*)

c Input the feed composition -------------------------------------------
1     write(*,600)' Number of components in the feed: '
      read(*,*)n
      write(*,*)'Input the feed composition in mole fractions ...'

      z(n)=1.
      do 10 i=1, n-1
        write(*,601)' z(', i, '): '
2       read(*,*)z(i)
c       Only non-negative quantities make physical sense
        if(z(i) .lt. 0.)then
          write(*,*)'Negative number!  Re-enter the last number.'
          goto 2
        endif
c       Calculate the balance for the last component
        z(n)=z(n)-z(i)
        write(*,601)' K(', i, '): '
        read(*,*)xK(i)
c       Read the last value
        if(i .eq. n-1)then
          write(*,602)' z(', n, '): ', z(n)
          write(*,601)' K(', n, '): '
          read(*,*)xK(n)
        endif
10    continue
      if(z(n) .lt. 0.)then
        write(*,*)'Invalid composition!  Start over again.'
        goto 1
      endif

c Call a routine to solve the flash distillation equation
c     V = RTBIS(function, lo_guess, hi_guess, accuracy)
      V = RTBIS(f, 0., 1., 1.e-6)

c Output header --------------------------------------------------------
      write(*,*)
      write(*,*)'Vapor fraction (V):  ', V
      write(*,*)'Liquid fraction (L): ', 1.-V
      write(*,*)
      write(*,*)'------------------------------------------------------'
      write(*,*)'                    Composition                       '
      write(*,*)'------------------------------------------------------'
      write(*,*)'Component      Feed           Vapor          Liquid   '
      write(*,*)'------------------------------------------------------'

c Calculate the liquid phase composition & print -----------------------
      do 20 i=1, n
        x(i) = z(i) / (1.+(xK(i)-1.)*V)
        y(i) = xk(i)*x(i)
        xtotal = xtotal + x(i)
        ytotal = ytotal + y(i)
        ztotal = ztotal + z(i)
        write(*,650)i, z(i), x(i), y(i)
20    continue

c Table footer ---------------------------------------------------------
      write(*,*)'------------------------------------------------------'
      write(*,651)' Total    ', ztotal, xtotal, ytotal

600   format(a\)
601   format(a, i1, a\)
602   format(a, i1, a, f9.7)
650   format(i6,4x,1p,3e15.7)
651   format(a,1p,3e15.7)
      end



c-----------------------------------------------------------------------
      function f(V)
c-----------------------------------------------------------------------
c Specify the function to be satisfied for flash distillation
c-----------------------------------------------------------------------
      common /cblock/n, z(10), xK(10)

      f=0.
      do 10 i=1, n
        f = f + (xK(i)-1.)*z(i) / (1.+(xK(i)-1.)*V)
10    continue
      end