[Photonics-commits] r337 - in pkg: SPPr/man Sparameters/man arrays/man baptMisc/man ebeamR/man mie/inst scatteringR/man spectroR/man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Nov 18 23:47:26 CET 2008


Author: baptiste
Date: 2008-11-18 23:47:26 +0100 (Tue, 18 Nov 2008)
New Revision: 337

Removed:
   pkg/baptMisc/man/cutIntervals.Rd
Modified:
   pkg/SPPr/man/SPPr-package.Rd
   pkg/SPPr/man/objF.Rd
   pkg/Sparameters/man/effectiveMedium.Rd
   pkg/arrays/man/plotSubSample.T.Rd
   pkg/ebeamR/man/ebeamR-package.Rd
   pkg/ebeamR/man/writePositionList.Rd
   pkg/mie/inst/.Rhistory
   pkg/scatteringR/man/Mie.Rd
   pkg/scatteringR/man/alpha_s.Rd
   pkg/scatteringR/man/array.factor.Rd
   pkg/scatteringR/man/runAdda.spectrum.Rd
   pkg/spectroR/man/assignFiles.Rd
   pkg/spectroR/man/exampleDF.Rd
   pkg/spectroR/man/extractSpectra.Rd
   pkg/spectroR/man/extractTspectra.Rd
   pkg/spectroR/man/locatePeaks.Rd
   pkg/spectroR/man/spectroR-package.Rd
Log:
man pages cleanup

Modified: pkg/SPPr/man/SPPr-package.Rd
===================================================================
--- pkg/SPPr/man/SPPr-package.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/SPPr/man/SPPr-package.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -31,8 +31,7 @@
 
 }
 \examples{
-clr()
-graphics.off()
+
 library(SPPr)
 
 data(scanSPP)

Modified: pkg/SPPr/man/objF.Rd
===================================================================
--- pkg/SPPr/man/objF.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/SPPr/man/objF.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -24,8 +24,7 @@
 }
 \seealso{}
 \examples{
-clr()
-graphics.off()
+	
 parameters <- list(layer2fit=1,epsr=c(-12),epsi=c(1.2),thick=c(40),
 cel=3e8,lambda=632.8e-9,nPrism=1.46,anglePrism=pi/3,nWater=1,
 zeroShift=-0.4,normalisation=23,pIni=c(-12,1.2,40,1))

Modified: pkg/Sparameters/man/effectiveMedium.Rd
===================================================================
--- pkg/Sparameters/man/effectiveMedium.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/Sparameters/man/effectiveMedium.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -37,7 +37,6 @@
 \seealso{}
 \examples{
 library(Sparameters)
-clr()
 
 data(exampleData)
 attach(exampleData)

Modified: pkg/arrays/man/plotSubSample.T.Rd
===================================================================
--- pkg/arrays/man/plotSubSample.T.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/arrays/man/plotSubSample.T.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -30,22 +30,22 @@
 }
 \seealso{\code{\link{plotSubSample.E}} for scaled extinction}
 \examples{
-library(arrays)
-library(constants)
-library(RColorBrewer)
-clr()
-#?arrays
-data(sample1)
-data(sample2)
-my.cols<-palette(c(brewer.pal(9,"Set1"),"black")[-c(6,8)])
-names(sample1)
-names(sample2)
-levels(sample2$disorder)
-cond.1 <- quote(size \%in\% c(1,3) & (polarisation == "long") & dose == "2" & pitch == 540) 
-cond.2 <- quote(disorder == 0 & (polarisation == "long") & dose == "1" & pitch == 540) 
-par() -> b.quiet # default
-par(bty="n",mfrow=c(2,1),mar=par("mar")+0.4)
-plotSubSample.E(sample1,cond.1,col=my.cols,xlim=c(480,920),lwd=2 )
-plotSubSample.E(sample2,cond.2,xlim=c(480,920),lwd=2 )
-par() -> b.quiet
+# library(arrays)
+# library(constants)
+# library(RColorBrewer)
+# 
+# #?arrays
+# data(sample1)
+# data(sample2)
+# my.cols<-palette(c(brewer.pal(9,"Set1"),"black")[-c(6,8)])
+# names(sample1)
+# names(sample2)
+# levels(sample2$disorder)
+# cond.1 <- quote(size \%in\% c(1,3) & (polarisation == "long") & dose == "2" & pitch == 540) 
+# cond.2 <- quote(disorder == 0 & (polarisation == "long") & dose == "1" & pitch == 540) 
+# par() -> b.quiet # default
+# par(bty="n",mfrow=c(2,1),mar=par("mar")+0.4)
+# plotSubSample.E(sample1,cond.1,col=my.cols,xlim=c(480,920),lwd=2 )
+# plotSubSample.E(sample2,cond.2,xlim=c(480,920),lwd=2 )
+# par() -> b.quiet
 }

Deleted: pkg/baptMisc/man/cutIntervals.Rd
===================================================================
--- pkg/baptMisc/man/cutIntervals.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/baptMisc/man/cutIntervals.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -1,29 +0,0 @@
-\name{cutIntervals}
-\alias{cutIntervals}
-%- Also NEED an '\alias' for EACH other topic documented here.
-\title{ splits the data into intervals}
-\description{
-  alternative to cut that returns numeric values
-}
-\usage{
-cutIntervals(x, ...)
-}
-\arguments{
-  \item{x}{vector}
-  \item{\dots}{arguments passed to cut(x, ...)}
-}
-\details{
-	calls cut(x, lebels=NULL, ...) and extracts the factor levels in numeric form
-}
-\value{matrix: intervals
-}
-\references{R-help mailing list,  cut}
-\author{baptiste Auguie}
-\note{ 
-	\section{Warning}{labels argument is not allowed,  use cut() directly in this case}
-}
-\seealso{\code{\link{cut}}, \code{\link{hist}}}
-\examples{
-	# not run
-# cutIntervals(1:5, 3)
-}

Modified: pkg/ebeamR/man/ebeamR-package.Rd
===================================================================
--- pkg/ebeamR/man/ebeamR-package.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/ebeamR/man/ebeamR-package.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -31,8 +31,6 @@
 
 }
 \examples{
-#require(constants)
-#clr() # clear workspace
 require(ebeamR)
 
 data(parameters)

Modified: pkg/ebeamR/man/writePositionList.Rd
===================================================================
--- pkg/ebeamR/man/writePositionList.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/ebeamR/man/writePositionList.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -27,8 +27,6 @@
 }
 \seealso{\code{\link{printParticle}} }
 \examples{
-#require(constants)
-#clr() # clear workspace
 require(ebeamR)
 
 data(parameters)

Modified: pkg/mie/inst/.Rhistory
===================================================================
--- pkg/mie/inst/.Rhistory	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/mie/inst/.Rhistory	2008-11-18 22:47:26 UTC (rev 337)
@@ -1,3 +1,7 @@
+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))
@@ -32,6 +36,7 @@
 help()
 sessionInfo
 sessionInfo()
+library(grid)
 #
 #
 #
@@ -76,7 +81,6 @@
 ?spectra
 library(constants)
 spectra()->test
-test
 library(spectra)
 spectra2excel(test,"test.txt")
 getWavelength(spectra())
@@ -239,6 +243,7 @@
 12753, 12996, 13057, 13149), class = "Date"), class = "zoo"))#
 #
 plot(d)
+d
 library(ifs)
 ?ifs
 #
@@ -386,7 +391,6 @@
 #
 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)#
 #
@@ -425,1856 +429,1685 @@
 grid.rect(gp = gpar(col = rgb(43/255, 140/255, 190/255)))#
 print(p2, newpage=F)#
 upViewport()
+?grid.layout
 #
-d = data.frame(a=1, b=2)#
-a = c("a", "b")#
-z = a
-d
-a
-z
 #
-subset(d, select=z)
-?subset
+library(mie)#
+library(grid)#
+data(AuJChristy)#
+plot.material()$fine.sample->material
 #
-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)#
+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))#
+eps1 <- 1.5^2#
+eps2 <- material$epsilon
+eps2
 #
 #
+kz <- matrix(sqrt(eps1*k0^2 - kx^2) + 0i, nrow = N.k, ncol=N.t, byrow=F)  #
+epsilon <- cbind(rep(eps1, length=length(eps2)), eps2, eps1)
 #
 #
-lambda <- 1.2#
+eps1 <- 1.5^2#
+eps2 <- material$epsilon#
+lambda <- material$lambda#
+k0 <- 2*pi/lambda#
+N.t <- 3#
 #
+kz <- matrix(sqrt(eps1*k0^2 - kx^2) + 0i, nrow = N.k, ncol=N.t, byrow=F)  #
+epsilon <- cbind(rep(eps1, length=length(eps2)), eps2, eps1)
+#
+kx <- k0*sin(pi/4)#
+kz <- matrix(sqrt(eps1*k0^2 - kx^2) + 0i, nrow = N.k, ncol=N.t, byrow=F)  #
+epsilon <- cbind(rep(eps1, length=length(eps2)), eps2, eps1)
+#
 k0 <- 2*pi/lambda#
