c ######################################################################
c #
c #                   TUNGSTEN PROPERTIES PACKAGE
c #                   -----------------------------
c #
c # Contains functions for the calculation of the thermo-physical
c # properties of pure Tungsten (W)
c #
c ######################################################################

c ######################################################################
      real function dW(T)
c ######################################################################
c #
c # Density of Manganese
c #
c # Range: 0 <= T <= 3400 K
c #
c #                        References
c #                        ----------
c # http://en.wikipedia.org/wiki/Tungsten
c #
c # variable          I/O    meaning                           units
c # --------------------------------------------------------------------
c #   T               x      absolute temperature              K
c #   dW                x    density                           Kg/m**3
c #
c #
c # Author : L.Bottura at Cryosoft
c # Version: 1.0   March 2013
c #
c ######################################################################
      implicit none
c * external variables
      real     T
c * fit variables      

c * local variables

c *      
      dW = 19250.0
c *

      return
      end
c ######################################################################
      real function cW(T)
c ######################################################################
c #
c # Specific heat of Tungsten
c #
c # Range: 1 <= T <= 3400 K
c #
c #                        References
c #                        ----------
c # G.K. White, S.J. Collocott, Heat Capacity of Reference Materials: 
c # Cu and W, J. Phys. Chem. Ref. Data, 13, No 4, 1251, 1984
c #
c # variable          I/O    meaning                           units
c # --------------------------------------------------------------------
c #   T               x      absolute temperature              K
c #   cW                x    specific heat                     J/Kg K
c #
c # Author : L.Bottura at Cryosoft
c # Version: 1.0   March 2013
c #
c ######################################################################
      implicit none
c * external variables
      real     T
c * fit variables
      real     T0
      real     Tmin,Tmax
      real     a1,a3,a5,a7
      real     AA,BB,CC,DD,a,b,c,d,na,nb,nc,nd
      data     T0  /  30.88000301 / 
      data     a1  /  0.005562986 / ,  a3  /  0.177964e-3 / , 
     &         a5  /  9.94444e-08 / ,  a7  / -1.13203e-11 /
      data     AA  / -0.211376751 / ,  BB  /  0.0742226861 / ,
     &         CC  / -0.00017876104 / ,  DD  /  2.92E-10    /
      data     a   /  22.12680831 / ,  b   /  52.03300092 / ,  
     &         c   /  69.91811762 / ,  d   /  3854.004399 /
      data     na  /  2.179684475 / ,  nb  /  2.310745633 / , 
     &         nc  /  3.427955022 / ,  nd  /  3.562898947 /
      data     Tmin / 1.0/, Tmax / 3400.0/
      save    
c * local variables   
      real     TT

c *
      TT=T
      TT=max(TT,Tmin)
      TT=min(TT,Tmax)
      
      if (TT.le.T0) then
         cW = a1*TT + a3*TT**3 + a5*TT**5 + a7*TT**7            
      elseif(TT.gt.T0) then
         cW = AA*TT**na /(1+TT/a)**na + BB*TT**nb / (1+TT/b)**nb +
     &        CC*TT**nc /(1+TT/c)**nc + DD*TT**nd / (1+TT/d)**nd
      endif

c *
      return
      end
c ######################################################################
      real function kW(T,B,RRR)
c ######################################################################
c #
c # Thermal conductivity of Tungsten
c #
c # Range: 2 <= T <= 3600 K, no B dependence, no RRR dependence
c #
c #                        References
c #                        ----------
c # C.Y. Ho, R.W. Powell, P.E.Liley, Thernal Conductivity of the 
c # Elements, J. Phys. Chem. Ref. Data, 1 (2), 1972
c #
c # variable          I/O    meaning                           units
c # --------------------------------------------------------------------
c #   T               x      absolute temperature              K
c #   B               x      magnetic field                    T
c #   RRR             x      residual resistivity ratio        -
c #   kW                x    thermal conductivity              W/m K
c #
c #
c # Author : L.Bottura at Cryosoft
c # Version: 2.0   December 2020
c #
c ######################################################################
      implicit none
c * external variables
      real     T,B,RRR
c * fit variables      
      real     Tmin,Tmax
      real     k0,alpha,beta,n,m
      real     Tc
      data     k0    / 108.0894507 / , alpha / 1515.96341   /,
     &         beta  / 0.000142973 /
      data     n     /  1.00409115 / , m     / 3.598429945 /                         
      data     Tc    / 27.8641308 / 
      data     Tmin / 1.0/, Tmax / 3600.0/
      save
c * local variables
      real     TT,k0exp,knum,kden,k1,k2
c *
      TT=T
      TT=min(TT,Tmax)
      TT=max(TT,Tmin) 

c * catch overflow in exponential
      k0exp = max(min((TT-Tc)/Tc,30.0),-30.0)

c *
c      knum = alpha*TT**n
c      kden = 1.0+beta*TT**m
c      write(6,*) 'TT ',TT,' knum ',knum,' kden ',kden
c      k1   = knum/kden
c      write(6,*) 'k1 ',k1
c      k2   = k0 * (1.0-exp(-k0exp))
c      write(6,*) ' k2 ',k2
      kW = (alpha*TT**n) / (1.0+beta*TT**m) + k0 * (1.0-exp(-k0exp))
     
c *
      return
      end
c ######################################################################
      real function rW(T,B,RRR)
