#### R CODE to "Biophysical models unravel associations between glucocorticoids and thermogenic metabolic rates across avian species"
#### J.G. Rubalcaba & B. Jimeno

#### The function endoTe() computes thermogenic metabolic rate as a function of ambient temperature, solar radiation, and wind velocity, based on the 
# furry ellipsoid model by Porter and Kearney (2009, PNAS). See manuscript & supplementary information for details (Rubalcaba & Jimeno 2021, Funct Ecol)


endoTe <- function(Tc,     # Core temperature (ºC)
                   M,      # Body mass (g)
                   BMR=NA, # Basal metabolic rate (W) - if NA, the function estimates BMR from body mass
                   shape,  # Shape factor (1: spherical, >1: ellipsoidal)
                   F_a,    # View factor (exposition to solar radiation, 0-1)
                   rf,     # Insulation layer depth (m)
                   abs,    # Absorbance to short-wave radiation (0-1)
                   emis,   # Emissivity of long-wave radiation (0-1)
                   Ta,     # Air temperature (ºC)
                   v,      # Wind speed (ms-1)
                   RH,     # Relative humidity (%)
                   S)      # Solar radiation (direct + scattered) (Wm-2)
{
  ## Basal Metabolic Rate ##
  if(is.na(BMR)){
    BMR_O2 = 6.7141 * M^0.6452 ## Fristoe et al. 2015; Someveille et al. 2018 ## Mammals: exp(1.3621)*M^0.7016 (pantheria)
    BMR = BMR_O2 / 172 # (W) 1W = 172 mLO2 h-1
  }
 
  ## Ellipsoid geometry (Porter and Kearney 2009) ## 
  k = 3
  b = (k*(M*1e-6/(4*pi*shape)))^(1/3)
  c = b
  a = shape*b
  S2 = a^2*b^2*c^2 / (a^2*b^2 + a^2*c^2 + b^2*c^2) 
  V = 4/3*pi*a*b*c # Body volume (m3)
  eccentric = sqrt(a^2-c^2)/a
  Area = 2*pi*b^2 + 2*pi*((a*b)/eccentric)*asin(eccentric) # Total surface area (m2)
  
  ## Insulation layer ## 
  kins = 0.027  # Thermal conductivity feathers (W m-1 ºC-1) - Porter and Kearney 2009
  ao = a + rf   # inner radium + feather layer (m)
  bo = b + rf
  co = c + rf
  Vo = 4/3*pi*ao*bo*co  # External volume (m3)
  ecc_out = sqrt(ao^2-co^2)/ao
  Ao = 2*pi*bo^2 + 2*pi*((ao*bo)/ecc_out)*asin(ecc_out) # External surface (m2)
  
  ## Thermal radiation ## 
  sigma = 5.670367e-8  # Stefan-Boltzmann constant (W m-2 K-4)
  R = 4 * Ao * sigma * emis * (Ta+273)^3 # Thermal radiation coefficient
  
  ## Convection ## 
  L = V^(1/3)   # Characteristic length (m) (Mitchell 1976)
  kb = 0.5 + 6.14*b + 0.439 # Thermal conductivity body (W m^-1 ºC-1) Porter and Kearney (2009)
  kf =  1.5207e-11*(Ta+273)^3 - 4.8574e-08*(Ta+273)^2 + 1.0184e-04*(Ta+273) - 3.9333e-04 # Thermal conductivity of air (W m-1 ºC-1)
  Pd = 101325   # Atmosferic pressure (Pa) 
  rho_da = Pd / (287.05 * (Ta+273))  # Density of dry air (kg m-3)
  Pv = exp(77.3450 + 0.0057 * (Ta+273) - 7235 / (Ta+273)) / (Ta+273)^8.2  # Water vapor pressure (Pa)
  RH = RH/100 # ransform relative humidity % to proportion 0-1
  x = RH * 0.62198 * Pv / Pd # Humidity ratio (mass_air_water / mas_dry_air)
  rho = rho_da * (1 + x) / (1 + 1.609 * x) # Density of moist air (Kg m-3)
  mu = 1.458e-6 * (Ta+273)^(3/2) / ((Ta+273) + 110.4)  # Air dynamic viscosity (Pa s)
  cp = 1000*(-3.46607e-20*Ta^6+9.121847e-17*Ta^5+1.079642e-13*Ta^4-5.714484e-10*Ta^3+5.773354e-07*Ta^2+8.976385e-06*Ta+1.00528623891845) # Specific heat capacity air (J Kg-1 ºC-1)
  
  Re = rho * v * L / mu     # Reynolds number
  Pr = mu * cp / kf         # Prandtl number
  Nu_forced = 2.0 + 0.6 * Re^(1/2) * (Pr)^(1/3) # Nusselt number forced convection
  
  Ta[which(Ta > Tc)] <- Tc # correction to avoid negative values in Ta - Tc below (let Gr be 0 when Ta > Tc)
  Gr = rho^2 * 1/(Ta+273) * 9.8 * L^3 * (Tc - Ta) / mu^2  # Grashof number
  Nu_free = 2.0 + 0.6 * Gr^(1/4) * Pr^(1/3) # Nusselt number free convection
  
  Nu_total = (Nu_free^3 + Nu_forced^3)^(1/3) 
  hc = Nu_total * kf/L # Convection coefficient
  
  ### Operative temperature ##
  Te = Ta + F_a * abs * S / (R + hc)
  
  ## Metabolic Heat Production ## 
  Rbody = S2 / (2*V*kb) # Body resistance to heat conduction
  Rconv = 1 / (hc*Ao)   # Convective heat transfer resistance
  Rir = 1 / R           # IR loss resistance
  Rins = (bo - b) / (Ao * kins) # Resistance of insulation layer 
  Rtot = Rbody + Rins + (Rconv * Rir) / (Rconv + Rir) # Total resistance
  
  Q_gen = (Tc - Te) / Rtot # Thermogenic heat (W)

  MR = Q_gen
  MR[MR < BMR] = BMR # Whole organism metabolic rate

  ## Model Output dataframe ##
  output <- data.frame(MR,  # Modeled thermogenic metabolic rate
                       BMR, # Estimated (or user-specified) BMR
                       Area,# Modeled skin surface area
                       V,   # body volume
                       Te,  # Modeled operative temperature
                       Ta,  # user-specified meteorological conditions (air temperature, wind velocity, relative humidity, and solar radiation)
                       v, 
                       RH, 
                       S)
  return(output)
}