+N.t <- 3#
+N.k <-length(k0)#
+kx <- k0*sin(pi/4)#
 #
-eps <- permittivity(wav=lambda)#
-n <- Re(sqrt(eps))#
-k <- Im(sqrt(eps))#
+kz <- matrix(sqrt(eps1*k0^2 - kx^2) + 0i, nrow = N.k, ncol=N.t, byrow=F)  #
+epsilon <- cbind(rep(eps1, N.k), eps2, eps1)
+epsilon
+epsilon[, 2]
+kz/epsilon
 #
-ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
-cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+eps1 <- 1.5^2#
+eps3 <- 1.5^2#
+eps2 <- material$epsilon#
+lambda <- material$lambda#
+k0 <- 2*pi/lambda#
+N.t <- 3#
+N.k <-length(k0)#
+kx <- k0*sin(pi/4)#
 #
-scale <- 2#
-N <- 100#
-xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+kz <- matrix(sqrt(eps1*k0^2 - kx^2) + 0i, nrow = N.k, ncol=N.t, byrow=F)  #
+epsilon <- cbind(rep(eps1, N.k), eps2, eps3)#
 #
-cart2polar <- function(xy){#
-r <- sqrt(xy[, 1]^2+xy[, 2]^2)	#
-angle <- atan(xy[, 2]/xy[, 1])#
-data.frame(r=r, angle=angle)#
+kze <- kz/epsilon
+epsilon%*%k0^2
+k0^2
+outer(epsilon, k0^2)
+keps <- epsilon * matrix(k0^2, nrow = N.k, ncol=N.t, byrow=F)
+keps
+lambda
+#
+foo <- function(kx, d=0.1){#
+	kz <- matrix(sqrt(keps - kx^2 + 0i), nrow = N.k, ncol=N.t, byrow=F)  #
+	kze <- kz / epsilon#
+res <- (kze[, 2] + kze[, 1])*(kze[, 3] + kze[, 1]) + (kze[, 2] - kze[, 1])*(kze[, 3] - kze[, 1])*exp(2i*d*kz[, 1])#
+res#
 }#
 #
-rphi <- cart2polar(xy)#
-remove <- which(rphi$r <= radius)#
+foo(k0*sin(pi/4))
 #
-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))))#
+foo(k0*sin(pi/4))#
 #
-levelplot(z~x*y, xy, cuts = 100,aspect="xy",  xlab=expression("distance / "*mu*m),#
-          ylab=expression("distance / "*mu*m), #
-          colorkey = T, region = TRUE)
+objF <- function(p){#
+	Mod(foo(p[1]+1i*p[2]))#
+}#
 #
+p0 <-  k0*sin(pi/4)#
+optim(p0, objF)->res
 #
-radius <- 0.2#
-with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+permittivity(0.5)
+permittivity
 #
 #
+permittivity(500)
+library(baptMisc)#
+clr()#
 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))
 #
-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#
+eps1 <- 1.5^2#
+eps3 <- 1.5^2#
+eps2 <- material$epsilon#
+lambda <- material$lambda#
 k0 <- 2*pi/lambda#
+N.t <- 3#
+N.k <-length(k0)#
+kx <- k0*sin(pi/4)#
 #
-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)#
+permittivity(500)
 #
-scale <- 2#
-N <- 100#
-xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+permittivity(500)
 #
-cart2polar <- function(xy){#
-r <- sqrt(xy[, 1]^2+xy[, 2]^2)	#
-angle <- atan(xy[, 2]/xy[, 1])#
-data.frame(r=r, angle=angle)#
-}#
+data(AuJChristy)#
+plot.material()$fine.sample->material#
 #
-rphi <- cart2polar(xy)#
-remove <- which(rphi$r <= radius)#
+eps1 <- 1.5^2#
+eps3 <- 1.5^2#
+eps2 <- material$epsilon#
+lambda <- material$lambda#
+k0 <- 2*pi/lambda#
+N.t <- 3#
+N.k <-length(k0)#
+kx <- k0*sin(pi/4)#
 #
-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)
+permittivity(500)
+setwd('/Users/baptiste/Documents/R/photonics/pkg/mie/R')
+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#
+}
 #
-xy$z <- 0#
-xy$z[-remove] <- c(apply(my.field, 2, function(x) Re(x%*%Conj(x))))#
+permittivity(w=500)
 #
-xy$z[remove] <- NA#
+foo <- function(kx, lambda=0.5, d=0.1){#
+	k0 <- 2*pi/lambda#
+	eps2 <- permittivity(lambda*1e3)#
+	epsilon <- c(1.5^2, eps2, 1.5^2)#
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)#
+	kze <- kz / epsilon#
+	res <- (kze[2] + kze[1])*(kze[3] + kze[1]) + (kze[2] - kze[1])*(kze[3] - kze[1]) * exp(2i*d*kz[, 1])#
+	res#
+}#
 #
-levelplot(z~x*y, xy, cuts = 100,aspect="xy",  xlab=expression("distance / "*mu*m),#
-          ylab=expression("distance / "*mu*m), #
-          colorkey = T, region = TRUE)
+foo(k0*sin(pi/4)[1])
 #
 #
-library(mie)#
-data(AuJChristy)#
-plot.material()$fine.sample->material#
+foo <- function(kx, lambda=0.5, d=0.1){#
+	k0 <- 2*pi/lambda#
+	eps2 <- permittivity(w=lambda*1e3)#
+	epsilon <- c(1.5^2, eps2, 1.5^2)#
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)#
+	kze <- kz / epsilon#
+	res <- (kze[2] + kze[1])*(kze[3] + kze[1]) + (kze[2] - kze[1])*(kze[3] - kze[1]) * exp(2i*d*kz[, 1])#
+	res#
+}#
 #
-radius <- 0.2#
-with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))#
+foo(k0*sin(pi/4)[1])
 #
+foo <- function(kx, lambda=0.5, d=0.1){#
+	k0 <- 2*pi/lambda#
+	eps2 <- permittivity(w=lambda*1e3)#
+	epsilon <- c(1.5^2, eps2, 1.5^2)#
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)#
+	kze <- kz / epsilon#
+	res <- (kze[2] + kze[1])*(kze[3] + kze[1]) + (kze[2] - kze[1])*(kze[3] - kze[1]) * exp(2i*d*kz[1])#
+	res#
+}#
 #
+foo(k0*sin(pi/4)[1])
+lambda=0.5
+	k0 <- 2*pi/lambda#
+	eps2 <- permittivity(w=lambda*1e3)
+	epsilon <- c(1.5^2, eps2, 1.5^2)#
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)#
+	kze <- kz / epsilon
+k0
+kx
+kx=1
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)
+kz
+	kze <- kz / epsilon
+k0*sin(pi/4)[1]
 #
+foo <- function(kx=1, lambda=0.5, d=0.1){#
+	k0 <- 2*pi/lambda#
+	eps2 <- permittivity(w=lambda*1e3)#
+	epsilon <- c(1.5^2, eps2, 1.5^2)#
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)#
+	kze <- kz / epsilon#
+	res <- (kze[2] + kze[1])*(kze[3] + kze[1]) + (kze[2] - kze[1])*(kze[3] - kze[1]) * exp(2i*d*kz[1])#
+	res#
+}#
 #
-lambda <- 1.2#
-lambda <- 0.5#
-k0 <- 2*pi/lambda#
+foo(k0*sin(pi/4)[1])
 #
-eps <- permittivity(wav=lambda)#
-n <- Re(sqrt(eps))#
-k <- Im(sqrt(eps))#
+p0 <-  k0*sin(pi/4)#
+optim(p0, objF)->res
 #
-ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
-cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+objF <- function(p){#
+	Mod(foo(kx=p[1]+1i*p[2]))#
+}#
 #
-scale <- 2#
-N <- 40#
-xy <- expand.grid(x=seq(-scale*radius, scale*radius, length=N), y=seq(-scale*radius, scale*radius, length=N))#
+p0 <-  k0*sin(pi/4)#
+optim(p0, objF)->res
 #
-cart2polar <- function(xy){#
-r <- sqrt(xy[, 1]^2+xy[, 2]^2)	#
-angle <- atan(xy[, 2]/xy[, 1])#
-data.frame(r=r, angle=angle)#
+p0 <-  c(k0*sin(pi/4), 0)#
+optim(p0, objF)->res
+res
+#
+objF <- function(p, lambda){#
+	Mod(foo(kx=p[1]+1i*p[2], lambda=lambda))#
 }#
 #
-rphi <- cart2polar(xy)#
-remove <- which(rphi$r <= radius)#
+p0 <-  c(k0*sin(pi/4), 0)
+optim(p0, objF)$par->res
+?optim
+material$lambda
 #
-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))#
+fooOptim <- function(lambda){#
+optim(p0, objF, lambda=lambda)$par->res#
+}#
 #
