# based on SW_Viscosity.m

SWviscosity = function(Temp,uTemp,Sal,uSal) {
  # SW_Viscosity    Dynamic viscosity of seawater
  #=========================================================================
  # USAGE:  mu = SW_Viscosity(Temp,uTemp,Sal,uSal)
  #
  # DESCRIPTION:
  #   Dynamic viscosity of seawater at atmospheric pressure (0.1 MPa) using
  #   Eq. (22) given in [1] which best fit the data of [2], [3] and [4].
  #   The pure water viscosity equation is a best fit to the data of [5].
  #   Values at temperature higher than the normal boiling temperature
  #   are calculated at the saturation pressure.
  #
  # INPUT:
  #   Temp  = temperature
  #   uTemp = temperature unit
  #        'C'  : [degree Celsius] (ITS-90)
  #        'K'  : [Kelvin]
  #        'F'  : [degree Fahrenheit]
  #        'R'  : [Rankine]
  #   Sal  = salinity
  #   uSal = salinity unit
  #        'ppt': [g/kg]  (reference-composition salinity)
  #        'ppm': [mg/kg] (in parts per million)
  #        'w'  : [kg/kg] (mass fraction)
  #        '%'  : [kg/kg] (in parts per hundred)
  #
  #   Note: Temp and Sal must have the same dimensions
  #
  # OUTPUT:
  #   mu = dynamic viscosity [kg/m-s]
  #
  #   Note: mu will have the same dimensions as Temp and Sal
  #
  # VALIDITY: 0 < Temp < 180 C and 0 < Sal < 150 g/kg;
  #
  # ACCURACY: 1.5%
  #
  # REVISION HISTORY:
  #   2009-12-18: Mostafa H. Sharqawy (mhamed@mit.edu), MIT
  #               - Initial version
  #   2012-06-06: Karan H. Mistry (mistry@mit.edu), MIT
  #               - Allow Temp,Sal input in various units
  #               - Allow Temp,Sal to be matrices of any size
  #   2014-12-15: Daniel Neumann (daniel.neumann@hzg.de), Helmholtz-Zentrum Geesthacht
  #               - converted MATLAB code to R
  #
  # DISCLAIMER:
  #   This software is provided "as is" without warranty of any kind.
  #   See the file sw_copy.m for conditions of use and licence.
  #
  # REFERENCES:
  #   [1] M. H. Sharqawy, J. H. Lienhard V, and S. M. Zubair, Desalination
  #       and Water Treatment, 16, 354-380, 2010. (http://web.mit.edu/seawater/)
  #   [2] B. M. Fabuss, A. Korosi, and D. F. Othmer, J., Chem. Eng. Data 14(2), 192, 1969.
  #   [3] J. D. Isdale, C. M. Spence, and J. S. Tudhope, Desalination, 10(4), 319 - 328, 1972
  #   [4] F. J. Millero, The Sea, Vol. 5, 3  80, John Wiley, New York, 1974
  #   [5] IAPWS release on the viscosity of ordinary water substance 2008
  #=========================================================================
  
  ## CHECK INPUT ARGUMENTS ----
  
  # CHECK THAT S&T HAVE SAME SHAPE
  if (length(Sal) != length(Temp)) {
    stop('check_stp: Sal & Temp must have same dimensions');
  }
  
  # CONVERT TEMPERATURE INPUT TO °C
  if (tolower(uTemp) == 'k') {
    Temp = Temp - 273.15
  } else if (tolower(uTemp) == 'f') {
    Temp = 5/9*(Temp-32)
  } else if (tolower(uTemp) == 'r') {
    Temp = 5/9*(Temp-491.67)
  } else if (tolower(uTemp) != 'c') {
    stop("Not a recognized temperature unit. Please use 'C', 'K', 'F', or 'R'");
  }
  
  # CONVERT SALINITY TO PPT
  if ( tolower(uSal)%in%c('ppm', 'mg/kg') ) {
    Sal = Sal/1000
  } else if ( tolower(uSal)%in%c('w', 'kg/kg') ) {
    Sal = Sal*1000;
  } else if ( tolower(uSal) == '%' ) {
    Sal = Sal*10;
  } else if ( !(tolower(uSal)%in%c('ppt','g/kg')) ) {
    stop("Not a recognized salinity unit. Please use 'ppt', 'ppm', 'w', or '%'")
  } 
  
  # CHECK THAT Sal & Temp ARE WITHIN THE FUNCTION RANGE
  if ( sum((Temp < 0) + (Temp > 180), na.rm = TRUE) > 0) {
    warning('At least one temperature value is out of range for density function 0 < Temp < 180 C. 
            Temperature is set to 0 or 180 C.')
    Temp[Temp > 180] = 180
    Temp[Temp < 0] = 0
  }
  
  if ( sum((Sal < 0) + (Sal > 150), na.rm = TRUE) > 0 ) {
    warning('At least one salinity value is out of range for density function 0 < Sal < 160 g/kg')
  }
  
  ## BEGIN ----
  sal = Sal/1000
  
  a = c(1.5700386464E-01,
        6.4992620050E+01,
        -9.1296496657E+01,
        4.2844324477E-05,
        1.5409136040E+00,
        1.9981117208E-02,
        -9.5203865864E-05,
        7.9739318223E+00,
        -7.5614568881E-02,
        4.7237011074E-04)
  
  mu_w = a[4] + 1 / ( a[1] * (Temp + a[2] )^2 + a[3] )
  
  A  = a[5] + a[6] * T + a[7] * Temp^2
  B  = a[8] + a[9] * T + a[10]* Temp^2
  mu = mu_w * ( 1 + A * sal + B * sal^2 )
  
}