################ Zebra Finch (Briga & Verhulst 2017)

Tc_Tgut = 40.3  # Core temperature (ºC) (Whittow 1986 in Avian Physiology)
M = 16          # Body mass (g)
BMR_Tgut = 0.22 # Basal metabolic rate (Briga & Verhulst 2017)

Ta <- seq(5, 40, length.out = 100) # Range of chamber temperatures (ºC)
out <- endoTe(Tc=Tc_Tgut, M, BMR=BMR_Tgut, shape=1.1, F_a=0.5, rf=0.005, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)
out_fluff <- endoTe(Tc=Tc_Tgut, M, BMR=BMR_Tgut, shape=1.5, F_a=0.5, rf=0.01, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)

plot(MR ~ Ta, out, ylab = "Metabolic rate (W)", xlab = "Air Temperature (ºC)", ylim=c(0,1),
     cex.axis=1.2, cex.lab=1.5, type='l', lwd=1, lty=2)
lines(MR ~ Ta, out_fluff)

################ House Sparrow (Hudson and Kimzey 1966, nocturnal MR Houston population)

Tc_Pdom = 38.6
M_Pdom = 33
BMR_Pdom = 2.34 * M_Pdom / 172

Ta <- seq(-5, 40, length.out = 100) 
out <- endoTe(Tc=Tc_Pdom, M=M_Pdom, BMR=BMR_Pdom, shape=1.5, F_a=0.5, rf=0.005, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)
out_fluff <- endoTe(Tc=Tc_Pdom, M=M_Pdom, BMR=BMR_Pdom, shape=1.5, F_a=0.5, rf=0.01, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)