-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))))#
+sapply(material$lambda, fooOptim)
 #
+plot.material(fine=10)$fine.sample->material
+material
 #
+test <- sapply(material$lambda, fooOptim)
+matplot(material$lambda,test)
+matplot(material$lambda,t(test))
 #
-p12 <- levelplot(z12~x*y, xy, cuts = 100,aspect="xy",  xlab=expression("distance / "*mu*m),#
-          ylab=expression("distance / "*mu*m), #
-          colorkey = T, region = TRUE)#
+test[1, ]->kxr#
+test[2, ]->kxi#
 #
-p065 <- levelplot(z065~x*y, xy, cuts = 100,aspect="xy",  xlab=expression("distance / "*mu*m),#
-          ylab=expression("distance / "*mu*m), #
-          colorkey = T, region = TRUE)#
+par(mfrow=c(2, 1))#
+plot(kxr, 2*pi/material$lambda)#
+plot(material$lambda, kxi)
 #
-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
+par(mfrow=c(1, 2))#
+plot(kxr, 2*pi/material$lambda)#
+plot(material$lambda, kxi)
 #
-myLay <- grid.layout(3, 1, widths=unit(1, "inches"),#
-                        heights=unit(0.25, "npc"),#
-                        just=just)#
+foo <- function(kx=1, lambda=0.5, d=0.1){#
+	k0 <- 2*pi/lambda#
+	eps2 <- permittivity(w=lambda*1e3)#
+	epsilon <- c(1.5^2, eps2, 1.0^2)#
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)#
+	kze <- kz / epsilon#
+	res <- (kze[2] + kze[1])*(kze[3] + kze[1]) + (kze[2] - kze[1])*(kze[3] - kze[1]) * exp(2i*d*kz[1])#
+	res#
+}#
 #
-pushViewport(viewport(layout=myLay)))
+foo(k0*sin(pi/4)[1])#
 #
-myLay <- grid.layout(3, 1, widths=unit(1, "inches"),#
-                        heights=unit(0.25, "npc"))#
+objF <- function(p, lambda){#
+	Mod(foo(kx=p[1]+1i*p[2], lambda=lambda))#
+}#
 #
-pushViewport(viewport(layout=myLay)))
+p0 <-  c(k0*sin(pi/4), 0)#
+fooOptim <- function(lambda){#
+optim(p0, objF, lambda=lambda)$par->res#
+}#
 #
-myLay <- grid.layout(3, 1, widths=unit(1, "inches"),#
-                        heights=unit(0.25, "npc"))#
+test <- sapply(material$lambda, fooOptim)#
 #
-pushViewport(viewport(layout=myLay))
-myLay <- grid.layout(3, 1, widths=unit(1, "inches"),#
-                        heights=unit(0.25, "npc"))#
+test[1, ]->kxr#
+test[2, ]->kxi#
 #
-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"))#
+par(mfrow=c(1, 2))#
+plot(kxr, 2*pi/material$lambda)#
+plot(material$lambda, kxi)
+library(baptMisc)#
+clr()#
+library(mie)#
+library(grid)#
+data(AuJChristy)#
+plot.material(fine=50)$fine.sample->material#
 #
-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"))#
+eps1 <- 1.5^2#
+eps3 <- 1.5^2#
+eps2 <- material$epsilon#
+lambda <- material$lambda#
+k0 <- 2*pi/lambda#
+N.t <- 3#
+N.k <-length(k0)#
+kx <- k0*sin(pi/4)#
 #
-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"))#
+permittivity(w=500)#
 #
-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"))#
+epsilon <- cbind(rep(eps1, N.k), eps2, eps3)#
+keps <- epsilon * matrix(k0^2, nrow = N.k, ncol=N.t, byrow=F)#
 #
-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()
+kz <- matrix(sqrt(eps1*k0^2 - kx^2) + 0i, nrow = N.k, ncol=N.t, byrow=F)  #
 #
-library(lattice)#
-myLay <- grid.layout(1, 3, widths=unit(2, "inches"),#
-                        heights=unit(2, "inches"))#
+kze <- kz/epsilon#
 #
-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)
+foo <- function(kx=1, lambda=0.5, d=0.1){#
+	k0 <- 2*pi/lambda#
+	eps2 <- permittivity(w=lambda*1e3)#
+	epsilon <- c(1.5^2, eps2, 1.0^2)#
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)#
+	kze <- kz / epsilon#
+	res <- (kze[2] + kze[1])*(kze[3] + kze[1]) + (kze[2] - kze[1])*(kze[3] - kze[1]) * exp(2i*d*kz[1])#
+	res#
+}#
 #
-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)))))#
+foo(k0*sin(pi/4)[1])#
 #
-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)))#
-})
+objF <- function(p, lambda){#
+	Mod(foo(kx=p[1]+1i*p[2], lambda=lambda))#
+}#
 #
-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)))#
-})
+p0 <-  c(k0*sin(pi/4), 0)#
+fooOptim <- function(lambda){#
+optim(p0, objF, lambda=lambda)$par->res#
+}#
 #
+test <- sapply(material$lambda, fooOptim)#
 #
+test[1, ]->kxr#
+test[2, ]->kxi#
+#
+par(mfrow=c(1, 2))#
+plot(kxr, 2*pi/material$lambda)#
+plot(material$lambda, kxi)
+library(baptMisc)#
+clr()#
 library(mie)#
 library(grid)#
 data(AuJChristy)#
-plot.material()$fine.sample->material#
+plot.material(fine=10)$fine.sample->material#
 #
-radius <- 0.2#
-with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))#
+eps1 <- 1.5^2#
+eps3 <- 1.5^2#
+eps2 <- material$epsilon#
+lambda <- material$lambda#
+k0 <- 2*pi/lambda#
+N.t <- 3#
+N.k <-length(k0)#
+kx <- k0*sin(pi/4)#
 #
 #
+permittivity(w=500)#
 #
 #
-lambda <- 1.2#
-lambda <- 0.5#
-k0 <- 2*pi/lambda#
+epsilon <- cbind(rep(eps1, N.k), eps2, eps3)#
+keps <- epsilon * matrix(k0^2, nrow = N.k, ncol=N.t, byrow=F)#
 #
-eps <- permittivity(wav=lambda)#
-n <- Re(sqrt(eps))#
-k <- Im(sqrt(eps))#
+kz <- matrix(sqrt(eps1*k0^2 - kx^2) + 0i, nrow = N.k, ncol=N.t, byrow=F)  #
 #
-ab <- Mie.ab(m=n+1i*k, x=k0 * radius)#
-cd <- Mie.cd(m=n+1i*k, x=k0 * radius)#
+kze <- kz/epsilon#
 #
-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)#
+foo <- function(kx=1, lambda=0.5, d=0.1){#
+	k0 <- 2*pi/lambda#
+	eps2 <- permittivity(w=lambda*1e3)#
+	epsilon <- c(1.5^2, eps2, 1.0^2)#
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)#
+	kze <- kz / epsilon#
+	res <- (kze[2] + kze[1])*(kze[3] + kze[1]) + (kze[2] - kze[1])*(kze[3] - kze[1]) * exp(2i*d*kz[1])#
+	res#
+}#
 #
-p065 <- levelplot(z065~x*y, xy, cuts = 100,aspect="xy",  xlab=expression("distance / "*mu*m),#
-          ylab=expression("distance / "*mu*m), #
-          colorkey = T, region = TRUE)#
+foo(k0*sin(pi/4)[1])#
 #
-p05 <- levelplot(z05~x*y, xy, cuts = 100,aspect="xy",  xlab=expression("distance / "*mu*m),#
-          ylab=expression("distance / "*mu*m), #
-          colorkey = T, region = TRUE)#
+objF <- function(p, lambda){#
+	Mod(foo(kx=p[1]+1i*p[2], lambda=lambda))#
+}#
 #
-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(...){#
+fooOptim <- function(lambda){#
 	#
-			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)))#
-})
+optim(c(2*pi/lambda, 0), objF, lambda=lambda)$par->res#
+}#
 #
-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)))#
-})
+test <- sapply(material$lambda, fooOptim)#
 #
-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
+test[1, ]->kxr#
+test[2, ]->kxi#
 #
-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)))#
-})
+par(mfrow=c(1, 2))#
+plot(kxr, 2*pi/material$lambda, t="b")#
+plot(material$lambda, kxi, t="b")
+library(baptMisc)#
+clr()#
+library(mie)#
+library(grid)#
+data(AuJChristy)#
+plot.material(fine=20)$fine.sample->material#
 #
-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)))#
-})
+eps1 <- 1.5^2#
+eps3 <- 1.5^2#
+eps2 <- material$epsilon#
+lambda <- material$lambda#
+k0 <- 2*pi/lambda#
+N.t <- 3#
+N.k <-length(k0)#
+kx <- k0*sin(pi/4)#
 #
-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)))#
-})
+permittivity(w=500)#
 #
