program simanj7

      implicit real*8 (a-h, o-z)
      parameter (n = 11, nd = 3000)
c      parameter (n = 16, 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 = 'jl1207b1'
      fil3 = 'jl1207b2'
      
      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'
      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)

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, ', 
     +           'Qb, 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)
         
         call fit (n, x, 10, 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(6)
         w1 = v(7)
         w2 = v(8)
         q1 = v(9)
         q2 = v(10)
         
               
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)
      ! 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 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)
         
         call fit (n, x, 10, 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)
      
      parameter (n = 11, neps = 4)
c      parameter (n = 16, 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 = 804
      iseed2 = 77
      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) = .461d0
      ub(1) = 10.d0
      
       x(2) = .235d0
      ub(2) = 1.d0
      
c      lb(3) = -5.d0
       x(3) = .090d0
      ub(3) = 10.d0
      
      lb(4) = -1.d0
       x(4) = .131d0
      ub(4) = 1.d0
      
      lb(5) = .0001d0
       x(5) = .0209d0
      ub(5) = 1.d0
      
      lb(6) = -1.d0
       x(6) = -.0704d0
      ub(6) = 1.d0
      
c      lb(7) = -5.d0
       x(7) = .776d0
      ub(7) = 2.d0
      
      lb(8) = -.5d0
       x(8) = .51d0
      ub(8) = 1.5d0
      
c      lb(9) = .5d0
       x(9) = .213d0
      ub(9) = 2.d0
      
      lb(10) = .1d0
       x(10) = .306d0
      ub(10) = .5d0
      
       x(11) = .768d0
      ub(11) = 1.5d0
      
c      lb(13) = .5d0
c      x(13) = 2.5d0
c      ub(13) = 5.d0
      
c      lb(14) = -.5d0
c      x(14) = 0.d0
c      ub(14) = .5d0

c      lb(15) = 1.d0
c      x(15) = 20.d0
c      ub(15) = 40.d0
      