plot(MR ~ Ta, out, lwd=1, lty=2, type='l', ylab = "Metabolic rate (W)", xlab = "Air Temperature (ºC)", pch=20, cex=1, ylim=c(0,1.5),
     cex.axis=1.2, cex.lab=1.5)
lines(MR ~ Ta, out_fluff)

################ Alaskan Ptarmigan (West 1972)
# Lagopus lagopus

Mlagopus = 564.5
Tclagoupus = 39.7
BMRlagopus = 699 / 172 # summer

Ta <- seq(-60, 40, length.out = 100)
out <- endoTe(Tc=Tclagoupus, M=Mlagopus, BMR=BMRlagopus, shape=1.5, F_a=0.5, rf=0.005, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)
out_fluff <- endoTe(Tc=Tclagoupus, M=Mlagopus, BMR=BMRlagopus, shape=1.5, F_a=0.5, rf=0.01, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)

plot(MR ~ Ta, out, lwd=1, lty=2, type='l', ylab = "Metabolic rate (W)", xlab = "Air Temperature (ºC)", 
     pch=20, cex=1, ylim=c(0,12), cex.axis=1.2, cex.lab=1.5)
lines(MR ~ Ta, out_fluff)

BMRlagopus = 625 / 172 # winter

out <- endoTe(Tc=Tclagoupus, M=Mlagopus, BMR=BMRlagopus, shape=1.5, F_a=0.5, rf=0.01, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)
out_fluff <- endoTe(Tc=Tclagoupus, M=Mlagopus, BMR=BMRlagopus, shape=1.5, F_a=0.5, rf=0.02, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)

plot(MR ~ Ta, out, lwd=1, lty=2, type='l', ylab = "Metabolic rate (W)", xlab = "Air Temperature (ºC)", 
     pch=20, cex=1, ylim=c(0,12), cex.axis=1.2, cex.lab=1.5)
lines(MR ~ Ta, out_fluff)

# Lagopus mutus
Mmutus = 432
Tcmutus = 41
BMRmutus = 628 / 172

Ta <- seq(-60, 40, length.out = 100)
out <- endoTe(Tc=Tcmutus, M=Mmutus, BMR=BMRmutus, shape=1.5, F_a=0.5, rf=0.01, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)
out_fluff <- endoTe(Tc=Tcmutus, M=Mmutus, BMR=BMRmutus, shape=1.5, F_a=0.5, rf=0.02, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)

plot(MR ~ Ta, out, lwd=1, lty=2, type='l', ylab = "Metabolic rate (W)", xlab = "Air Temperature (ºC)", pch=20, cex=1, ylim=c(0,12), cex.axis=1.2, cex.lab=1.5)
lines(MR ~ Ta, out_fluff)

################ Rodents Liu et al (2004)
# Clethrionomys rufocanus

Tc_Crufo <- 38.2
M_Crufo <- 23.05
BMR_Crufo <- M_Crufo * 3.62 / 172

Ta <- seq(5, 40, length.out = 100)
out <- endoTe(Tc=Tc_Crufo, M=M_Crufo, BMR=BMR_Crufo, shape=1.5, F_a=0.5, rf=0.002, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)

out$MR <- out$MR 

plot(MR ~ Ta, out, ylab = "Metabolic rate (W)", xlab = "Air Temperature (ºC)", pch=20, cex=1, ylim=c(0,1.5), cex.axis=1.2, cex.lab=1.5)

# Apodemus speciosus
Tc_Aspe <- 36.9
M_Aspe <- 28.52
BMR_Aspe <- M_Aspe * 2.69 / 172

Ta <- seq(5, 40, length.out = 100)
out <- endoTe(Tc=Tc_Aspe, M=M_Aspe, BMR=BMR_Aspe, shape=2.5, F_a=0.5, rf=0.003, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)

