program simani5

      implicit real*8 (a-h, o-z)
      parameter (n = 12, nd = 3000)
c      parameter (n = 13, nd = 3000)
      character*8 fil1, fil2, fil3
      real*8 mu, ms
      dimension x(n), vm(n), v(n)
      common /r1/ d1, d2, d3, d4, d5, d6, d7, d8, qp(nd), hb(nd), 
     + rw(nd), ms(nd), mu(nd), heh(nd), q(nd), he(nd), hh(nd), hu(nd), 
     + b(nd), th(nd), fi(nd), stdpi 
      common /i1/ np, istd
      
      fil1 = 'weirdat1'
      fil2 = 'jn1407a1'
      fil3 = 'jn1407a2'
      
      open (1, file=fil1//'.prn')
      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)
      qp(i) = q(i)/(b(i)*dsqrt(gr*he(i)**3))
      hb(i) = he(i)/b(i)
      heh(i) = hh(i)/he(i)
      rw(i) = he(i)*he(i)/(he(i)*he(i)+b(i)*hh(i))
      mu(i) = 1.d0/dtan(th(i))
      ms(i) = dtan(fi(i))
      if (qp(i) .le. 0.0d0) qp(i) = 1.d-10
      if (hb(i) .le. 0.0d0) hb(i) = 1.d-10
      if (heh(i) .le. 0.0d0) heh(i) = 1.d-10
      if (rw(i) .le. 0.0d0) rw(i) = 1.d-10
      if (ms(i) .le. 0.0d0) ms(i) = 1.d-10
      if (mu(i) .le. 0.0d0) mu(i) = 1.d-10
      i = i + 1
      go to 20
  100 continue
      np = i-1
      istd = 0
      
      call siman (x, rel, vm)

      write (2, *) 'd1, d2, d3, d4, d5, d6, d7, d8, d9'
      write (2, '(6(g14.5,1h,))') (x(i), i=1,n)
      write (2, *) 'v1, v2, v3, v4, v5, v6, v7, v8, v9'
      write (2, '(6(g14.5,1h,))') (vm(i), i=1,n)

      istd = 1
      call fcn (np, x, rel)

      write (3, *) 'Hb, R, Mu, Ms, He, Hh, Hu, B, Theta, Fi, ', 
     +               'Q, Qh, Error, Qpi, Qpih, Rel-Pi-Err'
      
      dummy = 0.d0
      do i = 1, np
         v(1) = hb(i)
         v(2) = heh(i)
         v(3) = rw(i)
         v(4) = ms(i)
         v(5) = mu(i)
         
         call fit (n, x, 5, v, qpih, dummy)
         
         qh = qpih*b(i)*dsqrt(gr*he(i)**3)
         qerr = qh-q(i)
         relpi = dabs(qpih-qp(i))/dsqrt(dabs(qpih*qp(i)))
         if (qpih .le. 0.d0) then
            relpi = 2.d0 + relpi
         endif
      
         write (3, '(16(g14.5,1h,))') 
     +   hb(i), rw(i), mu(i), ms(i), he(i), hh(i), hu(i), b(i), th(i), 
     +   fi(i), q(i), qh, qerr, qp(i), qpih, relpi

      end do

      close(1)
      close(2)
      close(3)
      
      end
      

      subroutine fit (n, d, m, v, q, rel)
      ! fitd0
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      
c      c1 = d(1) + rw*(d(2) + rw*d(3))
c      c1 = d(1) + rw*(d(2) + rw*rw*d(3))
      c1 = d(1) + rw*rw*(d(2) + rw*rw*d(3))
      if (c1 .lt. 0.0d0) then
         rel = rel - 10.d0*c1 
      endif
      c2 = 1.d0+(rw/d(6))**d(7)
      c3 = d(4) +d(5)/c1
      if (c3 .lt. 0.d0) then
         rel = rel - 10.d0*c3
      endif
      
      c4 = 1.d0+(rw/d(10))**d(11)
