[R] code optimization tips

[baptiste Auguié]

Hi,

Being new to R I'm asking for some advice on how to optimize the  
performance of the following piece of code:


> alpha_c <- function(lambda=600e-9,alpha_s=1e-14,N=400,spacing=1e-7){
>
> k<-2*pi/lambda
> ri<-c(0,0) # particle at the origin
> x<-c(-N:N)
> positions <- function(N) {
>    reps <- 2*N+1
>    matrix(c(rep(-N:N, each = reps), rep(-N:N, times = reps)),
>           nrow = 2, byrow = TRUE)
> }
> rj<-positions(N)*spacing # all positions in the 2N x 2N array
> rj<-rj[1:2,-((dim(rj)[2]-1)/2+1)] # remove rj=(0,0)
>
> mod<-function(x){sqrt(x[1]^2+x[2]^2)} # modulus
>
> sij <-function(rj){
> rij=mod(rj-ri)
> cos_ij=rj[1]/rij
> sin_ij=rj[2]/rij
>
> A<-(1-1i*k*rij)*(3*cos_ij^2-1)*exp(1i*k*rij)/(rij^3)
> B<-k^2*sin_ij^2*exp(1i*k*rij)/rij
>
> sij<-A+B
> }
>
> s_ij<-apply(rj,2,sij)
> S<-sum(s_ij)
> alpha_s/(1-alpha_s*S)
> }
> alpha_c()


This function is to be called for a few tens of values of lambda in a  
'for' loop, and possibly a couple of different N and spacing (their  
magnitude is typically around the default one).

This can be a bit slow ––– not that I would expect otherwise --- and  
I wonder if there is something I could do to optimize it (vectorize  
with respect to the lambda parameter?, change the units of the  
problem to deal with numbers closer to unity?,...)

Best regards,

baptiste