[Photonics-commits] r328 - in pkg: . baptMisc/R mie mie/R mie/data mie/inst mie/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Nov 15 13:36:10 CET 2008
Author: baptiste
Date: 2008-11-15 13:36:10 +0100 (Sat, 15 Nov 2008)
New Revision: 328
Added:
pkg/baptMisc/R/google.r
pkg/mie/
pkg/mie/DESCRIPTION
pkg/mie/R/
pkg/mie/R/.Rhistory
pkg/mie/R/Eint.r
pkg/mie/R/Escat.r
pkg/mie/R/Mie.R
pkg/mie/R/Mie.S.r
pkg/mie/R/Mie.ab.R
pkg/mie/R/Mie.cd.r
pkg/mie/R/Mie.coated.R
pkg/mie/R/Mie.coated.ab.R
pkg/mie/R/Mie.coated.spectra.R
pkg/mie/R/Mie.pt.r
pkg/mie/R/Mie.shell.r
pkg/mie/R/Mie.spectra.R
pkg/mie/R/hankel.r
pkg/mie/R/permittivity.r
pkg/mie/R/plot.material.R
pkg/mie/data/
pkg/mie/data/AuJChristy.rda
pkg/mie/inst/
pkg/mie/inst/.Rhistory
pkg/mie/inst/MieNearFieldPlot.png
pkg/mie/inst/MieNearFieldPlot2.png
pkg/mie/inst/Rplot.pdf
pkg/mie/inst/fieldPlot.pdf
pkg/mie/inst/mieNearField.rda
pkg/mie/inst/test.r
pkg/mie/inst/test2.r
pkg/mie/man/
pkg/mie/man/Mie.Rd
pkg/mie/testMat.r
Log:
mie first import
Added: pkg/baptMisc/R/google.r
===================================================================
--- pkg/baptMisc/R/google.r (rev 0)
+++ pkg/baptMisc/R/google.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,11 @@
+google <- function (string, rsite=FALSE, rarchive=FALSE){
+ ronly <- NULL
+ if (rsite & rarchive) stop("can't have both!")
+ if(rsite) ronly <- "sitesearch=r-project.org&"
+ RURL = "http://www.google.com/search"
+ if(rarchive) RURL <- "http://www.googlesyndicatedsearch.com/u/newcastlemaths"
+ RSearchURL = paste(RURL, "?", ronly, "q=",
+ string, sep = "")
+ browseURL(RSearchURL)
+ return(invisible(0))
+}
Added: pkg/mie/DESCRIPTION
===================================================================
--- pkg/mie/DESCRIPTION (rev 0)
+++ pkg/mie/DESCRIPTION 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,13 @@
+Package: mie
+Type: Package
+Title: Scattering related functions
+Version: 1.0
+Date: 2008-03-09
+Author: Baptiste Auguie
+Maintainer: baptiste <ba208 at ex.ac.uk>
+Description: Mie scattering for a sphere and a coated sphere.
+Depends: R (>= 1.8.0)
+License: GPL (>=2)
+URL: http://r-forge.r-project.org/projects/photonics/
+Packaged: Mon Mar 24 14:16:06 2008; baptiste
+LazyLoad: yes
\ No newline at end of file
Added: pkg/mie/R/.Rhistory
===================================================================
--- pkg/mie/R/.Rhistory (rev 0)
+++ pkg/mie/R/.Rhistory 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,475 @@
+library(EricksenLeslie)
+demo(pac=EricksenLeslie)
+demo(pac="EricksenLeslie")
+demo(pac="EricksenLeslie",flow_hom)
+#
+data(permittivitiesLuxpop) #
+plotMaterial("Cr",all=FALSE,range=c(0.5,0.8))
+#
+data(permittivitiesLuxpop) #
+plot.material("Cr",all=FALSE,range=c(0.5,0.8))
+?constants
+#
+data(luxpop) #
+plot.material("Cr",all=FALSE,range=c(0.5,0.8))
+baseUrl="http://www.luxpop.com/Material/"#
+material="Cr.nk"#
+readLines(url(paste(baseUrl,material,sep=""), open = "", blocking = TRUE, encoding = getOption("encoding")))->file#
+grep(";",file,value=TRUE)->comment#
+comment#
+read.table(textConnection(file),skip=length(comment))->values#
+head(values)#
+names(values)<-c("l","n","k")#
+attach(values)#
+nk2epsilon(list(n=n,k=k))->valuesEpsilon#
+valuesEpsilon$lambda<-l/1e4 #
+detach(values)#
+attach(valuesEpsilon)#
+#
+graphics.off()#
+par(mfrow=c(2,2),mar=c(5, 5, 1, 1),bty="n",pty="m")#
+cutePlot(lambda,epsr)#
+cutePlot(lambda,epsi,yl=expression(paste(epsilon[i])))#
+cutePlot(lambda,values$n,yl=expression(paste("n")))#
+cutePlot(lambda,values$k,yl=expression(paste("k")))
+version()
+help()
+sessionInfo
+sessionInfo()
+library(grid)
+#
+#
+#
+#
+#
+squareGrob <- grob(height=unit(.1, "npc"),#
+ cl="squareGrob")#
+#
+#
+#
+#
+widthDetails.squareGrob <- function(x) {#
+ convertUnit(unit(.1, "npc"),#
+ "npc",#
+ "y", "dimension",#
+ "x", "dimension")#
+}#
+#
+#
+#
+lay1 <- grid.layout(10, 2,#
+ widths=unit.c(unit(1, "null"),#
+ grobWidth(squareGrob)),#
+ heights=unit(.1, "npc"))#
+#
+... Once that squareGrob is defined, you can use it anywhere, so that's#
+a one-off cost. The following code shows the resulting layout ...#
+#
+#
+pushViewport(viewport(layout=lay1))#
+#
+pushViewport(viewport(layout.pos.col=1))#
+grid.rect(gp=gpar(fill="grey"))#
+popViewport()#
+#
+for (i in 1:10) {#
+ pushViewport(viewport(layout.pos.col=2,#
+ layout.pos.row=i))#
+ grid.rect(gp=gpar(fill="light grey"))#
+ popViewport()#
+}
+?spectra
+library(constants)
+spectra()->test
+test
+library(spectra)
+spectra2excel(test,"test.txt")
+getWavelength(spectra())
+source("http://r.research.att.com/survey.R")
+system("mate /var/folders/yF/yFY+OcggEhuLb4YeNRUrPU+++TI/-Tmp-//RtmptFKYkp/downloaded_packages")
+system("open /var/folders/yF/yFY+OcggEhuLb4YeNRUrPU+++TI/-Tmp-//RtmptFKYkp/downloaded_packages")
+library(baptMisc)
+library(SPPr)
+?SPPr
+graphics.off()#
+library(SPPr)#
+#
+data(scanSPP)#
+data(parameters)#
+head(scanSPP)#
+#
+parameters$zeroShift <- -0.4#
+parameters$thickness[3] <- 0#
+#
+dataNorm <- data2norm(scanSPP,parameters)#
+head(dataNorm)#
+#
+plotData(dataNorm, int=TRUE, col=1, ylim=c(0, 1))#
+criticalEdge(parameters)#
+#
+layer2fit <- 1#
+pIni <- c(-18,0.7, 45e-9, 1.05)#
+#
+plotModel(pIni, int=TRUE, col=2)#
+#
+result <- optim(pIni, objF)#
+#
+yFin <- fitModel(result$par)#
+yErr<-abs(yFin-dataNorm$N)#
+#
+plotModel(result$par,col=3)
+result
+graphics.off()#
+library(SPPr)#
+#
+data(scanSPP)#
+data(parameters)#
+head(scanSPP)#
+#
+parameters$zeroShift <- -0.4#
+parameters$thickness[3] <- 0#
+#
+dataNorm <- data2norm(scanSPP,parameters)#
+head(dataNorm)#
+#
+plotData(dataNorm, int=TRUE, col=1, ylim=c(0, 1))#
+criticalEdge(parameters)#
+#
+layer2fit <- 1#
+pIni <- c(-18,0.7, 45, 1.05)#
+#
+plotModel(pIni, int=TRUE, col=2)#
+#
+result <- optim(pIni, objF)#
+#
+yFin <- fitModel(result$par)#
+yErr<-abs(yFin-dataNorm$N)#
+#
+plotModel(result$par,col=3)
+#
+data(scanSPP)#
+head(scanSPP)
+x.axis <- seq(from=0.5, to=4.5, length.out=13112)
+which.min(abs(x.axis - 3.2))
+?lattice
+?xyplot
+brewer.pal
+xyplot
+library(fortune)
+library(fortunes)
+fortune()
+library(spatstat)
+? spatstat
+citation()
+citation("spatstat")
+#
+#
+xrange=c(-15,15)#
+yrange=c(0,16)#
+plot(0,xlim=xrange,ylim=yrange,type='n')#
+#
+#
+yr=seq(yrange[1],yrange[2],len=50)#
+offsetFn=function(y){2*sin(0+y/3)}#
+offset=offsetFn(yr)#
+leftE = function(y){-10-offsetFn(y)}#
+rightE = function(y){10+offsetFn(y)}#
+#
+#
+xp=c(leftE(yr),rev(rightE(yr)))#
+yp=c(yr,rev(yr))#
+polygon(xp,yp,col="#
+#
+#
+#
+#
+h=9#
+#
+#
+xt=seq(0,rightE(h),len=100)#
+yt=log(1+log(1+log(xt+1)))#
+yt=yt-min(yt)#
+yt=h*yt/max(yt)#
+#
+x=c(leftE(h),rightE(h),rev(xt),-xt)#
+y=c(h,h,rev(yt),yt)#
+polygon(x,y,col="red",border=NA)
+#
+library(Hmisc)#
+sexyNo<-function(){#
+ for(i in (1:20)){#
+ S<-paste(makeNstr(" ",i),"0",makeNstr(" ",(40-2*i)),"<--1\n",sep="")#
+ cat(S); Sys.sleep(1)#
+ }#
+ cat(paste(makeNstr(" ",21),"0~1\n",sep="")); Sys.sleep(2)#
+ for(i in (1:9)){#
+ cat(paste(makeNstr(" ",22),"|\n",sep="")); Sys.sleep(2)#
+ }#
+ cat(" ***!!! 0.1 !!!***\n")#
+}#
+#
+sexyNo()
+plotCI
+?plotCI
+load("/Users/baptiste/Documents/PhD/DraftsPapers/disorderArrays/analysis/diagModel.rda")
+library(ggplot2)
+River.Mile <- c(215 ,202, 198, 190, 185, 179, 148, 119, 61)#
+Cu <- rnorm(9)#
+Fe <- rnorm(9)#
+Mg <- rnorm(9)#
+Ti <- rnorm(9)#
+Ir <- rnorm(9)#
+r <- data.frame(River.Mile, Cu, Fe, Mg, Ti, Ir)#
+#
+z <- melt.data.frame(r, id.var="River.Mile")#
+#
+#
+qplot(River.Mile, value, facets=(variable~.), data=z)
+#
+xyplot(value~River.Mile | variable, data=z)
+?plotmatrix
+data(mtcars)
+#
+plotmatrix(mtcars[, 1:3])
+plotmatrix
+#
+plotmatrix( data=z)
+library(zoo)#
+d<-(structure(c(1.39981554315924, 0.89196314359498, 0.407816250252697,#
+0.823496839063978, 1.14429021220358, 1.23971035967413, 0.960868900583432,#
+0.927685306209829, 1.22072345292821, 0.249842897450642, 1.00879641624694,#
+0.925372139878243, 0.317259909172362, 0.382677149697482), index =#
+structure(c(11808,#
+11869, 11961, 11992, 12084, 12173, 12265, 12418, 12600, 12631,#
+12753, 12996, 13057, 13149), class = "Date"), class = "zoo"))#
+#
+plot(d)
+d
+library(ifs)
+?ifs
+#
+#
+#
+#
+#
+IFS.est <- IFS(y)#
+xx <- IFS.est$x#
+tt <- IFS.est$y#
+#
+ss <- pbeta(xx,.5,.1)#
+#
+#
+#
+ #
+par(mfrow=c(2,1)) #
+ #
+plot(ecdf(y),xlim=c(0,1),main="IFS estimator versus EDF")#
+lines(xx,ss,col="blue")#
+lines(xx,tt,col="red")
+#
+y<-rbeta(50,.5,.1)#
+#
+#
+#
+#
+IFS.est <- IFS(y)#
+xx <- IFS.est$x#
+tt <- IFS.est$y#
+#
+ss <- pbeta(xx,.5,.1)#
+#
+#
+#
+ #
+par(mfrow=c(2,1)) #
+ #
+plot(ecdf(y),xlim=c(0,1),main="IFS estimator versus EDF")#
+lines(xx,ss,col="blue")#
+lines(xx,tt,col="red")
+a=-0.7;b=-0.4 #
+Limits=c(-2,2)#
+MaxIter=60#
+cl=colours()#
+Step=seq(Limits[1],Limits[2],by=0.01)#
+PointsMatrix=array(0,dim=c(length(Step)*length(Step),3))#
+a1=0#
+plot(0,0,xlim=Limits,ylim=Limits,col="white")#
+#
+for(x in Step)#
+{#
+ for(y in Step)#
+ {#
+ n=0#
+ DIST=0#
+ x1=x;y1=y #
+ while(n<MaxIter & DIST<4)#
+ {#
+ newx=x1^2-y1^2+a#
+ newy=2*x1*y1+b#
+ DIST=newx^2+newy^2#
+ x1=newx;y1=newy#
+ n=n+1#
+ }#
+ if(DIST<4) colour=24 else colour=n*10#
+ #
+ a1=a1+1#
+ PointsMatrix[a1,]=c(x,y,colour)#
+ }#
+}#
+#
+plot(PointsMatrix[,1], PointsMatrix[,2], xlim=Limits, ylim=Limits, col=cl[PointsMatrix[,3]], pch=".")
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+p <- xyplot(y~x, data=df)#
+pushViewport(viewport(x=0.3, width = 0.5, height = 0.8, angle = 0,#
+name = "topvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+p <- xyplot(y~x, data=df)#
+pushViewport(viewport(x=0.3, width = 0.5, height = 0.8, angle = 0,#
+name = "topvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+grid.text("some text elsewhere", x=0.8, gp = gpar(col = rgb(43/255,#
+140/255, 190/255)))
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+df2 <- data.frame(x=rnorm(100), y=runif(x))#
+#
+p <- xyplot(y~x, data=df)#
+p2 <- xyplot(y~x, data=df2)#
+pushViewport(viewport(x=0.25, width = 0.5, height = 0.8, angle = 0,#
+name = "leftvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+pushViewport(viewport(x=0.75, width = 0.5, height = 0.8, angle = 0,#
+name = "rightvp"))#
+print(p2, newpage=F)
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+df2 <- data.frame(x <- rnorm(100), y=runif(x))#
+#
+p <- xyplot(y~x, data=df)#
+p2 <- xyplot(y~x, data=df2)#
+pushViewport(viewport(x=0.25, width = 0.5, height = 0.8, angle = 0,#
+name = "leftvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+pushViewport(viewport(x=0.75, width = 0.5, height = 0.8, angle = 0,#
+name = "rightvp"))#
+print(p2, newpage=F)
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+df2 <- data.frame(x <- rnorm(100), y=runif(x))#
+#
+p <- xyplot(y~x, data=df)#
+p2 <- xyplot(y~x, data=df2)#
+pushViewport(viewport(x=0.25, width = 0.5, height = 0.8, angle = 0,#
+name = "leftvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+pushViewport(viewport(x=0.75, width = 0.5, height = 0.8, angle = 0,#
+name = "rightvp"))#
+#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p2, newpage=F)
+setwd('/Users/baptiste/Documents/R')
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+df2 <- data.frame(x <- rnorm(100), y=runif(x))#
+#
+p <- xyplot(y~x, data=df)#
+p2 <- xyplot(y~x, data=df2)#
+pushViewport(viewport(x=0.25, width = 0.5, height = 0.8, angle = 0,#
+name = "leftvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+pushViewport(viewport(x=0.75, width = 0.5, height = 0.8, angle = 0,#
+name = "rightvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p2, newpage=F)#
+upViewport()
+#
+library(grid)#
+library(ggplot2)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+df2 <- data.frame(x <- rnorm(100), y=runif(x))#
+#
+p <- qplot(x,y, data=df)#
+p2 <- qplot(x,y, data=df2)#
+#
+pushViewport(viewport(x=0.25, width = 0.5, height = 0.8, angle = 0,#
+name = "leftvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+pushViewport(viewport(x=0.75, width = 0.5, height = 0.8, angle = 0,#
+name = "rightvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p2, newpage=F)#
+upViewport()
+?grid.layout
+setwd('/Users/baptiste/Documents/R/photonics/pkg/mie/R')
+`Mie.ab` <-#
+function(m = 1,x = 1){#
+#
+ z <- m * x#
+ nmax <- round(Re(2+x+4*x^(1/3)))#
+ nmx <- round(max(nmax,Mod(z)) + 16)#
+ n <- (1:nmax)#
+ nu <- (n + 0.5)#
+#
+ sx <- sqrt(0.5 * pi * x)#
+ px <- sx * besselJ(x,nu)#
+ p1x <- c(sin(x), px[1:(nmax-1)])#
+#
+ chx <- -sx * besselY(x,nu)#
+ ch1x <- c(cos(x), chx[1:(nmax-1)])#
+#
+ gsx <- px-1i*chx#
+ gs1x <- p1x-1i*ch1x#
+#
+ dnx <- rep(0i,nmx)#
+ #
+ #
+ for (j in seq(from=nmx,to=2,by=-1) ){ #
+ dnx[j-1] <- j/z-1/(dnx[j]+j/z)#
+ }#
+#
+ dn <- dnx[n] #
+ da <- dn/m + n/x#
+ db <- m * dn + n/x#
+#
+ an <- (da * px - p1x) / (da * gsx - gs1x)#
+ bn <- (db * px - p1x) / (db * gsx - gs1x)#
+#
+ result <- cbind(an, bn)#
+ return(as.data.frame(result))#
+}
+besselJ(1:3,1)
+besselJ(1:3,1:3)
+besselJ(1:3,rep(1:3,each=3))
+besselJ(1:3,2)
Added: pkg/mie/R/Eint.r
===================================================================
--- pkg/mie/R/Eint.r (rev 0)
+++ pkg/mie/R/Eint.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,25 @@
+Eint <- function(cd, k, r, phi=0, theta, nmax = length(cd$cn), pt=Mie.pt(cos(theta),nmax)){
+
+ kr <- k*r
+ j <- besselJ(kr, 0:(nmax+1))
+ Es.r <- 0
+ Es.theta <- 0
+ Es.phi <- 0
+
+ for(n in 1:nmax){
+ cn <- cd$cn[n]
+ dn <- cd$dn[n]
+ pn <- pt[n, 1]
+ taun <- pt[n, 2]
+ En <- (1i)^n* (2*n+1)/(n*(n+1))
+ Es.r.n <- En*(-1i*dn*cos(phi)*n*(n+1)*sin(theta)*pn*j[n+1]/kr)
+ Es.theta.n <- En*(cn*cos(phi)*pn*j[n+1] - 1i*dn*cos(phi)*taun * ((n*j[n] - (n+1)*j[n+2])/(2*n+1) + j[n+1]/kr))
+ Es.phi.n <- En*(-cn*sin(phi)*taun*j[n+1] + 1i*dn*sin(phi)*pn * ((n*j[n] - (n+1)*j[n+2])/(2*n+1) + j[n+1]/kr))
+ Es.r <- Es.r + Es.r.n
+ Es.theta <- Es.theta + Es.theta.n
+ Es.phi <- Es.phi + Es.phi.n
+ }
+
+ Esca <- c(Es.r=Es.r, Es.theta=Es.theta, Es.phi=Es.phi)
+ Esca
+}
Added: pkg/mie/R/Escat.r
===================================================================
--- pkg/mie/R/Escat.r (rev 0)
+++ pkg/mie/R/Escat.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,26 @@
+Escat <- function(ab, k, r, phi=0, theta, nmax = length(ab$an), pt=Mie.pt(cos(theta),nmax)){
+
+ kr <- k*r
+ hn <- hankel(0:(nmax+1), kr)
+
+ Es.r <- 0
+ Es.theta <- 0
+ Es.phi <- 0
+
+ for(n in 1:nmax){
+ an <- ab$an[n]
+ bn <- ab$bn[n]
+ pn <- pt[n, 1]
+ taun <- pt[n, 2]
+ En <- (1i)^n* (2*n+1)/(n*(n+1))
+ Es.r.n <- En*(1i*an*cos(phi)*n*(n+1)*sin(theta)*pn*hn[n+1]/kr)
+ Es.theta.n <- En*(-bn*cos(phi)*pn*hn[n+1] + 1i*an*cos(phi)*taun * ((n*hn[n] - (n+1)*hn[n+2])/(2*n+1) + hn[n+1]/kr))
+ Es.phi.n <- En*(bn*sin(phi)*taun*hn[n+1] - 1i*an*sin(phi)*pn * ((n*hn[n] - (n+1)*hn[n+2])/(2*n+1) + hn[n+1]/kr))
+ Es.r <- Es.r + Es.r.n
+ Es.theta <- Es.theta + Es.theta.n
+ Es.phi <- Es.phi + Es.phi.n
+ }
+
+ Esca <- c(Es.r=Es.r, Es.theta=Es.theta, Es.phi=Es.phi)
+ Esca
+}
Added: pkg/mie/R/Mie.R
===================================================================
--- pkg/mie/R/Mie.R (rev 0)
+++ pkg/mie/R/Mie.R 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,49 @@
+`Mie` <-
+function(m, x){
+
+if (x == 0){ # To avoid a singularity at x=0
+ result <- list(qext = 0, qsca = 0, qabs = 0)
+ return(as.data.frame(result))
+}
+ # This is the normal situation
+
+ nmax <- round(2 + x + 4 * x^(1/3))
+ n1 <- nmax-1
+
+ n <- (1:nmax)
+ cn <- 2*n+1
+ c1n <- n * (n+2) / (n+1)
+ c2n <- cn / n / (n+1)
+
+ x2 <- x^2
+
+ ab <- Mie.ab(m,x)
+
+ anp <- Re(ab$an)
+ anpp <- Im(ab$an)
+
+ bnp <- Re(ab$bn)
+ bnpp <- Im(ab$bn)
+
+ g1 <- matrix(0,nrow=nmax,ncol=4) #[0; 0; 0; 0]; # displaced numbers used for
+
+ g1[1:n1,] <- cbind(anp[2:nmax], anpp[2:nmax], bnp[2:nmax], bnpp[2:nmax])
+
+ dn <- cn * (anp+bnp)
+
+ q <- sum(dn)
+
+ qext <- 2 * q / x2
+
+ en <- cn * (anp^2 + anpp^2 + bnp^2 + bnpp^2)
+
+ q <- sum(en)
+
+ qsca <- 2 * q / x2
+
+ qabs <- qext - qsca
+
+result <- c(qext, qabs, qsca)
+return(result)
+}
+
Added: pkg/mie/R/Mie.S.r
===================================================================
--- pkg/mie/R/Mie.S.r (rev 0)
+++ pkg/mie/R/Mie.S.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,33 @@
+Mie.S <- function(m, x, u){
+
+
+# Computation of Mie Scattering functions S1 and S2
+# for complex refractive index m=m'+im",
+# size parameter x=k0*a, and u=cos(scattering angle),
+# where k0=vacuum wave number, a=sphere radius;
+# s. p. 111-114, Bohren and Huffman (1983) BEWI:TDD122
+# C. Mätzler, May 2002
+
+nmax <- round(2+x+4*x^(1/3))
+
+ab <- Mie.ab(m,x)
+
+an <- ab$an
+bn <- ab$bn
+
+pt <- Mie.pt(u,nmax)
+
+pin <- pt$p
+tin <- pt$tau
+
+n <- 1:nmax
+n2 <- (2*n+1)/(n*(n+1))
+
+pin <- n2*pin
+tin <- n2*tin
+
+S1 <- (an*pin+bn*tin)
+S2 <- (an*tin+bn*pin)
+
+as.data.frame(cbind(S1, S2))
+}
\ No newline at end of file
Added: pkg/mie/R/Mie.ab.R
===================================================================
--- pkg/mie/R/Mie.ab.R (rev 0)
+++ pkg/mie/R/Mie.ab.R 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,37 @@
+`Mie.ab` <-
+function(m = 1,x = 1){
+
+ z <- m * x
+ nmax <- round(Re(2+x+4*x^(1/3)))
+ nmx <- round(max(nmax,Mod(z)) + 16)
+ n <- (1:nmax)
+ nu <- (n + 0.5)
+
+ sx <- sqrt(0.5 * pi * x)
+ px <- sx * besselJ(x,nu)
+ p1x <- c(sin(x), px[1:(nmax-1)])
+
+ chx <- -sx * besselY(x,nu)
+ ch1x <- c(cos(x), chx[1:(nmax-1)])
+
+ gsx <- px-1i*chx
+ gs1x <- p1x-1i*ch1x
+
+ dnx <- rep(0i,nmx)
+
+ # Computation of Dn(z) according to (4.89) of B+H (1983)
+ for (j in seq(from=nmx,to=2,by=-1) ){
+ dnx[j-1] <- j/z-1/(dnx[j]+j/z)
+ }
+
+ dn <- dnx[n] # Dn(z), n=1 to nmax
+ da <- dn/m + n/x
+ db <- m * dn + n/x
+
+ an <- (da * px - p1x) / (da * gsx - gs1x)
+ bn <- (db * px - p1x) / (db * gsx - gs1x)
+
+ result <- cbind(an, bn)
+ return(as.data.frame(result))
+}
+
Added: pkg/mie/R/Mie.cd.r
===================================================================
--- pkg/mie/R/Mie.cd.r (rev 0)
+++ pkg/mie/R/Mie.cd.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,43 @@
+`Mie.cd` <-
+function(m = 1,x = 1){
+
+ z <- m * x
+ nmax <- round(Re(2+x+4*x^(1/3)))
+ nmx <- round(max(nmax,Mod(z)) + 16)
+ n <- (1:nmax)
+ nu <- (n + 0.5)
+ msq <- m^2
+
+ cnx <- rep(0+0i, nmx)
+ for (j in seq(nmx, 2, by=-1)){
+ cnx[j-1] <- j-z*z/(cnx[j]+j)
+ }
+ cnn <- cnx[n]
+
+ sqx <- sqrt(0.5*pi/x)
+ sqz <- sqrt(0.5*pi/z)
+ bx <- besselJ(x, nu)*sqx
+ bz <- 1/besselJ(x, nu)/sqz
+ yx <- besselY(x, nu)*sqx
+ hx <- bx + 1i*yx
+
+ b1x <- c(sin(x)/x, bx[1:(nmax-1)])
+ y1x <- c(-cos(x)/x, yx[1:(nmax-1)])
+ h1x <- b1x+1i*y1x
+ ###
+
+ ax <- x*b1x-n*bx
+ ahx <- x*h1x-n*hx
+
+ nenn1 <- msq*ahx-hx*cnn
+
+ nenn2 <- ahx-hx*cnn
+
+ cn <- bz*(bx*ahx-hx*ax)/nenn2
+
+ dn <- bz*m*(bx*ahx-hx*ax)/nenn1
+
+ result <- cbind(cn, dn)
+ return(as.data.frame(result))
+}
+
Added: pkg/mie/R/Mie.coated.R
===================================================================
--- pkg/mie/R/Mie.coated.R (rev 0)
+++ pkg/mie/R/Mie.coated.R 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,56 @@
+`Mie.coated` <-
+function(m1, m2, x, y){
+
+if (x==y) { # To reduce computing time
+ result <- Mie(m1,y)
+ return(result)
+}
+
+if (x==0) { # To avoid a singularity at x=0
+ result <- Mie(m2,y)
+ return(result)
+}
+
+if (m1==m2) { # To reduce computing time
+ result <- Mie(m1,y)
+ return(result)
+}
+
+# This is now the normal situation
+
+ nmax <- round(2 + y + 4 * y^(1/3))
+ n1 <- nmax-1
+
+ n <- (1:nmax)
+ cn <- 2*n+1
+ c1n <- n * (n+2) / (n+1)
+ c2n <- cn / n / (n+1)
+
+ y2 <- y^2
+
+ ab <- Mie.coated.ab(m1, m2, x, y)
+
+ anp <- Re(ab$an)
+ anpp <- Im(ab$an)
+
+ bnp <- Re(ab$bn)
+ bnpp <- Im(ab$bn)
+ g1 <- matrix(0,nrow=nmax,ncol=4)
+
+ g1[1:n1,] <- cbind(anp[2:nmax], anpp[2:nmax], bnp[2:nmax], bnpp[2:nmax])
+
+ dn <- cn * (anp+bnp)
+ q <- sum(dn)
+
+ qext <- 2 * q / y2
+
+ en <- cn * (anp^2 + anpp^2 + bnp^2 + bnpp^2)
+
+ q <- sum(en)
+ qsca <- 2 * q / y2
+ qabs <- qext - qsca
+
+ result <- c(qext, qabs, qsca)
+ return(result)
+}
+
Added: pkg/mie/R/Mie.coated.ab.R
===================================================================
--- pkg/mie/R/Mie.coated.ab.R (rev 0)
+++ pkg/mie/R/Mie.coated.ab.R 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,100 @@
+`Mie.coated.ab` <-
+function(m1 = 1, m2 =1, x = 1, y=1){
+
+m <- m2 / m1
+
+# z <- m * x
+# nmax <- round(Re(2+x+4*x^(1/3)))
+# nmx <- round(max(nmax,Mod(z)) + 16)
+# n <- (1:nmax)
+# nu <- (n + 0.5)
+
+
+u <- m1 * x
+v <- m2 * x
+w <- m2 * y # The arguments of Bessel Functions
+
+nmax <- round(Re(2+y+4*y^(1/3))) # The various nmax values
+mx <- max(Mod(m1*y),Mod(m2*y))
+nmx <- round(max(nmax,mx) + 16)
+nmax1 <- nmax-1
+n <- (1:nmax)
+
+# Computation of Dn(z), z=u,v,w according to (4.89) of B+H (1983)
+
+z <- u
+dnx <- rep(0i,nmx)
+for (j in seq(from=nmx,to=2,by=-1) ){
+ dnx[j-1] <- j/z-1/(dnx[j]+j/z)
+}
+dnu <- dnx[n]
+
+z <- v
+for (j in seq(from=nmx,to=2,by=-1) ){
+ dnx[j-1] <- j/z-1/(dnx[j]+j/z)
+}
+dnv <- dnx[n]
+
+z <- w
+for (j in seq(from=nmx,to=2,by=-1) ){
+ dnx[j-1] <- j/z-1/(dnx[j]+j/z)
+}
+dnw <- dnx[n]
+
+# Computation of Psi, Chi and Gsi Functions and their derivatives
+
+nu <- (n+0.5)
+
+sv <- sqrt(0.5*pi*v)
+pv <- sv * besselJ(v, nu)
+
+sw <- sqrt(0.5*pi*w)
+pw <- sw*besselJ(w, nu)
+
+sy <- sqrt(0.5*pi*y)
+py <- sy*besselJ(y, nu)
+
+p1y <- c(sin(y), py[1:(nmax-1)])
+chv <- -sv * besselY(v,nu)
+chw <- -sw * besselY(w,nu)
+chy <- -sy * besselY(y,nu)
+
+ch1y <- c(cos(y), chy[1:(nmax-1)])
+
+gsy <- py-1i*chy
+gs1y <- p1y-1i*ch1y
+
+# Computation of U, V, F Functions, avoiding products of Riccati-Bessel Fcts.
+
+uu <- m*dnu-dnv
+vv <- dnu/m-dnv
+
+fv <- pv/chv
+fw <- pw/chw
+
+ku1 <- uu*fv/pw
+kv1 <- vv*fv/pw
+
+ku2 <- uu*(pw-chw*fv) + (pw/pv)/chv
+kv2 <- vv*(pw-chw*fv)+(pw/pv)/chv
+
+dns1 <- ku1/ku2
+gns1 <- kv1/kv2
+
+# Computation of Dn_Schlange, Gn_Schlange
+
+dns <- dns1+dnw
+gns <- gns1+dnw
+
+a1 <- dns/m2+n/y
+b1 <- m2*gns+n/y
+
+# an and bn
+
+an <- (py*a1-p1y)/(gsy*a1-gs1y)
+
+bn <- (py*b1-p1y)/(gsy*b1-gs1y)
+result <- cbind(an, bn)
+return(as.data.frame(result))
+}
+
Added: pkg/mie/R/Mie.coated.spectra.R
===================================================================
--- pkg/mie/R/Mie.coated.spectra.R (rev 0)
+++ pkg/mie/R/Mie.coated.spectra.R 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,31 @@
+`Mie.coated.spectra` <-
+function(lambda, core.n, inner.radius = 0.04 ,
+ layer.thickness = 0.005, layer.n = 1.5,
+ medium = 1.333){
+
+ if (length(layer.n) == 1 & length(core.n) == 1) {
+ stop("both layer.n and core.epsilon are scalars, run Mie.coated instead")
+ }
+ if (length(layer.n) == 1 ) {
+ layer.n <- rep(layer.n, length = length(core.n)) # case a scalar is given (non-dispersive)
+ }
+ if (length(core.n) == 1 ) {
+ core.n <- rep(core.n, length = length(layer.n)) # case a scalar is given (non-dispersive)
+ }
+
+ outer.radius <- inner.radius + layer.thickness
+ area <- pi * outer.radius^2
+ k0 <- 2 * pi / lambda
+ k <- k0 * medium
+
+ x <- k * inner.radius
+ y <- k * outer.radius
+
+ m1 <- core.n / medium
+ m2 <- layer.n / medium
+
+ efficiencies <- mapply(Mie.coated, m1=m1, m2 = m2, x = x, y = y)
+
+ return(area*t(efficiencies)) # cross sections, this is a matrix that can be used with mapply to sweep over params
+}
+
Added: pkg/mie/R/Mie.pt.r
===================================================================
--- pkg/mie/R/Mie.pt.r (rev 0)
+++ pkg/mie/R/Mie.pt.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,26 @@
+Mie.pt <- function (u,nmax){
+#
+# pi_n and tau_n, -1 <= u <= 1, n1 integer from 1 to nmax
+# angular functions used in Mie Theory
+# Bohren and Huffman (1983), p. 94 - 95
+
+p <- rep(1, 3*nmax)
+tau <- rep(u, 3*nmax)
+p[2] <- 3*u
+tau[2] <- 3*cos(2*acos(u))
+
+for (n1 in (3:nmax)){
+
+ p1 <- (2*n1-1)/(n1-1)*p[n1-1]*u
+ p2 <- n1/(n1-1)*p[n1-2]
+ p[n1] <- p1-p2
+
+ t1 <- n1*u*p[n1]
+ t2 <- (n1+1)*p[n1-1]
+ tau[n1] <- t1-t2
+
+}
+
+as.data.frame(cbind(p, tau))
+}
+
Added: pkg/mie/R/Mie.shell.r
===================================================================
--- pkg/mie/R/Mie.shell.r (rev 0)
+++ pkg/mie/R/Mie.shell.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,69 @@
+# library(constants)
+# library(baptMisc)
+# library(scatteringR)
+# clr()
+# data(AuJChristy)
+#
+# cel <- 3e8
+#
+# plot.material("AuJChristy", range=c(0.4, 1.2), plot=F)$fine.sample->gold
+# lambda <- gold$lambda
+# epsilon1 <- 5.44
+# epsilon2 <- gold$epsilon
+# epsilon3 <- 1.78
+#
+# k0 <- 2*pi/lambda #* sqrt(epsilon3)
+# k <- 2*pi/lambda * sqrt(epsilon3)
+#
+# omega <- k0 * cel
+#
+# m12 <- (epsilon1 / epsilon3)
+# m22 <- (epsilon2 / epsilon3)
+# m1 <- sqrt(m12)
+# m2 <- sqrt(m22)
+# ma <- (m22 - m12) / (3*m22 + 2* m12)
+# mb <- (m22 - m12) / (1 + 2*m1)
+#
+#
+# a <- 20.0e-3 # inner radius
+# b <- 25e-3 # outer radius
+#
+# x <- a * omega / cel
+# y <- b * omega / cel
+#
+# a1a <- 80*(m12 + 2*m22) - 8*(m12 - m22)*(m12 + 10 * m22) *x^2 +
+# 2*m22*(2*m12^2 + 7 * m12*m22 - 10*m22^2)*x^4 -
+# m12 * m22^2*(3*m12 - 4*m22)*x^6
+#
+# a1b <- (m12 - m22) * (m12*m22*x^4 - 5*(m12 + m22)*x^2 + 50) *x^3
+#
+# a1ay <- (m22 - 1) * y^3 * (m22*y^4 - 5*(m22+1) + 50)
+#
+# a1by <- 80*(1+m22) + 8*(m22 - 1)*(10*m22+1)*y^2 -
+# 2*m22*(10*m22^2 - 7 * m22 - 2)*y^4 + m2^4*(4*m22 - 3)*y^6
+#
+# a1cy <- -60i*(m22 + 2) + (m22 - 1)*y^2 * (6i*(10+m22) - 40*y) -
+# m22 * y^4* (3i*(m22 - 4) - 4*(m22 - 2)*y)
+#
+# a1dy <- 24i*(m22 - 1)*(2 + (m22 + 1)*y^2) +
+# 16*(2*m22 + 1)*y^3 -
+# 6i*m22*(m2^4 + 5*m22 -2)*m22 * y^4 + 8*m22 *
+# (2*m22 - 1)*y^5 + 3*m2^4 *(m22 - 3)*y^6*1i -
+# 2*m2^4*(2*m22 - 3)*y^7
+#
+# a1 <- (4/5 * y^3 * a1a * a1ay + a1b * a1by) / (a1a * y^3 * a1cy - 4 * a1b * a1dy)
+# a2 <- 1i/15*(y^5*(1-m22) + ma*x^5 * (2 + 3* m22)) / ((3 + 2*m22) + 6*ma*x^5 /y^5 *(1-m22))
+# b1 <- 1i/45*(y^5*(1-m22) +3*m1*mb*x^5) / (1 - 1/6 * y^2*(m22 + 1))
+#
+# term.1 <- (12/5 * y^3 * a1a * a1ay + a1b * a1by) / (a1a * y^3 * a1cy - 4 * a1b * a1dy)
+# term.2 <- 5i/15*(y^5*(1-m22) + ma*x^5*(2+3*m22)) / ((3 + 2*m22) + 6*ma*x^5 /y^5*(1-m22))
+# term.3 <- 1i/15*(y^5*(1-m22) +3*m1*mb*x^5) / (1 - 1/6 * y^2*(m22 + 1))
+#
+# qsca <- 2*pi/k^2 *((term.1) + (term.2) + (term.3))
+#
+# qsca2 <- 2*pi/k^2 *(3*(abs(a1)^2 + abs(b1)^2) + 5*abs(a2)^2)
+# qext <- 2*pi/k^2 * Re(3*((a1) + (b1)) + 5*(a2))
+#
+# plot(lambda, qsca2)
+# plot(lambda, qext)
+# # lines(lambda, qsca)
\ No newline at end of file
Added: pkg/mie/R/Mie.spectra.R
===================================================================
--- pkg/mie/R/Mie.spectra.R (rev 0)
+++ pkg/mie/R/Mie.spectra.R 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,17 @@
+`Mie.spectra` <-
+function( lambda, epsilon, radius = 0.04 , medium = 1.333, efficiency = FALSE ){
+
+ area <- pi * (radius)^2
+ k0 <- 2 * pi / lambda
+ k <- k0 * medium
+ x <- k * radius
+ m <- sqrt(epsilon) / medium
+
+ efficiencies <- mapply(Mie, m = m, x = x)
+
+ if(efficiency == TRUE)
+ return(t(efficiencies)) # um^2
+ else
+ return(area*t(efficiencies)) # cross sections
+}
+
Added: pkg/mie/R/hankel.r
===================================================================
--- pkg/mie/R/hankel.r (rev 0)
+++ pkg/mie/R/hankel.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,9 @@
+hankel <- function(n, x){
+ sqx <- sqrt(0.5*pi/x)
+ nu <- (n+0.5)
+ bx <- besselJ(x,nu)*sqx
+ yx <- besselY(x,nu)*sqx
+
+ hx <- bx+1i*yx
+ hx
+}
\ No newline at end of file
Added: pkg/mie/R/permittivity.r
===================================================================
--- pkg/mie/R/permittivity.r (rev 0)
+++ pkg/mie/R/permittivity.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,5 @@
+permittivity <- function(material=AuJChristy, wavelength=632.8, plot=F, ...){
+ plot.material(material, plot=plot, ...)->tmp
+ predict(tmp$spline.real,wavelength*1e-9)$y +
+ 1i * predict(tmp$spline.imag,wavelength*1e-9)$y
+}
\ No newline at end of file
Added: pkg/mie/R/plot.material.R
===================================================================
--- pkg/mie/R/plot.material.R (rev 0)
+++ pkg/mie/R/plot.material.R 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,51 @@
+plot.material <- function (material = AuJChristy, all = FALSE, plot = FALSE,
+ range.lambda = c(0.4, 1.2), fine.n = 200 )
+{
+ lambda <- material$lambda * 1e+06
+ if (all == FALSE) {
+ limits <- c(which.min(abs(lambda - range.lambda[1])),
+ which.min(abs(lambda - range.lambda[2])))
+ if (limits[1] == limits[2])
+ stop("too few unique points")
+ lambda <- lambda[limits[1]:limits[2]]
+ epsr <- material$epsr[limits[1]:limits[2]]
+ epsi <- material$epsi[limits[1]:limits[2]]
+ }
+ else {
+ epsr <- material$epsr
+ epsi <- material$epsi
+ }
+
+ raw.data <- data.frame(lambda = lambda, epsilon = epsr + 1i* epsi, nk = sqrt(epsr + 1i* epsi))
+ with(raw.data,smooth.spline(lambda,Re(epsilon)))->spline.epsr
+ with(raw.data,smooth.spline(lambda,Im(epsilon)))->spline.epsi
+
+ lambda.fine <- seq(from=range(raw.data$lambda)[1],to=range(raw.data$lambda)[2],length=fine.n)
+ epsilon.fine <- predict(spline.epsr,lambda.fine)$y + 1i * predict(spline.epsi,lambda.fine)$y
+
+ fine.sample <- data.frame(lambda= lambda.fine,epsilon = epsilon.fine, nk = sqrt(epsilon.fine))
+
+if (plot == TRUE){
+ layout(rbind(c(1, 1), c(2, 3), c(4, 5)), widths = c(1, 1),
+ heights = c(0.8, 1, 1))
+ par(mar = c(1, 5, 4, 2), pty = "m", xpd = NA)
+ plot(c(-1, 1), c(-1, 1), t = "n", bty = "n", xaxt = "n",
+ yaxt = "n", xlab = "", ylab = "")
+ com <- gsub(";", "", material$comment)
+ text(0, 0, paste(com[-c(length(com) - 1, length(com))], collapse = "\n"),
+ cex = 1.2)
+
+ par(mar = c(5, 5, 1, 1), pty = "m", xpd = FALSE)
+ with(raw.data,plot(lambda, Re(epsilon), t = "p", ylab =~epsilon[r], xlab=~"wavelength / "*mu*m))
+ with(fine.sample,lines(lambda,Re(epsilon)))
+ with(raw.data,plot(lambda, Im(epsilon), t = "p", ylab =~epsilon[i], xlab=~"wavelength / "*mu*m))
+ with(fine.sample,lines(lambda,Im(epsilon)))
+ with(raw.data,plot(lambda, Re(nk), t = "p", ylab = "n", xlab=~"wavelength / "*mu*m))
+ with(fine.sample,lines(lambda,Re(nk)))
+ with(raw.data,plot(lambda, Im(nk), t = "p", ylab = "k", xlab=~"wavelength / "*mu*m))
+ with(fine.sample,lines(lambda,Im(nk)))
+}
+
+ invisible(list(raw.data = raw.data, fine.sample = fine.sample,
+ spline.real=spline.epsr, spline.imag=spline.epsi))
+}
Property changes on: pkg/mie/R/plot.material.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/mie/data/AuJChristy.rda
===================================================================
(Binary files differ)
Property changes on: pkg/mie/data/AuJChristy.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/mie/inst/.Rhistory
===================================================================
--- pkg/mie/inst/.Rhistory (rev 0)
+++ pkg/mie/inst/.Rhistory 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,2280 @@
+#
+data(permittivitiesLuxpop) #
+plotMaterial("Cr",all=FALSE,range=c(0.5,0.8))
+#
+data(permittivitiesLuxpop) #
+plot.material("Cr",all=FALSE,range=c(0.5,0.8))
+?constants
+#
+data(luxpop) #
+plot.material("Cr",all=FALSE,range=c(0.5,0.8))
+baseUrl="http://www.luxpop.com/Material/"#
+material="Cr.nk"#
+readLines(url(paste(baseUrl,material,sep=""), open = "", blocking = TRUE, encoding = getOption("encoding")))->file#
+grep(";",file,value=TRUE)->comment#
+comment#
+read.table(textConnection(file),skip=length(comment))->values#
+head(values)#
+names(values)<-c("l","n","k")#
+attach(values)#
+nk2epsilon(list(n=n,k=k))->valuesEpsilon#
+valuesEpsilon$lambda<-l/1e4 #
+detach(values)#
+attach(valuesEpsilon)#
+#
+graphics.off()#
+par(mfrow=c(2,2),mar=c(5, 5, 1, 1),bty="n",pty="m")#
+cutePlot(lambda,epsr)#
+cutePlot(lambda,epsi,yl=expression(paste(epsilon[i])))#
+cutePlot(lambda,values$n,yl=expression(paste("n")))#
+cutePlot(lambda,values$k,yl=expression(paste("k")))
+version()
+help()
+sessionInfo
+sessionInfo()
+#
+#
+#
+#
+#
+squareGrob <- grob(height=unit(.1, "npc"),#
+ cl="squareGrob")#
+#
+#
+#
+#
+widthDetails.squareGrob <- function(x) {#
+ convertUnit(unit(.1, "npc"),#
+ "npc",#
+ "y", "dimension",#
+ "x", "dimension")#
+}#
+#
+#
+#
+lay1 <- grid.layout(10, 2,#
+ widths=unit.c(unit(1, "null"),#
+ grobWidth(squareGrob)),#
+ heights=unit(.1, "npc"))#
+#
+... Once that squareGrob is defined, you can use it anywhere, so that's#
+a one-off cost. The following code shows the resulting layout ...#
+#
+#
+pushViewport(viewport(layout=lay1))#
+#
+pushViewport(viewport(layout.pos.col=1))#
+grid.rect(gp=gpar(fill="grey"))#
+popViewport()#
+#
+for (i in 1:10) {#
+ pushViewport(viewport(layout.pos.col=2,#
+ layout.pos.row=i))#
+ grid.rect(gp=gpar(fill="light grey"))#
+ popViewport()#
+}
+?spectra
+library(constants)
+spectra()->test
+test
+library(spectra)
+spectra2excel(test,"test.txt")
+getWavelength(spectra())
+source("http://r.research.att.com/survey.R")
+system("mate /var/folders/yF/yFY+OcggEhuLb4YeNRUrPU+++TI/-Tmp-//RtmptFKYkp/downloaded_packages")
+system("open /var/folders/yF/yFY+OcggEhuLb4YeNRUrPU+++TI/-Tmp-//RtmptFKYkp/downloaded_packages")
+library(baptMisc)
+library(SPPr)
+?SPPr
+graphics.off()#
+library(SPPr)#
+#
+data(scanSPP)#
+data(parameters)#
+head(scanSPP)#
+#
+parameters$zeroShift <- -0.4#
+parameters$thickness[3] <- 0#
+#
+dataNorm <- data2norm(scanSPP,parameters)#
+head(dataNorm)#
+#
+plotData(dataNorm, int=TRUE, col=1, ylim=c(0, 1))#
+criticalEdge(parameters)#
+#
+layer2fit <- 1#
+pIni <- c(-18,0.7, 45e-9, 1.05)#
+#
+plotModel(pIni, int=TRUE, col=2)#
+#
+result <- optim(pIni, objF)#
+#
+yFin <- fitModel(result$par)#
+yErr<-abs(yFin-dataNorm$N)#
+#
+plotModel(result$par,col=3)
+result
+graphics.off()#
+library(SPPr)#
+#
+data(scanSPP)#
+data(parameters)#
+head(scanSPP)#
+#
+parameters$zeroShift <- -0.4#
+parameters$thickness[3] <- 0#
+#
+dataNorm <- data2norm(scanSPP,parameters)#
+head(dataNorm)#
+#
+plotData(dataNorm, int=TRUE, col=1, ylim=c(0, 1))#
+criticalEdge(parameters)#
+#
+layer2fit <- 1#
+pIni <- c(-18,0.7, 45, 1.05)#
+#
+plotModel(pIni, int=TRUE, col=2)#
+#
+result <- optim(pIni, objF)#
+#
+yFin <- fitModel(result$par)#
+yErr<-abs(yFin-dataNorm$N)#
+#
+plotModel(result$par,col=3)
+#
+data(scanSPP)#
+head(scanSPP)
+x.axis <- seq(from=0.5, to=4.5, length.out=13112)
+which.min(abs(x.axis - 3.2))
+?lattice
+?xyplot
+brewer.pal
+xyplot
+library(fortune)
+library(fortunes)
+fortune()
+library(spatstat)
+? spatstat
+citation()
+citation("spatstat")
+#
+#
+xrange=c(-15,15)#
+yrange=c(0,16)#
+plot(0,xlim=xrange,ylim=yrange,type='n')#
+#
+#
+yr=seq(yrange[1],yrange[2],len=50)#
+offsetFn=function(y){2*sin(0+y/3)}#
+offset=offsetFn(yr)#
+leftE = function(y){-10-offsetFn(y)}#
+rightE = function(y){10+offsetFn(y)}#
+#
+#
+xp=c(leftE(yr),rev(rightE(yr)))#
+yp=c(yr,rev(yr))#
+polygon(xp,yp,col="#
+#
+#
+#
+#
+h=9#
+#
+#
+xt=seq(0,rightE(h),len=100)#
+yt=log(1+log(1+log(xt+1)))#
+yt=yt-min(yt)#
+yt=h*yt/max(yt)#
+#
+x=c(leftE(h),rightE(h),rev(xt),-xt)#
+y=c(h,h,rev(yt),yt)#
+polygon(x,y,col="red",border=NA)
+#
+library(Hmisc)#
+sexyNo<-function(){#
+ for(i in (1:20)){#
+ S<-paste(makeNstr(" ",i),"0",makeNstr(" ",(40-2*i)),"<--1\n",sep="")#
+ cat(S); Sys.sleep(1)#
+ }#
+ cat(paste(makeNstr(" ",21),"0~1\n",sep="")); Sys.sleep(2)#
+ for(i in (1:9)){#
+ cat(paste(makeNstr(" ",22),"|\n",sep="")); Sys.sleep(2)#
+ }#
+ cat(" ***!!! 0.1 !!!***\n")#
+}#
+#
+sexyNo()
+plotCI
+?plotCI
+load("/Users/baptiste/Documents/PhD/DraftsPapers/disorderArrays/analysis/diagModel.rda")
+library(ggplot2)
+River.Mile <- c(215 ,202, 198, 190, 185, 179, 148, 119, 61)#
+Cu <- rnorm(9)#
+Fe <- rnorm(9)#
+Mg <- rnorm(9)#
+Ti <- rnorm(9)#
+Ir <- rnorm(9)#
+r <- data.frame(River.Mile, Cu, Fe, Mg, Ti, Ir)#
+#
+z <- melt.data.frame(r, id.var="River.Mile")#
+#
+#
+qplot(River.Mile, value, facets=(variable~.), data=z)
+#
+xyplot(value~River.Mile | variable, data=z)
+?plotmatrix
+data(mtcars)
+#
+plotmatrix(mtcars[, 1:3])
+plotmatrix
+#
+plotmatrix( data=z)
+library(zoo)#
+d<-(structure(c(1.39981554315924, 0.89196314359498, 0.407816250252697,#
+0.823496839063978, 1.14429021220358, 1.23971035967413, 0.960868900583432,#
+0.927685306209829, 1.22072345292821, 0.249842897450642, 1.00879641624694,#
+0.925372139878243, 0.317259909172362, 0.382677149697482), index =#
+structure(c(11808,#
+11869, 11961, 11992, 12084, 12173, 12265, 12418, 12600, 12631,#
+12753, 12996, 13057, 13149), class = "Date"), class = "zoo"))#
+#
+plot(d)
+library(ifs)
+?ifs
+#
+#
+#
+#
+#
+IFS.est <- IFS(y)#
+xx <- IFS.est$x#
+tt <- IFS.est$y#
+#
+ss <- pbeta(xx,.5,.1)#
+#
+#
+#
+ #
+par(mfrow=c(2,1)) #
+ #
+plot(ecdf(y),xlim=c(0,1),main="IFS estimator versus EDF")#
+lines(xx,ss,col="blue")#
+lines(xx,tt,col="red")
+#
+y<-rbeta(50,.5,.1)#
+#
+#
+#
+#
+IFS.est <- IFS(y)#
+xx <- IFS.est$x#
+tt <- IFS.est$y#
+#
+ss <- pbeta(xx,.5,.1)#
+#
+#
+#
+ #
+par(mfrow=c(2,1)) #
+ #
+plot(ecdf(y),xlim=c(0,1),main="IFS estimator versus EDF")#
+lines(xx,ss,col="blue")#
+lines(xx,tt,col="red")
+a=-0.7;b=-0.4 #
+Limits=c(-2,2)#
+MaxIter=60#
+cl=colours()#
+Step=seq(Limits[1],Limits[2],by=0.01)#
+PointsMatrix=array(0,dim=c(length(Step)*length(Step),3))#
+a1=0#
+plot(0,0,xlim=Limits,ylim=Limits,col="white")#
+#
+for(x in Step)#
+{#
+ for(y in Step)#
+ {#
+ n=0#
+ DIST=0#
+ x1=x;y1=y #
+ while(n<MaxIter & DIST<4)#
+ {#
+ newx=x1^2-y1^2+a#
+ newy=2*x1*y1+b#
+ DIST=newx^2+newy^2#
+ x1=newx;y1=newy#
+ n=n+1#
+ }#
+ if(DIST<4) colour=24 else colour=n*10#
+ #
+ a1=a1+1#
+ PointsMatrix[a1,]=c(x,y,colour)#
+ }#
+}#
+#
+plot(PointsMatrix[,1], PointsMatrix[,2], xlim=Limits, ylim=Limits, col=cl[PointsMatrix[,3]], pch=".")
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+p <- xyplot(y~x, data=df)#
+pushViewport(viewport(x=0.3, width = 0.5, height = 0.8, angle = 0,#
+name = "topvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+p <- xyplot(y~x, data=df)#
+pushViewport(viewport(x=0.3, width = 0.5, height = 0.8, angle = 0,#
+name = "topvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+grid.text("some text elsewhere", x=0.8, gp = gpar(col = rgb(43/255,#
+140/255, 190/255)))
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+df2 <- data.frame(x=rnorm(100), y=runif(x))#
+#
+p <- xyplot(y~x, data=df)#
+p2 <- xyplot(y~x, data=df2)#
+pushViewport(viewport(x=0.25, width = 0.5, height = 0.8, angle = 0,#
+name = "leftvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+pushViewport(viewport(x=0.75, width = 0.5, height = 0.8, angle = 0,#
+name = "rightvp"))#
+print(p2, newpage=F)
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+df2 <- data.frame(x <- rnorm(100), y=runif(x))#
+#
+p <- xyplot(y~x, data=df)#
+p2 <- xyplot(y~x, data=df2)#
+pushViewport(viewport(x=0.25, width = 0.5, height = 0.8, angle = 0,#
+name = "leftvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+pushViewport(viewport(x=0.75, width = 0.5, height = 0.8, angle = 0,#
+name = "rightvp"))#
+print(p2, newpage=F)
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+df2 <- data.frame(x <- rnorm(100), y=runif(x))#
+#
+p <- xyplot(y~x, data=df)#
+p2 <- xyplot(y~x, data=df2)#
+pushViewport(viewport(x=0.25, width = 0.5, height = 0.8, angle = 0,#
+name = "leftvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+pushViewport(viewport(x=0.75, width = 0.5, height = 0.8, angle = 0,#
+name = "rightvp"))#
+#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p2, newpage=F)
+setwd('/Users/baptiste/Documents/R')
+library(grid)#
+library(lattice)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+df2 <- data.frame(x <- rnorm(100), y=runif(x))#
+#
+p <- xyplot(y~x, data=df)#
+p2 <- xyplot(y~x, data=df2)#
+pushViewport(viewport(x=0.25, width = 0.5, height = 0.8, angle = 0,#
+name = "leftvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+pushViewport(viewport(x=0.75, width = 0.5, height = 0.8, angle = 0,#
+name = "rightvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p2, newpage=F)#
+upViewport()
+#
+library(grid)#
+library(ggplot2)#
+#
+df <- data.frame(x=rnorm(100), y=rnorm(100))#
+df2 <- data.frame(x <- rnorm(100), y=runif(x))#
+#
+p <- qplot(x,y, data=df)#
+p2 <- qplot(x,y, data=df2)#
+#
+pushViewport(viewport(x=0.25, width = 0.5, height = 0.8, angle = 0,#
+name = "leftvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p, newpage=F)#
+upViewport()#
+pushViewport(viewport(x=0.75, width = 0.5, height = 0.8, angle = 0,#
+name = "rightvp"))#
+grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+print(p2, newpage=F)#
+upViewport()
+#
+d = data.frame(a=1, b=2)#
+a = c("a", "b")#
+z = a
+d
+a
+z
+#
+subset(d, select=z)
+?subset
+#
+subset(d, select=a)
+#
+subset(d, select=b)
+#
+subset(d, select="a")
+aa =a
+#
+subset(d, select=aa)
+factorial <- function(n) {#
+ #
+ if (!is.loaded(symbol.For('factorial_R_wrapper'))) { dyn.load('factorial.so') } #
+ #
+ #
+ #
+ #
+ #
+ #
+ returned_data = .Fortran('factorial_r_wrapper', n=as.integer(n), result=integer(1)) #
+ #
+ return(returned_data$result) #
+ }
+factorial <- function(n) {#
+ #
+ if (!is.loaded(symbol.For('factorial_R_wrapper'))) { dyn.load('factorial.so') } #
+ #
+ #
+ #
+ #
+ #
+ #
+ returned_data = .Fortran('factorial_r_wrapper', n=as.integer(n), result=integer(1)) #
+ #
+ return(returned_data$result) #
+ }#
+factorial(1)
+?Fortran
+?.Fortran
+#
+#
+library(mie)#
+data(AuJChristy)#
+plot.material()$fine.sample->material#
+#
+radius <- 0.2#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))#
+#
+#
+#
+#
+lambda <- 1.2#
+#
+k0 <- 2*pi/lambda#
+#
+eps <- permittivity(wav=lambda)#
+n <- Re(sqrt(eps))#
+k <- Im(sqrt(eps))#
+#
+ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+#
+scale <- 2#
+N <- 100#
+xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+#
+cart2polar <- function(xy){#
+r <- sqrt(xy[, 1]^2+xy[, 2]^2) #
+angle <- atan(xy[, 2]/xy[, 1])#
+data.frame(r=r, angle=angle)#
+}#
+#
+rphi <- cart2polar(xy)#
+remove <- which(rphi$r <= radius)#
+#
+my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+int.field <- mapply(Eint, r=rphi$r[remove], theta=rphi$angle[remove], MoreArgs=list(cd=cd,k=n*k0, phi=0)) #
+xy$z <- 0#
+xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+xy$z[remove] <- c(apply(int.field, 2, function(x) Re(x%*%Conj(x))))#
+#
+levelplot(z~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)
+#
+#
+radius <- 0.2#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+#
+#
+library(mie)#
+data(AuJChristy)#
+plot.material()$fine.sample->material#
+#
+radius <- 0.2#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+#
+levelplot(z~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)
+#
+lambda <- 1.2#
+lambda <- 0.5#
+k0 <- 2*pi/lambda#
+#
+eps <- permittivity(wav=lambda)#
+n <- Re(sqrt(eps))#
+k <- Im(sqrt(eps))#
+#
+ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+#
+scale <- 2#
+N <- 100#
+xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+#
+cart2polar <- function(xy){#
+r <- sqrt(xy[, 1]^2+xy[, 2]^2) #
+angle <- atan(xy[, 2]/xy[, 1])#
+data.frame(r=r, angle=angle)#
+}#
+#
+rphi <- cart2polar(xy)#
+remove <- which(rphi$r <= radius)#
+#
+my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+int.field <- mapply(Eint, r=rphi$r[remove], theta=rphi$angle[remove], MoreArgs=list(cd=cd,k=n*k0, phi=0)) #
+xy$z <- 0#
+xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+xy$z[remove] <- c(apply(int.field, 2, function(x) Re(x%*%Conj(x))))#
+#
+levelplot(z~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)
+#
+xy$z <- 0#
+xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+#
+xy$z[remove] <- NA#
+#
+levelplot(z~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)
+#
+#
+library(mie)#
+data(AuJChristy)#
+plot.material()$fine.sample->material#
+#
+radius <- 0.2#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))#
+#
+#
+#
+#
+lambda <- 1.2#
+lambda <- 0.5#
+k0 <- 2*pi/lambda#
+#
+eps <- permittivity(wav=lambda)#
+n <- Re(sqrt(eps))#
+k <- Im(sqrt(eps))#
+#
+ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+#
+scale <- 2#
+N <- 40#
+xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+#
+cart2polar <- function(xy){#
+r <- sqrt(xy[, 1]^2+xy[, 2]^2) #
+angle <- atan(xy[, 2]/xy[, 1])#
+data.frame(r=r, angle=angle)#
+}#
+#
+rphi <- cart2polar(xy)#
+remove <- which(rphi$r <= radius)#
+#
+my.field12 <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=2*pi/1.2, phi=0))#
+my.field065 <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=2*pi/0.65, phi=0))#
+my.field05 <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=2*pi/0.5, phi=0))#
+#
+xy$z12 <- 0#
+xy$z05 <- 0#
+xy$z065 <- 0#
+xy$z12[-remove] <- c(apply(my.field12, 2, function(x) Re(x%*%Conj(x))))#
+xy$z05[-remove] <- c(apply(my.field05, 2, function(x) Re(x%*%Conj(x))))#
+xy$z065[-remove] <- c(apply(my.field065, 2, function(x) Re(x%*%Conj(x))))#
+#
+#
+#
+p12 <- levelplot(z12~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)#
+#
+p065 <- levelplot(z065~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)#
+#
+p05 <- levelplot(z05~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)
+p05
+library(grid)
+?grid.layout
+#
+myLay <- grid.layout(3, 1, widths=unit(1, "inches"),#
+ heights=unit(0.25, "npc"),#
+ just=just)#
+#
+pushViewport(viewport(layout=myLay)))
+#
+myLay <- grid.layout(3, 1, widths=unit(1, "inches"),#
+ heights=unit(0.25, "npc"))#
+#
+pushViewport(viewport(layout=myLay)))
+#
+myLay <- grid.layout(3, 1, widths=unit(1, "inches"),#
+ heights=unit(0.25, "npc"))#
+#
+pushViewport(viewport(layout=myLay))
+myLay <- grid.layout(3, 1, widths=unit(1, "inches"),#
+ heights=unit(0.25, "npc"))#
+#
+pushViewport(viewport(layout=myLay))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+pushViewport(viewport(layout.pos.col=1, name="v11"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ #
+ upViewport()
+myLay <- grid.layout(3, 1, widths=unit(1, "inches"),#
+ heights=unit(0.25, "npc"))#
+#
+pushViewport(viewport(layout=myLay))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+pushViewport(viewport(layout.pos.col=1, name="v11"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ #
+ upViewport()#
+ upViewport()
+library(lattice)#
+myLay <- grid.layout(3, 1, widths=unit(1, "inches"),#
+ heights=unit(0.25, "npc"))#
+#
+pushViewport(viewport(layout=myLay))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+pushViewport(viewport(layout.pos.col=1, name="v11"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p05, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=2, name="v12"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p065, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=3, name="v13"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p12, newpage=F)#
+ upViewport()
+#
+library(lattice)#
+myLay <- grid.layout(3, 3, widths=unit(1, "inches"),#
+ heights=unit(0.25, "npc"))#
+#
+pushViewport(viewport(layout=myLay))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+pushViewport(viewport(layout.pos.col=1, name="v11"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p05, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=2, name="v12"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p065, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=3, name="v13"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p12, newpage=F)#
+ upViewport()
+#
+library(lattice)#
+myLay <- grid.layout(1, 3, widths=unit(1, "inches"),#
+ heights=unit(0.25, "npc"))#
+#
+pushViewport(viewport(layout=myLay))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+pushViewport(viewport(layout.pos.col=1, name="v11"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p05, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=2, name="v12"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p065, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=3, name="v13"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p12, newpage=F)#
+ upViewport()
+#
+library(lattice)#
+myLay <- grid.layout(1, 3, widths=unit(2, "inches"),#
+ heights=unit(2, "inches"))#
+#
+pushViewport(viewport(layout=myLay))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+pushViewport(viewport(layout.pos.col=1, name="v11"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p05, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=2, name="v12"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p065, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=3, name="v13"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p12, newpage=F)#
+ upViewport()
+#
+library(lattice)#
+myLay <- grid.layout(1, 3, widths=unit(5, "inches"),#
+ heights=unit(2, "inches"))#
+#
+pushViewport(viewport(layout=myLay))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+pushViewport(viewport(layout.pos.col=1, name="v11"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p05, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=2, name="v12"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p065, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=3, name="v13"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p12, newpage=F)#
+ upViewport()
+#
+library(lattice)#
+myLay <- grid.layout(1, 3, widths=unit(3, "inches"),#
+ heights=unit(2, "inches"))#
+#
+pushViewport(viewport(layout=myLay))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+pushViewport(viewport(layout.pos.col=1, name="v11"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p05, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=2, name="v12"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p065, newpage=F)#
+ upViewport()#
+pushViewport(viewport(layout.pos.col=3, name="v13"))#
+ grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
+ print(p12, newpage=F)#
+ upViewport()
+#
+mdf <- with(xy, #
+data.frame(x=x, y=y, z=c(z05, z065, z12), wavelength=factor(rep(c(0.5, 0.65, 1.2), each=length(x)))))
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy",scales="free", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy",scales=list(z="free"), xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)
+levelplot
+?levelplot
+mdf <- with(xy, #
+data.frame(x=x, y=y, z=c(z05/max(z05), z065/max(z065), z12/max(z12)), wavelength=factor(rep(c(0.5, 0.65, 1.2), each=length(x)))))#
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ panel.points(0, 0, cex=10)#
+})
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ panel.points(0, 0, cex=10, fill="black")#
+})
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ panel.points(0, 0, cex=10, pch=21, fill="black")#
+})
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ panel.points(0, 0, cex=12.5, pch=21, fill="black")#
+})
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = 0.2, gp=gpar(col=NA, fill=grey(0.5)))#
+})
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = 0.25, gp=gpar(col=NA, fill=grey(0.5)))#
+})
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "npc"), gp=gpar(col=NA, fill=grey(0.5)))#
+})
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col=NA, fill=grey(0.5)))#
+})
+#
+#
+library(mie)#
+library(grid)#
+data(AuJChristy)#
+plot.material()$fine.sample->material#
+#
+radius <- 0.2#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))#
+#
+#
+#
+#
+lambda <- 1.2#
+lambda <- 0.5#
+k0 <- 2*pi/lambda#
+#
+eps <- permittivity(wav=lambda)#
+n <- Re(sqrt(eps))#
+k <- Im(sqrt(eps))#
+#
+ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+#
+scale <- 2#
+N <- 200#
+xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+#
+cart2polar <- function(xy){#
+r <- sqrt(xy[, 1]^2+xy[, 2]^2) #
+angle <- atan(xy[, 2]/xy[, 1])#
+data.frame(r=r, angle=angle)#
+}#
+#
+rphi <- cart2polar(xy)#
+remove <- which(rphi$r <= radius)#
+#
+my.field12 <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=2*pi/1.2, phi=0))#
+my.field065 <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=2*pi/0.65, phi=0))#
+my.field05 <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=2*pi/0.5, phi=0))#
+#
+xy$z12 <- 0#
+xy$z05 <- 0#
+xy$z065 <- 0#
+xy$z12[-remove] <- c(apply(my.field12, 2, function(x) Re(x%*%Conj(x))))#
+xy$z05[-remove] <- c(apply(my.field05, 2, function(x) Re(x%*%Conj(x))))#
+xy$z065[-remove] <- c(apply(my.field065, 2, function(x) Re(x%*%Conj(x))))#
+#
+#
+#
+p12 <- levelplot(z12~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)#
+#
+p065 <- levelplot(z065~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)#
+#
+p05 <- levelplot(z05~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE)#
+#
+mdf <- with(xy, #
+data.frame(x=x, y=y, z=c(z05/max(z05), z065/max(z065), z12/max(z12)), wavelength=factor(rep(c(0.5, 0.65, 1.2), each=length(x)))))#
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col=NA, fill=grey(0.5)))#
+})
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), between=list(x=0.2), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), between=list(x=0.2), strip=strip.custom(var=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})
+?strip.default
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), between=list(x=0.2), strip=strip.custom(var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ strip.names = c(T, TRUE),#
+ ylab=expression("distance / "*mu*m), between=list(x=0.2), strip=strip.custom(var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), between=list(x=0.2), strip=strip.custom(strip.names = c(T, TRUE),var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})
+#
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), between=list(x=0.2), strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})
+#
+pdf("fieldPlot.pdf")#
+print(field.plt)#
+dev.off()
+#
+field.plot <- levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+ ylab=expression("distance / "*mu*m), between=list(x=0.2), strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+pdf("fieldPlot.pdf")#
+print(field.plot)#
+dev.off()
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2), layout=c(3, 1), #
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})
+?pdf
+#
+quartz(file="fieldPlot.pdf", height=3)#
+print(field.plot)#
+dev.off()
+#
+quartz(type="pdf", file="fieldPlot.pdf", height=3)#
+print(field.plot)#
+dev.off()
+quartz(type="native", file="fieldPlot.pdf", height=3)#
+print(field.plot)
+levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2), layout=c(1, 3), #
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})
+quartz(type="native", file="fieldPlot.pdf", height=5)#
+print(field.plot)
+#
+field.plot <- levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2), layout=c(3, 1), #
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})
+#
+quartz(type="native", file="fieldPlot.pdf", height=5)#
+print(field.plot)
+#
+quartz(type="png", file="MieNearFieldPlot.png", height=5)#
+#
+print(field.plot)#
+dev.off()
+#
+quartz(type="png", file="MieNearFieldPlot.png", height=5, width=7)#
+#
+print(field.plot)#
+dev.off()
+#
+#
+quartz(type="png", file="MieNearFieldPlot.png", height=5, width=7, bg="white", dpi=200)#
+#
+print(field.plot)#
+dev.off()
+#
+quartz(type="png", file="MieNearFieldPlot.png", height=4, width=7, bg="white", dpi=200)#
+#
+print(field.plot)#
+dev.off()
+#
+radius <- 0.1#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+#
+radius <- 0.2#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=T)$fine.sample->material
+#
+data(AgPalik)#
+plot.material("AgPalik", plot=T, range=c(0.3, 0.9))$fine.sample->material
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=T, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.2#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=T, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.05#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.5), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.05#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.0), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.05#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.5), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material()$fine.sample->material#
+#
+radius <- 0.05#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.5), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy")$fine.sample->material#
+#
+radius <- 0.05#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.5), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F)$fine.sample->material#
+#
+radius <- 0.05#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.5), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.05#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.5), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.05#
+#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.5), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.04#
+#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.1), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.04#
+#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.6), t="l", ylab=~sigma))
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.04#
+#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.5), t="l", ylab=~sigma))
+#
+with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
+ matplot(lambda, mdata, t="l", ylab=~sigma)#
+})
+head(mdata)
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.04#
+#
+#
+#
+#
+mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
+ matplot(lambda, mdata, t="l", ylab=~sigma)#
+})#
+#
+head(mdata)
+mdata
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.04#
+#
+#
+#
+#
+mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
+ matplot(lambda, mdata, t="l", ylab=~sigma);#
+ mdata#
+})#
+#
+head(mdata)
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.04#
+#
+#
+#
+#
+mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
+ matplot(lambda, mdata, t="l", ylab=~sigma);#
+ cbind(lambda, mdata)#
+})#
+#
+head(mdata)
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.04#
+#
+#
+#
+#
+mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
+ matplot(lambda, mdata, t="l", ylab=~sigma);#
+ cbind(lambda, mdata)#
+})#
+#
+names(mdata) <- c("wavelength", "cext", "cabs", "csat")#
+head(mdata)
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.04#
+#
+#
+#
+#
+mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
+ matplot(lambda, mdata, t="l", ylab=~sigma);#
+ as.data.frame(cbind(lambda, mdata))#
+})#
+#
+names(mdata) <- c("wavelength", "cext", "cabs", "csat")#
+head(mdata)
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+plot.material("AuJChristy", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+radius <- 0.04#
+#
+#
+#
+#
+mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
+ matplot(lambda, mdata, t="l", ylab=~sigma);#
+ as.data.frame(cbind(lambda, mdata))#
+})#
+#
+names(mdata) <- c("wavelength", "cext", "cabs", "csat")#
+head(mdata)#
+#
+export2excel(mdata, file="chris.txt")
+system("head chris.txt")
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+#
+radius <- 0.04#
+#
+#
+#
+#
+mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
+ matplot(lambda, mdata, t="l", ylab=~sigma);#
+ as.data.frame(cbind(lambda, mdata))#
+})#
+#
+names(mdata) <- c("wavelength", "cext", "cabs", "csat")#
+head(mdata)#
+#
+export2excel(mdata, file="chris.txt")
+library(constants)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+#
+radius <- 0.04#
+#
+mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
+ matplot(lambda, mdata, t="l", ylab=~sigma);#
+ as.data.frame(cbind(lambda, mdata))#
+})#
+#
+names(mdata) <- c("wavelength", "cext", "cabs", "csat")#
+head(mdata)#
+#
+export2excel(mdata, file="chris.txt")
+setwd('/Users/baptiste/Documents/R/testFortran')
+library(mie)#
+data(AgPalik)#
+plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+#
+#
+radius <- 0.04#
+#
+mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
+ matplot(lambda, mdata, t="l", ylab=~sigma);#
+ as.data.frame(cbind(lambda, mdata))#
+})#
+#
+names(mdata) <- c("wavelength", "cext", "cabs", "csat")#
+head(mdata)#
+#
+export2excel(mdata, file="chris.txt")
+radius <- 0.05#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+#
+library(mie)#
+library(grid)#
+data(AuJChristy)#
+plot.material()$fine.sample->material#
+#
+radius <- 0.2#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))#
+radius <- 0.05#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+#
+detach(package:constants)
+plot.material
+#
+plot.material()$fine.sample->material#
+#
+radius <- 0.2#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))#
+radius <- 0.05#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+radius <- 0.1#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+#
+oneField <- function(lambda=0.5, radius=0.1, N=10, scale=2){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+ #
+ return(xy)#
+}#
+#
+oneField()
+rep(c(0.5, 0.65, 1.2), 9)
+rep(c(0.05, 0.1, 0.2), each=3)
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)
+allFields
+oneField <- function(lambda=0.5, radius=0.1, N=10, scale=2){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+ #
+ return(cbind(xy, radius=radius, lambda=lambda)#
+}#
+#
+oneField()
+#
+oneField <- function(lambda=0.5, radius=0.1, N=10, scale=2){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+ #
+ return(cbind(xy, radius=radius, lambda=lambda))#
+}#
+#
+oneField()
+oneField <- function(lambda=0.5, radius=0.1, N=5, scale=2){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+ #
+ as.data.frame(cbind(xy, radius=radius, lambda=lambda))#
+}#
+#
+oneField()#
+#
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+mdf <- do.call(rbind, allFields)
+mdf
+#
+field.plot <- levelplot(z~x*y|wavelength*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2), layout=c(3, 1), #
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2), layout=c(3, 1), #
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+#
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2), layout=c(3, 2), #
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2), layout=c(3, 2), #
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+mdf <- do.call(rbind, allFields)#
+#
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2), layout=c(3, 2), #
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+#
+field.plot <- levelplot(z~x*y|lambda, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2), layout=c(3, 2), #
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+head(mdf)
+str(mdf)
+mdf$lambda
+mdf$radius
+#
+#
+field.plot <- levelplot(z~x*y|lambda, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+#
+oneField <- function(lambda=0.5, radius=0.1, N=5, scale=2){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+ xy$z <- xy$z/max(xy$z)#
+ as.data.frame(cbind(xy, radius=radius, lambda=lambda))#
+}#
+#
+oneField()#
+#
+#
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+mdf <- do.call(rbind, allFields)#
+#
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+#
+#
+oneField <- function(lambda=0.5, radius=0.1, N=10, scale=2){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+ xy$z <- xy$z/max(xy$z)#
+ as.data.frame(cbind(xy, radius=radius, lambda=lambda))#
+}#
+#
+oneField()#
+#
+#
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+mdf <- do.call(rbind, allFields)#
+#
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+#
+#
+oneField <- function(lambda=0.5, radius=0.1, N=10, region=1.0){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+ xy$z <- xy$z/max(xy$z)#
+ as.data.frame(cbind(xy, radius=radius, lambda=lambda))#
+}#
+#
+oneField()#
+#
+#
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+mdf <- do.call(rbind, allFields)#
+#
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+#
+oneField <- function(lambda=0.5, radius=0.1, N=10, region=1.0){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+ xy$z <- xy$z/max(xy$z)#
+ as.data.frame(cbind(xy, radius=radius, lambda=lambda))#
+}#
+#
+oneField()
+#
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+mdf <- do.call(rbind, allFields)
+#
+oneField(lambda=1)
+#
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)
+remove
+xy$z[-remove]
+xy
+dim(xy)
+lambda=0.5
+radius=0.1
+ N=10
+region=1.0
+#
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)
+#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))
+ my.field
+ xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))
+ xy$z <- xy$z/max(xy$z)
+#
+ as.data.frame(cbind(xy, radius=radius, lambda=lambda))
+#
+oneField(radius=0.2)
+#
+oneField <- function(lambda=0.5, radius=0.1, N=10, region=1.0){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+ xy$z <- xy$z/max(xy$z)#
+ as.data.frame(cbind(xy, radius=radius, lambda=lambda))#
+}#
+#
+oneField(radius=0.5)
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+mdf <- do.call(rbind, allFields)
+#
+oneField(radius=1)
+#
+oneField(radius=0.5)
+#
+oneField(radius=0.05)
+oneField <- function(lambda=0.5, radius=0.1, N=10, region=1.0){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ try(xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x)))) )#
+ xy$z <- xy$z/max(xy$z)#
+ as.data.frame(cbind(xy, radius=radius, lambda=lambda))#
+}#
+#
+oneField(radius=0.5)
+#
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+#
+mdf <- do.call(rbind, allFields)
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+#
+oneField <- function(lambda=0.5, radius=0.1, N=5, region=1.0){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ try(xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x)))) )#
+ xy$z <- xy$z/max(xy$z)#
+ as.data.frame(cbind(xy, radius=radius, lambda=lambda))#
+}#
+#
+oneField(radius=0.5)#
+#
+#
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+#
+mdf <- do.call(rbind, allFields)#
+#
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+oneField <- function(lambda=0.5, radius=0.1, N=5, region=1.0){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r <= radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ try(xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x)))) )#
+ xy$z <- xy$z/max(xy$z)#
+ as.data.frame(cbind(xy, radius=factor(radius), lambda=factor(lambda)))#
+}#
+#
+oneField(radius=0.5)#
+#
+#
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+#
+mdf <- do.call(rbind, allFields)#
+#
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))#
+#
+#
+cart2polar <- function(xy){#
+r <- sqrt(xy[, 1]^2+xy[, 2]^2) #
+angle <- atan(xy[, 2]/xy[, 1])#
+data.frame(r=r, angle=angle)#
+}#
+#
+oneField <- function(lambda=0.5, radius=0.1, N=5, region=1.0){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r < radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ try(xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x)))) )#
+ xy$z <- xy$z/max(xy$z)#
+ as.data.frame(cbind(xy, radius=factor(radius), lambda=factor(lambda)))#
+}#
+#
+oneField(radius=0.5)
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+#
+mdf <- do.call(rbind, allFields)#
+#
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+#
+oneField <- function(lambda=0.5, radius=0.1, N=20, region=1.0){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r < radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ try(xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x)))) )#
+ xy$z <- xy$z/max(xy$z)#
+ as.data.frame(cbind(xy, radius=factor(radius), lambda=factor(lambda)))#
+}#
+#
+#
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+#
+mdf <- do.call(rbind, allFields)#
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ grid.circle(name="my.circle", r = unit(0.05, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+field.plot
+library(latticeExtra)#
+useOuterStrips(field.plot)
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(...){#
+ #
+ panel.levelplot(...)#
+ #
+ my.rad <- as.numeric(unique(xy$radius[subscripts]))#
+ grid.circle(name="my.circle", r = unit(my.rad, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+library(latticeExtra)#
+useOuterStrips(field.plot)
+#
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(..., subscripts){#
+ #
+ panel.levelplot(...)#
+ #
+ my.rad <- as.numeric(unique(xy$radius[subscripts]))#
+ grid.circle(name="my.circle", r = unit(my.rad, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+library(latticeExtra)#
+useOuterStrips(field.plot)
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(..., subscripts){#
+ #
+ panel.levelplot(..., subscripts=subscripts)#
+ #
+ my.rad <- as.numeric(unique(xy$radius[subscripts]))#
+ grid.circle(name="my.circle", r = unit(my.rad, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+library(latticeExtra)#
+useOuterStrips(field.plot)
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(..., subscripts){#
+ #
+ panel.levelplot(..., subscripts=subscripts)#
+ #
+ my.rad <- as.numeric(unique(xy$radius[subscripts]))#
+ print(my.rad)#
+ grid.circle(name="my.circle", r = unit(my.rad, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+library(latticeExtra)#
+useOuterStrips(field.plot)
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(..., subscripts){#
+ #
+ panel.levelplot(..., subscripts=subscripts)#
+ #
+ my.rad <- as.numeric(unique(mdf$radius[subscripts]))#
+ print(my.rad)#
+ grid.circle(name="my.circle", r = unit(my.rad, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+library(latticeExtra)#
+useOuterStrips(field.plot)
+fact2numeric
+fact2num
+factor2num
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(..., subscripts){#
+ #
+ panel.levelplot(..., subscripts=subscripts)#
+ #
+ my.rad <- factor2num(unique(mdf$radius[subscripts]))#
+ print(my.rad)#
+ grid.circle(name="my.circle", r = unit(my.rad, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+library(latticeExtra)#
+useOuterStrips(field.plot)
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(y=0.2, x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(..., subscripts){#
+ #
+ panel.levelplot(..., subscripts=subscripts)#
+ #
+ my.rad <- factor2num(unique(mdf$radius[subscripts]))#
+ print(my.rad)#
+ grid.circle(name="my.circle", r = unit(my.rad, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+library(latticeExtra)#
+useOuterStrips(field.plot)
+#
+oneField <- function(lambda=0.5, radius=0.1, N=150, region=1.0){#
+ #
+ k0 <- 2*pi/lambda#
+ #
+ eps <- permittivity(wav=lambda)#
+ n <- Re(sqrt(eps))#
+ k <- Im(sqrt(eps))#
+ #
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+ #
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))#
+ xy$z <- 0#
+ rphi <- cart2polar(xy)#
+ remove <- which(rphi$r < radius)#
+ #
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))#
+ try(xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x)))) )#
+ xy$z <- xy$z/max(xy$z)#
+ as.data.frame(cbind(xy, radius=factor(radius), lambda=factor(lambda)))#
+}#
+#
+#
+#
+#
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)#
+#
+mdf <- do.call(rbind, allFields)#
+#
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy", #
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m), #
+ between=list(y=0.2, x=0.2),#
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)), #
+ colorkey = T, region = TRUE, panel=function(..., subscripts){#
+ #
+ panel.levelplot(..., subscripts=subscripts)#
+ #
+ my.rad <- factor2num(unique(mdf$radius[subscripts]))#
+ print(my.rad)#
+ grid.circle(name="my.circle", r = unit(my.rad, "native"), gp=gpar(col="black", fill=grey(0.5)))#
+})#
+#
+library(latticeExtra)#
+useOuterStrips(field.plot)
+#
+save("mieNearField.rda", field.plot, mdf)
+#
+save(file="mieNearField.rda", field.plot, mdf)
+setwd('/Users/baptiste/Documents/R/photonics/pkg/mie/inst')
+#
+quartz(type="png", file="MieNearFieldPlot2.png", height=7, width=7, bg="white", dpi=200)#
+#
+print(useOuterStrips(field.plot))#
+dev.off()
Added: pkg/mie/inst/MieNearFieldPlot.png
===================================================================
(Binary files differ)
Property changes on: pkg/mie/inst/MieNearFieldPlot.png
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/mie/inst/MieNearFieldPlot2.png
===================================================================
(Binary files differ)
Property changes on: pkg/mie/inst/MieNearFieldPlot2.png
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/mie/inst/Rplot.pdf
===================================================================
(Binary files differ)
Property changes on: pkg/mie/inst/Rplot.pdf
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/mie/inst/fieldPlot.pdf
===================================================================
(Binary files differ)
Property changes on: pkg/mie/inst/fieldPlot.pdf
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/mie/inst/mieNearField.rda
===================================================================
(Binary files differ)
Property changes on: pkg/mie/inst/mieNearField.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: pkg/mie/inst/test.r
===================================================================
--- pkg/mie/inst/test.r (rev 0)
+++ pkg/mie/inst/test.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,84 @@
+# library(baptMisc)
+# clr()
+library(mie)
+library(grid)
+data(AuJChristy)
+plot.material()$fine.sample->material
+
+radius <- 0.2
+# 1.2um is dipolar
+# 0.5um is octupolar
+# 0.65 is quadrupolar
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+
+# 1.2um is dipolar
+# 0.5um is octupolar
+# 0.65 is quadrupolar
+
+lambda <- 1.2
+lambda <- 0.5
+k0 <- 2*pi/lambda
+
+eps <- permittivity(wav=lambda)
+n <- Re(sqrt(eps))
+k <- Im(sqrt(eps))
+
+ab <- Mie.ab(m=n+1i*k, x=k0 * radius)
+cd <- Mie.cd(m=n+1i*k, x=k0 * radius)
+
+scale <- 2
+N <- 200
+xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))
+
+cart2polar <- function(xy){
+r <- sqrt(xy[, 1]^2+xy[, 2]^2)
+angle <- atan(xy[, 2]/xy[, 1])
+data.frame(r=r, angle=angle)
+}
+
+rphi <- cart2polar(xy)
+remove <- which(rphi$r <= radius)
+
+my.field12 <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=2*pi/1.2, phi=0))
+my.field065 <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=2*pi/0.65, phi=0))
+my.field05 <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=2*pi/0.5, phi=0))
+# int.field <- mapply(Eint, r=rphi$r[remove], theta=rphi$angle[remove], MoreArgs=list(cd=cd,k=n*k0, phi=0)) # dubious, is it n+1i*k?
+xy$z12 <- 0
+xy$z05 <- 0
+xy$z065 <- 0
+xy$z12[-remove] <- c(apply(my.field12, 2, function(x) Re(x%*%Conj(x))))
+xy$z05[-remove] <- c(apply(my.field05, 2, function(x) Re(x%*%Conj(x))))
+xy$z065[-remove] <- c(apply(my.field065, 2, function(x) Re(x%*%Conj(x))))
+# xy$z[remove] <- c(apply(int.field, 2, function(x) Re(x%*%Conj(x))))
+# xy$z[remove] <- NA#c(apply(int.field, 2, function(x) Re(x%*%Conj(x))))
+
+p12 <- levelplot(z12~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),
+ ylab=expression("distance / "*mu*m),
+ colorkey = T, region = TRUE)
+
+p065 <- levelplot(z065~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),
+ ylab=expression("distance / "*mu*m),
+ colorkey = T, region = TRUE)
+
+p05 <- levelplot(z05~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),
+ ylab=expression("distance / "*mu*m),
+ colorkey = T, region = TRUE)
+
+mdf <- with(xy,
+data.frame(x=x, y=y, z=c(z05/max(z05), z065/max(z065), z12/max(z12)), wavelength=factor(rep(c(0.5, 0.65, 1.2), each=length(x)))))
+
+field.plot <- levelplot(z~x*y|wavelength, mdf, cuts = 100,aspect="xy",
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m),
+ between=list(x=0.2), layout=c(3, 1),
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)),
+ colorkey = T, region = TRUE, panel=function(...){
+
+ panel.levelplot(...)
+ # panel.points(0, 0, cex=12.5, pch=21, fill="black")
+ grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))
+})
+
+quartz(type="png", file="MieNearFieldPlot2.png", height=4, width=7, bg="white", dpi=200)
+# quartz(type="native", file="fieldPlot.png", height=5)
+print(field.plot)
+dev.off()
Added: pkg/mie/inst/test2.r
===================================================================
--- pkg/mie/inst/test2.r (rev 0)
+++ pkg/mie/inst/test2.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,81 @@
+# library(baptMisc)
+# clr()
+detach(package:constants)
+library(mie)
+library(grid)
+data(AuJChristy)
+plot.material()$fine.sample->material
+
+radius <- 0.2
+# 1.2um is dipolar
+# 0.5um is octupolar
+# 0.65 is quadrupolar
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+
+radius <- 0.1
+# 0.5um is octupolar
+# 0.65 is quadrupolar
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+
+
+cart2polar <- function(xy){
+r <- sqrt(xy[, 1]^2+xy[, 2]^2)
+angle <- atan(xy[, 2]/xy[, 1])
+data.frame(r=r, angle=angle)
+}
+
+oneField <- function(lambda=0.5, radius=0.1, N=150, region=1.0){
+
+ k0 <- 2*pi/lambda
+
+ eps <- permittivity(wav=lambda)
+ n <- Re(sqrt(eps))
+ k <- Im(sqrt(eps))
+
+ ab <- Mie.ab(m=n+1i*k, x=k0 * radius)
+ cd <- Mie.cd(m=n+1i*k, x=k0 * radius)
+
+ xy <- expand.grid(x=seq(-region/2, region/2, length=N), y=seq(-region/2,region/2, length=N))
+ xy$z <- 0
+ rphi <- cart2polar(xy)
+ remove <- which(rphi$r < radius)
+
+ my.field <- mapply(Escat, r=rphi$r[-remove], theta=rphi$angle[-remove], MoreArgs=list(ab=ab,k=k0, phi=0))
+ try(xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x)))) )
+ xy$z <- xy$z/max(xy$z)
+ as.data.frame(cbind(xy, radius=factor(radius), lambda=factor(lambda)))
+}
+
+# oneField(radius=0.5)
+
+
+allFields <- mapply(oneField, lambda=rep(c(0.5, 0.65, 1.2), 3), radius=rep(c(0.05, 0.1, 0.2), each=3), SIMPLIFY=F)
+
+mdf <- do.call(rbind, allFields)
+
+field.plot <- levelplot(z~x*y|lambda*radius, mdf, cuts = 100,aspect="xy",
+ xlab=expression("distance / "*mu*m), ylab=expression("distance / "*mu*m),
+ between=list(y=0.2, x=0.2),
+ strip=strip.custom(strip.names = c(T, TRUE),sep="=", var.name=expression(lambda)),
+ colorkey = T, region = TRUE, panel=function(..., subscripts){
+
+ panel.levelplot(..., subscripts=subscripts)
+ # panel.points(0, 0, cex=12.5, pch=21, fill="black")
+ my.rad <- factor2num(unique(mdf$radius[subscripts]))
+ print(my.rad)
+ grid.circle(name="my.circle", r = unit(my.rad, "native"), gp=gpar(col="black", fill=grey(0.5)))
+})
+
+library(latticeExtra)
+
+#
+
+# save(file="mieNearField.rda", field.plot, mdf)
+#
+
+
+quartz(type="png", file="MieNearFieldPlot2.png", height=7, width=7, bg="white", dpi=200)
+# quartz(type="native", file="fieldPlot.png", height=5)
+print(useOuterStrips(field.plot))
+dev.off()
+
Added: pkg/mie/man/Mie.Rd
===================================================================
--- pkg/mie/man/Mie.Rd (rev 0)
+++ pkg/mie/man/Mie.Rd 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,79 @@
+\name{Mie}
+\alias{Mie}
+\alias{Mie.ab}
+\alias{Mie.coated.ab}
+\alias{Mie.coated.spectra}
+\alias{Mie.spectra}
+
+\title{Mie scattering}
+\description{
+Computation of Mie Efficiencies / cross-sections
+
+}
+\usage{
+Mie(m, x)
+Mie.ab(m = 1, x = 1)
+Mie.coated.ab(m1 = 1, m2 = 1, x = 1, y = 1)
+Mie.ab(m = 1, x = 1)
+Mie.spectra(lambda, epsilon, radius = 0.04, medium = 1.333, efficiency = FALSE)
+Mie.coated.spectra(lambda, epsilon, inner.radius = 0.04,
+ outer.radius = 0.045, shell = 1.5, medium = 1.333)
+}
+
+\arguments{
+ \item{m}{complex refractive index relative to the ambient medium }
+ \item{x}{size parameter in the ambient medium}
+ \item{lambda}{vector, wavelengths [um]}
+ \item{epsilon}{(complex) vector, permittivity values}
+ \item{(inner|outer.)radius}{scalar, [um] }
+ \item{medium}{scalar, incident medium refractive index}
+ \item{efficiency}{boolean, TRUE for efficiencies, FALSE for cross sections }
+}
+\details{
+Mie.spectra and Mie.coated.spectra are the high level functions that should be used for most applications. See references.
+}
+\value{
+data.frame,
+\item{qext }{extinction}
+\item{qsca }{scattering}
+\item{qabs }{absorption}
+}
+\references{Bohren and Huffman (1983)
+BEWI:TDD122, p. 103,119-122,477
+C. Matzler, May 2002, revised July 2002.}
+\author{baptiste Auguie}
+\note{
+efficiency = TRUE is not defined for coated spheres, its practical definition not being very clear.
+}
+\seealso{\code{\link{Mie.spectra}}, \code{\link{Mie.ab}}}
+\examples{
+library(scatteringR)
+
+library(constants)
+clr()
+graphics.off()
+require(RColorBrewer)
+palette(brewer.pal(8,"Set1"))
+
+data(goldData)
+attach(goldData)
+subset<- sort(sample(length(goldData[[1]]),100))
+goldData <-as.data.frame(goldData)[subset,]
+attach(goldData)
+
+XsectionsBare <- Mie.spectra(lambda = lambda, radius =0.05,
+ epsilon = epsilon,
+ medium= 1.333)
+
+XsectionsLayer <- Mie.coated.spectra(lambda = lambda,
+ epsilon = epsilon,
+ inner.radius = 0.05,
+ outer.radius = 0.055,
+ shell =1.5 ,
+ medium= 1.333)
+
+Xsections <- cbind(XsectionsBare, XsectionsLayer)
+matplot(lambda, Xsections,t="l",lwd=2,lty=rep(1:2,each=3),col=rep(1:3,2),xlab = expression(paste("wavelength [nm]")), ylab = expression(paste(sigma, " [",mu,m^2, "]")))
+legend("topright",fill=c(1:3),c("extinction","scattering","absorption"))
+}
+
Added: pkg/mie/testMat.r
===================================================================
--- pkg/mie/testMat.r (rev 0)
+++ pkg/mie/testMat.r 2008-11-15 12:36:10 UTC (rev 328)
@@ -0,0 +1,3 @@
+
+testmat1 <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16), nrow=4)
+testmat2 <- matrix(c(1,2,3,5,5,6,7,8,9,10,11,12,13,14,15,16), nrow=4)
\ No newline at end of file
More information about the Photonics-commits
mailing list