comfort_R: Models and Equations for Human Comfort Research

Documented in calcHbExUnsteady

#' @title Human Body Exergy Consumption Rate using Unsteady State Method
#' @description \code{calcHbExUnsteady} Function calculates the human body exergy consumPtion rate using unsteady state method based on a series of environmental variables.
#' @aliases calcHbExUnsteady
#' @aliases hbexunsteady
#' @aliases exunsteady
#' @usage calcHbExUnsteady(ta, tr, rh, vel, clo, met, tao, rho, frad = 0.7, 
#' eps = 0.95, ic = 1.085, ht = 171, wt = 70, tcr = 37, tsk = 36, basMet = 58.2,
#' warmUp = 60, cdil = 100, sigmatr = 0.25, dateTime)
#' @param ta a numeric value presenting air temperature in [degree C]
#' @param tr a numeric value presenting mean radiant temperature in [degree C]
#' @param vel a numeric value presenting air velocity in [m/s]
#' @param rh a numeric value presenting relative humidity [\%]
#' @param clo a numeric value presenting clothing insulation level in [clo]
#' @param met a numeric value presenting metabolic rate in [met]
#' @param tao a numeric vector presenting outdoor air temperature in [degree C]
#' @param rho numeric vector presenting outdoor relative humidity [\%]
#' @param frad a numeric vector presenting the fraction of body exposed to radiation 0.7(for seating), 0.73(for standing) [-]
#' @param eps a numeric vector presenting emissivity [-]
#' @param ic a numeric vector presenting permeability of clothing: 1.084 (average permeability), 0.4 (low permeability)
#' @param ht a numeric vector presenting body height in [cm]
#' @param wt a numeric vector presenting body weight in [kg]
#' @param tcr a numeric vector presenting initial value for core temperature in [degree C]
#' @param tsk a numeric vector presenting initial value for skin temperature in [degree C]
#' @param basMet a numeric vector presenting basal metabolic rate in [met]
#' @param warmUp a numeric vector presenting length of warm up period, i.e. number of times, loop is running for HBx calculation.
#' @param cdil a numeric vector presenting value for cdil in 2-node model of Gagge
#' @param sigmatr a numeric vector presenting value for cdil in 2-node model of Gagge.
#' @param dateTime a POsIxct vector of the times of measurement.
#' @details This function requires vectors of data including the corresponding time stamp. In case the time between two measurements is more than a minute, intermediate values are interpolated.
#' @returns
#' Returns a data.frame with the following columns.
#' Exergy input
#' \item{xInmetu }{Exergy input through metabolism}
#' \item{xInmetwcu }{Label warm/ cold for exergy input through metabolism}
#' \item{xInAIRwcu }{Exergy input through inhaled humid air}
#' \item{xInAIRwcwcu }{Label warm/ cold for exergy input through inhaled humid air}
#' \item{xInAIRwdu }{Exergy input through inhaled dry air}
#' \item{xInAIRwdwdu }{Label wet/ dry for exergy input through inhaled dry air}
#' \item{xInLUNGwcu }{Exergy input through water lung}
#' \item{xInLUNGwcwcu }{Label warm/ cold for exergy input through water lung}
#' \item{xInLUNGwdu }{Exergy input through water lung}
#' \item{xInLUNGwdwdu }{Label wet/ dry for exergy input through water lung}
#' \item{xInsheLLwcu }{Exergy input through water from sweat}
#' \item{xInsheLLwcwcu }{Label warm/ cold for exergy input through water from sweat}
#' \item{xInsheLLwdu }{Exergy input through water from sweat}
#' \item{xInsheLLwdwdu }{Label wet/ dry for exergy input through water from sweat}
#' \item{xInraDu }{Exergy input through radiation}
#' \item{xInraDwcu }{Label warm/ cold for exergy input through radiation}
#' \item{xIntotaLu }{total exergy input}
#' Exergy output
#' \item{xoutstorecoreu }{Exergy stored in core}
#' \item{xoutstoreshelu }{Exergy stored in shell}
#' \item{xoutaIRwcu }{Exergy output through exhaled humid air}
#' \item{xoutaIRwcwcu }{Label warm/ cold for exergy output through exhaled humid air}
#' \item{xoutaIRwdu }{Exergy output through exhaled dry air}
#' \item{xoutaIRwdwdu }{Label wet/ dry for exergy output through exhaled dry air}
#' \item{xoutswEATwcu }{Exergy output through water vapour from sweat}
#' \item{xoutswEATwcwcu }{Label warm/ cold for exergy output through water vapour from sweat}
#' \item{xoutswEATwdu }{Exergy output through water vapour from sweat}
#' \item{xoutswEATwdwdu }{Label wet/ dry for exergy output through water vapour from sweat}
#' \item{xoutraDu }{Exergy output through radiation}
#'\item{xoutraDwcu }{Label warm/ cold for exergy output through radiation}
#' \item{xoutCONVu }{Exergy output through convection}
#' \item{xoutCONVwcu }{Label warm/ cold for exergy output through convection}
#' \item{xouttotaLu }{total exergy output}
#' Exergy balance
#' \item{xconsu }{total exergy consumPtion}
#' Additional values
#' \item{tsku }{Calculated skin temperature}
#' \item{tcru }{Calculated core temperature}
#' \item{wu }{Calculated skin wettedness}
#' @examples
#' ## Define environmental parameters
#' ta <- seq(20,25,.1)
#' tr <- ta
#' rh <- rep(50, length(ta))
#' vel <- rep(.1, length(ta))
#' clo <- rep(.8, length(ta))
#' met <- rep(1.2, length(ta))
#' tao <- rep(5, length(ta))
#' rho <- rep(80, length(ta))
#' dateTime <- as.POSIXct(seq(0,by=60,length.out=length(ta)), origin="1970-01-01")
#' ## Calculation of human body exergy consumPtion rate
#' calcHbExUnsteady(ta, tr, rh, vel, clo, met, tao, rho, dateTime = dateTime)$xconsu
#' @author This function is based on a VBA code developed by masanori Shukuya. transformation of VBA-code and Excel procedures into R syntax by Marcel Schweiker.
#' @seealso see also \code{\link{calcComfInd}}
#' @references 
#' Schweiker, Kolarik, Dovjak & Shukuya (2016) <doi:10.1016/j.enbuild.2016.01.002>
#' 
#' Shukuya (2015) Calculation of human body-core and skin-layer temperatures under 
#' unsteady-state conditions-for unsteady-state human-body exergy analysis-, 
#' internal report of exergy-research group, Tech. rep.
#' @note According to Gagge's paper (1973), the value of 'cdil' may vary between 75 and 225 and 'sigma-tr' between 0.25 and 0.75. There is a note in the appendix of his paper saying two things:  1) whatever the values taken for cdil and sigma-tr, there must be no significant change in resulting thermal equilibrium. But, the values taken for cdil and sigmaTr do affect time to equilibrium. According to the analysis of schweiker et al. (2015), the values of 100 and .25 lead to the best fit of calculated and observed skin temperature.
#' @export