c      if (dabs(c3) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - c3*c3) + 10.d0
c         return
c      endif
      c5 = d(8) +d(9)/c3
      if (c5 .lt. 0.d0) c5 = 0.d0
      
      q = c1 + c3*hb*ms
     +      + c5*( hb*hb*(ms*ms+mu*mu) )**d(12) 
      
      return
      end            
      

      subroutine fcn (n, x, rel)

      implicit real*8 (a-h, o-z)
      dimension x(n), v(20)
      real*8 mu, ms, merr
      parameter (nd = 3000)
      common /r1/ d1, d2, d3, d4, d5, d6, d7, d8, qp(nd), hb(nd), 
     + rw(nd), ms(nd), mu(nd), heh(nd), q(nd), he(nd), hh(nd), hu(nd), 
     + b(nd), th(nd), fi(nd), stdpi 
      common /i1/ np, istd


      gr = 32.174d0
      mrerr = 0.d0
      rel = 0.d0
      relpi = 0.d0
      stdpi = 0.d0
      rms = 0.d0
      std = 0.d0
      dummy = 0.d0

      do i = 1, np
         v(1) = hb(i)
         v(2) = heh(i)
         v(3) = rw(i)
         v(4) = ms(i)
         v(5) = mu(i)
         
         call fit (n, x, 5, v, qpih, rel)
         
         qh = qpih*b(i)*dsqrt(gr*he(i)**3)

      if (istd .eq. 1) then 
         rms1 = qh-q(i)
         rms = rms + rms1*rms1
         merr = merr + rms1
         relpi1 = dabs(qpih-qp(i))/dsqrt(dabs(qpih*qp(i)))
         relpi = relpi + relpi1
         stdpi = stdpi + relpi1*relpi1
         std = std + rms1
      end if

         rel1 = dabs(qh-q(i))
         if (rel1 .gt. 1.d0) rel1 = rel1*rel1
         rel = rel + rel1
      
      end do
      
      rel = rel + 0.1d0*dabs(x(6)-1.d0) + 0.1d0*dabs(x(7)-1.d0)
      xn = dble(np)
      rel = rel/xn

      if (istd .eq. 1) then 
         merr = merr/xn
         stdpi = dsqrt((xn*stdpi-relpi*relpi)/(xn*(xn-1)))
         std = dsqrt((xn*rms-std*std)/(xn*(xn-1)))
         relpi2 = relpi/xn
         print *, 'np   mean-error   std-error   relative err   ',
     +          'std rel err'
         print *, np, merr, std, relpi2, stdpi 
         write (2, *) 'np   mean-error   std-error   relative err',   
     +          '   std rel err'
         write (2, '(i6,1h,, 4(g14.5,1h,))') np, merr, std, relpi2, 
     +                                       stdpi 
      end if
      
      return
      end


      subroutine siman (xopt, fopt, vm)
      
      implicit real*8 (a-h, o-z)
      
      parameter (n = 12, neps = 4)
c      parameter (n = 13, 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  maxf

      external fcn 

c  set input values of the input/output parameters.
      t = 100.d0
      maxf = .false.
      eps = 1.0d-8
      rt = .9d0
      iseed1 = 787
      iseed2 = 400
      ns = 21
      nt = 5
      maxevl = 800000
      iprint = 1
      do i = 1, n
         c(i) = 1.618d0
      end do

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

      x(1) = .463d0
      x(2) = 0.0001d0
      x(3) = 0.0001d0
      x(4) = .325d0
      x(5) = -.163d0
      x(6) = .514d0
      x(7) = 8.68d0
      x(8) = -.108d0
      x(9) = .73d0
      x(10) = .508d0
      x(11) = 2.52d0
      x(12) = .384d0
      
      lb(1) = .00001d0
      lb(6) = 0.00001d0
      lb(10) = 0.00001d0
      lb(12) = 0.00001d0

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

      write(*,1000) n, maxf, 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,maxf,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,maxf,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(n), lb(n), ub(n), c(n), vm(n), fstar(n),
     1                  xopt(n), xp(n), t, eps, rt, fopt
      integer  nacp(n), n, ns, nt, neps, nacc, maxevl, iprint,
     1         nobds, ier, nfcnev, iseed1, iseed2
      logical  maxf

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.
      iseed1 = 1
      iseed2 = 2
      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. maxf) f = -f
      nfcnev = nfcnev + 1
      fopt = f
      fstar(1) = f
      if(iprint .ge. 1) call prt2(maxf,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(maxf,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. maxf) fp = -fp
               nfcnev = nfcnev + 1
               if(iprint .ge. 3) call prt4(maxf,n,xp,x,fp,f)

c  if too many function evaluations occur, terminate the algorithm.
               if(nfcnev .ge. maxevl) then
                  call prt5
                  if (.not. maxf) 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(maxf)
                     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(maxf)
                  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)
c            if (ratio .gt. .6) then
c               vm(i) = vm(i)*(1. + c(i)*(ratio - .6)/.4)
c            else if (ratio .lt. .4) then
c               vm(i) = vm(i)/(1. + c(i)*((.4 - ratio)/.4))
c            end if
            if (ratio .gt. .618) then
