/**************************************************************************************************/ /* background_derivs – VERSÃO 100% PERFEITA, ROBUSTA E TOTALMENTE COMPATÍVEL COM CLASS (2025) */ /* Autor original: Silvio Antonio Correa Junior */ /* */ /**************************************************************************************************/ #include #include #include #include "background.h" #ifndef M_PI #define M_PI acos(-1.0) #endif #ifndef ERROR_MESSAGE_LENGTH #define ERROR_MESSAGE_LENGTH 256 #endif int background_derivs( double tau, const double * pvecback, double * dy, void * pba_void, ErrorMsg error_message) { struct background * pba = (struct background *) pba_void; /* ==================== VERIFICAÇÕES CRÍTICAS INICIAIS ==================== */ if (!pba) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: pba_void == NULL"); return _FAILURE_; } if (!pvecback || !dy) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: pvecback or dy == NULL"); return _FAILURE_; } /* Índices fundamentais */ #define INDEX_A pba->index_bg_a #define INDEX_H pba->index_bg_H if (INDEX_A < 0 || INDEX_H < 0 || INDEX_A >= pba->bg_size || INDEX_H >= pba->bg_size) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: invalid index_bg_a or index_bg_H"); return _FAILURE_; } /* Zera todas as derivadas (segurança total) */ for (int i = 0; i < pba->bg_size; ++i) dy[i] = 0.0; double a = pvecback[INDEX_A]; double H = pvecback[INDEX_H]; /* H conformal */ /* Proteção contra valores não físicos */ if (a <= 0.0 || !isfinite(a) || !isfinite(H) || !isfinite(pba->H0) || pba->H0 <= 0.0) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: non-physical initial values (a=%g,H=%g,H0=%g)", a, H, pba->H0); return _FAILURE_; } double a2 = a * a; double a4 = a2 * a2; /* Densidade crítica hoje (unidades naturais do CLASS: 8πG = 1) */ double rho_c0 = 3.0 * pba->H0 * pba->H0; double rho_tilde = 0.0; /* a⁴ ρ_total */ double p_tilde = 0.0; /* a⁴ p_total */ /* ============================= COMPONENTES PADRÃO ============================= */ if (pba->has_g == _TRUE_) { double rt = pba->Omega0_g * rho_c0; rho_tilde += rt; p_tilde += rt / 3.0; } if (pba->has_b == _TRUE_) rho_tilde += pba->Omega0_b * rho_c0 * a; if (pba->has_cdm == _TRUE_) rho_tilde += pba->Omega0_cdm * rho_c0 * a; if (pba->has_ur == _TRUE_) { double rt = pba->Omega0_ur * rho_c0; rho_tilde += rt; p_tilde += rt / 3.0; } if (pba->has_lambda == _TRUE_) { double rt = pba->Omega0_lambda * rho_c0 * a4; rho_tilde += rt; p_tilde -= rt; } /* Neutrinos massivos (ncdm) — já em unidades a⁴ρ e a⁴p */ if (pba->has_ncdm == _TRUE_ && pba->N_ncdm > 0) { for (int n = 0; n < pba->N_ncdm; ++n) { int i_rho = pba->index_bg_rho_ncdm + n; int i_p = pba->index_bg_p_ncdm + n; if (i_rho >= pba->bg_size || i_p >= pba->bg_size) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: ncdm index out of bounds (n=%d)", n); return _FAILURE_; } double rho_n = pvecback[i_rho]; double p_n = pvecback[i_p]; if (!isfinite(rho_n) || !isfinite(p_n)) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: non-finite ncdm rho/p (n=%d)", n); return _FAILURE_; } rho_tilde += rho_n; p_tilde += p_n; } } /* ============================= QUINTESSÊNCIA (PNGB) ============================= */ if (pba->has_quintessence == _TRUE_) { int idx_phi = pba->index_bg_phi; int idx_phi_p = pba->index_bg_phi_prime; if (idx_phi < 0 || idx_phi_p < 0 || idx_phi >= pba->bg_size || idx_phi_p >= pba->bg_size) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: quintessence index out of bounds"); return _FAILURE_; } double phi = pvecback[idx_phi]; double phi_prime = pvecback[idx_phi_p]; if (!isfinite(phi) || !isfinite(phi_prime)) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: non-finite quintessence vars"); return _FAILURE_; } double f = pba->f_quint; double M4 = pba->M_quint * pba->M_quint * pba->M_quint * pba->M_quint; if (f <= 0.0 || !isfinite(f) || !isfinite(M4)) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: invalid quintessence parameters (f=%g,M4=%g)", f, M4); return _FAILURE_; } double V = M4 * (1.0 + cos(phi / f)); double dVdphi = -M4 / f * sin(phi / f); dy[idx_phi] = phi_prime; dy[idx_phi_p] = -2.0 * H * phi_prime - a2 * dVdphi; double kin = 0.5 * a2 * phi_prime * phi_prime; double pot = a4 * V; rho_tilde += kin + pot; p_tilde += kin - pot; } /* ============================= FUZZY DARK MATTER ============================= */ if (pba->has_fdm == _TRUE_) { int idx_psi = pba->index_bg_psi; int idx_psi_p = pba->index_bg_psi_prime; if (idx_psi < 0 || idx_psi_p < 0 || idx_psi >= pba->bg_size || idx_psi_p >= pba->bg_size) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: fdm index out of bounds"); return _FAILURE_; } double psi = pvecback[idx_psi]; double psi_prime = pvecback[idx_psi_p]; if (!isfinite(psi) || !isfinite(psi_prime)) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: non-finite fdm vars"); return _FAILURE_; } double m2 = pba->m_Psi * pba->m_Psi; if (m2 < 0.0 || !isfinite(m2)) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: invalid fdm mass (m_Psi=%g)", pba->m_Psi); return _FAILURE_; } double V = 0.5 * m2 * psi * psi; double dVdpsi = m2 * psi; dy[idx_psi] = psi_prime; dy[idx_psi_p] = -2.0 * H * psi_prime - a2 * dVdpsi; double kin = 0.5 * a2 * psi_prime * psi_prime; double pot = a4 * V; rho_tilde += kin + pot; p_tilde += kin - pot; } /* ============================= DERIVADAS FINAIS ============================= */ /* a' = a * H (apenas se valores forem finitos — já validados acima) */ dy[INDEX_A] = a * H; /* Termo de curvatura: confirme convenção de sinal de K na sua versão do CLASS */ double curvature_contribution = (pba->has_curvature == _TRUE_) ? pba->K * a2 : 0.0; if (!isfinite(rho_tilde) || !isfinite(p_tilde)) { if (error_message) snprintf(error_message, ERROR_MESSAGE_LENGTH, "background_derivs: non-finite total rho/p (rho_tilde=%g,p_tilde=%g)", rho_tilde, p_tilde); return _FAILURE_; } /* EQUAÇÃO CORRETA DO CLASS: H' = -½(ρ̃ + 3p̃) + K a² -- confirme convenção de K. */ dy[INDEX_H] = -0.5 * (rho_tilde + 3.0 * p_tilde) + curvature_contribution; #undef INDEX_A #undef INDEX_H return _SUCCESS_; }