c      lb(16) = 1.d0
c      x(16) = 20.1d0
c      ub(16) = 40.d0
c      x(17) = 0.d0
c      x(18) = 0.d0
      
      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            
      

      subroutine fitd8 (n, d, m, v, q, rel)
      ! fitd8
      
      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)
      
      q = d(1)
      
      c1 = 1.d0 + rw*d(4)
      if (dabs(c1) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c1))
         c1 = 1.d-6
      endif
      c1 = (d(3) + rw*d(5))/c1
      if (c1 .lt. 0.d0) then
         rel = rel - 100.d0*c1
         c1 = 0.d0
      end if
      
      c2 = 1.d0+rw*d(7)
      if (dabs(c2) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c2))
         c2 = 1.d-6
      end if
      c2 = (d(6) + rw*d(8))/c2
      if (c2 .lt. 0.d0) then
         rel = rel - 100.d0*c2
         c2 = 0.d0
      end if
      v(6) = c1
      v(7) = c2
      
      q = q + c1*hb*ms
     +      + c2*( hb*hb*(ms*ms+mu*mu) )**d(2) 
      
      return
      end            
      

      subroutine fitd9 (n, d, m, v, q, rel)
      ! fitd9
      
      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)
      
      q = d(1)
      
      c1 = 1.d0 + rw*(d(4)+rw*d(9))
      if (dabs(c1) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c1))
         c1 = 1.d-6
      endif
      c1 = (d(3) + rw*d(5))/c1
      if (c1 .lt. 0.d0) then
         rel = rel - 100.d0*c1
         c1 = 0.d0
      end if
      
      c2 = 1.d0+rw*(d(7)+rw*d(10))
      if (dabs(c2) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c2))
         c2 = 1.d-6
      end if
      c2 = (d(6) + rw*d(8))/c2
      if (c2 .lt. 0.d0) then
         rel = rel - 100.d0*c2
         c2 = 0.d0
      end if
      v(6) = c1
      v(7) = c2
      
      q = q + c1*hb*ms
     +      + c2*( hb*hb*(ms*ms+mu*mu) )**d(2) 
      
      return
      end            
      

      subroutine fite0 (n, d, m, v, q, rel)
      ! fite0
      
      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)
      
      q = d(1)
      
      c1 = 1.d0 + rw*(d(4)+rw*d(9))
      if (dabs(c1) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c1))
         c1 = 1.d-6
      endif
      c1 = (d(3) + rw*(d(5)+rw*d(11)))/c1
      if (c1 .lt. 0.d0) then
         rel = rel - 100.d0*c1
         c1 = 0.d0
      end if
      
      c2 = 1.d0+rw*(d(7)+rw*d(10))
      if (dabs(c2) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c2))
         c2 = 1.d-6
      end if
      c2 = (d(6) + rw*(d(8)+rw*d(12)))/c2
      if (c2 .lt. 0.d0) then
         rel = rel - 100.d0*c2
         c2 = 0.d0
      end if
      v(6) = c1
      v(7) = c2
      
      q = q + c1*hb*ms
     +      + c2*( hb*hb*(ms*ms+mu*mu) )**d(2) 
      
      return
      end            
      

      subroutine fite1 (n, d, m, v, q, rel)
      ! fite1
      
      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)
      
      q = d(1)
      
      c1 = 1.d0 + rw*(d(4)+rw*(d(9)+rw*d(13)))
      if (dabs(c1) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c1))
         c1 = 1.d-6
      endif
      c1 = (d(3) + rw*(d(5)+rw*d(11)))/c1
      if (c1 .lt. 0.d0) then
         rel = rel - 100.d0*c1
         c1 = 0.d0
      end if
      
      c2 = 1.d0+rw*(d(7)+rw*(d(10)+rw*d(14)))
      if (dabs(c2) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c2))
         c2 = 1.d-6
      end if
      c2 = (d(6) + rw*(d(8)+rw*d(12)))/c2
      if (c2 .lt. 0.d0) then
         rel = rel - 100.d0*c2
         c2 = 0.d0
      end if
      v(6) = c1
      v(7) = c2
      
      q = q + c1*hb*ms
     +      + c2*( hb*hb*(ms*ms+mu*mu) )**d(2) 
      
      return
      end            
      

      subroutine fite2 (n, d, m, v, q, rel)
      ! fite2
      
      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)
      
      q = d(1)
      
      c1 = 1.d0 + rw*(d(4)+rw*(d(9)+rw*d(13)))
      if (dabs(c1) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c1))
         c1 = 1.d-6
      endif
      c1 = (d(3) + rw*(d(5)+rw*(d(11)+rw*d(15))))/c1
      if (c1 .lt. 0.d0) then
         rel = rel - 100.d0*c1
         c1 = 0.d0
      end if
      
      c2 = 1.d0+rw*(d(7)+rw*(d(10)+rw*d(14)))
      if (dabs(c2) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c2))
         c2 = 1.d-6
      end if
      c2 = (d(6) + rw*(d(8)+rw*(d(12)+rw*d(16))))/c2
      if (c2 .lt. 0.d0) then
         rel = rel - 100.d0*c2
         c2 = 0.d0
      end if
      v(6) = c1
      v(7) = c2
      
      q = q + c1*hb*ms
     +      + c2*( hb*hb*(ms*ms+mu*mu) )**d(2) 
      
      return
      end            
      

      subroutine fite3 (n, d, m, v, q, rel)
      ! fite3
      
      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)
      
      q = d(1)
      
      c1 = 1.d0 + rw*( d(4)+rw*( d(9)+rw*( d(13)+rw*d(17))))
      if (dabs(c1) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c1))
         c1 = 1.d-6
      endif
      c1 = (d(3) + rw*(d(5)+rw*(d(11)+rw*d(15))))/c1
      if (c1 .lt. 0.d0) then
         rel = rel - 100.d0*c1
         c1 = 0.d0
      end if
      
      c2 = 1.d0+rw*(d(7)+rw*(d(10)+rw*(d(14)+rw*d(18))))
      if (dabs(c2) .lt. 1.d-6) then
         rel = rel + 1.d8*(1.d-6-dabs(c2))
         c2 = 1.d-6
      end if
      c2 = (d(6) + rw*( d(8)+rw*( d(12)+rw*d(16))))/c2
      if (c2 .lt. 0.d0) then
         rel = rel - 100.d0*c2
         c2 = 0.d0
      end if
      v(6) = c1
      v(7) = c2
      
      q = q + c1*hb*ms
     +      + c2*( hb*hb*(ms*ms+mu*mu) )**d(2) 
      
      return
      end            
      

      subroutine fitf0 (n, d, m, v, q, rel)
      ! fitf0
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      rw = v(2)
      ms = v(3)
      mu = v(4)
      ssu = ms*ms+mu*mu      
      q = d(1)
      
      w1 = d(3) + d(4)*rw + d(5)*dsqrt(rw)
      if (w1 .lt. 0.d0) then
         rel = rel - 100.d0*w1
         w1 = 1.d-10
      end if
      
      w2 = 1.d0 + (w1/d(8))**d(9)
      w2 = d(6) + d(7)/w2
