Cálculo das fases da lua no R

A influência do ciclo lunar sobre os organismos marinhos é grande e, como consequência, sobre a atividade pesqueira. Muitas vezes ao analisar dados queremos relacionar a abundância das espécies ou a CPUE das pescarias à fase lunar.

Por esta razão, necessitei de uma função em R que indicasse a fase da lua a partir de uma data. Achei em http://www.paulsadowski.com/wsh/moonphase.htm um código em Visual Basic para este cálculo e o portei para R.
As fases da lua ficaram como:
nova -> crescente concava -> quarto crescente -> crescente convexa -> cheia -> minguante convexa -> quarto minguante -> minguante concava -> nova.

Verifiquei as respostas dadas pela função com o programa LunaBar e em um site Islâmico (http://islam.com.pt/), na seção Fases da Lua. Descobri que o calendário Islâmico é baseado no ciclo lunar e por isso eles oferecem a informação de maneira exata. Os resultados que obtive com a função foram bem próximos aos calculados no site e no LunaBar.

Acredito que a função pode ser utilizada sem inconvenientes, no entanto não posso garantir a exatidão dos cálculos.
A função pode ser facilmente alterada caso se deseje apenas a idade ou a fase da lua.
Achei interessante a utilização da função recode, do pacote car, que deve estar carregado.

require(car)

lua(2011,8,31)

Idade da lua (dias): 2
Fase da lua: crescente concava
Distância (raio): 56.67
Latitude Elíptica (graus): -4.8
Longitude Elíptica (graus): 190.11 

lua<- function(Y,M,D) {
# http://www.paulsadowski.com/wsh/moonphase.htm

# necessita do pacote cars
P2<-2*3.14159
YY<-Y-as.integer((12-M)/10)
MM=M+9
if (MM>=12) MM<-MM-12
K1=as.integer(365.25*(YY+4712))
K2=as.integer(30.6*MM+.5)
K3=as.integer(as.integer((YY/100)+49)*.75)-38

# J é a data às 12h UT do dia em questão

J<-K1+K2+D+59
if(J>2299160) J<-J-K3
# Calcula a fase sinódica
V<-(J-2451550.1)/29.530588853
V<-V-as.integer(V)
if(V<0) V<-V+1
IP<-V
# Idade da Lua em dias
AG<-IP*29.53
IP<-IP*P2 # Converte fase em radianos
# Calcula a distância
V<-(J-2451562.2)/27.55454988
V<-V-as.integer(V)
if(V<0) V<-V+1
DP<-V
DP<-DP*P2 # Converte em radianos

DI<-60.4-3.3*cos(DP)-.6*cos(2*IP-DP)-.5*cos(2*IP)
# Calcula a Latitude
V<-(J-2451565.2)/27.212220817
V<-V-as.integer(V)
if(V<0) V<-V+1
NP<-V
NP<-NP*P2 # Converte em radianos
LA<-5.1*sin(NP)
# Calcula a Longitude
V<-(J-2451555.8)/27.321582241
# Normaliza valores para o intervalo de 0 a 1
V<-V-as.integer(V)
if(V<0) V<-V+1
RP<-V
LO<-360*RP+6.3*sin(DP)+1.3*sin(2*IP-DP)+.7*sin(2*IP)
# fases em inglês
# http://home.hiwaay.net/~krcool/Astro/moon/moonphase/
Phase<-c("nova","crescente concava","quarto crescente","crescente convexa","cheia","minguante convexa","quarto minguante","minguante concava")
ThisPhase<-recode(as.integer(AG),"c(0,29)=1;c(1,2,3,4,5,6)=2;c(7)=3;c(8,9,10,11,12,13)=4;c(14)=5;c(15,16,17,18,19,20,21)=6;c(22)=7;c(23,24,25,26,27,28)=8")
message("Idade da lua (dias): ", as.integer(AG))
message("Fase da lua: ",Phase[ThisPhase])
message("Distância (raio): ",round(DI,2))
message("Latitude Elíptica (graus): ",round(LA,2))
message("Longitude Elíptica (graus): ",round(LO,2))
}