MATLAB style for loops in R & loop fails to fill each column in R

Question

I am in the process of converting the Simplified Psychrometric Chart MATLAB code from http://www.ohio.edu/mechanical/thermo/Applied/Chapt.7_11/Psychro_chart/psychro.html into R, but I am stuck on a couple of operations.

1) I am trying to implement the following MATLAB code in R:

# MATLAB code
# both pg and patm are defined earlier in the program

for phi = 0.1:0.1:0.4, % phi = relative humidity 10% - 40%
   w = 622*phi*pg./(patm-phi*pg);
   plot(t,w)
end

This is my attempt at translating the above code into R:

patm <- 101.325

phi <- as.matrix(seq(0.1, 0.4, 0.1))

pg <- matrix(c(0.61165, 0.65709, 0.70599, 0.75808, 0.81355, 0.87258, 0.93536, 1.00210, 1.07300, 1.14830, 
1.22820), nrow = 11, ncol = 1) 

w <- data.frame(matrix(nrow = nrow(pg), ncol = length(phi)))

for (i in ncol(w))
{
w[i] <- 622 * phi[[i]] * pg/(patm - phi[[i]] * pg)
}

Instead of filling each column, only the last column has the result:

 # X1   X2  X3  X4
 # 1    NA  NA  NA  1.505520
 # 2    NA  NA  NA  1.617658
 # 3    NA  NA  NA  1.738379
 # 4    NA  NA  NA  1.867026
 # 5    NA  NA  NA  2.004080
 # 6    NA  NA  NA  2.149996
 # 7    NA  NA  NA  2.305256
 # 8    NA  NA  NA  2.470394
 # 9    NA  NA  NA  2.645922
 # 10   NA  NA  NA  2.832450
 # 11   NA  NA  NA  3.030496

 dput(w)
 structure(list(X1 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
 NA), X2 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), X3 = c(NA, 
 NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), X4 = structure(c(1.50552046025963, 
 1.61765774182314, 1.73837871399276, 1.8670263620807, 2.00408001696641, 
 2.14999576929037, 2.30525601486727, 2.47039440308445, 2.64592183222691, 
 2.83245044300499, 3.03049575082621), .Dim = c(11L, 1L))), .Names = c("X1", 
 "X2", "X3", "X4"), row.names = c(NA, -11L), class = "data.frame")

How should the R code be revised to get the desired result?

The following code is what I would like to obtain (desired result).

#          X1                  X2                   X3              X4
# 1      0.3756981  0.7518503   1.128458    1.505520
# 2      0.4036271  0.8077785   1.212455    1.617658
# 3      0.4336856  0.8679764   1.302874    1.738379
# 4      0.4657082  0.9321142   1.399220    1.867026
# 5      0.4998122  1.0004283   1.501850    2.004080
# 6      0.5361091  1.0731432   1.611105    2.149996
# 7      0.5747165  1.1504960   1.727342    2.305256
# 8      0.6157644  1.2327491   1.850958    2.470394
# 9      0.6593768  1.3201530   1.982333    2.645922
# 10    0.7057024   1.4130080   2.121922    2.832450
# 11    0.7548656   1.5115656   2.270107    3.030496

2) Can the following MATLAB code be done in R using base plot or ggplot2?

# MATLAB code
# t1, tv,, and wg1 are defined previously in the program

for i = 1:7,
   plot([t1(i),tv0(i)], [wg1(i),0],'g-')
end

This is how the MATLAB plot call could look in R without the for loop:

plot(x1, y1, x2, y2)

I want to thank you in advance.

Answer 1

I wrote the following code based on the comments, which uses ggplot2. It works, but I was wondering if there is a way to simplify the code for use in ggplot2?

Thank you.

patm <- 101.325 # standard atmosphere (kPa)
rair <- 0.287 # gas constant of air (kJ/kg.K)

tpg <- fread("water_sat_t_p_1.txt") # Saturation Properties for Water — 1°C
Temperature Increments
tpg <- setDF(tpg)
tpg <- tpg[1:51, 1:2]
names(tpg)[1] <- "Temperature (°C)"
tpg[, 2] <- tpg[,2] * 1e3
names(tpg)[1] <- "Pressure (kPa)"

