program simanj8

      implicit real*8 (a-h, o-z)
c      parameter (n = 11, nd = 3000)
      parameter (n = 22, nd = 3000)
      character*8 fil1, fil2, fil3
      real*8 mu, ms
      dimension x(n), vm(n), v(20)
      common /r1/ 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 = 'weirdat6'
      fil2 = 'jl1707b1'
      fil3 = 'jl1707b2'
      
      open (1, file=fil1//'.prn')
      open (2, file=fil2//'.prn')
      open (3, file=fil3//'.prn')
      read (1, *)
      gr = 32.174d0
      rgr = dsqrt(gr)

   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))
c      hb(i) = he(i)/b(i)
c      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))
      if (mu(i) .lt. 1.d-4) mu(i) = 0.d0
      ms(i) = dtan(fi(i))
      if (ms(i) .lt. 1.d-4) ms(i) = 0.d0
      if (qp(i) .le. 0.0d0) qp(i) = 1.d-10
c      if (hb(i) .le. 0.0d0) hb(i) = 1.d-10
c      if (heh(i) .le. 0.0d0) heh(i) = 1.d-10
      if (rw(i) .le. 0.0d0) rw(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'
c      write (2, '(6(g14.5,1h,))') (x(i), i=1,n)
      write (2, '(6(g17.8,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)

c      write (3, *) 'Hb, R, W1, W2, Sc1, Sc2, Sc3, Mu, Ms, He, Hh, ', 
      write (3, *) 'Mu, Ms, He, Hh, Hu, B, R, W1, W2, ', 
c     +           'Qb, Q1, Q2, Q, Qh, Error, AbsErr, Qpi, Qpih, ',
     +           'Qr, Q1, Q2, Q, Qh, Error, AbsErr, Qpi, Qpih, ',
     +           'Rel-Pi-Err, Theta, Fi'
      
      dummy = 0.d0
      
      do i = 1, np
         v(1) = he(i)
         v(2) = b(i)
         v(3) = rw(i)
         v(4) = ms(i)
         v(5) = mu(i)
         v(6) = hh(i)
         v(7) = hu(i)
         
         call fit (n, x, 12, v, qh, dummy)         
         qpih = qh/(b(i)*dsqrt(gr*he(i)**3))
         qerr = qh-q(i)
         abserr = dabs(qerr)
         relpi = dabs(qpih-qp(i))/dsqrt(dabs(qpih*qp(i)))
         if (qpih .le. 0.d0) then
            relpi = 2.d0 + relpi
         endif
         qb = v(8)
         w1 = v(9)
         w2 = v(10)
         q1 = v(11)
         q2 = v(12)
         
               
c         write (3, '(21(g14.5,1h,))') 
c     +   hb(i), rw(i), v(5), v(6), v(7), v(8), v(9), mu(i), ms(i),  
         write (3, '(21(g14.5,1h,))') 
     +   mu(i), ms(i), he(i), hh(i), hu(i), b(i), rw(i), w1, w2,  
     +   qb, q1, q2, q(i), qh, qerr, abserr,
     +   qp(i), qpih, relpi, th(i), fi(i)

      end do

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

      subroutine fit (n, d, m, v, q, rel)
      ! fith2
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      hh = v(6)
      hu = v(7)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+dexp((rw-d(8))/d(9))
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*dsqrt(sg)

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      q = qb + q1*he + q2*he**d(10)
      
      fe = d(11) +d(12)*dexp(-ms) +d(13)*dexp(-mu) +d(14)*dexp(-hh)
     +          +d(15)*dexp(-hu) +d(16)*dexp(-rw)
      if (fe .le. 0.d0) then
      	rel = rel - 1.d6*fe
      	fe = 0.d0
      endif
      
      fc = 1.d0 +d(17)*(hh**d(18)) +d(19)*(rw**d(20)) 
     +      +d(21)*(sg**d(22))
      if (fc .le. 0.d0) then
      	rel = rel - 1.d6*fc
      	fe = 0.d0
      endif
      
      q = q*5.672213d0*(he**fe)*fc

c      v(8) = qb
      v(9) = w1
      v(10) = w2
c      v(11) = q1
c      v(12) = q2
      v(11) = q1*he
      v(12) = q2*he
      v(8) = qb + v(11)
      
      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/ 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

      if (x(7)*x(9) .lt. 0.d0) rel = -100.*x(7)*x(9)
      
      do i = 1, np
         v(1) = he(i)
         v(2) = b(i)
         v(3) = rw(i)
         v(4) = ms(i)
         v(5) = mu(i)
         v(6) = hh(i)
         v(7) = hu(i)
         
         call fit (n, x, 12, v, qh, rel)
         
         qpih = qh/(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
      
      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)
      
c      parameter (n = 11, neps = 4)
      parameter (n = 22, 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 = 551
      iseed2 = 100
      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) = 0.d0
         ub(i) = 1.d0
      end do

      lb(1) = 1.d-8
       x(1) = .553d0
      ub(1) = 10.d0
      
       x(2) = .143d0
      ub(2) = 1.d0
      
c      lb(3) = -5.d0
       x(3) = 2.39d0
      ub(3) = 10.d0
      
      lb(4) = -1.d0
       x(4) = -.365d0
      ub(4) = 1.d0
      
      lb(5) = .0001d0
       x(5) = .181d0
      ub(5) = 1.d0
      
      lb(6) = -1.d0
       x(6) = .00297d0
      ub(6) = 1.d0
      
c      lb(7) = -5.d0
       x(7) = .899d0
      ub(7) = 2.d0
      
      lb(8) = -.5d0
       x(8) = .261d0
      ub(8) = 1.5d0
      
c      lb(9) = .5d0
       x(9) = .184d0
      ub(9) = 2.d0
      
      x(10) = 1.25d0
      ub(10) = 2.d0
      
      do i = 11, n
         lb(i) = -2.d0
         x(i) = 0.d0
         ub(i) = 2.d0
      end do
      
      x(11) = 1.49d0
      x(12) = -.00836
      x(13) = .0749
      x(14) = -.0733
      x(15) = .0584
      x(16) = -0.0391
      x(17) = .00256
      x(18) = 9.95
      ub(18) = 20.
      x(19) = -.186
      x(20) = 3.71
      ub(20) = 20.
      x(21) = -.00161
      x(22) = .933
      ub(22) = 5.
      
      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 = 499
      iseed2 = 293
      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 fitg0 (n, d, m, v, q, rel)
      ! fitg0
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w1 = 1.d0+(rw/d(4))**d(5)
c      if (dabs(w1) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w1*w1) + 10.d0
c         return
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      
      w2 = 1.d0+(rw/d(8))**d(9)
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*(sg**d(10))
      
      q = qb + q1*he + q2*(he**d(11))*(b**(1.d0-d(11)))
      q = q*5.672213d0*he*dsqrt(he)

      v(6) = qb
      v(7) = w1
      v(8) = w2
      v(9) = q1
      v(10) = q2
      
      return
      end            
      

      subroutine fitg1 (n, d, m, v, q, rel)
      ! fitg1
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+(rw/d(8))**d(9)
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*(sg**d(10))

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      
      q = qb + q1*he + q2*(he**d(11))
      q = q*5.672213d0*he*dsqrt(he)

      v(6) = qb
      v(7) = w1
      v(8) = w2
      v(9) = q1
      v(10) = q2
      
      return
      end            
      

      subroutine fitg2 (n, d, m, v, q, rel)
      ! fitg2
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+dexp((rw-d(8))/d(9))
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*(sg**d(10))

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      
      q = qb + q1*he + q2*(he**d(11))
      q = q*5.672213d0*he*dsqrt(he)

      v(6) = qb
      v(7) = w1
      v(8) = w2
      v(9) = q1
      v(10) = q2
      
      return
      end            
      

      subroutine fitg3 (n, d, m, v, q, rel)
      ! fitg3
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      hh = v(6)
      hu = v(7)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+dexp((rw-d(8))/d(9))
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*(sg**d(10))

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      
      f = d(12) +d(13)*dexp(-ms) +d(14)*dexp(-mu) +d(15)*dexp(-hh)
     +          +d(16)*dexp(-hu)
      
      q = qb + q1*he + q2*(he**d(11))
      q = q*5.672213d0*(he**f)

      v(8) = qb
      v(9) = w1
      v(10) = w2
      v(11) = q1
      v(12) = q2
      
      return
      end            
      

      subroutine fitg4 (n, d, m, v, q, rel)
      ! fitg4
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      hh = v(6)
      hu = v(7)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+dexp((rw-d(8))/d(9))
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*(sg**d(10))

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      
      f = d(12) +d(13)*dexp(-ms) +d(14)*dexp(-mu) +d(15)*dexp(-hh)
     +          +d(16)*dexp(-hu)
      
      q = qb + q1*he + q2*(he**d(11))
      q = q*5.672213d0*(he**f)*(1.d0+d(17)*hh)

      v(8) = qb
      v(9) = w1
      v(10) = w2
      v(11) = q1
      v(12) = q2
      
      return
      end            
      

      subroutine fitg5 (n, d, m, v, q, rel)
      ! fitg5
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      hh = v(6)
      hu = v(7)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+dexp((rw-d(8))/d(9))
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*(sg**d(10))

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      
      f = d(12) +d(13)*dexp(-ms) +d(14)*dexp(-mu) +d(15)*dexp(-hh)
     +          +d(16)*dexp(-hu)
      
      q = qb + q1*he + q2*(he**d(11))
      q = q*5.672213d0*(he**f)*(1.d0+d(17)*hh+d(18)*rw)

      v(8) = qb
      v(9) = w1
      v(10) = w2
      v(11) = q1
      v(12) = q2
      
      return
      end            
      

      subroutine fitg6 (n, d, m, v, q, rel)
      ! fitg6
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      hh = v(6)
      hu = v(7)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+dexp((rw-d(8))/d(9))
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*(sg**d(10))

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      
      f = d(12) +d(13)*dexp(-ms) +d(14)*dexp(-mu) +d(15)*dexp(-hh)
     +          +d(16)*dexp(-hu)
      
      q = qb + q1*he + q2*(he**d(11))
      q = q*5.672213d0*(he**f)*(1.d0+d(17)*hh+d(18)*rw+d(19)*ms)

      v(8) = qb
      v(9) = w1
      v(10) = w2
      v(11) = q1
      v(12) = q2
      
      return
      end            
      

      subroutine fitg7 (n, d, m, v, q, rel)
      ! fitg7
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      hh = v(6)
      hu = v(7)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+dexp((rw-d(8))/d(9))
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*(sg**d(10))

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      
      f = d(12) +d(13)*dexp(-ms) +d(14)*dexp(-mu) +d(15)*dexp(-hh)
     +          +d(16)*dexp(-hu)
      
      q = qb + q1*he + q2*(he**d(11))
      q = q*5.672213d0*(he**f)*(1.d0+d(17)*hh+d(18)*rw+d(19)*ms*ms)

      v(8) = qb
      v(9) = w1
      v(10) = w2
      v(11) = q1
      v(12) = q2
      
      return
      end            
      

      subroutine fitg8 (n, d, m, v, q, rel)
      ! fitg8
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      hh = v(6)
      hu = v(7)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+dexp((rw-d(8))/d(9))
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*dsqrt(sg)

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      
      fe = d(10) +d(11)*dexp(-ms) +d(12)*dexp(-mu) +d(13)*dexp(-hh)
     +          +d(14)*dexp(-hu)
      if (fe .le. 0.d0) then
        rel = rel - 1.d6*fe
        fe = 0.d0
      endif
      
      fc = 1.d0 +d(15)*hh +d(16)*rw +d(17)*hu
      if (fc .le. 0.d0) then
        rel = rel - 1.d6*fc
        fe = 0.d0
      endif
      
      q = qb + (q1 + q2)*he
      q = q*5.672213d0*(he**fe)*fc

c      v(8) = qb
      v(9) = w1
      v(10) = w2
c      v(11) = q1
c      v(12) = q2
      v(11) = q1*he
      v(12) = q2*he
      v(8) = qb + v(11)
      
      return
      end            
      

      subroutine fitg9 (n, d, m, v, q, rel)
      ! fitg9
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      hh = v(6)
      hu = v(7)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+dexp((rw-d(8))/d(9))
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*dsqrt(sg)

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      q = qb + q1*he + q2*he**d(10)
      
      fe = d(11) +d(12)*dexp(-ms) +d(13)*dexp(-mu) +d(14)*dexp(-hh)
     +          +d(15)*dexp(-hu) +d(16)*dexp(-rw)
      if (fe .le. 0.d0) then
        rel = rel - 1.d6*fe
        fe = 0.d0
      endif
      
      fc = 1.d0 +d(17)*hh +d(18)*rw +d(19)*hu
      if (fc .le. 0.d0) then
        rel = rel - 1.d6*fc
        fe = 0.d0
      endif
      
      q = q*5.672213d0*(he**fe)*fc

c      v(8) = qb
      v(9) = w1
      v(10) = w2
c      v(11) = q1
c      v(12) = q2
      v(11) = q1*he
      v(12) = q2*he
      v(8) = qb + v(11)
      
      return
      end            
      

      subroutine fith0 (n, d, m, v, q, rel)
      ! fith0
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      hh = v(6)
      hu = v(7)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+dexp((rw-d(8))/d(9))
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*dsqrt(sg)

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      q = qb + q1*he + q2*he**d(10)
      
      fe = d(11) +d(12)*dexp(-ms) +d(13)*dexp(-mu) +d(14)*dexp(-hh)
     +          +d(15)*dexp(-hu) +d(16)*dexp(-rw)
      if (fe .le. 0.d0) then
        rel = rel - 1.d6*fe
        fe = 0.d0
      endif
      
      fc = 1.d0 +d(17)*hh +d(18)*rw +d(19)*mu
      if (fc .le. 0.d0) then
        rel = rel - 1.d6*fc
        fe = 0.d0
      endif
      
      q = q*5.672213d0*(he**fe)*fc

c      v(8) = qb
      v(9) = w1
      v(10) = w2
c      v(11) = q1
c      v(12) = q2
      v(11) = q1*he
      v(12) = q2*he
      v(8) = qb + v(11)
      
      return
      end            
      

      subroutine fith1 (n, d, m, v, q, rel)
      ! fith1
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      he = v(1)
      b =  v(2)
      rw = v(3)
      ms = v(4)
      mu = v(5)
      hh = v(6)
      hu = v(7)
      sg = ms*ms+mu*mu
            
      qb = d(1)*b
      
      w2 = 1.d0+dexp((rw-d(8))/d(9))
c      if (dabs(w2) .lt. 0.01d0) then
c         rel = rel + 1.d8*(1.d-4 - w2*w2) + 10.d0
c         return
c      endif
      w2 = d(6) +d(7)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = w2*dsqrt(sg)

      w1 = 1.d0+dexp((w2-d(4))/d(5))
c      if (w1 .lt. 1.d0) then
c         rel = rel + 1.d4*(1.d0 - w1)
c      endif
      w1 = d(2) +d(3)/w1
      if (w1 .lt. 0.d0) then
         rel = rel - 1.d6*w1
         w1 = 0.d0
      endif
      q1 = w1*ms
      q = qb + q1*he + q2*he**d(10)
      
      fe = d(11) +d(12)*dexp(-ms) +d(13)*dexp(-mu) +d(14)*dexp(-hh)
     +          +d(15)*dexp(-hu) +d(16)*dexp(-rw)
      if (fe .le. 0.d0) then
        rel = rel - 1.d6*fe
        fe = 0.d0
      endif
      
      fc = 1.d0 +d(17)*hh +d(18)*rw +d(19)*dsqrt(sg)
      if (fc .le. 0.d0) then
        rel = rel - 1.d6*fc
        fe = 0.d0
      endif
      
      q = q*5.672213d0*(he**fe)*fc

c      v(8) = qb
      v(9) = w1
      v(10) = w2
c      v(11) = q1
c      v(12) = q2
      v(11) = q1*he
      v(12) = q2*he
      v(8) = qb + v(11)
      
      return
      end