Dear list,
This is quite a specific question requiring the package orthopolynom.
This package provides a nice implementation of the Legendre
polynomials, however I need the associated Legendre polynomial which
can be readily expressed in terms of the mth order derivative of the
corresponding Legendre polynomial. (For the curious, I'm trying to
calculate spherical harmonics [*]).
Because legendre.polynomials(l) returns a list of Legendre polynomials
of degree 0 to l, I'd like to make use of the whole list of them at a
time rather than wasting this information. For a given degree "l" I
therefore have a list of l+1 polynomials. For each of these I want to
compute l+1 derivatives, from m= 0 to m=l. The last step is to
evaluate all of these polynomials with a vector argument and return a
list of data.frames.
I've come up with the following hack but it's really ugly,
require(orthopolynom)
md <- function(.p, m=2){
test <- list()
if(.p==0) pl.list <- rep(as.polylist(.p), m+1) else {
pl.list <- as.polylist(.p)
for(n in seq(1, m+1)){
pl.list[[n+1]] <- deriv(pl.list[[n]])
}
}
rev(pl.list) # ascending order
}
l <- 3 # example
theta <- seq(0, pi, length= 10) # the variable to evaluate the
polynomials at
Pl <- as.polylist(legendre.polynomials(l))
Plm <- lapply(seq_along(Pl), function(ind) md(Pl[[ind]], ind-1))
Plm.theta <- lapply(seq_along(Plm), function(ind) # treat each order l
sapply(seq_along(Plm[[ind]]), function(ind2) # treat each order m
(-1)^ind2 *(1-cos(theta)^2)^(ind2/2) * as.function(Plm[[ind]]
[[ind2]])( cos(theta)) )) # evaluate the expression in theta
I tried (unsuccessfully) to get inspiration from Recall() but since I
want to store the intermediate derivatives it doesn't seem very
suitable anyway.
Any advice is welcome!
[*] http://en.wikipedia.org/wiki/Spherical_harmonic
_____________________________
Baptiste AuguiƩ
School of Physics
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK
Phone: +44 1392 264187
http://newton.ex.ac.uk/research/emag
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.