R/array.factor.R

array.factor <- function(N.part = 4, start = 0.6, end = 0.8, N.lambda = 2^7, pitch = 0.5*1.46){
	
	tlambda<-seq(from = start, to = end, length = N.lambda)
	tk <- 2*pi/tlambda
	
	rj <- t(expand.grid(-N.part:N.part, -N.part:N.part) * pitch)
	rj <- rj[1:2,  -((dim(rj)[2]-1)/2 + 1)] # remove rj=(0,0)

	moduli <- apply(rj, 2, function(x) sqrt(x %*% x))
	cis <- rj / matrix(rep(moduli, each = 2), nrow=2)
	cosines <- matrix(rep(cis[1, ], each = N.lambda), nrow = N.lambda)
	sines <- matrix(rep(cis[2, ], each = N.lambda), nrow = N.lambda)
	rij <- matrix(rep(moduli, each = N.lambda), nrow = N.lambda)

	exponentials <- outer(tk, moduli, function(x, y) { exp(1i*x*y) })
	kr <- outer(tk, moduli)

	A <- (1 - 1i * kr) * (3 * cosines^2 - 1) * exponentials / rij^3
	B <- tk^2 * sines^2 * exponentials / rij

	S <- rowSums(A+B)
	list(S=S, A=A, B=B, lambda=tlambda)
}

Try the array package in your browser

Any scripts or data that you put into this service are public.

array documentation built on May 2, 2019, 6:07 p.m.