setwd("/Users/amaral/Tracker") New_Horizons<-read.table("New_Horizons.txt", header=TRUE) Plutão<-read.table("Plutão.txt", header=TRUE) MU69<-read.table("2014 MU69.txt", header=TRUE) #plot(New_Horizons$x,New_Horizons$y, type="n", xlab="abscissa", ylab="ordenada") #grid() #lines(New_Horizons$x,New_Horizons$y) #lines(Plutão$x,Plutão$y) #lines(MU69$x,MU69$y) #abline(v=0, lty=3) #abline(h=0, lty=3) #box(lwd=3) ## 1) Limpar dados (tirar linhas com NA em x ou y) Plutão_limpo <- subset(Plutão, !is.na(x) & !is.na(y)) ## 2) Coordenadas polares no referencial com o Sol na origem Plutão_limpo$r <- sqrt(Plutão_limpo$x^2 + Plutão_limpo$y^2) Plutão_limpo$theta <- atan2(Plutão_limpo$y, Plutão_limpo$x) # em radianos ## 3) Ajustar 1/r = A + B cos(theta) + C sin(theta) Plutão_limpo$inv_r <- 1 / Plutão_limpo$r Plutão_limpo$c <- cos(Plutão_limpo$theta) Plutão_limpo$s <- sin(Plutão_limpo$theta) fit <- lm(inv_r ~ c + s, data = Plutão_limpo) A <- coef(fit)[["(Intercept)"]] B <- coef(fit)[["c"]] C <- coef(fit)[["s"]] ## 4) Recuperar parâmetros orbitais geométricos p <- 1 / A e <- sqrt(B^2 + C^2) / A omega <- atan2(C, B) # argumento do periastro (rad) a <- p / (1 - e^2) # semieixo maior da elipse cat("p =", p, "\n") cat("e =", e, "\n") cat("ω =", omega, "rad\n") cat("a =", a, "\n") ## 5) Gerar a órbita completa (0 a 2π, por exemplo) theta_grid <- seq(0, 2*pi, length.out = 1000) r_grid <- p / (1 + e * cos(theta_grid - omega)) x_grid <- r_grid * cos(theta_grid) y_grid <- r_grid * sin(theta_grid) orbit <- data.frame(theta = theta_grid, x = x_grid, y = y_grid) ## --- NOVO: calcular limites para incluir tudo --- x_all <- c(Plutão_limpo$x, orbit$x) y_all <- c(Plutão_limpo$y, orbit$y) ## 6) Desenhar o arco observado + órbita extrapolada plot(Plutão_limpo$x, Plutão_limpo$y, asp = 1, pch = 16, xlab = "x", ylab = "y", xlim = range(x_all), ylim = range(y_all), main = "órbitas keplerianas ajustadas", type="n") #lines(orbit$x, orbit$y, col = "red", lwd = 2) #points(0, 0, pch = 3, col = "blue", cex = 1.2) # foco (Sol) #legend("topleft", # legend = c("dados (arco)", "órbita ajustada", "Sol"), # col = c("black", "red", "blue"), # pch = c(16, NA, 3), # lty = c(NA, 1, NA), # bty = "n") ## 0) MU69: data.frame com colunas t, x, y ## MU69 <- read.table("MU69.txt", header = TRUE) # se estiveres a ler de ficheiro ## 1) Limpar dados (tirar linhas com NA em x ou y) MU69_limpo <- subset(MU69, !is.na(x) & !is.na(y)) ## 2) Coordenadas polares (Sol na origem) MU69_limpo$r <- sqrt(MU69_limpo$x^2 + MU69_limpo$y^2) MU69_limpo$theta <- atan2(MU69_limpo$y, MU69_limpo$x) # radianos ## 3) Ajustar 1/r = A + B cos(theta) + C sin(theta) MU69_limpo$inv_r <- 1 / MU69_limpo$r MU69_limpo$c <- cos(MU69_limpo$theta) MU69_limpo$s <- sin(MU69_limpo$theta) fit_MU69 <- lm(inv_r ~ c + s, data = MU69_limpo) coef(fit_MU69) A_MU <- coef(fit_MU69)[["(Intercept)"]] B_MU <- coef(fit_MU69)[["c"]] C_MU <- coef(fit_MU69)[["s"]] ## 4) Parâmetros orbitais geométricos p_MU <- 1 / A_MU e_MU <- sqrt(B_MU^2 + C_MU^2) / A_MU omega_MU <- atan2(C_MU, B_MU) # argumento do periastro em rad a_MU <- p_MU / (1 - e_MU^2) # semieixo maior cat("MU69:\n") cat("p =", p_MU, "\n") cat("e =", e_MU, "\n") cat("ω =", omega_MU, "rad\n") cat("a =", a_MU, "\n") ## (Se quiseres confirmar: para estes dados deve dar algo perto de ## p ≈ 3.26e2, e ≈ 0.47, a ≈ 4.18e2, em unidades dos teus x,y.) ## 5) Gerar a órbita completa (0 a 2π) theta_grid_MU <- seq(0, 2*pi, length.out = 1000) r_grid_MU <- p_MU / (1 + e_MU * cos(theta_grid_MU - omega_MU)) x_grid_MU <- r_grid_MU * cos(theta_grid_MU) y_grid_MU <- r_grid_MU * sin(theta_grid_MU) orbita_MU <- data.frame(theta = theta_grid_MU, x = x_grid_MU, y = y_grid_MU) ## 6) Definir limites do gráfico com base em dados + órbita x_all_MU <- c(MU69_limpo$x, orbita_MU$x) y_all_MU <- c(MU69_limpo$y, orbita_MU$y) ## 7) Desenhar arco observado + órbita completa #plot(MU69_limpo$x, MU69_limpo$y, # asp = 1, pch = 16, # xlab = "x", ylab = "y", # xlim = range(x_all_MU), # ylim = range(y_all_MU), # main = "MU69: arco observado e órbita kepleriana ajustada") #lines(orbita_MU$x, orbita_MU$y, col = "red", lwd = 2) #points(0, 0, pch = 3, col = "blue", cex = 1.2) # Sol (foco) #legend("topleft", # legend = c("dados (arco)", "órbita ajustada", "Sol"), # col = c("black", "red", "blue"), # pch = c(16, NA, 3), # lty = c(NA, 1, NA), # bty = "n") ## Semieixo maior de Plutão em UA (valor astronómico) a_Pl_au <- 39.48 # podes ajustar para o valor que preferires ## Escala usada nos teus dados (unidades do gráfico por UA) escala_ua <- a / a_Pl_au # 'a' é o semieixo maior de Plutão em unidades do teu gráfico escala_ua ## Parâmetros de Neptuno em UA a_Nep_au <- 30.07 e_Nep <- 0.0086 ## Semieixo maior de Neptuno nas tuas unidades a_Nep <- a_Nep_au * escala_ua ## Semieixo menor b_Nep <- a_Nep * sqrt(1 - e_Nep^2) ## Parâmetro: anomalia excêntrica E de 0 a 2π E_Nep <- seq(0, 2*pi, length.out = 1000) ## Elipse kepleriana com o Sol num foco (origem) x_Nep <- a_Nep * cos(E_Nep) - a_Nep * e_Nep # desloca o centro para pôr o foco na origem y_Nep <- b_Nep * sin(E_Nep) Neptuno_orbita <- data.frame(x = x_Nep, y = y_Nep) ## Pontos de intercepção entre New Horizons e a órbita de Neptuno ## Elipse de Neptuno já definida a <- a_Nep b <- b_Nep e <- e_Nep ## Coeficientes da recta ajustada à New Horizons ajuste <- lm(New_Horizons$y[50:197] ~ New_Horizons$x[50:197]) c_int <- coef(ajuste)[1] # intercepto m_slope<- coef(ajuste)[2] # declive m <- m_slope c <- c_int ## Coeficientes do polinómio Ax^2 + Bx + C = 0 A <- a^2 * m^2 + b^2 B <- 2 * a^2 * c * m + 2 * a * b^2 * e C <- a^2 * b^2 * e^2 - a^2 * b^2 + a^2 * c^2 disc <- B^2 - 4 * A * C if (disc < 0) { stop("Sem intersecção real entre a recta e a órbita elíptica de Neptuno.") } x1 <- (-B + sqrt(disc)) / (2 * A) x2 <- (-B - sqrt(disc)) / (2 * A) y1 <- m * x1 + c y2 <- m * x2 + c pontos_intersecção <- rbind( c(x = x1, y = y1), c(x = x2, y = y2) ) pontos_intersecção pontos_intersecção <- data.frame(x=pontos_intersecção[,1], y=pontos_intersecção[,2]) ## Gráfico completo com pontos de intercepção grid() abline(v=0, lty=3) abline(h=0, lty=3) abline(lm(New_Horizons$y[50:197]~New_Horizons$x[50:197]), lty=3, lwd=0.6) lines(orbit$x, orbit$y, col = "cyan", lwd = 2) lines(orbita_MU$x, orbita_MU$y, col = "green", lwd = 2) lines(Neptuno_orbita$x, Neptuno_orbita$y, col = "violet", lwd = 2) points(0, 0, pch = 3, col = "olivedrab1", cex = 1.2) # Sol (foco) lines(Plutão$x,Plutão$y) lines(MU69$x,MU69$y) lines(New_Horizons$x,New_Horizons$y, col="red") points(c(0,0), pch=21, bg="yellow", cex=0.4) points(pontos_intersecção, pch=21, bg="red", cex=0.4) box(lwd=3) ## Fazer a GIF animada com Neptuno das efemérides (Swiss Ephemeris) library(ggplot2) library(gganimate) library(dplyr) library(gifski) library(swephR) ## 0) Inicializar a Swiss Ephemeris swe_set_ephe_path(NULL) # usa caminho default (Moshier ou ficheiros swephRdata) data(SE) # garante que o objecto SE está carregado ## ATENÇÃO: usar só flags que existem em swephR ## heliocêntrico + Moshier iflag <- SE$FLG_MOSEPH + SE$FLG_HELCTR ## Função utilitária: converter t do Tracker em dia juliano ## Aqui estou a supor que t está em anos a partir de 2006-01-19. ## Se a escala de t for outra, ajusta 'escala_t_para_dias'. t_to_jd <- function(t, year0 = 2006, month0 = 1, day0 = 19, hour0 = 0, escala_t_para_dias = 365.25) { jd0 <- swe_julday(year0, month0, day0, hour0, SE$GREG_CAL) jd0 + t * escala_t_para_dias } ## 1) Preparar trajectórias com etiqueta do corpo (dados do Tracker) Plutao_mov <- Plutão[, c("t","x","y")] Plutao_mov$corpo <- "Plutão" MU69_mov <- MU69[, c("t","x","y")] MU69_mov$corpo <- "Arrokoth" NH_mov <- New_Horizons[, c("t","x","y")] NH_mov$corpo <- "New Horizons" traj_all <- bind_rows(Plutao_mov, MU69_mov, NH_mov) |> filter(!is.na(x), !is.na(y), !is.na(t)) |> arrange(t) ## 2) Intervalo comum de t a todos os três (sincronização real) t_min <- max(tapply(traj_all$t, traj_all$corpo, min)) t_max <- min(tapply(traj_all$t, traj_all$corpo, max)) traj_all <- traj_all |> filter(t >= t_min, t <= t_max) ## 3) Criar um corpo extra para o Sol, com t em todo o intervalo t_seq <- sort(unique(traj_all$t)) if (length(t_seq) == 0) { stop("Não há valores de t comuns aos três corpos depois do filtro.") } Sol_mov <- data.frame( t = t_seq, x = 0, y = 0, corpo = "Sol" ) ## 3b) Criar trajectória de Neptuno a partir das efemérides (Swiss Ephemeris) ## mapear t -> dia juliano (ajusta escala_t_para_dias se for necessário) jd_seq <- t_to_jd(t_seq, year0 = 2006, month0 = 1, day0 = 19, hour0 = 0, escala_t_para_dias = 365.25) ## pré-alocar vectores de Neptuno nas tuas unidades nep_x <- numeric(length(jd_seq)) nep_y <- numeric(length(jd_seq)) for (i in seq_along(jd_seq)) { jd <- jd_seq[i] res <- swe_calc_ut(jd, SE$NEPTUNE, iflag) if (res$return < 0) { stop(paste("Erro em swe_calc_ut para Neptuno (i =", i, ", jd =", jd, "):", res$serr)) } lon_deg <- res$xx[1] # longitude eclíptica heliocêntrica r_au <- res$xx[3] # distância em UA lon_rad <- lon_deg * pi/180 x_au <- r_au * cos(lon_rad) y_au <- r_au * sin(lon_rad) nep_x[i] <- x_au * escala_ua nep_y[i] <- y_au * escala_ua } Neptuno_mov <- data.frame( t = t_seq, x = nep_x, y = nep_y, corpo = "Neptuno" ) ## 4) Juntar tudo: dados Tracker + Sol + Neptuno traj_all2 <- bind_rows(traj_all, Sol_mov, Neptuno_mov) ## 5) Órbitas teóricas de fundo bg_orbits <- bind_rows( cbind(orbit[, c("x","y")], tipo = "Órbita Plutão"), cbind(orbita_MU[, c("x","y")], tipo = "Órbita Arrokoth"), cbind(Neptuno_orbita[, c("x","y")], tipo = "Órbita Neptuno") ) ## 6) Limites incluindo tudo (dados + órbitas) x_lim <- range(c(traj_all2$x, bg_orbits$x), na.rm = TRUE) y_lim <- range(c(traj_all2$y, bg_orbits$y), na.rm = TRUE) ## 7) Construir a GIF p <- ggplot() + # órbitas teóricas de fundo geom_path(data = bg_orbits, aes(x, y, group = tipo, colour = tipo), alpha = 0.3, linewidth = 0.5) + # trajectórias dos corpos em movimento (Plutão, Arrokoth, New Horizons, Neptuno) geom_path(data = traj_all2 |> filter(corpo != "Sol"), aes(x, y, group = corpo, colour = corpo), alpha = 0.6) + # pontos instantâneos (inclui todos os corpos, Sol incluído) geom_point(data = traj_all2, aes(x, y, colour = corpo), size = 2) + scale_colour_manual( values = c( # corpos "Plutão" = "deepskyblue3", "Arrokoth" = "darkgreen", "Neptuno" = "purple", "New Horizons" = "red2", "Sol" = "gold", # órbitas de fundo "Órbita Plutão" = "deepskyblue3", "Órbita Arrokoth" = "darkgreen", "Órbita Neptuno" = "purple" ) ) + coord_equal(xlim = x_lim, ylim = y_lim) + theme_minimal() + theme(legend.position = "bottom") + labs(x = "x", y = "y", title = "t = {round(frame_along, 3)}") + transition_reveal(t) ## 8) Gerar e gravar GIF anim_save("NewHorizons_planetas_com_Neptuno_ephem.gif", animation = p, fps = 20, duration = 12, width = 800, height = 800, renderer = gifski_renderer()) ## Para converter a GIF para um MPEG4 #ffmpeg -i ~/Tracker/NewHorizons_planetas_com_Neptuno_ephem.gif -movflags faststart -pix_fmt yuv420p -vf "scale=trunc(iw/2)*2:trunc(ih/2)*2" New_Horizons2a.mp4