Deconstructed Beetles Original Code
Authors removed for double blind review
16/02/2022
Introduction
In this file we provide the code used for the analysis of the data from the paper “Deconstructed beetles: green composite materials with potential for passive cooling due to high near-infrared reflectivity”.
In our manuscript, we describe the optical properties and underlying structures that produce very high near-infrared reflectance in three species of green scarab beetles Xylonichus eucalypti, Anoplognathus prasinus and Paraschizognathus olivaceus. Although the chiral structures in scarabs have been well described in previous literature, we have limited knowledge of the natural mechanisms to manipulate near-infrared light (NIR; 700 – 1700 nm). In addition, previous studies in light manipulation by the beetle elytra have been focused mainly on the cuticle. Thus, our aim was to explore if other elements in the elytra can contribute to their characteristic visible and NIR reflectance.
Three green scarabs with High NIR reflectance
The original .Rmd file and raw data can be accessed here: https://bit.ly/3bmwGC4 (This link contains information that may disclose author’s identity)
Setting up
Libraries
library(ggplot2)
library(wesanderson)
library(gridExtra)
library(rptR)
library(cowplot)
library(pavo)
library(tidyverse)
library(dplyr)
library(kableExtra)
library(graphics)
library(matrixStats)Data Sets
Reflect<-read.csv("../Data/1_HemisphericalReflectance.csv") %>%
as.rspec(.) %>%
procspec(.,fixneg = "zero")
Trans<-read.csv("../Data/2_Transmittance.csv") %>%
filter (wl <= 1700 & wl >= 400) %>%
as.rspec(.) %>%
procspec(.,fixneg = "zero")
Transmicro <- read.csv("../Data/3_TransmittanceMicroScale.csv")
Cylinders <- read.csv("../Data/4_XylonichusCylinderDiameters.csv")
Centroids <- read.csv("../Data/5_XylonichusCentroidDistance.csv")
QuasiOrder <- read.csv("../Data/6_PrumModelResults.csv")
Ordered <- read.csv("../Data/7_PhotonicCrystalModelResults.csv")
Scenarios <- read.csv("../Data/8_PhotonicCrystalScenarios.csv")Optical properties (macro)
Comparison between the optical properties at a macro scale for the three beetles. (Figure 2 in the manuscript)
P_Prsi<-ggplot() +
geom_ribbon(data=Trans,
mapping = aes(x=wl,y=100-prsi3d,
ymin=100,ymax=100-prsi3d),
fill="#b7b7b7", colour="black", alpha=0.5)+ # Transmittance
geom_ribbon(data=Trans,
mapping = aes(x=wl,y=100-prsi3d,
ymin=Reflect$prsi3,ymax=100-prsi3d),
colour="black", fill="#2b8fff", alpha=0.5)+ # Absorbance
geom_area(data=Reflect, aes(x=wl, y=prsi3),
fill="#f8f813", colour="black", alpha=0.5) + # Reflectance
scale_x_continuous(limits=c(400, 1700),
expand = c(0, 0))+
scale_y_continuous(limits=c(0, 100), breaks = c(50,100),
expand = c(0, 0),
sec.axis = sec_axis(~100-.,breaks = c(100,50))) +
theme_classic()+
theme(axis.title.x = element_blank(),
axis.text.x = element_text(color="white"))+
theme(axis.text.y = element_text(size =18)) +
ylab( " " )+
theme(plot.margin = unit(c(0.6,0.4,0.1,0.4), "cm"))+
geom_vline(xintercept = 700, linetype="dotted",
color = "black", size=1)
P_Xyln<-ggplot() +
geom_ribbon(data=Trans,
mapping = aes(x=wl,y=100-xyln3d,
ymin=100,ymax=100-xyln3d),
fill="#b7b7b7", colour="black", alpha=0.5)+ # Transmittance
geom_ribbon(data=Trans,
mapping = aes(x=wl,y=100-xyln3d,
ymin=Reflect$xyln3,ymax=100-xyln3d),
colour="black", fill="#2b8fff", alpha=0.5)+ # Absorbance
geom_area(data=Reflect, aes(x=wl, y=xyln3),
fill="#f8f813", colour="black", alpha=0.5) + # Reflectance
scale_x_continuous(limits=c(400, 1700),
expand = c(0, 0))+
scale_y_continuous(limits=c(0, 100), breaks = c(50,100),
expand = c(0, 0),
sec.axis = sec_axis(~100-.,breaks = c(100,50))) +
theme_classic()+
theme(axis.title.x = element_blank(),
axis.text.x = element_text(color="white"))+
theme(axis.text.y = element_text(size =18)) +
ylab( " " )+
theme(plot.margin = unit(c(0.15,0.4,0.15,0.4), "cm"))+
geom_vline(xintercept = 700, linetype="dotted",
color = "black", size=1)
P_Ocul<-ggplot() +
geom_ribbon(data=Trans,
mapping = aes(x=wl,y=100-oliv2d,
ymin=100,ymax=100-oliv2d),
fill="#b7b7b7", colour="black", alpha=0.5)+ # Transmittance
geom_ribbon(data=Trans,
mapping = aes(x=wl,y=100-oliv2d,
ymin=Reflect$oliv2,ymax=100-oliv2d),
colour="black", fill="#2b8fff", alpha=0.5)+ # Absorbance
geom_area(data=Reflect, aes(x=wl, y=oliv2),
fill="#f8f813", colour="black", alpha=0.5) + # Reflectance
scale_x_continuous(limits=c(400, 1700),
expand = c(0, 0))+
scale_y_continuous(limits=c(0, 100), breaks = c(50,100),
expand = c(0, 0),
sec.axis = sec_axis(~100-.,breaks = c(100,50))) +
theme_classic()+
theme(axis.title.x = element_blank())+
theme(axis.text.y = element_text(size =18)) +
theme(axis.text.x = element_text(size =18)) +
ylab( " " )+
theme(plot.margin = unit(c(0.1,0.4,0.3,0.4), "cm"))+
geom_vline(xintercept = 700, linetype="dotted",
color = "black", size=1)# Align the three plots:
plot_grid(P_Prsi, P_Xyln, P_Ocul,
align = "v", nrow=3, ncol=1, rel_widths = c(1, 1))
Yellow = Reflectance grey = Transmittance (inverted axis) blue = Absorbance
Elytra Architecture
We observed the cross section of the beetle elytra with SEM (Figure 3). Images were analyzed with the software ImageJ. Results from these measurements are in Data/OtherSoftware and in table 1 of the manuscript.
Scatterometry
Reflectance
We measured the reflectanc eof the composite material (cuticle + white underlay) in two configurations (dorsal and ventral). Data analyzed with a costumized code in MatLab. Raw files in the folder Data/OtherSoftware. (Figure 4 in the manuscript top and middle row)
Transmittance
The three beetles have very similar transmittance profile in their cuticle. (Figure 4 bottom row)
# Prepare for plot
Transmicro<-
Transmicro %>%
mutate(across(spp, factor, levels=c("A.pras","X.euca","P.ocul")))# plot
ggplot(data=Transmicro, aes(wl,Transm,spp))+
geom_line()+
facet_wrap(~spp)+
theme_bw()+
theme(
strip.background = element_blank(),
strip.text.x = element_blank()
)
Nanostructures
Morphological analysis: Statistical Descriptors of the particles in the white underlay of Xylonichus eucalypti
Calculations
Cylinder diameters
Cylinders <-
Cylinders %>%
mutate(MajorAxisnm = MajorAxis * 1000, # from microns to nm
MinorAxisnm = MinorAxis * 1000) # from microns to nm
# Prepare for the plot
HistAxis<-
Cylinders %>%
dplyr::select (MajorAxisnm, MinorAxisnm ) %>%
gather (Axis, Length)
# Plot
ggplot(HistAxis, aes(x=Length, fill=Axis, color=Axis))+
geom_histogram(fill="gray", alpha=0.3,
position="identity", binwidth=25)+
theme_bw()+
theme(legend.position="top")+
ylab("Counts")+
xlab("Cylinder Diameter (nm)")+
theme(legend.position = c(0.8,0.5))+
theme(legend.title = element_blank())+
theme(legend.text = element_text(size=7))+
scale_color_manual(labels = c("Major axis", "Minor axis"),
values=c("#ffd781", "#93c7ff")) +
theme(legend.key.size = unit(0.4, 'cm'))
Span of the diameters:
quantile(Cylinders$MajorAxisnm, c(0.05,0.95))## 5% 95%
## 603.25 1027.00quantile(Cylinders$MinorAxisnm, c(0.05,0.95))## 5% 95%
## 289.7 646.0Centre-to-centre distances
The software provides the centroid distances as the x and y position in the image for each particle, in micrometers. To calculate the centre to centre distance, it is necesary to use the formula for euclidian distances
Centroids <-
Centroids %>%
# Calculate euclidean distance:
mutate(Euclidean = sqrt(((x1-x)^2)+((y2-y)^2))) %>%
# convert from microns to nm
mutate(Euclideannm = Euclidean * 1000)
# Plot:
ggplot(Centroids, aes(x=Euclideannm))+
geom_histogram( color="#ffd781" ,fill="#ffd781", alpha=0.3,
position="identity",binwidth=50)+
theme_bw()+
theme(legend.position="top")+
ylab("Counts")+
xlab("Centre - Centre Distance (nm)")+
theme(legend.position = c(0.8,0.5))+
theme(legend.title = element_blank())+
theme(legend.text = element_text(size=6))
Circularity
This parameter indicates how circular are the particles. When the circularity is closer to 1, this means that the relationship between area and perimeter is close to 1, so the particle has a regular perimeter
# Plot:
ggplot(Cylinders, aes(x=circ.)) +
geom_histogram( color="#0e0928" ,fill="#93c7ff", alpha=0.5,
position="identity")+
theme_bw()+
theme(legend.position="top")+
ylab("Counts")+
xlab("Circularity (4π*area/perimeter^2)")+
scale_x_continuous(breaks=seq(0,1,0.2))+
theme(legend.position = c(0.8,0.5))+
theme(legend.title = element_blank())+
theme(legend.text = element_text(size=6))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Aspect ratio
Major Axis/minor Axis
If this value is far away from 1, it means that the elipse is flattened. if close to 1, it means that the ellipse is close to a circle.
# Plot:
ggplot(Cylinders, aes(x=AR)) +
geom_histogram( color="#0e0928" ,fill="#93c7ff", alpha=0.5,
position="identity")+
theme_bw()+
theme(legend.position="top")+
ylab("Counts")+
xlab("Aspect Ratio (Major Axis/Minor Axis")+
theme(legend.position = c(0.8,0.5))+
theme(legend.title = element_blank())+
theme(legend.text = element_text(size=6))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Angle
This parameter indicates varies between 0-180 degrees, and it is the angle between the primary axis and a line parallel to the x-axis of the image.
# Plot:
ggplot(Cylinders, aes(x=angle)) +
geom_histogram( color="#0e0928" ,fill="#93c7ff", alpha=0.5,
position="identity")+
theme_bw()+
theme(legend.position="top")+
ylab("Counts")+
xlab("Angle (ellipse axis to the image x axis)")+
scale_x_continuous(breaks=seq(0,180,20))+
theme(legend.position = c(0.8,0.5))+
theme(legend.title = element_blank())+
theme(legend.text = element_text(size=6))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Correlations
We evaluated the correlation between the centre to centre distance and the diameter
First we need to create unique codes for the particles including the TEM number, the rectangle number and the ROI number.
In this way, we can combine the two data frames, cylinders diameter and centre to centre distance. Thus, for each particle we will obtain its diameter and its centre-to-centre distance to one of its adjacent particles.
Cylinders2 <-
Cylinders %>%
unite("ROICode",c("section","rectangle","ROI")) %>%
dplyr::select(ROICode,MajorAxisnm,MinorAxisnm)
Centroids2 <-
Centroids %>%
mutate(label2 = substr(label,16,18)) %>%
unite("ROICode",c("label2","ROI"),remove=FALSE) %>%
unite("ROICompCode",c("label2","ROIcomp"),remove=FALSE) %>%
dplyr::select(ROICode,ROICompCode,Euclideannm) %>%
# order:
arrange(ROICompCode) %>%
#remove particles that do not have diameter info:
filter(ROICompCode!="1_4_17")
StructStats <-
Cylinders2 %>%
# select the ROIs for which we have the centroid distances:
filter (ROICode %in% (Centroids2$ROICompCode)) %>%
# order:
arrange(ROICode) %>%
# Join data frames: cylinders with Centroids
bind_cols(.,Centroids2) %>%
dplyr::select(-ROICode...1) %>%
rename(AdjROI = ROICode...4) %>%
dplyr::select(4,1,2,3,5) %>%
filter(MajorAxisnm<1400)
head(StructStats)## ROICompCode MajorAxisnm MinorAxisnm AdjROI Euclideannm
## 1 1_1_11 862 550 1_1_10 962.2333
## 2 1_1_14 765 592 1_1_10 719.4616
## 3 1_1_25 711 596 1_1_30 891.4034
## 4 1_1_27 732 560 1_1_30 610.0467
## 5 1_1_28 728 562 1_1_30 984.9152
## 6 1_1_29 782 521 1_1_30 921.2481In this data frame, the first column is the label of the particles (i.e. transversal section of a cylinder in the white underlay of X. eucalypti) that we are studying. Second column and third column are the diameter. Since the particles are an ellipse, we kept the major and minor axis. Column 4 is the the label of a particle adjacent to the particle in column 1. Column 5 is the centre to centre distance between these two particles.
Now we can test if the diameter of the particle correlates to the centre to centre distance.
cor.test(StructStats$MajorAxisnm,StructStats$Euclideannm)##
## Pearson's product-moment correlation
##
## data: StructStats$MajorAxisnm and StructStats$Euclideannm
## t = 0.57344, df = 84, p-value = 0.5679
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1514332 0.2707390
## sample estimates:
## cor
## 0.06244572cor.test(StructStats$MinorAxisnm,StructStats$Euclideannm)##
## Pearson's product-moment correlation
##
## data: StructStats$MinorAxisnm and StructStats$Euclideannm
## t = -1.1363, df = 84, p-value = 0.2591
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.32640536 0.09121409
## sample estimates:
## cor
## -0.1230391ggplot(StructStats, aes(x = Euclideannm))+
geom_point(aes(x = Euclideannm,y = MajorAxisnm),
size=3, alpha=0.5,pch=21,
col="black",fill="#93c7ff")+
geom_point(aes(x = Euclideannm,y = MinorAxisnm),
size=3, alpha=0.5,pch=21,
col="black",fill="#f8f813")+
theme_bw()+
xlab("Centre-centre distance (nm)")+
ylab("Diameter of particles (nm)")+
scale_x_continuous(breaks=seq(500,1500,200))+
scale_y_continuous(breaks=seq(400,1400,200))+
geom_hline(yintercept=mean(StructStats$MajorAxisnm),
linetype="dashed", color = "#0e0928")+
geom_hline(yintercept=mean(StructStats$MinorAxisnm),
linetype="dashed", color = "#abab0d")
Blue = major axis and yellow = minor axis
We can also test if the sum of the dimeter/2 of two adjacent particles explain its centre to centre distance?
SubCylind <-
Cylinders2 %>%
filter (ROICode %in% (StructStats$AdjROI))
StructStats$MajAxisAdj <- SubCylind$MajorAxisnm [factor(StructStats$AdjROI)]
StructStats$MinAxisAdj <- SubCylind$MinorAxisnm [factor(StructStats$AdjROI)]
StructStats <-
StructStats %>%
mutate(SumMajAxis = MajorAxisnm/2 + MajAxisAdj/2) %>%
mutate(SumMinAxis = MinorAxisnm/2 + MinAxisAdj/2) ModEucAxis <- lm (StructStats$Euclideannm ~ StructStats$SumMajAxis + StructStats$SumMinAxis)
summary(ModEucAxis)##
## Call:
## lm(formula = StructStats$Euclideannm ~ StructStats$SumMajAxis +
## StructStats$SumMinAxis)
##
## Residuals:
## Min 1Q Median 3Q Max
## -390.10 -131.58 -12.37 89.04 661.33
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 729.3293 271.7935 2.683 0.00879 **
## StructStats$SumMajAxis 0.2762 0.2300 1.201 0.23318
## StructStats$SumMinAxis -0.1664 0.4076 -0.408 0.68412
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 203 on 83 degrees of freedom
## Multiple R-squared: 0.01799, Adjusted R-squared: -0.005678
## F-statistic: 0.7601 on 2 and 83 DF, p-value: 0.4709ggplot(StructStats, aes(x = Euclideannm))+
geom_point(aes(x = Euclideannm,y = SumMajAxis),
size=3, alpha=0.5,pch=21,
col="black",fill="#93c7ff")+
geom_point(aes(x = Euclideannm,y = SumMinAxis),
size=3, alpha=0.5,pch=21,
col="black",fill="#f8f813")+
theme_bw()+
xlab("Centre-centre distance (nm)")+
ylab("Sum of r for adjacent particles (nm)")+
scale_x_continuous(breaks=seq(500,1500,200))+
geom_hline(yintercept=mean(StructStats$SumMajAxis),
linetype="dashed", color = "#0e0928")+
geom_hline(yintercept=mean(StructStats$SumMinAxis),
linetype="dashed", color = "#abab0d")
Results
Statistical descriptors of the cylinders in the X. eucalypti white underlay structures (Table 2)
maxA<- round(c(
mean(Cylinders$MajorAxisnm),
sd(Cylinders$MajorAxisnm),
quantile(Cylinders$MajorAxisnm, c(0.25,0.75)),
length(Cylinders$MajorAxisnm)),2)
minA <- round(c(
mean(Cylinders$MinorAxisnm),
sd(Cylinders$MinorAxisnm),
quantile(Cylinders$MinorAxisnm, c(0.25,0.75)),
length(Cylinders$MinorAxisnm)),2)
cent <- round(c(
mean(Centroids$Euclideannm),
sd(Centroids$Euclideannm),
quantile(Centroids$Euclideannm, c(0.25,0.75)),
length(Centroids$Euclideannm)),2)
Ecir <- round(c(
mean(Cylinders$circ.),
sd(Cylinders$circ.),
quantile(Cylinders$circ., c(0.25,0.75)),
length(Cylinders$circ.)),2)
ARat <- round(c(
mean(Cylinders$AR),
sd(Cylinders$AR),
quantile(Cylinders$AR, c(0.25,0.75)),
length(Cylinders$AR)),2)
# For the angles change the reference system
# Either the elypse is aligned with the horizontal
# Or it is aligned with the normal of the elytron
# angle should take values only between 0 and 90
Angs1 <- Cylinders$angle[Cylinders$angle<90] #less than 90
Angs2 <- 180-(Cylinders$angle[Cylinders$angle>90])# >90, invert
Angs3 <- c(Angs1,Angs2) # join them again
Angl <- round(c(
mean(Angs3),
sd(Angs3),
quantile(Angs3, c(0.25,0.75)),
length(Angs3)),2)
Table 2 in the manuscript:
| Parameter | Mean | sd | 25% quantile | 75% quantile | n |
|---|---|---|---|---|---|
| Ellipse Major axis | 807.4 | 168.62 | 741.25 | 883 | 738 |
| Ellipse Minor axis | 510.87 | 109.37 | 462 | 581 | 738 |
| Aspect Ratio | 1.63 | 0.46 | 1.34 | 1.8 | 738 |
| Centre-Centre | 867.69 | 200.47 | 732.28 | 966.54 | 88 |
| Circularity | 0.83 | 0.11 | 0.79 | 0.9 | 738 |
| Angle | 20.85 | 14.86 | 10.04 | 29.06 | 737 |
The particles are ellipses in which the aspect ratio is ~1.6 which means the major axis is ~1.6 times the minor axis.
There is no correlation between the major/minor axis of the ellipses and the centre to centre distance. This is because the particles are not perfectly circular and the spacing between ellipses is determined by the chitin matrix (which varies considerably).
The particles are often at angles. Their major axis can vary between 0 and 30 degrees from the horizontal of the image. This means that if we model the cross section of the cylinders as ellipses, the major axis of those ellipses will often align with the ventral surface of the elytron and almost never with the normal to the elytron.
Procedure
This is a short description of the procedure used in Fiji/ImajeJ to extract the centre to centre distance and the diameter of the transversal section of the cylinders from the TEM images
Here we refer as “particle” to the circle correspondent to a transversal section of a cylinder in the structure. Cylinders are filled with air and surrounded by a matrix of chitin.
One rectangle of one TEM section
Then, we used the plugin Analyze particles in Fiji/Image J to obtain statistical descriptors for the shapes of the transversal sections of the cylinders.
To do this we convert the image to black and white, then apply a threshold:
Example of a threshold to identify particles
Then, the plugin assigns a number to each particle, fits an ellipse to each particle and outputs the basic stats.
Example of the numbers assigned to each particle in the area (each particle is also refered as the region of interes, ROI)
Note: The procedure also includes two extra steps to remove particles that are too small and probably just artifacts (not shown here for brevity).
Models
Photonic crystal
We simplified the structure into parallel cylinders with circular cross-section, and analysed its interaction with light of different wavelengths and polarizations using the Finite Element Method (FEM) as implemented in COMSOL Multiphysics
We compared models in which the cylinders transversal section is a circle of diameter 700, 820 and 1000 nm in two polarizations: incident light polarized parallel (continuous) and perpendicular (dashed) to the cylinder axis
OrderModel <-
Ordered %>%
gather (mod, Reflectance, -wl) %>%
mutate (Polarization = substr(mod,1,7),
Cylinder = substr(mod,8,13))
PolarizA <-
OrderModel %>%
filter (Polarization == "Perpen")
PolarizB <-
OrderModel %>%
filter (Polarization == "Parall")
LINES <- c("ParallP0700nm" = "solid",
"ParallP0820nm" = "solid",
"ParallP1000nm" = "solid",
"PerpenP0700nm" = "dashed",
"PerpenP0820nm" = "dashed",
"PerpenP1000nm" = "dashed")
COLORS <- c("ParallP0700nm" = "#c4a922",
"ParallP0820nm" = "black",
"ParallP1000nm" = "grey",
"PerpenP0700nm" = "#c4a922",
"PerpenP0820nm" = "black",
"PerpenP1000nm" = "grey")
ggplot(OrderModel, aes(x=wl, y = Reflectance,
linetype = factor(mod), color=mod))+
geom_line(size=0.8)+
theme_bw() +
scale_linetype_manual(values = LINES) +
scale_color_manual(values = COLORS )+
theme(legend.position="none")+
xlab("Wavelength (nm)") +
ylab("Reflectance (%)") +
xlim(400,1700) +
theme(axis.title.y = element_text(
margin = margin(t = 0, r = 10, b = 0, l = 0)))
Different Scenarios
In the data set scenarios the names of the columns consist on the description of each scenario as follow:
First 4 characters: Reflectance or transmittance. Since we did not include any loses in the models, these two parameters are inverse.
Following 4 characters: A + 3 numbers. This is the radius of the particle in nm (i.e. of the transversal section of a cylinder), assuming the particles are circles (not ellipses)
Following 4 characters: P + 3 numbers. This is the centre to centre distance between two particles in nm.
Following 4 characters: Indicates the polarization, X polarization is perpendicular to the cylinders and z polarization is parallel to them.
The final letter refers to the light incidence, whether the measurement was done normal (N) or off normal (O). The off normal measurement was done at 30 degrees angle.
Setting up
We used different data frames to explore the effect of the variability on each parameter
Particle Size
# Focus on different particle size keeping the other variables constant
Particle <-
Scenarios %>%
dplyr::select(wl, contains ("Refl")) %>% # select reflectance
dplyr::select(wl, contains ("lN")) %>% # normal incidence
dplyr::select(wl, contains ("P868")) # At fixed distance 868 nm
# Re-arrange: each column is one variable, each row one observation
toplotParticle <-
Particle %>%
gather (name, Reflectance, -wl) %>%
mutate (polarization = substr(name,13,16)) %>%
mutate (sizeQ = substr(name,5,8)) %>%
dplyr::select (1,4,5,3)
# X polarization
tplotParx <-
toplotParticle %>%
dplyr::filter (polarization == "Xpol") %>%
dplyr::select (wl,sizeQ,Reflectance)
#z polarization
tplotParz <-
toplotParticle %>%
dplyr::filter (polarization == "Zpol") %>%
dplyr::select (1,3,4)Spacing
# Focus on different spacing keeping the other variables constant
Spacing <-
Scenarios %>%
dplyr::select(wl, contains ("Refl")) %>% # select reflectance
dplyr::select(wl, contains ("lN")) %>% # normal incidence
dplyr::select(wl, contains ("A330")) # At fixed diameter 330 nm
# Re-arrange: each column is one variable, each row one observation
tplotSpac <-
Spacing %>%
gather (name, Reflectance, -wl) %>%
mutate (polarization = substr(name,13,16)) %>%
mutate (spacingQ = substr(name,9,12)) %>%
dplyr::select (1,4,5,3)
# Consider only the z-polarization for visualization purposes
tplotSpacz <-
tplotSpac %>%
filter (polarization == "Zpol") %>%
dplyr::select (1,3,4)Polarization
# Focus on different polarizations keeping the other variables constant
tplotpol <-
tplotSpac %>% #normal incidence and fixed diameter 330 nm
filter (spacingQ == "P868") %>% # spacing at 868 nm
dplyr::select (1,2,4)Angle
# Focus on different angle of incidence keeping the other variables constant
AngleR <-
Scenarios %>%
dplyr::select(wl, contains ("Refl")) %>% # select reflectance
dplyr::select(wl, contains ("A330")) %>% # radius 330 nm
dplyr::select(wl, contains ("P868")) %>% # At fixed distance 868 nm
dplyr::select(wl, contains("Zpol")) # only z polarization
# Re-arrange: each column is one variable, each row one observation
toplotAngle <-
AngleR %>%
gather (name, Reflectance, -wl) %>%
mutate (Angle = substr(name,17,17)) %>%
dplyr::select (1,3,4) Average
toplotAvg <-
Scenarios %>%
dplyr::select(wl, contains ("Refl")) %>% # select reflectance
dplyr::select(wl, contains ("lN")) %>% # normal incidence
mutate(AvgReflxpol =
rowMeans(.[,c(2,4,6,8,10)])) %>%
mutate(AvgReflzpol =
rowMeans(.[,c(3,5,7,9,11)])) %>%
dplyr::select(wl,AvgReflxpol, AvgReflzpol)4 pannels
Here we display the effect of the different scenarios: Particle size, centre-to-centre distance, polarization and angle (Figure 5 in the manuscript).
We first create each of the 4 figures separately:
Particle
S1 <- ggplot(data= tplotParz[tplotParz$sizeQ!="A330",], aes(x=wl,y=Reflectance)) +
geom_line(aes(col=sizeQ, group=sizeQ, linetype= sizeQ))+
scale_color_manual(values=c("#93c7ff","black"),
labels=c("301nm ","366nm ")) +
scale_linetype_manual(values=c(1,1),
labels=c("301nm ","366nm ")) +
theme_test() +
theme(legend.position = c(0.86, 0.7),
legend.title= element_blank(),
legend.key.size = unit(0.2, 'cm'),
legend.spacing.y = unit(0, "mm"),
legend.background = element_blank(),
legend.text=element_text(size=6.5),
legend.key.height = unit(0.3, 'cm'), #change legend key height
legend.key.width = unit(0.4, 'cm')) +#change legend key width
scale_y_continuous(breaks=c(0,0.5,1)) +
theme(axis.title.x=element_blank())+
scale_x_continuous(breaks=c(400,800, 1200,1600))+
theme(axis.title.y = element_text(margin=margin(r=10)))+
ylab("Reflectance (fraction)") +
theme(plot.margin = unit(c(0.5,0.1,0.05,0.3), "cm"))Spacing
S2 <- ggplot(data= tplotSpacz[tplotSpacz$spacingQ!="P868",], aes(x=wl,y=Reflectance)) +
geom_line(aes(col=spacingQ, group=spacingQ, linetype=spacingQ))+
scale_color_manual(values=c("#93c7ff","black"),
labels=c("732nm ","967nm ")) +
scale_linetype_manual(values=c(1,1),
labels=c("732nm ","967nm ")) +
theme_test() +
theme(legend.position = c(0.86, 0.7),
legend.title= element_blank(),
legend.key.size = unit(0.2, 'cm'),
legend.spacing.y = unit(0, "mm"),
legend.background = element_blank(),
legend.text=element_text(size=6.5),
legend.key.height = unit(0.3, 'cm'), #change legend key height
legend.key.width = unit(0.4, 'cm')) + #change legend key width
scale_y_continuous(breaks=c(0,0.5,1)) +
theme(axis.title.x=element_blank(),
axis.title.y=element_blank())+
scale_x_continuous(breaks=c(400,800, 1200,1600)) +
theme(plot.margin = unit(c(0.5,0.1,0.05,0.1), "cm"))Polarizations
S3 <- ggplot(data= tplotpol, aes(x=wl,y=Reflectance)) +
geom_line(aes(col=polarization,
group=polarization,
linetype=polarization)) +
scale_color_manual(values=c("#93c7ff","black"),
labels=c("x pol ","z pol ")) +
scale_linetype_manual(values=c(1,1),
labels=c("x pol ","z pol ")) +
theme_test() +
theme(legend.position = c(0.86, 0.7),
legend.title= element_blank(),
legend.key.size = unit(0.2, 'cm'),
legend.spacing.y = unit(0, "mm"),
legend.background = element_blank(),
legend.text=element_text(size=6.5),
legend.key.height = unit(0.3, 'cm'), #change legend key height
legend.key.width = unit(0.4, 'cm'))+#change legend key width
scale_y_continuous(breaks=c(0,0.5,1)) +
theme(axis.title.y = element_text(margin=margin(r=10)))+
scale_x_continuous(breaks=c(400,800, 1200,1600))+
xlab("wavelength (nm)")+
ylab("Reflectance (fraction)")+
theme(plot.margin = unit(c(0.1,0.1,0.05,0.3), "cm"))Angle
S4 <- ggplot(data= toplotAngle, aes(x=wl,y=Reflectance))+
geom_line(aes(col=Angle,
group=Angle,
linetype=Angle)) +
scale_color_manual(values=c("black","#93c7ff"),
labels=c("Norm ","30 deg ")) +
scale_linetype_manual(values=c(1,1),
labels=c("Norm ","30 deg ")) +
theme_test() +
theme(legend.position = c(0.86, 0.7),
legend.title= element_blank(),
legend.key.size = unit(0.2, 'cm'),
legend.spacing.y = unit(0, "mm"),
legend.background = element_blank(),
legend.text=element_text(size=6.5),
legend.key.height = unit(0.3, 'cm'), #change legend key height
legend.key.width = unit(0.4, 'cm')) + #change legend key width
scale_y_continuous(breaks=c(0,0.5,1)) +
theme(axis.title.y=element_blank())+
scale_x_continuous(breaks=c(400,800, 1200,1600)) +
xlab("wavelength (nm)") +
theme(plot.margin = unit(c(0.1,0.1,0.05,0.1), "cm"))Then join them together to compare:
plot_grid(S1, S2, S3, S4,
align = "hv", nrow=2, ncol=2, rel_widths = c(1,1),
scale = 1, vjust=0.5)
Combined scenarios
Multiple scenarios can occur at the same time.
We wanted to recreate the general trend for the additive effects of these scenarios. Thus, we obtained weighted averages for different polarization and different angle. (Figure 8 in the supplementary materials)
This means we averaged all the models with normal incidence and different particle size and spacing by polarization. And we obtained a data frame with the two polarizations in the two different angles of incidence
# Average per polarization and angle
WSc <-
Scenarios %>%
dplyr::select(wl, contains ("Refl")) %>% # select reflectance
mutate(xpolNorm = rowMeans(
dplyr::select(.,ends_with ("XpolN")), na.rm = TRUE)) %>% # average x pol Norm
mutate(zpolNorm = rowMeans(
dplyr::select(.,ends_with ("ZpolN")), na.rm = TRUE)) %>% # average z pol Norm
rename(zpol30 = ReflA330P868ZpolO) %>%
rename(xpol30 = ReflA330P868XpolO) %>%
dplyr::select(1,12,13,14,15)
head(WSc)## wl xpol30 zpol30 xpolNorm zpolNorm
## 1 350 0.89183 0.78633 0.2474110 0.383458
## 2 352 0.87658 0.78852 0.3282038 0.371852
## 3 354 0.68220 0.80385 0.4052974 0.356248
## 4 356 0.45578 0.85418 0.4529960 0.361314
## 5 358 0.70678 0.84597 0.4458740 0.330380
## 6 360 0.65343 0.79169 0.3805180 0.304974# Arrange for plotting
toplotWSc <-
WSc %>%
gather (Lightcondition, Reflectance, -wl)
# Find the weighted mean
Avvline <-
WSc %>%
mutate(All = rowMeans(WSc[2:5])) %>%
mutate(xpol = rowMeans(WSc[c(2,4)])) %>%
mutate(zpol = rowMeans(WSc[c(3,5)])) %>%
dplyr::select(1,6,7,8)
head(Avvline)## wl All xpol zpol
## 1 350 0.5772572 0.5696205 0.584894
## 2 352 0.5912890 0.6023919 0.580186
## 3 354 0.5618988 0.5437487 0.580049
## 4 356 0.5310675 0.4543880 0.607747
## 5 358 0.5822510 0.5763270 0.588175
## 6 360 0.5326530 0.5169740 0.548332ggplot()+
geom_line(data= toplotWSc, aes(x=wl,y=Reflectance,
group= Lightcondition, color=Lightcondition,
linetype=Lightcondition)) +
scale_linetype_manual(values=c("longdash","solid","longdash", "solid"))+
scale_color_manual(values=c("#93c7ff","#93c7ff","grey","grey"))+
# geom_line(data=Avvline, aes(x=wl, y=zpol), lwd=1.3, col="black")+
# geom_line(data=Avvline, aes(x=wl, y=xpol), lwd=1.3, col="blue")+
geom_line(data=Avvline, aes(x=wl, y=All), lwd=1.3, col="black")+
theme_bw()+
ylab ("Standardized Reflectance") +
xlab("Wavelength (nm)") 
Quasi-ordered model
The plot shows the probability of reflectance of light inciding perpendicular to the cylinder axis. (Supplementary materials 9)
PlotM1 <-
QuasiOrder %>%
mutate(., MeanS1 = rowMeans(dplyr::select(., starts_with("Xyln1."))),
., MeanS2 = rowMeans(dplyr::select(., starts_with("Xyln2."))),
., MeanS3 = rowMeans(dplyr::select(., starts_with("Xyln3."))),
., MeanS4 = rowMeans(dplyr::select(., starts_with("Xyln4."))),
., MeanS5 = rowMeans(dplyr::select(., starts_with("Xyln5."))),
., MeanS6 = rowMeans(dplyr::select(., starts_with("Xyln6."))),
., MeanS7 = rowMeans(dplyr::select(., starts_with("Xyln7."))),
., MeanS8 = rowMeans(dplyr::select(., starts_with("Xyln8."))),
., MeanS9 = rowMeans(dplyr::select(., starts_with("Xyln9.")))) %>%
dplyr::select(wl,contains("Mean"))
PlotM1Means <-
PlotM1 %>%
mutate(Xyl.sd = rowSds(as.matrix(PlotM1[,2:10]))) %>%
mutate(Xyl.mean = rowMeans2(as.matrix(PlotM1[,2:10]))) %>%
dplyr::select(wl,Xyl.sd,Xyl.mean)
ggplot(PlotM1Means, aes(x=wl)) +
geom_ribbon(aes(y = Xyl.mean*100,
ymin = Xyl.mean*100 - Xyl.sd*100,
ymax = Xyl.mean*100 + Xyl.sd*100,
fill = "#DCDCDC"), alpha = .4) +
geom_line(aes(x=wl, y=Xyl.mean*100),
size=1, colour="black")+
scale_fill_manual(values=c("#DCDCDC"))+
theme_bw() +
ylab("Probablity of Reflectance x 100")+
xlab("wavelength (nm)")+
theme(legend.position="none")+
xlim(400,1700) +
theme(axis.title.y = element_text(
margin = margin(t = 0, r = 10, b = 0, l = 0)))## Warning: Removed 14 row(s) containing missing values (geom_path).