program scaletesta1

      implicit real*8 (a-h, o-z)
      parameter (n = 9, nd = 2500)
      character*8 fil1, fil2
      real*8 mu, ms
      
      fil1 = 'scaletst'
      fil2 = 'my1607a3'
      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
      d1 = 3.4978
      d2 = -0.024658
      d3 = -0.3336
      d4 = 0.67084
      d5 = -0.34901
      d6 = -0.22005
c      d7 = -6.0681
c      d8 = 0.38644
c      d9 = 0.64499
      i = 0

   20 read(1, *, end=100) q, he, hh, th, hu, b, fi
      qpi = q/dsqrt(gr*he**5)
      hb = he/b
      w = he/(hh+he)
      wm = 1.d0-w
      mu = 1.d0/dtan(th)
      ms = dtan(fi)
      if (ms .le. 0.0) ms = 1.d-10
      
      qpih = d1*((he/(hu+he))**d2)*((he/(hh+he))**d3)*(hb**(-d4))
     +         *(dcos(fi)**d5)*(dsin(th)**d6)
     
      qh = qpih*dsqrt(gr*he**5)
      if (qh .le. 0.d0) then
         qh = 1.d-10
         qph = 1.d-10
         rel = rel+2.d0
         relpi = relpi+2.d0
      else
         rel1 = dabs(qh-q)/dsqrt(dabs(qh*q))
         rel = rel + rel1
         relpi1 = dabs(qpih-qpi)/dsqrt(dabs(qpih*qpi))
         relpi = relpi + relpi1
      end if
      std = std + rel1*rel1
      stdpi = stdpi + relpi1*relpi1
      write (2, '(12g14.5, f7.4, 2g14.5, f7.4)') 
     +   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