c-----------------------------------------------------------------------
c Find the fugacity of nitrogen
c-----------------------------------------------------------------------
c Instructor: Nam Sun Wang
c-----------------------------------------------------------------------
c Programming Note:
c   Call QDAGS of IMSL to integrate in FUNCTION FUGACITY
c-----------------------------------------------------------------------

c Pass all species-dependent constants through the common block.
      common /cblock/A0, a, B0, b, c, R, Temp

c The values for nitrogen are (from Example 3.10 of Riggs) -------------
      A0 = 136.2315   ! in kPa m^6/kmole^2
      a  = 0.02617    ! in m^3/kmole
      B0 = 0.05046    ! in m^3/kmole
      b  = -0.00691   ! in m^3/kmole
      c  = 4.2E4      ! in m^3*K^3/kmole
      R  = 8.3144     ! Gas law constant in kPa*m^3/(kmole-K)=Joule/(mole-K)

      print *, 'Find the fugacity of nitrogen.'
      print *, 'Enter pressure (kPa): '
      read *, P
      print *, 'Enter temperature (Kelvin): '
      read *, T
      print *, 'The fugacity (in kPa) is: ', fugacity(P,T)
      end


c-----------------------------------------------------------------------
      function fugacity(P,T)
c-----------------------------------------------------------------------
c Return the fugacity based on the Beattie-Bridgemen Equation of State
c     fugacity ... fugacity of a gas in kPa (output)
c     P        ... pressure in kPa (1atm = 101.325kPa) (input)
c     T        ... temperature in Kelvin (input)
c-----------------------------------------------------------------------
c Fugacity is calculated from the following equation.
c                   ôP
c  ln(fugacity/P) = ³  (p*V/(R*T)-1)/p dp
c                   õ0
c  ln(fugacity/P) = integral (from 0 to P) (p*V/(R*T)-1)/p dp
c-----------------------------------------------------------------------
c Programming Note:
c   Call QDAGS of IMSL to integrate
c-----------------------------------------------------------------------

      external fugafcn
      common /cblock/A0, a, B0, b, c, R, Temp

c Cannot include the function argument in COMMON; so store it in another variable
      Temp = T

c Numerically evaluate the integral first, then find the fugacity ------
c   call QDAGS(func, a, b, errabs, errrel, s, errest)
c   where func ... name of the function to integrate,
c                  must be in the form func(x)
c         a    ... lower integration limit (input)
c         b    ... upper integration limit (input)
c         errabs.. absolute error tolerance (input)
c         errrel.. relative error tolerance (input)
c         s    ... estimate of the integral (output)
c         errest ... estimated error (output)
      call QDAGS(fugafcn, 0., P, 1.E-4, 0., xintegral, errest)
      fugacity = exp(xintegral)*P
      end


c-----------------------------------------------------------------------
      function fugafcn(P)
c-----------------------------------------------------------------------
c Specify the integrand function in fugacity calculation
c
c Fugacity is calculated from the following equation.
c  ln(fugacity/P) = integral (from 0 to P) (p*V/(R*T)-1)/p dp
c-----------------------------------------------------------------------
      common /cblock/A0, a, B0, b, c, R, T

c Create an output vector whose size matches the input
      fugafcn = ( P*V(P,T) / (R*T) - 1.) / P
      end


c-----------------------------------------------------------------------
      function V(P,T)
c-----------------------------------------------------------------------
c Return the specific volume of an nonideal gas
c calculated from Beattie-Bridgemen Equation of State
c     V ... specific volume in m^3/kmole=liter/mole (output)
c     P ... pressure in kPa (1atm = 101.325kPa) (input)
c     T ... temperature in Kelvin (input)
c All the necessary constants are available in a common block
c-----------------------------------------------------------------------
c Beattie-Bridgeman Equation of State (EOS)
c       R*T*(1-e)          A
c   P = ---------*(V+B) - ----
c         V**2            V**2
c where
c     A = A0*(1-a/V)
c     B = B0*(1-b/V)
c     e = c/(V*T**3)
c Multiply both sides of the EOS by V**4 gives a 4th degree polynomial
c     P*V**4 = R*T*(V-c/T**3)*(V**2+B0*(V-b)) - V*A0*(V-a)
c Collect terms ...
c  0 = P*V**4 - R*T*V**3 + (R*c/T**2-R*T*B0+A0)*V**2 + ...
c       + (R*c*B0/T**2+R*T*B0*b-A0*a)*V - R*c*B0*b/T**2
c-----------------------------------------------------------------------
c Programming Note:
c   Call of ZPORC of IMSL to find the polynomial roots.
c-----------------------------------------------------------------------

c Declare variables ----------------------------------------------------
      parameter (n=4)
      real coeff(n+1)
      complex roots(n)
c All species-dependent constants comes through the common block.
      common /cblock/A0, a, B0, b, c, R, Temp

c Make up the coefficient vector ---------------------------------------
      coeff(5) = P                               ! V**4 term
      coeff(4) = -R*T                            ! V**3 term
      coeff(3) = R*c/T**2 - R*T*B0+A0            ! V**2 term
      coeff(2) = R*c*B0/T**2 + R*T*B0*b - A0*a   ! V term
      coeff(1) = -R*c*B0*b/T**2                  ! constant term

c Call a routine to find the roots -------------------------------------
c     CALL ZRHQR(coeff, n, rootr, rooti)   ! from Numerical Recipes
      CALL ZPORC (n, coeff, roots)

c Return a valid root; pick the root near the ideal gas law. -----------
c The 4th root seems to be the right one
      v = roots(4)

      end