This is an OpenPGP/MIME signed message (RFC 2440 and 3156) --------------enig807B2312A20EAF60129FDDFA Content-Type: text/plain; charset=ISO-8859-1 Content-Transfer-Encoding: quoted-printable Hi, I am writing this email, because I am not sure if the issue I have discovered is a bug or not. For a few days I have been fiddling around with a small program that calculates the reflectance of multilayer dielectric mirrors (eg. bragg mirrors). It does not really matter what I did, what matters is that I could not do that using R. I had some sample data which I used to test if my R program works correctly (if it fits the data it works). This sample data was for two different incident angles with respect to the normal of the hypothetical bragg mirror (0=B0 and 70=B0). The plots for 0=B0 matched perfectly but those for 70=B0 were quite a bit off. After trying out all sorts of things for several days, I tried the same algorithm in octave (matlab clone, but open source). There, I got a perfect match for both 0=B0 and 70=B0. I almost could not believe it, R must have calculated wrong. The R version I use at home at the moment is 2.2.1 on gentoo linux. I also tried version 2.1.0 at my university (debian stable), both had the same result. I uploaded everything to reproduce this to http://n.ethz.ch/student/homartin/r-vs-octave/ There, you can also find the .pdf file of the comparison of the reference, the R output and the one from octave. The calculations are quite complex, nevertheless, I would be very happy to read your reply concerning this problem. Thanks in advance, best Regards, Martin --=20 Better to use medicines at the outset than at the last moment. --------------enig807B2312A20EAF60129FDDFA Content-Type: application/pgp-signature; name="signature.asc" Content-Description: OpenPGP digital signature Content-Disposition: attachment; filename="signature.asc" -----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.2.2 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org iD8DBQFEj2Yw899ysxFmsp8RArgrAJwLYHCy0vf5m6/iAQwtt+/+pOFPWACg/Q13 6xkVN0cpZjeKAeGWj7qKdOM=hv4P -----END PGP SIGNATURE----- --------------enig807B2312A20EAF60129FDDFA--
> Hi, > > I am writing this email, because I am not sure if the issue I have > discovered is a bug or not. > > For a few days I have been fiddling around with a small program that > calculates the reflectance of multilayer dielectric mirrors (eg. bragg > mirrors). It does not really matter what I did, what matters is that I > could not do that using R. > I had some sample data which I used to test if my R program works > correctly (if it fits the data it works). This sample data was for two > different incident angles with respect to the normal of the hypothetical > bragg mirror (0=B0 and 70=B0). The plots for 0=B0 matched perfectly but t> hose > for 70=B0 were quite a bit off. After trying out all sorts of things for > several days, I tried the same algorithm in octave (matlab clone, but > open source). There, I got a perfect match for both 0=B0 and 70=B0. I alm> ost > could not believe it, R must have calculated wrong. > > The R version I use at home at the moment is 2.2.1 on gentoo linux. I > also tried version 2.1.0 at my university (debian stable), both had the > same result. > > I uploaded everything to reproduce this to > > http://n.ethz.ch/student/homartin/r-vs-octave/ > > There, you can also find the .pdf file of the comparison of the > reference, the R output and the one from octave. > > The calculations are quite complex, nevertheless, I would be very happy > to read your reply concerning this problem. > > Thanks in advance, best Regards, > Martin > >I get only the r-devel digest, so I could'nt respond the direct way (this message will not show up in the correct thread, probably), but: there is nothing wrong with the R results AFAICS: I translated your octave script one-to-one to R (see attachment). this produces (within floating point prec. (double)) the same results. I propose you go through this and your own R-script (which honestly was unreadable to me :-)) with browser() or debug() and compare results for Ms1, Ms2, Msr on the way to localise the differences in computation. if you uncomment the last few lines, source("Rs.R") will give you vectors which you can directly compare to the octave results (note that the leading 499 zeros are still there, since I mimicked your octave script) so: really no bug, right? joerg -------------- next part -------------- An embedded and charset-unspecified text was scrubbed... Name: Rs.R Url: https://stat.ethz.ch/pipermail/r-devel/attachments/20060614/6cc27ce0/attachment.pl
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/