saji/gotelem
1
0
Fork 0
Code Issues 7 Actions Packages Releases 2 Activity
gotelem/py/notebooks/test_solar.ipynb
saji c772d6e95f some python crap
2023-07-10 18:10:21 -05:00

325 lines
8.8 KiB
Plaintext
Raw Permalink Blame History

{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2023-06-21T00:28:49.748311944Z",
"start_time": "2023-06-21T00:28:49.744946948Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/saji/Documents/Code/buildroot/gotelem/py/pytelem/optimus.py:50: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
" @jit\n",
"/home/saji/Documents/Code/buildroot/gotelem/py/pytelem/optimus.py:65: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
" @jit\n",
"/home/saji/Documents/Code/buildroot/gotelem/py/pytelem/optimus.py:79: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
" @jit\n",
"/home/saji/Documents/Code/buildroot/gotelem/py/pytelem/optimus.py:479: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
" @jit\n"
]
}
],
"source": [
"\n",
"from pytelem.optimus import *\n",
"from jax import jit"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2023-06-21T00:28:50.052452543Z",
"start_time": "2023-06-21T00:28:50.050159245Z"
},
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'optim' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[2], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m ffast \u001b[39m=\u001b[39m jit(optim\u001b[39m.\u001b[39msolar_position)\n",
"\u001b[0;31mNameError\u001b[0m: name 'optim' is not defined"
]
}
],
"source": [
"ffast = jit(optim.solar_position)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"ExecuteTime": {
"end_time": "2023-06-21T00:29:01.075229307Z",
"start_time": "2023-06-21T00:29:01.042571064Z"
},
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": [
"timestamp = np.array(1687306901)\n",
"timestamps = np.array([1687306901, 1687306902, 1687306903, 1687306906])"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false,
"jupyter": {
"outputs_hidden": false
}
},
"outputs": [],
"source": [
"\n",
"jd = timestamps / 86400.0 + 2440587.5\n",
"jde = jd + DELTA_T / 86400.0\n",
"jce = (jde - 2451545) / 36525\n",
"\n",
"jme = jce / 10\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"\n",
"# todo: make more elegant?\n",
"# TODO: vectorize later? it's kinda complex\n",
"\n",
"# heliocentric longitude\n",
"l_rad = np.zeros_like(timestamp)\n",
"for idx, vec in enumerate(HELIO_L):\n",
" l_rad = l_rad + helio_vector(vec, jme) * jme ** idx\n",
"\n",
"l_rad = l_rad / 10e8\n",
"l_deg = np.rad2deg(l_rad) % 360\n",
"\n",
"# heliocentric latitude\n",
"b_rad = np.zeros_like(timestamp)\n",
"for idx, vec in enumerate(HELIO_B):\n",
" b_rad = b_rad + helio_vector(vec, jme) * jme ** idx\n",
"b_rad = b_rad / 10e8\n",
"b_deg = b_rad % 360\n",
"\n",
"# heliocentric radius\n",
"r_rad = np.zeros_like(timestamp)\n",
"for idx, vec in enumerate(HELIO_R):\n",
" r_rad = r_rad + helio_vector(vec, jme) * jme ** idx\n",
"r_rad = r_rad / 10e8\n",
"r_deg = r_rad % 360\n",
"\n",
"theta = (l_deg + 180) % 360\n",
"beta = -1 * b_deg\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(4, 5)\n"
]
}
],
"source": [
"\n",
"def cubic_poly(a,b,c,d):\n",
" return a + b * jce + c * jce ** 2 + (jce ** 3) / d\n",
"X0 = cubic_poly(297.85036, 445267.111480, -0.0019142, 189474)\n",
"X1 = cubic_poly(357.52772, 35999.050340, -0.0001603, -300000)\n",
"X2 = cubic_poly(134.96298, 477198.867398, 0.0086972, 56250)\n",
"X3 = cubic_poly(93.27191, 483202.017538, -0.0036825, 327270)\n",
"X4 = cubic_poly(125.04452, 1934.136261, 0.0020708, 450000)\n",
"\n",
"X = np.vstack([X0, X1, X2, X3, X4]).T\n",
"print(X.shape)\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"\n",
"nut = NUTATION_ABCD_ARRAY\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(4, 63)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(nut[:,0] + jce[..., np.newaxis] * nut[:,1])"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(4, 5)\n",
"(63, 5)\n",
"d psi shape (4, 63)\n",
"[-0.00333032 -0.00333032 -0.00333032 -0.00333032]\n"
]
}
],
"source": [
"\n",
"print(X.shape)\n",
"print(NUTATION_YTERM_ARRAY.shape)\n",
"d_psi = (nut[:,0] + jce[..., np.newaxis] * nut[:,1]) * np.sin(np.sum(X[:, np.newaxis, :] * NUTATION_YTERM_ARRAY[np.newaxis, ...], axis=-1))\n",
"\n",
"print(f\"d psi shape {d_psi.shape}\")\n",
"\n",
"nut_long = np.sum(d_psi, axis=-1) / 36,000,000\n",
"\n",
"print(nut_long)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(4, 11)\n"
]
}
],
"source": [
"u = timestamps[:, np.newaxis] ** np.arange(0,11)\n",
"print(u.shape)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"\n",
"epsilon_0 = np.array(\n",
" [[\n",
" 84381.448,\n",
" -4680.93,\n",
" 1.55,\n",
" 1999.25,\n",
" -51.38,\n",
" -249.67,\n",
" -39.05,\n",
" 7.12,\n",
" 27.87,\n",
" 5.79,\n",
" 2.45,\n",
" ]]\n",
" )\n",
"res = u * epsilon_0\n"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Array([-8.6115017e+12, -5.7323689e+12, -6.4693841e+12, -1.1346544e+13], dtype=float32)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res.sum(axis=-1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "telemstrategy",
"language": "python",
"name": "telemstrategy"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Reference in a new issue View git blame Copy permalink