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')