-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)))#
-})#
+epsilon <- cbind(rep(eps1, N.k), eps2, eps3)#
+keps <- epsilon * matrix(k0^2, nrow = N.k, ncol=N.t, byrow=F)#
 #
-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
+kz <- matrix(sqrt(eps1*k0^2 - kx^2) + 0i, nrow = N.k, ncol=N.t, byrow=F)  #
 #
-quartz(file="fieldPlot.pdf", height=3)#
-print(field.plot)#
-dev.off()
+kze <- kz/epsilon#
 #
-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()
 #
+foo <- function(kx=1, lambda=0.5, d=0.1){#
+	k0 <- 2*pi/lambda#
+	eps2 <- permittivity(w=lambda*1e3)#
+	epsilon <- c(1.5^2, eps2, 1.0^2)#
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)#
+	kze <- kz / epsilon#
+	res <- (kze[2] + kze[1])*(kze[3] + kze[1]) + (kze[2] - kze[1])*(kze[3] - kze[1]) * exp(2i*d*kz[1])#
+	res#
+}#
 #
-quartz(type="png", file="MieNearFieldPlot.png", height=5, width=7, bg="white", dpi=200)#
+foo(k0*sin(pi/4)[1])#
 #
-print(field.plot)#
-dev.off()
+objF <- function(p, lambda){#
+	Mod(foo(kx=p[1]+1i*p[2], lambda=lambda))#
+}#
 #
-quartz(type="png", file="MieNearFieldPlot.png", height=4, width=7, bg="white", dpi=200)#
 #
-print(field.plot)#
-dev.off()
+fooOptim <- function(lambda){#
+	#
+optim(c(2*pi/lambda, 0), objF, lambda=lambda)$par->res#
+}#
 #
-radius <- 0.1#
-with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+test <- sapply(material$lambda, fooOptim)#
 #
-radius <- 0.2#
+test[1, ]->kxr#
+test[2, ]->kxi#
 #
+par(mfrow=c(1, 2))#
+plot(kxr, 2*pi/material$lambda, t="b")#
+plot(material$lambda, kxi, t="b")
+test
+library(baptMisc)#
+clr()#
+library(mie)#
+library(grid)#
+data(AuJChristy)#
+plot.material(fine=20)$fine.sample->material#
 #
+eps1 <- 1.5^2#
+eps3 <- 1.5^2#
+eps2 <- material$epsilon#
+lambda <- material$lambda#
+k0 <- 2*pi/lambda#
+N.t <- 3#
+N.k <-length(k0)#
+kx <- k0*sin(pi/4)#
 #
-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#
+permittivity(w=500)#
 #
-radius <- 0.2#
 #
+epsilon <- cbind(rep(eps1, N.k), eps2, eps3)#
+keps <- epsilon * matrix(k0^2, nrow = N.k, ncol=N.t, byrow=F)#
 #
+kz <- matrix(sqrt(eps1*k0^2 - kx^2) + 0i, nrow = N.k, ncol=N.t, byrow=F)  #
 #
-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#
+kze <- kz/epsilon#
 #
-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#
+foo <- function(kx=1, lambda=0.5, d=0.05){#
+	k0 <- 2*pi/lambda#
+	eps2 <- permittivity(w=lambda*1e3)#
+	epsilon <- c(1.5^2, eps2, 1.0^2)#
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)#
+	kze <- kz / epsilon#
+	res <- (kze[2] + kze[1])*(kze[3] + kze[1]) + (kze[2] - kze[1])*(kze[3] - kze[1]) * exp(2i*d*kz[1])#
+	res#
+}#
 #
-radius <- 0.05#
+foo(k0*sin(pi/4)[1])#
 #
+objF <- function(p, lambda){#
+	Mod(foo(kx=p[1]+1i*p[2], lambda=lambda))#
+}#
 #
 #
-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#
+fooOptim <- function(lambda){#
+	#
+optim(c(2*pi/lambda, 0), objF, lambda=lambda)$par->res#
+}#
 #
-radius <- 0.05#
+test <- sapply(material$lambda, fooOptim)#
 #
+test[1, ]->kxr#
+test[2, ]->kxi#
 #
+par(mfrow=c(1, 2))#
+plot(kxr, 2*pi/material$lambda, t="b")#
+plot(material$lambda, kxi, t="b")
+library(baptMisc)#
+clr()#
+library(mie)#
+library(grid)#
+data(AuJChristy)#
+plot.material(fine=20)$fine.sample->material#
 #
-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#
+eps1 <- 1.5^2#
+eps3 <- 1.5^2#
+eps2 <- material$epsilon#
+lambda <- material$lambda#
+k0 <- 2*pi/lambda#
+N.t <- 3#
+N.k <-length(k0)#
+kx <- k0*sin(pi/4)#
 #
-radius <- 0.05#
 #
+permittivity(w=500)#
 #
 #
-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#
+epsilon <- cbind(rep(eps1, N.k), eps2, eps3)#
+keps <- epsilon * matrix(k0^2, nrow = N.k, ncol=N.t, byrow=F)#
 #
-radius <- 0.05#
+kz <- matrix(sqrt(eps1*k0^2 - kx^2) + 0i, nrow = N.k, ncol=N.t, byrow=F)  #
 #
+kze <- kz/epsilon#
 #
 #
-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#
 #
+foo <- function(kx=1, lambda=0.5, d=0.05){#
+	k0 <- 2*pi/lambda#
+	eps2 <- permittivity(w=lambda*1e3)#
+	epsilon <- c(1.5^2, eps2, 1.0^2)#
+	kz <- sqrt(epsilon*k0^2 - kx^2 + 0i)#
+	kze <- kz / epsilon#
+	res <- (kze[2] + kze[1])*(kze[3] + kze[1]) + (kze[2] - kze[1])*(kze[3] - kze[1]) * exp(2i*d*kz[2])#
+	res#
+}#
 #
+foo(k0*sin(pi/4)[1])#
 #
+objF <- function(p, lambda){#
+	Mod(foo(kx=p[1]+1i*p[2], lambda=lambda))#
+}#
 #
-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#
+fooOptim <- function(lambda){#
+	#
+optim(c(2*pi/lambda, 0), objF, lambda=lambda)$par->res#
+}#
 #
+test <- sapply(material$lambda, fooOptim)#
 #
+test[1, ]->kxr#
+test[2, ]->kxi#
 #
+par(mfrow=c(1, 2))#
+plot(kxr, 2*pi/material$lambda, t="b")#
+plot(material$lambda, kxi, t="b")
+setwd('/Users/baptiste/Documents/R')
 #
-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#
+objF <- function(p, lambda){#
 #
-radius <- 0.04#
+	res <- Mod(foo(kx=p[1]+1i*p[2], lambda=lambda))#
+	if(p[1]^2+p[2]^2 > (2*pi/lambda)^2) res <- res* 1e9	#
+res#
+}#
 #
 #
+fooOptim <- function(lambda){#
+	#
+optim(c(2*pi/lambda, 0), objF, lambda=lambda)$par->res#
+}#
 #
+test <- sapply(material$lambda, fooOptim)#
 #
-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#
+test[1, ]->kxr#
+test[2, ]->kxi#
 #
-radius <- 0.04#
+par(mfrow=c(1, 2))#
+plot(kxr, 2*pi/material$lambda, t="b")#
+plot(material$lambda, kxi, t="b")
+load("angle520S0-5.rda")#
+sizes.spectra <- pitch520#
+load("angle550S0-10.rda")#
+sizes.spectra <- pitch550
+load("angle520S0-5.rda")#
+sizes.spectra -> pitch520#
+load("angle550S0-10.rda")#
+sizes.spectra -> pitch550
 #
+wavelength <- pitch550 at wavelength#
+intensity <- as.data.frame(cbind(pitch520 at intensity, pitch550 at intensity))#
+labels <- as.data.frame(rbind(pitch520 at intensity, pitch550 at intensity))#
+labels$pitch <- factor(rep(c(0.52, 0.55), each=length(pitch520 at intensity)))
 #
 #
+wavelength <- pitch550 at wavelength#
+intensity <- as.data.frame(cbind(pitch520 at intensity, pitch550 at intensity))#
+labels <- as.data.frame(rbind(pitch520 at labels, pitch550 at labels))#
+labels$pitch <- factor(rep(c(0.52, 0.55), each=length(pitch520 at intensity)))
 #
-with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1.5), t="l", ylab=~sigma))
+wavelength <- pitch550 at wavelength#
+intensity <- as.data.frame(cbind(pitch520 at intensity, pitch550 at intensity))#
+labels <- as.data.frame(rbind(pitch520 at labels, pitch550 at labels))#
+labels$pitch <- factor(rep(c(0.52, 0.55), each=length(pitch520 at labels)))
 #
-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#
+labels$pitch <- factor(rep(c(0.52, 0.55), each=ncol(pitch520 at labels)))
+nrow(pitch520 at labels)
 #
