Shape transition error
76 views
Skip to first unread message
Shane Morris
Apr 14, 2026, 9:37:54 AMApr 14
to NicheMapR
Hi Mike,
I've noticed an error that occurs when an endotherm (as an ellipsoid) transitions from a shape (SHAPE_B) where their percentage of subcutaneous fat is enough to have a non-zero fat thickness (output FAT_THICK) to a shape where their fat thickness becomes 0. It does not occur when subcutaneous fat is not present (SUBQFAT=0). It appears to be an issue with the calculation of the flesh conductivity (K_FLESH). The flesh conductivity does not change when the SHAPE_B parameter changes. Could be something to do with how the calculation of whether FAT_THICK is zero or non-zero occurs in the GEOM_ENDO function which is after the skin conductivity value (AK1_LAST) is calculated in SOLVENDO.f? Below is the code to explain what I mean. It is for a hibernating brown bear.
Cheers,
Shane.
# load libraries
library(NicheMapR); library(tidyverse)
# Input variables for a hibernating bear
TAs <- seq(-30, 30, 1)
AMASS <- 135
Q10_torpor <- 1
QBASAL <- 0.3509*AMASS
TC <- 37.5
TC_MAX <- 40
TC_MIN <- 30
TC_INC <- 0.01
# The variables the issue is related to
FATPCT <- 20
SUBQFAT <- 1
SHAPE_B <- 1.3
SHAPE_B_MAX <- 6
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out1 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out1$TAs <- TAs
# Plot the metabolic rate
ggplot(endo.out1)+ geom_line(aes(x=TAs, y=enbal.QGEN))+ theme_classic()
# Why this peak at 7 degrees?
ggplot(endo.out1)+ geom_line(aes(x=TAs, y=treg.SHAPE_B))+ theme_classic()
ggplot(endo.out1)+ geom_line(aes(x=TAs, y=morph.AREA))+ theme_classic()
# Also why does the surface area decrease when posture changes? The surface area of a more spherical object should be lower than one that is more ellipsoidal. Is it to do with the fur flattening, thus decreasing the surface area?
# This blip in metabolism is tied to subcutaneous fat. Here I will plot a model with no subcutaneous fat, one with 50% extra fat and one with 50% less fat
SUBQFAT <- 0
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out2 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out2$TAs <- TAs
# This blip is clearly tied to the fat percentage. It dramatically rises with more fat and dramatically decreases with
SUBQFAT <- 1
FATPCT <- 30
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out3 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out3$TAs <- TAs
FATPCT <- 10
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out4 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out4$TAs <- TAs
endo.out2 %>% mutate(FAT="None")
# Plot the metabolic rate
(endo.out1 %>% mutate(FAT="Normal") %>% bind_rows(endo.out2 %>% mutate(FAT="None"), endo.out3 %>% mutate(FAT="50% extra"), endo.out4 %>% mutate(FAT="50% less"))) %>%
ggplot()+ geom_line(aes(x=TAs, y=enbal.QGEN, group=FAT, colour=FAT), linewidth=2)+ theme_classic()
# It appears to be due to the transition from ellipsoid ratios 4.5 to 4.6
# If I go back to the original fat variables and keep the shape always at 4 or 5 it does not happen
FATPCT <- 20.1
SHAPE_B <- 4.5
SHAPE_B_MAX <- 4.5
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out5 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out5$TAs <- TAs
SHAPE_B <- 4.5
SHAPE_B_MAX <- 4.6
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out6 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out6$TAs <- TAs
SHAPE_B <- 4.6
SHAPE_B_MAX <- 4.6
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out7 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out7$TAs <- TAs
(endo.out5 %>% mutate(SHAPE="Just 4.5") %>% bind_rows(endo.out6 %>% mutate(SHAPE="4.5 to 4.6"), endo.out7 %>% mutate(SHAPE="Just 4.6"))) %>%
ggplot()+ geom_line(aes(x=TAs, y=enbal.QGEN, group=SHAPE, colour=SHAPE), linewidth=2)+ theme_classic()
# Plotting all variables to look at patterns that may explain
for(i in 1:ncol(endo.out5)){
min5 <- min(endo.out5[, i], na.rm=TRUE)
max5 <- max(endo.out5[, i], na.rm=TRUE)
plot(TAs, endo.out5[, i], type="l", main=paste0("Variable ", colnames(endo.out5)[i]), ylim=c(min5*0.5, max5*1.5))
lines(TAs, endo.out6[, i], col="blue")
lines(TAs, endo.out7[, i], col="red")
abline(v=-4, col="grey")
abline(v=11, col="grey")
}
# From looking at the variables. The only one I can see that is a candidate is the K_FLESH
(endo.out5 %>% mutate(SHAPE="Just 4.5") %>% bind_rows(endo.out6 %>% mutate(SHAPE="4.5 to 4.6"), endo.out7 %>% mutate(SHAPE="Just 4.6"))) %>%
ggplot()+ geom_line(aes(x=TAs, y=treg.K_FLESH, group=SHAPE, colour=SHAPE), linewidth=2)+ theme_classic() +
geom_vline(xintercept=-4, colour="grey")+ geom_vline(xintercept=11, colour="grey")
# it appears to be caused by the delay in the increase in flesh conductivity
# Can see here that the transition in SHAPE_B occurs but the skin conductivity doesn't
endo.out6 %>%
ggplot()+ geom_line(aes(x=TAs, y=treg.SHAPE_B), linewidth=2)+ theme_classic() +
geom_vline(xintercept=-4, colour="grey")+ geom_vline(xintercept=11, colour="grey")
library(NicheMapR); library(tidyverse)
# Input variables for a hibernating bear
TAs <- seq(-30, 30, 1)
AMASS <- 135
Q10_torpor <- 1
QBASAL <- 0.3509*AMASS
TC <- 37.5
TC_MAX <- 40
TC_MIN <- 30
TC_INC <- 0.01
# The variables the issue is related to
FATPCT <- 20
SUBQFAT <- 1
SHAPE_B <- 1.3
SHAPE_B_MAX <- 6
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out1 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out1$TAs <- TAs
# Plot the metabolic rate
ggplot(endo.out1)+ geom_line(aes(x=TAs, y=enbal.QGEN))+ theme_classic()
# Why this peak at 7 degrees?
ggplot(endo.out1)+ geom_line(aes(x=TAs, y=treg.SHAPE_B))+ theme_classic()
ggplot(endo.out1)+ geom_line(aes(x=TAs, y=morph.AREA))+ theme_classic()
# Also why does the surface area decrease when posture changes? The surface area of a more spherical object should be lower than one that is more ellipsoidal. Is it to do with the fur flattening, thus decreasing the surface area?
# This blip in metabolism is tied to subcutaneous fat. Here I will plot a model with no subcutaneous fat, one with 50% extra fat and one with 50% less fat
SUBQFAT <- 0
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out2 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out2$TAs <- TAs
# This blip is clearly tied to the fat percentage. It dramatically rises with more fat and dramatically decreases with
SUBQFAT <- 1
FATPCT <- 30
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out3 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out3$TAs <- TAs
FATPCT <- 10
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out4 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out4$TAs <- TAs
endo.out2 %>% mutate(FAT="None")
# Plot the metabolic rate
(endo.out1 %>% mutate(FAT="Normal") %>% bind_rows(endo.out2 %>% mutate(FAT="None"), endo.out3 %>% mutate(FAT="50% extra"), endo.out4 %>% mutate(FAT="50% less"))) %>%
ggplot()+ geom_line(aes(x=TAs, y=enbal.QGEN, group=FAT, colour=FAT), linewidth=2)+ theme_classic()
# It appears to be due to the transition from ellipsoid ratios 4.5 to 4.6
# If I go back to the original fat variables and keep the shape always at 4 or 5 it does not happen
FATPCT <- 20.1
SHAPE_B <- 4.5
SHAPE_B_MAX <- 4.5
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out5 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out5$TAs <- TAs
SHAPE_B <- 4.5
SHAPE_B_MAX <- 4.6
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out6 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out6$TAs <- TAs
SHAPE_B <- 4.6
SHAPE_B_MAX <- 4.6
ptm <- proc.time()
endo.out <- lapply(1:length(TAs), function(x) {
endoR(TA = TAs[x], VEL = 0.1, RH = 5, QSOLR = 0,
AMASS = AMASS, Q10 = Q10_torpor, QBASAL = QBASAL,
TC = TC, TC_MAX = TC_MAX, TC_MIN = TC_MIN,
SHAPE_B = SHAPE_B,SHAPE_B_MAX = SHAPE_B_MAX,
FATPCT = FATPCT, SUBQFAT = SUBQFAT,
TORPOR = 1,
SHAPE=4, SAMODE=2, DHAIRD = 0.000076, DHAIRV = 0.000076, LHAIRD = 0.0759,
LHAIRV = 0.0828,ZFURD = 0.0559,ZFURV = 0.0628, RHOD = 8.5e+6, RHOV = 8.5e+6*0.75,
REFLD = 0.430, REFLV = 0.29, PZFUR = 0.1, PVEN = 0.25,UNCURL = 0.1,
PCTWET = 0.1, PCTWET_MAX = 10,PCTWET_INC = 0.25, DELTAR = 5,EXTREF = 25,
PANT_INC = 0.1, PANT_MAX = 10,PANT_MULT = 1,AK1 = 1.5,AK2 = 0.23,AK1_MAX = 2.8,
AK1_INC = 0.1,ANDENS = 1033.6, RQ = 0.7 )
})
proc.time() - ptm
# combine data frame and add ambient temperature
endo.out7 <- do.call("rbind", lapply(endo.out, data.frame))
endo.out7$TAs <- TAs
(endo.out5 %>% mutate(SHAPE="Just 4.5") %>% bind_rows(endo.out6 %>% mutate(SHAPE="4.5 to 4.6"), endo.out7 %>% mutate(SHAPE="Just 4.6"))) %>%
ggplot()+ geom_line(aes(x=TAs, y=enbal.QGEN, group=SHAPE, colour=SHAPE), linewidth=2)+ theme_classic()
# Plotting all variables to look at patterns that may explain
for(i in 1:ncol(endo.out5)){
min5 <- min(endo.out5[, i], na.rm=TRUE)
max5 <- max(endo.out5[, i], na.rm=TRUE)
plot(TAs, endo.out5[, i], type="l", main=paste0("Variable ", colnames(endo.out5)[i]), ylim=c(min5*0.5, max5*1.5))
lines(TAs, endo.out6[, i], col="blue")
lines(TAs, endo.out7[, i], col="red")
abline(v=-4, col="grey")
abline(v=11, col="grey")
}
# From looking at the variables. The only one I can see that is a candidate is the K_FLESH
(endo.out5 %>% mutate(SHAPE="Just 4.5") %>% bind_rows(endo.out6 %>% mutate(SHAPE="4.5 to 4.6"), endo.out7 %>% mutate(SHAPE="Just 4.6"))) %>%
ggplot()+ geom_line(aes(x=TAs, y=treg.K_FLESH, group=SHAPE, colour=SHAPE), linewidth=2)+ theme_classic() +
geom_vline(xintercept=-4, colour="grey")+ geom_vline(xintercept=11, colour="grey")
# it appears to be caused by the delay in the increase in flesh conductivity
# Can see here that the transition in SHAPE_B occurs but the skin conductivity doesn't
endo.out6 %>%
ggplot()+ geom_line(aes(x=TAs, y=treg.SHAPE_B), linewidth=2)+ theme_classic() +
geom_vline(xintercept=-4, colour="grey")+ geom_vline(xintercept=11, colour="grey")
NicheMapR
Apr 16, 2026, 4:28:15 AMApr 16
to NicheMapR
Hi Shane,
Interesting issue - the calculation of the fat layer for an ellipsoid is quite a complex one. I tried the same code for a cylinder and the issue didn't happen. I'll dig into it soon and let you know.
All the best,
Mike
Mike
Shane Morris
Apr 21, 2026, 5:04:02 AMApr 21
to NicheMapR
Hi Mike,
Thanks, that would be great. Yes, it is much more complex than I initially thought. I was not aware that the fat thickness changed so much as the animal "stretches out".
Cheers,
Shane.
Shane Morris
Jun 11, 2026, 7:47:55 AMJun 11
to NicheMapR
Hi Mike,
Have you had a chance to look at this? Is there anyway to simply include a minimum fat thickness to avoid the fat from disappearing? As you can the fat thickness should likely be set to a minimum of 0.02 - 0.025.
Cheers,
Shane.
NicheMapR
Jun 16, 2026, 7:43:30 AMJun 16
to NicheMapR
Hi Shane,
Sorry for the delay in responding - I found a work-around (with the help of Claude code) - now there's an iteration using Newton's method to solve for the fat thickness that doesn't miss the solution. It's pushed to github.
All the best,
Mike
p.s. with the fat layer, for the HomoTherm human model I realised that the subcutaneous fat is bypassed as peripheral blood flow increases and I implemented that for the human model but not yet for the endotherm model, so this means that having a fat layer will make the animal unrealistically vulnerable to high temeprature
Shane Morris
Jun 17, 2026, 6:16:49 AMJun 17
to NicheMapR
Hi Mike,
Thanks for clearing that up. One additional thing I noticed that is unrelated to the SHAPE issue is that in the endoR you describe the PANT_MULT as the multiplier on basal metabolic rate at maximum panting level but in the FORTRAN code you calculate it as
PANT_COST=((PANT-1.)/(PANT_MAX+1e-6-1.)*(PANT_MULT-1.)*QBASAL_REF)
which means the cost is incurred once PANT exceeds 1. Was this intentional and the description is wrong or vice versa? Just thought I should flag it in case it was not what was intended.
Thanks again for sorting the SHAPE issue,
Shane
NicheMapR
Jun 18, 2026, 4:16:04 PM (13 days ago) Jun 18
to NicheMapR
Hi Shane,
It's what is intended in that by the time PANT has reached PANT_MAX, the "basal" (really, "minimum") metabolic rate has risen by a factor of PANT_MULT, but it rises linearly along the way. E.g.:
If QBASAL_REF is 10, PANT_MAX is 5 and PANT_MULT is 1.1, at PANT of 2, PANT_COST is 0.25 and the new QBASAL is 10.25. If PANT is PANT_MAX then PANT_COST is 1 and QBASAL is 11, so it's been multiplied by 1.1 = PANT_MULT.
One issue is that as the core temperature rises, the QBASAL_REF goes up with the Q10 effect so that needs to be kept in mind.
All the best,
Mike
Mike
Shane Morris
Jun 19, 2026, 9:27:20 AM (13 days ago) Jun 19
to NicheMapR
Hi Mike,
That makes sense. Thanks again.
All the best,
Shane.
Reply all
Reply to author
Forward
0 new messages