program scaletesta2

      implicit real*8 (a-h, o-z)
      parameter (n = 13, nd = 2)
      character*8 fil1, fil2
      real*8 mu, ms, d(n), v(10)
      
      fil1 = 'scaleta0'
      fil2 = 'my2507a3'
      open (1, file=fil1//'.prn')
      open (2, file=fil2//'.prn')
      write (2, *) 'Hb,  W,  Mu,  Ms,  He,  Hh,  Hu,  B,  Theta,  Fi, ', 
     +               'Q,  Qh,  Rel-Error,  Qpi,  Qpih,  Rel-Pi-Err'

      read (1, *)
      gr = 32.174d0
      rel = 0.d0
      relpi = 0.d0
      std = 0.d0
      stdpi = 0.d0
      i = 0
      d(1) = 8.5073d0
      d(2) = 0.65821d0
      d(3) = -.062313d0
      d(4) = 0.011479d0
      d(5) = 0.0093485d0
      d(6) = -0.25833d0
      d(7) = 5.089d0
      d(8) = 0.1026d0
      d(9) = -0.0078405d0
      d(10) = -1.1749d0
      d(11) = -6.7532d0
      d(12) = -0.15133d0
      d(13) = 1.0727d0

   20 read(1, *, end=100) q, he, hh, th, hu, b, fi
      qpi = q/(b*dsqrt(gr*he**3))
      hb = he/b
c      w = he/(hh+he)
      w = he*he/(he*he+b*hh)
      if (w .ge. 1.d0) w = 1.d0 - 1.d-10
c      wm = 1.d0-w
      mu = 1.d0/dtan(th)
      ms = dtan(fi)
      if (hb .le. 0.0) hb = 1.d-10
      if (w .le. 0.0) w = 1.d-10
      if (ms .le. 0.0) ms = 1.d-10
      if (mu .le. 0.0) mu = 1.d-10
      v(1) = hb
      v(2) = w
      v(3) = mu
      v(4) = ms
      
      call fit (n, d, 4, v, qpih)
     
      qh = qpih*b*dsqrt(gr*he**3)
      relpi1 = dabs(qpih-qpi)/dsqrt(dabs(qpih*qpi))
      rel1 = dabs(qh-q)/dsqrt(dabs(qh*q))
      if (qh .le. 0.d0) then
         rel1 = 2.d0 + rel1*rel1
         relpi = 2.d0 +relpi1*relpi1
      endif
      relpi = relpi + relpi1
      stdpi = stdpi + relpi1*relpi1
      rel = rel + rel1
      std = std + rel1*rel1
      
      write (2, '(16(g14.5,1h,))') 
     +   hb, w, mu, ms, he, hh, hu, b, th, 
     +   fi, q, qh, rel1, qpi, qpih, relpi1
      
      i = i + 1
      go to 20
      
  100 continue
      xn = dble(i)
      std = dsqrt((xn*std-rel*rel)/(xn*(xn-1)))
      stdpi = dsqrt((xn*stdpi-relpi*relpi)/(xn*(xn-1)))
      rel = rel/xn
      relpi = relpi/xn
      write (2, *) i, relpi, stdpi, rel, std
      print *, i, relpi, stdpi, rel, std
      
      close(1)
      close(2)
      pause

      end
      
      
      
            

      subroutine fit (n, d, m, v, q)
      
      implicit real*8 (a-h, o-z)
      dimension d(n), v(m)
      real*8 mu, ms
      
      hb = v(1)
      w = v(2)
c      wm = 1.d0 - w
      mu = v(3)
      ms = v(4)
      qa8 = 0.47175 +hb*( 0.20017 +w*(-2.0328+w*(8.0383+w*(-9.115
     +      +w*3.2259))) )*ms 
     + +(hb**0.62922)*( 0.49922 +w*(1.5502+w*(-7.1752+w*(7.7394
     +      +w*(-2.5628)))) )
     +  *((ms*ms+mu*mu)**0.36655)
      q = ( d(1) + d(2)*(ms**d(3))*(mu**d(4))/
     +                   ( 1.d0+d(5)*(ms**d(6))*(mu**d(7)) ) )*qa8
      q = ( d(8) + d(9)*(hb**d(10))/( 1.d0+d(11)*(hb**d(12)) ) )*q
      q = q**d(13)
      
      return
      end            


      subroutine fitb1 (n, d, m, v, q)
      
      implicit real*8 (a-h, o-z)
      dimension d(n), v(m)
      real*8 mu, ms
      
      hb = v(1)
      w = v(2)
c      wm = 1.d0 - w
      mu = v(3)
      ms = v(4)
      qa8 = 0.47175 +hb*( 0.20017 +w*(-2.0328+w*(8.0383+w*(-9.115
     +      +w*3.2259))) )*ms 
     + +(hb**0.62922)*( 0.49922 +w*(1.5502+w*(-7.1752+w*(7.7394
     +      +w*(-2.5628)))) )
     +  *((ms*ms+mu*mu)**0.36655)
      q = ( d(1) + d(2)*(ms**d(3))*(mu**d(4))/
     +                   ( 1.d0+d(5)*(ms**d(6))*(mu**d(7)) ) )*qa8
      q = ( d(8) + d(9)*(hb**d(10))/( 1.d0+d(11)*(hb**d(12)) ) )*q
      q = q**d(13)
      
      return
      end            


      subroutine fita8 (n, d, m, v, q)
      
      implicit real*8 (a-h, o-z)
      dimension d(n), v(m)
      real*8 mu, ms
      
      hb = v(1)
      w = v(2)
c      wm = 1.d0 - w
      mu = v(3)
      ms = v(4)
      q = d(1) +hb*( d(2) +w*(d(3)+w*(d(4)+w*(d(5)+w*d(6)))) )*ms 
     + +(hb**d(7))*( d(8) +w*(d(9)+w*(d(10)+w*(d(11)+w*d(12)))) )
     +  *((ms*ms+mu*mu)**d(13))
      
      return
      end            


      subroutine fit01 (n, d, m, v, q)
      
      implicit real*8 (a-h, o-z)
      dimension d(n), v(m)
      real*8 mu, ms
      
      hb = v(1)
      w = v(2)
      wm = 1.d0 - w
      mu = v(3)
      ms = v(4)
      q = d(1) +d(2)*(hb**d(3))*(w**d(4))*(ms**d(5)) 
     +         +d(6)*(hb**d(7))*(wm**d(8))*(mu**d(9))
      
      return
      end