c-----------------------------------------------------------------------
c Regression of Antoine coefficients by minimizing sse.
c Estimate the Antoine constants of water by fitting the given steam data
c to the equation
c   Psat = exp( A-B/(T+C))
c where Psat is the saturation pressure in psia and T is temperature in øF.
c Programming Note: Use POWELL of Numerical Recipes (p411)
c                   Link with NR.LIB
c Instructor: Nam Sun Wang
c-----------------------------------------------------------------------

c Declare arrays -------------------------------------------------------
      parameter (nmax=1000, m=3)
      real param(m), xi(m,m)
      dimension T(nmax), P(nmax)

c Extend data T and P to external subroutines
      common /cblock/ T, P, ipt

c Program header -------------------------------------------------------
      print *, 'Nonlinear regression of Antoine coefficients of water'
      print *, '  Psat = exp(A-B/(T+C))'

c Read data from a file ------------------------------------------------
      open(1,file='steam.prn',status='old')
      do i=1, nmax
19      read(1,*,err=19,end=21) T(i), P(i)
      enddo
21    close(1)
      ipt=i-1
      print *, ' Number of points read = ', ipt

c Initial guess --------------------------------------------------------
      param(1) =    15.
      param(2) = 10000.
      param(3) =   460.

c Call a routine from Numerical Recipes to minimize sse(x) -------------
c Usage: call powell(p, xi, m, mp, ftol, iter, fret)
c   where
c   func   ... a fixed function name (which is not part of the argument)
c              POWELL calls func(p)
c              p ... independent variable
c   p(m)      ... initial starting point (input); best point (output)
c   xi(mp,mp) ... matrix containing the initial set of directions (input)
c   m         ... dimension; number of variables (input)
c   mp        ... physical row dimension of xi as declared
c   ftol      ... fractional tolerance in the function value (input)
c   iter      ... number of iterations taken
c   fret      ... value of f at p

c Set the initial directions to be unit vectors
      do i=1, m
        xi(i,i)=1.
      enddo

      call POWELL(param, xi, m, m, 1.e-4, iter, sse)

c Print the coefficients -----------------------------------------------
      A = param(1)
      B = param(2)
      C = param(3)
      print *, 'The coefficients are: '
      print *, '  A = ', A
      print *, '  B = ', B
      print *, '  C = ', C
      print *, 'The value of sse is: ', sse

      end



c-----------------------------------------------------------------------
      function func(param)
c-----------------------------------------------------------------------
c Evaluate the function (which is sse) to be minimized
c Instructor: Nam Sun Wang
c ----------------------------------------------------------------------

      parameter (nmax=1000)
      real param(1)
      dimension T(nmax), P(nmax)

c Extend data T and P to external subroutines
      common /cblock/ T, P, ipt

c Use my own variables
      A = param(1)
      B = param(2)
      C = param(3)

c Calculate sse (sum of squared error)
      sse = 0.
      do i=1, ipt
        error = P(i) - exp(A-B/(T(i)+C))
        sse   = sse + error*error
      enddo

c Make an assignment to the function name
      func = sse

      end