1

皆さん、次の Visual Basic コードを R に翻訳したいと思います。

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function WetBulb(T As Double, WDes As Double, PAtm As Double)
' Function to calculate wet-bulb temperature from dry-bulb
' and humidity ratio
Dim Wsat As Double
Dim TWBOld As Double
Dim WOld As Double
Dim TWBNew As Double
Dim TWB As Double
Dim WStar As Double
Dim W As Double
Dim slope As Double
Wsat = HumRatRH(T, RHMax, PAtm)
TWBOld = T
WOld = Wsat
TWBNew = TWBOld - 1
Do
    TWB = TWBNew
    WStar = HumRatRH(TWB, RHMax, PAtm)
    W = ((HfgRef - (CpWat - CpVap) * TWB) * WStar - CpAir * (T - TWB)) / (HfgRef + CpVap * T - CpWat * TWB)
    slope = (W - WOld) / (TWB - TWBOld)
    TWBNew = TWB - (W - WDes) / slope
    If Abs(W - WDes) < Abs(WOld - WDes) Then
        WOld = W
        TWBOld = TWB
    End If
Loop Until Abs((TWBNew - TWB) / TWB) < tolRel
WetBulb = TWB
End Function

私が抱えている問題は、ループにベクトルが含まれていることです。そのため、このループと if ステートメントを何らかの形でベクトル化する必要があります。以下は私の試みですが、ベクトル化する必要がある2つのうちの1つだけをベクトル化すると思います。スニペットが実行されるように、必要な関数と定数をすべて含めました。機能は一番下にあります。また、正解のスニペット テスト コードも含めました。

どんな助けでも大歓迎です。

# Constants independent of unit system
NMol   = 0.62198       # ratio of molecular weights, Mvap/MAir
RHMax  = 1             # maximum relative humidity, 1 or 100 (if percent)
tolRel = 0.000001      # relative error tolerance for iteration

# Constants for English Units
# Note: constants currently configured for PAtm in atmospheres
HfgRef = 1061          # heat of vaporization at 0C, Btu/hr.lbm.F
CpVap = 0.444          # specific heat of water vapor, Btu/hr.lbm.F
CpWat = 1              # specific heat of liquid water, Btu/hr.lbm.F
CpAir = 0.24           # specific heat of dry air, Btu/hr.lbm.F
RAir = 0.02521         # gas constant for air, (user pressure).ft3/lbm.R
kPaMult = 101.325      # multiplier to get kPascals from user pressure
TAbs = 459.67          # add to user temperature to get absolute temp
TKelMult = 0.555556    # multiplier to get Kelvin from user temp
TAmb = 70              # typical temperature in user units (initial value)
#####################################################################
SatPress <- function(TArg) {

# Define constants for vapor pressure correlations
C1  = -5674.5359
C2  = -0.51523058
C3  = -0.009677843
C4  = 0.00000062215701
C5  = 2.0747825E-09
C6  = -9.484024E-13
C7  = 4.1635019
C8  = -5800.2206
C9  = -5.516256
C10 = -0.048640239
C11 = 0.000041764768
C12 = -0.000000014452093
C13 = 6.5459673

T = (TArg + TAbs) * TKelMult
# Use different correlations for pressure over ice or water
    kPa.lo = exp(C1 / T + C2 + T * C3 + T * T * (C4 + T * (C5 + C6 * T)) + C7 * log(T))
    kPa.hi = exp(C8 / T + C9 + T * (C10 + T * (C11 + T * C12)) + C13 * log(T))
kPa = ifelse(T < 273.15, kPa.lo, kPa.hi)
SatPress = kPa / kPaMult
return(SatPress)

}
#####################################################################

HumRatRH = function(T,RH,PAtm) {
# function to calculate humidity ratio from temperature
# and relative humidity
pw = SatPress(T) * RH / RHMax
HumRatRH = NMol * pw / (PAtm - pw)
return(HumRatRH)
}
#####################################################################
WetBulb = function(T, WDes,PAtm) {
# Function to calculate wet-bulb temperature from dry-bulb
# and humidity ratio
Wsat = HumRatRH(T, RHMax, PAtm)
TWBOld = T
WOld = Wsat
TWBNew = TWBOld - 1
iterate.TWB = function(x) {
    repeat {
    TWB = TWBNew
    WStar = HumRatRH(TWB, RHMax, PAtm)
    W = ((HfgRef - (CpWat - CpVap) * TWB) * WStar - CpAir * (T - TWB)) / (HfgRef + CpVap * T - CpWat * TWB)
    slope = (W - WOld) / (TWB - TWBOld)
    TWBNew = TWB - (W - x) / slope
    TWBOld=ifelse(abs(W - x) < abs(WOld - x),TWB,TWBOld) # update TWBOld first
    WOld=ifelse(abs(W - x) < abs(WOld - x),w,WOld)       # then update WOld
    if (abs((TWBNew - TWB) / TWB) < tolRel) break()
    }
    return(TWB)
}
WetBulb = sapply(WDes, iterate.TWB)
return(WetBulb)
}

#####################################################################

temp = c(80,55,100)
w = c(0.011,0.009,0.016)
PAtm = 0.8187308
WetBulb(temp,w,PAtm)

# The correct answer:
# 62.95381538   51.3986312   74.02877887
4

1 に答える 1

2

関数をベクトル化する最も簡単な方法fは、 を使用することVectorizeです。デフォルトでは、fすべての引数に関してベクトル化します。この場合、3 つの引数のうち 2 つだけをベクトル化したいので、これを で指定しますvectorize.args。

WetBulb <- Vectorize(WetBulb, vectorize.args=c("T", "WDes"))

sapply(また、内部の を削除することもできますWetBulb。) これは必ずしもベクトル化を取得する最も効率的な方法ではありません (基本的にmapply呼び出しの構文糖衣です) が、確かに最も単純です。

于 2013-07-22T16:44:55.907 に答える