Re: [Rd] Bug or not? (PR#8977)

From: Jean lobry <lobry_at_biomserv.univ-lyon1.fr>
Date: Wed 14 Jun 2006 - 17:55:34 GMT

I think you are computing your matrix to power N+1 instead of N. This code seems to work for me:

Rs <- function(lambda, theta, N){

   n0 <-  3.66
   n1 <-  2.4
   n2 <-  1.45
   nn <-  1.00

   lambda_ref <- 1000

   d1 <- lambda_ref / 4 / n1
   d2 <- lambda_ref / 4 / n2

   theta0 <- asin( nn / n0 * sin( theta ) )
   theta1 <- asin( nn / n1 * sin( theta ) )
   theta2 <- asin( nn / n2 * sin( theta ) )

   etha_s0 <- -n0 * cos( theta0 )
   etha_s1 <- n1 * cos( theta1 )

   etha_s2 <- n2 * cos( theta2 )
   etha_sn <- nn * cos( theta )

   delta1 <- 2 * pi / lambda * n1 * d1 * cos( theta1 )    delta2 <- 2 * pi / lambda * n2 * d2 * cos( theta2 )

   Ms1 <- matrix( c( cos( delta1 ) , 1i * etha_s1 * sin( delta1 ) , 1i / etha_s1 * sin( delta1 ) , cos( delta1 ) ), 2 , 2 )

   Ms2 <- matrix( c( cos( delta2 ) , 1i * etha_s2 * sin( delta2 ) , 1i / etha_s2 * sin( delta2 ) , cos( delta2 ) ), 2 , 2 )

   Mstmp <- Ms2 %*% Ms1
   Msr <- Mstmp
   for(i in 1:(N-1)) Msr <- Msr %*% Mstmp

   Rs <- ( abs( ( etha_sn * ( Msr[1, 1] - etha_s0 * Msr[1, 2] ) - ( Msr[2 , 1] - etha_s0 * Msr[2, 2] ) ) / ( etha_sn * ( Msr[1, 1] - etha_s0 * Msr[1, 2] ) + ( Msr[2, 1] - etha_s0 * Msr[2, 2] ) ) ) )^2

   return(Rs)
}

lambda_range <- 500:1500
s0 <- sapply(lambda_range, Rs, theta = 0, N = 5) s75 <- sapply(lambda_range, Rs, theta = 75*pi/180, N = 5)

ref <- read.table("Mathcad2.txt", header=FALSE, col.names=c('wl','a0', 'a75')) o0deg<-scan("o0deg.dat", "numeric", sep=" ", skip=5) o75deg<-scan("o75deg.dat", "numeric", sep=" ", skip=5) o0deg<-o0deg[-1]
o75deg<-o75deg[-1]

pdf(file="comparison.pdf", width=11.8, height=8.3)

par(mar = c(3.5,3.5,4,3.5))

plot(ref$wl,ref$a0,ylim=c(0,
1),type="l",col="1",lty=2,xlab="",ylab="",axes=FALSE, lwd = 1.5)

lines(lambda_range,s0,type="l",col="2", lty=2, lwd = 2)
lines(lambda_range,o0deg,type="l",col="3", lty=2)
lines(ref$wl,ref$a75,type="l",col="1",lty=3, lwd = 1.5)
lines(lambda_range,s75,col="2", lty=3, lwd = 2)
lines(lambda_range,o75deg,col="3", lty=3)

axis(1, at=seq(lambda_min,lambda_max,20)) axis(2)

mtext("wavelength [nm]", side=1, line=2)
mtext("reflection", side=2, line=2)
mtext(paste("bragg mirror, s-polarized: central wavelength ", 
lambda_ref," nm, ", N, " pairs of layers", seq=""), side=3, line=2.5, cex=1.2)
mtext(paste("n0=", n0 (1), "; n1=", n1(1),"; n2=", n2(1),"; n3=", nn(1)),side=3,line=1.5, cex=0.7)
legend("topleft", c("0A","75A"), lty=2:3) legend("topright", c("Reference","R", "Octave"), col=1:3, lty=1)

dev.off()

-- 
Jean R. Lobry            (lobry@biomserv.univ-lyon1.fr)
Laboratoire BBE-CNRS-UMR-5558, Univ. C. Bernard - LYON I,
43 Bd 11/11/1918, F-69622 VILLEURBANNE CEDEX, FRANCE
allo  : +33 472 43 12 87     fax    : +33 472 43 13 88
http://pbil.univ-lyon1.fr/members/lobry/

______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Thu Jun 15 05:07:27 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Wed 14 Jun 2006 - 22:22:26 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.