plot(MR ~ Ta, out, ylab = "Metabolic rate (W)", xlab = "Air Temperature (ºC)", pch=20, cex=1, ylim=c(0,1.5), cex.axis=1.2, cex.lab=1.5)


# Apodemus agrarius
Tc_Aagra <- 36.1
M_Aagra <- 24.37
BMR_Aagra <- M_Aagra * 3.29 / 172

Ta <- seq(5, 40, length.out = 100)
out <- endoTe(Tc=Tc_Aagra, M=M_Aagra, BMR=BMR_Aagra, shape=2.5, F_a=0.5, rf=0.002, abs=0.3, emis=0.95, Ta, v=0.001, RH=50, S=0)

plot(MR ~ Ta, out, ylab = "Metabolic rate (W)", xlab = "Air Temperature (ºC)", pch=20, cex=1, ylim=c(0,1.5), cex.axis=1.2, cex.lab=1.5)

################ Solar radiation I: Wolf & Walsberg 1996 (Auriparus flaviceps)

Tc = 40 
F_a = 1
v <- seq(0, 3, length.out = 100)
  
out_sun_fluff <- endoTe(Tc, M=7, BMR=NA, shape=1.5, F_a, rf=0.01, abs=0.3, emis=0.95, Ta=15, v, RH=50, S=500)
out_sun <- endoTe(Tc, M=7, BMR=NA, shape=1.5, F_a, rf=0.005, abs=0.3, emis=0.95, Ta=15, v, RH=50, S=500)
out_nosun_fluff <- endoTe(Tc, M=7, BMR=NA, shape=1.5, F_a, rf=0.01, abs=0.3, emis=0.95, Ta=15, v, RH=50, S=0)
out_nosun <- endoTe(Tc, M=7, BMR=NA, shape=1.5, F_a, rf=0.005, abs=0.3, emis=0.95, Ta=15, v, RH=50, S=0)

plot(MR ~ v, out_sun, ylab = "Metabolic rate (W)", xlab = "Wind speed (m / s)", ylim =c(0.1, 0.5), type="l", lty=2, cex.axis=1.2, cex.lab=1.5)
lines(MR ~ v, out_sun_fluff)
lines(MR ~ v, out_nosun, lty=2)
lines(MR ~ v, out_nosun_fluff)

################ Solar radiation II: Wolf & Walsberg 1995 (Ground squirrels)

m_Slat = 202
m_Ssat = 332

v <- seq(0, 4, length.out = 100)

out_Slateralis <- endoTe(Tc=37, M=m_Slat, BMR=NA, shape=1.5, F_a=0.6, rf=0.001, abs=0.3, emis=0.95, Ta=15, v, RH=50, S=780)
out_Slateralisnosun <- endoTe(Tc=37, M=m_Slat, BMR=NA, shape=1.5, F_a=0.6, rf=0.001, abs=0.3, emis=0.95, Ta=15, v, RH=50, S=0)

plot(MR ~ v, out_Slateralis, ylab = "Metabolic rate (W)", xlab = "Wind speed (m / s)", type="l", ylim =c(0, 5), cex.axis=1.2, cex.lab=1.5)
lines(MR ~ v, out_Slateralisnosun, pch=20)

# 

out_Ssaturatus <- endoTe(Tc=37, M=m_Ssat, BMR=NA, shape=1.5, F_a=0.33, rf=0.0015, abs=0.3, emis=0.95, Ta=15, v, RH=50, S=780)
out_Ssaturatusnosun <- endoTe(Tc=37, M=m_Ssat,BMR=NA,shape=1.5, F_a=0.33, rf=0.0015, abs=0.3, emis=0.95, Ta=15, v, RH=50, S=0)

plot(MR ~ v, out_Ssaturatus, ylab = "Metabolic rate (W)", xlab = "Wind speed (m / s)", type="l", ylim =c(0, 8), cex.axis=1.2, cex.lab=1.5)
lines(MR ~ v, out_Ssaturatusnosun, pch=20)
##