-radius <- 0.04#
+labels$pitch <- factor(rep(c(0.52, 0.55), each=nrow(pitch520 at labels)))
+pitch520 at labels
 #
+labels$pitch <- factor(rep(0.52, nrow(pitch520 at labels)), rep(0.55, nrow(pitch550 at labels)))
 #
+labels$pitch <- factor(rep(0.52, ncol(pitch520 at labels)), rep(0.55, ncol(pitch550 at labels)))
 #
+both <- spectra(wavelength=wavelength, intensity=intensity, labels=labels, setup=pitch550 at setup)
+load("angle520S0-5.rda")#
+sizes.spectra -> pitch520#
+load("angle550S0-10.rda")#
+sizes.spectra -> pitch550#
 #
-mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
-	matplot(lambda, mdata, t="l", ylab=~sigma)#
-})#
+library(spectra)#
+library(lattice)#
+ls()#
 #
-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#
+wavelength <- pitch550 at wavelength#
+intensity <- as.data.frame(cbind(pitch520 at intensity, pitch550 at intensity))#
+labels <- as.data.frame(rbind(pitch520 at labels, pitch550 at labels))#
+labels$pitch <- factor(rep(0.52, ncol(pitch520 at labels)), rep(0.55, ncol(pitch550 at labels)))#
 #
-radius <- 0.04#
+both <- spectra(wavelength=wavelength, intensity=intensity, labels=labels, setup=pitch550 at setup)
+both
 #
+cext.melt <- as.data.frame(both)#
+diff <- 0.52 * 1.5#
+dark.col <- colorRampPalette(c(brewer.pal(8, "Set1")[1], brewer.pal(8, "Set1")[2]))#
+pastel.col <- colorRampPalette(c(brewer.pal(8, "Pastel1")[1], brewer.pal(8, "Pastel1")[2]))#
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+				panel.segments(diff, 0, diff, 1, lty=2, col="grey30")#
+				panel.segments(diff/sqrt(2), 0, diff/sqrt(2), 0.1, lty=2, col="grey30")#
+	},#
+	  		panel.groups = function(..., group.number) { #
 #
+	  				 panel.xyplot(...)#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
+labels
 #
+labels$pitch <- factor(rep(0.52, ncol(pitch520 at intensity)), rep(0.55, ncol(pitch550 at intensity)))
+ ncol(pitch520 at intensity)
+ncol(pitch550 at intensity)
 #
-mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
-	matplot(lambda, mdata, t="l", ylab=~sigma);#
-	mdata#
-})#
+labels$pitch <- factor(c(rep(0.52, ncol(pitch520 at intensity)), rep(0.55, ncol(pitch550 at intensity))))
 #
-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#
+wavelength <- pitch550 at wavelength#
+intensity <- as.data.frame(cbind(pitch520 at intensity, pitch550 at intensity))#
+labels <- as.data.frame(rbind(pitch520 at labels, pitch550 at labels))#
+labels$pitch <- factor(c(rep(0.52, ncol(pitch520 at intensity)), rep(0.55, ncol(pitch550 at intensity))))#
 #
-radius <- 0.04#
+both <- spectra(wavelength=wavelength, intensity=intensity, labels=labels, setup=pitch550 at setup)#
 #
+cext.melt <- as.data.frame(both)#
+diff <- 0.52 * 1.5#
+dark.col <- colorRampPalette(c(brewer.pal(8, "Set1")[1], brewer.pal(8, "Set1")[2]))#
+pastel.col <- colorRampPalette(c(brewer.pal(8, "Pastel1")[1], brewer.pal(8, "Pastel1")[2]))#
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+				panel.segments(diff, 0, diff, 1, lty=2, col="grey30")#
+				panel.segments(diff/sqrt(2), 0, diff/sqrt(2), 0.1, lty=2, col="grey30")#
+	},#
+	  		panel.groups = function(..., group.number) { #
 #
+	  				 panel.xyplot(...)#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
 #
+xyplot(value~wavelength|pitch, data=cext.melt, subset= theta <= 10, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+				panel.segments(diff, 0, diff, 1, lty=2, col="grey30")#
+				panel.segments(diff/sqrt(2), 0, diff/sqrt(2), 0.1, lty=2, col="grey30")#
+	},#
+	  		panel.groups = function(..., group.number) { #
 #
-mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
-	matplot(lambda, mdata, t="l", ylab=~sigma);#
-	cbind(lambda, mdata)#
-})#
+	  				 panel.xyplot(...)#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, subset= theta <= 3, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+				panel.segments(diff, 0, diff, 1, lty=2, col="grey30")#
+				panel.segments(diff/sqrt(2), 0, diff/sqrt(2), 0.1, lty=2, col="grey30")#
+	},#
+	  		panel.groups = function(..., group.number) { #
 #
-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#
+	  				 panel.xyplot(...)#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
 #
-radius <- 0.04#
+xyplot(value~wavelength|pitch, data=cext.melt, subset= factor2num(theta) <= 3, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+				panel.segments(diff, 0, diff, 1, lty=2, col="grey30")#
+				panel.segments(diff/sqrt(2), 0, diff/sqrt(2), 0.1, lty=2, col="grey30")#
+	},#
+	  		panel.groups = function(..., group.number) { #
 #
+	  				 panel.xyplot(...)#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )
+factor2num(cext.melt$theta)
 #
+xyplot(value~wavelength|pitch, data=cext.melt, subset= factor2num(theta) ==0, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+				panel.segments(diff, 0, diff, 1, lty=2, col="grey30")#
+				panel.segments(diff/sqrt(2), 0, diff/sqrt(2), 0.1, lty=2, col="grey30")#
+	},#
+	  		panel.groups = function(..., group.number) { #
 #
+	  				 panel.xyplot(...)#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )
 #
-mdata <- with(material, {mdata <- Mie.spectra(lambda, epsilon, radius=radius, medium=1.5);#
-	matplot(lambda, mdata, t="l", ylab=~sigma);#
-	cbind(lambda, mdata)#
-})#
+xyplot(value~wavelength|pitch, data=cext.melt, subset= factor2num(theta) <=0.1, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+				panel.segments(diff, 0, diff, 1, lty=2, col="grey30")#
+				panel.segments(diff/sqrt(2), 0, diff/sqrt(2), 0.1, lty=2, col="grey30")#
+	},#
+	  		panel.groups = function(..., group.number) { #
 #
-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#
+	  				 panel.xyplot(...)#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )
+subset(cex.melt,  factor2num(theta) <=0.1)
 #
-radius <- 0.04#
+cext.melt <- subset(cex.melt,  factor2num(cext.melt$theta) <=0.1)
+subset
 #
+cext.melt <- subset(both,  factor2num(theta) <=0.1)
 #
+cext.melt <- subset(both,  factor2num(theta) <=0.1)#
+diff <- 0.52 * 1.5#
+dark.col <- colorRampPalette(c(brewer.pal(8, "Set1")[1], brewer.pal(8, "Set1")[2]))#
+pastel.col <- colorRampPalette(c(brewer.pal(8, "Pastel1")[1], brewer.pal(8, "Pastel1")[2]))#
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+				panel.segments(diff, 0, diff, 1, lty=2, col="grey30")#
+				panel.segments(diff/sqrt(2), 0, diff/sqrt(2), 0.1, lty=2, col="grey30")#
+	},#
+	  		panel.groups = function(..., group.number) { #
 #
+	  				 panel.xyplot(...)#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )
 #
-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))#
-})#
+cext.melt <- as.data.frame(subset(both,  factor2num(theta) <=0.1))#
+diff <- 0.52 * 1.5#
+dark.col <- colorRampPalette(c(brewer.pal(8, "Set1")[1], brewer.pal(8, "Set1")[2]))#
+pastel.col <- colorRampPalette(c(brewer.pal(8, "Pastel1")[1], brewer.pal(8, "Pastel1")[2]))#
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+				panel.segments(diff, 0, diff, 1, lty=2, col="grey30")#
+				panel.segments(diff/sqrt(2), 0, diff/sqrt(2), 0.1, lty=2, col="grey30")#
+	},#
+	  		panel.groups = function(..., group.number) { #
 #
-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#
+	  				 panel.xyplot(...)#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
 #
-radius <- 0.04#
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+			#
+	},#
+	  		panel.groups = function(..., group.number) { #
 #
+	  				 panel.xyplot(...)#
+						panel.segments(diff[group.number], 0, diff[group.number], 1, lty=2, col="grey30")#
+						panel.segments(diff[group.number]/sqrt(2), 0, diff[group.number]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )
 #
+both <- spectra(wavelength=wavelength, intensity=intensity, labels=labels, setup=pitch550 at setup)#
 #
+cext.melt <- as.data.frame(both)#
+cext.melt <- as.data.frame(subset(both,  factor2num(theta) <=0.1))#
+diff <- c(0.52, 0.55) * 1.5#
+dark.col <- colorRampPalette(c(brewer.pal(8, "Set1")[1], brewer.pal(8, "Set1")[2]))#
+pastel.col <- colorRampPalette(c(brewer.pal(8, "Pastel1")[1], brewer.pal(8, "Pastel1")[2]))#
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+			#
+	},#
+	  		panel.groups = function(..., group.number) { #
+	  				 panel.xyplot(...)#
+						panel.segments(diff[group.number], 0, diff[group.number], 1, lty=2, col="grey30")#
+						panel.segments(diff[group.number]/sqrt(2), 0, diff[group.number]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
 #
-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))#
-})#
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+			#
+	},#
+	  		panel.groups = function(..., group.number) { #
+	  				 panel.xyplot(...)#
+	print(group.number)#
+						panel.segments(diff[group.number], 0, diff[group.number], 1, lty=2, col="grey30")#
+						panel.segments(diff[group.number]/sqrt(2), 0, diff[group.number]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+			#
+	},#
+	  		panel.groups = function(..., group.number) { #
+	  				 panel.xyplot(...)#
+	#
+	ind <- panel.number()#
+						panel.segments(diff[ind], 0, diff[ind], 1, lty=2, col="grey30")#
+						panel.segments(diff[ind]/sqrt(2), 0, diff[ind]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
 #
-names(mdata) <- c("wavelength", "cext", "cabs", "csat")#
-head(mdata)#
-#
-export2excel(mdata, file="chris.txt")
-system("head chris.txt")
+lambda <- 0.5#
+k0 <- 2*pi/lambda#
 library(constants)#
-data(AgPalik)#
-plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+data(Constants)#
+omega <- Constants$cel * k0#
 #
+limits <- 1#
+N <- 300#
 #
-radius <- 0.04#
+xy <- expand.grid(x=seq(-limits, limits, length=N), y=seq(-limits, limits, 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)#
 #
+kr <- rphi$r * k0#
+time <- 100#
+phase <- kr - omega*time#
 #
-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))#
-})#
+E.theta <- ifelse(kr==0, 0, cos(rphi$angle)*exp(1i*phase)/kr)#
 #