t <- data.frame(tpg[, 1]) # Temperature (°C)
pg <- data.frame(tpg[, 2]) # saturation vapor pressure (kPa)
wg <- 622 * pg/(patm - pg) # saturation specific humidity

phi <- seq(0.1, 0.4, 0.1) # phi = relative humidity 10% - 40%
w <- data.frame(matrix(nrow = nrow(pg), ncol = length(phi)))

for (i in 1:ncol(w))
{
w[i] <- 622 * phi[[i]] * pg/(patm - phi[[i]] * pg)
}

phi1 <- seq(0.6, 0.8, 0.2) # phi = 60%, 80%
w1 <- data.frame(matrix(nrow = nrow(pg), ncol = length(phi1)))

for (i in 1:ncol(w1))
{
w1[i] <- 622 * phi1[[i]] * pg/(patm - phi1[[i]] * pg)
}


# specific volume and enthalpy/wet-bulb Temperature
tpg1 <- fread("water_sat_t_p_5.txt") # Saturation Properties for Water — 5°C 
Temperature Increments
tpg1 <- setDF(tpg1)
tpg1 <- tpg1[1:7, 1:2]
names(tpg1)[1] <- "Temperature (°C)"
tpg1[, 2] <- tpg1[,2] * 1e3
names(tpg1)[1] <- "Pressure (kPa)"

t1 <- data.frame(tpg1[, 1]) # Temperature (°C)
pg1 <- data.frame(tpg1[, 2]) # saturation vapor pressure (kPa)

wg1 <- data.frame(622 * pg1/(patm - pg1)) # saturation specific humidity
# specific volume of dry air (cubic m/kg dry air) (green)
vol <- rair * (t1 + 273.15)/(patm - pg1) # specific vol at saturation
tv0 <- patm * vol/rair - 273.15 # air temperature at zero humidity

zerocol <- data.frame(rep_len(0, nrow(t1)))


# wet bulb temperature (also enthalpy) lines (red)
h <- t1 + 2.5 * wg1 # enthalpy (kJ/kg dry air) (displayed)
t0 <- h # temperature at zero humidity for enthalpy h


# enthalpy axis and enthalpy lines (black)
h <- seq(10, 110, 10) # enthalpy (kJ/kg-dry-air)

t0a <- h # temperature at zero humidity
t1a <- (h - 12.5)/3.5 # temperature on enthalpy axis
w2 <- t1a + 5 # specific humidity on enthalpy axis

zerocol1 <- rep_len(0, length(t1a))


# plotting with ggplot2
forplot1 <- data.frame(t, wg) # values to be plotted
x1 <- forplot1[,1] # x values
y1 <- forplot1[,2] # y values

forplot1 <- data.frame(forplot1, w)
y2 <- forplot1[, 3] # y values
y3 <- forplot1[, 4] # y values
y4 <- forplot1[, 5] # y values
y5 <- forplot1[, 6] # y values

forplot1 <- data.frame(forplot1, w1)
y6 <- forplot1[, 7] # y values
y7 <- forplot1[, 8] # y values


forplot2 <- data.frame(t1, tv0, wg1, zerocol)
yy1 <- forplot2[1, 1] # x values
yy2 <- forplot2[1, 2] # y values
yy3 <- forplot2[1, 3] # x values
yy4 <- forplot2[1, 4] # y values

yy5 <- forplot2[2, 1] # x values
yy6 <- forplot2[2, 2] # y values
yy7 <- forplot2[2, 3] # x values
yy8 <- forplot2[2, 4] # y values

yy9 <- forplot2[3, 1] # x values
yy10 <- forplot2[3, 2] # y values
yy11 <- forplot2[3, 3] # x values
yy12 <- forplot2[3, 4] # y values

yy13 <- forplot2[4, 1] # x values
yy14 <- forplot2[4, 2] # y values
yy15 <- forplot2[4, 3] # x values
yy16 <- forplot2[4, 4] # y values

yy17 <- forplot2[5, 1] # x values
yy18 <- forplot2[5, 2] # y values
yy19 <- forplot2[5, 3] # x values
yy20 <- forplot2[5, 4] # y values

