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 ([EMAIL PROTECTED])
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
