# http://www.srh.noaa.gov/epz/?n=wxcalc

# -----------------------------------------------------------------------------
# http://www.srh.noaa.gov/images/epz/wxcalc/tempConvert.pdf
# -----------------------------------------------------------------------------

Tf2Tc <- function(x) {
	return( (5/9) * (x - 32) )
}

Tf2Tk <- function(x) {
	return( (5/9) * (x - 32) + 273.15 )
}

Tc2Tf <- function(x) {
	return( (9/5) * (x + 32) )
}

Tc2Tk <- function(x) {
	return( x + 273.15 )
}

Tk2Tf <- function(x) {
	return( (9/5) * (x - 273.15) ) + 32
}

Tk2Tc <- function(x) {
	return( x - 273.15 )
}

Tf2Tr <- function(x) { # R is Rankine
	return( x + 459.69 )
}

Tr2Tf <- function(x) { # R is Rankine
	return( x - 459.69 )
}

Tc2Tr <- function(x) { # R is Rankine
	return( (9/5) * (x + 32 ) + 459.69 )
}

Tr2Tc <- function(x) { # R is Rankine
	return( (5/9) * (x - 459.69) - 32 )
}

# -----------------------------------------------------------------------------
#http://www.srh.noaa.gov/images/epz/wxcalc/pressureConversion.pdf
# -----------------------------------------------------------------------------

Pinhg2Pmmhg <- function(x) {
	return( 25.4 * x )
}

Pmmhg2Pinhg <- function(x) {
	# return( x / 25.4 )
	return( 0.03937008 * x )
}

Pinhg2Pmb <- function(x) {
	return( 33.8639 * x )
}

Pmb2Pinhg <- function(x) {
	# return( x / 33.8639 )
	return( 0.0295300 * x )
}

Pinhg2Pkpa <- function(x) {
	return( 33.8639 * x  / 10)
}

Pkpa2Pinhg <- function(x) {
	return( 0.295300 * x )
}

Pinhg2Ppsi <- function(x) {
	return( 0.491154 * x )
}

Ppsi2Pinhg <- function(x) {
	return( 0.491154 * x )
}

Pmmhg2Pmb <- function(x) {
	return( 1.333224 * x )
}

Pmb2Pmmhg <- function(x) {
	return( 0.750062 * x )
}

Pkpa2Pmmhg <- function(x) {
	return( 7.50062 * x )
}

Pmmhg2Ppsi <- function(x) {
	return( 0.0193368 * x )
}

Ppsi2Pmmhg <- function(x) {
	return( 51.7149 * x )
}

Pmb2Pkpa <- function(x) {
	return( x / 10 )
}

Pkpa2Pmb <- function(x) {
	return( 10 * x )
}

Ppsi2Pmb <- function(x) {
	return( 68.9476 * x )
}

Pkpa2Psi <- function(x) {
	return( 0.145038 * x )
}

Ppsi2Pkpa <- function(x) {
	return( 6.89476 * x )
}

Pmmhg2Patm <- function(x) {
	return( 760 * x )
}

Patm2Pmmgh <- function(x) {
	return( 0.0013157894 * x )
}

Pinhg2Patm <- function(x) {
	return( 0.033421 * x )
}

Patm2Pinhg <- function(x) {
	return( 29.9213 * x )
}

Ppsi2Patm <- function(x) {
	return( 0.0680457 * x )
}

Patm2Psi <- function(x) {
	return( 14.6960 * x )
}

Pmb2Patm <- function(x) {
	return( 0.009869233 * x )
}

Patm2Pmb <- function(x) {
	return( 1013.250 * x )
}

Pkpa2Patm <- function(x) {
	return( 0.009869233 * x )
}

Patm2Pkpa <- function(x) {
	return( 101.3250 * x )
}

# -----------------------------------------------------------------------------
# http://www.srh.noaa.gov/images/epz/wxcalc/windConversion.pdf
# -----------------------------------------------------------------------------

# miles per hour to knots
Wmph2Wknts <- function(x) {
	return( 0.8689762 * x )
}

# knots to miles per hour
Wknts2Wmph <- function(x) {
	return( 1.1507794 * x )
}

# miles per hour to meters per second
Wmph2Wmps <- function(x) {
	return( 0.44704 * x )
}

# meters per second to miles per hour
Wmps2Wmph <- function(x) {
	return( 0.44704 * x )
}

# miles per hour to feet per second
Wmph2Wftps <- function(x) {
	return( 1.46667 * x )
}