c               vm(i) = vm(i)*(1. + (c(i)-1.)*(ratio - .618)/.382)
               vm(i) = vm(i)*(1. + 1.618*(ratio - .618))
            else if (ratio .lt. .382) then
c               vm(i) = vm(i)/(1. + (c(i)-1.)*(.382 - ratio)/.382)
               vm(i) = vm(i)/(1. + 1.618*(.382 - ratio))
            end if
            if (vm(i) .gt. (ub(i)-lb(i))) then
c               vm(i) = ub(i) - lb(i)
c               vm(i) = (0.75d0+0.2499*ranmar())*(ub(i) - lb(i)) 
               vm(i) = (0.618d0+0.382*ranmar())*(ub(i) - lb(i)) 
c               vm(i) = 0.618d0*(ub(i) - lb(i)) 
c               vm(i) = (0.309d0+0.382*ranmar())*(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(maxf,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. maxf) 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()
      real*8 u(97), c, cd, cm
      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
         if (uni .gt. 1.0) uni = 0.9999
         ranmar = uni
      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(maxf,n,x,f)

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

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

      return
      end

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

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

      write(*,'(''  '')')
      call prtvec(x,n,'current x')
      if (maxf) 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(maxf,n,xp,x,fp,f)

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

      write(*,'(''  '')')
      call prtvec(x,n,'current x')
      if (maxf) 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(maxf)

      logical  maxf

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

      return
      end

      subroutine prt7(maxf)

      logical  maxf

      if (maxf) 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(maxf,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  maxf

      totmov = nup + ndown + nrej

      write(*,'(/,
     1  '' intermediate results before next temperature reduction'',/)')
      write(*,'(''  current temperature:            '',g12.5)') t
      if (maxf) 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
      

      subroutine fitb8 (n, d, m, v, q, rel)
c      n = 8      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 c(20)
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      w = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1) 
      c(1) = d(2)+w**d(3)
      if (c(1) .le. 0.d0) then
         c(1) = 1.d-8
         rel = rel + 100.
      end if
      c(2) = d(4)+w**d(5)
      if (c(2) .le. 0.d0) then
         c(2) = 1.d-8
         rel = rel + 100.
      end if
      q = q + c(1)*ms + c(2)*((1.d0-w)**d(6))*((ms*ms+mu*mu)**d(7))
      
      return
      end            
      

      subroutine fitb9 (n, d, m, v, q, rel)
c      n = 8      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
c      real*8 c(20)
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      w = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1) 
      c1 = w-d(3)
      if (c1 .le. 0.d0) then
         c1 = 1.d-8
         rel = rel + 100.
      else
         c1 = c1**d(4)
      end if
      c2 = d(5)+w**d(6)
      if (c2 .le. 0.d0) then
         c2 = 1.d-8
         rel = rel + 100.
      end if
      q = q + d(2)*c1*ms + c2*((1.d0-w)**d(7))*((ms*ms+mu*mu)**d(8))
      
      return
      end            
      

      subroutine fitc0 (n, d, m, v, q, rel)
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
c      real*8 c(20)
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      w = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1) 
      t1 = 1.-(d(3)-d(4))**40
      t2 = 1.-(d(6)-d(7))**40
      if ((t1 .gt. 0.d0) .or. (t2 .gt. 0.d0)) then
         rel = rel + 1000.
         return
      end if
      c1 = 1.d0 
     +   +( d(3)*exprep(-d(4)*w) - d(4)*exprep(-d(3)*w) )/(d(4)-d(3))
      if (c1 .lt. 0.d0) then
         c1 = 1.d-8 - c1
         rel = rel + 100.
      end if
      c2 =  exprep( -d(6)*(w+d(8)) ) -  exprep( -d(7)*(w+d(8)) )       
      c2 = c2 /(d(7)-d(6))
      if (c2 .lt. 0.d0) then
         c2 = 1.d-8 -c2
         rel = rel + 100.
      end if
      c3 = w-d(8)
      if (c3 .lt. 0.d0) then
         c3 = 1.d-8 - c3
         rel = rel + 100.
      end if      
      q = q + d(2)*c1*hb*ms 
     +    + d(5)*c2*((hb*hb*(ms*ms+mu*mu))**d(9))
      
      return
      end            
      

      subroutine fitc1 (n, d, m, v, q, rel)
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
c      real*8 c(20)
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      r = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1) 
c      t1 = 1.-(d(3)-d(4))**40
c      t2 = 1.-(d(6)-d(7))**40
c      if ((t1 .gt. 0.d0) .or. (t2 .gt. 0.d0)) then
c         rel = rel + 1000.
c         return
c      end if
      c1 = d(2) +r*(d(3) +r*(d(4) +r*d(5))) 
      if (c1 .lt. 0.d0) then
         c1 = 1.d-8 
         rel = rel + 1000.*c1
      end if
      c2 = exprep( -d(7)*(r+d(8)) ) -  exprep( -(1.d0+d(7))*(r+d(8)) )       
      if (c2 .lt. 0.d0) then
         c2 = 1.d-8
         rel = rel + 1000.*c2
      end if
      c3 = r-d(8)
      if (c3 .lt. 0.d0) then
         c3 = 1.d-8
         rel = rel + 1000.*c3
      end if      
      q = q + c1*hb*ms 
     +    + d(6)*c2*((hb*hb*(ms*ms+mu*mu))**d(9))
      
      return
      end            
      

      subroutine fitc2 (n, d, m, v, q, rel)
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      r = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1)
       
      t1 = 100.-( 1.d4*(d(6)-d(4))  )**2
      t2 = 100.-( 1.d4*(d(11)-d(9)) )**2
      if (t1 .gt. 0.d0) then
         rel = rel + 10.d0*t1
         return
      end if
      if (t2 .gt. 0.d0) then
         rel = rel + 10.d0*t2
         return
      end if
      rel = rel + 10.d0*(d(7)-1.d0)**2
      
      c1 = r+d(5)