-names(mdata) <- c("wavelength", "cext", "cabs", "csat")#
-head(mdata)#
+Int <- Re(E.theta)#
 #
-export2excel(mdata, file="chris.txt")
-library(constants)#
-data(AgPalik)#
-plot.material("AgPalik", plot=F, range=c(0.3, 0.9))$fine.sample->material#
+xy$z <- Int#
 #
+contourplot(z~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
-radius <- 0.04#
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
+	panel.levelplot(x, y, z, contour=T, region=F, at=pretty(c(range(xy$z)[1],range(xy$z)[2]/20), 10),#
+					col=grey(0.3),lty=3, ...)#
+	#
+},           colorkey = F, region = T)
 #
-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))#
-})#
+lambda <- 0.5#
+k0 <- 2*pi/lambda#
+library(constants)#
+data(Constants)#
+omega <- Constants$cel * k0#
 #
-names(mdata) <- c("wavelength", "cext", "cabs", "csat")#
-head(mdata)#
+limits <- 10#
+N <- 200#
 #
-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#
+xy <- expand.grid(x=seq(-limits, limits, length=N), y=seq(-limits, limits, 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)#
+}#
 #
-radius <- 0.04#
+rphi <- cart2polar(xy)#
 #
-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))#
-})#
+kr <- rphi$r * k0#
+time <- 100#
+phase <- kr - omega*time#
 #
-names(mdata) <- c("wavelength", "cext", "cabs", "csat")#
-head(mdata)#
+E.theta <- ifelse(kr==0, 0, cos(rphi$angle)*exp(1i*phase)/kr)#
 #
-export2excel(mdata, file="chris.txt")
-radius <- 0.05#
+Int <- Re(E.theta)#
 #
+xy$z <- Int#
 #
+contourplot(z~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
-with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
+	panel.levelplot(x, y, z, contour=T, region=F, at=pretty(c(range(xy$z)[1],range(xy$z)[2]/20), 10),#
+					col=grey(0.3),lty=3, ...)#
+	#
+},           colorkey = F, region = T)
 #
-library(mie)#
-library(grid)#
-data(AuJChristy)#
-plot.material()$fine.sample->material#
+p1 <- contourplot(z~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
-radius <- 0.2#
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
+	#
+},           colorkey = F, region = T)#
 #
+p1
+show.settings()
 #
+trellis.par.get("regions")
 #
-with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))#
-radius <- 0.05#
+update(p1, par.settings=list(regions=list(col=c(grey(0.1), grey(0.5)))))
 #
+update(p1, par.settings=list(regions=list(col=c(grey(seq(0, 1, len=20))))))
 #
+update(p1, par.settings=list(regions=list(col=c(grey(seq(0, 1, len=100))))))
 #
-with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.5)))
 #
-detach(package:constants)
-plot.material
 #
-plot.material()$fine.sample->material#
+p1 <- contourplot(z~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
-radius <- 0.2#
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
+	#
+},           colorkey = F, region = T)#
 #
 #
 #
-with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))#
-radius <- 0.05#
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.5)))
 #
+p1 <- contourplot(z~x*y, rbind(xy, shifty()), cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
+	#
+},           colorkey = F, region = T)#
 #
-with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
-radius <- 0.1#
 #
 #
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.5)))
 #
-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)#
-}#
+p1 <- contourplot(z~x*y, rbind(xy, shifty(xy)), cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
-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){#
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
 	#
-	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)#
-}#
+},           colorkey = F, region = T)
 #
-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))#
+shifty <- function(xy, dy=1){#
+	with(xy, y=y-dy)#
 }#
 #
-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))#
-}#
+rbind(xy, shifty())->xy2#
 #
-oneField()#
 #
+p1 <- contourplot(z~x*y, xy2, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
+	#
+},           colorkey = F, region = T)#
 #
-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)))#
-})#
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.5)))
 #
-field.plot
+rbind(xy, shifty(xy))->xy2#
 #
 #
-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(...){#
+p1 <- contourplot(z~x*y, xy2, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
+#
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
 	#
-			panel.levelplot(...)#
-			#
-			grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
-})#
+},           colorkey = F, region = T)#
 #
-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
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.5)))
+xy
+with(xy, y=y-dy)
 #
+shifty <- function(xy, dy=1){#
+	xy$y <- xy$y-dy)#
+}#
 #
-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)#
+rbind(xy, shifty(xy))->xy2
 #
+shifty <- function(xy, dy=1){#
+	xy$y <- xy$y-dy#
+}#
 #
-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)))#
-})#
+rbind(xy, shifty(xy))->xy2
 #
-field.plot
+p1 <- contourplot(z~x*y, xy2, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
-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(x, y, z, contour=F, region=T, at,col,  ...)#
 	#
-			panel.levelplot(...)#
-			#
-			grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
-})#
+},           colorkey = F, region = T)#
 #
-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)))#
-})#
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.5)))
 #
-field.plot
 #
-oneField <- function(lambda=0.5, radius=0.1, N=5, scale=2){#
+p1 <- contourplot(z~x*y, xy2, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
+#
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
 	#
-	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))#
-}#
+},           colorkey = F, region = T)#
 #
-oneField()#
 #
 #
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.1)))
 #
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.8)))
 #
-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)#
 #
+shifty <- function(xy, dy=1){#
+	xy$y <- xy$y-dy#
+	xy#
+}#
 #
-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)))#
-})#
+rbind(xy, shifty(xy))->xy2#
 #
-field.plot
 #
+p1 <- contourplot(z~x*y, xy2, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
-oneField <- function(lambda=0.5, radius=0.1, N=10, scale=2){#
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
 	#
-	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))#
-}#
+},           colorkey = F, region = T)#
 #
-oneField()#
 #
 #
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.8)))
 #
 #
-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)#
+p1 <- contourplot(z~x*y, xy2, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),xlim=c(-10, 10), ylim=c(-10, 10), #
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
-#
-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(x, y, z, contour=F, region=T, at,col,  ...)#
 	#
-			panel.levelplot(...)#
-			#
-			grid.circle(name="my.circle", r = unit(0.2, "native"), gp=gpar(col="black", fill=grey(0.5)))#
-})#
+},           colorkey = F, region = T)#
 #
-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))#
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.8)))
+repeat
+?repeat
+}}}
+?
+}
+?for
+})
+)
+#
+shifty <- function(xy, dy=1, N=1){#
+do.call(rbind, lapply(seq_along(1:N), function(ii) {xy$y <- xy$y-ii*dy ; xy}))#
 }#
+shifty(xy)
 #
-oneField()#
+shifty <- function(xy, dy=1, N=1){#
+rbind(xy, do.call(rbind, lapply(seq_along(1:N), function(ii) {xy$y <- xy$y-ii*dy ; xy})))#
+}#
 #
+xy2 <- shifty(xy, N=3)#
 #
+p1 <- contourplot(z~x*y, xy2, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),xlim=c(-10, 10), ylim=c(-10, 10), #
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
+	#
+},           colorkey = F, region = T)#
 #