# feet per second to miles per hour
Wftps2Wmph <- function(x) {
	return( 0.681818 * x )
}

# miles per hour to kilometers per hour
Wmph2Wkmph <- function(x) {
	return( 1.609344 * x )
}

# kilometers per hour to miles per hour
Wkmph2Wmph <- function(x) {
	return( 0.621371 *  x )
}

# knots to meters per second
Wknts2Wmps <- function(x) {
	return( 0.5144444 * x )
}

# meters per second to knots
Wmps2Wknts <- function(x) {
	return( 1.9438445 * x )
}

# knots to feet per second
Wknts2Wftps  <- function(x) {
	return( 1.6878099 * x )
}

# feet per second to knots
Wftps2Wknts <- function(x) {
	return( 0.5924838 * x )
}

# knots to kilometers per hour
# ??? mislabeled in source material ???
Wknts2Wkmph <- function(x) {
	return( 1.852 * x )
}

# kilometers per hour to knots
# ??? mislabeled in source material ???
Wkmph2Wknts <- function(x) {
	return( 1.852 * x )
}

# feet per second to meters per second
Wftps2Wmps <- function(x) {
	return( 0.3048 * x )
}

# meters per second to feet per second
Wmps2Wftps <- function(x) {
	return( 13.28084 * x )
}

# meters per second to kilometers per hour
Wmps2Wkmph <- function(x) {
	return( 3.6 * x )
}

# kilometers per hour to meters per second
Wkmph2Wmps <- function(x) {
	return( 1.852 * x )
}

# feet per second to kilometers per hour
Wftps2Wkmph <- function(x) {
	return( 1.09728 * x )
}

# kilometers per hour to feet per second
Wkmph2Wftps <- function(x) {
	return( 1.852 * x )
}

# -----------------------------------------------------------------------------
# http://www.srh.noaa.gov/images/epz/wxcalc/vaporPressure.pdf
# -----------------------------------------------------------------------------

# actual vapor pressure
getVP <- function(T) { # T = temperature (C)
	return ( 6.11*10^((7.5*T)/(237.7+T)) )
}

# relative humidity in millibars
getRH <- function(T,Td) {
	#T = temperature (C)
	#Td = dew point temperature (C)
	e = getVP(Td) # saturation vapor pressure
	es = getVP(T) # actual vapor pressure
	return ( 100 * e / es ) # millibars
}

# -----------------------------------------------------------------------------
#http://www.srh.noaa.gov/images/epz/wxcalc/rhTdFromWetBulb.pdf
# -----------------------------------------------------------------------------

getVPwb <- function(T) { # T = temperature (C)
	return ( 6.112*exp((17.67*T)/(243.5+T)) )
	# http://www.fao.org/docrep/X0490E/x0490e07.htm
	# return ( 6.108*exp((17.27*T)/(237.3+T)) )
}

# get relative humidity from wet bulb temperature
getRHwb <- function(T,Tw,P) {
	# T = temperature (C)
	# Tw = wet bulb temperature (C)
	# P = station pressure (mb)
	
	es = getVPwb(T) # saturation vapor pressure
	ew = getVPwb(Tw) # wet-bulb vapor pressure
	e = ew - P * (T - Tw) * 0.00066 * (1 + (0.00115 * Tw) )
	return ( 100 * e / es )
}
	
# get dew point temperature from wet bulb temperature
getTDwb <- function(T,Tw,P) {
	# T = temperature (C)
	# Tw = wet bulb temperature (C)
	# P = station pressure (mb)
	
	es = getVPSwb(T) # saturation vapor pressure
	ew = getVPWwb(Tw) # wet-bulb vapor pressure
	e = ew - P * (T - Tw) * 0.00066 * (1 + (0.00115 * Tw) )
	return ( (243.5 * log(e/6.112)) / (17.67 - ln(e/6.112)) )
}

# -----------------------------------------------------------------------------
#http://www.srh.noaa.gov/images/epz/wxcalc/densityAltitude.pdf
# -----------------------------------------------------------------------------

getDensityAltitude <- function(T,Td,P) {
	# T = temperature (C)
	# Td = dew point temperature (C)
	# P = station pressure (mb)	
	
	T = Tc2Tk(T) # convert to Kelvin
	e = getVP(Td) # vaporization temperature
	Tv = T / (1 - ( (e/P) * (1-0.622) ) ) # virtual temperature
	Tv = Tk2Tr(Tv) # convert T from Kelvin to Rankine
	P = Pmb2Pinhg(P) # convert P from mb to in Hg
	return( 145366 * (1 - ((17.326*Pinhg)/Tv)^0.235) )
}

