# based on SW_Density.m

SWdensity = function(Temp,uTemp,Sal,uSal) {
  # SW_Density    Density of seawater
  #=========================================================================
  # USAGE:  rho = SW_Density(Temp,uTemp,Sal,uSal)
  #
  # DESCRIPTION:
  #   Density of seawater at atmospheric pressure (0.1 MPa) using Eq. (8)
  #   given by [1] which best fit the data of [2] and [3]. The pure water
  #   density equation is a best fit to the data of [4]. 
  #   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:
  #   rho = density [kg/m^3]
  #
  #   Note: rho will have the same dimensions as Temp and Sal
  #
  # VALIDITY: 0 < Temp < 180 C; 0 < Sal < 160 g/kg;
  # 
  # ACCURACY: 0.1%
  # 
  # 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 license.
  # 
  # 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] Isdale, and Morris, Desalination, 10(4), 329, 1972.
  #   [3] Millero and Poisson, Deep-Sea Research, 28A (6), 625, 1981
  #   [4]  IAPWS release on the Thermodynamic properties of ordinary water substance, 1996. 
  #=========================================================================
  
  ## CHECK INPUT ARGUMENTS ----
  
  # CHECK THAT Sal AND Temp 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 > 160), 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(9.9992293295E+02,
        2.0341179217E-02,
        -6.1624591598E-03,
        2.2614664708E-05,
        -4.6570659168E-08)
  
  b = c(8.0200240891E+02,
        -2.0005183488E+00,
        1.6771024982E-02,
        -3.0600536746E-05,
        -1.6132224742E-05)
  
  rho_w = a[1] + a[2] * Temp + a[3] * Temp^2 + a[4] * Temp^3 + a[5] * Temp^4
  D_rho = b[1] * sal + b[2] * sal * Temp + b[3] * sal * Temp^2 + b[4] * sal * Temp^3 + b[5] * sal^2 * Temp^2
  rho   = rho_w + D_rho
  
  return(rho)
}