diff --git a/symcelmech/README.md b/symcelmech/README.md new file mode 100644 index 0000000..5336ba0 --- /dev/null +++ b/symcelmech/README.md @@ -0,0 +1,218 @@ +# symcelmech + +An open-source Maxima package for analytical celestial mechanics based on Lie transformations. + +`symcelmech` automates the construction of high-order perturbation theories for Hamiltonian systems using the Hori-Deprit method (Lie series). It provides a complete symbolic pipeline — from potential derivation through orbital averaging to generating function computation — supporting both closed-form (exact in eccentricity) and series expansion modes. + +## Features + +- Automated Hori-Deprit/Mersman recursive Lie triangle up to arbitrary order +- Dual-track integration: closed-form (true anomaly domain) and Poisson series (mean anomaly domain) +- Perturbing potential derivation: zonal harmonics (J_n), tesseral/sectorial harmonics (C_nm, S_nm), third-body perturbations, and solar radiation pressure (SRP) +- Hansen coefficient engine for exact orbital averaging +- Memoized Poisson bracket computation with antisymmetry exploitation +- Mean-to-osculating element transformation via Lie series corrections +- Deprit's elimination of the parallax +- Lagrange, Gauss, and Hamilton planetary equations (classical and non-singular equinoctial) +- Delaunay-Poincare canonical transformations +- Elliptic expansion series (equation of the center, r/a, cos f, sin f, cos E) +- Physical constants database (Earth, Moon, Mercury, Apophis, Didymos) + +## Requirements + +- [Maxima](https://maxima.sourceforge.io/) (tested with 5.47+) +- Maxima packages: `orthopoly`, `vect` (bundled with standard Maxima distributions) + +## Installation + +### Option 1: Add to Maxima search path (recommended) + +Clone the repository: + +```bash +git clone https://github.com//symcelmech.git +``` + +Add the package directory to your `maxima-init.mac` (created automatically in `~/.maxima/` on Linux/macOS or `%USERPROFILE%/maxima/` on Windows): + +```maxima +/* Add this line to your maxima-init.mac */ +file_search_maxima: cons("/path/to/symcelmech/", file_search_maxima)$ +``` + +Then, from any Maxima session: + +```maxima +load("symcelmech")$ +``` + +### Option 2: Install via maxima-packages + +`symcelmech` follows the [maxima-packages](https://github.com/maxima-project-on-github/maxima-packages) convention and can be installed from the community repository (pending acceptance): + +```maxima +/* After maxima-packages is set up */ +load("symcelmech/symcelmech.mac")$ +``` + +### Option 3: Direct load + +If you prefer not to modify `maxima-init.mac`, load the entry point with the full path: + +```maxima +load("/path/to/symcelmech/symcelmech.mac")$ +``` + +### Option 4: Load individual modules + +```maxima +load("/path/to/symcelmech/derivepotential.mac")$ +load("/path/to/symcelmech/planetaryequations.mac")$ +``` + +Internal dependencies are resolved automatically. + +## Running Tests + +The package includes a regression test suite compatible with Maxima's `run_testsuite`: + +```maxima +batch("rtest_symcelmech.mac", test)$ +``` + +## Package Structure + +| Module | Description | +|---|---| +| `lietransformations.mac` | Core Hori-Deprit solver (`solve_hori_deprit_recurrence`) and Lie series corrections (`apply_lie_corrections`) | +| `analyticalfunctions.mac` | Poisson bracket engines (series and closed-form), generating function quadrature, trigonometric expansion tools, elimination of the parallax | +| `derivepotential.mac` | Perturbing potential derivation: zonal, tesseral, third-body, and SRP | +| `average.mac` | Secular/periodic separation, Hansen coefficient computation, orbital averaging | +| `celmechseries.mac` | Elliptic expansion series and eccentricity truncation tools | +| `algebraicutils.mac` | Memoization, expression sanitization, coordinate substitution rules | +| `planetaryequations.mac` | Lagrange, Gauss, and Hamilton equations; Delaunay-Poincare transformations | +| `astrodata.mac` | Physical and orbital constants for supported celestial bodies | + +### Dependency Graph + +``` +lietransformations.mac +├── analyticalfunctions.mac +├── average.mac +│ └── analyticalfunctions.mac +├── celmechseries.mac +└── algebraicutils.mac + └── analyticalfunctions.mac + +derivepotential.mac +├── analyticalfunctions.mac +└── celmechseries.mac + +planetaryequations.mac (standalone) +astrodata.mac (standalone) +``` + +## Quick Start + +### Example: J2 Frozen Orbit Theory (Closed-Form) + +```maxima +/* Load the pipeline */ +load("lietransformations.mac")$ +load("derivepotential.mac")$ +load("astrodata.mac")$ + +/* Set up canonical dependencies for closed-form mode */ +setup_canonical_dependencies()$ + +/* Derive the J2 zonal potential */ +U2: derive_zonal_potential(2, mu, r, f, a, e, i, g, R, true)$ + +/* Solve the Hori-Deprit recurrence to second order */ +vars: [l, g, h, L, G, H]$ +result: solve_hori_deprit_recurrence( + 2, /* max perturbation order */ + [H0, U2, 0], /* Hamiltonian list [F0, F1, F2] */ + vars, /* canonical variables */ + f, /* averaging variable */ + n_mean, /* fundamental frequency */ + 2, /* max order for generating functions */ + true /* closed-form mode */ +)$ + +/* Extract transformed Hamiltonians and generating functions */ +f_star_list: result[1]$ +s_gen_list: result[2]$ +``` + +### Example: Lagrange Planetary Equations + +```maxima +load("planetaryequations.mac")$ + +/* Compute rates of change under a disturbing function R */ +eqs: lagrange_planetary_equations(R_pot, a, e, i, h, g, M, n)$ +``` + +## Companion Package + +`symcelmech` is designed to work alongside [CelestialMechanics.jl](https://github.com/guilhermedeolivpaes/CelestialMechanics.jl), a high-performance numerical toolkit in Julia for orbital propagation and validation of analytical theories. Together they form a hybrid symbolic-numerical ecosystem for mission design and perturbation analysis. + +## Citation + +If you use `symcelmech` in your research, please cite the accompanying paper (under review): + +```bibtex +@article{paes2026symcelmech, + author = {Paes, Guilherme de Oliveira and Berton, Lilian and Vilhena de Moraes, Rodolpho}, + title = {A Hybrid Symbolic-Numerical Framework for Artificial Satellite + Theory and Dynamics using Maxima and Julia}, + year = {2026}, + note = {Submitted to Springer Nature} +} +``` + +## Acknowledgements + +Parts of the documentation, code comments, test scripts, and some code optimizations were developed with assistance from AI language models (Claude/Anthropic and Gemini/Google) and reviewed by the author. All scientific content, mathematical formulations, and architectural decisions are the author's own work. + +## Author + +**Guilherme de Oliveira Paes** + +- Mathematician +- MSc in Orbital Dynamics — National Institute for Space Research (INPE) +- PhD candidate in Computer Science — Federal University of São Paulo +- Contact: oliveira.guilherme1643@gmail.com / guilherme.paes@unifesp.br + +## Repository Structure + +``` +symcelmech/ +├── symcelmech.mac # entry point (load this) +├── analyticalfunctions.mac # Poisson brackets, quadrature, trigonometric tools +├── algebraicutils.mac # memoization, sanitization, substitution rules +├── astrodata.mac # physical constants database +├── average.mac # Hansen coefficients, orbital averaging +├── celmechseries.mac # elliptic expansion series +├── derivepotential.mac # perturbing potential derivation +├── lietransformations.mac # Hori-Deprit solver, Lie corrections +├── planetaryequations.mac # Lagrange, Gauss, Hamilton equations +├── rtest_symcelmech.mac # regression test suite +├── README.md +├── LICENSE +└── .gitignore +``` + +## Submitting to maxima-packages + +This repository is structured for submission to the official [maxima-packages](https://github.com/maxima-project-on-github/maxima-packages) community repository. To submit: + +1. Fork `maxima-project-on-github/maxima-packages` +2. Create a `symcelmech/` subdirectory in your fork +3. Copy all `.mac` files and the `rtest_symcelmech.mac` test script +4. Open a pull request + +## License + +This project is licensed under the [GNU General Public License v2.0](LICENSE). \ No newline at end of file diff --git a/symcelmech/algebraicutils.mac b/symcelmech/algebraicutils.mac new file mode 100644 index 0000000..84a11b9 --- /dev/null +++ b/symcelmech/algebraicutils.mac @@ -0,0 +1,467 @@ +/* algebraicutils.mac */ + +/* turn off reset warnings to load the package cleanly */ +suppress_trace_warnings: true$ + +/* adds the directory of this source file (.mac) to the maxima + search path (file_search_maxima). + this ensures that the interpreter can locate other internal + packages (such as series_celestial_mechanics.mac) residing + in the same folder, regardless of the worksheet's location. */ +file_search_maxima: cons(pathname_directory(load_pathname), file_search_maxima)$ + +load("analyticalfunctions.mac")$ + +/* ============================================================================ + FUNCTION: setup_pb_cache + Description: + initializes and manages the memoization system for poisson bracket calculations. + + Parameters: + vlist: the list of canonical variables to be used by the poisson bracket engine. + is_closed: routing flag for the poisson bracket engine. + ============================================================================ */ +setup_pb_cache(vlist, is_closed) := ( + kill(_pb_hash, pb_memo), + + /* funmake creates the function structure by injecting the content of vlist and the flag in a closed form */ + define(_pb_hash[f1, f2], funmake('poisson_bracket, [f1, f2, vlist, is_closed])), + + /* logical wrapper for managing antisymmetry and serial formats */ + pb_memo(A, B) := block( + [gA, gB], + /* protection against void terms */ + if A = 0 or B = 0 then return(0), + + /* converts to a general representation to avoid breaking the hash */ + gA : ratdisrep(A), + gB : ratdisrep(B), + + if gA = gB then return(0), + + if orderlessp(gA, gB) then + _pb_hash[gA, gB] + else + -_pb_hash[gB, gA] + ) +)$ + +/* memoization in maxima, using square brackets instead of parentheses, transforms the function into an array function + or a table lookup function. unlike a regular function that reprocesses the code every time, + a function with memoization works like a notebook: the first time you request a calculation, maxima calculates and + notes the result; in subsequent times, it simply reads what it noted. */ + +/* ============================================================================ + FUNCTION: deatomize_parameters + Description: + returns a list of substitutions for orbital and physical parameters. + + Parameters: + none + ============================================================================ */ +deatomize_parameters() :=block( + return( + [ + cos_i=cos(i), + sin_i=sin(i), + eta^2=1-e^2, + eta=sqrt(1-e^2), + eta_3^2=1-e_3^2, + eta_3=sqrt(1-e_3^2), + eta_o^2=1-e_o^2, + eta_o=sqrt(1-e_o^2), + k_o=mu_o*beta, + theta=omega_rot*t, + w_rot=omega_rot, + n=sqrt(mu/(a^3)), + p=a*(1-e^2) + ] + ) +)$ + +/* in lisp family languages like maxima, the block function always automatically returns the result of the last evaluated row. + this means we dont even need to use the return keyword, and we can simplify the entire function to a single line without declaring local variables. */ + +/* ============================================================================ + FUNCTION: deatomize_zonal_harmonics + Description: + returns a list of substitutions for zonal harmonics expansion. + + Parameters: + n_0: initial index of the expansion. + n_f: final index of the expansion. + ============================================================================ */ +deatomize_zonal_harmonics(n_0, n_f) :=block( + makelist( + concat('k_j, n) = concat('j, n)*R^n*mu, + n, + n_0, + n_f + ) +)$ + +/* flatten to transform the sublists into a single flat list that the subst function can read. + to integrate s afterwards, simply put it into a list: + [ + concat('k_c, n, m) = concat('c, n, m) * R^n * mu, + concat('k_s, n, m) = concat('s, n, m) * R^n * mu + ] */ + +/* ============================================================================ + FUNCTION: deatomize_tesseral_harmonics + Description: + returns a flat list of substitutions for tesseral harmonics. + + Parameters: + n_0: initial index of the expansion. + n_f: final index of the expansion. + ============================================================================ */ +deatomize_tesseral_harmonics(n_0, n_f) := block( + flatten( + makelist( + makelist( + concat('k_c, n, m) = concat('c, n, m) * R^n * mu, + m, 1, n + ), + n, n_0, n_f + ) + ) +)$ + + +/* ============================================================================ + FUNCTION: subst_for_julia + Description: + provides predefined sets of symbolic substitutions to bridge the gap between + analytical derivations in maxima and numerical implementation in julia. + + Parameters: + name_string: a string identifying the desired substitution set. + ============================================================================ */ +subst_for_julia(name_string) := block( + [name_up], + name_up: supcase(name_string), + + /* the string must be entirely uppercase in the function. */ + + if name_up = "SRP" then ( + return([ + a_o = a_sun, + e_o = e_sun, + i_o = i_sun, + mu_o = mu_sun, + l_o = l_sun, + g_o = g_sun, + h_o = h_sun, + eta_o = eta_sun, + lambda_o = lambda_sun, + n_o = sqrt(mu_sun/(a_sun^3)) + ]) + ) elseif name_up = "DELTOKEP" then ( + return([ + L = sqrt(mu*a), + G = sqrt(mu*a) * sqrt(1 - e^2), + H = sqrt(mu*a) * sqrt(1 - e^2) * cos_i + ]) + ) elseif name_up = "KEPTODEL" then ( + return([ + a = L^2/mu, + e = sqrt(1 - (G/L)^2), + cos_i = H/G, + sin_i = sqrt(1 - (H/G)^2) + ]) + ) elseif name_up = "TIME" then ( + return([ + l_3 = n_3 * t + l_3_0, + l_sun = n_sun * t + l_sun_0 + ]) + )else ( + print("error: replacement for ", name_string, " not found in the database."), + return([]) + ) +)$ + + +/* ============================================================================ + FUNCTION: coeff_of_term + Description: + extracts the cumulative numerical coefficient from a symbolic term. + + Parameters: + u: the symbolic term to be analyzed. + ============================================================================ */ +coeff_of_term(u) := block( + [c: 1], + if numberp(u) then return(u), + if not(mapatom(u)) and op(u) = "*" then ( + for f in args(u) do ( + if numberp(f) then c: c * f + else if floatnump(f) then c: c * f + ) + ) else if not(numberp(u)) then ( + /* if it is neither a number nor a product, the coefficient is 1. */ + c: 1 + ), + return(c) +)$ + +/* ============================================================================ + FUNCTION: chop_series + Description: + a precision filter used to control expression swell in high order series + expansions by removing terms smaller than a specified tolerance. + + Parameters: + expr: the symbolic series or expression to be filtered. + tol: the numerical threshold for truncation. + ============================================================================ */ +chop_series(expr, tol) := block( + [temp: expr, c], + if mapatom(expr) then ( + if numberp(expr) and abs(expr) < tol then return(0) else return(expr) + ), + + if op(expr) = "+" then ( + temp: 0, + for termo in args(expr) do ( + c: coeff_of_term(termo), + /* it only filters if it is a comparable number */ + if numberp(c) then ( + if abs(c) > tol then temp: temp + termo + ) else ( + temp: temp + termo /* keep going if you cannot determine the number */ + ) + ) + ) else ( + c: coeff_of_term(expr), + if numberp(c) and abs(c) < tol then temp: 0 + ), + return(temp) +)$ + +/* ============================================================================ + FUNCTION: sanitize_hamiltonian_expression_old + Description: + sanitizes, simplifies, and standardizes symbolic hamiltonian expressions. + + Parameters: + expr: the symbolic expression to be processed. + is_closed: flag for applying advanced simplifications for closed form theories. + ecc_order: order for taylor expansion in eccentricity, or false. + ============================================================================ */ +sanitize_hamiltonian_expression_old(expr, is_closed, ecc_order) := block( + [temp: expr, tol: 1e-15], + + + temp : chop_series(temp, tol), + + if is_closed then ( + /* kill sin squared plus cos squared and group common denominators. */ + temp : trigsimp(ratsimp(temp)), + + /* factor polynomials */ + temp : factor(temp), + + temp : ratsubst(eta^2, 1-e^2, temp), + temp : ratsubst(eta, G/L, temp) + ), + + + temp : ratsubst(G*sin_i, sqrt(G^2 - H^2), temp), + temp : ratsubst(L*e, sqrt(L^2 - G^2), temp), + temp : ratsubst(L*e, sqrt(L-G)*sqrt(L+G), temp), + + /* unifies inclination and semi axis. */ + temp : subst([ + cos(i) = cos_i, + sin(i)^2 = (1 - cos_i^2), + cos(2*i) = (2*cos_i^2 - 1), + cos(4*i) = (8*cos_i^2*(cos_i^2 - 1) + 1), + a = L^2/mu, + H = G*cos_i + ], temp), + + /* delaunay and s_gen identities */ + temp : ratsubst(1 - cos_i^2, sin_i^2, temp), + temp : ratsubst(e, sqrt(1-G^2/L^2), temp), + temp : ratsubst(sin_i, sqrt(1-H^2/G^2), temp), + temp : ratsubst(L^2*e^2, L^2 - G^2, temp), + temp : ratsubst(L^2*(1-e^2), G^2, temp), + + /* taylor expansion */ + if not(is_closed) and ecc_order # false then ( + /* the ratdisrep here guarantees that taylor does not work with cre garbage. */ + temp : ratdisrep(taylor(temp, e, 0, ecc_order)) + ), + + /* returning a general representation for the hash table will not cause an error */ + return(ratdisrep(temp)) +)$ + +/* ============================================================================ + FUNCTION: sanitize_hamiltonian_expression + Description: + sanitizes, simplifies, and standardizes symbolic hamiltonian expressions. + + Parameters: + expr: the symbolic expression to be processed. + is_closed: flag for applying advanced simplifications for closed form theories. + ecc_order: order for taylor expansion in eccentricity, or false. + ============================================================================ */ +sanitize_hamiltonian_expression(expr, is_closed, ecc_order) := block( + [temp: expr, tol: 1e-15], + + temp : chop_series(temp, tol), + + if is_closed then ( + temp : trigsimp(ratsimp(temp)), + temp : factor(temp), + temp : ratsubst(eta^2, 1-e^2, temp), + temp : ratsubst(eta, G/L, temp), + + /* closed form specific delaunay identities */ + temp : ratsubst(G*sin_i, sqrt(G^2 - H^2), temp), + temp : ratsubst(L*e, sqrt(L^2 - G^2), temp), + temp : ratsubst(L*e, sqrt(L-G)*sqrt(L+G), temp), + + temp : subst([ + cos(i) = cos_i, + sin(i)^2 = (1 - cos_i^2), + cos(2*i) = (2*cos_i^2 - 1), + cos(4*i) = (8*cos_i^2*(cos_i^2 - 1) + 1), + a = L^2/mu, + H = G*cos_i + ], temp), + + temp : ratsubst(1 - cos_i^2, sin_i^2, temp), + temp : ratsubst(e, sqrt(1-G^2/L^2), temp), + temp : ratsubst(sin_i, sqrt(1-H^2/G^2), temp), + temp : ratsubst(L^2*e^2, L^2 - G^2, temp), + temp : ratsubst(L^2*(1-e^2), G^2, temp) + + ) else ( + /* Series mode: minimal sanitization - no substitutions. + Variable mapping (a, H, cos(i)) is applied once per order + at the final f_star[n] level, not on every intermediate B_{n,k}. */ + if ecc_order # false then + temp : ratdisrep(taylor(temp, e, 0, ecc_order)) + ), + + return(ratdisrep(temp)) +)$ + +/* ============================================================================ + FUNCTION: clean_pb + Description: + simplifies the results of poisson bracket operations by enforcing the + fundamental kinematic relationship between eccentricity and eta. + + Parameters: + pb: the raw output from a poisson bracket calculation. + ============================================================================ */ +clean_pb(pb) := block( + [pb_simp, pb_eta, e_s, eta_s], + + e_s : gensym("e_"), + eta_s : gensym("eta_"), + + pb_simp : subst([e = e_s, eta = eta_s], pb), + pb_simp : trigsimp(ratsimp(pb_simp)), + pb_eta : ratsubst(eta_s^2, 1 - e_s^2, factor(pb_simp)), + + return(subst([e_s = e, eta_s = eta], pb_eta)) +)$ + +/* ============================================================================ + FUNCTION: get_keplerian_to_delaunay_rules + Description: + returns a list of inverse relations from keplerian to delaunay elements. + + Parameters: + mu: the gravitational parameter. + ============================================================================ */ +get_keplerian_to_delaunay_rules(mu):=block( + + /* declares local variable */ + [delaunay_inv], + + /* inverse relations */ + assume(G>0, L>0), + + delaunay_inv: [ + a = L^2/mu, + e^2 = 1 - (G/L)^2, + e = sqrt(1 - (G/L)^2), + cos(i) = H/G, + sin(i) = sqrt(1 - (H/G)^2), + cos(i)^2 = (H/G)^2, + sin(i)^2 = 1 - (H/G)^2, + /* multiples in case trigexpand does not break everything */ + cos(2*i) = 2*(H/G)^2 - 1, + sin(2*i) = 2*(H/G)*sqrt(1 - H^2/G^2), + cos(3*i) = 4*(H/G)^3 - 3*(H/G), + cos(4*i) = 8*(H/G)^4 - 8*(H/G)^2 + 1, + sin(3*i) = 3*sqrt(1 - H^2/G^2) - 4*(sqrt(1 - H^2/G^2))^3, + omega = g, + Omega = h, + M = l, + M_3 = l_3, + Omega_3 = h_3, + omega_3 = g_3, + M_o = L_o, + omega_o = g_o, + Omega_o = h_o + ], + + return(delaunay_inv) +)$ + +/* ============================================================================================================= */ +/* ============================================================================================================= */ +/* ============================================================================================================= */ + +/* ============================================================================ + FUNCTION: get_delaunay_to_kepler_rules + Description: + returns a list of mapping rules from delaunay variables to keplerian elements. + + Parameters: + mu_val: the gravitational parameter. + d_list: the list of delaunay elements. + k_list: the list of keplerian elements. + ============================================================================ */ +get_delaunay_to_kepler_rules(mu_val, d_list, k_list) := block( + [L, G, H, l, g, h, a, e, i, M, w, W, local_rules], + + /* unpacks the lists */ + [L, G, H, l, g, h] : d_list, + [a, e, i, M, w, W] : k_list, + + local_rules : [ + /* canonical momenta to keplerian */ + L^8 = mu_val^4 * a^4, + L^6 = mu_val^3 * a^3, + L^4 = mu_val^2 * a^2, + L^2 = mu_val * a, + L = sqrt(mu_val * a), + + G^4 = mu_val^2 * a^2 * (1 - e^2)^2, + G^2 = mu_val * a * (1 - e^2), + G = sqrt(mu_val * a * (1 - e^2)), + + H^2 = mu_val * a * (1 - e^2) * cos(i)^2, + H = sqrt(mu_val * a * (1 - e^2)) * cos(i), + + /* mixed relations */ + G^2/L^2 = (1 - e^2), + H^2/G^2 = cos(i)^2, + sqrt(1 - G^2/L^2) = e, + + /* canonical angles to kepler */ + l = M, + g = w, + h = W + ], + + return(local_rules) +)$ diff --git a/symcelmech/analyticalfunctions.mac b/symcelmech/analyticalfunctions.mac new file mode 100644 index 0000000..5f7519d --- /dev/null +++ b/symcelmech/analyticalfunctions.mac @@ -0,0 +1,635 @@ +/* analyticalfunctions.mac */ + +/* turn off reset warnings to load the package cleanly */ +redefine_warn: false$ + +file_search_maxima: cons(pathname_directory(load_pathname), file_search_maxima)$ + +kill(a, e, i, f, l, g, h, L, G, H)$ + + +/* ============================================================================ + FUNCTION: setup_canonical_dependencies + Description:configures the symbolic differentiation engine of maxima by + * establishing functional dependencies and custom derivative rules gradefs for + * orbital mechanics. this allows the poisson bracket engine to automatically apply + * the chain rule across different coordinate systems e.g. from delaunay + * actions to keplerian geometry. + * key configurations + * 1. delaunay kepler mapping defines derivatives of eccentricity e + * semi major axis a and inclination i with respect to delaunay + * actions l g h. + * 2. orbital geometry dna establishes the exact partial derivatives for + * the radius r and true anomaly f. + * 3. closed form jacobian implements the exact transformation between the + * true anomaly f and mean anomaly l including the df de variation + * at constant l. + * 4. homogeneous formalism sets up time dependent rules for third body + * anomalies l_3 l_o and planetary rotation theta relative to the + * independent variable tau. + * technical note + * this setup is essential for closed form theories where derivatives + * involving the true anomaly must be handled without series expansion. + Parameters: none + * ============================================================================ */ +setup_canonical_dependencies() := ( + depends(a, L), depends(e, [L, G]), depends(i, [G, H]), + depends(r, [a, e, f]), depends(sin_i, [G, H]), + depends(cos_i, [G, H]), depends(eta, [L, G]), + depends(f, [l, e]), + + /* derivative of e with respect to l */ + gradef(e, L, (1 - e^2) / (L * e)), + + /* derivative of e with respect to g */ + gradef(e, G, -G / (L^2 * e)), + + /* derivative of sin i with respect to g */ + gradef(sin_i, G, ( H^2) / (G^3 * sin_i)), + /* derivative of sin i with respect to h */ + gradef(sin_i, H, (-H) / (G^2 * sin_i)), + + /* derivative of cos i with respect to g */ + gradef(cos_i, G, -H/(G^2)), + /* derivative of cos i with respect to h */ + gradef(cos_i, H, 1/G), + + /* delaunay to kepler */ + /* derivative of a with respect to l */ + gradef(a, L, 2*a/L), + /* derivative of i with respect to g */ + gradef(i, G, cot(i)/G), + /* derivative of i with respect to h */ + gradef(i, H, -1/(G * sin(i))), + + /* ----------------------------- */ + gradef(r, a, r/a), + /* derivative of r with respect to e */ + gradef(r, e, -r * (cos(f) + e) / (1 - e^2)), + /* derivative of r with respect to f */ + gradef(r, f, (r^2 * e * sin(f)) / (a * (1 - e^2))), + + /* ----------------------------- */ + /* derivative of eta with respect to l */ + gradef(eta, L, (-G^2) / (L^3 * eta)), + /* derivative of eta with respect to g */ + gradef(eta, G, G / (L^2 * eta)), + + /* ----------------------------- */ + /* rules for homogeneous formalism */ + /* truncated mode in eccentricity */ + /* derivative of l_3 with respect to t */ + gradef(l_3, tau, n_3), + /* derivative of l_o with respect to t */ + gradef(l_o, tau, n_o), + /* derivative of l with respect to t */ + gradef(l, tau, n), + + /* for tesseral harmonics in the inertial frame of reference */ + gradef(theta, tau, w_rot), /* theta equals w_rot tau */ + + print("canonical dependencies activated.") +)$ + +/* ============================================================================ + FUNCTION: teardown_canonical_dependencies + Description:cleans the symbolic environment by removing all functional + * dependencies previously established by setup canonical dependencies. + * usage should be called at the end of a solver run or when switching between + * different dynamical models to prevent unwanted derivative cross talk or + * memory overhead in the symbolic tree. + Parameters: none + * ============================================================================ */ +teardown_canonical_dependencies() := ( + + remove(f, dependency), remove(r, dependency), + remove(a, dependency), remove(e, dependency), + remove(eta, dependency), remove(i, dependency), + remove(sin_i, dependency), remove(cos_i, dependency), + + + remove(l_3, dependency), remove(l_o, dependency), + remove(l, dependency), remove(theta, dependency), + + print("canonical dependencies deactivated.") +)$ + + +/* ============================================================================ + FUNCTION: pb_series + Description:an optimized version of the series poisson bracket engine. it adopts the + * linear extraction strategy from the closed form engine to improve + * performance in high order expansions. by decoupling symbolic differentiation + * from algebraic summation it reduces the computational overhead on the + * expression tree. + Parameters: + * a: symbolic series to be bracketed + * b: symbolic series to be bracketed + * vars: canonical variables q p + * ============================================================================ */ +pb_series(A, B, Vars) := block( + [ + n_dims: length(Vars)/2, + Aq: [], Ap: [], Bq: [], Bp: [], + PB: 0, i + ], + + if A = B then return(0), + + /* 1. mass extraction effect pb closed */ + /* we extract all derivatives before starting any algebra */ + for i: 1 thru n_dims do ( + Aq: endcons(diff(A, Vars[i]), Aq), + Ap: endcons(diff(A, Vars[i + n_dims]), Ap), + Bq: endcons(diff(B, Vars[i]), Bq), + Bp: endcons(diff(B, Vars[i + n_dims]), Bp) + ), + + /* 2. linearized sum */ + /* we use ratexpand to keep the polynomial series structure clean */ + for i: 1 thru n_dims do ( + PB: PB + Aq[i]*Bp[i] - Ap[i]*Bq[i] + ), + + return(ratexpand(PB)) +)$ + + +/* ============================================================================ + FUNCTION: pb_closed_form + Description:a high performance poisson bracket engine optimized for closed form orbital + * theories. unlike the standard series pb this function treats the radius r + * and true anomaly f as independent symbolic entities manually applying + * the analytical chain rule to map them back to delaunay actions. + Parameters: + * a: symbolic expressions involving orbital parameters + * b: symbolic expressions involving orbital parameters + * vars: canonical variable list + * ============================================================================ */ +pb_closed_form(A, B, Vars) := block( + [ + PB: 0, + Ar, Af, Aa, Ae, Aeta, Ai, Ag, Ah, AG, AH, + Br, Bf, Ba, Be, Beta, Bi, Bg, Bh, BG, BH, + df_dl, dr_dl, dr_dL, dr_dG, df_dL, df_dG, + da_dL, de_dL, de_dG, deta_dL, deta_dG, di_dG, di_dH, + C_rf, C_ra, C_re, C_reta, C_fa, C_fe, C_feta, + C_rG, C_fG, C_eG, C_etaG, C_iG, C_iH, + extended: is(length(Vars) = 8) + ], + + if A = B then return(0), + + /* 2. explicit derivatives independent */ + Ar: diff(A, r), Af: diff(A, f), Aa: diff(A, a), Ae: diff(A, e), + Aeta: diff(A, eta), Ai: diff(A, i) + diff(A, sin_i)*cos_i - diff(A, cos_i)*sin_i, + Ag: diff(A, g), Ah: diff(A, h), AG: diff(A, G), AH: diff(A, H), + + Br: diff(B, r), Bf: diff(B, f), Ba: diff(B, a), Be: diff(B, e), + Beta: diff(B, eta), Bi: diff(B, i) + diff(B, sin_i)*cos_i - diff(B, cos_i)*sin_i, + Bg: diff(B, g), Bh: diff(B, h), BG: diff(B, G), BH: diff(B, H), + + /* 3. monomial kinematics */ + /* we eliminate the binomial 1 e cos f using the orbit equation. + purely monomial denominators avoid gcd explosion in ram */ + + df_dl : (a^2 * eta) / r^2, + dr_dl : (a * e * sin(f)) / eta, + + da_dL : 2*a/L, + de_dL : eta^2 / (L * e), + deta_dL : -eta / L, + + de_dG : -eta / (L * e), + deta_dG : 1 / L, + + di_dG : cos_i / (L * eta * sin_i), + di_dH : -1 / (L * eta * sin_i), + + /* radial and angular derivatives in l and g keeping l constant */ + dr_dL : 2*r/L - ((a - r) * eta^2) / (L * e^2), + dr_dG : ((a - r) * eta) / (L * e^2), + + df_dL : (sin(f) / (L * e)) * (1 + (a * eta^2) / r), + df_dG : -(sin(f) / (L * e * eta)) * (1 + (a * eta^2) / r), + + df_de_exact : (sin(f) / (1 - e^2)) * (2 + e * cos(f)), + dr_de_exact : -a * cos(f), + + + /* 4. pre expanded fundamental scalars */ + C_rf : ratexpand( dr_dl * df_dL - dr_dL * df_dl ), + C_ra : ratexpand( dr_dl * da_dL ), + C_re : ratexpand( dr_dl * de_dL ), + C_reta : ratexpand( dr_dl * deta_dL ), + C_fa : ratexpand( df_dl * da_dL ), + C_fe : ratexpand( df_dl * de_dL ), + C_feta : ratexpand( df_dl * deta_dL ), + + C_rG : ratexpand(dr_dG), + C_fG : ratexpand(df_dG), + C_eG : ratexpand(de_dG), + C_etaG : ratexpand(deta_dG), + C_iG : ratexpand(di_dG), + C_iH : ratexpand(di_dH), + + /* 5. the factored master equation */ + PB : ratexpand( + /* block l l */ + C_rf * (Ar*Bf - Af*Br) + + C_ra * (Ar*Ba - Aa*Br) + + C_re * (Ar*Be - Ae*Br) + + C_reta * (Ar*Beta - Aeta*Br) + + C_fa * (Af*Ba - Aa*Bf) + + C_fe * (Af*Be - Ae*Bf) + + C_feta * (Af*Beta - Aeta*Bf) + + /* block g g */ + + (Ag*BG - AG*Bg) + + C_rG * (Ag*Br - Ar*Bg) + + C_fG * (Ag*Bf - Af*Bg) + + C_eG * (Ag*Be - Ae*Bg) + + C_etaG * (Ag*Beta - Aeta*Bg) + + C_iG * (Ag*Bi - Ai*Bg) + + /* block h h */ + + (Ah*BH - AH*Bh) + + C_iH * (Ah*Bi - Ai*Bh) + ), + + if extended then block( + [Atau, AT, Btau, BT, Al3, Bl3, Ath, Bth, + tau_var: Vars[4], T_var: Vars[8]], + + /* pure partial derivatives */ + Al3: diff(A, l_3), Bl3: diff(B, l_3), + Ath: diff(A, theta), Bth: diff(B, theta), + AT: diff(A, T_var), BT: diff(B, T_var), + + /* chain rule d dtau equals n_3 d dl_3 plus w_rot d dtheta */ + Atau: n_3*Al3 + w_rot*Ath, + Btau: n_3*Bl3 + w_rot*Bth, + + PB: PB + ratexpand(Atau*BT - AT*Btau), + + PB : ratsubst(df_de_exact, 'diff(f, e), PB), + PB : ratsubst(dr_de_exact, 'diff(r, e), PB) + ), + + return(PB) +)$ + + +/* ============================================================================ + FUNCTION: poisson_bracket + Description:top level router for poisson bracket calculations. decouples the lie + * solver logic from the specific kinematic implementation. + Parameters: + * a: expressions to be bracketed + * b: expressions to be bracketed + * vars: canonical variables + * is_closed_form: routing flag. if true triggers pb closed form + * if false triggers pb series. + * ============================================================================ */ +poisson_bracket(A, B, Vars, is_closed_form) := block( + if is_closed_form then + return(pb_closed_form(A, B, Vars)) + else + return(pb_series(A, B, Vars)) +)$ + + +/* ============================================================================ + FUNCTION: generating_function_quadrature + Description:performs a high performance symbolic integration quadrature of periodic + * hamiltonian terms to compute the generating function s. instead of using + * the general purpose integrate engine it employs a frequency extraction + * method in the complex exponential domain which is faster and more reliable + * for trigonometric series. + Parameters: + * expr: the periodic part of the hamiltonian to be integrated + * var: the integration variable e.g. f or l + * ============================================================================ */ +generating_function_quadrature(expr, var) := block( + [temp, int_val, lista_termos, e_s, a_s, eta_s], + + e_s : gensym("e_"), a_s : gensym("a_"), eta_s : gensym("eta_"), + temp : subst([e = e_s, a = a_s, eta = eta_s], expr), + + temp : if freeof(%i, temp) then expand(exponentialize(temp)) else temp, + if (temp = 0) then return(0), + + lista_termos : if not atom(temp) and op(temp) = "+" then args(temp) else [temp], + + int_val : apply("+", map(lambda([u], + block([freq], + freq : ratsimp(diff(u, var) / u), + if freq = 0 then 0 else u / freq + ) + ), lista_termos)), + + temp : realpart(rectform(int_val)), + + /* restores */ + temp : subst([e_s = e, a_s = a, eta_s = eta], temp), + return(temp) +)$ + + + + +/* defines the pattern for any expression */ +matchdeclare(aa, all)$ + +/* ============================================================================ + FUNCTION: celmec_trigexpand + Description:a specialized trigonometric expansion tool tailored for celestial mechanics. + * it targets a specific variable usually an anomaly like f or l and + * decomposes multiple angle trigonometric terms into powers of sines and + * cosines while applying addition formulas to separate fast and slow + * arguments. + Parameters: + * expr: the symbolic expression to be expanded + * v_target: the target variable e.g. true anomaly f whose + * multiple angles will be expanded into powers + * ============================================================================ */ +celmec_trigexpand(expr, v_target) := block( + [rule_c, rule_s], + + /* activates the trigonometric expansion flags locally */ + local(trigexpandplus, trigexpandtimes), + trigexpandplus: true, + trigexpandtimes: true, + + /* cosine rule cos a plus n f */ + rule_c: defrule(rc, cos(aa), + block([c, rest, cv, sv], + c: coeff(aa, v_target), + rest: expand(aa - c*v_target), + /* expands cos c f and sin c f into powers */ + cv: trigexpand(cos(c*v_target)), + sv: trigexpand(sin(c*v_target)), + cos(rest)*cv - sin(rest)*sv + ) + ), + + /* sin rule sin a plus n f */ + rule_s: defrule(rs, sin(aa), + block([c, rest, cv, sv], + c: coeff(aa, v_target), + rest: expand(aa - c*v_target), + /* expands cos c f and sin c f into powers */ + cv: trigexpand(cos(c*v_target)), + sv: trigexpand(sin(c*v_target)), + sin(rest)*cv + cos(rest)*sv + ) + ), + + apply1(expr, rc, rs) +)$ + +/* ============================================================================ + * FUNCTION: eliminate_parallax + * ============================================================================ + * Description: + * Implements Deprit's elimination of the parallax (Deprit 1981) reformulated + * in Delaunay variables following Lara, San-Juan & Lopez-Ochoa (2014, 2015). + * + * The transformation removes all short-period terms from the Hamiltonian + * EXCEPT those proportional to 1/r^2, by exploiting the identity: + * 1/r^m = (1/r^2) * [(1 + e*cos(f))^(m-2)] / (a*eta^2)^(m-2), m > 2 + * + * Unlike a standard Delaunay normalization (which averages over l or f), + * the parallax elimination only removes the EXPLICIT appearance of f, + * preserving the implicit dependence through r. The resulting Hamiltonian + * H_{0,1} still contains short-period terms via the factor 1/r^2. + * + * Reference: + * - Deprit, A. (1981). Celest. Mech. 24(2), 111-153. + * - Lara, M., San-Juan, J.F., Lopez-Ochoa, L.M. (2015). Celest. Mech. Dyn. Astron. + * - Coffey, S. & Deprit, A. (1982). J. Guidance 5(4), 366-371. + * + * Parameters: + * U_raw (expr): The raw perturbing potential in terms of (r, f, g, i). + * Must be proportional to 1/r^(n+1) for degree n. + * f_var (symbol): The true anomaly variable (typically 'f'). + * r_var (symbol): The radius variable (typically 'r'). + * a_var (symbol): The semi-major axis variable. + * e_var (symbol): The eccentricity variable. + * eta_var (symbol): The eccentricity function sqrt(1-e^2). + * n_mean (expr) : The Keplerian mean motion (mu^2/L^3). + * + * Returns: + * A list [H_star, W_parallax] where: + * H_star : The simplified Hamiltonian (terms free of explicit f, + * still containing 1/r^2). + * W_parallax : The generating function of the parallax transformation. + * + * Usage Example: + * >> U2: derive_zonal_potential(2, mu, r, f, a, e, i, g, R, true); + * >> [H_par, W_par]: eliminate_parallax(U2, f, r, a, e, eta, mu^2/L^3); + * ============================================================================ */ +eliminate_parallax(U_raw, f_var, r_var, a_var, e_var, eta_var, n_mean) := block( + [ + U_reduced, U_expanded, terms_list, H_star, H_per, + W_par, term, jac_df_dl + ], + + /* --------------------------------------------------------------- + * Step 1: Reduce all powers of 1/r to 1/r^2 using the identity + * 1/r^m = (1/r^2) * (1 + e*cos(f))^(m-2) / (a*eta^2)^(m-2) + * + * Strategy: multiply U_raw by r^2, substitute r = a*eta^2/(1+e*cos(f)), + * expand trigonometrically, then divide back by r^2. + * This converts all (1/r^m) factors into (1/r^2) * polynomial_in_cos_f. + * --------------------------------------------------------------- */ + + /* Multiply by r^2 to clear the 1/r^2 factor */ + U_reduced: ratsimp(U_raw * r_var^2), + + /* Now substitute r = a*eta^2/(1+e*cos(f)) in the REMAINING r dependence. + After multiplying by r^2, the original 1/r^3 becomes 1/r, + the 1/r^4 becomes 1/r^2 -> 1, etc. + So the remaining r factors get replaced. */ + U_reduced: subst(r_var = a_var * eta_var^2 / (1 + e_var*cos(f_var)), U_reduced), + U_reduced: ratsimp(U_reduced), + + /* trigreduce to linearize into cos(kf + mg) form and simplification */ + U_reduced: trigreduce(U_reduced), + U_reduced: ratexpand(U_reduced), + + /* --------------------------------------------------------------- + * Step 2: Separate into parallax-free and periodic parts. + * + * H_star gets terms FREE of the explicit true anomaly f. + * H_per gets terms containing f explicitly. + * + * IMPORTANT: This is NOT an average. H_star retains the 1/r^2 + * factor (which will be restored), preserving implicit f-dependence. + * --------------------------------------------------------------- */ + terms_list: args(U_reduced), + H_star: 0, + H_per: 0, + + for term in terms_list do ( + if freeof(f_var, term) then + H_star: H_star + term + else + H_per: H_per + term + ), + + /* Restore the 1/r^2 factor */ + H_star: H_star / r_var^2, + H_per: H_per / r_var^2, + + /* Simplify */ + H_star: ratsimp(H_star), + H_per: ratsimp(H_per), + + /* --------------------------------------------------------------- + * Step 3: Compute the generating function W by quadrature. + * + * From the homological equation: n * dW/dl = H_per + * Using the relation: a^2 * eta * dl = r^2 * df (Eq. 21 of Lara) + * We get: dW/df = (H_per * r^2) / (n * a^2 * eta) + * + * But H_per already has the 1/r^2 factor, so: + * dW/df = H_per * r^2 / (n * a^2 * eta) + * = [terms_periodic / r^2] * r^2 / (n * a^2 * eta) + * = terms_periodic / (n * a^2 * eta) + * + * The integrand is a polynomial in cos(f) and sin(f), so the + * quadrature is elementary. + * --------------------------------------------------------------- */ + + /* Form the integrand: H_per * r^2 / (n_mean * a^2 * eta) + = (periodic_terms/r^2) * r^2 / (n_mean * a^2 * eta) + = periodic_terms / (n_mean * a^2 * eta) */ + + /* Use the generating_function_quadrature on the periodic part + after multiplying by the Jacobian factor */ + W_par: -generating_function_quadrature(H_per * r_var^2, f_var) / (n_mean * a_var^2 * eta_var), + + /* Sanitize */ + W_par: sanitize_hamiltonian_expression(W_par, true, false), + H_star: sanitize_hamiltonian_expression(H_star, true, false), + + print("=== Parallax Elimination Complete ==="), + print(" H_star terms:", nterms(H_star)), + print(" W_par terms:", nterms(W_par)), + print(" Contains r:", not freeof(r_var, H_star)), + print(" Free of f:", freeof(f_var, H_star)), + + return([H_star, W_par]) +)$ + +/* ============================================================================ + * keplerian geometry closed form + * ============================================================================ + Description:a collection of exact geometric relations for the two body problem. + * these functions facilitate the analytical mapping between the different + * anomalies mean eccentric and true and the radial distance without + * relying on series expansions. + * ============================================================================ */ + + +/* ============================================================================ + FUNCTION: sin_f_exact + Description:computes the sine of the true anomaly f as a function of the + * eccentric anomaly e and eccentricity e. + Parameters: + * f_symb: symbol for true anomaly + * e_symb: symbol for eccentricity + * e_symb_ecc: symbol for eccentric anomaly + * ============================================================================ */ +sin_f_exact(f_symb, e_symb, E_symb) := block( + [sin_f], + + sin_f: ((1-e_symb^2)^(1/2)/(1-e_symb*cos(E_symb)))*sin(E_symb), + + return(sin_f) +)$ + +/* ============================================================================ + FUNCTION: cos_f_exact + Description:computes the cosine of the true anomaly f as a function of the + * eccentric anomaly e and eccentricity e. + Parameters: + * f_symb: symbol for true anomaly + * e_symb: symbol for eccentricity + * e_symb_ecc: symbol for eccentric anomaly + * ============================================================================ */ +cos_f_exact(f_symb, e_symb, E_symb) := block( + + [cos_f], + + cos_f: (cos(E_symb) - e_symb)/(1-e_symb*cos(E_symb)), + + return(cos_f) +)$ + +/* ============================================================================ + FUNCTION: r_over_a_exact + Description:computes the normalized radial distance r a as a function of the + * eccentric anomaly e and eccentricity e. + Parameters: + * e_symb: symbol for eccentricity + * e_symb_ecc: symbol for eccentric anomaly + * ============================================================================ */ +r_over_a_exact(e_symb, E_symb) := block( + + [r_a], + + r_a: (1 - e_symb*cos(E_symb)), + + return(r_a) +)$ + +/* ============================================================================ + FUNCTION: r_two_body_exact + Description:computes the exact radial distance r in terms of the true anomaly f. + Parameters: + * e_symb: symbol for eccentricity + * a_symb: symbol for semi major axis + * f_symb: symbol for true anomaly + * ============================================================================ */ +r_two_body_exact(e_symb, a_symb, f_symb) := block( + [r], + + r: a_symb*(1-e_symb^2)/(1+e_symb*cos(f_symb)), + + return(r) +)$ + +/* ============================================================================ + FUNCTION: dm_f + Description:computes the exact jacobian dm df representing the ratio between the + * mean anomaly increment and the true anomaly increment. + Parameters: + * e_symb: symbol for eccentricity + * f_symb: symbol for true anomaly + * ============================================================================ */ +dM_f(e_symb, f_symb) := block( + [dM], + + dM: (1 - e_symb^2)^(3/2) / (1 + e_symb * cos(f_symb))^2, + + return(dM) +)$ + +/* ============================================================================ + FUNCTION: dm_avg_exact + Description:computes the exact jacobian dm de representing the ratio between the + * mean anomaly increment and the eccentric anomaly increment derived + * directly from keplers equation m equals e minus e sin e. + Parameters: + * e_symb: symbol for eccentricity + * e_symb_ecc: symbol for eccentric anomaly + * ============================================================================ */ +dM_avg_exact(e_symb, E_symb) := block( + [dM], + + dM: (1-e_symb*cos(E_symb)), + + return(dM) +)$ diff --git a/symcelmech/astrodata.mac b/symcelmech/astrodata.mac new file mode 100644 index 0000000..695280b --- /dev/null +++ b/symcelmech/astrodata.mac @@ -0,0 +1,134 @@ +/* astrodata.mac */ + +/* author: Guilherme de Oliveira Paes */ +/* contact: guilherme.paes@unifesp.br / oliveira.guilherme1643@gmail.com */ +/* 19/04/2026 */ +/* version: 1.0 */ + +suppress_trace_warnings: true$ + + +file_search_maxima: cons(pathname_directory(load_pathname), file_search_maxima)$ + +/* ============================================================================ + * FUNCTION: get_astro_params + * ============================================================================ + * Description: + * Acts as a centralized database of physical, gravitational, and orbital + * constants for various celestial bodies. It provides the necessary parameters + * for building complex Hamiltonians (Geopotential, Third-Body, SRP, etc.) + * with high-fidelity numerical values. + * + * Parameters: + * astro_name (string): The name of the celestial body (e.g., "Earth", + * "Moon", "Mercury", "Apophis", "Didymos"). + * The input is case-insensitive. + * + * Returns: + * A list of symbolic equations [parameter = value, ...] ready for use with + * 'subst'. Returns an empty list and prints an error if the body is not + * in the database. + * + * Standard Units: + * - Gravitational Parameter (mu): km^3/s^2 + * - Mean Radius (R): km + * - Rotation Rate (omega_rot): rad/s + * - Distances (d_AU / a_sun): AU or km (as specified in the return keys) + * - Angles (i): rad + * - Harmonics (j_n, c_nm, s_nm): Dimensionless + * + * Supported Bodies: + * 1. Planets/Moons: EARTH (WGS84/EGM), MOON (LP150Q), MERCURY (MESSENGER). + * 2. Asteroids: APOPHIS (Shape model-derived), DIDYMOS (DART mission specs). + * 3. Orbital Systems: DIDYMOS_SYSTEM (Sun-relative orbital elements). + * ============================================================================ */ +get_astro_params(astro_name):= block( + [name_up], + name_up: supcase(astro_name), + + /* + mu = km^3/s^2; R = km omega_rot = rad/s; d_AU = AU, a_body = km, i = rad + */ + + /* the string must be entirely uppercase in the function */ + if name_up = "EARTH" then ( + return([ + mu=398600.4418, + R=6378.137, + a_sun_AU=1.0, + j2=0.00108263, + j3=-2.5327e-6, + j4=-1.6196e-6, + J5=-2.276e-7, + J6=5.4068e-7, + c22=1.5744e-6, + s22=-0.9038e-6, + omega_rot=7.2921151467e-5 + ]) + ) elseif name_up = "MOON" then ( + return([ + mu=4.90280012e3, R=1738.1, d_AU=1.0, a_earth_km=384400.0, e_earth=0.0549, + + j2=2.032365e-4, j3=8.585355e-6, j4=-9.860033e-6, j5=8.393511e-7, j6=-1.331149e-5, + j7=-2.205653e-5, j8=-9.434727e-6, j9=1.511956e-5, j10=3.863970e-6, j11=5.130603e-6, j12=9.103736e-6, + + c22=2.235491e-5, c31=2.849956e-5, c32=4.847437e-6, c33=1.712582e-6, + s22=1.535185e-8, s31=5.953103e-6, s32=1.689177e-6, s33=-2.564964e-7, + c41=-5.0879936e-7, c42=7.84175859e-8, c43=5.92099402e-8, c44=-3.98407411e-9, + s41=-4.49144872e-7, s42=1.48177868e-7, s43=-1.20077667e-8, s44=6.52571425e-9, + + omega_rot=2.6617e-6 + ]) + ) elseif name_up = "MERCURY" then ( + return([ + mu=22032.09, + R=2439.7, + j2=6.0e-5, + j3=1.188e-5, + j4=1.95e-5, + c22=1.24973e-5, + d_AU=0.387, + e_sun=0.206, + a_sun=5.789437596090001e7 + ]) + ) elseif name_up = "APOPHIS" then ( + return([ + mu = 2.86e-9, + R = 0.17, + j2 = 0.1344360544, + j3 = -0.0303639272, + j4 = -0.0463008733, + c22 = 0.0472123118, + c42 = -0.0045027695, + omega_rot = 5.711155929300816e-5 + ]) + ) elseif name_up = "DIDYMOS" then ( + return([ + mu = 3.489324040e-8, + R = 0.39, + j2 = 8.35e-2, + j3 = -2.27e-2, + j4 = -5.5e-3, + c22 = 1.82e-3, + c31=-6.22e-3, + s31=-5.94e-3, + c32=-7.10e-5, + s32=1.01e-3, + c33=8.47e-5, + s33=-3.07e-4, + c42=-2.06e-5, + c44=5.91e-7, + omega_rot = 0.0007725088531821048 + ]) + ) elseif name_up = "DIDYMOS_SYSTEM" then ( + return([ + d_AU=1.642564399038672, + a_sun=2.4572406756e8, + e_sun=0.3832284321571334, + i_sun=0.059587 + ]) + )else ( + print("Error: Body '", astro_name, "' not found in the database."), + return([]) + ) +)$ \ No newline at end of file diff --git a/symcelmech/average.mac b/symcelmech/average.mac new file mode 100644 index 0000000..5a155d1 --- /dev/null +++ b/symcelmech/average.mac @@ -0,0 +1,313 @@ +/* average.mac */ + +suppress_trace_warnings: true$ +redefine_warn: false$ + +file_search_maxima: cons(pathname_directory(load_pathname), file_search_maxima)$ + + +load("analyticalfunctions.mac")$ + +assume(L > G, G > H, G > 0, L > 0)$ +assume(e > 0, e < 1, e3 > 0, e3 < 1)$ + +/* ============================================================================ + FUNCTION: leaf_count + Description:measures the structural complexity and symbolic weight of an expression + * by recursively counting the number of atoms leaves in its expression tree. + * while nterms only counts the number of summands leaf count provides a + * metric for the total computational cost and memory footprint of the terms. + * usage note critical for monitoring expression swell during high order lie + * transformations allowing the user to decide when to apply chop series + * or further simplifications. + Parameters: + * expr: the symbolic expression to be analyzed + * ============================================================================ */ +leaf_count(expr) := if atom(expr) then 1 else 1 + apply("+", map(leaf_count, args(expr)))$ + +/* ============================================================================ + FUNCTION: map_average + Description:acts as a symbolic filter to segregate a hamiltonian into its secular + * mean and periodic fluctuating components relative to a specific + * averaging variable for example true anomaly f or mean anomaly l. + * internal logic automatically handles single expressions sums or lists. + * iterates through each summand to perform a clean freeof check. + * this separation is the precursor to the homologous equation solver + * defining which terms contribute to the transformed energy and which + * to the generating function. + Parameters: + * ham_input: the hamiltonian expression to be processed + * var_avg: the variable of integration or averaging + * ============================================================================ */ +map_average(ham_input, var_avg) := block( + [ham_list, total_secular: 0, total_periodic: 0, res], + + /* converts the input into a list to iterate term by term */ + ham_list : if listp(ham_input) then ham_input + else if not atom(ham_input) and op(ham_input) = "+" then args(ham_input) + else [ham_input], + + /* iterates over each term of the hamiltonian to separate secular components from periodic ones according to the model closed or series */ + for h in ham_list do ( + + /* here h is already the term prepared by the prepare potential poisson function */ + if freeof(var_avg, h) then ( + total_secular : total_secular + h + ) else ( + total_periodic : total_periodic + h + ) + ), + + [total_secular, total_periodic] +)$ + + +/* ============================================================================ + * hansen coefficients and orbital averaging hughes method + Description:implements the calculation of hansen coefficients x0 n m e and the + * associated orbital averaging technique for terms of the form r n exp i m f. + * this module allows for exact analytical averaging over the true anomaly f + * by mapping trigonometric terms directly to eccentricity functions. + * ============================================================================ */ + +kill(hansen_x0_pos, hansen_x0_pos_calc)$ + +/* ============================================================================ + FUNCTION: hansen_x0_pos_calc + Description:recursively computes the hansen coefficients x0 n m e for positive + * powers of the radial distance n greater than 0. + * mathematical reference uses the recursion relations established by hughes 1977 + * and modified for stability in symbolic environments. + Parameters: + * n: power of the radial distance r n + * m: frequency of the true anomaly m f + * ecc: eccentricity e + * ============================================================================ */ +hansen_x0_pos_calc[n, m, ecc] := block( + [local_m : abs(m)], + + if n = 1 and local_m = 0 then return(expand(1 + ecc^2/2)), + if n = 2 and local_m = 0 then return(expand(1 + 3*ecc^2/2)), + if n = 1 and local_m = 1 then return(expand(-3*ecc/2)), + if n = 2 and local_m = 1 then return(expand(-2*ecc*(1 + ecc^2/4))), + + if local_m > 1 then ( + return(expand( + (1/(ecc*(n - local_m + 2))) * (ecc*(n + local_m) * hansen_x0_pos_calc[n, local_m - 2, ecc] + + 2*(local_m - 1) * hansen_x0_pos_calc[n, local_m - 1, ecc]) + )) + ), + + if n > 2 then ( + return(expand( + ((2*n + 1)/(n + 1)) * hansen_x0_pos_calc[n - 1, local_m, ecc] - + ((n - local_m)*(n + local_m))/(n*(n + 1)) * (1 - ecc^2) * hansen_x0_pos_calc[n - 2, local_m, ecc] + )) + ), + + return(0) +)$ + +/* ============================================================================ + FUNCTION: hansen_x0_pos + Description:alias for hansen x0 pos calc. recursively computes the hansen + * coefficients x0 n m e for positive powers of the radial distance. + Parameters: + * n: power of the radial distance r n + * m: frequency of the true anomaly m f + * ecc: eccentricity e + * ============================================================================ */ +hansen_x0_pos[n, m, ecc] := hansen_x0_pos_calc[n, m, ecc]$ + + +/* ============================================================================ + FUNCTION: hansen_x0_neg + Description:computes the hansen coefficients x0 minus n plus 1 m e for negative powers of + * the radial distance n less than 0. this is critical for averaging gravitational + * potentials 1 r n. method employs a closed form analytical binomial expansion. + Parameters: + * n: negative power magnitude actually passed as positive related to n plus 1 + * m: frequency + * e: eccentricity + * ============================================================================ */ +hansen_x0_neg(n, m, e) := block( + [p, coef, k, m_abs, j, limit], + p : n - 1, + m_abs : abs(m), + + /* physical nullity rule */ + if m_abs > p then return(0), + + /* analytical binomial expansion */ + /* k starts at m abs and jumps by 2 until p */ + limit : floor((p - m_abs)/2), + coef : sum( + binomial(p, m_abs + 2*j) * (e^(m_abs + 2*j) / 2^(m_abs + 2*j)) * binomial(m_abs + 2*j, j), + j, 0, limit + ), + + return(expand(radcan(coef / (1 - e^2)^(n - 1/2)))) +)$ + + + +/* ============================================================================ + FUNCTION: get_arg_safe + Description:safely takes the argument of a function cos or sin within an expression. + Parameters: + * ex: the expression + * fn: the function to search for + * ============================================================================ */ +get_arg_safe(ex, fn) := catch( + scanmap(lambda([u], if not(atom(u)) and op(u) = fn then throw(part(u, 1)) else u), ex), + ex +)$ + +/* ============================================================================ + FUNCTION: process_hansen_term + Description:core pattern matching unit that identifies the power of r and the + * trigonometric frequency of f within a single term replacing them + * with the corresponding hansen averaged coefficient. + Parameters: + * termo: a single summand from the hamiltonian + * vf: the true anomaly variable f + * eccs: clean symbol for eccentricity + * as: clean symbol for semi major axis + * ============================================================================ */ +process_hansen_term(termo, vf, eccS, aS) := block( + [ang, mm, pp, n_exp, coef_puro, f_fn], + + /* 1. identifies the trigonometric function and extracts the angle */ + if not freeof(cos, termo) then f_fn : cos + else if not freeof(sin, termo) then f_fn : sin + else f_fn : false, + + if f_fn # false then ( + ang : get_arg_safe(termo, f_fn), + mm : coeff(ang, vf), /* m from mathematica */ + pp : ratexpand(ang - mm*vf) /* p from mathematica where your g is */ + ) else ( + mm : 0, pp : 0 + ), + + /* 2. hunts the power of r */ + n_exp : hipow(termo, r), + if n_exp = 0 then n_exp : lopow(termo, r), + + /* 3. isolates the coefficient cc from mathematica */ + coef_puro : termo, + if f_fn # false then coef_puro : subst(1, f_fn(ang), coef_puro), + if n_exp # 0 then coef_puro : subst(1, r^n_exp, coef_puro), + + /* 4. reconstructs using hansen preserving pp */ + if n_exp # 0 then ( + return(coef_puro * aS^n_exp * (if n_exp > 0 then hansen_x0_pos[n_exp, mm, eccS] + else hansen_x0_neg(abs(n_exp)-1, mm, eccS)) * (if f_fn = cos then cos(pp) else if f_fn = sin then sin(pp) else 1)) + ) else return(0) /* term without r n in the average is garbage */ +)$ + +/* ============================================================================ + FUNCTION: apply_hansen_hughes + Description:high level iterator that applies the hansen hughes averaging process over + * a sum of terms. it linearizes the expression using trigreduce before + * mapping the process hansen term function across all summands. + Parameters: + * expr: the expression to be averaged + * varf: the true anomaly variable + * eccs: the eccentricity variable + * as: the semi major axis variable + * ============================================================================ */ +apply_hansen_hughes(expr, varF, eccS, aS) := block( + [temp, lista], + /* ensures there is no rat or cre and expands to free the terms */ + temp : totaldisrep(expr), + temp : expand(trigreduce(trigexpand(temp))), + + if atom(temp) or op(temp) # "+" then lista : [temp] else lista : args(temp), + + return(apply("+", map(lambda([u], process_hansen_term(u, varF, eccS, aS)), lista))) +)$ + +/* ============================================================================ + FUNCTION: hansen_third_body_average + Description:specialized top level function for third body potential averaging. + * it manages the jacobian weight usually r 2 a 2 eta performs + * the hansen based quadrature and splits the resulting hamiltonian + * into its secular and periodic components. + Parameters: + * u3b_raw: raw third body potential + * var_f: true anomaly variable + * jacobian_eta: jacobian weight for the integration + * ============================================================================ */ +hansen_third_body_average(U3b_raw, var_f, jacobian_eta) := block( + [U3b_jac, weighted_term, secular_part, periodic_part, e_sym, a_sym], + + /* 0. clean symbols for the hansen engine */ + e_sym : gensym("e_"), + a_sym : gensym("a_"), + + /* 1. weighting and trigonometric linearization */ + U3b_jac : U3b_raw * jacobian_eta, + + /* protects e and a during trigreduce and expand */ + weighted_term : subst([e = e_sym, a = a_sym], U3b_jac), + weighted_term : trigreduce(weighted_term), + weighted_term : expand(weighted_term), + + /* 2. hansen engine */ + secular_part : apply_hansen_hughes(weighted_term, var_f, e_sym, a_sym), + + /* 3. periodic part */ + periodic_part : expand(weighted_term - secular_part), + + /* 4. restores the original symbols */ + secular_part : subst([e_sym = e, a_sym = a], secular_part), + periodic_part : subst([e_sym = e, a_sym = a], periodic_part), + + return([secular_part, periodic_part]) +)$ + +/* ============================================================================ + FUNCTION: hansen_srp_average + Description:extracts the secular average and the periodic part of the solar radiation + * pressure srp potential exactly without eccentricity expansion. + * computational strategy keeps the symbols in keplerian form a e i to optimize the hansen + * engine and preserve compatibility with the kinematics of pb closed form. + * the conversion to delaunay actions l g h is delegated to the post processor + * sanitize hamiltonian expression. + Parameters: + * u_srp: the solar radiation pressure potential + * var_f: the true anomaly variable + * ============================================================================ */ +hansen_srp_average(U_srp, var_f) := block( + [U_exp, hamSecular, hamPeriodic], + + /* 1. trigonometric linearization in the keplerian domain + * expands and reduces the combined cosines generated in derive srp potential + * maintaining the integrity of the pure symbols. + */ + U_exp : trigreduce(U_srp), + U_exp : expand(U_exp), + + /* 2. secular extraction via hansen coefficients hughes + * since the radial power of the srp is n equals 1 the apply hansen hughes function + * via hansen x0 pos will generate exact and finite polynomials in e. + * there is no need for a jacobian here as the geometry is not altered. + */ + hamSecular : apply_hansen_hughes(U_exp, var_f, e, a), + + /* 3. periodic part extraction + * the periodic part is simply the difference between the expanded osculating + * potential and its secular average. this part will serve as input for the + * integration of the first order generating function. + */ + hamPeriodic : expand(U_exp - hamSecular), + + /* returns the list in the format that solve hori deprit recurrence + * expects to receive in f input list secular periodic + */ + return([hamSecular, hamPeriodic]) +)$ + +/* turns on confirmation to return to normal. */ +suppress_trace_warnings: false$ diff --git a/symcelmech/celmechseries.mac b/symcelmech/celmechseries.mac new file mode 100644 index 0000000..c4ab09b --- /dev/null +++ b/symcelmech/celmechseries.mac @@ -0,0 +1,316 @@ +/* celmechseries.mac */ + +suppress_trace_warnings: true$ + +/* --- global package settings --- */ + +/* declares that the variables belong to the set of real numbers */ +declare([e, E, M], real)$ + +/* adds restrictions on the range of values for integration and simplification */ +assume(e >= 0, e < 1, E >= 0, E <= 2*%pi, M >= 0, M <= 2*%pi)$ + + +/* ============================================================================ + FUNCTION: _get_e_series + Description:computes the series expansion of the eccentric anomaly (e) as a + * function of the mean anomaly (m) and eccentricity (e). + Parameters: + * order_e (int): truncation order for the eccentricity (e). + * m_symb (sym): symbol for mean anomaly. + * e_symb (sym): symbol for eccentricity. + * e_symb_ecc (sym): symbol for eccentric anomaly. + * ============================================================================ */ +_get_E_series(order_e, M_symb, e_symb, E_symb) := block([E_approx, n], + /* start with approximation e approx m */ + E_approx: M_symb, + + /* iterate to get higher order terms in e */ + for n: 1 thru order_e do ( + E_approx: M_symb + e_symb*sin(E_approx) + ), + + /* expand in taylor series */ + return(taylor(E_approx, e_symb, 0, order_e)) +)$ + +/* ============================================================================ + FUNCTION: get_cos_f_series + Description:generates the series expansion for cos(f) in terms of m and e. + Parameters: + * order_e (int): truncation order for the eccentricity. + * m_symb (sym): symbol for mean anomaly. + * e_symb (sym): symbol for eccentricity. + * e_symb_ecc (sym): symbol for eccentric anomaly. + * ============================================================================ */ +get_cos_f_series(order_e, M_symb, e_symb, E_symb) := block( + /* declare all local vars */ + [E_series, cos_f_expr, result], + + /* get series of e */ + E_series: _get_E_series(order_e, M_symb, e_symb, E_symb), + + /* expression for cos(f) */ + cos_f_expr: (cos(E_symb) - e_symb)/(1 - e_symb*cos(E_symb)), + + /* substitute e with the series and expand */ + result: subst(E_symb = E_series, cos_f_expr), + result: taylor(result, e_symb, 0, order_e), + result: trigreduce(result), /* converts powers to multiple angles */ + /* instead of appearing cos(m)^2, for example, it will appear in terms of cos(2m)... */ + return(result) +)$ + +/* ============================================================================ + FUNCTION: get_sin_f_series + Description:generates the series expansion for sin(f) in terms of m and e. + Parameters: + * order_e (int): truncation order for the eccentricity. + * m_symb (sym): symbol for mean anomaly. + * e_symb (sym): symbol for eccentricity. + * e_symb_ecc (sym): symbol for eccentric anomaly. + * ============================================================================ */ +get_sin_f_series(order_e, M_symb, e_symb, E_symb) := block( + /* declare all local vars */ + [E_series, sin_f_expr, result], + + /* get series of e */ + E_series: _get_E_series(order_e, M_symb, e_symb, E_symb), + + /* expression for sin(f) */ + sin_f_expr: sqrt(1 - e_symb^2)*sin(E_symb)/(1 - e_symb*cos(E_symb)), + + /* substitute e with the series and expand */ + result: subst(E_symb = E_series, sin_f_expr), + result: taylor(result, e_symb, 0, order_e), + result: trigreduce(result), /* converts powers to multiple angles */ + /* instead of appearing cos(m)^2, for example, it will appear in terms of cos(2m)... */ + + return(result) +)$ + +/* ============================================================================ + FUNCTION: get_r_over_a_kepler_series + Description:computes the series expansion of the normalized radial distance (r/a). + Parameters: + * order_e (int): truncation order for the eccentricity. + * m_symb (sym): symbol for mean anomaly. + * e_symb (sym): symbol for eccentricity. + * e_symb_ecc (sym): symbol for eccentric anomaly. + * ============================================================================ */ +get_r_over_a_kepler_series(order_e, M_symb, e_symb, E_symb) := block( + /* declare all local vars */ + [E_series, r_a_expr], + + /* 1. get the series for e(m, e) */ + E_series: _get_E_series(order_e, M_symb, e_symb, E_symb), + + /* 2. define the exact equation */ + r_a_expr: 1 - e_symb*cos(E_symb), + + /* 3. substitute the series of e in the equation */ + r_a_expr_subst: subst(E_symb = E_series, r_a_expr), + + /* 4. expand the final result in series of e */ + result: taylor(r_a_expr_subst, e_symb, 0, order_e), + return(trigreduce(result)) +)$ + +/* ============================================================================ + FUNCTION: get_cos_e_series + Description:computes the series expansion of cos(e). + Parameters: + * order_e (int): truncation order for the eccentricity. + * m_symb (sym): symbol for mean anomaly. + * e_symb (sym): symbol for eccentricity. + * e_symb_ecc (sym): symbol for eccentric anomaly. + * ============================================================================ */ + get_cos_E_series(order_e, M_symb, e_symb, E_symb) := block( + /* declare all local vars */ + [r_a_Series, cos_E_expr, result], + + /* 1. get the series for r/a in order (order_e+1) */ + r_a_series: get_r_over_a_kepler_series(order_e+1, M_symb, e_symb, E_symb), + + /* 2. apply the algebraic expression cos(e) = (1 - r/a) / e */ + cos_E_expr: (1 - r_a_series) / e_symb, + + /* 3. re-expand the result to the final order order_e */ + result: taylor(cos_E_expr, e_symb, 0, order_e), + + /* 4. simplify trigonometric functions (e.g., cos(m)^2 -> cos(2m)) */ + return(trigreduce(result)) +)$ + +/* ============================================================================ + FUNCTION: get_f_series + Description:computes the series for the true anomaly (f). + Parameters: + * order_e (int): truncation order for the eccentricity. + * m_symb (sym): symbol for mean anomaly. + * e_symb (sym): symbol for eccentricity. + * e_symb_ecc (sym): symbol for eccentric anomaly. + * ============================================================================ */ +get_f_series(order_e, M_symb, e_symb, E_symb) := block( + [sin_f_series, cos_f_series, f_expr, result], + + /* 1. get the series of sin(f) and cos(f) in terms of m and e */ + sin_f_series: get_sin_f_series(order_e, M_symb, e_symb, E_symb), + cos_f_series: get_cos_f_series(order_e, M_symb, e_symb, E_symb), + + /* 2. calculate f using atan2 */ + /* the series of f(m, e). the order 0 term of this is m_symb. */ + f_expr: atan2(sin_f_series, cos_f_series), + + /* 3. re-expand and simplify */ + /* we only apply taylor() to the f series to ensure + the correct truncation at order_e. */ + result: taylor(f_expr, e_symb, 0, order_e), + + return(trigreduce(result)) +)$ + +/* ============================================================================ + FUNCTION: get_f_minus_m_series + Description:computes the series for the equation of the center (f - m). + Parameters: + * order_e (int): truncation order for the eccentricity. + * m_symb (sym): symbol for mean anomaly. + * e_symb (sym): symbol for eccentricity. + * e_symb_ecc (sym): symbol for eccentric anomaly. + * ============================================================================ */ +get_f_minus_M_series(order_e, M_symb, e_symb, E_symb) := block( + [sin_f_series, cos_f_series, f_expr, f_M_expr, result], + + /* 1. get the series of sin(f) and cos(f) in terms of m and e */ + sin_f_series: get_sin_f_series(order_e, M_symb, e_symb, E_symb), + cos_f_series: get_cos_f_series(order_e, M_symb, e_symb, E_symb), + + /* 2. calculate f using atan2 */ + /* the series of f(m, e) */ + f_expr: atan2(sin_f_series, cos_f_series), + + /* 3. calculate the equation of the center (f - m) */ + /* the order 0 (e=0) of f_expr is m_symb, so they cancel out + the resulting series starts at o(e^1). */ + f_M_expr: f_expr - M_symb, + + /* 4. re-expand and simplify */ + /* we use taylor() again to ensure that the expression + (resulting from the subtraction) is correctly truncated at order_e. */ + result: taylor(f_M_expr, e_symb, 0, order_e), + + /* trigreduce() converts powers (e.g., sin(m)^2) into multiple + angles (e.g., cos(2m)), if any. */ + return(trigreduce(result)) +)$ + +/* ============================================================================ + FUNCTION: truncate_ecc_delaunay + Description:specialized truncation tool for expressions using delaunay variables (l, g). + Parameters: + * expr (expr): the expression to be truncated. + * max_ord (int): maximum power of eccentricity to preserve. + * g_var (sym): the delaunay g action symbol. + * l_var (sym): the delaunay l action symbol. + * ============================================================================ */ +truncate_ecc_delaunay(expr, max_ord, G_var, L_var) := block( + [temp, eps, result], + + /* 1. forces maxima to treat eps as positive to avoid abs */ + assume(eps > 0), + + /* 2. introduces the dummy parameter eps */ + temp: subst(G_var = L_var * sqrt(1 - eps^2), expr), + + /* 3. simplification and expansion */ + temp: taylor(temp, eps, 0, max_ord), + result: ratdisrep(temp), + + /* 4. cleanup: removes the assumption so as not to affect the rest of the system */ + forget(eps > 0), + + /* 5. returns to delaunay and expands */ + result: subst(eps = sqrt(1 - G_var^2/L_var^2), result), + return(result) +)$ + +/* ============================================================================ + FUNCTION: truncate_ecc_keplerian + Description:truncates expressions in eccentricity (e) while properly handling true anomaly dependencies. + Parameters: + * expr (expr): the expression to be truncated. + * max_ord (int): maximum power of eccentricity to preserve. + * ============================================================================ */ +truncate_ecc_keplerian(expr, max_ord) := block( + [temp], + + /* 1. Only substitute r if present (series mode should be freeof r) */ + if freeof(r, expr) then + temp: expr + else + temp: subst(r = r_two_body_exact(e, a, f), expr), + + /* 2. Taylor expansion in e */ + temp: taylor(temp, e, 0, max_ord), + + /* 3. Convert to common polynomial */ + temp: ratdisrep(temp), + + return(temp) +)$ + +/* ============================================================================ + FUNCTION: truncate_closed_form + Description:truncates expressions in closed-form by sweeping the polynomial. + Parameters: + * expr (expr): the expression to be truncated. + * max_ord (int): maximum power of eccentricity to preserve. + * ============================================================================ */ +truncate_closed_form(expr, max_ord) := block( + [res, term_list, final_expr: 0], + + /* 1. protects the exact kinematics of the closed form */ + res: subst([eta = ETA_DUMMY, sin_i = SINI_DUMMY, cos_i = COSI_DUMMY], expr), + + /* 2. forces total expansion (ratexpand) to separate the terms */ + res: ratexpand(res), + + /* 3. transforms the giant expression into an iterable list */ + if op(res) = "+" then + term_list: args(res) + else + term_list: [res], + + /* 4. sweeps and truncates */ + for term in term_list do ( + if hipow(term, e) <= max_ord then + final_expr: final_expr + term + ), + + /* 5. returns the geometric protection */ + final_expr: subst([ETA_DUMMY = eta, SINI_DUMMY = sin_i, COSI_DUMMY = cos_i], final_expr), + + /* 6. repacks the equation (drops the leaf_count) */ + final_expr: ratsimp(final_expr), + + return(final_expr) +)$ + +/* ============================================================================ + FUNCTION: shows_available_functions + Description:prints the list of available functions in this module. + Parameters: none + * ============================================================================ */ +shows_available_functions() := block( + print("available functions:"), + print(" - get_cos_f_series(order_e, M_symb, e_symb, E_symb)"), + print(" - get_sin_f_series(order_e, M_symb, e_symb, E_symb)"), + print(" - get_r_over_a_kepler_series(order_e, M_symb, e_symb, E_symb)"), + print(" - get_cos_E_series(order_e, M_symb, e_symb, E_symb)"), + print(" - get_f_series(order_e, M_symb, e_symb, E_symb)"), + print(" - get_f_minus_M_series(order_e, M_symb, e_symb, E_symb)") +)$ + +/* turns confirmation back to normal */ +suppress_trace_warnings: false$ diff --git a/symcelmech/derivepotential.mac b/symcelmech/derivepotential.mac new file mode 100644 index 0000000..dd7e52e --- /dev/null +++ b/symcelmech/derivepotential.mac @@ -0,0 +1,544 @@ +/* derivepotential.mac */ + +suppress_trace_warnings: true$ + +/* global package configurations */ + +/* adds the directory of this own file mac to the search path file search maxima. + this ensures it can find other packages such as series celestial mechanics mac + that are in the same folder regardless of where the worksheet is. */ +file_search_maxima: cons(pathname_directory(load_pathname), file_search_maxima)$ + +/* now a simple load works because maxima will look in the directory we just added. */ +load("analyticalfunctions.mac")$ /* for celmec trigexpand */ +load("celmechseries.mac")$ + +/* loads the vect package for gradient jacobian calculations */ +load("vect")$ + +/* loads orthopoly for pnm legendre p with 3 arguments */ +load("orthopoly")$ + +/* removes the added path so as not to pollute the global environment permanently. */ +file_search_maxima: rest(file_search_maxima)$ + +/* ============================================================================ + FUNCTION: derive_zonal_potential + Description:computes the potential for zonal harmonics m=0 of degree n in + * keplerian elements. matches the classical physics convention where + * acceleration a = -grad u. + Parameters: + * n: degree of the harmonic + * mu: gravitational parameter + * r_var: radius variable scalar + * f: true anomaly + * a: semi major axis + * e: eccentricity + * i: inclination + * omega: argument of periapsis + * r: equatorial radius of the central body + * compact_vars: boolean flag. if true compacts constants into a single kn term. + * ============================================================================ */ +derive_zonal_potential(n, mu, r_var, f, a, e, i, omega, R, compact_vars) := block( + [jn, sin_phi, pn, U, k_n], + + /* generates the jn symbol e.g. j2 j3 */ + jn: concat('j, n), + + /* sine of latitude sin phi = sin i * sin f + omega */ + sin_phi: sin(i) * sin(f + omega), + + /* legendre p requires the orthopoly package to be loaded */ + pn: legendre_p(n, sin_phi), + + /* generates the potential. u = mu * jn * r^n / r^n+1 * pn sin phi. + note the potential is positive here because we use the physics convention + for potential energy v = -u astro. since jn = -cn0 the signs cancel out */ + U: mu * jn * (R^n / r_var^(n + 1)) * pn, + + U: ratsubst(cos_i, cos(i), U), + U: ratsubst(sin_i, sin(i), U), + + /* compacting physical constants */ + if compact_vars then ( + /* generates the kn name e.g. k2 k3 */ + k_n: concat('k_j, n), + + /* substitutes in the potential */ + U: ratsubst(k_n, jn * R^n * mu, U) + ), + + return(U) +)$ + +/* defines r expr algebraic expression and r simplification symbol. + these are global so the functions can use them. */ +r_expr: sqrt(x^2 + y^2 + z^2)$ +r2_expr: x^2 + y^2 + z^2$ + +/* ============================================================================ + FUNCTION: derive_cartesian_zonal_harmonic + Description:computes the potential and acceleration for zonal harmonics m=0 + * of degree n in cartesian coordinates. + Parameters: + * n: degree of the harmonic + * sep_coord: boolean flag. if true returns acceleration components ax ay az. + * if false returns the scalar potential u. + * ============================================================================ */ +derive_cartesian_zonal_harmonic(n, sep_coord) := block( + [ + jn, sin_theta, Pn_cart, U, U_expanded, U_simpl, + U_grad, ax_simpl, ay_simpl, az_simpl + ], + + /* defines the jn symbol e.g. j2 j3 */ + jn: concat('j, n), + + sin_theta: z / r_expr, + + /* gets the legendre polynomial e.g. 3*u^2-1 / 2 for n=2 */ + Pn_cart: legendre_p(n, sin_theta), + + /* generates the potential. u = mu * jn * r^n / r^n+1 * pn sin phi. + note the potential is positive here because we use the physics convention + for potential energy v = -u astro. since jn = -cn0 the signs cancel out */ + U: (mu * R^n * jn) / (r_expr^(n+1)) * Pn_cart, + + /* converts u into a purely algebraic expression */ + /* this will expand r expr and simplify the polynomial */ + U_expanded: trigexpand(U), + + U_simpl: ratsubst(r, r_expr, U_expanded), + + if sep_coord then ( + + /* calculates the gradient now over a purely algebraic expression */ + U_grad: jacobian([U_expanded], [x, y, z])[1], + + /* remembering that a = -grad u */ + ax_simpl: ratsubst(r, r_expr, -U_grad[1]), + ay_simpl: ratsubst(r, r_expr, -U_grad[2]), + az_simpl: ratsubst(r, r_expr, -U_grad[3]), + + /* returns the factored equations */ + return ([ + subst((x^2+y^2+z^2)=r^2, factor(ax_simpl)), + subst((x^2+y^2+z^2)=r^2, factor(ay_simpl)), + subst((x^2+y^2+z^2)=r^2, factor(az_simpl)) + ]) + ) else ( + return(factor(U_simpl)) + ) +)$ + + +/* ============================================================================ + FUNCTION: derive_tesseral_potential + Description:derives the tesseral potential cnm in terms of the true anomaly. + Parameters: + * n: degree of the harmonic + * m: order of the harmonic + * mu: gravitational parameter + * r: radius variable scalar + * f: true anomaly + * a: semi major axis + * e: eccentricity + * i: inclination + * omega: argument of periapsis + * omega: longitude of the ascending node + * r: equatorial radius of the central body + * compact_vars: boolean flag. if true compacts constants into a single kc term. + * frame: boolean flag. if true uses inertial frame. if false uses rotating frame. + * include_snm: boolean flag. if true includes snm terms. + * ============================================================================ */ +derive_tesseral_potential(n, m, mu, r, f, a, e, i, omega, Omega, R, compact_vars, frame, include_snm) := block( + [u, X_bf, Y_bf, Z_bf, p_nm_term, multipoles, z_var, u_pot, c_nm, s_nm, h_rel], + + /* frame = true is inertial and frame = false is rotating */ + + /* 1. definition of the longitude angle relative to the body */ + /* in the inertial frame we subtract the rotation h - theta. in the rotating + frame h is already the relative coordinate. theta = w rot * t */ + h_rel : if frame then h - theta else h, + + /* 2. coefficients cnm and snm are constant in the body frame */ + c_nm: concat('C, n, m), + s_nm: concat('S, n, m), + + /* zeros snm before heavy algebra begins */ + if not include_snm then ( + s_nm: 0 + ), + + /* 3. orbital geometry */ + u: f + omega, + + /* 4. projection in body fixed directions */ + X_bf: cos(u)*cos(h_rel) - sin(u)*cos(i)*sin(h_rel), + Y_bf: cos(u)*sin(h_rel) + sin(u)*cos(i)*cos(h_rel), + Z_bf: sin(u)*sin(i), + + /* 5. latitude term legendre derivative */ + p_nm_term: diff(legendre_p(n, z_var), z_var, m), + p_nm_term: subst(z_var = Z_bf, p_nm_term), + + /* 6. multipole expansion via complex newton binomial */ + /* equivalent to complexexpand re c - is * x + iy ^ m in mathematica. + this generates the cos m lambda and sin m lambda terms in a polynomial form. */ + multipoles: realpart(expand((c_nm - %i*s_nm) * (X_bf + %i*Y_bf)^m)), + + /* 7. assembly of the hamiltonian */ + u_pot: -mu * (R^n / r^(n + 1)) * p_nm_term * multipoles, + + /* 8. preparation for the lie deprit solver */ + u_pot: ratsubst(cos_i, cos(i), u_pot), + u_pot: ratsubst(sin_i, sin(i), u_pot), + + /* compaction of physical constants */ + if compact_vars then ( + /* generates the name k cnm e.g. k c22 k c31 */ + k_c: concat('k_c, n, m), + u_pot: ratsubst(k_c, c_nm * R^n * mu, u_pot), + + /* only substitutes s nm if the include snm flag is true */ + if include_snm then ( + /* generates the name k snm e.g. k s22 k s31 */ + k_s: concat('k_s, n, m), + u_pot: ratsubst(k_s, s_nm * R^n * mu, u_pot) + ), + + u_pot: factor(u_pot) + ), + + return(u_pot) +)$ + +/* ============================================================================ + FUNCTION: derive_cartesian_inertial_tesseral_sectorial + Description:computes the gravitational potential or acceleration gradient for + * tesseral and sectorial harmonics in a cartesian inertial frame. + Parameters: + * n: degree of the harmonic + * m: order of the harmonic + * t: time variable + * omega_rot: rotation rate of the body mean motion + * sep_coord: boolean flag. if true returns acceleration components ax ay az. + * if false returns the scalar potential u. + * sep_cnm_snm: boolean flag. if true separates the output into cnm and snm parts. + * ============================================================================ */ +derive_cartesian_inertial_tesseral_sectorial(n, m, t, omega_rot, sep_coord, sep_cnm_snm) := block( + [ + cnm, snm, sin_theta, + lambda_inercial, theta_body, lambda_rot, + Pnm_u, Pnm_cart, U, U_expanded, + U_grad, U_simpl, + ax_simpl, ay_simpl, az_simpl, + ax_cnm_simpl, ay_cnm_simpl, az_cnm_simpl, + ax_snm_simpl, ay_snm_simpl, az_snm_simpl, + U_cnm, U_snm + ], + + /* 1. define cnm and snm symbols */ + cnm: concat('c, n, m), + snm: concat('s, n, m), + + /* 2. define inertial spherical coordinates */ + sin_theta: z / r_expr, + lambda_inercial: atan2(y, x), + theta_body: omega_rot * t, /* omega rot is the mean motion n */ + lambda_rot: lambda_inercial - theta_body, /* this is the longitude in the rotating frame */ + + /* 3. define the potential u in symbolic form */ + Pnm_u: assoc_legendre_p(n, m, u), + Pnm_cart: ratsubst(sin_theta, u, Pnm_u), + + /* 3c. build the formula for u using the rotating longitude lambda rot */ + U: -(mu * R^n) / (r_expr^(n+1)) * Pnm_cart * (cnm * cos(m*lambda_rot) + snm * sin(m*lambda_rot)), + + /* 4. expand u + this will now generate explicit terms of sin omega rot t and cos omega rot t */ + U_expanded: trigexpand(U), + + /* 5. calculate the gradient grad u in inertial coordinates x y z */ + U_grad: jacobian([U_expanded], [x, y, z])[1], + + /* 6. substitute x^2+y^2+z^2 with r^2 for a clean output */ + U_simpl: ratsubst(r, r_expr, U_expanded), + + if sep_coord then ( + + if sep_cnm_snm then ( + /* we extract components first to apply ratcoef safely to scalars */ + ax_cnm_simpl: ratsubst(r, r_expr, cnm * ratcoef(-U_grad[1], cnm)), + ay_cnm_simpl: ratsubst(r, r_expr, cnm * ratcoef(-U_grad[2], cnm)), + az_cnm_simpl: ratsubst(r, r_expr, cnm * ratcoef(-U_grad[3], cnm)), + + ax_snm_simpl: ratsubst(r, r_expr, snm * ratcoef(-U_grad[1], snm)), + ay_snm_simpl: ratsubst(r, r_expr, snm * ratcoef(-U_grad[2], snm)), + az_snm_simpl: ratsubst(r, r_expr, snm * ratcoef(-U_grad[3], snm)), + + /* 7. return the factored equations */ + /* i used trigexpand again because in step 4 + maxima cannot expand atan2 when m > 1 + leaving atan2 explicit in the equations */ + return ([ + factor(trigexpand(ax_cnm_simpl)), + factor(trigexpand(ay_cnm_simpl)), + factor(trigexpand(az_cnm_simpl)), + factor(trigexpand(ax_snm_simpl)), + factor(trigexpand(ay_snm_simpl)), + factor(trigexpand(az_snm_simpl)) + ]) + + ) else ( + ax_simpl: ratsubst(r, r_expr, -U_grad[1]), + ay_simpl: ratsubst(r, r_expr, -U_grad[2]), + az_simpl: ratsubst(r, r_expr, -U_grad[3]), + + /* 7. return the factored equations */ + /* i used trigexpand again because in step 4 + maxima cannot expand atan2 when m > 1 + leaving atan2 explicit in the equations */ + return ([ + factor(trigexpand(ax_simpl)), + factor(trigexpand(ay_simpl)), + factor(trigexpand(az_simpl)) + ]) + ) + + ) else ( + if sep_cnm_snm then ( + U_cnm: cnm * ratcoef(U_simpl, cnm), + U_snm: snm * ratcoef(U_simpl, snm), + return([factor(trigexpand(U_cnm)), factor(trigexpand(U_snm))]) + ) else ( + return(factor(trigexpand(U_simpl))) + ) + ) +)$ + + +/* ============================================================================ + FUNCTION: get_inertial_vector + Description:computes the inertial position vector given orbital elements. + Parameters: + * f: true anomaly + * omega: argument of periapsis + * omega: longitude of the ascending node + * i: inclination + * ============================================================================ */ +get_inertial_vector(f, omega, Omega, i) := block( + [u, cu, su, cOmega, sOmega, ci, si], + u: f + omega, + cu: cos(u), su: sin(u), + cOmega: cos(Omega), sOmega: sin(Omega), + ci: cos(i), si: sin(i), + + [ + cOmega * cu - sOmega * su * ci, + sOmega * cu + cOmega * su * ci, + su * si + ] +)$ + + +/* ============================================================================ + FUNCTION: get_general_third_body_factors + Description:computes the geometric factors for the third body perturbation. + Parameters: + * n: degree of the harmonic + * mupert: gravitational parameter of the perturbing body + * rpert: radial distance of the perturbing body + * fpert: true anomaly of the perturbing body + * omegapert: argument of periapsis of the perturbing body + * omegapert: longitude of the ascending node of the perturbing body + * ipert: inclination of the perturbing body + * rsat: radial distance of the satellite + * fsat: true anomaly of the satellite + * omegasat: argument of periapsis of the satellite + * omegasat: longitude of the ascending node of the satellite + * isat: inclination of the satellite + * ============================================================================ */ +get_general_third_body_factors(n, muPert, rPert, fPert, omegaPert, OmegaPert, iPert, + rSat, fSat, omegaSat, OmegaSat, iSat) := block( + [posSat, posPert, cosPsi, Pn], + + /* 1. unit vectors */ + posSat: get_inertial_vector(fSat, omegaSat, OmegaSat, iSat), + posPert: get_inertial_vector(fPert, omegaPert, OmegaPert, iPert), + + /* 2. dot product to find cos psi */ + /* in maxima the . operator performs the dot product of lists */ + cosPsi: posSat . posPert, + + /* keeps in sum of cosines essential for analytical integration */ + cosPsi: trigreduce(cosPsi), + + /* legendre p comes from the orthopoly package loaded in the header */ + Pn: legendre_p(n, cosPsi), + + [ + muPert * rPert^-(n + 1), /* c factor */ + rSat^n, /* r factor */ + Pn /* p factor */ + ] +)$ + + +/* ============================================================================ + FUNCTION: derive_general_third_body + Description:computes the complete third body perturbing potential. + Parameters: + * orderecc: order of eccentricity not directly used here but kept for signature + * n: degree of the harmonic + * mu3: gravitational parameter of the perturbing body + * r3: radial distance of the perturbing body + * f3: true anomaly of the perturbing body + * a3: semi major axis of the perturbing body + * e3: eccentricity of the perturbing body + * i3: inclination of the perturbing body + * omega3: argument of periapsis of the perturbing body + * omega3: longitude of the ascending node of the perturbing body + * r: radial distance of the satellite + * f: true anomaly of the satellite + * a: semi major axis of the satellite + * e: eccentricity of the satellite + * i: inclination of the satellite + * omega: argument of periapsis of the satellite + * omega: longitude of the ascending node of the satellite + * ============================================================================ */ +derive_general_third_body(orderEcc, n, mu3, r3, f3, a3, e3, i3, omega3, Omega3, r, f, a, e, i, omega, Omega) := block( + [factors, Cfac, rFac, Pfac, U], + + print("1. calculating vector geometry symbolic in f and r..."), + + factors: get_general_third_body_factors( + n, mu3, r3, f3, omega3, Omega3, i3, + r, f, omega, Omega, i + ), + + Cfac: factors[1], + rFac: factors[2], + Pfac: factors[3], + + U: -Cfac * rFac * Pfac, + + print("2. process completed."), + return(U) +)$ + +/* ============================================================================ + FUNCTION: derive_srp_potential + Description:computes the solar radiation pressure srp perturbing potential. + Parameters: + * r: radial distance of the satellite + * a: semi major axis of the satellite + * e: eccentricity of the satellite + * i: inclination of the satellite + * f: true anomaly of the satellite + * g: argument of periapsis of the satellite + * h: longitude of the ascending node of the satellite + * musun: gravitational parameter of the sun + * beta: srp coefficient + * rsun: radial distance of the sun + * asun: semi major axis of the sun + * esun: eccentricity of the sun + * isun: inclination of the sun + * fsun: true anomaly of the sun + * gsun: argument of periapsis of the sun + * hsun: longitude of the ascending node of the sun + * ============================================================================ */ +derive_srp_potential(r, a, e, i, f, g, h, muSun, beta, rSun, aSun, eSun, iSun, fSun, gSun, hSun) := block( + [COSI, SENI, A, B, C, D, EE, cosPsi, U], + + /* orbital geometry */ + COSI : cos(i), + SENI : sin(i), + + A : 1/4 * (1 + COSI) * (1 - cos(iSun)), + B : 1/4 * (1 - COSI) * (1 + cos(iSun)), + C : 1/4 * (1 + COSI) * (1 + cos(iSun)), + D : 1/4 * (1 - COSI) * (1 - cos(iSun)), + EE : 1/2 * SENI * sin(iSun), + + print("1. assembling cos psi..."), + + /* exact cos psi */ + cosPsi : + C * cos(f + g + h - fSun - gSun - hSun) + + A * cos(f + g + h + fSun + gSun - hSun) + + B * cos(f + g - h + fSun + gSun + hSun) + + D * cos(f + g - h - fSun - gSun + hSun) + + EE * (cos(f + g - fSun - gSun) - cos(f + g + fSun + gSun)), + + print("2. assembling the potential..."), + /* exact srp potential */ + U : beta * muSun * (r / rSun^2) * cosPsi, + + U: ratsubst(cos_i, cos(i), U), + U: ratsubst(sin_i, sin(i), U), + + + print("3. process completed."), + return(U) +)$ + +/* ============================================================================ + FUNCTION: prepare_potential_poisson + Description:prepares the potential for poisson series expansion. + Parameters: + * expr: original potential expression + * order_e: truncation order for eccentricity + * m_symb: symbol for mean anomaly + * e_symb: symbol for eccentricity + * e_symb: symbol for eccentric anomaly + * f_symb: symbol for true anomaly + * r_symb: symbol for radial distance + * r_expr: expression to substitute for r_symb + * ============================================================================ */ +prepare_potential_poisson(expr, order_e, M_symb, e_symb, E_symb, f_symb, r_symb, r_expr) := block( + [U, temp, cos_f, sin_f], + + U: subst(r_symb=r_expr, expr), + + cos_f: get_cos_f_series(order_e, M_symb, e_symb, E_symb), + sin_f: get_sin_f_series(order_e, M_symb, e_symb, E_symb), + + /* 1. replaces r and f with the series in l and e */ + temp: subst([cos(f_symb)=cos_f, sin(f_symb)=sin_f], celmec_trigexpand(U, f_symb)), + + /* 2. taylor series to truncate eccentricity */ + temp: ratdisrep(taylor(temp, e_symb, 0, order_e)), + + /* 3. trigonometric linearization */ + temp: trigreduce(temp), + + /* 4. final organization */ + return(ratexpand(temp)) +)$ + +/* ============================================================================ + FUNCTION: prepare_potential_closedform + Description:prepares a perturbing potential for closed form orbital averaging. + * the process involves substituting the radial distance r with its conic + * expression multiplying by the jacobian of the transformation one dimensional + * jacobian matrix and applying algebraic and trigonometric simplifications. + * the output is a linearized expansion that facilitates the separation of + * secular and periodic terms. + Parameters: + * expr: the original perturbing potential expression. + * jac: the jacobian of the transformation e.g. dl df. + * r_sym: symbol for radial distance. + * r_expr: the analytical expression for r radial distance. + * ============================================================================ */ +prepare_potential_closedform(expr, jac, r_sym, r_expr) := block( + [expr_r, expr_simp, expr_trig, expr_int, ci_s, si_s], + + ci_s : gensym("ci_"), si_s : gensym("si_"), + expr_r : subst([r_sym = r_expr, cos_i = ci_s, sin_i = si_s], expr * jac), + expr_simp : ratsimp(expr_r), + expr_trig : trigreduce(expr_simp), + expr_int : ratexpand(expr_trig), + + return(subst([ci_s = cos_i, si_s = sin_i], expr_int)) +)$ diff --git a/symcelmech/lietransformations.mac b/symcelmech/lietransformations.mac new file mode 100644 index 0000000..1556f0e --- /dev/null +++ b/symcelmech/lietransformations.mac @@ -0,0 +1,424 @@ + +/* lietransformations.mac */ + +suppress_trace_warnings: true$ +redefine_warn: false$ + +/* captures the directory where this file (.mac) resides. */ +file_search_maxima: cons(pathname_directory(load_pathname), file_search_maxima)$ + +/* loads modules that will be used */ +load("analyticalfunctions.mac")$ +load("average.mac")$ +load("celmechseries.mac")$ +load("algebraicutils.mac")$ + + +/* ============================================================================ + * FUNCTION: compute_B + * ============================================================================ + * Description: + * Computes the intermediate terms (B_{n,k}) of the recursive Lie triangle + * according to the Hori-Deprit/Mersman formalism. This function handles the + * cross-coupling between previous Hamiltonian terms and generating functions + * using Poisson Brackets. + * + * Mathematical Reference: + * Implements Equation 11 and Equation 12 from Mersman (1970), "A Unified + * Treatment of the Theory of Hamiltonian Systems and Lie Series." + * + * Parameters: + * nn (int): Current row of the Lie triangle (order of perturbation). + * kk (int): Current column of the Lie triangle. + * vars_list (list): Canonical variables for the Poisson Brackets. + * max_order_s (int): Maximum order allowed for generating functions. + * closed_form (bool): If true, performs exact algebraic simplification. + * If false, operates in series expansion mode. + * order_trunc_e (int/bool): Order of eccentricity truncation (if in series mode). + * + * Returns: + * The computed symbolic expression for the term B[nn, kk], properly + * sanitized and simplified according to the integration mode. + * + * Internal Logic: + * - k=0: Boundary condition. Returns the initial transformed term (f_star[nn]). + * - k>0: Evaluates the recursive summation of Poisson Brackets: + * B[n, k] = (1/k) * sum_{j=1}^{n-k+1} {B[n-j, k-1], S[j]} + * - Employs 'pb_memo' for high-speed Poisson Bracket retrieval. + * - Handles 'trigreduce' and 'rat' conversion to maintain series structure + * when in expansion mode. + * ============================================================================ */ +compute_B(nn, kk, vars_list, max_order_s, closed_form, order_trunc_e) := block( + + [jj, sum_val: 0, temp_list, reduced_list, temp_F, temp_S, new_term, temp_args], + + if kk = 0 then return(f_star[nn]), /* boundary condition for column zero (equation 12, Mersman, 1970) */ + + for jj:1 thru min(nn - kk + 1, max_order_s) do ( + if closed_form then ( + new_term: pb_memo(B_cache[nn-jj, kk-1], s_gen[jj]), + sum_val: clean_pb(sum_val + new_term) + + ) else ( + /* recurrence for intermediate terms (equation 11, Mersman, 1970) */ + /* called via memoization: ignores vars_list because it's already fixed in the cache. */ + new_term : pb_memo(B_cache[nn-jj, kk-1], s_gen[jj]), + + if order_trunc_e # false then ( + new_term : truncate_ecc_keplerian(new_term, order_trunc_e) + ), + + sum_val : sum_val + new_term + ) + ), /* end of for loop */ + + + if closed_form then ( + new_term : clean_pb(sum_val / kk) + ) else ( + /* for expansion in e, the rat() function helps maintain the series structure. */ + /*new_term : trigreduce(rat(sum_val / kk))*/ + new_term : ratexpand(sum_val / kk), + if new_term = 0 then + new_term : 0 + else if atom(new_term) or not(op(new_term) = "+") then + new_term : trigreduce(new_term) + else + new_term : map(trigreduce, new_term) + ), + + return(sanitize_hamiltonian_expression(new_term, closed_form, order_trunc_e)) +)$ + + +/* ============================================================================ + * FUNCTION: solve_hori_deprit_recurrence + * ============================================================================ + * Description: + * Core analytical perturbation engine based on the Lie Transformation Method + * (Hori-Deprit/Mersman formalism). It constructs the recursive Lie triangle + * to compute transformed Hamiltonians (secular/mean) and the associated + * generating functions of the system up to an arbitrary order. + * + * Integration Architectures (Dual-Track): + * 1. Closed-Form (closed_form = true): Performs exact integration in the true + * anomaly domain by substituting the orbit equation. It natively incorporates + * Deprit's Parallax Elimination by appending the Equation of the Center + * (f - l) to the first-order generating function (S1), effectively + * shielding the system from singularities at e = 0. + * 2. Series Expansion (closed_form = false): Applies standard orbital averaging + * theory in Delaunay variables (averaging over mean anomaly 'l'). Ideal for + * handling pre-processed external perturbations, such as Third-Body + * effects expanded via Hansen coefficients. + * + * Parameters: + * max_order (int): Maximum order of the Lie expansion (e.g., 2, 3, 4). + * f_input_list (list): List of original Hamiltonians [F0, F1, F2...]. + * Supports sub-lists [secular, periodic] if a term + * has been pre-processed. + * vars_list (list): List of canonical variables (e.g., [l, g, h, L, G, H]). + * avg_var (sym): Integration variable for averaging (e.g., 'f' or 'l'). + * n_avg (sym): Fundamental frequency (e.g., mean motion 'n_mean'). + * max_order_s (int): Maximum order for the explicit calculation of + * generating functions. + * closed_form (bool): If true, solves in closed form (exact in e). + * If false, uses Taylor/Poisson series expansion. + * [opt_args] (opts): 1. order_e: Truncation order for eccentricity. + * 2. trunc_gen: (bool) Truncate generating functions? + * 3. order_trunc_gen: Truncation order for gen. functions. + * + * Returns: + * A list containing two arrays: [list(f_star), list(s_gen)] + * [1] f_star: The transformed Hamiltonians F_n* (Secular/Mean Energies). + * [2] s_gen: The Generating Functions S_n. + * + * Architectural Notes: + * - Employs aggressive memory management (Garbage Collection via 'gc()' and + * 'kill(labels)') to prevent Lisp Heap exhaustion during the expansion of + * complex symbolic trees in higher orders (>= 2). + * - In Closed-Form (Order 1), the function automatically restores time-geometry + * canonicity by introducing the Equation of the Center (f - l) into S1, + * consistent with classical literature (Deprit, 1981). + * - Designed for "Cascading Transformations": the output of one call can be + * seamlessly injected into another to perform successive averages (e.g., + * eliminating 'f', then 'f_3', then 'h'). + * ============================================================================ */ +solve_hori_deprit_recurrence(max_order, f_input_list, vars_list, avg_var, n_avg, max_order_s, closed_form, [opt_args]) := block( + [n, k, j, entry_term, g_brackets, parts, periodic_total, temp_to_avg, n_term_f, n_term_s, + order_e, trunc_gen, f_phys], + + /* + [opt_args], because Maxima only allows one optional argument. To use more than one, a single + name is used to capture all the extra arguments and then distribute them within the block. + */ + + /* extracting optional arguments from the 'opt_args' list */ + order_e: if length(opt_args) >= 1 then opt_args[1] else false, + trunc_gen: if length(opt_args) >= 2 then opt_args[2] else false, + order_trunc_gen: if length(opt_args) >= 3 then opt_args[3] else false, + + /* clear local/global cache */ + kill(f_star, s_gen, B_cache, f_phys), + + /* prepare the dependencies ONLY if it is not a closed (or global) format. */ + if not(closed_form) then ( + setup_canonical_dependencies() + ), + + /* calls analyticalutils memoization */ + setup_pb_cache(vars_list, closed_form), + + /* initial condition: order 0 (H0) */ + f_star[0]: if listp(f_input_list[1]) then f_input_list[1][1] else f_input_list[1], + f_phys[0]: f_star[0], + B_cache[0, 0]: f_star[0], /* definition of the basis of the recursive triangle (Equation 12, Mersman, 1970) */ + + if closed_form then ( + r_expr: r_two_body_exact(e,a,f), + r_expr_eta: ratsubst(eta^2, 1-e^2, r_expr) + ), + + print("=== Starting Hori-Deprit Method ==="), + if closed_form then ( + print(" --- Closed Form ---"), + print(" - Average over:", avg_var) + + ) else ( + print("--- Series Expansion ---"), + print(" - Ecc Truncation order:", order_e), + print(" - Gen. function truncation order:", order_trunc_gen) + ), + + for n: 1 thru max_order do ( + print("Processing order ", n, "..."), + + /* collection of the entry term (Fn) */ + entry_term: f_input_list[n+1], + + /* accumulation of brackets (Lie's series contribution) */ + g_brackets: 0, + + /* contribution k >= 2 (higher-order curvature) */ + for k: 2 thru n do ( + B_cache[n, k]: compute_B(n, k, vars_list, max_order_s, closed_form, order_e), + g_brackets: g_brackets - B_cache[n, k] /* accumulation of the Lie series to compose the nth order (equation 13, Mersman, 1970) */ + ), + + /* contribution k = 1 (interaction between Hamiltonians and previous generatrices) */ + if n > 1 then ( + for j: 1 thru min(n - 1, max_order_s) do ( + /* for closed_form, reconstruct the Lie base with the original physicist. */ + ham_lie : f_star[n-j], + g_brackets: g_brackets - pb_memo(ham_lie, s_gen[j]) + ), + + /* expands the total sum of brackets only once */ + g_brackets: ratexpand(g_brackets) /* ratexpand is necessary to avoid zeroing Fn*, n > 1 */ + ), + + /* average and secular/periodic separation */ + block([entrada_sec, entrada_per], + /* collects the input term (Fn) - checks if the user passed [sec, per] or a single term */ + entrada_sec : if listp(entry_term) then entry_term[1] else entry_term, + entrada_per : if listp(entry_term) then entry_term[2] else 0, + + /* g_brackets already contain the Lie couplings calculated in compute_B */ + temp_to_avg: entrada_sec + g_brackets, + + /* --------------------------------------------------------------------------------- + * AUTO-BYPASS OF PRE-SECULAR POTENTIALS (Ex: Third Body via Hansen) + * If 'temp_to_avg' is already secularized with respect to 'r' (does not possess the primary geometric variable), the prepare_potential_closedform routine will not execute + * no replacement and will act only as a safe pass-through, maintaining the integrity of terms from distant disturbing bodies (r_3) + * --------------------------------------------------------------------------------- */ + + temp_to_avg: prepare_potential_closedform(temp_to_avg, 1, r, r_expr_eta), + /* average: works for series (l) or closed form (f)*/ + parts : map_average(temp_to_avg, avg_var), + + f_star[n] : parts[1], + f_star[n] : sanitize_hamiltonian_expression(f_star[n], closed_form, order_e), + + /* the total periodic part is the sum of the original perturbation plus Lie couplings */ + periodic_total: entrada_per + parts[2], + + /* generation of s_n (generating function) */ + if n <= max_order_s then ( + /* if there is no periodic part, the generator is zero */ + if periodic_total = 0 then ( + s_gen[n]: 0 + + ) else ( + /* quadrature of the periodic part */ + s_gen[n]: -generating_function_quadrature(periodic_total, avg_var) / n_avg, + + /*s_gen[n]: clean_atan2(s_gen[n]),*/ + s_gen[n]: sanitize_hamiltonian_expression(s_gen[n], closed_form, order_trunc_gen), + + /* eccentricity truncation */ + if trunc_gen then ( + if closed_form then + s_gen[n]: truncate_closed_form(s_gen[n], order_trunc_gen) + else + s_gen[n]: truncate_ecc_keplerian(s_gen[n], order_trunc_gen) + ) + ) + ) else ( + s_gen[n]: 0 + ) + ), + + /* Apply variable mapping once per order (series mode only) */ + if not(closed_form) then ( + f_star[n]: subst([ + cos(i) = cos_i, + sin(i)^2 = (1 - cos_i^2), + a = L^2/mu, + H = G*cos_i + ], f_star[n]), + f_star[n]: ratsimp(f_star[n]), + + if n <= max_order_s then + s_gen[n]: subst([ + cos(i) = cos_i, + sin(i)^2 = (1 - cos_i^2), + a = L^2/mu, + H = G*cos_i + ], s_gen[n]) + ), + + /* f_star statistics */ + n_term_f : nterms(f_star[n]), + print(" > Statistics f_star[", n, "]: ", n_term_f, " terms"), + + B_cache[n, 0]: f_star[n], /* updating the zero column with the transformed Hamiltonian (Equation 12, Mersman, 1970) */ + + /* s_gen statistics */ + n_term_s : nterms(s_gen[n]), + print(" > Statistics s_gen[", n, "]: ", n_term_s, " terms"), + + /* cache for B(n,1) */ + if n < max_order then ( + temp_sum: 0, + for j: 1 thru min(n, max_order_s) do ( + ham_lie : f_star[n-j], + temp_sum: temp_sum + pb_memo(ham_lie, s_gen[j]), + + B_cache[n, 1]: sanitize_hamiltonian_expression(temp_sum, closed_form, order_e) + ) + ), + + print(" -> f_star[", n, "] ok"), + + /* aggressive temporary memory cleanup */ + kill(labels), /* remove the history of %i and %o that Maxima keeps */ + if n > 1 then gc() /* activates the Lisp Garbage Collector */ + + /* + kill(labels): Maxima stores each intermediate step in memory so you can call them with %. In a 3rd order loop, + this stores thousands of useless expressions. + + gc(): Maxima runs on Lisp. Lisp doesn't always clear memory immediately when a variable is discarded. The gc() function forces this + cleanup, which usually "empties" the RAM usage you see in the task manager. + */ + ), + + if not(closed_form) then ( + teardown_canonical_dependencies(), + + /* Post-processing: eliminate G from secular Hamiltonians + to remove spurious 1/e^2 singularities */ + for n: 1 thru max_order do ( + f_star[n]: subst(G = L*sqrt(1-e^2), f_star[n]), + f_star[n]: taylor(f_star[n], e, 0, order_e), + f_star[n]: ratdisrep(f_star[n]), + f_star[n]: ratsimp(f_star[n]), + s_gen[n]: trigsimp(s_gen[n]), + s_gen[n]: ratsimp(s_gen[n]) + ) + ), + + print("=== Algorithm completed ==="), + return([listarray(f_star), listarray(s_gen)]) +)$ + + +/* ============================================================================ + * FUNCTION: apply_lie_corrections + * ============================================================================ + * Description: + * Computes the periodic corrections (short and long-period variations) for + * orbital elements using the Lie series operator. This function maps mean + * elements to osculating elements (direct transform) or vice-versa (inverse + * transform) by evaluating the Lie series displacement. + * + * Parameters: + * mean_element (expr/list): The symbolic element(s) to be transformed + * (e.g., a, e, i, g, h, l). + * s_gen_list (list): The list of generating functions [S1, S2, ...] + * produced by the Hori-Deprit solver. + * vars_list (list): Canonical variables for the Poisson Brackets. + * closed_form (bool): If true, uses the closed-form Poisson engine. + * [inverse_opt] (bool): Optional. If true, performs the inverse + * transformation (osculating -> mean) by negating the + * generating functions. + * + * Returns: + * The total periodic displacement (Delta) for the given element(s). + * The osculating element is then X_osc = X_mean + Delta. + * + * Internal Logic: + * - Implements the Lie series recurrence: + * X_n = (1/n) * sum_{j=1}^{n} {X_{n-j}, S_j} + * - Starts with X_0 as the identity term. + * - Sums terms from order 1 to N to obtain the total correction. + * - Supports mapping over lists of elements for batch processing. + * ============================================================================ */ +apply_lie_corrections(mean_element, s_gen_list, vars_list, closed_form, [inverse_opt]) := block( + [internal_transform, is_inverse], + + is_inverse: if length(inverse_opt) > 0 then inverse_opt[1] else false, + + internal_transform(expr) := block( + [n, k, max_order, s_list, x_triangle, total_sum], + local(x_triangle), + max_order: length(s_gen_list), + + s_list: if is_inverse then map("-", s_gen_list) else s_gen_list, + + /* Initializes column zero (k=0) of the triangle for variables */ + x_triangle[0,0] : expr, + for n:1 thru max_order do ( + x_triangle[n,0] : 0 + ), + + /* Populate the Lie Triangle using Mersman/Deprit recursion */ + for n:1 thru max_order do ( + for k:1 thru n do ( + x_triangle[n,k]: 0, + for j:1 thru min(n - k + 1, max_order) do ( + /* Direct use of the Poisson parenthesis. Evaluate whether using compute_B has any advantages in the future. */ + x_triangle[n,k]: x_triangle[n,k] + poisson_bracket(x_triangle[n-j, k-1], s_list[j], vars_list, closed_form) + ), + x_triangle[n,k]: x_triangle[n,k] / k, + + if closed_form then + x_triangle[n,k] /*: ratsimp(x_triangle[n,k])*/ + ) + ), + + /* The total correction is the sum of all the elements of the triangle */ + total_sum: 0, + for n:1 thru max_order do ( + for k:1 thru n do ( + total_sum: total_sum + x_triangle[n,k] + ) + ), + return(total_sum) + ), + + if listp(mean_element) then + map(internal_transform, mean_element) + else + internal_transform(mean_element) +)$ + diff --git a/symcelmech/planetaryequations.mac b/symcelmech/planetaryequations.mac new file mode 100644 index 0000000..c4615e0 --- /dev/null +++ b/symcelmech/planetaryequations.mac @@ -0,0 +1,497 @@ +/* planetaryequations.mac */ + +suppress_trace_warnings: true$ + +kill(a, e, i, f, l, g, h, L, G, H)$ + +/*--- global package configurations --- */ + +/* -------------------------------------------------------------------------- + * implementation note: global variables vs. _symb arguments + * -------------------------------------------------------------------------- + * + * the 'declare' and 'assume' lines above define the properties + * (e.g., 'real', 'a>0') for the *standard global symbols* ('a', 'e', 'i', etc.). + * + * the functions of the package (like 'lagrange_planetary_equations') use + * the '_symb' convention (e.g., 'a_symb', 'e_symb') in their *arguments* + * to be modular and work with any symbol. + * + * when to use: + * when calling the function with standard symbols, like: + * >> lagrange_planetary_equations(r_j2, a, e, i, omega, omega, m, n) + * + * maxima passes the symbol 'a' (which it knows is > 0) to the argument 'a_symb'. + * + * this global configuration is correct and should not be changed to '_symb'. + * + * -------------------------------------------------------------------------- */ + +/* declares that variables belong to the set of real numbers. */ + +declare([a, e, i, Omega, omega, M, n], real)$ + +/* adds restrictions on the range of values for integration and simplification. */ + +assume(a>0, e>=0, e<1, i>=0, i<=%pi)$ + +depends([e], [L, G])$ +depends(sin_i, [G, H])$ +depends(cos_i, [G, H])$ + +/* d(e)/dl */ +gradef(e, L, (1 - e^2) / (L * e))$ + +/* d(e)/dg */ +gradef(e, G, -G / (L^2 * e))$ + +/* d(sin_i)/dg */ +gradef(sin_i, G, ( H^2) / (G^3 * sin_i))$ +/* d(sin_i)/dh */ +gradef(sin_i, H, (-H) / (G^2 * sin_i))$ + +/* d(cos_i)/dg */ +gradef(cos_i, G, -H/(G^2))$ +/* d(cos_i)/dh */ +gradef(cos_i, H, 1/G)$ + +/* -------------------------------------------------------------------------- + FUNCTION: hamilton_equations(f, momentavars, coordvars, t) + * -------------------------------------------------------------------------- + * + Description: + * extracts hamiltons canonical equations from a hamiltonian f, + * automatically generating the rates of change of generalized coordinates + * and conjugate momenta. + * + * the function is completely generic and independent of the chosen set of + * canonical variables: delaunay (l,g,h,l,g,h), poincare, whittaker, + * or any other canonical pair (q,p). + * + * arguments: + * f - [symbolic expression] the hamiltonian of the system. + * momentavars - [list] list of conjugate momenta [p1, p2, ..., pn]. + * coordvars - [list] list of generalized coordinates [q1, q2, ..., qn]. + * t - [symbol] independent variable (time). + * + * return: + * [list] list with 2n differential equations in the form: + * [dq1/dt = df/dp1, ..., dqn/dt = df/dpn, + * dp1/dt = -df/dq1, ..., dpn/dt = -df/dqn] + * + * usage examples: + * delaunay: hamilton_equations(h, [l,g,h], [l,g,h], t) + * poincare: hamilton_equations(h_p, [x1,x2,x3], [y1,y2,y3], t) + * extended: hamilton_equations(h, [l,g,h,t], [l,g,h,tau], t) + * + * reference: + * + * + * + * -------------------------------------------------------------------------- */ +hamilton_equations(F, momentaVars, coordVars, t) := block( + [n, eqs_q, eqs_p], + n : length(coordVars), + + /* 1. coordinates: dq/dt = dh/dp */ + eqs_q : makelist('diff(coordVars[i], t) = diff(F, momentaVars[i]), i, 1, n), + + /* 2. momenta: dp/dt = -dh/dq */ + eqs_p : makelist('diff(momentaVars[i], t) = -diff(F, coordVars[i]), i, 1, n), + + return(append(eqs_q, eqs_p)) +)$ + +/* -------------------------------------------------------------------------- + FUNCTION: delaunay_to_poincare(expr) + * -------------------------------------------------------------------------- + * + Description: + * applies the canonical transformation from delaunay (l,g,h,l,g,h) to the + * non-singular poincare variables (y1,y2,y3,x1,x2,x3), removing + * singularities at e=0 and i=0. + * + * the transformation is defined by: + * + * x1 = l, y1 = l + g + h, + * x2 = sqrt(2) * sqrt(l-g) * cos(g+h), y2 = -sqrt(2) * sqrt(l-g) * sin(g+h), + * x3 = sqrt(2) * sqrt(g-h) * cos(h), y3 = -sqrt(2) * sqrt(g-h) * sin(h). + * + * arguments: + * expr - [symbolic expression] expression in delaunay variables + * (l, g, h, l, g, h). + * + * return: + * [expression] the same expression converted to poincare variables + * (x1, x2, x3, y1, y2, y3). + * + * reference: + * - ferraz-mello, s. (2007). canonical perturbation theories: degenerate + * systems and resonance. springer. eq. 7.7. + * - lara, m. (2010). hamiltonian perturbation solutions for spacecraft orbit prediction. eq. 6.25. + * - morbidelli, a. (2002). modern celestial mechanics-aspects of solar system dynamics. eqs. 1.76 and 1.78 + * + * -------------------------------------------------------------------------- */ +delaunay_to_poincare(expr) := block( + [rules, result], + + /* inverse transformation: (x,y) -> (l,g,h,l,g,h) + solving (6.25) for delaunay variables: + + l = x1 + g = x1 - (x2^2 + y2^2)/2 + h = x1 - (x2^2 + y2^2)/2 - (x3^2 + y3^2)/2 + + g + h = atan2(-y2, x2) + h = atan2(-y3, x3) + g = atan2(-y2, x2) - atan2(-y3, x3) + l = y1 - atan2(-y2, x2) + */ + rules : [ + L = X1, + G = X1 - (X2^2 + y2^2)/2, + H = X1 - (X2^2 + y2^2)/2 - (X3^2 + y3^2)/2, + l = y1 - atan2(-y2, X2), + g = atan2(-y2, X2) - atan2(-y3, X3), + h = atan2(-y3, X3) + ], + + result : subst(rules, expr), + return(result) +)$ + + +/* -------------------------------------------------------------------------- + FUNCTION: poincare_to_delaunay(expr) + * -------------------------------------------------------------------------- + * + Description: + * applies the direct transformation, from poincare (y1,y2,y3,x1,x2,x3) to + * delaunay (l,g,h,l,g,h). + * + * arguments: + * expr - [expression] expression in poincare variables. + * + * return: + * [expression] converted to delaunay variables. + * + * reference: + * - ferraz-mello, s. (2007). canonical perturbation theories: degenerate + * systems and resonance. springer. eq. 7.7. + * - lara, m. (2010). hamiltonian perturbation solutions for spacecraft orbit prediction. eq. 6.25. + * - morbidelli, a. (2002). modern celestial mechanics-aspects of solar system dynamics. eqs. 1.76 and 1.78 + * + * -------------------------------------------------------------------------- */ +poincare_to_delaunay(expr) := block( + [rules, result], + + /* direct transformation: (l,g,h,l,g,h) -> (x,y) according to eq. 6.25 */ + rules : [ + X1 = L, + X2 = sqrt(2) * sqrt(L - G) * cos(g + h), + X3 = sqrt(2) * sqrt(G - H) * cos(h), + y1 = l + g + h, + y2 = -sqrt(2) * sqrt(L - G) * sin(g + h), + y3 = -sqrt(2) * sqrt(G - H) * sin(h) + ], + + result : subst(rules, expr), + return(result) +)$ + + + +/* + -------------------------------------------------------------------------- + * implementation note: the use of "_symb" + * -------------------------------------------------------------------------- + * + * it is essential that this function receives all symbols of the + * orbital elements (a_symb, e_symb, etc.) as arguments. + * + * reason: + * lagrange equations depend on partial derivatives like `diff(r, e_symb)`. + * if we used the global symbol 'e', the function would fail silently + * when calculating perturbations that depend on a different symbol (like 'e_3' + * of a third body). + * + * example: + * if r depends on 'e_3' and the function calls `diff(r, e)`, the result + * will be zero, because 'r' does not depend on the global symbol 'e'. + * + * by passing all symbols as arguments, we ensure that the + * coefficients (e.g., 1/n_symb*a_symb) and the derivatives (e.g., diff(r, m_symb)) + * use the correct symbols for any calculation. + -------------------------------------------------------------------------- +*/ + +/* -------------------------------------------------------------------------- + FUNCTION: lagrange_planetary_equations(r, a_symb, e_symb, i_symb, h_symb, g_symb, m_symb, n_symb) + * -------------------------------------------------------------------------- + * + Description: + * calculates lagranges planetary equations, which provide the rate of + * change (time derivative) of the six classical orbital elements + * (a, e, i, omega, omega, m) under the influence of a + * disturbing function r. + * + * this function is modular and expects all symbols used + * in the orbital elements to be passed as arguments, to ensure + * that the partial derivatives (diff) are calculated correctly. + * + * arguments: + * r - [symbolic expression] the disturbing function (potential). + * a_symb - [symbol] the symbol for the semi-major axis. + * e_symb - [symbol] the symbol for the eccentricity. + * i_symb - [symbol] the symbol for the inclination. + * h_symb - [symbol] the symbol for the longitude of the ascending node. + * g_symb - [symbol] the symbol for the argument of perigee. + * m_symb - [symbol] the symbol for the mean anomaly. + * n_symb - [symbol] the symbol for the mean motion. + * + * return: + * [list] a list containing the six symbolic time derivatives: + * [da/dt, de/dt, di/dt, domega/dt, domega/dt, dm/dt] + * + * -------------------------------------------------------------------------- */ +lagrange_planetary_equations(R_pot, a_symb, e_symb, i_symb, h_symb, g_symb, M_symb, n_symb) := block( + + /* declares local variables */ + [dadt, dedt, didt, dOmegadt, domegadt, dMdt], + + dadt: 'diff(a, t) = (2/(n_symb*a_symb)) * diff(R_pot, M_symb), + + dedt: 'diff(e, t) = -(sqrt(1-e_symb^2))/(n_symb*a_symb^2*e_symb)*diff(R_pot, g_symb) + + (1-e_symb^2)/(n_symb*a_symb^2*e_symb)*diff(R_pot, M_symb), + + didt: 'diff(i, t) = -1/(n_symb*a_symb^2*sqrt(1-e_symb^2)*sin(i_symb))*diff(R_pot, h_symb) + + cos(i_symb)/(n_symb*a_symb^2*sqrt(1-e_symb^2)*sin(i_symb))*diff(R_pot, g_symb), + + dOmegadt: 'diff(Omega, t) = 1/(n_symb*a_symb^2*sqrt(1-e_symb^2)*sin(i_symb))*diff(R_pot, i_symb), + + domegadt: 'diff(omega, t) = (sqrt(1-e_symb^2))/(n_symb*a_symb^2*e_symb)*diff(R_pot, e_symb) + - cos(i_symb)/(n_symb*a_symb^2*sqrt(1-e_symb^2)*sin(i_symb))*diff(R_pot, i_symb), + + dMdt: 'diff(M, t) = n_symb - 2/(n_symb*a_symb)*diff(R_pot, a_symb) + - (1-e_symb^2)/(n_symb*a_symb^2*e_symb)*diff(R_pot, e_symb), + + return([dadt, dedt, didt, dOmegadt, domegadt, dMdt]) +)$ + + +/* -------------------------------------------------------------------------- + FUNCTION: lagrange_planetary_equations_nonsingular + * -------------------------------------------------------------------------- + * + Description: + * calculates lagranges planetary equations in non-singular variables + * (modified equinoctial elements), which remove the singularities + * present in classical equations for circular (e = 0) and + * equatorial (i = 0) orbits. + * + * the set of variables used is: + * a : semi-major axis + * h1 : e * sin(omega + omega) + * k1 : e * cos(omega + omega) + * p1 : tan(i/2) * sin(omega) + * q1 : tan(i/2) * cos(omega) + * lambda : m + omega + omega (mean longitude) + * + * arguments: + * r_pot - [symbolic expression] the disturbing function. + * must be expressed in terms of equinoctial elements + * (a_s, h1_s, k1_s, p1_s, q1_s, lambda_s). + * a_s - [symbol] semi-major axis. + * h1_s - [symbol] e*sin(varpi), where varpi = omega + omega. + * k1_s - [symbol] e*cos(varpi). + * p1_s - [symbol] tan(i/2)*sin(omega). + * q1_s - [symbol] tan(i/2)*cos(omega). + * lambda_s - [symbol] mean longitude (m + omega + omega). + * n_s - [symbol] mean motion. + * + * return: + * [list] [da/dt, dh1/dt, dk1/dt, dp1/dt, dq1/dt, dlambda/dt] + * + * references: + * - broucke, r. a. and cefola, p. j. (1972). on the equinoctial orbit + * elements. celestial mechanics, 5(3), 303-310. + * - walker, m. j. h., ireland, b. and owens, j. (1985). a set of modified + * equinoctial orbit elements. celestial mechanics, 36, 409-419. + * - battin, r. h. (1999). an introduction to the methods of astrodynamics. + * aiaa education series. sec. 10.3. + * + * -------------------------------------------------------------------------- */ +lagrange_planetary_equations_nonsingular(R_pot, a_s, h1_s, k1_s, p1_s, q1_s, lambda_s, n_s) := block( + [dadt, dh1dt, dk1dt, dp1dt, dq1dt, dlambdadt, + e2_s, eta_s, A_s, B_s, C_s], + + /* abbreviations */ + e2_s : h1_s^2 + k1_s^2, /* e^2 */ + eta_s : sqrt(1 - e2_s), /* sqrt(1 - e^2) */ + A_s : n_s * a_s^2, /* n * a^2 */ + B_s : 1 + p1_s^2 + q1_s^2, /* 1 + tan(i/2)^2 = 2/(1+cos(i)) = 2*sec^2(i/2) ... no: = 1 + p^2 + q^2 */ + + /* 1. semi-major axis */ + dadt : 'diff(a, t) = (2 / (n_s * a_s)) * diff(R_pot, lambda_s), + + /* 2. h1 = e*sin(varpi) */ + dh1dt : 'diff(h1, t) = -(eta_s / A_s) * diff(R_pot, k1_s) + + (eta_s * k1_s / (A_s * (1 + eta_s))) * diff(R_pot, lambda_s) + - (h1_s / (2 * A_s * eta_s)) * (p1_s * diff(R_pot, p1_s) + q1_s * diff(R_pot, q1_s)), + + /* 3. k1 = e*cos(varpi) */ + dk1dt : 'diff(k1, t) = (eta_s / A_s) * diff(R_pot, h1_s) + + (eta_s * h1_s / (A_s * (1 + eta_s))) * diff(R_pot, lambda_s) + + (k1_s / (2 * A_s * eta_s)) * (p1_s * diff(R_pot, p1_s) + q1_s * diff(R_pot, q1_s)), + + /* 4. p1 = tan(i/2)*sin(omega) */ + dp1dt : 'diff(p1, t) = (B_s / (2 * A_s * eta_s)) * diff(R_pot, q1_s), + + /* 5. q1 = tan(i/2)*cos(omega) */ + dq1dt : 'diff(q1, t) = -(B_s / (2 * A_s * eta_s)) * diff(R_pot, p1_s), + + /* 6. mean longitude (lambda = m + omega + omega) */ + dlambdadt : 'diff(lambda, t) = n_s + - (2 / (n_s * a_s)) * diff(R_pot, a_s) + - (eta_s / A_s * (1 + eta_s)) + * (h1_s * diff(R_pot, h1_s) + k1_s * diff(R_pot, k1_s)) + + (1 / (2 * A_s * eta_s)) + * (p1_s * diff(R_pot, p1_s) + q1_s * diff(R_pot, q1_s)), + + return([dadt, dh1dt, dk1dt, dp1dt, dq1dt, dlambdadt]) +)$ + + +/* -------------------------------------------------------------------------- + FUNCTION: gauss_planetary_equations + * -------------------------------------------------------------------------- + * + Description: + * calculates gauss planetary equations in rsw form + * (radial-transverse-normal), which provide the rate of change of the six + * classical orbital elements (a, e, i, omega, omega, m) under the influence + * of an arbitrary perturbing force decomposed into the three orbital components. + * + * unlike lagrange equations, which require a scalar disturbing + * function r (potential), gauss equations accept arbitrary perturbing + * forces, including non-conservative forces (e.g., drag, propulsion). + * + * arguments: + * a_dr - [expression] radial component of the perturbing acceleration + * (positive outwards, in the direction of the position vector). + * a_dt - [expression] transverse component of the perturbing acceleration + * (positive in the direction of motion, perpendicular to a_dr + * in the orbital plane). + * a_dh - [expression] normal component of the perturbing acceleration + * (positive in the direction of the orbital angular momentum). + * a_symb - [symbol] semi-major axis. + * e_symb - [symbol] eccentricity. + * i_symb - [symbol] inclination. + * h_symb - [symbol] longitude of the ascending node (omega). + * g_symb - [symbol] argument of perigee (omega). + * f_symb - [symbol] true anomaly. + * n_symb - [symbol] mean motion. + * + * return: + * [list] [da/dt, de/dt, di/dt, domega/dt, domega/dt, dm/dt] + * + * reference: + * battin, r. h. (1999). an introduction to the mathematics and methods of astrodynamics. + * aiaa education series. eq. 10.41. + * + * -------------------------------------------------------------------------- */ +gauss_planetary_equations(a_dr, a_dt, a_dh, a_symb, e_symb, i_symb, + h_symb, g_symb, f_symb, n_symb) := block( + [dadt, dedt, didt, dOmegadt, domegadt, dMdt, + p_s, r_s, hm_s, b_s, theta_s], + + /* orbital abbreviations */ + p_s : a_symb * (1 - e_symb^2), + r_s : p_s / (1 + e_symb * cos(f_symb)), + hm_s : n_symb * a_symb^2 * sqrt(1 - e_symb^2), + b_s : a_symb * sqrt(1 - e_symb^2), + theta_s : f_symb + g_symb, + + /* 1. semi-major axis (da/dt) */ + dadt : 'diff(a, t) = (2 * a_symb^2 / hm_s) + * (e_symb * sin(f_symb) * a_dr + (p_s / r_s) * a_dt), + + /* 2. eccentricity (de/dt) */ + dedt : 'diff(e, t) = (1 / hm_s) + * (p_s * sin(f_symb) * a_dr + + ((p_s + r_s) * cos(f_symb) + r_s * e_symb) * a_dt), + + /* 3. inclination (di/dt) */ + didt : 'diff(i, t) = (r_s * cos(theta_s) / hm_s) * a_dh, + + /* 4. longitude of the ascending node (domega/dt) */ + dOmegadt : 'diff(Omega, t) = (r_s * sin(theta_s)) + / (hm_s * sin(i_symb)) * a_dh, + + /* 5. argument of perigee (domega/dt) */ + domegadt : 'diff(omega, t) = (1 / (hm_s * e_symb)) + * (-p_s * cos(f_symb) * a_dr + (p_s + r_s) * sin(f_symb) * a_dt) + - (r_s * sin(theta_s) * cos(i_symb)) + / (hm_s * sin(i_symb)) * a_dh, + + /* 6. mean anomaly (dm/dt) */ + dMdt : 'diff(M, t) = n_symb + + (b_s / (a_symb * hm_s * e_symb)) + * ((p_s * cos(f_symb) - 2 * r_s * e_symb) * a_dr + - (p_s + r_s) * sin(f_symb) * a_dt), + + return([dadt, dedt, didt, dOmegadt, domegadt, dMdt]) +)$ + +/* ============================================================================ + FUNCTION: prepare_planetary_equations + Description: + prepares the perturbing potential list for planetary equations. + + Parameters: + r_pot_list: list of perturbing potential expressions. + order_e: truncation order for eccentricity. + r_exp: conic expression for the radial distance. + a_symb: symbol for semi-major axis. + e_symb: symbol for eccentricity. + i_symb: symbol for inclination. + f_symb: symbol for true anomaly. + m_symb: symbol for mean anomaly. + e_symb_ecc: symbol for eccentric anomaly. + r_symb: symbol for radial distance. + + Returns: + the final prepared potential series substituting short period terms. + ============================================================================ */ +prepare_planetary_equations(R_pot_list, order_e, r_exp, a_symb, e_symb, i_symb, f_symb, M_symb, E_symb, r_symb) := block( + [R_sum, R_sum_trig, R_zonal_no_r, R_zonal_exp, cos_f, sin_f, R_zonal_series], + + /* 1. sums the potentials in the list */ + R_sum: apply("+", R_pot_list), + + /* 2. replaces the atoms sin_i/cos_i with the correct functional symbol */ + R_sum_trig: subst([sin_i=sin(i_symb), cos_i=cos(i_symb)], R_sum), + + /* 3. replaces the radius r with the conic expression (r_exp) */ + R_zonal_no_r: subst(r_symb=r_exp, R_sum_trig), + + /* 4. expands trigonometric terms mixing f_symb with other angles (e.g., g) */ + R_zonal_exp: celmec_trigexpand(R_zonal_no_r, f_symb), + + /* 5. generates the kepler series for the short period */ + cos_f: get_cos_f_series(order_e, M_symb, e_symb, E_symb), + sin_f: get_sin_f_series(order_e, M_symb, e_symb, E_symb), + + /* 6. replaces the series in the expanded potential */ + R_zonal_series: subst([cos(f_symb)=cos_f, sin(f_symb)=sin_f], R_zonal_exp), + + return(R_zonal_series) +)$ + +shows_available_functions() :=block( + print("available functions:"), + print(" - lagrange_planetary_equations(R, a_symb, e_symb, i_symb, h_symb, g_symb, M_symb, n_symb)") +)$ + +/* turns confirmation back to normal */ +suppress_trace_warnings: false$ diff --git a/symcelmech/rtest_symcelmech.mac b/symcelmech/rtest_symcelmech.mac new file mode 100644 index 0000000..c74a424 --- /dev/null +++ b/symcelmech/rtest_symcelmech.mac @@ -0,0 +1,207 @@ +/* rtest_symcelmech.mac */ +/* regression test suite for the symcelmech package */ + +(load("symcelmech.mac"), 'done); +'done$ + +/* ============================================================================ + test group 1: elliptic expansions (celmechseries) + ============================================================================ */ + +/* cos(f) series at order 2 should recover the classical expansion */ +subst(e = 0, ratdisrep(get_cos_f_series(1, M, e, E))); +cos(M)$ + +/* sin(f) series at order 0 is simply sin(M) */ +ratdisrep(get_sin_f_series(0, M, e, E)); +sin(M)$ + +/* r/a series at order 0 is 1 */ +subst(e = 0, ratdisrep(get_r_over_a_kepler_series(1, M, e, E))); +1$ + +/* ============================================================================ + test group 2: analytical geometry (analyticalfunctions) + ============================================================================ */ + +/* exact radial distance at f=0 gives periapsis */ +r_two_body_exact(e, a, 0); +a*(1 - e^2)/(1 + e)$ + +/* exact radial distance at f=pi gives apoapsis */ +r_two_body_exact(e, a, %pi); +a*(1 - e^2)/(1 - e)$ + +/* jacobian dM/df at e=0 is 1 (circular orbit) */ +subst(e = 0, dM_f(e, f)); +1$ + +/* ============================================================================ + test group 3: potential derivation (derivepotential) + ============================================================================ */ + +/* J2 zonal potential should have the correct structure */ +(U2: derive_zonal_potential(2, mu, r, f, a, e, i, g, R, false), + freeof(j3, U2)); +true$ + +/* J2 potential should contain j2 */ +(U2: derive_zonal_potential(2, mu, r, f, a, e, i, g, R, false), + not freeof(j2, U2)); +true$ + +/* J3 zonal potential should contain j3 */ +(U3: derive_zonal_potential(3, mu, r, f, a, e, i, g, R, false), + not freeof(j3, U3)); +true$ + +/* compact mode should introduce k_j2 */ +(U2c: derive_zonal_potential(2, mu, r, f, a, e, i, g, R, true), + not freeof(k_j2, U2c)); +true$ + +/* inertial vector should be a unit vector (norm = 1) at e=0, circular orbit */ +(v: get_inertial_vector(f, 0, 0, 0), + trigsimp(v[1]^2 + v[2]^2 + v[3]^2)); +1$ + +/* ============================================================================ + test group 4: averaging and hansen coefficients (average) + ============================================================================ */ + +/* hansen X0(1,0) = 1 + e^2/2 */ +hansen_x0_pos_calc[1, 0, e]; +1 + e^2/2$ + +/* hansen X0(2,0) = 1 + 3*e^2/2 */ +hansen_x0_pos_calc[2, 0, e]; +1 + 3*e^2/2$ + +/* hansen X0(1,1) = -3*e/2 */ +hansen_x0_pos_calc[1, 1, e]; +-3*e/2$ + +/* map_average: constant expression is entirely secular */ +map_average(42, f); +[42, 0]$ + +/* map_average: cos(f) is entirely periodic */ +map_average(cos(f), f); +[0, cos(f)]$ + +/* map_average: sum with mixed terms */ +map_average(A + B*cos(f), f); +[A, B*cos(f)]$ + +/* ============================================================================ + test group 5: algebraic utilities (algebraicutils) + ============================================================================ */ + +/* chop_series removes small coefficients */ +chop_series(1e-20*x + 3*y, 1e-15); +3*y$ + +/* chop_series preserves large coefficients */ +chop_series(5*x + 3*y, 1e-15); +5*x + 3*y$ + +/* coeff_of_term extracts the number from a product */ +coeff_of_term(7*x*y); +7$ + +/* coeff_of_term of a pure number */ +coeff_of_term(42); +42$ + +/* deatomize_parameters returns a list */ +listp(deatomize_parameters()); +true$ + +/* ============================================================================ + test group 6: poisson brackets (analyticalfunctions) + ============================================================================ */ + +/* pb_series: bracket of a function with itself is zero */ +pb_series(L^2, L^2, [l, g, h, L, G, H]); +0$ + +/* pb_series: {l, L} = 1 (canonical pair) */ +pb_series(l, L, [l, g, h, L, G, H]); +1$ + +/* pb_series: {g, G} = 1 */ +pb_series(g, G, [l, g, h, L, G, H]); +1$ + +/* pb_series: {h, H} = 1 */ +pb_series(h, H, [l, g, h, L, G, H]); +1$ + +/* pb_series: antisymmetry {L, l} = -1 */ +pb_series(L, l, [l, g, h, L, G, H]); +-1$ + +/* pb_series: non-conjugate pair is zero */ +pb_series(l, G, [l, g, h, L, G, H]); +0$ + +/* ============================================================================ + test group 7: planetary equations (planetaryequations) + ============================================================================ */ + +/* hamilton_equations returns 2n equations */ +(eqs: hamilton_equations(H_test, [L, G, H], [l, g, h], t), + length(eqs)); +6$ + +/* lagrange_planetary_equations returns 6 equations */ +(eqs: lagrange_planetary_equations(R_pot, a, e, i, h, g, M, n), + length(eqs)); +6$ + +/* gauss_planetary_equations returns 6 equations */ +(eqs: gauss_planetary_equations(a_dr, a_dt, a_dh, a, e, i, h, g, f, n), + length(eqs)); +6$ + +/* ============================================================================ + test group 8: astrodata + ============================================================================ */ + +/* earth parameters should be a non-empty list */ +(params: get_astro_params("Earth"), + listp(params) and length(params) > 0); +true$ + +/* moon parameters should contain j2 */ +(params: get_astro_params("Moon"), + not freeof(j2, params)); +true$ + +/* case insensitivity */ +(params: get_astro_params("earth"), + listp(params) and length(params) > 0); +true$ + +/* unknown body returns empty list */ +get_astro_params("Pluto"); +[]$ + +/* ============================================================================ + test group 9: generating function quadrature + ============================================================================ */ + +/* integral of cos(f) should be sin(f) */ +generating_function_quadrature(cos(f), f); +sin(f)$ + +/* integral of sin(f) should be -cos(f) */ +generating_function_quadrature(sin(f), f); +-cos(f)$ + +/* integral of cos(2*f) should be sin(2*f)/2 */ +generating_function_quadrature(cos(2*f), f); +sin(2*f)/2$ + +'done; +'done$ diff --git a/symcelmech/symcelmech.mac b/symcelmech/symcelmech.mac new file mode 100644 index 0000000..f2bf426 --- /dev/null +++ b/symcelmech/symcelmech.mac @@ -0,0 +1,19 @@ +/* symcelmech.mac -- entry point for the symcelmech package */ + +/* ensure the package directory is in the search path */ +if not member(pathname_directory(load_pathname), file_search_maxima) then + file_search_maxima: cons(pathname_directory(load_pathname), file_search_maxima)$ + +/* core modules */ +load("astrodata.mac")$ +load("celmechseries.mac")$ +load("analyticalfunctions.mac")$ +load("algebraicutils.mac")$ +load("average.mac")$ +load("derivepotential.mac")$ +load("lietransformations.mac")$ +load("planetaryequations.mac")$ + +symcelmech_version : "1.0.0"$ + +print(concat("symcelmech ", symcelmech_version, " loaded."))$