# -----------------------------------------------------------------------------
# http://www.srh.noaa.gov/images/epz/wxcalc/heatIndex.pdf
# -----------------------------------------------------------------------------

getHeatIndex <- function(T,Td) {
	# T = temperature (C)
	# Td = dew point temperature (C)
	rh = getRH(T,Td)
	Hidx = ( -42.379 
		+ (2.04901523 * T)
		+ (10.14333127 * rh)
		- (0.22475541 * T * rh)
		- (6.83783 * 10^(-3) * T^2)
		- (5.481717 * 10^(-2) * rh^2)
		+ (1.22874 * 10^(-3) * T^2 * rh)
		+ (8.5282 * 10^(-4) * T * rh^2)
		- (1.99 * 10^(0-6) * T^2 * rh^2) )
		
	return (Hidx)
}

# -----------------------------------------------------------------------------
# http://www.srh.noaa.gov/images/epz/wxcalc/mixingRatio.pdf
# -----------------------------------------------------------------------------

getMixingRatio <- function(T,P) {
	# T = temperature (C)
	# P = station pressure (mb)
	e = getVP(T)
	return( 621.97 * ( e / (P-e) ) )
}

# -----------------------------------------------------------------------------
# http://www.srh.noaa.gov/images/epz/wxcalc/pressureAltitude.pdf
# -----------------------------------------------------------------------------

# pressure altitude in meters
getPressureAltitude <- function(P) {
	# P = pressure (mb)
	return ( 0.3048 * 145366.45 * (1 - (P/1013.25)^0.190284) )
}

# -----------------------------------------------------------------------------
# http://www.srh.noaa.gov/images/epz/wxcalc/speedOfSound.pdf
# -----------------------------------------------------------------------------

# speed of sound in knots
getSpeedOfSound <- function(T) { # temperature (C)
	T = Tc2Tk(T) # convert to Kelvin
	return ( 643.855 * (T/273.15)^0.5 )
}

# -----------------------------------------------------------------------------
# http://www.srh.noaa.gov/images/epz/wxcalc/stationPressure.pdf
# -----------------------------------------------------------------------------

getStationPressure <- function(H,P) {
	# H = elevation in meters
	# P = pressure in mb
	P = Pmb2Pinhg(P) # convert to inches Hg
	Pstn = ( P * ( (288-0.0065*H)/288 )^5.2561 )
	Pstn = Pinhg2Pmb(Pstn) # convert back to mb
	return(Pstn)
}

# -----------------------------------------------------------------------------
# http://www.srh.noaa.gov/images/epz/wxcalc/virtualTemperature.pdf
# -----------------------------------------------------------------------------

getVirtualTemperature <- function(T,Td,P) {
	# T = temperature (C)
	# Td = dew point temperature (C)
	# P = station pressure (mb)
	
	Tv = (T + 273.15) / (1 - 0.379*( ( 6.11 * 10^( (7.5*Td) / (237.7+Td) ) ) / P ) )
	return (Tv)	
}

# -----------------------------------------------------------------------------
# http://www.srh.noaa.gov/images/epz/wxcalc/wetBulbTdFromRh.pdf
# -----------------------------------------------------------------------------

getDewpointTemperature <- function(T,RH,P) {
	# T = temperature (C)
	# RH = relative humidity
	# P = station pressure (mb)

	es = getVP(T) # saturated vapor pressure
	return ( (237.7 * log ((es * RH)/611)) / (7.5 - log ((es * RH)/611)) )
}

# use Skew-T diagram
# getWetbulbTemperature

# -----------------------------------------------------------------------------
# http://www.srh.noaa.gov/images/epz/wxcalc/windChill.pdf
# -----------------------------------------------------------------------------


getWindChill <- function(T,W) {
	# T = temperature (C)
	# W = wind speed (mph)
	
	wc = ( 35.74
		+ (0.6215 * T)
		- (35.75 * W^0.16)
		+ (0.4275 * T * W^0.16) )
	return (wc)
}

# watts per meter squared (W/m^2)
getWindChillwps <- function(T,W) {
	# T = temperature (C)
	# W = wind speed (meters per second)
	
	wc = ( (33-T) * 
		(12.1452
		+ 11.6222 * sqrt(W)
		- 1.16222 * W) )
	return (wc)
}