c      if (dabs(c2) .lt. 1.d-6) then
c         rel = rel + 1.d8*(1.d-6-dabs(c2))
c         c2 = 1.d-6
c      end if
      sc = 1.d0 + d(10)*dexp(-d(11)*ms) + d(12)*dexp(-d(13)*ssu)
      if (sc .lt. 0.d0) then
         rel = rel - 100.d0*sc
         sc = 0.d0
      end if
      v(5) = w1
      v(6) = w2
      v(7) = sc
      
      q = q + w1*hb*ms
     +      + w2*( hb*hb*ssu )**d(2)
      q = sc*q 
      
      return
      end            
      

      subroutine fitf1 (n, d, m, v, q, rel)
      ! fitf1
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      rw = v(2)
      ms = v(3)
      mu = v(4)
      ssu = ms*ms+mu*mu      
      q = d(1)
      
      w1 = (d(3) + d(4)*rw*rw)/(1.d0 + d(5)*rw*rw)
      if (w1 .lt. 0.d0) then
         rel = rel - 100.d0*w1
         w1 = 1.d-10
      end if
      
      w2 = 1.d0 + (w1/d(8))**d(9)
      w2 = d(6) + d(7)/w2
c      if (dabs(c2) .lt. 1.d-6) then
c         rel = rel + 1.d8*(1.d-6-dabs(c2))
c         c2 = 1.d-6
c      end if
      sc = 1.d0 + d(10)*dexp(-d(11)*ms) + d(12)*dexp(-d(13)*ssu)
      if (sc .lt. 0.d0) then
         rel = rel - 100.d0*sc
         sc = 0.d0
      end if
      v(5) = w1
      v(6) = w2
      v(7) = sc
      
      q = q + w1*hb*ms
     +      + w2*( hb*hb*ssu )**d(2)
      q = sc*q 
      
      return
      end            
      

      subroutine fitf2 (n, d, m, v, q, rel)
      ! fitf2
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      rw = v(2)
      ms = v(3)
      mu = v(4)
      ssu = ms*ms+mu*mu      
      q = d(1)
      
      w1 = d(3) + d(4)*(rw**d(5))
      if (w1 .lt. 0.d0) then
         rel = rel - 100.d0*w1
         w1 = 1.d-10
      end if
      
      w2 = 1.d0 + (w1/d(8))**d(9)
      w2 = d(6) + d(7)/w2
      if (w2 .lt. 0.d0) then
         rel = rel - 100.d0*w2
         w2 = 1.d-10
      end if
      sc = 1.d0 + d(10)*dexp(-d(11)*ms) + d(12)*dexp(-d(13)*ssu)
      if (sc .lt. 0.d0) then
         rel = rel - 100.d0*sc
         sc = 0.d0
      end if
      v(5) = w1
      v(6) = w2
      v(7) = sc
      
      q = q + w1*hb*ms
     +      + w2*( hb*hb*ssu )**d(2)
      q = sc*q 
      
      return
      end            
      

      subroutine fitf3 (n, d, m, v, q, rel)
      ! fitf3
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      rw = v(2)
      ms = v(3)
      mu = v(4)
      ssu = ms*ms+mu*mu      
      q = d(1)
      
      w1 = d(3) + d(4)*(rw**d(5))
      if (w1 .lt. 0.d0) then
         rel = rel - 100.d0*w1
         w1 = 1.d-10
      end if
      
      w2 = 1.d0 + (w1/d(8))**d(9)
      w2 = d(6) + d(7)/w2
      if (w2 .lt. 0.d0) then
         rel = rel - 100.d0*w2
         w2 = 1.d-10
      end if
      sc = 1.d0 + d(10)*dexp(-d(11)*ms) 
     +           + d(12)*dexp(-d(13)*(ssu**d(14)))
      if (sc .lt. 0.d0) then
         rel = rel - 100.d0*sc
         sc = 0.d0
      end if
      v(5) = w1
      v(6) = w2
      v(7) = sc
      
      q = q + w1*hb*ms
     +      + w2*( hb*hb*ssu )**d(2)
      q = sc*q 
      
      return
      end            
      

      subroutine fitf4 (n, d, m, v, q, rel)
      ! fitf4
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      rw = v(2)
      ms = v(3)
      mu = v(4)
      ss = ms*ms+mu*mu      
      
      sc1 = 1.d0 + d(1)*dexp(-d(2)*ms) + d(3)*dexp(-d(4)*ss)
      if (sc1 .lt. 0.d0) then
         rel = rel - 100.d0*sc1
         sc1 = 0.d0
      end if
      
      w1 = d(5) + d(6)*(rw**d(7))
      if (w1 .lt. 0.d0) then
         rel = rel - 100.d0*w1
         w1 = 1.d-10
      end if
      
      sc2 = 1.d0 + d(8)*dexp(-d(9)*ms)
      if (sc2 .lt. 0.d0) then
         rel = rel - 100.d0*sc2
         sc2 = 0.d0
      end if
      
      w2 = 1.d0 + (w1*4.02107d0)**5.2024d0
      w2 = -0.17565d0 + .78212d0/w2
      if (w2 .lt. 0.d0) then
         rel = rel - 100.d0*w2
         w2 = 1.d-10
      end if
      
      sc3 = 1.d0 + d(10)*dexp(-d(11)*(ss**d(12)))
      if (sc3 .lt. 0.d0) then
         rel = rel - 100.d0*sc3
         sc3 = 0.d0
      end if
      
      v(5) = w1
      v(6) = w2
      v(7) = sc1
      v(8) = sc2
      v(9) = sc3
      
      q = 0.4674d0*sc1 + w1*hb*ms*sc2 +w2*(hb*hb*(ss**0.36714d0))*sc3
      
      return
      end            
      

      subroutine fitf5 (n, d, m, v, q, rel)
      ! fitf5
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      rw = v(2)
      ms = v(3)
      mu = v(4)
      ss = ms*ms+mu*mu      
      q = d(1)
      
      w1 = d(3) + d(4)*(rw**d(5))
      if (w1 .lt. 0.d0) then
         rel = rel - 100.d0*w1
         w1 = 1.d-10
      end if
      
      w2 = 1.d0 + (w1/d(8))**d(9)
      w2 = d(6) + d(7)/w2
      if (w2 .lt. 0.d0) then
         rel = rel - 100.d0*w2
         w2 = 1.d-10
      end if
      sc = 1.d0 + d(10)*dexp(-d(11)*ms) + d(12)*dexp(-d(13)*ss)
     +          + d(14)*dexp(-d(15)*ms-d(16)*ss)
      if (sc .lt. 0.d0) then
         rel = rel - 100.d0*sc
         sc = 0.d0
      end if
      v(5) = w1
      v(6) = w2
      v(7) = sc
      
      q = q + w1*hb*ms
     +      + w2*( hb*hb*ss )**d(2)
      q = sc*q 
      
      return
      end            
      

      subroutine fitf6 (n, d, m, v, q, rel)
      ! fitf6
      
      implicit real*8 (a-h, o-z)
      real*8 d(n), v(m) 
      real*8 mu, ms
      
      hb = v(1)
      rw = v(2)
      ms = v(3)
      mu = v(4)
      ss = ms*ms+mu*mu      
      q = d(1)
      
      w1 = d(3) + d(4)*(rw**d(5))
      if (w1 .lt. 0.d0) then
         rel = rel - 100.d0*w1
         w1 = 1.d-10
      end if
      if (d(4)*d(5) .lt. 0.d0) rel = rel - 100.d0*d(4)*d(5)
      
      w2 = 1.d0 + (w1/d(8))**d(9)
      w2 = d(6) + d(7)/w2
      if (w2 .lt. 0.d0) then
         rel = rel - 100.d0*w2
         w2 = 1.d-10
      end if
      sc = 1.d0 + d(10)*dexp(-d(11)*ms) + d(12)*dexp(-d(13)*ss)
     +          + d(14)*dexp(-d(15)*ms-d(16)*ss)
      if (sc .lt. 0.d0) then
         rel = rel - 100.d0*sc
         sc = 0.d0
      end if
      v(5) = w1
      v(6) = w2
      v(7) = sc
      
      q = q + w1*hb*ms
     +      + w2*( hb*hb*ss )**d(2)
      q = sc*q 
      
      return
      end            
      

      subroutine fitf7 (n, d, m, v, q, rel)
      ! fitf7
      
      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*he*5.6772213d0*dsqrt(he)
      
      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**d(6))*(he**d(7))
      
      w2 = 1.d0+(rw/d(10))**d(11)
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(8) +d(9)/w2
      if (w2 .lt. 0.d0) w2 = 0.d0
      q2 = 5.6772213d0*(sg**d(12))*(he**d(13))
      
      q = qb + w1*q1 + w2*q2

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

      subroutine fitf8 (n, d, m, v, q, rel)
      ! fitf8
      
      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 + q2)*he)*5.672213d0*he*dsqrt(he)

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

      subroutine fitf9 (n, d, m, v, q, rel)
      ! fitf9
      
      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))
      
c      q = qb + q1*he + q2*(he**d(11))*(b**(1.d0-d(11)))
      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 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