calcHbExUnsteady <- function(ta, tr, rh, vel, clo, met, tao, rho, frad = .7, eps = .95,  ic = 1.085, ht=171, wt=70, tcr=37, tsk=36, basMet=58.2, warmUp=60, cdil=100, sigmatr=.25, dateTime){

  ## definition of output variables
  ## Exergy input
  xInmetu <- xInmetwcu <- xInAIRwcu <- xInAIRwcwcu <- xInAIRwdu <- xInAIRwdwdu <- xInLUNGwcu <- xInLUNGwcwcu <- xInLUNGwdu <- xInLUNGwdwdu <- xInsheLLwcu <- xInsheLLwcwcu <- xInsheLLwdu <- xInsheLLwdwdu <- xInraDu <- xInraDwcu <- xIntotaLu <- NA

  ## Exergy output
  xoutstorecoreu <- xoutstoreshelu <- xoutaIRwcu <- xoutaIRwcwcu <- xoutaIRwdu <- xoutaIRwdwdu <- xoutswEATwcu <- xoutswEATwcwcu <- xoutswEATwdu <- xoutswEATwdwdu <- xoutraDu <- xoutraDwcu <- xoutCONVu <- xoutCONVwcu <- xouttotaLu <- NA

  ## balance and additional variables
  xconsu <- tsku <- tcru <- wu <- NA

  tko <- 273.15
  tskSet <- 33.7
  tcrSet <- 36.8
  lr <- 16.5 * 10  ^  (-3)
  csw <- 170
  mAir <- 28.97 * 0.001;
  mWater <- 18.015 * 0.001;
  rGas <- 8.31446 #[J/(molK)]
  Row <- 1000; Roa <- 1.2#[kg/m3]
  cpBody <- 3490; cpa <- 1005; cpv <- 1846; cpw <- 4186 #[J/(kgK)]
  PO <- 101325#[N/m2<-J/m3]

  ## radiation area of human body - taking account for covered parts by other parts of body (e.g. inner parts of arm) - coefficients taken from Fanger
  aBody <- 0.008883 * ht  ^  0.663 * wt  ^  0.444
  rMA <- wt / aBody

  ## A series of calculation timestep by timestep
  ##  adjust here for input of individual conditions in each line

  runs <- length(ta)

  i <- 1

  icl <- clo[i] # clothing insulation [clo]
  va <- vel[i] # Indoor air velocity [m/s]
  metrate <- met[i] # metabolic rate [met]
  tmr <- tr[i] # mean radiant temp [degree C]
  tair <- ta[i] # air temp [degree C]
  phia <- rh[i] # relative humidity indoors [%]
  toEnv <- tao[i] # outdoor temp [degree C]
  phiaEnv <- rho[i] # relative humidity outdoors [%]

  ## warm-up period - see also discussion in paper
  for (j in 1:warmUp){
    # extract single value from time series data
    dtCal <- 60
    qmet <- metaTherm(metrate, basMet)
    fcl <- 1 + 0.15 * icl
    rcl <- 0.155 * icl

    TKa <- tair + tko; tkmr <- tmr + tko

    tkoEnv <- toEnv + tko
    pvEnv <- pVapor(toEnv, phiaEnv) # function see Book Shukuya p. 268

    hr <- 6.13 * frad * eps # radiative heat transfer coefficient
    hc <- hcvG(va, metrate, basMet) # convective heat transfer coefficient
    top <- (hr * tmr + hc * tair) / (hr + hc)
    ffcl <- 1 / (1 + 0.155 * icl * fcl * (hr + hc));
    fpcl <- 1 / (1 + 0.155 * icl * fcl * hc / ic) # related to evaporation
    cl <- 1 / (1 + 0.155 * icl * fcl * hc / ic)
    im <- hc * fpcl / ((hr + hc) * ffcl)

    iclStar <- 0.6; fclStar <- 1 + 0.15 * iclStar
    hcStar <- hcvG(0.1, 1, basMet)
    ffclStar <- 1 / (1 + 0.155 * iclStar * fclStar * (hr + hcStar))
    fpclStar <- 1 / (1 + 0.155 * iclStar * fclStar * hcStar / ic)
    imStar <- hcStar * fpclStar / ((hr + hcStar) * ffclStar);
    cres <- 0.0014 * (basMet * metrate) * (34 - tair) # heat transfer to environment by convection in relation to respiration (empirical equation by Fanger)

    pva <- pVapor(tair, phia)
    eres <- 0.0173 * (basMet * metrate) * (5.87 - pva / 1000) # heat removed by evaporation in relation to respiration (empirical equation by Fanger)


    qShiv <- mshiv(tcrSet, tskSet, tcr, tsk)
    vblS <- vblCdilStr(cdil, sigmatr, tcrSet, tskSet, tcr, tsk)# blood flow rate
    # next line is core of Gagge model. If blood flow becomes lower, than skin layer gets more dominant
    alfaSk <- 0.0418 + 0.745 / (vblS + 0.585); # factor to adjust for difference in dominance
    kS <- 5.28 + 1.163 * vblS # conductance between core and skin layer
    Qcr <- (1 - alfaSk) * rMA * cpBody; Qsk <- alfaSk * rMA * cpBody # heat capacity of core and skin layer

    tcrN <- (1 - dtCal * kS / Qcr) * tcr + dtCal / Qcr * (qmet + qShiv - (cres + eres) + kS * tsk) # core temperature at time step before step calculated
    psks <- pVapor(tsk, 100)# saturated water vapour pressure at skin surface calculated with pure water , i.e. adaptive processes such as less salt in sweat might be put heresee also p. 281 of book Shukuya
    emax <- fpcl * (fcl * lr * hc) * (psks - pva) # max rate of water dispersion when the whole skin surface is wet
    tb <- alfaSk * tsk + (1 - alfaSk) * tcr; # average body temperature using calculated tsk and tcr
    tbSet <- alfaSk * tsk + (1 - alfaSk) * tcrSet # average body temperature using set-point values for tsk and tcr
    ersw <- csw * (tb - tbSet) * exp((tsk - tskSet) / 10.7) # amount of sweat secretion

    w <- 0.06 + 0.94 * ersw / emax
    if (w < 0.06){
      w <- 0.06
    }
    if (1 < w){
      w <- 1
    }

    DTQ <- dtCal / Qsk
    tskN <- (1 - DTQ * kS - DTQ * fcl * ffcl * (hr + hc)) * tsk - DTQ * w * fcl * lr * hc * fpcl * psks + DTQ * (kS * tcr + fcl * (hr + hc) * ffcl * top + w * fcl * lr * hc * fpcl * pva) # tsk in next step
    tclN <- ((1 / rcl) * tskN + fcl * hr * tmr + fcl * hc * tair) / (1 / rcl + fcl * hr + fcl * hc)

    ## Exergy balance

    ## Thermal exergy generation by metabolism

    tkcrN <- tcrN + tko; tkskN <- tskN + tko; tkclN <- tclN + tko
    metaEnergy <- qmet + qShiv
    xMet <- metaEnergy * (1 - tkoEnv / tkcrN);
    xMetwc <- wcXCheck(tkcrN, tkoEnv)

    ## Inhaled humid air

    Vin <- 1.2 * 10  ^  (-6) * metaEnergy # Volume of air intake [V/s]
    cpav <- cpa * (mAir / (rGas * TKa)) * (PO - pva) + cpv * (mWater / (rGas * TKa)) * pva
    xInhaleWc <- Vin * wcEx(cpav, TKa, tkoEnv);
    xInhaleWcwc <- wcXCheck(TKa, tkoEnv)

    xInhaleWd <- Vin * wdEx(TKa, tkoEnv, pva, pvEnv);
    xInhaleWdwd <- wdXCheck(pva, pvEnv)

    ## Liquid water generated in the core by metabolism to be dispersed into the exhaled air

    VwCore <- Vin * (0.029 - 0.049 * 10  ^  (-4) * pva)
    xLwCoreWc <- VwCore * Roa * wcEx(cpw, tkcrN, tkoEnv); xLwCoreWcwc <- wcXCheck(tkcrN, tkoEnv)

    pvs_env <- pVapor(toEnv, 100)
    xLwCoreWet <- VwCore * Roa * rGas / mWater * tkoEnv * log(pvs_env / pvEnv)
    xLwCoreWetwd <- "wet"

    ## Liquid water generated in the shell by metabolism to be secreted as sweat

    vwShellRow <- w * emax / (2450 * 1000)
    xLwShellWc <- vwShellRow * wcEx(cpw, tkskN, tkoEnv); xLwShellWcwc <- wcXCheck(tkskN, tkoEnv)

    xLwShellWd <- vwShellRow * wdExLw(tkoEnv, pvs_env, pva, pvEnv);
    xLwShellWdwd <- wdXCheck(pva, pvEnv)

    ## radiant exergy input

    xInRad <- fcl * hr * (tkmr - tkoEnv)  ^  2 / (tkmr + tkoEnv);
    xInRadwc <- wcXCheck(tkmr, tkoEnv)

    ## total exergy input

    xIntotal <- xMet + xInhaleWc + xInhaleWd + xLwCoreWc + xLwCoreWet + xLwShellWc + xLwShellWd + xInRad

    ## Exergy stored

    xStCore <- Qcr * (1 - tkoEnv / tkcrN) * (tcrN - tcr) / dtCal
    xStShell <- Qsk * (1 - tkoEnv / tkskN) * (tskN - tsk) / dtCal

    ## Exhaled humid air

    pvssCr <- pVapor(tcrN, 100)
    cpav <- cpa * (mAir / (rGas * tkcrN)) * (PO - pvssCr) + cpv * (mWater / (rGas * tkcrN)) * pvssCr
    xExhaleWc <- Vin * wcEx(cpav, tkcrN, tkoEnv);
    xExhaleWcwc <- wcXCheck(tkcrN, tkoEnv)

    xExhaleWd <- Vin * wdEx(tkcrN, tkoEnv, pvssCr, pvEnv); xExhaleWdwd <- wdXCheck(pvssCr, pvEnv)

    ## water vapor originating from the sweat and humid air containing the evaporated sweat

    xSweatWc <- vwShellRow * wcEx(cpv, tkclN, tkoEnv);
    xSweatWcwc <- wcXCheck(tkclN, tkoEnv)

    xSweatWd <- vwShellRow * wdExLw(tkoEnv, pva, pva, pvEnv);
    xSweatWdwd <- wdXCheck(pva, pvEnv)

    ## radiant exergy output

    xOutRad <- fcl * hr * (tkclN - tkoEnv)  ^  2 / (tkclN + tkoEnv);
    xOutRadwc <- wcXCheck(tkclN, tkoEnv)

    ## Exergy transfer by convection

    xOutConv <- fcl * hc * (tkclN - TKa) * (1. - tkoEnv / tkclN);
    xOutConvwc <- wcXCheck(tkclN, tkoEnv)

    xouttotal <- xStCore + xStShell + xExhaleWc + xExhaleWd + xSweatWc + xSweatWd + xOutRad + xOutConv

    xConsumption <- xIntotal - xouttotal

    #etStar <- calcet(top, tair, phia, w, im, 50, im);

    tcr <- tcrN
    tsk <- tskN

  }

  ## Output values
  ## Exergy input
  xInmetu[i] <- xMet #metabolism
  xInmetwcu[i] <- xMetwc #metabolism warm/cold
  xInAIRwcu[i] <- xInhaleWc	# inhaled humid air
  xInAIRwcwcu[i] <- xInhaleWcwc
  xInAIRwdu[i] <- xInhaleWd #
  xInAIRwdwdu[i] <- xInhaleWdwd # wet/dry
  xInLUNGwcu[i] <- xLwCoreWc # water lung
  xInLUNGwcwcu[i] <- xLwCoreWcwc
  xInLUNGwdu[i] <- xLwCoreWet
  xInLUNGwdwdu[i] <- xLwCoreWetwd
  xInsheLLwcu[i] <- xLwShellWc # water sweat
  xInsheLLwcwcu[i] <- xLwShellWcwc
  xInsheLLwdu[i] <- xLwShellWd
  xInsheLLwdwdu[i] <- xLwShellWdwd
  xInraDu[i] <- xInRad # radiation in
  xInraDwcu[i] <- xInRadwc
  xIntotaLu[i] <- xIntotal # totaL exergy input

  ## Exergy output
  xoutstorecoreu[i] <- xStCore # exergy stored
  xoutstoreshelu[i] <- xStShell
  xoutaIRwcu[i] <- xExhaleWc # exhaled humid air
  xoutaIRwcwcu[i] <- xExhaleWcwc
  xoutaIRwdu[i] <- xExhaleWd
  xoutaIRwdwdu[i] <- xExhaleWdwd
  xoutswEATwcu[i] <- xSweatWc # water vapour
  xoutswEATwcwcu[i] <- xSweatWcwc
  xoutswEATwdu[i] <- xSweatWd
  xoutswEATwdwdu[i] <- xSweatWdwd
  xoutraDu[i] <- xOutRad # radiation out
  xoutraDwcu[i] <- xOutRadwc
  xoutCONVu[i] <- xOutConv # convection
  xoutCONVwcu[i] <- xOutConvwc
  xouttotaLu[i] <- xouttotal # totaL exergy out

  ## balance
  xconsu[i] <- xConsumption # exergy consumPtion total
  tsku[i] <- tsk
  tcru[i] <- tcr
  wu[i] <- w

  ## from here start dynamic series of data from row 2 of dataset
  for (i in 2:runs){

    ## Calc Dtcal dependend on timediff between rows + approx if more than 20 seconds time difference
    dtCal <-  as.numeric(difftime(dateTime[i], dateTime[i-1], units="secs")) # in seconds

    if(dtCal > 20){ # if time difference is greater than 20 seconds (leading to the calcs getting unstable), intermediate values or calculated

      nRunsInter <- trunc(dtCal/20-1, 0) # calc number of minutes in between interval
      dtCal <- dtCal / (nRunsInter+1)

      for (k in 1: nRunsInter) {

        ## get linear approximated values for intermediate times
        icl <- clo[i-1] + k*((clo[i] - clo[i-1])/(nRunsInter+1)) # clothing insulation [clo]
        va <- vel[i-1] + k*((vel[i] - vel[i-1])/(nRunsInter+1)) # Indoor air velocity [m/s]
        metrate <- met[i-1] + k*((met[i] - met[i-1])/(nRunsInter+1)) # metabolic rate [met]
        tmr <- tr[i-1] + k*((tr[i] - tr[i-1])/(nRunsInter+1)) # mean radiant temp [degree C]
        tair <- ta[i-1] + k*((ta[i] - ta[i-1])/(nRunsInter+1)) # air temp [degree C]
        phia <- rh[i-1] + k*((rh[i] - rh[i-1])/(nRunsInter+1)) # relative humidity indoors [%]
        toEnv <- tao[i-1] + k*((tao[i] - tao[i-1])/(nRunsInter+1)) # outdoor temp [degree C]
        phiaEnv <- rho[i-1] + k*((rho[i] - rho[i-1])/(nRunsInter+1)) # relative humidity outdoors [%]

        qmet <- metaTherm(metrate, basMet)
        fcl <- 1 + 0.15 * icl
        rcl <- 0.155 * icl

        TKa <- tair + tko; tkmr <- tmr + tko

        tkoEnv <- toEnv + tko
        pvEnv <- pVapor(toEnv, phiaEnv) #--------function see Book Shukuya p. 268

        hr <- 6.13 * frad * eps # radiative heat transfer coefficient
        hc <- hcvG(va, metrate, basMet) # convective heat transfer coefficient
        top <- (hr * tmr + hc * tair) / (hr + hc)
        ffcl <- 1 / (1 + 0.155 * icl * fcl * (hr + hc));
        fpcl <- 1 / (1 + 0.155 * icl * fcl * hc / ic)
        cl <- 1 / (1 + 0.155 * icl * fcl * hc / ic)
        im <- hc * fpcl / ((hr + hc) * ffcl)

        iclStar <- 0.6; fclStar <- 1 + 0.15 * iclStar
        hcStar <- hcvG(0.1, 1, basMet)
        ffclStar <- 1 / (1 + 0.155 * iclStar * fclStar * (hr + hcStar))
        fpclStar <- 1 / (1 + 0.155 * iclStar * fclStar * hcStar / ic)
        imStar <- hcStar * fpclStar / ((hr + hcStar) * ffclStar);
        cres <- 0.0014 * (basMet * metrate) * (34 - tair) # heat transfer to environment by convection in relation to respiration (empirical equation by Fanger)

        pva <- pVapor(tair, phia)
        eres <- 0.0173 * (basMet * metrate) * (5.87 - pva / 1000) # heat removed by evaporation in relation to respiration (empirical equation by Fanger)

        qShiv <- mshiv(tcrSet, tskSet, tcr, tsk)
        vblS <- vblCdilStr(cdil, sigmatr, tcrSet, tskSet, tcr, tsk)# blood flow rate
        ## next line is core of Gagge model. If blood flow becomes lower, than skin layer gets more dominant
        alfaSk <- 0.0418 + 0.745 / (vblS + 0.585); # factor to adjust for difference in dominance
        kS <- 5.28 + 1.163 * vblS # conductance between core and skin layer
        Qcr <- (1 - alfaSk) * rMA * cpBody; Qsk <- alfaSk * rMA * cpBody # heat capacity of core and skin layer

        tcrN <- (1 - dtCal * kS / Qcr) * tcr + dtCal / Qcr * (qmet + qShiv - (cres + eres) + kS * tsk) # core temperature at time step before step calculated
        psks <- pVapor(tsk, 100)# saturated water vapour pressure at skin surface calculated with pure water , i.e. adaptive processes such as less salt in sweat might be put heresee also p. 281 of book Shukuya
        emax <- fpcl * (fcl * lr * hc) * (psks - pva) # max rate of water dispersion when the whole skin surface is wet
        tb <- alfaSk * tsk + (1 - alfaSk) * tcr; # average body temperature using calculated tsk and tcr
        tbSet <- alfaSk * tsk + (1 - alfaSk) * tcrSet # average body temperature using set-point values for tsk and tcr
        ersw <- csw * (tb - tbSet) * exp((tsk - tskSet) / 10.7) # amount of sweat secretion

        w <- 0.06 + 0.94 * ersw / emax
        if (w < 0.06){
          w <- 0.06
        }
        if (1 < w){
          w <- 1
        }

        DTQ <- dtCal / Qsk
        tskN <- (1 - DTQ * kS - DTQ * fcl * ffcl * (hr + hc)) * tsk - DTQ * w * fcl * lr * hc * fpcl * psks + DTQ * (kS * tcr + fcl * (hr + hc) * ffcl * top + w * fcl * lr * hc * fpcl * pva) # tsk in next step
        tclN <- ((1 / rcl) * tskN + fcl * hr * tmr + fcl * hc * tair) / (1 / rcl + fcl * hr + fcl * hc)

        #etStar <- calcet(top, tair, phia, w, im, 50, im);

        tcr <- tcrN
        tsk <- tskN

      } # end for k times in between interval

    } # end if Dtcal >20

    ## extract single value from time series data
    icl <- clo[i] # clothing insulation [clo]
    va <- vel[i] # Indoor air velocity [m/s]
    metrate <- met[i] # metabolic rate [met]
    tmr <- tr[i] # mean radiant temp [degree C]
    tair <- ta[i] # air temp [degree C]
    phia <- rh[i] # relative humidity indoors [%]
    toEnv <- tao[i] # outdoor temp [degree C]
    phiaEnv <- rho[i] # relative humidity outdoors [%]

    qmet <- metaTherm(metrate, basMet)
    fcl <- 1 + 0.15 * icl
    rcl <- 0.155 * icl

    TKa <- tair + tko; tkmr <- tmr + tko

    tkoEnv <- toEnv + tko
    pvEnv <- pVapor(toEnv, phiaEnv) # function see Book Shukuya p. 268

    hr <- 6.13 * frad * eps # radiative heat transfer coefficient
    hc <- hcvG(va, metrate, basMet) # convective heat transfer coefficient #
    top <- (hr * tmr + hc * tair) / (hr + hc)
    ffcl <- 1 / (1 + 0.155 * icl * fcl * (hr + hc));
    fpcl <- 1 / (1 + 0.155 * icl * fcl * hc / ic) # related to evaporation
    cl <- 1 / (1 + 0.155 * icl * fcl * hc / ic)
    im <- hc * fpcl / ((hr + hc) * ffcl)

    iclStar <- 0.6; fclStar <- 1 + 0.15 * iclStar
    hcStar <- hcvG(0.1, 1, basMet)
    ffclStar <- 1 / (1 + 0.155 * iclStar * fclStar * (hr + hcStar))
    fpclStar <- 1 / (1 + 0.155 * iclStar * fclStar * hcStar / ic)
    imStar <- hcStar * fpclStar / ((hr + hcStar) * ffclStar);
    cres <- 0.0014 * (basMet * metrate) * (34 - tair) # heat transfer to environment by convection in relation to respiration (empirical equation by Fanger)

    pva <- pVapor(tair, phia)
    eres <- 0.0173 * (basMet * metrate) * (5.87 - pva / 1000) # heat removed by evaporation in relation to respiration (empirical equation by Fanger)

    qShiv <- mshiv(tcrSet, tskSet, tcr, tsk)
    vblS <- vblCdilStr(cdil, sigmatr, tcrSet, tskSet, tcr, tsk)# blood flow rate
    # next line is core of Gagge model. If blood flow becomes lower, than skin layer gets more dominant
    alfaSk <- 0.0418 + 0.745 / (vblS + 0.585); # factor to adjust for difference in dominance
    kS <- 5.28 + 1.163 * vblS # conductance between core and skin layer
    Qcr <- (1 - alfaSk) * rMA * cpBody; Qsk <- alfaSk * rMA * cpBody # heat capacity of core and skin layer

    tcrN <- (1 - dtCal * kS / Qcr) * tcr + dtCal / Qcr * (qmet + qShiv - (cres + eres) + kS * tsk) # core temperature at time step before step calculated
    psks <- pVapor(tsk, 100)# saturated water vapour pressure at skin surface calculated with pure water , i.e. adaptive processes such as less salt in sweat might be put heresee also p. 281 of book Shukuya
    emax <- fpcl * (fcl * lr * hc) * (psks - pva) # max rate of water dispersion when the whole skin surface is wet
    tb <- alfaSk * tsk + (1 - alfaSk) * tcr; # average body temperature using calculated tsk and tcr
    tbSet <- alfaSk * tsk + (1 - alfaSk) * tcrSet # average body temperature using set-point values for tsk and tcr
    ersw <- csw * (tb - tbSet) * exp((tsk - tskSet) / 10.7) # amount of sweat secretion

    w <- 0.06 + 0.94 * ersw / emax
    if (w < 0.06){
      w <- 0.06
    }
    if (1 < w){
      w <- 1
    }

    DTQ <- dtCal / Qsk
    tskN <- (1 - DTQ * kS - DTQ * fcl * ffcl * (hr + hc)) * tsk - DTQ * w * fcl * lr * hc * fpcl * psks + DTQ * (kS * tcr + fcl * (hr + hc) * ffcl * top + w * fcl * lr * hc * fpcl * pva) # tsk in next step
    tclN <- ((1 / rcl) * tskN + fcl * hr * tmr + fcl * hc * tair) / (1 / rcl + fcl * hr + fcl * hc)

    ## Exergy balance

    ## Thermal exergy generation by metabolism

    tkcrN <- tcrN + tko; tkskN <- tskN + tko; tkclN <- tclN + tko
    metaEnergy <- qmet + qShiv
    xMet <- metaEnergy * (1 - tkoEnv / tkcrN);
    xwc <- wcXCheck(tkcrN, tkoEnv)

    ## Inhaled humid air

    Vin <- 1.2 * 10  ^  (-6) * metaEnergy # Volume of air intake [V/s]
    cpav <- cpa * (mAir / (rGas * TKa)) * (PO - pva) + cpv * (mWater / (rGas * TKa)) * pva
    xInhaleWc <- Vin * wcEx(cpav, TKa, tkoEnv); xwc <- wcXCheck(TKa, tkoEnv)

    xInhaleWd <- Vin * wdEx(TKa, tkoEnv, pva, pvEnv); xwd <- wdXCheck(pva, pvEnv)

    ## Liquid water generated in the core by metabolism to be dispersed into the exhaled air

    VwCore <- Vin * (0.029 - 0.049 * 10  ^  (-4) * pva)
    xLwCoreWc <- VwCore * Roa * wcEx(cpw, tkcrN, tkoEnv); xwc <- wcXCheck(tkcrN, tkoEnv)

    pvs_env <- pVapor(toEnv, 100)
    xLwCoreWet <- VwCore * Roa * rGas / mWater * tkoEnv * log(pvs_env / pvEnv)

    ## Liquid water generated in the shell by metabolism to be secreted as sweat

    vwShellRow <- w * emax / (2450 * 1000)
    xLwShellWc <- vwShellRow * wcEx(cpw, tkskN, tkoEnv); xwc <- wcXCheck(tkskN, tkoEnv)

    xLwShellWd <- vwShellRow * wdExLw(tkoEnv, pvs_env, pva, pvEnv); xwd <- wdXCheck(pva, pvEnv)

    ## radiant exergy input

    xInRad <- fcl * hr * (tkmr - tkoEnv)  ^  2 / (tkmr + tkoEnv); xwc <- wcXCheck(tkmr, tkoEnv)

    ## total exergy input

    xIntotal <- xMet + xInhaleWc + xInhaleWd + xLwCoreWc + xLwCoreWet + xLwShellWc + xLwShellWd + xInRad

    ## Exergy stored

    xStCore <- Qcr * (1 - tkoEnv / tkcrN) * (tcrN - tcr) / dtCal
    xStShell <- Qsk * (1 - tkoEnv / tkskN) * (tskN - tsk) / dtCal

    ## Exhaled humid air

    pvssCr <- pVapor(tcrN, 100)
    cpav <- cpa * (mAir / (rGas * tkcrN)) * (PO - pvssCr) + cpv * (mWater / (rGas * tkcrN)) * pvssCr
    xExhaleWc <- Vin * wcEx(cpav, tkcrN, tkoEnv); xwc <- wcXCheck(tkcrN, tkoEnv)

    xExhaleWd <- Vin * wdEx(tkcrN, tkoEnv, pvssCr, pvEnv); xwd <- wdXCheck(pvssCr, pvEnv)

    ## water vapor originating from the sweat and humid air containing the evaporated sweat

    xSweatWc <- vwShellRow * wcEx(cpv, tkclN, tkoEnv); xwc <- wcXCheck(tkclN, tkoEnv)

    xSweatWd <- vwShellRow * wdExLw(tkoEnv, pva, pva, pvEnv); xwd <- wdXCheck(pva, pvEnv)

    ## radiant exergy output

    xOutRad <- fcl * hr * (tkclN - tkoEnv)  ^  2 / (tkclN + tkoEnv); xwc <- wcXCheck(tkclN, tkoEnv)

    ## Exergy transfer by convection

    xOutConv <- fcl * hc * (tkclN - TKa) * (1. - tkoEnv / tkclN); xwc <- wcXCheck(tkclN, tkoEnv)

    xouttotal <- xStCore + xStShell + xExhaleWc + xExhaleWd + xSweatWc + xSweatWd + xOutRad + xOutConv

    xConsumption <- xIntotal - xouttotal

    #etStar <- calcet(top, tair, phia, w, im, 50, im);

    tcr <- tcrN
    tsk <- tskN

    ## Output values
    ## Exergy input
    xInmetu[i] <- signif(xMet, 4) #metabolism
    xInmetwcu[i] <- xMetwc #metabolism warm/cold
    xInAIRwcu[i] <- signif(xInhaleWc, 4)	# inhaled humid air
    xInAIRwcwcu[i] <- xInhaleWcwc
    xInAIRwdu[i] <- signif(xInhaleWd, 4) #
    xInAIRwdwdu[i] <- xInhaleWdwd # wet/dry
    xInLUNGwcu[i] <- signif(xLwCoreWc, 4) # water lung
    xInLUNGwcwcu[i] <- xLwCoreWcwc
    xInLUNGwdu[i] <- signif(xLwCoreWet, 4)
    xInLUNGwdwdu[i] <- xLwCoreWetwd
    xInsheLLwcu[i] <- signif(xLwShellWc, 4) # water sweat
    xInsheLLwcwcu[i] <- xLwShellWcwc
    xInsheLLwdu[i] <- signif(xLwShellWd, 4)
    xInsheLLwdwdu[i] <- xLwShellWdwd
    xInraDu[i] <- signif(xInRad, 4) # radiation in
    xInraDwcu[i] <- xInRadwc
    xIntotaLu[i] <- signif(xIntotal, 4) # totaL exergy input

    ## Exergy output
    xoutstorecoreu[i] <- signif(xStCore, 4) # exergy stored
    xoutstoreshelu[i] <- signif(xStShell, 4)
    xoutaIRwcu[i] <- signif(xExhaleWc, 4) # exhaled humid air
    xoutaIRwcwcu[i] <- xExhaleWcwc
    xoutaIRwdu[i] <- signif(xExhaleWd, 4)
    xoutaIRwdwdu[i] <- xExhaleWdwd
    xoutswEATwcu[i] <- signif(xSweatWc, 4) # water vapour
    xoutswEATwcwcu[i] <- xSweatWcwc
    xoutswEATwdu[i] <- signif(xSweatWd, 4)
    xoutswEATwdwdu[i] <- xSweatWdwd
    xoutraDu[i] <- signif(xOutRad, 4) # radiation out
    xoutraDwcu[i] <- xOutRadwc
    xoutCONVu[i] <- signif(xOutConv, 4) # convection
    xoutCONVwcu[i] <- xOutConvwc
    xouttotaLu[i] <- signif(xouttotal, 4) # totaL exergy out

    ## balance
    xconsu[i] <- signif(xConsumption, 4) # exergy consumPtion total
    tsku[i] <- signif(tsk, 4)
    tcru[i] <- signif(tcr, 4)
    wu[i] <- signif(w, 4)

  }# End for

  #setStar <- calcet(top, tair, phia, w, im, 50, imStar);
  results <- data.frame(

    ## Output values
    #setStar,
    #etStar,
    ## Exergy input
    xInmetu, #metabolism
    xInmetwcu, #metabolism warm/cold
    xInAIRwcu, 	# inhaled humid air
    xInAIRwcwcu,
    xInAIRwdu, #
    xInAIRwdwdu, # wet/dry
    xInLUNGwcu, # water lung
    xInLUNGwcwcu,
    xInLUNGwdu,
    xInLUNGwdwdu,
    xInsheLLwcu, # water sweat
    xInsheLLwcwcu,
    xInsheLLwdu,
    xInsheLLwdwdu,
    xInraDu, # radiation in
    xInraDwcu,
    xIntotaLu, # totaL exergy input

    ## Exergy output
    xoutstorecoreu, # exergy stored core
    xoutstoreshelu, # exergy stored shell
    xoutaIRwcu, # exhaled humid air
    xoutaIRwcwcu,
    xoutaIRwdu,
    xoutaIRwdwdu,
    xoutswEATwcu, # water vapour
    xoutswEATwcwcu,
    xoutswEATwdu,
    xoutswEATwdwdu,
    xoutraDu, # radiation out
    xoutraDwcu,
    xoutCONVu, # convection
    xoutCONVwcu,
    xouttotaLu, # totaL exergy out

    ## balance
    xconsu, # exergy consumPtion total
    tsku, # skin temperature
    tcru, # core temperature
    wu, # skin wettedness
    stringsAsFactors=FALSE
  )

  results
} # end of main program