* sbf2flow.f:   Q+D to implement flow model derived in SBF II
*
* Syntax: sbf2flow x y z => vr sigmav vx vy vz
*
* Rev 1.0 - 2 July 1999 - John Tonry
*
      character*80 line
      if(iargc().lt.3) then
         write(6,*) 'Syntax: sbf2flow SGX SGY SGZ'
         stop
      end if

      call getarg(1,line)
      read(line,*) x
      call getarg(2,line)
      read(line,*) y
      call getarg(3,line)
      read(line,*) z

      call sbf2flow(x,y,z,vx,vy,vz,vr,vsig)

      write(6,1000) nint(vr), nint(vsig), nint(vx), nint(vy), nint(vz)
 1000 format(8i8)

      stop
      end
      
      subroutine sbf2flow(x,y,z,vx,vy,vz,vr,vsig)
*
* Rev 1.0 - 2 July 1999 - John Tonry
*
      ezp(x) = exp(amax1(-20.0,x))

      if(x.eq.0.and.y.eq.0.and.z.eq.0) x = 0.001
* Cosmology parameters
      h0 = 78.4
      wx = -55
      wy = 143
      wz =  -8
      omega = 0.2

* Common parameters
      rcore = 2
      thermal = 187

* VA parameters
      xv = -3.1
      yv = 16.6
      zv = -1.6
      deltav = 0.974
      gammav = 1.5
      rcutv = 11.7
      sigv = 650

* GA parameters
      xg = -36.7
      yg =  14.6
      zg = -17.1
      deltag = 0.781
      gammag = 2.0
      rcutg = 49.5
      sigg = 500

* Fornax parameters
      xf =  -1.9
      yf = -15.0
      zf = -13.4
      sigf = 235

* Quadrupole parameters
      rquad = 50
      qxx =  2.57
      qxy =  2.04
      qxz = -7.87
      qyz = -3.63
      qzz =  8.55

**************************************************

* Peculiar velocity
      vx = wx
      vy = wy
      vz = wz

      ivrw = nint((vx*x + vy*y + vz*z) / sqrt(x*x+y*y+z*z))

* Hubble flow contribution
      vx = vx + h0*x
      vy = vy + h0*y
      vz = vz + h0*z

      ivrh = nint(h0*sqrt(x*x+y*y+z*z))

* Virgo Attractor
      distv = amax1(0.001, sqrt(xv**2+yv**2+zv**2))
      d3 = (distv/rcore)*(distv/rcore)*(distv/rcore)
      delta0 = deltav*d3/(ezp(-distv/rcutv)*((1+d3)**(1-gammav/3)-1))

      ratv = amax1(0.001, sqrt((x-xv)**2+(y-yv)**2+(z-zv)**2))
      d3 = (ratv/rcore)*(ratv/rcore)*(ratv/rcore)
      delta = delta0/d3*(ezp(-ratv/rcutv)*((1+d3)**(1-gammav/3)-1))

      uinfall = h0/3 * ratv * omega**0.6 * delta * (1+delta)**-0.25
      
      ux = -uinfall * (x-xv)/ratv
      uy = -uinfall * (y-yv)/ratv
      uz = -uinfall * (z-zv)/ratv

      vx = vx + ux
      vy = vy + uy
      vz = vz + uz

      ivrv = nint((ux*x + uy*y + uz*z) / sqrt(x*x+y*y+z*z))

* Great Attractor
      distg = amax1(0.001, sqrt(xg**2+yg**2+zg**2))
      d3 = (distg/rcore)*(distg/rcore)*(distg/rcore)
      delta0 = deltag*d3/(ezp(-distg/rcutg)*((1+d3)**(1-gammag/3)-1))

      ratg = amax1(0.001, sqrt((x-xg)**2+(y-yg)**2+(z-zg)**2))
      d3 = (ratg/rcore)*(ratg/rcore)*(ratg/rcore)
      delta = delta0/d3*(ezp(-ratg/rcutg)*((1+d3)**(1-gammag/3)-1))

      uinfall = h0/3 * ratg * omega**0.6 * delta * (1+delta)**-0.25
      
      ux = -uinfall * (x-xg)/ratg
      uy = -uinfall * (y-yg)/ratg
      uz = -uinfall * (z-zg)/ratg

      vx = vx + ux
      vy = vy + uy
      vz = vz + uz

      ivrg = nint((ux*x + uy*y + uz*z) / sqrt(x*x+y*y+z*z))

* Quadrupole
      ratq = sqrt(x*x+y*y+z*z)
      cut = ezp(-0.5*ratq*ratq/(rquad*rquad))
      qyy = -qxx - qzz
      ux = cut*(x*qxx + y*qxy + z*qxz)
      uy = cut*(x*qxy + y*qyy + z*qyz)
      uz = cut*(x*qxz + y*qyz + z*qzz)

      vx = vx + ux
      vy = vy + uy
      vz = vz + uz

      ivrq = nint((ux*x + uy*y + uz*z) / sqrt(x*x+y*y+z*z))

* Radial component
      vr = (vx*x + vy*y + vz*z) / sqrt(x*x+y*y+z*z)

* Thermal velocity
      vsig = thermal*thermal
      vsig = vsig + sigv*sigv*ezp(-ratv*ratv/(rcore*rcore))
      vsig = vsig + sigg*sigg*ezp(-ratg*ratg/(rcore*rcore))
* Fornax
      ratf = amax1(0.001, sqrt((x-xf)**2+(y-yf)**2+(z-zf)**2))
      vsig = vsig + sigf*sigf*ezp(-ratf*ratf/(rcore*rcore))
      vsig = sqrt(vsig)

* Details (if desired)
C      write(6,1000) ivrh, ivrv, ivrq, ivrw, ivrg, nint(vr), nint(vsig)

      return
      end