---
title: "Using AeroSampleR to model aerosol sampling efficiency"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using AeroSampleR to model aerosol sampling efficiency}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
  markdown: 
    wrap: 72
---

```{=html}
<style>
body {
text-align: justify}
</style>
```

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## AeroSampleR Introduction

Air sampling systems are used in many applications that require the monitoring of hazardous airborne particles. When tubing is used to collect aerosol particles, some of the particles are lost along the way.  The efficiency of the system depends on the particle and tubing configuration. 

This version of AeroSampleR provides sampling efficiency for a limited set of system elements. Sampling systems always include a probe, after which are combinations of straight tubing and bends. While some systems include expansion or contraction elements, or sample splitters. These components are not covered in this version of AeroSampleR. The probe model is limited to a simple open ended pipe in still air. 

AeroSampleR relies on the concept of aerodynamic median activity diameter (AMAD), which accounts for particle density and shape, leaving equivalent spherical water droplets as the modeling targets. 

Efficiency functions are based predominantly on testing with aerosol particles through stainless steel tubing. The [Zhang](https://doi.org/10.1016/j.jaerosci.2012.05.007),
[McFarland](https://doi.org/10.1021/es960975c), and [Pui](https://doi.org/10.1080/02786828708959166) bend models are used in this package.

The aerosol transport models are based on tests on clean systems. This package is designed primarily for new tubing designs. If a system is not maintained clean -- and free of condensation, there can be no expectation that sampling efficiency models will be accurate. 

```{r echo=FALSE, message=FALSE, warning=FALSE}
library(AeroSampleR)
library(ggplot2)
library(dplyr)
library(flextable)
```

```{r echo=FALSE, message=FALSE, warning=FALSE, comment = ">"}
sys_df <- structure(list( el_num = c("1", "2", "3", "4"), 
    el_type = c("probe", "tube", "bend", "tube"), 
    length_cm = c(NA, 111.76, NA, 146.05), 
    angle_to_horiz = c(NA, 90, NA, 0), 
    orient = c("u", NA, NA, NA), 
    bend_angle = c(NA, NA, 90, NA), 
    bend_rad_cm = c(NA, NA, 12.7, NA)), 
    row.names = c(NA, -4L), 
    class = c("tbl_df", "tbl", "data.frame"))
cat("\n")
```

### example data
The first task in evaluating a system is to set up a table that includes all the "elements" with the following column headers:

- `el_num` sequential number of the element
- `el_type` starting with "probe", followed by "tube" and "bend" elements
- `length_cm` length of tubes in centimeters. Leave blank for probe and bends.
- `angle_to_horiz` degree of slope of straight tube elements. Leave blank for probes and bends.
- `orient` orientation of the probe. Options are "u" for up, "d" for down and "h" for horizontal
``bend_angle` how many degrees a bend turns the sample. Typically, 90 or 45. Leave blank for probes and tubes.
- `bend_rad_cm` the bend radius in centimeters. Leave blank for probes and tubes.

```{r echo=FALSE, message=FALSE, warning=FALSE, comment = ">"}
ft <- flextable(sys_df)
ft <- colformat_double(ft, digits = 0)
ft
```

### 1) Get the system data

You need to get this data into R and call it `sys_df`, the "system" data frame. There are many options on how to do this. Below are two examples. You will have to provide the path and the file name, but for the example, we'll show the file to be called system.txt or system.xlsx in `c:/work`.  Here are two:

a. Use base R, (`utils` package that is loaded with base R) to read a text file:

```{r echo=TRUE, eval=FALSE}
sys_df <- read.table(
  file = "c:/work/system.txt",
  header = TRUE
  )
```

b. Use the readxl package to read a spreadsheet of the 'xlsx' format:

```{r echo=TRUE, eval=FALSE}
sys_df <- readxl::read_xlsx(path = "c:/work/system.xlsx",
           sheet = "Sheet1", #default - update if needed
           range = "A1:G5", #put in entire range
           col_types = c("numeric",
                         "text",
                         "numeric",
                         "numeric",
                         "text",
                         "numeric",
                         "numeric")
                            )
```                            

### 2) Create particle distribution with `particle_dist()`

This function provides a logarithmic distribution of 1000 particle sizes and an additional set of discrete particles. By default, the logarithmically-distributed particles have an AMAD of 5 and a lognormal standard deviation of 2.5, consistent with ICRP 66. 

The discrete particles are 1, 5, and 10 micrometers AMAD.

```{r echo=TRUE, message=TRUE, warning=FALSE}
df <- particle_dist() #Default

```

```{r echo=FALSE, message=TRUE, warning=FALSE, fig.width=5, fig.height= 3}
df |> filter(dist == "log_norm") |> 
  ggplot(aes(D_p, dens)) + geom_point(color = "blue") +
  ggtitle("distribution of lognormal particle sizes")

df |> filter(dist == "log_norm") |> 
  mutate("activity" = D_p ^3 * dens)  |> 
  ggplot(aes(D_p, activity)) + geom_point(color = "blue") +
  ggtitle("relative activity by particle size",
          subtitle = "diameter cubed times density")

```


### 3) Set up the parameters for tube size, flow rate, temperature, and pressure. These parameters are not particle dependent and so can be kept in a small separate data frame.

-   D_tube is the inner diameter of the tube {cm}
-   Q_lpm is the flow rate of air {lpm}
-   T_C is the system pressure {Celsius}
-   P_kPa is the pressure of the system {kPa}

```{r echo=TRUE, message=TRUE, warning=FALSE}
# In this example the tubing wall is 1.65 mm thick. 
params <- set_params_1("D_tube" = 2.54 - (2 * 0.165), #1 inch tube diameter
                       "Q_lpm" = 2 * 28.3, #2 cfm converted to lpm
                       "T_C" = 25, 
                       "P_kPa" = 101.325)
```

Next, we compute the particle size-dependent parameters. These include factors for transport efficiency computation. - Cunningham Correction Factor {C_c} - terminal settling velocity {v_ts} - Particle Reynolds
number (tube) {Re_p} - Stokes number {Stk}

```{r echo=TRUE, message=TRUE, warning=FALSE}
df <- set_params_2(df, params)
```

At this point, our main particle distribution data frame has been modified with computed factors for use in the transport efficiency models, row by row.

### 4) Next, we compute the efficiency, element by element in transport order. We have only four elements in this example and we will evaluate them with `prob_eff()`, `tube_eff()`, `bend_eff()`, and lastly `tube_eff()` again. This will add columns to our particle data frame.

Calculate the efficiency of the probe via `prob_eff()` and add it to a new data frame. The orient argument sets the orientation of the probe. "u" means the probe is vertically upward. "d" is for a vertically downward facing probe. "h" is for a probe in a side configuration. The probe is in
the first row, so we use `[1]` to identify the orient parameter.

```{r echo=TRUE, message=TRUE, warning=FALSE}
df <- probe_eff(df, params, orient = sys_df$orient[1]) 
```

Calculate the efficiency of the first tube. Tube Efficiency is found using `tube_eff()` function. The length is given in cm (`length_cm`) and the angle from tube to horizontal orientation parameter
(`angle_to_horiz`) is specified here. All three parameters can be added to the above data frame, which will return a column for each distribution.

```{r echo=TRUE, message=TRUE, warning=FALSE}
df <- tube_eff(df, 
               params, 
               L = sys_df$length_cm[2] / 100, 
               angle_to_horiz = sys_df$angle_to_horiz[2], 
               elnum = sys_df$el_num[2])
```

Calculate the efficiency of the bend. Here, we'll take the Zhang model option. Bend Efficiency is found via the `bend_eff()` function and is where you will choose to use one of three different tube models {Zhang,
McFarland, or Pui}. The bend angle and element number are also listed in the function.

```{r echo=TRUE, message=TRUE, warning=FALSE}
df <- bend_eff(df, params, method = "Zhang", 
               bend_angle = sys_df$bend_angle[3],
               bend_radius = sys_df$bend_rad_cm[3] / 100, 
               elnum = sys_df$el_num[3])
```

Finally, we'll calculate transport efficiency through the last tube element.

```{r echo=TRUE, message=TRUE, warning=FALSE}
df <- tube_eff(df, params, L = sys_df$length_cm[2] / 100, 
               angle_to_horiz = sys_df$angle_to_horiz[4], 
               elnum = sys_df$el_num[4])
```

### At this point, the transport efficiencies have been built into the data frame. Let's have a look at the bottom few rows.
*It doesn't all fit horizontally, so match the top portion and the bottom portion by the row number.*
```{r echo=TRUE, message=TRUE, warning=FALSE}
tail(df)
```

### 5) Generate reports with `report_basic`, `report_plots` and `report_cum_plots`

The `report_basic` function provides total system efficiency for either all of the logarithmically distributed particles or all of the discrete particle sizes.

The `report_plots` function shows individual element efficiency. The `report_cum_plots` function shows cumulative efficiency through the system. This plot takes efficiency data from the rows of the data frame.
It therefore only works for individually selected particle sizes.

We'll show the parameter set first, so that the output message on the basic report, regarding units, makes sense.

```{r echo=TRUE, message=FALSE, warning=FALSE, fig.width=5, fig.height= 3, fig.align = 'center'}
params[, 7] <- formatC(params[, 7], digits = 2, format = "e")
params[, 8] <- formatC(params[, 8], digits = 2, format = "e")
params[, 11] <- formatC(params[, 11], digits = 2, format = "e")
params[, 3] <- formatC(params[, 3], digits = 4)
params[, 10] <- formatC(params[, 10], digits = 4)
ft <- flextable(params)
ft <- set_caption(ft, "system parameters")
ft
```

```{r echo=TRUE, message=FALSE, warning=FALSE, fig.width=5, fig.height= 3, fig.align = 'center'}

ft <- flextable(report_basic(df, params, "discrete"))
ft <- colformat_double(ft, digits = 3)
ft <- set_caption(ft, "results for discrete particle diameters")
ft

report_plots(df, "discrete")
report_cum_plots(df, 1)
report_cum_plots(df, 5)
report_cum_plots(df, 10)


ft <- flextable(report_basic(df, params, "log"))
ft <- colformat_double(ft, digits = 3)
ft <- set_caption(ft, "results for log distribution of particle diameters")
ft  

```

### 6) Optional extra reports for the logarithmically distributed particle set

The `report_log_mass` function provides details on every particle size.

Since there are 1000 data points, the full ouptput is probably not suitable for a typical report. The report provides the following columns of output: 

-  microns = the particle size in micrometers (microns is
shorter, but is considered supersed by micrometers) 
-  probs = relative probability of a particle being in this size bin 
-  bin_eff = the overall system efficiency for a particle of this size 
-  amb_mass = the probability of the particle multiplied by the mass of a spherical particle with the size given and density of 1 g per ml. This is the relative mass of this particle size in the ambient air being sampled.
-  sampled_mass = the relative mass that made it through the sampling system with this particle size 
-  bin_frac_lost = ambient mass in this bin minus the sampled mass in the bin, divided by the ambient mass
-  total_frac_lost = ambient mass in this bin minus the sampled mass in the bin, divided by the sum of the ambient mass

A random selection of ten rows of the 1000 rows are provided below:

```{r echo=TRUE, message=TRUE, warning=FALSE, fig.width=5, fig.height= 3, fig.align = 'center'}
df_log <- report_log_mass(df)[sort(sample(1:1000, 10)), ]
# need to make format changes so that flextable will show scientific notation
df_log[, 1] <- formatC(df_log[, 1], digits = 4)
df_log[, 2] <- formatC(df_log[, 2], digits = 2, format = "e")
df_log[, 3] <- formatC(df_log[, 3], digits = 2, format = "e")
df_log[, 4] <- formatC(df_log[, 4], digits = 2, format = "e")
df_log[, 5] <- formatC(df_log[, 5], digits = 2, format = "e")
df_log[, 6] <- formatC(df_log[, 6], digits = 2, format = "e")
df_log[, 7] <- formatC(df_log[, 7], digits = 2, format = "e")
ft <- flextable(df_log)
ft <- colformat_double(ft, digits = 3)
ft <- set_caption(ft, "results for random sample of 1000 particle diameters from the log set")
ft
```

The particle mass modeled in the ambient air and sampled through the air sampling system is shown with the function `report_plots`.

```{r echo=TRUE, message=TRUE, warning=FALSE, fig.width=5, fig.height= 3, fig.align = 'center'}
report_plots(df, "log")
```