矢量化迭代循環

Question

夥計們，我想下面的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語句。下面是我的嘗試，但我認爲我只是矢量化了我需要矢量化的兩個之一。我已經包含了所有必要的函數和常量，以便代碼片段可以運行。該功能位於底部。我還包含了一個帶有正確答案的代碼段測試代碼。

任何幫助都非常感謝。

# 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 

Answer 1

到vectorise功能f最簡單的方法是使用Vectorize。默認情況下，它對所有參數矢量化爲f。在這種情況下，您只想爲3個參數中的2個進行矢量化處理，因此您可以通過vectorize.args指定。

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

（你也可以刪除裏面WetBulb的sapply）。這並不一定是最有效的方式來獲得矢量化（它基本上是一個mapply調用語法糖），但它肯定是最簡單的。