# set working directory to source file location!

# make sure to install current SPOT version:
# install.packages("SPOT")
# then load SPOT with:
require(SPOT)

# load current version of cyclone analytical model functions
source('cyclone.R')

### Standard example from loeffler
funCyclone(c(1.26,2.5,.42,.65,.6,.2),fluid=list(Mu=1.85e-5,Ve=(50/36)/0.12,lambdag=1/200,Rhop=2000,Rhof=1.2,Croh=0.05))

### Data for the present case study:
## intervals of the particle size distribution
intervals <- c(0,1.,2.7,5.5,8.7,12.7,16.9,21.1,25.4,30.8,63)*1e-6 
## values for each interval
delta <- rep(0.1, length(intervals)-1)

h <- 160/1000
g <- read.csv("bart16eAllData.csv")

## Remove outliers:
g <- g[g$E < 100 & g$E > 80,]
g <- droplevels(g)

res <- c()
for (i in 1:17){
  df <- data.frame(g[i,1:7])
  be <- df$be/1000
  Da <- df$Da/1000
  Dt <- df$Dt/1000
  he <- df$he/1000
  ht <- df$ht/1000
  
  ### Process parameters:
  ## Inlet gas velocity # [m/s]
  ve <- 20

  ## Wall friction coefficient   
  lambdag <- 1/200
  ## viscosity
  mu <-  1.8e-5; # [Pa s] = [kg/(m s)] 
  ## gas density
  rhof <- 1.2 # [kg/m^3] 
  ## particle density
  rhop <- 2700 # [kg/m^3] 
  
  ## gas flow rate
  vp <- ve*be*he # [m^3 / s]
  
  ## raw gas concentration: 6g dust particles in 10 seconds
  croh <- 6*10^(-3)/(vp * 10)   # [kg / m^3] 
  
  ##  calculate collection efficiency E
  y <- funCyclone(x=c(Da,h,Dt,ht,he,be),
             fluid=list(Mu=mu,Ve=ve,lambdag=lambdag,Rhop=rhop,Rhof=rhof,Croh=croh),
             intervals=intervals,
             delta=delta)
  
  res <- c(res, -y[3] * 100) 
}
## print all results
print(res) 

## build linear model
X <- g[1:17,1:7]
df1 <- data.frame(X,res)
names(df1)
lm1 <- lm(res ~ he + be + Dt + ht  + Da , data=df1)
summary(lm1)

## refine linear model with stepwise AIC
require(MASS) # if MASS is not available: install.packages("MASS")
lm2 <- stepAIC(lm1, scope= list(upper = ~ he + be + Dt + ht + Da , lower = ~ 1), trace=TRUE)
## this model is shown in eq. 6 in the article:
summary(lm2) 
library("xtable")
xtable(lm2)


#######################
### Table 5: (M-A):

res <- c()
for (i in 1:3){  # Loeffler, Muschelknautz, Stairmand
  for(ht in c(0, 35/1000, 44/1000)){
    df <- data.frame(g[i,1:7])
    be <- df$be/1000
    Da <- df$Da/1000
    Dt <- df$Dt/1000
    he <- df$he/1000
    
    ### Process parameters:
    ## Inlet gas velocity # [m/s]
    ve <- 20
    
    ## Wall friction coefficient   
    lambdag <- 1/200
    ## viscosity
    mu <-  1.8e-5; # [Pa s] = [kg/(m s)] 
    ## gas density
    rhof <- 1.2 # [kg/m^3] 
    ## particle density
    rhop <- 2700 # [kg/m^3] 
    
    ## gas flow rate
    vp <- ve*be*he # [m^3 / s]
    
    ## raw gas concentration: 6g dust particles in 10 seconds
    croh <- 6*10^(-3)/(vp * 10)   # [kg / m^3] 
    
    ##  calculate collection efficiency E
    y <- funCyclone(x=c(Da,h,Dt,ht,he,be),
                    fluid=list(Mu=mu,Ve=ve,lambdag=lambdag,Rhop=rhop,Rhof=rhof,Croh=croh),
                    intervals=intervals,
                    delta=delta)
    
    res <- c(res, -y[3] * 100) #Ew.
  }
}
print(res) # Column (M-A)