gotelem/py/notebooks/test_solar.ipynb

325 lines
8.8 KiB
Plaintext
Raw Normal View History

2023-06-22 23:29:08 +00:00
{
"cells": [
{
"cell_type": "code",
2023-07-10 23:10:21 +00:00
"execution_count": 1,
2023-06-22 23:29:08 +00:00
"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": [
2023-07-10 23:10:21 +00:00
"/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"
2023-06-22 23:29:08 +00:00
]
}
],
"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)"
]
},
2023-07-10 23:10:21 +00:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
2023-06-22 23:29:08 +00:00
{
"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": {
2023-07-10 23:10:21 +00:00
"display_name": "telemstrategy",
2023-06-22 23:29:08 +00:00
"language": "python",
2023-07-10 23:10:21 +00:00
"name": "telemstrategy"
2023-06-22 23:29:08 +00:00
},
"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
}