program simang2

      implicit real*8 (a-h, o-z)
      parameter (n = 8, nd = 2500)
      external flushunit
      integer flushunit      
      character*8 fil1, fil2, fil3
      real*8 mu, ms, mui, msi
      dimension x(n), v(n)
      common /r1/ d1, d2, d3, d4, d5, d6, d7, d8, qp(nd), hb(nd), 
     + w(nd), ms(nd), mu(nd), a2(nd), q(nd), he(nd), hh(nd), hu(nd), 
     + b(nd), th(nd), fi(nd) 
      common /i1/ np
      
      fil1 = 'my16rawb'
      fil2 = 'my1707a1'
      fil3 = 'my1707a2'
      open (1, file=fil1//'.prn')
c      rewind 1
      open (2, file=fil2//'.prn')
      open (3, file=fil3//'.prn')
      read (1, *)
      gr = 32.174d0

   10 continue
      i = 1

   20 read (1, *, end=100) q(i), he(i), hh(i), th(i), hu(i), b(i), fi(i)
      b(i) = b(i)/12.d0
      qp(i) = q(i)/(b(i)*dsqrt(gr*he(i)**3))
      hb(i) = he(i)/b(i)
      w(i) = he(i)/(hh(i)+he(i))
      mu(i) = dcotan(th(i))
      ms(i) = dtan(fi(i))
      if (ms(i) .le. 0.0) ms(i) = 1.d-10
      i = i + 1
      go to 20
  100 continue
      np = i-1
      
      call siman (x, rel, v)

      write (2, *) 'd1, d2, d3, d4, d5, d6, d7, d8'
      write (2, '(10g14.5)') (x(i), i=1,n)
      write (2, *) 'relpi'
      write (2, '(g14.5)') rel
      write (2, *) 'v1, v2, v3, v4, v5, v6, v7, v8'
      write (2, '(10g14.5)') (v(i), i=1,n)

      np0 = 0
      rel = 0.d0
      relpi = 0.d0
      d1 = x(1)
      d2 = x(2)
      d3 = x(3)
      d4 = x(4)
      d5 = x(5)
      d6 = x(6)
      d7 = x(7)
      d8 = x(8)
      write (3, *) 'Hb, W, Mu, Ms, He, Hh, Hu, B, Theta, Fi, 
     +               Q, Qh, Rel-Error, Qpi, Qpih, Rel-Pi-Err'

      do i = 1, np
         
         qpi = qp(i)
         hbi = hb(i)
         wi = w(i)
         wim = 1.d0 - wi
         mui = mu(i)
         msi = ms(i)
         
         if (wi .eq. 0.0) wi = 1.d-10
         if (wim .eq. 0.0) wim = 1.d-10
         tt = (wi**d3)*(msi**d4) + d5*(wim**d6)*(mui**d7)
         if (tt .lt. 0.0) then
            tt = -((-tt)**d8)
         else
            tt = tt**d8
         end if
         qph = d1 + d2*hbi*tt 
         qh = qph*b(i)*dsqrt(gr*he(i)**3)/12.d0
         if (qh .le. 0.d0) then
            qh = 1.d-10
            qph = 1.d-10
            rel = rel + 2.d0
            relpi = relpi + 2.d0
         else
            rel1 = dabs(qh-q(i))/dsqrt(dabs(qh*q(i)))
            rel = rel + rel1
            relpi1 = dabs(qph-qpi)/dsqrt(dabs(qph*qpi))
            relpi = relpi + relpi1
         end if
         
         b(i) = 12.d0*b(i)
         write (3, '(12g14.5, f7.4, 2g14.5, f7.4)') 
     +   hbi, wi, mui, msi, he(i), hh(i), hu(i), b(i), th(i), 
     +   fi(i), q(i), qh, rel1, qpi, qph, relpi1
      end do

      rel = rel/dble(np)
      relpi = relpi/dble(np)
      write (2, *) 'np, rel, relpi'
      write (2, '(i5, 2g14.5)') np, rel, relpi
      
      call flushunit (2)
      call flushunit (3)
      close(1)
      close(2)
      close(3)

      end
      

      subroutine fcn (n, x, rel)

      implicit real*8 (a-h, o-z)
      dimension x(n)
      real*8 mu, ms, mui, msi
      parameter (nd=2500)
      common /r1/ d1, d2, d3, d4, d5, d6, d7, d8, qp(nd), hb(nd), 
     + w(nd), ms(nd), mu(nd), a2(nd), q(nd), he(nd), hh(nd), hu(nd), 
     + b(nd), th(nd), fi(nd) 
      common /i1/ np


      rel = 0.d0
      d1 = x(1)
      d2 = x(2)
      d3 = x(3)
      d4 = x(4)
      d5 = x(5)
      d6 = x(6)
      d7 = x(7)
      d8 = x(8)

      do i = 1, np
         
         qpi = qp(i)
         hbi = hb(i)
         wi = w(i)
         wim = 1.d0-wi
         mui = mu(i)
         msi = ms(i)
         
         if (wi .eq. 0.0) wi = 1.d-10
         if (wim .le. 0.0) wim = 1.d-10
         tt = (wi**d3)*(msi**d4) + d5*(wim**d6)*(mui**d7)
         if (tt .lt. 0.0) then
            tt = -((-tt)**d8)
         else
            tt = tt**d8
         end if
         qph = d1 + d2*hbi*tt 
         if (qph .le. 0.d0) then
            qph = 1.d-10
            rel = rel + 2.d0
         else
            rel = rel + dabs(qph-qpi)/dsqrt(dabs(qph*qpi))
         end if
      end do
      rel = rel/dble(np)
      
      return
      end


      subroutine siman (xopt, fopt, vm)
      
      implicit real*8 (a-h, o-z)
c      integer*2 iyy, imm, idd, HRS, MINS, SECS, HSECS
      parameter (n = 8, neps = 4)
      double precision  lb(n), ub(n), x(n), xopt(n), c(n), vm(n),
     1                  fstar(neps), xp(n), t, eps, rt, fopt

      integer  nacp(n), ns, nt, nfcnev, ier, iseed1, iseed2,
     1         maxevl, iprint, nacc, nobds

      logical  max

      external fcn 

c  set input values of the input/output parameters.
      t = 100.d0
      max = .false.
      eps = 1.0d-8
      rt = .9d0
      iseed1 = 1
      iseed2 = 2
      ns = 20
      nt = 5
      maxevl = 800 000
      iprint = 1
      do i = 1, n
         c(i) = 2.0d0
      end do

c  start values and limits
      do i = 1, n
         x(i) = 1.d0
         lb(i) = -5.d0
         ub(i) = 5.d0
      end do
      lb(1)=0.01d0

      
      do i = 1, n
c         vm(i) = ub(i) - lb(i)
         vm(i) = 0.75d0*(ub(i) - lb(i))
      end do

      write(*,1000) n, max, t, rt, eps, ns, nt, neps, maxevl, iprint,
     1              iseed1, iseed2

c      call prtvec(x,n,'starting values')
c      call prtvec(vm,n,'initial step length')
c      call prtvec(lb,n,'lower bound')
c      call prtvec(ub,n,'upper bound')
c      call prtvec(c,n,'c vector')
c      write(*,'(/,''  ****   end of driver routine output   ****''
c     1          /,''  ****   before call to sa.             ****'')')      

      call sa(n,x,max,rt,eps,ns,nt,neps,maxevl,lb,ub,c,iprint,iseed1,
     1        iseed2,t,vm,xopt,fopt,nacc,nfcnev,nobds,ier,
     2        fstar,xp,nacp)

      write(*,'(/,''  ****   results after sa   ****   '')')      
      call prtvec(xopt,n,'solution')
      call prtvec(vm,n,'final step length')
      write(*,1001) fopt, nfcnev, nacc, nobds, t, ier

1000  format(/,' simulated annealing example',/,
     1       /,' number of parameters: ',i3,'   maximazation: ',l5,
     2       /,' initial temp: ', g8.2, '   rt: ',g8.2, '   eps: ',g8.2,
     3       /,' ns: ',i3, '   nt: ',i2, '   neps: ',i2,
     4       /,' maxevl: ',i10, '   iprint: ',i1, '   iseed1: ',i4,
     5       '   iseed2: ',i4)
1001  format(/,' optimal function value: ',g20.13
     1       /,' number of function evaluations:     ',i10,
     2       /,' number of accepted evaluations:     ',i10,
     3       /,' number of out of bound evaluations: ',i10,
     4       /,' final temp: ', g20.13,'  ier: ', i3)

c      stop
      return
      end


      subroutine sa(n,x,max,rt,eps,ns,nt,neps,maxevl,lb,ub,c,iprint,
     1              iseed1,iseed2,t,vm,xopt,fopt,nacc,nfcnev,nobds,ier,
     2              fstar,xp,nacp)

      implicit real*8 (a-h, o-z)
      double precision  x(*), lb(*), ub(*), c(*), vm(*), fstar(*),
     1                  xopt(*), xp(*), t, eps, rt, fopt
      integer  nacp(*), n, ns, nt, neps, nacc, maxevl, iprint,
     1         nobds, ier, nfcnev, iseed1, iseed2
      logical  max

c  type all internal variables.
      double precision  f, fp, p, pp, ratio
      integer  nup, ndown, nrej, nnew, lnobds, h, i, j, m
      logical  quit

c  type all functions.
      double precision  exprep
      real  ranmar 

c  initialize the random number generator ranmar.
      call rmarin(iseed1,iseed2)

c  set initial values.
      nacc = 0
      nobds = 0
      nfcnev = 0
      ier = 99

      do 10, i = 1, n
         xopt(i) = x(i)
         nacp(i) = 0
10    continue

      do 20, i = 1, neps
         fstar(i) = 1.0d+20
20    continue 

c  if the initial temperature is not positive, notify the user and 
c  return to the calling routine.  
      if (t .le. 0.0) then
         write(*,'(/,''  the initial temperature is not positive. ''
     1             /,''  reset the variable t. ''/)')
         ier = 3
         return
      end if

c  if the initial value is out of bounds, notify the user and return
c  to the calling routine.
      do 30, i = 1, n
         if ((x(i) .gt. ub(i)) .or. (x(i) .lt. lb(i))) then
            call prt1
            ier = 2
            return
         end if
30    continue

c  evaluate the function with input x and return value as f.
      call fcn(n,x,f)

c  if the function is to be minimized, switch the sign of the function.
c  note that all intermediate and final output switches the sign back
c  to eliminate any possible confusion for the user.
      if(.not. max) f = -f
      nfcnev = nfcnev + 1
      fopt = f
      fstar(1) = f
      if(iprint .ge. 1) call prt2(max,n,x,f)

c  start the main loop. note that it terminates if (i) the algorithm
c  succesfully optimizes the function or (ii) there are too many
c  function evaluations (more than maxevl).
100   nup = 0
      nrej = 0
      nnew = 0
      ndown = 0
      lnobds = 0

      do 400, m = 1, nt
         do 300, j = 1, ns
            do 200, h = 1, n

c  generate xp, the trial value of x. note use of vm to choose xp.
               do 110, i = 1, n
                  if (i .eq. h) then
                     xp(i) = x(i) + (ranmar()*2.- 1.) * vm(i)
                  else
                     xp(i) = x(i)
                  end if

c  if xp is out of bounds, select a point in bounds for the trial.
                  if((xp(i) .lt. lb(i)) .or. (xp(i) .gt. ub(i))) then
                    xp(i) = lb(i) + (ub(i) - lb(i))*ranmar()
                    lnobds = lnobds + 1
                    nobds = nobds + 1
                    if(iprint .ge. 3) call prt3(max,n,xp,x,fp,f)
                  end if
110            continue

c  evaluate the function with the trial point xp and return as fp.
               call fcn(n,xp,fp)
               if(.not. max) fp = -fp
               nfcnev = nfcnev + 1
               if(iprint .ge. 3) call prt4(max,n,xp,x,fp,f)

c  if too many function evaluations occur, terminate the algorithm.
               if(nfcnev .ge. maxevl) then
                  call prt5
                  if (.not. max) fopt = -fopt
                  ier = 1
                  return
               end if

c  accept the new point if the function value increases.
               if(fp .ge. f) then
                  if(iprint .ge. 3) then
                     write(*,'(''  point accepted'')')
                  end if
                  do 120, i = 1, n
                     x(i) = xp(i)
120               continue
                  f = fp
                  nacc = nacc + 1
                  nacp(h) = nacp(h) + 1
                  nup = nup + 1

c  if greater than any other point, record as new optimum.
                  if (fp .gt. fopt) then
                  if(iprint .ge. 3) then
                     write(*,'(''  new optimum'')')
                  end if
                     do 130, i = 1, n
                        xopt(i) = xp(i)
130                  continue
                     fopt = fp
                     nnew = nnew + 1
                  end if

c  if the point is lower, use the metropolis criteria to decide on
c  acceptance or rejection.
               else
                  p = exprep((fp - f)/t)
                  pp = ranmar()
                  if (pp .lt. p) then
                     if(iprint .ge. 3) call prt6(max)
                     do 140, i = 1, n
                        x(i) = xp(i)
140                  continue
                     f = fp
                     nacc = nacc + 1
                     nacp(h) = nacp(h) + 1
                     ndown = ndown + 1
                  else
                     nrej = nrej + 1
                     if(iprint .ge. 3) call prt7(max)
                  end if
               end if

200         continue
300      continue

c  adjust vm so that approximately half of all evaluations are accepted.
         do 310, i = 1, n
            ratio = dfloat(nacp(i)) /dfloat(ns)
            if (ratio .gt. .6) then
               vm(i) = vm(i)*(1. + c(i)*(ratio - .6)/.4)
            else if (ratio .lt. .4) then
               vm(i) = vm(i)/(1. + c(i)*((.4 - ratio)/.4))
            end if
            if (vm(i) .gt. (ub(i)-lb(i))) then
               vm(i) = ub(i) - lb(i)
            end if
310      continue

         if(iprint .ge. 2) then
            call prt8(n,vm,xopt,x)
         end if

         do 320, i = 1, n
            nacp(i) = 0
320      continue

400   continue

      if(iprint .ge. 1) then
         call prt9(max,n,t,xopt,vm,fopt,nup,ndown,nrej,lnobds,nnew)
      end if

c  check termination criteria.
      quit = .false.
      fstar(1) = f
      if ((fopt - fstar(1)) .le. eps) quit = .true.
      do 410, i = 1, neps
         if (abs(f - fstar(i)) .gt. eps) quit = .false.
410   continue

c  terminate sa if appropriate.
      if (quit) then
         do 420, i = 1, n
            x(i) = xopt(i)
420      continue
         ier = 0
         if (.not. max) fopt = -fopt
         if(iprint .ge. 1) call prt10
         return
      end if

c  if termination criteria is not met, prepare for another loop.
      t = rt*t
      do 430, i = neps, 2, -1
         fstar(i) = fstar(i-1)
430   continue
      f = fopt
      do 440, i = 1, n
         x(i) = xopt(i)
440   continue

c  loop again.
      go to 100

      end

      function  exprep(rdum)
c  this function replaces exp to avoid under- and overflows and is
c  designed for Pentium 4 machines running Watcom F77 v11. it may be
c  necessary to modify it for other machines. note that the maximum
c  and minimum values of exprep are such that they have no effect on
c  the algorithm.

      double precision  rdum, exprep

      if (rdum .gt. 705.) then
         exprep = 1.50525 d+306
      else if (rdum .lt. -705.) then
         exprep = 6.6434 d-307
      else
         exprep = exp(rdum)
      end if

      return
      end

      subroutine rmarin(ij,kl)
c  this subroutine and the next function generate random numbers. see
c  the comments for sa for more information. the only changes from the
c  orginal code is that (1) the test to make sure that rmarin runs first
c  was taken out since sa assures that this is done (this test didn't
c  compile under ibm's vs fortran) and (2) typing ivec as integer was
c  taken out since ivec isn't used. with these exceptions, all following
c  lines are original.

c this is the initialization routine for the random number generator
c     ranmar()
c note: the seed variables can have values between:    0 <= ij <= 31328
c                                                      0 <= kl <= 30081
      implicit real*8 (a-h, o-z)
      real*8 u(97), c, cd, cm
      integer i97, j97
      common /raset1/ u, c, cd, cm, i97, j97
      if( ij .lt. 0  .or.  ij .gt. 31328  .or.
     *    kl .lt. 0  .or.  kl .gt. 30081 ) then
          print '(a)', ' the first random number seed must have a value
     *between 0 and 31328'
          print '(a)',' the second seed must have a value between 0 and
     *30081'
            stop
      endif
      i = mod(ij/177, 177) + 2
      j = mod(ij    , 177) + 2
      k = mod(kl/169, 178) + 1
      l = mod(kl,     169)
      do 2 ii = 1, 97
         s = 0.0
         t = 0.5
         do 3 jj = 1, 24
            m = mod(mod(i*j, 179)*k, 179)
            i = j
            j = k
            k = m
            l = mod(53*l+1, 169)
            if (mod(l*m, 64) .ge. 32) then
               s = s + t
            endif
            t = 0.5 * t
3        continue
         u(ii) = s
2     continue
      c = 362436.0 / 16777216.0
      cd = 7654321.0 / 16777216.0
      cm = 16777213.0 /16777216.0
      i97 = 97
      j97 = 33
      return
      end

      function ranmar()
c      real urand
c      integer iseed
      real*8 u(97), c, cd, cm
      integer i97, j97
      common /raset1/ u, c, cd, cm, i97, j97
         uni = u(i97) - u(j97)
         if( uni .lt. 0.0 ) uni = uni + 1.0
         u(i97) = uni
         i97 = i97 - 1
         if(i97 .eq. 0) i97 = 97
         j97 = j97 - 1
         if(j97 .eq. 0) j97 = 97
         c = c - cd
         if( c .lt. 0.0 ) c = c + cm
         uni = uni - c
         if( uni .lt. 0.0 ) uni = uni + 1.0
c         if ((uni .lt. 0.0) .or. (uni .gt. 1.0)) pause
         if (uni .gt. 1.) uni = 0.99999
         ranmar = uni
c         ranmar = urand(iseed)
      return
      end

      subroutine prt1
c  this subroutine prints intermediate output, as does prt2 through
c  prt10. note that if sa is minimizing the function, the sign of the
c  function value and the directions (up/down) are reversed in all
c  output to correspond with the actual function optimization. this
c  correction is because sa was written to maximize functions and
c  it minimizes by maximizing the negative a function.

      write(*,'(/,''  the starting value (x) is outside the bounds ''
     1          /,''  (lb and ub). execution terminated without any''
     2          /,''  optimization. respecify x, ub or lb so that  ''
     3          /,''  lb(i) .lt. x(i) .lt. ub(i), i = 1, n. ''/)')

      return
      end

      subroutine prt2(max,n,x,f)

      implicit real*8 (a-h, o-z)
      double precision  x(*), f
      integer  n
      logical  max

      write(*,'(''  '')')
      call prtvec(x,n,'initial x')
      if (max) then
         write(*,'(''  initial f: '',/, g25.18)') f
      else
         write(*,'(''  initial f: '',/, g25.18)') -f
      end if

      return
      end

      subroutine prt3(max,n,xp,x,fp,f)

      implicit real*8 (a-h, o-z)
      double precision  xp(*), x(*), fp, f
      integer  n
      logical  max

      write(*,'(''  '')')
      call prtvec(x,n,'current x')
      if (max) then
         write(*,'(''  current f: '',g25.18)') f
      else
         write(*,'(''  current f: '',g25.18)') -f
      end if
      call prtvec(xp,n,'trial x')
      write(*,'(''  point rejected since out of bounds'')')

      return
      end

      subroutine prt4(max,n,xp,x,fp,f)

      implicit real*8 (a-h, o-z)
      double precision  xp(*), x(*), fp, f
      integer  n
      logical  max

      write(*,'(''  '')')
      call prtvec(x,n,'current x')
      if (max) then
         write(*,'(''  current f: '',g25.18)') f
         call prtvec(xp,n,'trial x')
         write(*,'(''  resulting f: '',g25.18)') fp
      else
         write(*,'(''  current f: '',g25.18)') -f
         call prtvec(xp,n,'trial x')
         write(*,'(''  resulting f: '',g25.18)') -fp
      end if

      return
      end

      subroutine prt5

      write(*,'(/,''  too many function evaluations; consider ''
     1          /,''  increasing maxevl or eps, or decreasing ''
     2          /,''  nt or rt. these results are likely to be ''
     3          /,''  poor.'',/)')

      return
      end

      subroutine prt6(max)

      logical  max

      if (max) then
         write(*,'(''  though lower, point accepted'')')
      else
         write(*,'(''  though higher, point accepted'')')
      end if

      return
      end

      subroutine prt7(max)

      logical  max

      if (max) then
         write(*,'(''  lower point rejected'')')
      else
         write(*,'(''  higher point rejected'')')
      end if

      return
      end

      subroutine prt8(n,vm,xopt,x)

      implicit real*8 (a-h, o-z)
      double precision  vm(*), xopt(*), x(*)
      integer  n

      write(*,'(/,
     1  '' intermediate results after step length adjustment'',/)')
      call prtvec(vm,n,'new step length (vm)')
      call prtvec(xopt,n,'current optimal x')
      call prtvec(x,n,'current x')
      write(*,'('' '')')

      return
      end

      subroutine prt9(max,n,t,xopt,vm,fopt,nup,ndown,nrej,lnobds,nnew)

      implicit real*8 (a-h, o-z)
      double precision  xopt(*), vm(*), t, fopt
      integer  n, nup, ndown, nrej, lnobds, nnew, totmov
      logical  max

      totmov = nup + ndown + nrej

      write(*,'(/,
     1  '' intermediate results before next temperature reduction'',/)')
      write(*,'(''  current temperature:            '',g12.5)') t
      if (max) then
         write(*,'(''  max function value so far:  '',g25.18)') fopt
         write(*,'(''  total moves:                '',i8)') totmov
         write(*,'(''     uphill:                  '',i8)') nup
         write(*,'(''     accepted downhill:       '',i8)') ndown
         write(*,'(''     rejected downhill:       '',i8)') nrej
         write(*,'(''  out of bounds trials:       '',i8)') lnobds
         write(*,'(''  new maxima this temperature:'',i8)') nnew
      else
         write(*,'(''  min function value so far:  '',g25.18)') -fopt
         write(*,'(''  total moves:                '',i8)') totmov
         write(*,'(''     downhill:                '',i8)')  nup
         write(*,'(''     accepted uphill:         '',i8)')  ndown
         write(*,'(''     rejected uphill:         '',i8)')  nrej
         write(*,'(''  trials out of bounds:       '',i8)')  lnobds
         write(*,'(''  new minima this temperature:'',i8)')  nnew
      end if
      call prtvec(xopt,n,'current optimal x')
      call prtvec(vm,n,'step length (vm)')
      write(*,'('' '')')

      return
      end

      subroutine prt10

      write(*,'(/,''  sa achieved termination criteria. ier = 0. '',/)')

      return
      end

      subroutine prtvec(vector,ncols,name)
c  this subroutine prints the double precision vector named vector.
c  elements 1 thru ncols will be printed. name is a character variable
c  that describes vector. note that if name is given in the call to
c  prtvec, it must be enclosed in quotes. if there are more than 10
c  elements in vector, 10 elements will be printed on each line.

      integer ncols
      double precision vector(ncols)
      character *(*) name

      write(*,1001) name

      if (ncols .gt. 10) then
         lines = int(ncols/10.)

         do 100, i = 1, lines
            ll = 10*(i - 1)
            write(*,1000) (vector(j),j = 1+ll, 10+ll)
  100    continue

         write(*,1000) (vector(j),j = 11+ll, ncols)
      else
         write(*,1000) (vector(j),j = 1, ncols)
      end if

 1000 format( 10(g12.5,1x))
 1001 format(/,25x,a)

      return
      end