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
Rs <- function (lambda, theta, N) {
n0 = 3.66
n1 = 2.4
n2 = 1.45
nn = 1.00
lambda_ref = 1000
lambda_min = 500
lambda_max = 1500
increment = 5
ymin = 0
ymax = 1
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,byrow=T)
Ms2 = matrix(c(cos( delta2 ) , 1i / etha_s2 * sin( delta2 ) , 1i * etha_s2 *
sin( delta2) , cos( delta2 )) ,2,2,byrow=T)
dum = Ms2 %*% Ms1
Msr <- diag(2)
for (i in 1:N) Msr <- dum %*% Msr
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
Rs
}
##r0deg <- r75deg <- numeric(1001)
##for (i in 500:1500) r0deg[i]=Rs(i,0,5)
##for (i in 500:1500) r75deg[i]=Rs(i,75*pi/180,5)
##
##oct0 <- read.table('o0deg,dat')
##oct75 <- read.table('o75deg,dat')
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
