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

[Jean lobry]

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 at 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/