c ######################################################################
c #
c # Electrical resistivity of Tungsten
c #
c # Range: 1 <= T <= 3600 K,  0 <= B <= ?? T,  1.5 <= RRR <= 3000
c #
c #                        References
c #                        ----------
c # W.J. De Haas, J. De Nobel, The Thermal and The Electrical Resistance
c # of a Tungsten Single-crystal at Low Temperatures and in Magnetic 
c # Fields, Physica 5, 449, 1938
c #
c # variable          I/O    meaning                           units
c # --------------------------------------------------------------------
c #   T               x      absolute temperature              K
c #   B               x      magnetic field                    T
c #   RRR             x      residual resistivity ratio        -
c #   rW                x    resistivity                       Ohm m
c #
c # Author : L.Bottura at CryoSoft
c # Version: 1.0   March 2013
c #
c ######################################################################
      implicit none
c * external variables
      real     T,B,RRR
c * fit variables
      real     Tmin,Tmax
      parameter(Tmin   = 1,   Tmax   = 3600.0)
      real     Bmin,Bmax
      parameter(Bmin   = 0.0, Bmax   = 30.0)
      real     RRRmin,RRRmax
      parameter(RRRmin = 1.5, RRRmax = 3000.0)
      real     rho273
      parameter(rho273=4.50e-8)
      real     p1,p2,p3,p4,p5,p6,p7
      parameter(p1=2.6037e-17, p2=4.5661, p3=0.1974e11,
     &          p4=-1.2259, p5=10, p6=1, p7=0.333)
c * local variables   
      real     TT,BB,R,rhozero,arg,rhoi,rhoi0,rho0,magr
      real     magrW

c *
      TT=T
      TT=max(TT,Tmin)
      TT=min(TT,Tmax)

      BB=B
      BB=max(BB,Bmin)
      BB=min(BB,Bmax)

      R =RRR
      R =max(R,RRRmin)
      R =min(R,RRRmax)
	
c * resistivity at absolute zero
      rhozero = rho273/(R-1.0)

c * resistivity at the temperature t
      arg     = min((p5/TT)**p6,30.0)
      rhoi    = p1*TT**p2/(1.+p1*p3*TT**(p2+p4)*exp(-arg))
      rhoi0   = p7*rhoi*rhozero/(rhoi+rhozero)
      rho0    = rhozero+rhoi+rhoi0

c * transverse magneto-resistance factor
      magr    = magrW(TT,BB,R)

c * resistivity
      rW      = magr * rho0
	
c *
      return
      end
c ######################################################################
c #
c # Auxiliary functions and calculations
c #
c ######################################################################

c ######################################################################
      real function magrW(T,B,RRR)
c ######################################################################
c #
c # Magneto-resistivity factor of Tungsten, given as the ratio between
c # resistivity in transverse magnetic field B to resistivity at zero 
c # field
c #
c # Range: 1 <= T <= 3600 K,  0 <= B <= ?? T,  1.5 <= RRR <= 3000
c #
c #                        References
c #                        ----------
c # W.J. De Haas, J. De Nobel, The Thermal and The Electrical Resistance
c # of a Tungsten Single-crystal at Low Temperatures and in Magnetic 
c # Fields, Physica 5, 449, 1938
c #
c # variable          I/O    meaning                           units
c # --------------------------------------------------------------------
c #   T               x      absolute temperature              K
c #   B               x      magnetic field                    T
c #   RRR             x      residual resistivity ratio        -
c #   magrW             x    magneto-resistivity factor        -
c #
c # Author : L.Bottura at Cryosoft
c # Version: 1.0   March 2013
c #
c ######################################################################
      implicit none
c * external variables
      real     T,B,RRR
c * fit variables
      real     Tmin,Tmax
      parameter(Tmin   = 1,   Tmax   = 3600.0)
      real     Bmin,Bmax
      parameter(Bmin   = 0.0, Bmax   = 30.0)
      real     RRRmin,RRRmax
      parameter(RRRmin = 1.5, RRRmax = 3000.0)
      real     brrmin,brrmax
      parameter(brrmin = 0.0, brrmax = 40.0e3)
      real     rho273,rhorrr
      parameter(rho273=4.50e-8,rhorrr=0.0)
      real     p1,p2,p3,p4,p5,p6,p7
      parameter(p1=2.6037e-17, p2=4.5661, p3=0.1974e11,
     &          p4=-1.2259, p5=10, p6=1, p7=0.333)
      real     a1,a2,a3,a4
      parameter(a1=4.67476e-5, a2=1.9686018, a3=4.62332e-5, 
     &          a4=0.9407898)

c * local variables   
      real     TT,BB,R,rhozero,rhoice,arg,rhoi,rhoi0,rho0,brr,magr

c *
      TT=T
      TT=max(TT,Tmin)
      TT=min(TT,Tmax)

      BB=B
      BB=max(BB,Bmin)
      BB=min(BB,Bmax)

      R =RRR
      R =max(R,RRRmin)
      R =min(R,RRRmax)

c * resistivity at absolute zero
      rhozero = rho273/(R-1.0)

c * fit for the resistivity at the ice temperature 
      rhoice  = rho273 + rhorrr/R

c * resistivity at the temperature t
      arg     = min((p5/TT)**p6,30.0)
      rhoi    = p1*TT**p2/(1.+p1*p3*TT**(p2+p4)*exp(-arg))
      rhoi0   = p7*rhoi*rhozero/(rhoi+rhozero)
      rho0    = rhozero+rhoi+rhoi0

c * product of field and residual resistivity ratio
      brr     = BB*rhoice/rho0
      brr     = max(brr,brrmin)
      brr     = min(brr,brrmax)
	
c * fit for the transverse magneto-resistance increase
      if(brr .gt. 1.) then
         magr = a1*brr**a2/(1.0+a3*brr**a4)
      else
         magr = 0.0
      endif

c * transverse magneto-resistance factor
      magrW   = magr + 1.0
	
c *
      return
      end