-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)))#
-})#
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.8)))
 #
-field.plot
+p1 <- contourplot(z~x*y, xy2, cuts = 100, xlab=expression("distance / "*mu*m),xlim=c(-10, 10), ylim=c(-10, 10), #
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
-oneField <- function(lambda=0.5, radius=0.1, N=10, region=1.0){#
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
 	#
-	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))#
-}#
+},           colorkey = F, region = T)#
 #
-oneField()
 #
 #
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.8)))
 #
-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)
+xy2 <- shifty(xy, N=8)#
 #
-oneField(lambda=1)
+p1 <- contourplot(z~x*y, xy2, cuts = 100, xlab=expression("distance / "*mu*m),xlim=c(-10, 10), ylim=c(-10, 10), #
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
-	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#
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
 	#
-	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)
+},           colorkey = F, region = T)#
 #
-	#
-	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)
+update(p1, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.8)))
 #
-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)
+p2 <- contourplot(z~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
+	panel.levelplot(x, y, z, contour=T, region=F, at=pretty(c(range(xy$z)[1],range(xy$z)[2]/20), 10),#
+					col=grey(0.3),lty=3, ...)#
+	#
+},           colorkey = F, region = T)
+p2
 #
-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)
+lambda <- 0.5#
+k0 <- 2*pi/lambda#
+library(constants)#
+data(Constants)#
+omega <- Constants$cel * k0#
 #
-oneField(radius=1)
+limits <- 4#
+N <- 200#
 #
-oneField(radius=0.5)
+xy <- expand.grid(x=seq(-limits, limits, length=N), y=seq(-limits, limits, length=N))#
 #
-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))#
+cart2polar <- function(xy){#
+r <- sqrt(xy[, 1]^2+xy[, 2]^2)	#
+angle <- atan(xy[, 2]/xy[, 1])#
+data.frame(r=r, angle=angle)#
 }#
 #
-oneField(radius=0.5)
+rphi <- cart2polar(xy)#
 #
+kr <- rphi$r * k0#
+time <- 100#
+phase <- kr - omega*time#
 #
+E.theta <- ifelse(kr==0, 0, cos(rphi$angle)*exp(1i*phase)/kr)#
 #
-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)#
+Int <- Re(E.theta)#
 #
-mdf <- do.call(rbind, allFields)
+xy$z <- Int#
 #
-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))#
+shifty <- function(xy, dy=1, N=1){#
+rbind(xy, do.call(rbind, lapply(seq_along(1:N), function(ii) {xy$y <- xy$y-ii*dy ; xy})))#
 }#
 #
-oneField(radius=0.5)#
+xy2 <- shifty(xy, N=8)#
 #
+p1 <- contourplot(z~x*y, xy2, cuts = 100, xlab=expression("distance / "*mu*m),xlim=c(-10, 10), ylim=c(-10, 10), #
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
+	#
+},           colorkey = F, region = T)#
 #
 #
-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)#
+p2 <- contourplot(z~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
 #
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
+	panel.levelplot(x, y, z, contour=T, region=F, at=pretty(c(range(xy$z)[1],range(xy$z)[2]/20), 10),#
+					col=grey(0.3),lty=3, ...)#
+	#
+},           colorkey = F, region = T)#
 #
 #
+update(p2, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=0.8)))
 #
-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)#
+lambda <- 0.5#
+k0 <- 2*pi/lambda#
+library(constants)#
+data(Constants)#
+omega <- Constants$cel * k0#
 #
-mdf <- do.call(rbind, allFields)#
+limits <- 3#
+N <- 200#
 #
+xy <- expand.grid(x=seq(-limits, limits, length=N), y=seq(-limits, limits, length=N))#
 #
-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)))#
-}#
+rphi <- cart2polar(xy)#
 #
-oneField(radius=0.5)
+kr <- rphi$r * k0#
+time <- 100#
+phase <- kr - omega*time#
 #
+E.theta <- ifelse(kr==0, 0, cos(rphi$angle)*exp(1i*phase)/kr)#
 #
-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)#
+Int <- Re(E.theta)#
 #
-mdf <- do.call(rbind, allFields)#
+xy$z <- Int#
 #
 #
-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)))#
-})#
+shifty <- function(xy, dy=1, N=1){#
+rbind(xy, do.call(rbind, lapply(seq_along(1:N), function(ii) {xy$y <- xy$y-ii*dy ; xy})))#
+}#
 #
-field.plot
+xy2 <- shifty(xy, N=8)#
 #
-oneField <- function(lambda=0.5, radius=0.1, N=20, region=1.0){#
+p1 <- contourplot(z~x*y, xy2, cuts = 100, xlab=expression("distance / "*mu*m),xlim=c(-10, 10), ylim=c(-10, 10), #
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
+#
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
 	#
-	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)))#
-}#
+},           colorkey = F, region = T)#
 #
 #
 #
 #
-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(...){#
+p2 <- contourplot(z~x*y, xy, cuts = 100,aspect="xy", xlab=expression("distance / "*mu*m),#
+          ylab=expression("distance / "*mu*m),  pretty=T, labels=F, #
+panel=function(x, y, z, contour, region, at,col,lty,   ...){#
+#
+	panel.levelplot(x, y, z, contour=F, region=T, at,col,  ...)#
+	panel.levelplot(x, y, z, contour=T, region=F, at=pretty(c(range(xy$z)[1],range(xy$z)[2]/20), 10),#
+					col=grey(0.3),lty=3, ...)#
 	#
-			panel.levelplot(...)#
-			#
-			grid.circle(name="my.circle", r = unit(0.05, "native"), gp=gpar(col="black", fill=grey(0.5)))#
-})#
+},           colorkey = F, region = T)#
 #
-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)))#
-})#
+update(p2, par.settings=list(regions=list(col=grey(seq(0, 1, len=100)), alpha=1)))
+setwd('/Users/baptiste/Documents/PhD/modelling')
 #
-library(latticeExtra)#
-useOuterStrips(field.plot)
 #
+update(p2, par.settings=list(regions=list(col=grey(seq(0, 1, len=20)), alpha=1)))
+load("angle520S0-5.rda")#
+sizes.spectra -> pitch520#
+load("angle550S0-10.rda")#
+sizes.spectra -> pitch550#
 #
-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(spectra)#
+library(lattice)#
+ls()#
 #
-library(latticeExtra)#
-useOuterStrips(field.plot)
+wavelength <- pitch550 at wavelength#
+intensity <- as.data.frame(cbind(pitch520 at intensity, pitch550 at intensity))#
+labels <- as.data.frame(rbind(pitch520 at labels, pitch550 at labels))#
+labels$pitch <- factor(c(rep(0.52, ncol(pitch520 at intensity)), rep(0.55, ncol(pitch550 at intensity))))#
 #
