O R é um programa livre multiplataforma para análises estatísticas que pode ser baixado em seu site ou adicionado na lista de repositórios de máquinas linux. Suas possibilidades de aplicação em diversas áreas são praticamente ilimitadas.
Neste blog postarei o resultado de minha experiência em sua utilização nas áreas de dinâmica de populações de peixes, ciência pesqueira e ecologia.
As postagens deste blog se destina, além de mim mesmo, a iniciantes no R e alunos da minha área de atuação.
Aprendi muito em livros e nas listas de discussão R-help e a R_STAT, mas ainda tenho muito pela frente. Agradeço desde já qualquer contribuição.

quarta-feira, 31 de agosto de 2011

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

Um comentário:

  1. Montei esse algoritmo no JavaScript
    ano = ;
    mes = ;
    dia = ;
    a = 2*3.14159;
    b = ano-((12-mes)/10-((12-mes)/10)%1);
    c = (mes+8)%12+1;
    d = 365.25*(b+4712)-(365.25*(b+4712))%1;
    e = (30.6*c+0.5)-(30.6*c+0.5)%1;
    f = (b/100)+49-((b/100)+49)%1;
    g = f*0.75-(f*0.75)%1;
    h = g-38-(g-38)%1;
    if (d+e+dia+59 > 2299160) {
    i = d+e+dia+59-h;
    } else {
    i = d+e+dia+59;
    }
    j = (i-2451550.1)/29.530588853;
    if (j-(j-j%1) < 0) {
    k = j-(j-j%1)+1;
    } else {
    k = j-(j-j%1);
    }
    idadelua = k*29.530588853-(k*29.530588853)%1;
    if (idadelua == 0 || idadelua == 29) {
    faselua = "Nova";
    }
    else if (idadelua >= 0 && idadelua <= 6) {
    faselua = "Crescente Concava";
    }
    else if (idadelua == 7) {
    faselua = "4° Crescente";
    }
    else if (idadelua >= 8 && idadelua <= 13) {
    faselua = "Crescente Convexa";
    }
    else if (idadelua == 14) {
    faselua = "Cheia";
    }
    else if (idadelua >= 15 && idadelua <= 21) {
    faselua = "Minguante Convexa";
    }
    else if (idadelua == 22) {
    faselua = "4° Minguante";
    }
    else {
    faselua = "Minguante Concava";
    }

    ResponderExcluir