diff --git a/vignettes/EnergyIndicators.Rmd b/vignettes/EnergyIndicators.Rmd
index f4a1a04b722a6b250efadc015808dc5691811969..5a65c4759ee4932b76b791e1f3787fa6dcafaf25 100644
--- a/vignettes/EnergyIndicators.Rmd
+++ b/vignettes/EnergyIndicators.Rmd
@@ -30,24 +30,38 @@ Although wind turbines cannot extract all of the kinetic energy in the wind, and
 As an example, we simulate a time series of 1000 wind speed values from a Weibull distribution with scale factor of 6 and a shape factor of 2, which represent a sample of wind speed values obtained at a single location. The Weibull distribution is often assumed to fit observed wind speed values to a probability distribution function. Then, each instantaneous wind speed value is converted to its equivalent WPD.
 The `mean` and `sd` of the WPD can be employed to summarize the wind resource in that location. Otherwise, we can plot the histograms to see the full distribution of values:
 
-```{r, fig.width=7}
+```
 library(CSIndicators)
 set.seed(1)
-oldpar <- par(no.readonly = TRUE) 
 wind <- rweibull(n = 1000, shape = 2, scale = 6)
 WPD <- WindPowerDensity(wind)
 mean(WPD)
+```
+
+```
+## [1] 170.6205
+```
+
+```
 sd(WPD)
+```
+
+```
+## [1] 251.1349
+```
+
+```
 par(mfrow = c(1, 2))
 hist(wind, breaks = seq(0, 20))
 hist(WPD, breaks = seq(0, 4000, 200))
 ```
+![WPD](./Figures/WPD_histogram.png)
 
 As you can see the histogram of the WPD is highly skewed, even if the wind speed was only a little skewed!
 
 If not specified, an air density of 1.225 kg/m^3 is assumed. Otherwise, the parameter `ro` can be set to a fixed value (for instance the mean air density at the site elevation could be used), or a timeseries of density values measured at each time stamp can be used to obtain more accurate results.
 
-```{r}
+```
 WPD <- WindPowerDensity(wind, ro = 1.15)
 ```
 
@@ -61,16 +75,16 @@ Notice that power curves are intended to be used with 10-minutal steady wind spe
 
 Following on the previous example, we will compute now the CF that would be obtained from our sample of 1000 wind speed values when using a turbine of class IEC I, and compare it to the CF values for a class III:
 
-```{r, fig.width=7}
+```
 WCFI <- WindCapacityFactor(wind, IEC_class = "I")
 WCFIII <- WindCapacityFactor(wind, IEC_class = "III")
 par(mfrow = c(1, 3))
 hist(wind, breaks = seq(0, 20))
 hist(WCFI, breaks = seq(0, 1, 0.05), ylim = c(0, 500))
 hist(WCFIII, breaks = seq(0, 1, 0.05), ylim = c(0, 500))
-par(oldpar) 
 ```
 
+![WCF](./Figures/WCF_histogram.png)
 
 From the CF histograms we can see that, for this particular wind speed distribution, the IEC I turbine (designed for high winds) producess less energy than the IEC III turbine, which is more suitable for this range of wind speed values.