yy21 <- forplot2[6, 1] # x values
yy22 <- forplot2[6, 2] # y values
yy23 <- forplot2[6, 3] # x values
yy24 <- forplot2[6, 4] # y values

yy25 <- forplot2[7, 1] # x values
yy26 <- forplot2[7, 2] # y values
yy27 <- forplot2[7, 3] # x values
yy28 <- forplot2[7, 4] # y values



forplot3 <- data.frame(t1, t0, wg1, zerocol)
yyy1 <- forplot3[1, 1] # x values
yyy2 <- forplot3[1, 2] # y values
yyy3 <- forplot3[1, 3] # x values
yyy4 <- forplot3[1, 4] # y values

yyy5 <- forplot3[2, 1] # x values
yyy6 <- forplot3[2, 2] # y values
yyy7 <- forplot3[2, 3] # x values
yyy8 <- forplot3[2, 4] # y values

yyy9 <- forplot3[3, 1] # x values
yyy10 <- forplot3[3, 2] # y values
yyy11 <- forplot3[3, 3] # x values
yyy12 <- forplot3[3, 4] # y values

yyy13 <- forplot3[4, 1] # x values
yyy14 <- forplot3[4, 2] # y values
yyy15 <- forplot3[4, 3] # x values
yyy16 <- forplot3[4, 4] # y values

yyy17 <- forplot3[5, 1] # x values
yyy18 <- forplot3[5, 2] # y values
yyy19 <- forplot3[5, 3] # x values
yyy20 <- forplot3[5, 4] # y values

yyy21 <- forplot3[6, 1] # x values
yyy22 <- forplot3[6, 2] # y values
yyy23 <- forplot3[6, 3] # x values
yyy24 <- forplot3[6, 4] # y values



forplot4 <- data.frame(cbind(t0a, t1a, zerocol1, w2))
yyyy1 <- forplot4[1, 1] # x values
yyyy2 <- forplot4[1, 2] # y values
yyyy3 <- forplot4[1, 3] # x values
yyyy4 <- forplot4[1, 4] # y values

yyyy5 <- forplot4[2, 1] # x values
yyyy6 <- forplot4[2, 2] # y values
yyyy7 <- forplot4[2, 3] # x values
yyyy8 <- forplot4[2, 4] # y values

yyyy9 <- forplot4[3, 1] # x values
yyyy10 <- forplot4[3, 2] # y values
yyyy11 <- forplot4[3, 3] # x values
yyyy12 <- forplot4[3, 4] # y values

yyyy13 <- forplot4[4, 1] # x values
yyyy14 <- forplot4[4, 2] # y values
yyyy15 <- forplot4[4, 3] # x values
yyyy16 <- forplot4[4, 4] # y values

yyyy17 <- forplot4[5, 1] # x values
yyyy18 <- forplot4[5, 2] # y values
yyyy19 <- forplot4[5, 3] # x values
yyyy20 <- forplot4[5, 4] # y values

yyyy21 <- forplot4[6, 1] # x values
yyyy22 <- forplot4[6, 2] # y values
yyyy23 <- forplot4[6, 3] # x values
yyyy24 <- forplot4[6, 4] # y values

yyyy25 <- forplot4[7, 1] # x values
yyyy26 <- forplot4[7, 2] # y values
yyyy27 <- forplot4[7, 3] # x values
yyyy28 <- forplot4[7, 4] # y values

yyyy29 <- forplot4[8, 1] # x values
yyyy30 <- forplot4[8, 2] # y values
yyyy31 <- forplot4[8, 3] # x values
yyyy32 <- forplot4[8, 4] # y values

yyyy33 <- forplot4[9, 1] # x values
yyyy34 <- forplot4[9, 2] # y values
yyyy35 <- forplot4[9, 3] # x values
yyyy36 <- forplot4[9, 4] # y values

yyyy37 <- forplot4[10, 1] # x values
yyyy38 <- forplot4[10, 2] # y values
yyyy39 <- forplot4[10, 3] # x values
yyyy40 <- forplot4[10, 4] # y values