c      if (c1 .lt. 0.d0) rel = rel - 1000.d0*c1
      c2 = exprep( -d(4)*c1 ) -  exprep( -d(6)*c1 ) 
      c2 = d(2) + d(3)*c2/(d(6)-d(4))
      if (c2 .lt. 0.d0) rel = rel - 1000.*c2
      
      c3 = r + d(10)
c      if (c3 .lt. 0.d0) rel = rel - 1000.d0*c3
      c4 = exprep( -d(9)*c3 ) -  exprep( -d(11)*c3 ) 
      c4 = d(8)*c4/(d(11)-d(9))      
      if (c4 .lt. 0.d0) rel = rel - 1000.*c4
      
      q = q + c2*((hb*ms)**d(7)) + c4*((hb*hb*(ms*ms+mu*mu))**d(12))
      
      return
      end            
      

      subroutine fitc3 (n, d, m, v, q, rel)
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      r = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1)
       
      c1 = d(2) +r*(d(3) +r*(d(4) +r*d(5)))
      if (c1 .lt. 0.d0) rel = rel - 1000.*c1
      c2 = d(7) +r*(d(8) +r*(d(9) +r*d(10)))
      if (c2 .lt. 0.d0) rel = rel - 1000.*c2
      
      q = q + c1*( hb**d(6) )*ms 
     +      + c2*( hb**d(11) )*( (ms*ms+mu*mu)**d(12) )
      
      return
      end            
      

      subroutine fitc4 (n, d, m, v, q, rel)
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      r = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1)
       
      c1 = d(2) +r*(d(3) +r*(d(4) +r*(d(5) +r*d(6))))
      if (c1 .lt. 0.d0) rel = rel - 1.d6*c1
      c2 = d(7) +r*(d(8) +r*(d(9) +r*(d(10) +r*d(11))))
      if (c2 .lt. 0.d0) rel = rel - 1.d6*c2
      
      q = q + c1*hb*ms 
     +      + c2*( hb**d(12) )*( (ms*ms+mu*mu)**d(13) )
      
      return
      end            
      

      subroutine fitc5 (n, d, m, v, q, rel)
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      r = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1)
       
      c1 = d(2) +r*(d(3) +r*r*(d(4) +r*r*(d(5) +r*r*d(6))))
      if (c1 .lt. 0.d0) rel = rel - 1.d6*c1
      c2 = d(7) +r*(d(8) +r*r*(d(9) +r*r*(d(10) +r*r*d(11))))
      if (c2 .lt. 0.d0) rel = rel - 1.d6*c2
      
      q = q + c1*hb*ms 
     +      + c2*( hb**d(12) )*( (ms*ms+mu*mu)**d(13) )
      
      return
      end            
      

      subroutine fitc6 (n, d, m, v, q, rel)
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      r = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1)
       
      c1 = d(2) +r*r*(d(3) +r*r*(d(4) +r*r*(d(5) +r*r*d(6))))
      if (c1 .lt. 0.d0) rel = rel - 1.d6*c1
      c2 = d(7) +r*r*(d(8) +r*r*(d(9) +r*r*(d(10) +r*r*d(11))))
      if (c2 .lt. 0.d0) rel = rel - 1.d6*c2
      
      q = q + c1*hb*ms 
     +      + c2*( hb**d(12) )*( (ms*ms+mu*mu)**d(13) )
      
      return
      end
      
      
      subroutine fitc7 (n, d, m, v, q, rel)
      ! fitc7
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      r = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1)
       
      c1 = d(2) +d(3)/(1.d0+(r/d(4))**d(5))
      if (c1 .lt. 0.d0) rel = rel - 1.d6*c1
      c2 = d(8) +r*r*(d(9) +r*r*(d(10) +r*r*d(11)))
      if (c2 .lt. 0.d0) rel = rel - 1.d6*c2
      
      q = q + c1*(hb**d(6))*(ms**d(7))
     +      + c2*( hb**d(12) )*( (ms*ms+mu*mu)**d(13) )
      
      return
      end            
                  
      

      subroutine fitc8 (n, d, m, v, q, rel)
      ! fitc8
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      r = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1)
      
      c1 = 1.d0+(r/d(4))**d(5)
      if (dabs(c1) .lt. 0.01d0) then
         rel = rel + 1.d8*(1.d-4 - c1*c1) + 10.d0
         return
      endif
      c2 = d(2) +d(3)/c1
      if (c2 .lt. 0.d0) rel = rel - 1.d6*c2
      
      c3 = 1.d0+(r/d(10))**d(11)
      if (dabs(c3) .lt. 0.01d0) then
         rel = rel + 1.d8*(1.d-4 - c3*c3) + 10.d0
         return
      endif
      c4 = d(8) +d(9)/c3
      if (c4 .lt. 0.d0) rel = rel - 1.d6*c4
      
      q = q + c2*(hb**d(6))*(ms**d(7))
     +      + c4*( hb**d(12) )*( (ms*ms+mu*mu)**d(13) )
      
      return
      end            
      

      subroutine fitc9 (n, d, m, v, q, rel)
      ! fitc9
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      r = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1)
      
      c1 = 1.d0+(r/d(4))**d(5)
      if (dabs(c1) .lt. 0.01d0) then
         rel = rel + 1.d8*(1.d-4 - c1*c1) + 10.d0
         return
      endif
      c2 = d(2) +d(3)/c1
      if (c2 .lt. 0.d0) rel = rel - 1.d6*c2
      
      c3 = 1.d0+(r/d(10))**d(11)
      if (dabs(c3) .lt. 0.01d0) then
         rel = rel + 1.d8*(1.d-4 - c3*c3) + 10.d0
         return
      endif
      c4 = d(8) +d(9)/c3
      if (c4 .lt. 0.d0) c4 = 0.d0
      
      q = q + c2*(hb**d(6))*(ms**d(7))
     +      + c4*( hb**d(12) )*( (ms*ms+mu*mu)**d(13) )
      
      return
      end            
      

      subroutine fitd0 (n, d, m, v, q, rel)
      ! fitd0
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      r = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1)
      
      c1 = 1.d0+(r/d(4))**d(5)
