Skip to content

Instantly share code, notes, and snippets.

@nozma
Created December 8, 2012 14:10
Show Gist options
  • Save nozma/4240385 to your computer and use it in GitHub Desktop.
Save nozma/4240385 to your computer and use it in GitHub Desktop.
湿度計算セット
## --- 気温から飽和水蒸気圧を求める関数色々 ---
GofGra <- function(t){
## 入力:気温(℃)
## 出力:飽和水蒸気圧(hPa)
10^
(10.79574 * (1 - 273.16/(t + 273.15)) -
5.02800 * log((t + 273.15)/273.16, 10) +
1.50475 * 10^{-4} * {1-10^(-8.2969 * ((t + 273.15)/273.16 - 1))} +
0.42873 * 10^{-3} * {10^(4.76955*(1 - 273.16/(t + 273.15))) - 1} +
0.78614)
}
Okada <- function(t){
## 入力:気温(℃)
## 出力:飽和水蒸気圧(hPa)
exp(1.809378 +
0.07266115 * t +
(-3.003879 * 10^-4) * t^2 +
(1.181765 * 10^-6) * t^3 +
(-3.863083 * 10^-9) * t^4)
}
murray <- function(t){
## 入力:気温(℃)
## 出力:飽和水蒸気圧(hPa)
6.1078 * exp(17.2693882 * t /(t + 237.3))
}
tetens <- function(t){
## 入力:気温(℃)
## 出力:飽和水蒸気圧(hPa)
6.11 * 10^(7.5*t/(t + 237.3))
}
## --- 乾湿計から水蒸気圧を求める関数色々 ---
sprung <- function(td, tw){
## 入力:td(乾球温度、℃) tw(湿球温度、℃)
## 出力:水蒸気圧(hPa)
GofGra(tw) - 0.000662 * 1013.25 * (td - tw)
}
vp1 <- function(td, tw){ ## 風速1~1.5m/sの場合に使用
## 入力:td(乾球温度、℃) tw(湿球温度、℃)
## 出力:水蒸気圧(hPa)
GofGra(tw) - 0.0008 * 1013.25 * (td - tw) * (1 + (tw / 610))
}
vp2 <- function(td, tw){ ## 風速0.5m/s以下の場合に使用
## 入力:td(乾球温度、℃) tw(湿球温度、℃)
## 出力:水蒸気圧(hPa)
GofGra(tw) - 0.0012 * 1013.25 * (td - tw) * (1 + (tw / 610))
}
vp3 <- function(td, tw){ ## 通風なしの場合
## 入力:td(乾球温度、℃) tw(湿球温度、℃)
## 出力:水蒸気圧(hPa)
GofGra(tw) - 0.0008 * 1013.25 * (td - tw)
}
abs.hum <- function(t, U){
## 入力:t…気温(℃)、U…相対湿度(%)
## 出力:絶対湿度(g/m^3)
2.166740 * 10^2 * U * GofGra(t)/(100 * (t + 273.15))
}
## --- 各種湿度諸量を求める関数 ---
mix.r <- function(t, U, P = 1013.25){
## 入力:t…気温(℃)、U…相対湿度(%)
## :P…大気圧(hPa、省略時標準大気圧1013.25)
## 出力:混合比(kg/kg)
0.622 * (U/100 * GofGra(t)) / (P - U/100 * GofGra(t))
}
sup.hum <- function(t, U, P = 1013.25){
## 入力:t…気温(℃)、U…相対湿度(%)
## :P…大気圧(hPa、省略時標準大気圧1013.25)
## 出力:比湿(kg/kg)
0.622 * (U/100 * GofGra(t)) / (P - 0.378 * U/100 * GofGra(t))
}
mol.frac <- function(t, U, P = 1013.25){
## 入力:t…気温(℃)、U…相対湿度(%)
## :P…大気圧(hPa、省略時標準大気圧1013.25)
## 出力:モル分率(mol/mol)
U/100 * GofGra(t)/P
}
vpd <- function(t, U){
## 入力:t…気温(℃)、U…相対湿度(%)
## 出力:飽差(hPa)
GofGra(t) * (1 - U/100)
}
w.dif <- function(t, U){
## 入力:t…気温(℃)、U…相対湿度(%)
## 出力:飽差?(g/m^3)
abs.hum(t, 100) - abs.hum(t, U)
}
dew <- function(t, U){
## 入力:t…気温(℃)、U…相対湿度(%)
## 出力:露点(℃)
-(log(GofGra(t)*U/100/6.1078) * 237.3) /
(log(GofGra(t)*U/100/6.1078) - 17.2693882)
}
dew.dep <- function(t, U){
## 入力:t…気温(℃)、U…相対湿度(%)
## 出力:湿数(℃)
t - dew(t, U)
}
dewdep2rh <- function(t, dd){
## 入力:t…気温(℃)、dd…湿数(℃)
## 出力:相対湿度(%)
GofGra(t-dd)/GofGra(t) * 100
}
## --- その他の湿度諸量 ---
DI1 <- function(t, U){
## 入力:t…気温(℃)、U…相対湿度(%)
## 出力:不快指数
0.81 * t + U/100 * (0.99 * t - 14.3) + 46.3
}
DI2 <- function(t, td){
## 入力:t…気温(℃)、td…露点(℃)
## 出力:不快指数
0.99 * t + 0.36 * td + 41.5
}
DI3 <- function(t, tw){
## 入力:t…気温(℃)、tw…通風なしの湿球温度(℃)
## 出力:不快指数
0.72 * t * (t * tw) + 40.6
}
EU <- function(U, EU1 = NULL){
## 入力:U…相対湿度(%) ベクトル
## :EU1…前日実効湿度(%) 入力時は相対湿度は当日分のみで可
## 出力:実効湿度(%)
if(is.null(EU1)){ # 前日の実効湿度なし
(1 - 0.7) * (U[1] + sum(0.7^(1:(length(U)-1)) * U[-1]))
} else { # 前日の実効湿度あり
(1 - 0.7) * U[1] + 0.7 * EU1
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment