c Last time edited 27.12.2000 by Risto Hanninen
c
c     Subroutine that calculates the hydrodynamic parametes for 3He-A.
c     input:
c         p = pressure in bars
c         temp = temperature relative to the critical temperature
c     output:
c         rhop=rho_parallel is given in units of rho(hbar/2m_3)^2
c         lambdam=lambda_m is given in units of chi_m
c         other parameters are all given in units of rhop
c         rhos = rho_perp
c         rho0 = -(rho_parallel-rho_perp)
c         CC=C,C0,Ks,Kt,Kb,K1=K5,K2=K6
c     subroutines:
c         fermpar, included in this file
c         crossc,  in separate file

      subroutine hydrop(p,temp,rhos,rho0,CC,C0,Ks,Kt,Kb,K1,K2,
     &                  rhop,lambdam)

      implicit none
      double precision p,F1s,F0a,F1a,xi0,Tc0,meff
      double precision temp,gap,alpha,beta,gamma,fs,fa,be2,al4,YY
      double precision rhos,rhop,rho0,CC,C0,Ks,Kt,Kb,K1,K2,lambdam

      call fermpar(p,F1s,F0a,F1a,xi0,Tc0)
      call crossc(temp,gap,alpha,beta,gamma)
c     some test values (used by M.R. Williams, 1979)
c      F1s = 15.66
c      F1a = -0.54
c      F0a = -0.7375
      fs = F1s/3d0
      fa = F1a/3d0
      meff = 1d0+fs
      be2 = 2d0*beta
      al4 = 4d0*alpha

      if (temp.ge.1.) then
         rhop = 0.0
         rhos = 2.0
         rho0 = 1.0
         CC = 0.5
         C0 = 1.0
         Ks = 0.5
         Kt = 0.5
         Kb = 1.5
         K1 = 1.0
         K2 = 2.0
         lambdam = 0.0
         goto 200
      end if

      rhop = be2/(1d0+fs*(1d0-be2))
      rhos = al4/(1d0+fs*(1d0-al4))
      rho0 = rhos-rhop

      CC = beta/(1d0+fs*(1d0-al4))
      C0 = CC+0.5*rhop

      Ks = alpha/meff
      Kt = ((alpha+be2)/3d0+fs*beta**2/(1d0+fs*(1d0-be2)))/meff
      Kb = (beta+2*gamma/3d0+fs*beta**2/(1d0+fs*(1d0-al4)))/meff

      K1 = 2d0*beta*(1d0+fa)/((1d0+fa*(1d0-be2))*meff)
      K2 = 4d0*alpha*(1d0+fa)/((1d0+fa*(1d0-al4))*meff)

      YY = 1-2*(al4+beta)/3.0
      lambdam = (1-YY)/(1+F0a*YY)
      
      rhos = rhos/rhop
      rho0 = rho0/rhop
      CC = CC/rhop
      C0 = C0/rhop
      Ks = Ks/rhop
      Kt = Kt/rhop
      Kb = Kb/rhop
      K1 = K1/rhop
      K2 = K2/rhop

 200  continue
      write(*,'(A,F6.3,A)') 'Pressure: ',p,' bars'
      write(*,'(A,F7.4,A,F7.4,A,F7.4)') 
     &  'F1s  = ',F1s,' F0a = ',F0a,' F1a = ',F1a
      write(*,'(A,F7.4,A,F7.4,A)') 'meff = ',meff,' xi0 = ',xi0,' nm'
      write(*,*) 'Hydrodynamic parameters:'
      write(*,'(A,F7.5,A,F7.5)') 'rhop = ',rhop,' lamb = ',lambdam
      write(*,'(A,F7.5,A,F7.5)') 'rhos = ',rhos,' rho0 = ',rho0
      write(*,'(A,F7.5,A,F7.5)') 'CC   = ',CC  ,' C0   = ',C0
      write(*,'(A,F7.5,A,F7.5)') 'Ks   = ',Ks  ,' Kt   = ',Kt
      write(*,'(A,F7.5,A,F7.5)') 'K1   = ',K1  ,' K2   = ',K2
      write(*,'(A,F7.5)') 'Kb   = ',Kb

      end
      
c *************************************************************      

c     Calculates F1s, F1a, xi0 and Tc0 for pressure given  
c     Taken from Halperin, W.P., and E. Varoquax, Helium Three 
c     ed. W.P. Halperin, and L.P. Pitaevskii, Elsevier (1990). 
      subroutine fermpar(p,F1s,F0a,F1a,xi0,Tc0)

      implicit none

      double precision pi,p,hbar,kB,V,vF,c1,CnRT,mrtio,xi0,Tc0
      double precision F1s,F0a,F1a,F2a
      PARAMETER        (pi=3.141592653589793d0)

      hbar = 1.0545d-34
      kB = 1.3805d-23

c     molar volume
      V=36.82-1.185*p+0.08609*p**2-4.187d-3*p**3+1.071d-4*p**4
     &     -1.082d-6*p**5
c     Fermi velocity
      vF=59.03-1.967*p+0.08349*p**2-2.304d-3*p**3+2.553d-5*p**4
c     First sound
      c1=183.1+17.00*p-0.8940*p**2+0.03824*p**3-8.943d-4*p**4
     &     +8.484d-6*p**5
c     specific heat C/nRT
      CnRT=2.78+0.07151*p-1.678d-3*p**2+5.342d-5*p**3-6.160d-7*p**4
      mrtio=2.8+0.1292*p-3.188d-3*p**2+9.372d-5*p**3-1.030d-6*p**4
c     Tc0 in mK
      Tc0=0.9294+0.1387*p-6.930d-3*p**2+2.569d-4*p**3-5.725d-6*p**4
     &     +5.301d-8*p**5
      xi0=1000*hbar*vF/(2*pi*kB*Tc0)
c      xi0 in nanometers (nm)
      xi0=1.0d9*xi0

      F1s = 3.0*(mrtio-1.0)
      F0a = -7.007d-1-6.232d-3*p+2.057d-4*p**2-1.823d-6*p**3
      F1a = -5.678d-1-4.753d-2*p+1.791d-3*p**2-2.273d-5*p**3
      F2a = -8.774d-1+4.753d-2*p-6.399d-4*p**2

      end