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

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.00
quantile(Cylinders$MinorAxisnm, c(0.05,0.95))
##    5%   95% 
## 289.7 646.0

Centre-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.2481

In 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.06244572
cor.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.1230391
ggplot(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.4709
ggplot(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

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

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)

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.548332
ggplot()+
  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).