added tests

R/scattnlay.r

@@ -1,7 +1,7 @@
 Layer <- setClass("Layer",
                   representation(
                     m="complex",
-                    d="numeric"
+                    r="numeric"
                     ))
 
 Scatterer <- setClass("Scatterer",
@@ -13,7 +13,7 @@ Scatterer <- setClass("Scatterer",
 
 setMethod("show", signature(object="Layer"),function(object){
   cat("Layer Data \n")
-  cat("Diameter:", object@d," nm\n")
+  cat("Diameter:", 2*object@r," nm\n")
   cat("Epsilon:", object@m,"\n")
 })
 
@@ -24,12 +24,12 @@ setMethod("show", signature(object="Scatterer"),function(object){
 
 setGeneric("na<-",function(x,value) standardGeneric("na<-"))
 setGeneric("lambda<-",function(x,value) standardGeneric("lambda<-"))
-setGeneric("d<-",function(x,value) standardGeneric("d<-"))
+setGeneric("r<-",function(x,value) standardGeneric("r<-"))
 setGeneric("m<-",function(x,value) standardGeneric("m<-"))
 
 setReplaceMethod("na","Scatterer", function(x,value) {x@na <- value; validObject(x); x})
 setReplaceMethod("lambda","Scatterer", function(x,value) {x@lambda <- value; validObject(x); x})
-setReplaceMethod("d","Layer", function(x,value) {x@d <- value; validObject(x); x})
+setReplaceMethod("r","Layer", function(x,value) {x@r <- value; validObject(x); x})
 setReplaceMethod("m","Layer", function(x,value) {x@m <- value; validObject(x); x})
 
 setGeneric("scattnlay", function(object) {standardGeneric("scattnlay")})

oil_scat.R

@@ -0,0 +1,27 @@
+library(Rscattnlay)
+library(ggplot2)
+library(reshape2)
+
+S <- Scatterer()  # create a scatterer object
+
+dr <- Layer()     # and a layer representing the droplet
+na(S) <- 1.33    # set the ambient index
+lambda(S) <- 600
+
+m(dr) <- 1.431+0i;
+
+scatnlay_ar <- function(r){
+  r(dr) <- r
+  St <- S+dr;              # PUT THE STACK HERE
+  Q <- scattnlay(St);
+}
+
+dr_size <- c(1:100)
+
+output <- sapply(dr_size,scatnlay_ar)
+
+data_mat <- data.frame(t(output),dr_size,'oil')
+colnames(data_mat) <- c("Qext", "Qsca", "Qabs", "Qbk", "Qpr", "g", "Albedo", "nmax","Size","Layers")
+
+scat_plot <- ggplot(data=data_mat, aes(x=Size,y=Qsca)) 
+scat_plot + geom_line() + theme_minimal(base_size=22) + xlim(c(0,100)) + xlab("Radius (nm)")

test_2.R

@@ -0,0 +1,81 @@
+library(Rscattnlay)
+library(ggplot2)
+library(reshape2)
+
+S <- Scatterer()  # create a scatterer object
+
+np <- Layer()     # and a layer representing the np
+lipid <- Layer()  # and other layers
+inner_layer <- Layer()
+water <- Layer()
+
+na(S) <- 1.33     # set the ambient index
+r(np) <- 35       # and the size (radius) of the np in nm
+
+r(water) <- 40
+m(water) <- 1.33+0i
+
+r(inner_layer) <- 38 
+m(inner_layer) <- 1.33+0i
+
+r(lipid) <- 40    # and the size (radius) of the other layers
+m(lipid) <- 1.45+0i   #fixed value
+
+# load up silver data
+palik_ag_vis <- read.table("palik_ag_vis_hb.csv", header=TRUE)
+colnames(palik_ag_vis) <- c("lambda","n","k","eps_real","eps_imag")
+
+lambda_palik <- palik_ag_vis$lambda #convert to nm
+n_palik <- palik_ag_vis$n 
+k_palik <- palik_ag_vis$k
+
+#interpolation as palik data is sparse
+n = 512 #number of points
+spl_n <- approx(lambda_palik,n_palik,n=n)
+spl_k <- approx(lambda_palik,k_palik,n=n)
+lambda = spl_n$x
+
+scatnlay_ar <-function(k){
+  m(np) <- spl_n$y[k]+spl_k$y[k]*(0+1i);
+  lambda(S) <- lambda[k];
+  St <- S+np+water;              # PUT THE STACK HERE
+  Q <- scattnlay(St);
+}
+
+scatnlay_ar_water_lipid <-function(k){
+  m(np) <- spl_n$y[k]+spl_k$y[k]*(0+1i);
+  lambda(S) <- lambda[k];
+  St <- S+np+inner_layer+lipid;              # PUT THE STACK HERE
+  Q <- scattnlay(St);
+}
+
+output <- sapply(1:n,scatnlay_ar)
+output_lipid <- sapply(1:n,scatnlay_ar_water_lipid)
+data_mat_np <- data.frame(t(output),lambda,'water')
+data_mat_lipid <- data.frame(t(output_lipid),lambda,'lipid')
+
+colnames(data_mat_np) <- c("Qext", "Qsca", "Qabs", "Qbk", "Qpr", "g", "Albedo", "nmax","Lambda","Layers")
+colnames(data_mat_lipid) <- c("Qext", "Qsca", "Qabs", "Qbk", "Qpr", "g", "Albedo", "nmax","Lambda","Layers")
+
+data_mat <- rbind(data_mat_np,data_mat_lipid)
+
+ext_water <- subset(data_mat_np, Lambda < 700, select= c("Qext","Lambda"))
+ext_water$Lambda[which.max(ext_water$Qext)]
+
+ext_lipid <- subset(data_mat_lipid, Lambda < 700, select= c("Qext","Lambda"))
+ext_lipid$Lambda[which.max(ext_lipid$Qext)]
+
+dif_Qext <- data.frame(ext_lipid$Qext-ext_water$Qext) ## might be able to use aggregate here
+dif_spectra <- cbind(ext_lipid$Lambda,dif_Qext)
+colnames(dif_spectra) <- c("lambda","D_Qext")
+
+palik_ag_melt <- melt(subset(palik_ag_vis,select = c("lambda","n","k")), id=c("lambda"),variable.name = "type", 
+                      value.name = "Index")
+ext_plot <- ggplot(data=data_mat, aes(x=Lambda,y=Qext, group=as.factor(Layers), color=as.factor(Layers))) 
+ext_plot + geom_line() + theme_minimal(base_size=22) + xlim(c(350,550))+ ylim(c(0,10)) + labs(colour = "Experiment")  + xlab("Wavelength (nm)")
+
+scat_plot <- ggplot(data=data_mat, aes(x=Lambda,y=Qsca, group=as.factor(Layers), color=as.factor(Layers))) 
+scat_plot + geom_line() + theme_minimal(base_size=22) + xlim(c(350,550))+ ylim(c(0,7.5)) + labs(colour = "Experiment")  + xlab("Wavelength (nm)")
+
+dif_plot <- ggplot(data=dif_spectra, aes(x=lambda,y=D_Qext))
+dif_plot + geom_line() + theme_minimal(base_size=22) + xlim(c(350,550)) + xlab("Wavelength (nm)")

test_3.R

@@ -0,0 +1,53 @@
+#scattnlay()
+library(Rscattnlay)
+library(ggplot2)
+library(reshape2)
+
+S <- Scatterer()  # create a scatterer object
+
+np <- Layer()     # and a layer representing the np
+na(S) <- 1.33    # set the ambient index
+r(np) <- 55      # and the size (radius) of the np in nm
+
+# load up silver data
+palik_ag_vis <- read.table("palik_ag_vis_hb.csv", header=TRUE)
+colnames(palik_ag_vis) <- c("lambda","n","k","eps_real","eps_imag")
+
+lambda_palik <- palik_ag_vis$lambda #convert to nm
+n_palik <- palik_ag_vis$n 
+k_palik <- palik_ag_vis$k
+
+#interpolation as palik data is sparse
+n = 500#number of points
+spl_n <- approx(lambda_palik,n_palik,n=n)
+spl_k <- approx(lambda_palik,k_palik,n=n)
+lambda = spl_n$x
+
+scatnlay_ar <-function(k){
+  m(np) <- spl_n$y[k]+spl_k$y[k]*(0+1i);
+  lambda(S) <- lambda[k];
+  St <- S+np;              # PUT THE STACK HERE
+  Q <- scattnlay(St);
+}
+
+output <- sapply(1:n,scatnlay_ar)
+data_mat_np <- data.frame(t(output),lambda,'water')
+
+Data1 <- data.frame(read.table("/media/mbajb/data/Linux/Comsol Projects/Nottingham Simulations/Data8.txt",skip=5))
+colnames(Data1) <- c("Lambda","Frequency","Qscat")
+Data1$scal <- Data1$Qscat/max(Data1$Qscat)
+Data1$Lambda <- Data1$Lambda*1e9
+
+colnames(data_mat_np) <- c("Qext", "Qsca", "Qabs", "Qbk", "Qpr", "g", "Albedo", "nmax","Lambda","Layers")
+
+data_mat <- subset(data_mat_np,Lambda<700 & Lambda >350)
+#data_mat$Qsca <- data_mat$Qsca/max(data_mat$Qsca)
+
+palik_ag_melt <- melt(subset(palik_ag_vis,select = c("lambda","n","k")), id=c("lambda"),variable.name = "type", 
+                      value.name = "Index")
+
+scat_plot <- ggplot(data=data_mat, aes(x=Lambda,y=Qext)) 
+scat_plot + geom_line() + theme_minimal(base_size=22) + xlim(c(300,500)) + xlab("Wavelength (nm)")
+
+#scat_plot2 <- ggplot(data=Data1, aes(x=Lambda,y=scal)) 
+#scat_plot2 + geom_line() + theme_minimal(base_size=22) + xlim(c(350,700))+ ylim(c(0,1.2)) + xlab("Wavelength (nm)")