矢量化迭代循环

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调用语法糖），但它肯定是最简单的。