API tools faq
paste
Guest User

Earth orbit simulation.ipynb

a guest
Aug 1st, 2020
139
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 36.71 KB | None | 0 0
raw download clone embed print report
  1. {
  2. "cells": [
  3. {
  4. "cell_type": "markdown",
  5. "metadata": {},
  6. "source": [
  7. "## Simulating the Earth's Orbit"
  8. ]
  9. },
  10. {
  11. "cell_type": "markdown",
  12. "metadata": {},
  13. "source": [
  14. "This notebook collects the calculations done in order to answer the following question posted on math.stackexchange\n",
  15. "\n",
  16. "https://math.stackexchange.com/questions/3766767/is-there-a-simple-ish-function-for-modeling-seasonal-changes-to-day-night-durati\n",
  17. "\n",
  18. "The answer is [this one](https://math.stackexchange.com/questions/3766767/is-there-a-simple-ish-function-for-modeling-seasonal-changes-to-day-night-durati/3772816#3772816), by user \"JonathanZ supports MonicaC\", who is also the author of this notebook."
  19. ]
  20. },
  21. {
  22. "cell_type": "markdown",
  23. "metadata": {},
  24. "source": [
  25. "This notebook requires SymPy, the Python symbolic manipulation library. \n",
  26. "\n",
  27. "Also, a lot of the code in here was written to help produce LaTeX for the linked answer, and to create the Desmos graphs that are linked from the posted answer. You probably will need to read that answer to understand what is happening in this notebook."
  28. ]
  29. },
  30. {
  31. "cell_type": "code",
  32. "execution_count": 1,
  33. "metadata": {},
  34. "outputs": [],
  35. "source": [
  36. "%reset -f"
  37. ]
  38. },
  39. {
  40. "cell_type": "code",
  41. "execution_count": 2,
  42. "metadata": {},
  43. "outputs": [],
  44. "source": [
  45. "import sympy as sp"
  46. ]
  47. },
  48. {
  49. "cell_type": "markdown",
  50. "metadata": {},
  51. "source": [
  52. "## Create parameters"
  53. ]
  54. },
  55. {
  56. "cell_type": "code",
  57. "execution_count": 3,
  58. "metadata": {},
  59. "outputs": [],
  60. "source": [
  61. "# lattitude\n",
  62. "phi = sp.symbols('phi')"
  63. ]
  64. },
  65. {
  66. "cell_type": "code",
  67. "execution_count": 4,
  68. "metadata": {},
  69. "outputs": [],
  70. "source": [
  71. "# axial tilt of the Earth\n",
  72. "eps = sp.symbols('epsilon')"
  73. ]
  74. },
  75. {
  76. "cell_type": "code",
  77. "execution_count": 5,
  78. "metadata": {},
  79. "outputs": [],
  80. "source": [
  81. "# alpha is time of day in radians\n",
  82. "# i.e. alpha: 0 -> 2*pi = midnight to midnight\n",
  83. "#\n",
  84. "# beta is day of the year in radians\n",
  85. "# i.e. beta: 0 -> 2*pi = winter solstice to winter solstice\n",
  86. "alpha, beta = sp.symbols('alpha beta')"
  87. ]
  88. },
  89. {
  90. "cell_type": "markdown",
  91. "metadata": {},
  92. "source": [
  93. "## Direction vectors and Rotation matrices"
  94. ]
  95. },
  96. {
  97. "cell_type": "code",
  98. "execution_count": 6,
  99. "metadata": {},
  100. "outputs": [
  101. {
  102. "data": {
  103. "text/latex": [
  104. "$\\displaystyle \\left[\\begin{matrix}\\cos{\\left(\\phi \\right)}\\\\0\\\\\\sin{\\left(\\phi \\right)}\\end{matrix}\\right]$"
  105. ],
  106. "text/plain": [
  107. "Matrix([\n",
  108. "[cos(phi)],\n",
  109. "[ 0],\n",
  110. "[sin(phi)]])"
  111. ]
  112. },
  113. "execution_count": 6,
  114. "metadata": {},
  115. "output_type": "execute_result"
  116. }
  117. ],
  118. "source": [
  119. "# Surface normal to point at lattitude phi\n",
  120. "N = sp.Matrix([sp.cos(phi), 0, sp.sin(phi)])\n",
  121. "N\n",
  122. "\n"
  123. ]
  124. },
  125. {
  126. "cell_type": "code",
  127. "execution_count": 7,
  128. "metadata": {},
  129. "outputs": [
  130. {
  131. "name": "stdout",
  132. "output_type": "stream",
  133. "text": [
  134. "\\left[\\begin{matrix}\\cos{\\left(\\phi \\right)}\\\\0\\\\\\sin{\\left(\\phi \\right)}\\end{matrix}\\right]\n"
  135. ]
  136. }
  137. ],
  138. "source": [
  139. "print(sp.latex(N))"
  140. ]
  141. },
  142. {
  143. "cell_type": "code",
  144. "execution_count": 8,
  145. "metadata": {},
  146. "outputs": [
  147. {
  148. "data": {
  149. "text/latex": [
  150. "$\\displaystyle \\left[\\begin{matrix}\\cos{\\left(\\alpha \\right)} & \\sin{\\left(\\alpha \\right)} & 0\\\\- \\sin{\\left(\\alpha \\right)} & \\cos{\\left(\\alpha \\right)} & 0\\\\0 & 0 & 1\\end{matrix}\\right]$"
  151. ],
  152. "text/plain": [
  153. "Matrix([\n",
  154. "[ cos(alpha), sin(alpha), 0],\n",
  155. "[-sin(alpha), cos(alpha), 0],\n",
  156. "[ 0, 0, 1]])"
  157. ]
  158. },
  159. "execution_count": 8,
  160. "metadata": {},
  161. "output_type": "execute_result"
  162. }
  163. ],
  164. "source": [
  165. "# Rotation of the Earth throughout the day\n",
  166. "\n",
  167. "M_rot = sp.Matrix([\n",
  168. " [sp.cos(alpha), sp.sin(alpha), 0],\n",
  169. " [-sp.sin(alpha), sp.cos(alpha), 0],\n",
  170. " [0, 0, 1]\n",
  171. "])\n",
  172. "M_rot"
  173. ]
  174. },
  175. {
  176. "cell_type": "code",
  177. "execution_count": 9,
  178. "metadata": {},
  179. "outputs": [
  180. {
  181. "name": "stdout",
  182. "output_type": "stream",
  183. "text": [
  184. "\\left[\\begin{matrix}\\cos{\\left(\\alpha \\right)} & \\sin{\\left(\\alpha \\right)} & 0\\\\- \\sin{\\left(\\alpha \\right)} & \\cos{\\left(\\alpha \\right)} & 0\\\\0 & 0 & 1\\end{matrix}\\right]\n"
  185. ]
  186. }
  187. ],
  188. "source": [
  189. "print(sp.latex(M_rot))"
  190. ]
  191. },
  192. {
  193. "cell_type": "code",
  194. "execution_count": 10,
  195. "metadata": {},
  196. "outputs": [
  197. {
  198. "data": {
  199. "text/latex": [
  200. "$\\displaystyle \\left[\\begin{matrix}\\cos{\\left(\\epsilon \\right)} & 0 & \\sin{\\left(\\epsilon \\right)}\\\\0 & 1 & 0\\\\- \\sin{\\left(\\epsilon \\right)} & 0 & \\cos{\\left(\\epsilon \\right)}\\end{matrix}\\right]$"
  201. ],
  202. "text/plain": [
  203. "Matrix([\n",
  204. "[ cos(epsilon), 0, sin(epsilon)],\n",
  205. "[ 0, 1, 0],\n",
  206. "[-sin(epsilon), 0, cos(epsilon)]])"
  207. ]
  208. },
  209. "execution_count": 10,
  210. "metadata": {},
  211. "output_type": "execute_result"
  212. }
  213. ],
  214. "source": [
  215. "# Axial tilt of the Earth\n",
  216. "\n",
  217. "M_tilt = sp.Matrix([\n",
  218. " [sp.cos(eps),0,sp.sin(eps)],\n",
  219. " [0,1,0],\n",
  220. " [-sp.sin(eps),0,sp.cos(eps)]\n",
  221. "])\n",
  222. "M_tilt"
  223. ]
  224. },
  225. {
  226. "cell_type": "code",
  227. "execution_count": 11,
  228. "metadata": {},
  229. "outputs": [
  230. {
  231. "name": "stdout",
  232. "output_type": "stream",
  233. "text": [
  234. "\\left[\\begin{matrix}\\cos{\\left(\\epsilon \\right)} & 0 & \\sin{\\left(\\epsilon \\right)}\\\\0 & 1 & 0\\\\- \\sin{\\left(\\epsilon \\right)} & 0 & \\cos{\\left(\\epsilon \\right)}\\end{matrix}\\right]\n"
  235. ]
  236. }
  237. ],
  238. "source": [
  239. "print(sp.latex(M_tilt))"
  240. ]
  241. },
  242. {
  243. "cell_type": "code",
  244. "execution_count": 12,
  245. "metadata": {},
  246. "outputs": [
  247. {
  248. "data": {
  249. "text/latex": [
  250. "$\\displaystyle \\left[\\begin{matrix}\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} + \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}\\\\- \\sin{\\left(\\alpha \\right)} \\cos{\\left(\\phi \\right)}\\\\- \\sin{\\left(\\epsilon \\right)} \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\phi \\right)} + \\sin{\\left(\\phi \\right)} \\cos{\\left(\\epsilon \\right)}\\end{matrix}\\right]$"
  251. ],
  252. "text/plain": [
  253. "Matrix([\n",
  254. "[ sin(epsilon)*sin(phi) + cos(alpha)*cos(epsilon)*cos(phi)],\n",
  255. "[ -sin(alpha)*cos(phi)],\n",
  256. "[-sin(epsilon)*cos(alpha)*cos(phi) + sin(phi)*cos(epsilon)]])"
  257. ]
  258. },
  259. "execution_count": 12,
  260. "metadata": {},
  261. "output_type": "execute_result"
  262. }
  263. ],
  264. "source": [
  265. "N_ab = M_tilt * M_rot * N\n",
  266. "N_ab"
  267. ]
  268. },
  269. {
  270. "cell_type": "code",
  271. "execution_count": 13,
  272. "metadata": {},
  273. "outputs": [
  274. {
  275. "name": "stdout",
  276. "output_type": "stream",
  277. "text": [
  278. "\\left[\\begin{matrix}\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} + \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}\\\\- \\sin{\\left(\\alpha \\right)} \\cos{\\left(\\phi \\right)}\\\\- \\sin{\\left(\\epsilon \\right)} \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\phi \\right)} + \\sin{\\left(\\phi \\right)} \\cos{\\left(\\epsilon \\right)}\\end{matrix}\\right]\n"
  279. ]
  280. }
  281. ],
  282. "source": [
  283. "print(sp.latex(N_ab))"
  284. ]
  285. },
  286. {
  287. "cell_type": "code",
  288. "execution_count": 14,
  289. "metadata": {},
  290. "outputs": [
  291. {
  292. "data": {
  293. "text/latex": [
  294. "$\\displaystyle 1$"
  295. ],
  296. "text/plain": [
  297. "1"
  298. ]
  299. },
  300. "execution_count": 14,
  301. "metadata": {},
  302. "output_type": "execute_result"
  303. }
  304. ],
  305. "source": [
  306. "# Verifying unit length\n",
  307. "sp.simplify(N_ab.dot(N_ab))"
  308. ]
  309. },
  310. {
  311. "cell_type": "code",
  312. "execution_count": 15,
  313. "metadata": {},
  314. "outputs": [
  315. {
  316. "name": "stdout",
  317. "output_type": "stream",
  318. "text": [
  319. "\\left[\\begin{matrix}\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} + \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}\\\\- \\sin{\\left(\\alpha \\right)} \\cos{\\left(\\phi \\right)}\\\\- \\sin{\\left(\\epsilon \\right)} \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\phi \\right)} + \\sin{\\left(\\phi \\right)} \\cos{\\left(\\epsilon \\right)}\\end{matrix}\\right]\n"
  320. ]
  321. }
  322. ],
  323. "source": [
  324. "print(sp.latex(N_ab))"
  325. ]
  326. },
  327. {
  328. "cell_type": "code",
  329. "execution_count": 16,
  330. "metadata": {},
  331. "outputs": [
  332. {
  333. "data": {
  334. "text/latex": [
  335. "$\\displaystyle \\left[\\begin{matrix}- \\cos{\\left(\\beta \\right)}\\\\- \\sin{\\left(\\beta \\right)}\\\\0\\end{matrix}\\right]$"
  336. ],
  337. "text/plain": [
  338. "Matrix([\n",
  339. "[-cos(beta)],\n",
  340. "[-sin(beta)],\n",
  341. "[ 0]])"
  342. ]
  343. },
  344. "execution_count": 16,
  345. "metadata": {},
  346. "output_type": "execute_result"
  347. }
  348. ],
  349. "source": [
  350. "# Instead of putting the Earth in the appropriate place in its orbit\n",
  351. "# all we care about is where it is relative to the Sun, so we'll just \n",
  352. "# 'move' the Sun instead.\n",
  353. "\n",
  354. "sol_b = sp.Matrix([-sp.cos(beta), -sp.sin(beta), 0])\n",
  355. "sol_b"
  356. ]
  357. },
  358. {
  359. "cell_type": "code",
  360. "execution_count": 17,
  361. "metadata": {},
  362. "outputs": [
  363. {
  364. "name": "stdout",
  365. "output_type": "stream",
  366. "text": [
  367. "\\left[\\begin{matrix}- \\cos{\\left(\\beta \\right)}\\\\- \\sin{\\left(\\beta \\right)}\\\\0\\end{matrix}\\right]\n"
  368. ]
  369. }
  370. ],
  371. "source": [
  372. "print(sp.latex(sol_b))"
  373. ]
  374. },
  375. {
  376. "cell_type": "markdown",
  377. "metadata": {},
  378. "source": [
  379. "# Key result - Formula for the Solar Angle"
  380. ]
  381. },
  382. {
  383. "cell_type": "code",
  384. "execution_count": 18,
  385. "metadata": {},
  386. "outputs": [
  387. {
  388. "data": {
  389. "text/latex": [
  390. "$\\displaystyle \\sin{\\left(\\alpha \\right)} \\sin{\\left(\\beta \\right)} \\cos{\\left(\\phi \\right)} - \\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} \\cos{\\left(\\beta \\right)} - \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\beta \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}$"
  391. ],
  392. "text/plain": [
  393. "sin(alpha)*sin(beta)*cos(phi) - sin(epsilon)*sin(phi)*cos(beta) - cos(alpha)*cos(beta)*cos(epsilon)*cos(phi)"
  394. ]
  395. },
  396. "execution_count": 18,
  397. "metadata": {},
  398. "output_type": "execute_result"
  399. }
  400. ],
  401. "source": [
  402. "# SA_ = Solar Angle\n",
  403. "SA_cos = N_ab.dot(sol_b)\n",
  404. "sp.expand(SA_cos)"
  405. ]
  406. },
  407. {
  408. "cell_type": "code",
  409. "execution_count": 19,
  410. "metadata": {},
  411. "outputs": [
  412. {
  413. "name": "stdout",
  414. "output_type": "stream",
  415. "text": [
  416. "\\sin{\\left(\\alpha \\right)} \\sin{\\left(\\beta \\right)} \\cos{\\left(\\phi \\right)} - \\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} \\cos{\\left(\\beta \\right)} - \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\beta \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}\n"
  417. ]
  418. }
  419. ],
  420. "source": [
  421. "print(sp.latex(sp.expand(SA_cos)))"
  422. ]
  423. },
  424. {
  425. "cell_type": "code",
  426. "execution_count": 20,
  427. "metadata": {},
  428. "outputs": [],
  429. "source": [
  430. "def cos_to_horizon_deg(exp):\n",
  431. " return 90 - 180*sp.acos(exp)/sp.pi"
  432. ]
  433. },
  434. {
  435. "cell_type": "markdown",
  436. "metadata": {},
  437. "source": [
  438. "# Now create Desmos visualizations"
  439. ]
  440. },
  441. {
  442. "cell_type": "markdown",
  443. "metadata": {},
  444. "source": [
  445. "### First visualization with simple time implementation"
  446. ]
  447. },
  448. {
  449. "cell_type": "code",
  450. "execution_count": 21,
  451. "metadata": {},
  452. "outputs": [],
  453. "source": [
  454. "# Reserved for 'time'\n",
  455. "t = sp.symbols('t')"
  456. ]
  457. },
  458. {
  459. "cell_type": "code",
  460. "execution_count": 22,
  461. "metadata": {},
  462. "outputs": [
  463. {
  464. "name": "stdout",
  465. "output_type": "stream",
  466. "text": [
  467. "90 - 57.2956455309397*arccos(-(sin(0.0174533333333333*L)*sin(0.0174533333333333*p) + cos(0.0174533333333333*L)*cos(0.0174533333333333*p)*cos(6.2832*x/H))*cos(6.2832*x/(H*Y)) + sin(6.2832*x/H)*sin(6.2832*x/(H*Y))*cos(0.0174533333333333*L))\n"
  468. ]
  469. }
  470. ],
  471. "source": [
  472. "# Desmos graph 1\n",
  473. "\n",
  474. "# Single letter vars for Desmos\n",
  475. "L = sp.symbols('L') \n",
  476. "p = sp.symbols('p')\n",
  477. "H = sp.symbols('H') \n",
  478. "Y = sp.symbols('Y') \n",
  479. "x = sp.symbols('x') \n",
  480. "\n",
  481. "\n",
  482. "# 't' will be in units of 'hours'\n",
  483. "def desmosify1(e):\n",
  484. " e1 = e.subs([\n",
  485. " (eps, 2*sp.pi*p*sp.Rational(1,360)),\n",
  486. " (phi, sp.pi*L*sp.Rational(1,180)),\n",
  487. " (alpha, 2*sp.pi*t/H),\n",
  488. " (beta, 2*sp.pi * t/(H*Y))\n",
  489. " ])\n",
  490. " \n",
  491. " # Rename t (time) as x, as Desmos wants to plot y vs. x\n",
  492. " # Also, Desmos doesn't like my 'pi'\n",
  493. " e2 = e1.subs(t,x).subs(sp.pi, 3.1416)\n",
  494. " \n",
  495. " return str(e2).replace('acos', 'arccos')\n",
  496. "\n",
  497. "print(desmosify1(cos_to_horizon_deg(SA_cos)))"
  498. ]
  499. },
  500. {
  501. "cell_type": "markdown",
  502. "metadata": {},
  503. "source": [
  504. "### Visualizations 2 & 3 - Solar height over course of a day"
  505. ]
  506. },
  507. {
  508. "cell_type": "code",
  509. "execution_count": 23,
  510. "metadata": {},
  511. "outputs": [],
  512. "source": [
  513. "# Desmos graphs 2 & 3\n",
  514. "\n",
  515. "# beta is determined by a slider for day of yer\n",
  516. "# beta: 0 -> 2*pi to d: 0 -> 365\n",
  517. "\n",
  518. "# alpha is should already have been converted to 't' (in hours)\n",
  519. "# in the passed in expression.\n",
  520. "\n",
  521. "# Single letter vars for Desmos\n",
  522. "L = sp.symbols('L') \n",
  523. "d = sp.symbols('d')\n",
  524. "x = sp.symbols('x') # for Desmos\n",
  525. "\n",
  526. "def desmosify23(e):\n",
  527. " e1 = e.subs([\n",
  528. " (eps, 2*sp.pi*23.44*sp.Rational(1,360)),\n",
  529. " (phi, sp.pi*L*sp.Rational(1,180)),\n",
  530. " (beta, 2*sp.pi * d*sp.Rational(1,365))\n",
  531. " ])\n",
  532. " \n",
  533. " # Rename t (time) as x, as Desmos wants to plot y vs. x\n",
  534. " # Also, Desmos doesn't like my 'pi'\n",
  535. " e2 = e1.subs(t,x).subs(sp.pi, 3.1416)\n",
  536. " \n",
  537. " return str(e2).replace('acos', 'arccos')\n"
  538. ]
  539. },
  540. {
  541. "cell_type": "code",
  542. "execution_count": 24,
  543. "metadata": {},
  544. "outputs": [
  545. {
  546. "name": "stdout",
  547. "output_type": "stream",
  548. "text": [
  549. "90 - 57.2956455309397*arccos(-(0.397789385116828*sin(0.0174533333333333*L) + 0.917476759971813*cos(0.0174533333333333*L)*cos(0.2618*x))*cos(0.0172142465753425*d) + sin(0.0172142465753425*d)*sin(0.2618*x)*cos(0.0174533333333333*L))\n"
  550. ]
  551. }
  552. ],
  553. "source": [
  554. "SA_daily = SA_cos.subs([\n",
  555. " (alpha, 2*sp.pi*t*sp.Rational(1,24))\n",
  556. "])\n",
  557. "\n",
  558. "print(desmosify23(cos_to_horizon_deg(SA_daily)))"
  559. ]
  560. },
  561. {
  562. "cell_type": "code",
  563. "execution_count": 53,
  564. "metadata": {},
  565. "outputs": [
  566. {
  567. "data": {
  568. "text/latex": [
  569. "$\\displaystyle - \\left(\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} + \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}\\right) \\cos{\\left(\\beta \\right)} + \\sin{\\left(\\alpha \\right)} \\sin{\\left(\\beta \\right)} \\cos{\\left(\\phi \\right)}$"
  570. ],
  571. "text/plain": [
  572. "-(sin(epsilon)*sin(phi) + cos(alpha)*cos(epsilon)*cos(phi))*cos(beta) + sin(alpha)*sin(beta)*cos(phi)"
  573. ]
  574. },
  575. "execution_count": 53,
  576. "metadata": {},
  577. "output_type": "execute_result"
  578. }
  579. ],
  580. "source": [
  581. "SA_daily\n",
  582. "SA_cos"
  583. ]
  584. },
  585. {
  586. "cell_type": "code",
  587. "execution_count": 54,
  588. "metadata": {},
  589. "outputs": [
  590. {
  591. "name": "stdout",
  592. "output_type": "stream",
  593. "text": [
  594. "- \\left(\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} + \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}\\right) \\cos{\\left(\\beta \\right)} + \\sin{\\left(\\alpha \\right)} \\sin{\\left(\\beta \\right)} \\cos{\\left(\\phi \\right)}\n"
  595. ]
  596. }
  597. ],
  598. "source": [
  599. "print(sp.latex(SA_cos))"
  600. ]
  601. },
  602. {
  603. "cell_type": "code",
  604. "execution_count": 27,
  605. "metadata": {},
  606. "outputs": [
  607. {
  608. "name": "stdout",
  609. "output_type": "stream",
  610. "text": [
  611. "- \\frac{180 \\operatorname{acos}{\\left(- \\left(\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} + \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)} \\cos{\\left(\\frac{\\pi t}{12} \\right)}\\right) \\cos{\\left(\\beta \\right)} + \\sin{\\left(\\beta \\right)} \\sin{\\left(\\frac{\\pi t}{12} \\right)} \\cos{\\left(\\phi \\right)} \\right)}}{\\pi} + 90\n"
  612. ]
  613. }
  614. ],
  615. "source": [
  616. "print(sp.latex(cos_to_horizon_deg(SA_daily)))"
  617. ]
  618. },
  619. {
  620. "cell_type": "markdown",
  621. "metadata": {},
  622. "source": [
  623. "### Sidereal 'cheat'"
  624. ]
  625. },
  626. {
  627. "cell_type": "code",
  628. "execution_count": 28,
  629. "metadata": {},
  630. "outputs": [
  631. {
  632. "data": {
  633. "text/latex": [
  634. "$\\displaystyle - \\left(\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} + \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)} \\cos{\\left(\\beta - \\frac{\\pi t}{12} \\right)}\\right) \\cos{\\left(\\beta \\right)} - \\sin{\\left(\\beta \\right)} \\sin{\\left(\\beta - \\frac{\\pi t}{12} \\right)} \\cos{\\left(\\phi \\right)}$"
  635. ],
  636. "text/plain": [
  637. "-(sin(epsilon)*sin(phi) + cos(epsilon)*cos(phi)*cos(beta - pi*t/12))*cos(beta) - sin(beta)*sin(beta - pi*t/12)*cos(phi)"
  638. ]
  639. },
  640. "execution_count": 28,
  641. "metadata": {},
  642. "output_type": "execute_result"
  643. }
  644. ],
  645. "source": [
  646. "# Now with the Sidereal cheat\n",
  647. "\n",
  648. "N_ab_sidereal = M_tilt * M_rot.subs(alpha, alpha - beta) * N\n",
  649. "SA_cos_sidereal = N_ab_sidereal.dot(sol_b)\n",
  650. "\n",
  651. "SA_daily_sidereal = SA_cos_sidereal.subs([\n",
  652. " #(alpha, 2*sp.pi*t/24.0)\n",
  653. " (alpha, 2*sp.pi*t*sp.Rational(1,24))\n",
  654. "])\n",
  655. "\n",
  656. "SA_daily_sidereal"
  657. ]
  658. },
  659. {
  660. "cell_type": "code",
  661. "execution_count": 29,
  662. "metadata": {},
  663. "outputs": [
  664. {
  665. "name": "stdout",
  666. "output_type": "stream",
  667. "text": [
  668. "- \\left(\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} + \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)} \\cos{\\left(\\beta - \\frac{\\pi t}{12} \\right)}\\right) \\cos{\\left(\\beta \\right)} - \\sin{\\left(\\beta \\right)} \\sin{\\left(\\beta - \\frac{\\pi t}{12} \\right)} \\cos{\\left(\\phi \\right)}\n"
  669. ]
  670. }
  671. ],
  672. "source": [
  673. "print(sp.latex(SA_daily_sidereal))"
  674. ]
  675. },
  676. {
  677. "cell_type": "code",
  678. "execution_count": 30,
  679. "metadata": {},
  680. "outputs": [
  681. {
  682. "name": "stdout",
  683. "output_type": "stream",
  684. "text": [
  685. "90 - 57.2956455309397*arccos(-(0.397789385116828*sin(0.0174533333333333*L) + 0.917476759971813*cos(0.0174533333333333*L)*cos(0.0172142465753425*d - 0.2618*x))*cos(0.0172142465753425*d) - sin(0.0172142465753425*d)*sin(0.0172142465753425*d - 0.2618*x)*cos(0.0174533333333333*L))\n"
  686. ]
  687. }
  688. ],
  689. "source": [
  690. "print(desmosify23(cos_to_horizon_deg(SA_daily_sidereal)))"
  691. ]
  692. },
  693. {
  694. "cell_type": "code",
  695. "execution_count": 31,
  696. "metadata": {},
  697. "outputs": [
  698. {
  699. "name": "stdout",
  700. "output_type": "stream",
  701. "text": [
  702. "- \\frac{180 \\operatorname{acos}{\\left(- \\left(\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} + \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)} \\cos{\\left(\\beta - \\frac{\\pi t}{12} \\right)}\\right) \\cos{\\left(\\beta \\right)} - \\sin{\\left(\\beta \\right)} \\sin{\\left(\\beta - \\frac{\\pi t}{12} \\right)} \\cos{\\left(\\phi \\right)} \\right)}}{\\pi} + 90\n"
  703. ]
  704. }
  705. ],
  706. "source": [
  707. "print(sp.latex(cos_to_horizon_deg(SA_daily_sidereal)))"
  708. ]
  709. },
  710. {
  711. "cell_type": "markdown",
  712. "metadata": {},
  713. "source": [
  714. "## Formula for length of daylight"
  715. ]
  716. },
  717. {
  718. "cell_type": "code",
  719. "execution_count": 57,
  720. "metadata": {},
  721. "outputs": [
  722. {
  723. "name": "stdout",
  724. "output_type": "stream",
  725. "text": [
  726. "\\sin{\\left(\\alpha \\right)} \\sin{\\left(\\beta \\right)} \\cos{\\left(\\phi \\right)} - \\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} \\cos{\\left(\\beta \\right)} - \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\beta \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}\n"
  727. ]
  728. }
  729. ],
  730. "source": [
  731. "sp.expand(SA_cos)\n",
  732. "print(sp.latex(sp.expand(SA_cos)))"
  733. ]
  734. },
  735. {
  736. "cell_type": "code",
  737. "execution_count": 60,
  738. "metadata": {},
  739. "outputs": [
  740. {
  741. "data": {
  742. "text/plain": [
  743. "(2*atan((A - sqrt(A**2 + B**2 - C**2))/(B - C)),\n",
  744. " 2*atan((A + sqrt(A**2 + B**2 - C**2))/(B - C)))"
  745. ]
  746. },
  747. "execution_count": 60,
  748. "metadata": {},
  749. "output_type": "execute_result"
  750. }
  751. ],
  752. "source": [
  753. "A, B, C, w = sp.symbols('A B C w')\n",
  754. "\n",
  755. "# w is pi*t/12\n",
  756. "\n",
  757. "f = A*sp.sin(w) + B*sp.cos(w) + C\n",
  758. "# This call to sp.solve() is fairly slow. Give it 20-30 seconds to run\n",
  759. "(sunrise, sunset) = sp.solve(f, w)\n",
  760. "sunrise, sunset"
  761. ]
  762. },
  763. {
  764. "cell_type": "code",
  765. "execution_count": 34,
  766. "metadata": {},
  767. "outputs": [
  768. {
  769. "data": {
  770. "text/latex": [
  771. "$\\displaystyle 2 \\operatorname{atan}{\\left(\\frac{A - \\sqrt{A^{2} + B^{2} - C^{2}}}{B - C} \\right)}$"
  772. ],
  773. "text/plain": [
  774. "2*atan((A - sqrt(A**2 + B**2 - C**2))/(B - C))"
  775. ]
  776. },
  777. "execution_count": 34,
  778. "metadata": {},
  779. "output_type": "execute_result"
  780. }
  781. ],
  782. "source": [
  783. "sunrise"
  784. ]
  785. },
  786. {
  787. "cell_type": "code",
  788. "execution_count": 58,
  789. "metadata": {},
  790. "outputs": [
  791. {
  792. "name": "stdout",
  793. "output_type": "stream",
  794. "text": [
  795. "2 \\operatorname{atan}{\\left(\\frac{A - \\sqrt{A^{2} + B^{2} - C^{2}}}{B - C} \\right)}\n"
  796. ]
  797. }
  798. ],
  799. "source": [
  800. "print(sp.latex(sunrise))"
  801. ]
  802. },
  803. {
  804. "cell_type": "code",
  805. "execution_count": 59,
  806. "metadata": {},
  807. "outputs": [
  808. {
  809. "name": "stdout",
  810. "output_type": "stream",
  811. "text": [
  812. "2 \\operatorname{atan}{\\left(\\frac{A + \\sqrt{A^{2} + B^{2} - C^{2}}}{B - C} \\right)}\n"
  813. ]
  814. }
  815. ],
  816. "source": [
  817. "print(sp.latex(sunset))"
  818. ]
  819. },
  820. {
  821. "cell_type": "code",
  822. "execution_count": 35,
  823. "metadata": {},
  824. "outputs": [
  825. {
  826. "data": {
  827. "text/latex": [
  828. "$\\displaystyle 2 \\operatorname{atan}{\\left(\\frac{A + \\sqrt{A^{2} + B^{2} - C^{2}}}{B - C} \\right)}$"
  829. ],
  830. "text/plain": [
  831. "2*atan((A + sqrt(A**2 + B**2 - C**2))/(B - C))"
  832. ]
  833. },
  834. "execution_count": 35,
  835. "metadata": {},
  836. "output_type": "execute_result"
  837. }
  838. ],
  839. "source": [
  840. "sunset"
  841. ]
  842. },
  843. {
  844. "cell_type": "code",
  845. "execution_count": 40,
  846. "metadata": {},
  847. "outputs": [
  848. {
  849. "name": "stdout",
  850. "output_type": "stream",
  851. "text": [
  852. "2 \\operatorname{atan}{\\left(\\frac{A - \\sqrt{A^{2} + B^{2} - C^{2}}}{B - C} \\right)}\n"
  853. ]
  854. }
  855. ],
  856. "source": [
  857. "print(sp.latex(sunrise))"
  858. ]
  859. },
  860. {
  861. "cell_type": "code",
  862. "execution_count": 47,
  863. "metadata": {},
  864. "outputs": [],
  865. "source": [
  866. "#A, B, C, w = sp.symbols('A B C w')\n",
  867. "\n",
  868. "# w is pi*t/12\n",
  869. "\n",
  870. "#f = A*sp.sin(w) + B*sp.cos(w) + C\n",
  871. "# Fairly slow\n",
  872. "#(sunrise, sunset) = sp.solve(f, w)\n",
  873. "lod_abc = (sunset-sunrise)\n",
  874. "lod_abc2= 2*sp.pi - (sunrise-sunset)\n"
  875. ]
  876. },
  877. {
  878. "cell_type": "code",
  879. "execution_count": 48,
  880. "metadata": {},
  881. "outputs": [
  882. {
  883. "data": {
  884. "text/latex": [
  885. "$\\displaystyle - 2 \\operatorname{atan}{\\left(\\frac{- \\sqrt{\\sin^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\phi \\right)} - \\sin^{2}{\\left(\\epsilon \\right)} \\sin^{2}{\\left(\\phi \\right)} \\cos^{2}{\\left(\\beta \\right)} + \\cos^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\epsilon \\right)} \\cos^{2}{\\left(\\phi \\right)}} + \\sin{\\left(\\beta \\right)} \\cos{\\left(\\phi \\right)}}{\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} \\cos{\\left(\\beta \\right)} - \\cos{\\left(\\beta \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}} \\right)} + 2 \\operatorname{atan}{\\left(\\frac{\\sqrt{\\sin^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\phi \\right)} - \\sin^{2}{\\left(\\epsilon \\right)} \\sin^{2}{\\left(\\phi \\right)} \\cos^{2}{\\left(\\beta \\right)} + \\cos^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\epsilon \\right)} \\cos^{2}{\\left(\\phi \\right)}} + \\sin{\\left(\\beta \\right)} \\cos{\\left(\\phi \\right)}}{\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} \\cos{\\left(\\beta \\right)} - \\cos{\\left(\\beta \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}} \\right)}$"
  886. ],
  887. "text/plain": [
  888. "-2*atan((-sqrt(sin(beta)**2*cos(phi)**2 - sin(epsilon)**2*sin(phi)**2*cos(beta)**2 + cos(beta)**2*cos(epsilon)**2*cos(phi)**2) + sin(beta)*cos(phi))/(sin(epsilon)*sin(phi)*cos(beta) - cos(beta)*cos(epsilon)*cos(phi))) + 2*atan((sqrt(sin(beta)**2*cos(phi)**2 - sin(epsilon)**2*sin(phi)**2*cos(beta)**2 + cos(beta)**2*cos(epsilon)**2*cos(phi)**2) + sin(beta)*cos(phi))/(sin(epsilon)*sin(phi)*cos(beta) - cos(beta)*cos(epsilon)*cos(phi)))"
  889. ]
  890. },
  891. "execution_count": 48,
  892. "metadata": {},
  893. "output_type": "execute_result"
  894. }
  895. ],
  896. "source": [
  897. "# Replace A, B, and C with functions of our parameters\n",
  898. "\n",
  899. "lod = lod_abc.subs([\n",
  900. " (A, sp.sin(beta) * sp.cos(phi)),\n",
  901. " (B, -sp.cos(beta) * sp.cos(eps) * sp.cos(phi)),\n",
  902. " (C, -sp.sin(eps) * sp.sin(phi) * sp.cos(beta))\n",
  903. "])\n",
  904. "\n",
  905. "lod2 = lod_abc2.subs([\n",
  906. " (A, sp.sin(beta) * sp.cos(phi)),\n",
  907. " (B, -sp.cos(beta) * sp.cos(eps) * sp.cos(phi)),\n",
  908. " (C, -sp.sin(eps) * sp.sin(phi) * sp.cos(beta))\n",
  909. "])\n",
  910. "\n",
  911. "lod"
  912. ]
  913. },
  914. {
  915. "cell_type": "code",
  916. "execution_count": 52,
  917. "metadata": {},
  918. "outputs": [
  919. {
  920. "data": {
  921. "text/latex": [
  922. "$\\displaystyle - 2 \\operatorname{atan}{\\left(\\frac{- \\sqrt{\\sin^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\phi \\right)} - \\sin^{2}{\\left(\\epsilon \\right)} \\sin^{2}{\\left(\\phi \\right)} \\cos^{2}{\\left(\\beta \\right)} + \\cos^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\epsilon \\right)} \\cos^{2}{\\left(\\phi \\right)}} + \\sin{\\left(\\beta \\right)} \\cos{\\left(\\phi \\right)}}{\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} \\cos{\\left(\\beta \\right)} - \\cos{\\left(\\beta \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}} \\right)} + 2 \\operatorname{atan}{\\left(\\frac{\\sqrt{\\sin^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\phi \\right)} - \\sin^{2}{\\left(\\epsilon \\right)} \\sin^{2}{\\left(\\phi \\right)} \\cos^{2}{\\left(\\beta \\right)} + \\cos^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\epsilon \\right)} \\cos^{2}{\\left(\\phi \\right)}} + \\sin{\\left(\\beta \\right)} \\cos{\\left(\\phi \\right)}}{\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} \\cos{\\left(\\beta \\right)} - \\cos{\\left(\\beta \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}} \\right)} + 2 \\pi$"
  923. ],
  924. "text/plain": [
  925. "-2*atan((-sqrt(sin(beta)**2*cos(phi)**2 - sin(epsilon)**2*sin(phi)**2*cos(beta)**2 + cos(beta)**2*cos(epsilon)**2*cos(phi)**2) + sin(beta)*cos(phi))/(sin(epsilon)*sin(phi)*cos(beta) - cos(beta)*cos(epsilon)*cos(phi))) + 2*atan((sqrt(sin(beta)**2*cos(phi)**2 - sin(epsilon)**2*sin(phi)**2*cos(beta)**2 + cos(beta)**2*cos(epsilon)**2*cos(phi)**2) + sin(beta)*cos(phi))/(sin(epsilon)*sin(phi)*cos(beta) - cos(beta)*cos(epsilon)*cos(phi))) + 2*pi"
  926. ]
  927. },
  928. "execution_count": 52,
  929. "metadata": {},
  930. "output_type": "execute_result"
  931. }
  932. ],
  933. "source": [
  934. "lod2"
  935. ]
  936. },
  937. {
  938. "cell_type": "code",
  939. "execution_count": 49,
  940. "metadata": {},
  941. "outputs": [
  942. {
  943. "name": "stdout",
  944. "output_type": "stream",
  945. "text": [
  946. "- 2 \\operatorname{atan}{\\left(\\frac{- \\sqrt{\\sin^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\phi \\right)} - \\sin^{2}{\\left(\\epsilon \\right)} \\sin^{2}{\\left(\\phi \\right)} \\cos^{2}{\\left(\\beta \\right)} + \\cos^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\epsilon \\right)} \\cos^{2}{\\left(\\phi \\right)}} + \\sin{\\left(\\beta \\right)} \\cos{\\left(\\phi \\right)}}{\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} \\cos{\\left(\\beta \\right)} - \\cos{\\left(\\beta \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}} \\right)} + 2 \\operatorname{atan}{\\left(\\frac{\\sqrt{\\sin^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\phi \\right)} - \\sin^{2}{\\left(\\epsilon \\right)} \\sin^{2}{\\left(\\phi \\right)} \\cos^{2}{\\left(\\beta \\right)} + \\cos^{2}{\\left(\\beta \\right)} \\cos^{2}{\\left(\\epsilon \\right)} \\cos^{2}{\\left(\\phi \\right)}} + \\sin{\\left(\\beta \\right)} \\cos{\\left(\\phi \\right)}}{\\sin{\\left(\\epsilon \\right)} \\sin{\\left(\\phi \\right)} \\cos{\\left(\\beta \\right)} - \\cos{\\left(\\beta \\right)} \\cos{\\left(\\epsilon \\right)} \\cos{\\left(\\phi \\right)}} \\right)}\n"
  947. ]
  948. }
  949. ],
  950. "source": [
  951. "print(sp.latex(lod))"
  952. ]
  953. },
  954. {
  955. "cell_type": "markdown",
  956. "metadata": {},
  957. "source": [
  958. "https://www.desmos.com/calculator/e5ta8kejtt"
  959. ]
  960. },
  961. {
  962. "cell_type": "code",
  963. "execution_count": 50,
  964. "metadata": {},
  965. "outputs": [
  966. {
  967. "name": "stdout",
  968. "output_type": "stream",
  969. "text": [
  970. "- 7.63941940412529 \\operatorname{arctan}{\\left(\\frac{- \\sqrt{- \\sin^{2}{\\left(0.0174533333333333 L \\right)} \\sin^{2}{\\left(0.0174533333333333 p \\right)} \\cos^{2}{\\left(\\frac{6.2832 x}{Y} \\right)} + \\sin^{2}{\\left(\\frac{6.2832 x}{Y} \\right)} \\cos^{2}{\\left(0.0174533333333333 L \\right)} + \\cos^{2}{\\left(0.0174533333333333 L \\right)} \\cos^{2}{\\left(0.0174533333333333 p \\right)} \\cos^{2}{\\left(\\frac{6.2832 x}{Y} \\right)}} + \\sin{\\left(\\frac{6.2832 x}{Y} \\right)} \\cos{\\left(0.0174533333333333 L \\right)}}{\\sin{\\left(0.0174533333333333 L \\right)} \\sin{\\left(0.0174533333333333 p \\right)} \\cos{\\left(\\frac{6.2832 x}{Y} \\right)} - \\cos{\\left(0.0174533333333333 L \\right)} \\cos{\\left(0.0174533333333333 p \\right)} \\cos{\\left(\\frac{6.2832 x}{Y} \\right)}} \\right)} + 7.63941940412529 \\operatorname{arctan}{\\left(\\frac{\\sqrt{- \\sin^{2}{\\left(0.0174533333333333 L \\right)} \\sin^{2}{\\left(0.0174533333333333 p \\right)} \\cos^{2}{\\left(\\frac{6.2832 x}{Y} \\right)} + \\sin^{2}{\\left(\\frac{6.2832 x}{Y} \\right)} \\cos^{2}{\\left(0.0174533333333333 L \\right)} + \\cos^{2}{\\left(0.0174533333333333 L \\right)} \\cos^{2}{\\left(0.0174533333333333 p \\right)} \\cos^{2}{\\left(\\frac{6.2832 x}{Y} \\right)}} + \\sin{\\left(\\frac{6.2832 x}{Y} \\right)} \\cos{\\left(0.0174533333333333 L \\right)}}{\\sin{\\left(0.0174533333333333 L \\right)} \\sin{\\left(0.0174533333333333 p \\right)} \\cos{\\left(\\frac{6.2832 x}{Y} \\right)} - \\cos{\\left(0.0174533333333333 L \\right)} \\cos{\\left(0.0174533333333333 p \\right)} \\cos{\\left(\\frac{6.2832 x}{Y} \\right)}} \\right)}\n"
  971. ]
  972. }
  973. ],
  974. "source": [
  975. "# Desmosify lod\n",
  976. "L = sp.symbols('L') \n",
  977. "p = sp.symbols('p')\n",
  978. "Y = sp.symbols('Y') \n",
  979. "d = sp.symbols('d') \n",
  980. "x = sp.symbols('x') \n",
  981. "\n",
  982. "def desmosify_lod( expr ):\n",
  983. " t0 = expr*12/sp.pi # take care of the pi/12, i.e. rescale to 24 hour day\n",
  984. "\n",
  985. " t1 = t0.subs([\n",
  986. " (eps, 2*sp.pi*p*sp.Rational(1,360)),\n",
  987. " (phi, sp.pi*L*sp.Rational(1,180)),\n",
  988. " (beta, 2*sp.pi * d/Y)\n",
  989. " ])\n",
  990. "\n",
  991. " t2 = t1.subs(d,x).subs(sp.pi, 3.1416)\n",
  992. " t3 = sp.latex(t2).replace(\"**\", \"^\").replace(\"atan\", \"arctan\")\n",
  993. " return t3\n",
  994. "\n",
  995. "print(desmosify_lod(lod))"
  996. ]
  997. },
  998. {
  999. "cell_type": "code",
  1000. "execution_count": 51,
  1001. "metadata": {},
  1002. "outputs": [
  1003. {
  1004. "name": "stdout",
  1005. "output_type": "stream",
  1006. "text": [
  1007. "- 7.63941940412529 \\operatorname{arctan}{\\left(\\frac{- \\sqrt{- \\sin^{2}{\\left(0.0174533333333333 L \\right)} \\sin^{2}{\\left(0.0174533333333333 p \\right)} \\cos^{2}{\\left(\\frac{6.2832 x}{Y} \\right)} + \\sin^{2}{\\left(\\frac{6.2832 x}{Y} \\right)} \\cos^{2}{\\left(0.0174533333333333 L \\right)} + \\cos^{2}{\\left(0.0174533333333333 L \\right)} \\cos^{2}{\\left(0.0174533333333333 p \\right)} \\cos^{2}{\\left(\\frac{6.2832 x}{Y} \\right)}} + \\sin{\\left(\\frac{6.2832 x}{Y} \\right)} \\cos{\\left(0.0174533333333333 L \\right)}}{\\sin{\\left(0.0174533333333333 L \\right)} \\sin{\\left(0.0174533333333333 p \\right)} \\cos{\\left(\\frac{6.2832 x}{Y} \\right)} - \\cos{\\left(0.0174533333333333 L \\right)} \\cos{\\left(0.0174533333333333 p \\right)} \\cos{\\left(\\frac{6.2832 x}{Y} \\right)}} \\right)} + 7.63941940412529 \\operatorname{arctan}{\\left(\\frac{\\sqrt{- \\sin^{2}{\\left(0.0174533333333333 L \\right)} \\sin^{2}{\\left(0.0174533333333333 p \\right)} \\cos^{2}{\\left(\\frac{6.2832 x}{Y} \\right)} + \\sin^{2}{\\left(\\frac{6.2832 x}{Y} \\right)} \\cos^{2}{\\left(0.0174533333333333 L \\right)} + \\cos^{2}{\\left(0.0174533333333333 L \\right)} \\cos^{2}{\\left(0.0174533333333333 p \\right)} \\cos^{2}{\\left(\\frac{6.2832 x}{Y} \\right)}} + \\sin{\\left(\\frac{6.2832 x}{Y} \\right)} \\cos{\\left(0.0174533333333333 L \\right)}}{\\sin{\\left(0.0174533333333333 L \\right)} \\sin{\\left(0.0174533333333333 p \\right)} \\cos{\\left(\\frac{6.2832 x}{Y} \\right)} - \\cos{\\left(0.0174533333333333 L \\right)} \\cos{\\left(0.0174533333333333 p \\right)} \\cos{\\left(\\frac{6.2832 x}{Y} \\right)}} \\right)} + 24.0\n"
  1008. ]
  1009. }
  1010. ],
  1011. "source": [
  1012. "print(desmosify_lod(lod2))"
  1013. ]
  1014. },
  1015. {
  1016. "cell_type": "markdown",
  1017. "metadata": {},
  1018. "source": [
  1019. "## Work area"
  1020. ]
  1021. },
  1022. {
  1023. "cell_type": "code",
  1024. "execution_count": null,
  1025. "metadata": {},
  1026. "outputs": [],
  1027. "source": [
  1028. "w, r, s = sp.symbols('w r s')\n",
  1029. "\n",
  1030. "f = sp.sin(s)*sp.sin(w) + 4*sp.cos(w) + 2\n",
  1031. "\n",
  1032. "# Incredibly slow\n",
  1033. "#sp.solve(f, w, exclude=[r,s])"
  1034. ]
  1035. },
  1036. {
  1037. "cell_type": "code",
  1038. "execution_count": null,
  1039. "metadata": {},
  1040. "outputs": [],
  1041. "source": [
  1042. "#sp.solve( SA_daily, t, exclude=[eps, phi, beta])\n",
  1043. "SA_daily.subs([(t, 0)])"
  1044. ]
  1045. },
  1046. {
  1047. "cell_type": "code",
  1048. "execution_count": null,
  1049. "metadata": {},
  1050. "outputs": [],
  1051. "source": [
  1052. "# simplified version of lod\n",
  1053. "b = sp.symbols('b')\n",
  1054. "lod_simplest = lod.subs([\n",
  1055. " (eps, 3.14 * 23.5/180.0),\n",
  1056. " (phi, 3.14*42.36/180.0),\n",
  1057. " (beta, 6.28*x/365.0)\n",
  1058. "])\n",
  1059. "\n",
  1060. "print(sp.latex(lod_simplest).replace(\"atan\", \"arctan\"))"
  1061. ]
  1062. },
  1063. {
  1064. "cell_type": "code",
  1065. "execution_count": null,
  1066. "metadata": {},
  1067. "outputs": [],
  1068. "source": [
  1069. "# simplified version of lod2\n",
  1070. "b = sp.symbols('b')\n",
  1071. "lod_simplest2 = lod2.subs([\n",
  1072. " (eps, 3.14 * 23.5/180.0),\n",
  1073. " (phi, 3.14*42.36/180.0),\n",
  1074. " (beta, 6.28*x/365.0)\n",
  1075. "])\n",
  1076. "\n",
  1077. "print(sp.latex(lod_simplest2).replace(\"atan\", \"arctan\"))"
  1078. ]
  1079. }
  1080. ],
  1081. "metadata": {
  1082. "kernelspec": {
  1083. "display_name": "Python 3",
  1084. "language": "python",
  1085. "name": "python3"
  1086. },
  1087. "language_info": {
  1088. "codemirror_mode": {
  1089. "name": "ipython",
  1090. "version": 3
  1091. },
  1092. "file_extension": ".py",
  1093. "mimetype": "text/x-python",
  1094. "name": "python",
  1095. "nbconvert_exporter": "python",
  1096. "pygments_lexer": "ipython3",
  1097. "version": "3.6.10"
  1098. }
  1099. },
  1100. "nbformat": 4,
  1101. "nbformat_minor": 4
  1102. }
  1103.  
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!