# 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)
}

# -----------------------------------------------------------------------------
# http://www.fao.org/docrep/X0490E/x0490e07.htm
# -----------------------------------------------------------------------------

# get daily radiaion (top of atmosphere)
# http://www.fao.org/docrep/X0490E/x0490e07.htm
getRtoa <- function(lat,jday) {
	# jday is julian day of the year
	# Rn = net radiance per day # MJ / (day * m^2)
	
	# Gsc solar constant = 0.0820 MJ m-2 min-1,
	Gsc = 0.0820 # MJ/(min*m^2)
	# dr inverse relative distance Earth-Sun
	dr = 1 + 0.033*cos(2*pi*jday/365)
	# theta = latitude (radians)
	theta = (pi/180) * lat
	# delta = solar decimation (radians)
	delta = 0.409*sin((2*pi*jday/365 - 1.39))
	# w = sunset hour angle (radians)
	w = acos(-tan(theta) * tan(delta))
	Rn = (24 * 60 / pi) * Gsc * dr * ( 
			w*sin(theta)*sin(delta) 
			+ cos(theta)*cos(delta)*sin(w) )
	return(Rtoa)
}

getDaylightHours(lat,jday) {
	# theta = latitude (radians)
	theta = (pi/180) * lat
	# delta = solar decimation (radians)
	delta = 0.408*sin((2*pi*jday/365 - 1.39))
	# w = sunset hour ange (radians)
	w = acos(-tan(theta) * tan(delta))
	N = (24/pi)*w
	return(N)
}

# Rs solar or shortwave radiation [MJ m-2 day-1],
getRs(lat,jday){

	# n actual duration of sunshine [hour],
	# The actual duration of sunshine, n, is recorded with a Campbell Stokes sunshine recorder
	# n = ????
	
	# N maximum possible duration of sunshine or daylight hours [hour],
	N = getDayLightHours(lat,jday)
	# n/N relative sunshine duration [-],

	# RB:20120708: replace n with something ... based on cloudiness?
	n = N
	
	# Ra extraterrestrial radiation [MJ m-2 day-1],
	Ra = getRtoaDaily(lat,jday)

	# 'as' regression constant, expressing the fraction of extraterrestrial radiation 
	# reaching the earth on overcast days (n = 0),
	# as+bs fraction of extraterrestrial radiation reaching the earth on clear days (n = N).}
	
	# Where no actual solar radiation data are available and no calibration has been carried out 
	# for improved as and bs parameters, the values as = 0.25 and bs = 0.50 are recommended.
	as = 0.25
	bs = 0.50 
	
	Rs = (as+bs*(n/N))*Ra
	return(Rs)
}

#Rn = Rns - Rnl

# Clear-sky solar radiation (Rso) When calibrated values for as and bs are not available:
getRso <- function(Rtoa,z) {
	Rso = (0.75 + 2*10^-5*z)*Ra 
	return(Rso)
}

# Net solar or net shortwave radiation (Rns)
# The net shortwave radiation resulting from the balance 
# between incoming and reflected solar radiation is given by

getRns <- function(Rs) {
	# Rs the incoming solar radiation [MJ m-2 day-1].
	# albedo or canopy reflection coefficient, 
	# which is 0.23 for the hypothetical grass reference crop [dimensionless],
	Rns = (1-albedo)*Rs
	return(Rns)
}

# Net longwave radiation (Rnl)
getRnl <- function(T,Td,z,lat,jday) {
	Rs = getRs(lat,jday)
	Rso = getRso(getRtoa(lat,jday),z)
	e = getVPwb(Td)
	Tk = Tc2Tk(T)
	sigma_sb = 4.903
	Rnl = sigma_sb*Tk*(0.34-0.14*sqrt(e))*(1.35*Rs/Rso - 0.35)
	return(Rnl)
}

# Net radiation (Rn)
# The net radiation (Rn) is the difference between 
# the incoming net shortwave radiation (Rns) and 
# the outgoing net longwave radiation (Rnl):
getRn <- function(Rns,Rnl) {
	Rn = Rns - Rnl
	return(Rn)
}

# soil heat flux
getSoilHeatFlux <- function(Tnow,Tprev,Tnext) {

	# G soil heat flux [MJ m-2 day-1],
	# cs soil heat capacity [MJ m-3 °C-1],
	# Ti air temperature at time i [°C],
	# Ti-1 air temperature at time i-1 [°C],
	# delta_t length of time interval [day],
	# delta_z effective soil depth [m].
	
	# As the soil temperature lags air temperature, 
	# the average temperature for a period should be considered 
	# when assessing me daily soil heat flux, i.e., D_t should exceed one day. 
	# The depth of penetration of the temperature wave is determined 
	# by the length of the time interval. The effective soil depth, D_z, 
	# is only 0.10-0.20 m for a time interval of one or a few days 
	# but might be 2 m or more for monthly periods. 
	# The soil heat capacity is related to its mineral composition and water content.
	
	# G = cs * (Tj - Ti) * delta_z / delta_t
	
	# For day and ten-day periods:
	# G = 0
	
	# For monthly periods, where the month before and after are known
	# G = 0.07 (Tmonth_(i+1) - Tmonth_(i-1)) 
	
	# For monthly periods, where the month before is known but not month after
	# G = 0.14 (Tmonth_(i) - Tmonth_(i-1)) 
	
	if (Tnext == 999) {G = 0.014 * (Tnow - Tprev)}
	else {G = 0.07 * (Tnext - Tprev)}

	return(G)
}

# http://www.fao.org/docrep/X0490E/x0490e07.htm
# evapotranspiration process
getPMET <- function (T,Td,P,W,Rtoa) {

	# Atmospheric pressure (P)
	# given as P in the calling variables
	
	# Latent heat of vaporization (lambda)
	# As lambda_v varies only slightly over normal temperature ranges 
	# a single value of 2.45 MJ kg-1 is taken in the simplification 
	# of the FAO Penman-Monteith equation
	lambda_v = 2.45

	# Psychrometric constant (gamma)
	gamma = 0.0016286 * P / lambda_v # psychometric constant (kPA /K)

	# Air temperature
	# Tmean given as T in calling variables
	
	# Mean saturation vapour pressure (es)
	# calculated by getVPwb from T
	
	# Actual vapour pressure (ea) derived from dewpoint temperature
	# calculated by getVPwb from Td
	
	# delta_e = (es -e) = (1-rh)*es
	delta_e = getVPwb(T) - getVPwb(Td)
	
	# Net Radiation 
	# net incoming (bottom of atmos) - net outgoing
	Rns = getRns(getRs(lat,jday))
	Rnl = getRnl(T,Td,z,lat,jday)
	Rn = Rns - Rnl

	# Wind Speed (m/s)
	# given as W in calling variables

	PMET = ( 0.408*delta*(Rn - G) + gamma*900/(Tk)*W*(delta_e) )/(delta + gamma*(1+0.34*W))

	# http://www.fao.org/docrep/X0490E/x0490e07.htm
	#m = 4098*(0.6108*exp((17.27*T)/(T+237.3))) / (T+237.3)^2 # (kPa/C)
	#Emass = ( m*Rn + gamma*6.43*(1+0.536*U2)*delta_e ) / (lambda_v * (m + gamma) ) 
	#return(Emass)
	
	return(PMET)
}

get ET_alt <- function(Tmax,Tmin,lat,jday)
	Rtoa = getRtoa(lat,jday)
	ET_alt = 0.0023((Tmax+Tmin)/2 + 17.8)((Tmax - Tmin)^0.5)*Rtoa 
	return(ET_alt)
}