yyyy41 <- forplot4[11, 1] # x values
yyyy42 <- forplot4[11, 2] # y values
yyyy43 <- forplot4[11, 3] # x values
yyyy44 <- forplot4[11, 4] # y values

forplot5 <- data.frame(c(0, 25, 5, 30))
yyyyy1 <- forplot5[1,]
yyyyy2 <- forplot5[2,]
yyyyy3 <- forplot5[3,]
yyyyy4 <- forplot5[4,]


pdf(file = "Simple_Psychrometric_Chart.pdf") # create the pdf file with the 
# given name
dc <- ggplot(data = forplot1, aes(x = x1, y = y1)) 
+ geom_line(linetype = 2, colour = "red") 
+ labs(list(title = "Simplified Psychrometric Chart at 1 atm total pressure", 
x = "Dry Bulb Temperature (°C)")) 
+ ylab(expression(paste("Specific Humidity (", omega, ") grams 
moisture / kilogram dry air"))) 
+ theme_bw() + geom_line(aes(x = x1, y = y2), linetype = 1, colour = "blue") 
+ geom_line(aes(x = x1, y = y3), linetype = 1, colour = "blue") 
+ geom_line(aes(x = x1, y = y4), linetype = 1, colour = "blue") 
+ geom_line(aes(x = x1, y = y5), linetype = 1, colour = "blue") 
+ geom_line(aes(x = x1, y = y6), linetype = 1, colour = "blue") 
+ geom_line(aes(x = x1, y = y7), linetype = 1, colour = "blue") 
+ geom_line(data = forplot2[1,], aes(x = c(yy1, yy2), y = c(yy3, yy4)), 
linetype = 2, colour = "green") 
+ geom_line(data = forplot2[2,], aes(x = c(yy5, yy6), y = c(yy7, yy8)), 
linetype = 2, colour = "green")
+ geom_line(data = forplot2[3,], aes(x = c(yy9, yy10), y = c(yy11, yy12)), 
linetype = 2, colour = "green") 
+ geom_line(data = forplot2[4,], aes(x = c(yy13, yy14), y = c(yy15, yy16)), 
linetype = 2, colour = "green") 
+ geom_line(data = forplot2[5,], aes(x = c(yy17, yy18), y = c(yy19, yy20)), 
linetype = 2, colour = "green") 
+ geom_line(data = forplot2[6,], aes(x = c(yy21, yy22), y = c(yy23, yy24)), 
linetype = 2, colour = "green") 
+ geom_line(data = forplot2[7,], aes(x = c(yy25, yy26), y = c(yy27, yy28)), 
linetype = 2, colour = "green")
+ geom_line(data = forplot3[1,], aes(x = c(yyy1, yyy2), y = c(yyy3, yyy4)), 
linetype = 2, colour = "red") 
+ geom_line(data = forplot3[2,], aes(x = c(yyy5, yyy6), y = c(yyy7, yyy8)),
linetype = 2, colour = "red") 
+ geom_line(data = forplot3[3,], aes(x = c(yyy9, yyy10), y = c(yyy11, yyy12)), 
linetype = 2, colour = "red") 
+ geom_line(data = forplot3[4,], aes(x = c(yyy13, yyy14), y = c(yyy15, yyy16)), 
linetype = 2, colour = "red") 
+ geom_line(data = forplot3[5,], aes(x = c(yyy17, yyy18), y = c(yyy19, yyy20)), 
linetype = 2, colour = "red") 
+ geom_line(data = forplot3[6,], aes(x = c(yyy21, yyy22), y = c(yyy23, yyy24)), 
linetype = 2, colour = "red") 
+ geom_line(data = forplot4[1,], aes(x = c(yyyy1, yyyy2), y = c(yyyy3, yyyy4)), 
linetype = 2, colour = "black") 
+ geom_line(data = forplot4[2,], aes(x = c(yyyy5, yyyy6), y = c(yyyy7, yyyy8)), 
linetype = 2, colour = "black") 
+ geom_line(data = forplot4[3,], aes(x= c(yyyy9, yyyy10), 
y = c(yyyy11, yyyy12)), linetype = 2, colour = "black") 
+ geom_line(data = forplot4[4,], aes(x = c(yyyy13, yyyy14), 
y = c(yyyy15, yyyy16)), linetype = 2, colour = "black") 
+ geom_line(data = forplot4[5,], aes(x = c(yyyy17, yyyy18), 
y = c(yyyy19, yyyy20)), linetype = 2, colour = "black") 
+ geom_line(data = forplot4[6,], aes(x = c(yyyy21, yyyy22), 
y = c(yyyy23, yyyy24)), linetype = 2, colour = "black") 
+ geom_line(data = forplot4[7,], aes(x = c(yyyy25, yyyy26), 
y = c(yyyy27, yyyy28)), linetype = 2, colour = "black") 
+ geom_line(data = forplot4[8,], aes(x = c(yyyy29, yyyy30), 
y = c(yyyy31, yyyy32)), linetype = 2, colour = "black") 
+ geom_line(data = forplot4[9,], aes(x = c(yyyy33, yyyy34), 
y = c(yyyy35, yyyy36)), linetype = 2, colour = "black") 
+ geom_line(data = forplot4[10,], aes(x = c(yyyy37, yyyy38), 
y = c(yyyy39, yyyy40)), linetype = 2, colour = "black") 
+ geom_line(data = forplot4[11,], aes(x = c(yyyy41, yyyy42), 
y = c(yyyy43, yyyy44)), linetype = 2, colour = "black")
+ geom_line(data = forplot5, aes(x = c(yyyyy1, yyyyy2), y = c(yyyyy3, yyyyy4)), 
linetype = 2, colour = "black") 
+ coord_cartesian(xlim = c(0,50), ylim = c(0, 30)); dc
# use ggplot2 to create the specified plot where the x values are factor values 
# with levels, the color and the shapes are based on the R object Temperature2, 
# the size of each shaped object is 4 and the x axis has the labels from 
# 0 to 100,000
dev.off() # close the plotting device

Answer 2

Here is how you would port the for loop in an idiomatic way:

phi <- seq(0.1, 0.4, 0.1)
pg <- c(0.61165, 0.65709, 0.70599, 0.75808, 0.81355, 0.87258, 0.93536, 1.00210,
        1.07300, 1.14830, 1.22820) 
p <- rep(phi, each=length(pg))
matrix(622 * p * pg/(patm - p * pg), nrow=length(pg))

No explicit looping needed.

Answer 3

Here is the exact translation. Did not attempt attempt to improve the implementation in R.

# Plotting the Simplified Psychrometric Chart
# Izzi Urieli 2/27/2008, translated to R on Oct 7, 2014

setwd("set the path where the two data files are located")

tpg = read.table(file = "t_pg.txt")
t = tpg[,1]
pg = tpg[,2]
patm = 101.325
rair = 0.287
wg = 622*pg/(patm-pg)
plot(t,wg,type="l", col = "red", xlim=c(0, 50), ylim=c(5, 30), main = "Simplified Psychrometric Chart",
     xlab = "Dry Bulb Temperature (deg C)",
     ylab = "Specific Humidity (gm vap/kg dry air)",
     xaxs="i", yaxs="i")

for (phi in seq(0.1,0.4,0.1)) {
  w = 622*phi*pg/(patm-phi*pg)
  lines(t,w)
}

for (phi in seq(0.6,0.8,0.2)) {
  w = 622*phi*pg/(patm-phi*pg);
  lines(t,w)
}

tpg1 = read.table("t_pg1.txt")
t1 = tpg1[,1]
pg1 = tpg1[,2]
wg1 = 622*pg1/(patm-pg1)
vol = rair*(t1+273)/(patm-pg1)
tv0 = patm*vol/rair-273

for (i in 1:7) {
  lines(c(t1[i],tv0[i]),c(wg1[i],0), type = "l", col = "green")
}

h = t1 + 2.5*wg1
t0 = h

for (i in 1:6) {
  lines(c(t1[i],t0[i]),c(wg1[i],0), type = "l", col = "red")
}

for (h in seq(10, 110, 10)) {
t0 = h
t1 = (h - 12.5)/3.5
w1 = t1 + 5
lines(c(t0,t1),c(0,w1), type = "l", col = "black")
}

lines(c(0,25), c(5,30), type = "l", col = "black")
grid(10,10)