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

[Martin Hoffmann]

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"

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 254 bytes
Desc: OpenPGP digital signature
Url : https://stat.ethz.ch/pipermail/r-devel/attachments/20060614/4ad917a8/attachment-0001.bin