Joerg van den Hoff wrote: > >> 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
NO BUG! I found out what was wrong and it was all my fault!
The etha functions were the problem. They should be etha_sX(thetaX)
where thetaX is a function of theta. I made the etha_sX functions that
way and defined the etha_sX functions to be functions of theta, NOT thetaX.
But when I used these I called them with thetaX instead of only theta.
I would like to apologize for the "bug", I hope this didn't cost you too
much time. Rest assured I feel quite embarrassed now :-?.
I can imagine the script was hard to read, I am new to R and sometimes
it is still hard for me to get things done at all.
Thanks for your time and the great product of course ;-).
Best regards,
Martin
>
>
> ------------------------------------------------------------------------
>
> 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')
--
Sometimes I wonder if I'm in my right mind. Then it passes off and I'm
as intelligent as ever.
-- Samuel Beckett, "Endgame"
signature.asc
Description: OpenPGP digital signature
______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