c      if (dabs(c1) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - c1*c1) + 10.d0
c         return
c      endif
      c2 = d(2) +d(3)/c1
      if (c2 .lt. 0.d0) then
         rel = rel - 1.d6*c2
         c2 = 0.d0
      endif
      
      c3 = 1.d0+(r/d(8))**d(9)
c      if (dabs(c3) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - c3*c3) + 10.d0
c         return
c      endif
      c4 = d(6) +d(7)/c3
      if (c4 .lt. 0.d0) c4 = 0.d0
      
      q = q + c2*hb*ms
     +      + c4*( hb*hb*(ms*ms+mu*mu) )**d(10) 
      
      return
      end            
      

      subroutine fitd1 (n, d, m, v, q, rel)
      ! fitd1
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      heh = v(2)
      r = v(3)
      ms = v(4)
      mu = v(5)
      q = d(1)
      
      c1 = 1.d0+(r/d(4))**d(5)
      c2 = d(2) +d(3)/c1
      if (c2 .lt. 0.d0) then
         rel = rel - 1.d6*c2
         c2 = 0.d0
      endif
      
      c3 = 1.d0+(r/d(8))**d(9)
      c4 = d(6) +d(7)/c3
      if (c4 .lt. 0.d0) c4 = 0.d0
      
      q = q + c2*hb*ms
     +      + c4*( (hb**2.5d0)*(ms*ms+mu*mu) )**d(10) 
      
      return
      end