-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){#
+both <- spectra(wavelength=wavelength, intensity=intensity, labels=labels, setup=pitch550 at setup)#
+#
+cext.melt <- as.data.frame(both)#
+cext.melt <- as.data.frame(subset(both,  factor2num(theta) <=0.1))#
+diff <- c(0.52, 0.55) * 1.5#
+dark.col <- colorRampPalette(c(brewer.pal(8, "Set1")[1], brewer.pal(8, "Set1")[2]))#
+pastel.col <- colorRampPalette(c(brewer.pal(8, "Pastel1")[1], brewer.pal(8, "Pastel1")[2]))#
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+			#
+	},#
+	  		panel.groups = function(..., group.number) { #
+	  				 panel.xyplot(...)#
 	#
-			panel.levelplot(..., subscripts=subscripts)#
+	ind <- panel.number()#
+						panel.segments(diff[ind], 0, diff[ind], 1, lty=2, col="grey30")#
+						panel.segments(diff[ind]/sqrt(2), 0, diff[ind]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
+pastel.col <- colorRampPalette(c(brewer.pal(8, "Pastel1")[1], brewer.pal(8, "Pastel1")[2]))#
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), between=list(x=0.2),#
+scales=list(axs="i"), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
 			#
-			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.groups = function(..., group.number) { #
+	  				 panel.xyplot(...)#
 	#
-			panel.levelplot(..., subscripts=subscripts)#
+	ind <- panel.number()#
+						panel.segments(diff[ind], 0, diff[ind], 1, lty=2, col="grey30")#
+						panel.segments(diff[ind]/sqrt(2), 0, diff[ind]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30"), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), between=list(x=0.2),#
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
 			#
-			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.groups = function(..., group.number) { #
+	  				 panel.xyplot(...)#
 	#
-			panel.levelplot(..., subscripts=subscripts)#
+	ind <- panel.number()#
+						panel.segments(diff[ind], 0, diff[ind], 1, lty=2, col="grey30")#
+						panel.segments(diff[ind]/sqrt(2), 0, diff[ind]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30", x=list(axs="i")), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), between=list(x=0.2),#
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
 			#
-			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.groups = function(..., group.number) { #
+	  				 panel.xyplot(...)#
 	#
-			panel.levelplot(..., subscripts=subscripts)#
+	ind <- panel.number()#
+						panel.segments(diff[ind], 0, diff[ind], 1, lty=2, col="grey30")#
+						panel.segments(diff[ind]/sqrt(2), 0, diff[ind]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30", x=list(axs="i", at=seq(0.5, 1, by=0.1))), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), between=list(x=0.2),xlim=c(0.5, 1), #
+par.settings = list(	superpose.symbol = list(col=pastel.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
 			#
-			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)))#
-})#
+	},#
+	  		panel.groups = function(..., group.number) { #
+	  				 panel.xyplot(...)#
+	#
+	ind <- panel.number()#
+						panel.segments(diff[ind], 0, diff[ind], 1, lty=2, col="grey30")#
+						panel.segments(diff[ind]/sqrt(2), 0, diff[ind]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30", x=list(axs="i", at=seq(0.5, 1, by=0.1))), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p
 #
-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){#
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), between=list(x=0.2),xlim=c(0.5, 1), #
+par.settings = list(	superpose.symbol = list(col=dark.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
+			#
+	},#
+	  		panel.groups = function(..., group.number) { #
+	  				 panel.xyplot(...)#
 	#
-			panel.levelplot(..., subscripts=subscripts)#
+	ind <- panel.number()#
+						panel.segments(diff[ind], 0, diff[ind], 1, lty=2, col="grey30")#
+						panel.segments(diff[ind]/sqrt(2), 0, diff[ind]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30", x=list(axs="i", at=seq(0.5, 1, by=0.1))), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )
+p <- #
+xyplot(value~wavelength|pitch, data=cext.melt, groups=(theta), between=list(x=0.2),xlim=c(0.5, 1), #
+par.settings = list(	superpose.symbol = list(col=dark.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
 			#
-			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){#
+	},#
+	  		panel.groups = function(..., group.number) { #
+	  				 panel.xyplot(...)#
 	#
-	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)))#
-}#
+	ind <- panel.number()#
+						panel.segments(diff[ind], 0, diff[ind], 1, lty=2, col="grey30")#
+						panel.segments(diff[ind]/sqrt(2), 0, diff[ind]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30", x=list(axs="i", at=seq(0.5, 1, by=0.1))), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )#
+				#
+p#
 #
 #
+pdf("angleS.pdf")#
+print(p)#
+dev.off()
 #
+				#
+peak.ind <- which.max(df[, ncol(df)])#
+peaks <- sapply(df[peak.ind, 2:ncol(df)], "[" )#
+N <- (2*(1:15))#
+plot(N, peaks, ylab="peak intensity", xlab="Number of dipoles")
 #
-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)#
+pdf("angleS.pdf", width=6, height=5)#
+print(p)#
+dev.off()
+setwd('/Users/baptiste/Documents/R/photonics/pkg/cda/inst/examples')
 #
-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)#
+xyplot(value~wavelength, data=as.data.frame(both),subset=pitch==0.55,  groups=(theta), between=list(x=0.2),xlim=c(0.5, 1), #
+par.settings = list(	superpose.symbol = list(col=dark.col(length(unique(cext.melt$theta))), cex = 0.5), #
+superpose.line = list(col=dark.col(length(unique(cext.melt$theta))), lty = 1)#
+),#
+	type="l" , #
+	 		panel = function(...){ #
+	  			#
+	  	     	panel.grid(h = -1, v = -1, col = grey(0.9), lty = 1, lwd = 1)#
+	  			panel.superpose(..., lwd=2)#
 			#
-			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)))#
-})#
+	},#
+	  		panel.groups = function(..., group.number) { #
+	  				 panel.xyplot(...)#
+	#
+	ind <- panel.number()#
+						panel.segments(diff[ind], 0, diff[ind], 1, lty=2, col="grey30")#
+						panel.segments(diff[ind]/sqrt(2), 0, diff[ind]/sqrt(2), 0.1, lty=2, col="grey30")#
+	  		     }, 	xlab= expression("wavelength / "*mu*m),  ylab = expression(sigma[ext] / mu*m^2), #
+				      scales = list(col = "grey30", x=list(axs="i", at=seq(0.5, 1, by=0.1))), #
+					auto.key = list(space = "right", lines=F, points=T,text=paste(0:5),   cex.title = 1)#
+				  )
 #
-library(latticeExtra)#
-useOuterStrips(field.plot)
 #
-save("mieNearField.rda", field.plot, mdf)
+library(mie)#
+library(grid)#
+data(AuJChristy)#
+plot.material()$fine.sample->material#
 #
-save(file="mieNearField.rda", field.plot, mdf)
+radius <- 0.2#
+#
+#
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))
 setwd('/Users/baptiste/Documents/R/photonics/pkg/mie/inst')
 #
-quartz(type="png", file="MieNearFieldPlot2.png", height=7, width=7, bg="white", dpi=200)#
+radius <- 0.1#
 #
-print(useOuterStrips(field.plot))#
-dev.off()
+#
+with(material, matplot(lambda, Mie.spectra(lambda, epsilon, radius=radius, medium=1), t="l", ylab=~sigma))

Modified: pkg/scatteringR/man/Mie.Rd
===================================================================
--- pkg/scatteringR/man/Mie.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/scatteringR/man/Mie.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -50,8 +50,6 @@
 library(scatteringR)
 
 library(constants)
-clr()
-graphics.off()
 require(RColorBrewer)
 palette(brewer.pal(8,"Set1"))
 

Modified: pkg/scatteringR/man/alpha_s.Rd
===================================================================
--- pkg/scatteringR/man/alpha_s.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/scatteringR/man/alpha_s.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -27,8 +27,6 @@
 \examples{
 library(scatteringR)
 library(constants)
-clr()
-graphics.off()
 
 data(AuJChristy)
 plot.material(mat="AuJChristy",range=c(0.5,0.9),fine=300) -> gold

Modified: pkg/scatteringR/man/array.factor.Rd
===================================================================
--- pkg/scatteringR/man/array.factor.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/scatteringR/man/array.factor.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -32,7 +32,7 @@
 \examples{
 	library(baptMisc)
 	library(RColorBrewer)
-	clr()
+	
 	palette(brewer.pal(8, "Set1"))
 
 	t.extent <- c(10, 20, 100)

Modified: pkg/scatteringR/man/runAdda.spectrum.Rd
===================================================================
--- pkg/scatteringR/man/runAdda.spectrum.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/scatteringR/man/runAdda.spectrum.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -32,8 +32,6 @@
 \examples{
 library(scatteringR)
 library(constants)
-clr()
-
 data(luxpop)
 #### water ###
 plot.material("H2O_palik",range=c(0.4,2)) -> raw.data

Modified: pkg/spectroR/man/assignFiles.Rd
===================================================================
--- pkg/spectroR/man/assignFiles.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/spectroR/man/assignFiles.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -26,8 +26,8 @@
 \examples{
 library(spectroR)
 library(constants)
-clr()
 
+
 assignFiles()
 ls()
 

Modified: pkg/spectroR/man/exampleDF.Rd
===================================================================
--- pkg/spectroR/man/exampleDF.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/spectroR/man/exampleDF.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -27,9 +27,8 @@
 
 library(constants)
 library(spectroR)
-clr()
-graphics.off()
 
+
 require(RColorBrewer) 
 #####  graph customization #####
 palette(brewer.pal(8,"Set1")) 

Modified: pkg/spectroR/man/extractSpectra.Rd
===================================================================
--- pkg/spectroR/man/extractSpectra.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/spectroR/man/extractSpectra.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -35,9 +35,6 @@
 
 library(constants)
 library(spectroR)
-clr()
-graphics.off()
-
 require(RColorBrewer) 
 #####  graph customization #####
 palette(brewer.pal(8,"Set1")) 

Modified: pkg/spectroR/man/extractTspectra.Rd
===================================================================
--- pkg/spectroR/man/extractTspectra.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/spectroR/man/extractTspectra.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -34,8 +34,8 @@
 \examples{
 library(spectroR)
 library(constants)
-clr()
 
+
 assignFiles()
 ls()
 

Modified: pkg/spectroR/man/locatePeaks.Rd
===================================================================
--- pkg/spectroR/man/locatePeaks.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/spectroR/man/locatePeaks.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -31,8 +31,8 @@
 \examples{
 library(spectroR)
 library(constants)
-clr()
 
+
 assignFiles()
 ls()
 

Modified: pkg/spectroR/man/spectroR-package.Rd
===================================================================
--- pkg/spectroR/man/spectroR-package.Rd	2008-11-18 22:40:45 UTC (rev 336)
+++ pkg/spectroR/man/spectroR-package.Rd	2008-11-18 22:47:26 UTC (rev 337)
@@ -33,9 +33,8 @@
 
 library(constants)
 library(spectroR)
-clr()
-graphics.off()
 
+
 require(RColorBrewer) 
 #####  graph customization #####
 palette(brewer.pal(8,"Set1")) 



More information about the Photonics-commits mailing list