Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Fixed problem with planetocentric-specific parameters not being passe…
Browse files Browse the repository at this point in the history 
…d properly during an encounter. RMVS now works with oblateness.
  • Loading branch information
@daminton
daminton committed Jul 6, 2021
1 parent f8d52a1 commit 27dbf18
Show file tree
Hide file tree
Showing 3 changed files with 7 additions and 7 deletions.
4 changes: 2 additions & 2 deletions examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -482,7 +482,7 @@
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])\n",
"Coordinates:\n",
" id int64 4\n",
" * time (d) (time (d)) float64 0.0 1.0 2.0 3.0 4.0 ... 363.0 364.0 365.0 366.0</pre><div class='xr-wrap' hidden><div class='xr-header'><div class='xr-obj-type'>xarray.DataArray</div><div class='xr-array-name'>'px'</div><ul class='xr-dim-list'><li><span class='xr-has-index'>time (d)</span>: 367</li></ul></div><ul class='xr-sections'><li class='xr-section-item'><div class='xr-array-wrap'><input id='section-31768327-0ffd-429f-b0f5-581033c6d801' class='xr-array-in' type='checkbox' checked><label for='section-31768327-0ffd-429f-b0f5-581033c6d801' title='Show/hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-array-preview xr-preview'><span>0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0</span></div><div class='xr-array-data'><pre>array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" * time (d) (time (d)) float64 0.0 1.0 2.0 3.0 4.0 ... 363.0 364.0 365.0 366.0</pre><div class='xr-wrap' hidden><div class='xr-header'><div class='xr-obj-type'>xarray.DataArray</div><div class='xr-array-name'>'px'</div><ul class='xr-dim-list'><li><span class='xr-has-index'>time (d)</span>: 367</li></ul></div><ul class='xr-sections'><li class='xr-section-item'><div class='xr-array-wrap'><input id='section-fabb7d40-273a-4c06-a3fc-0399e88877ca' class='xr-array-in' type='checkbox' checked><label for='section-fabb7d40-273a-4c06-a3fc-0399e88877ca' title='Show/hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-array-preview xr-preview'><span>0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0</span></div><div class='xr-array-data'><pre>array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
Expand All @@ -503,7 +503,7 @@
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])</pre></div></div></li><li class='xr-section-item'><input id='section-a46df25f-b5e9-4e65-a77c-9b583e164594' class='xr-section-summary-in' type='checkbox' checked><label for='section-a46df25f-b5e9-4e65-a77c-9b583e164594' class='xr-section-summary' >Coordinates: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>id</span></div><div class='xr-var-dims'>()</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>4</div><input id='attrs-935ad6f8-92ab-4b29-8117-a352b238a841' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-935ad6f8-92ab-4b29-8117-a352b238a841' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-85807948-60a1-47e5-b247-9bf1e322b9a9' class='xr-var-data-in' type='checkbox'><label for='data-85807948-60a1-47e5-b247-9bf1e322b9a9' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array(4)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time (d)</span></div><div class='xr-var-dims'>(time (d))</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>0.0 1.0 2.0 ... 364.0 365.0 366.0</div><input id='attrs-967ca8e6-4c15-4896-8ec8-6e69d29dfd0c' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-967ca8e6-4c15-4896-8ec8-6e69d29dfd0c' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-35772f34-bf44-4cc6-9c73-0150e032e457' class='xr-var-data-in' type='checkbox'><label for='data-35772f34-bf44-4cc6-9c73-0150e032e457' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([ 0., 1., 2., ..., 364., 365., 366.])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-80315f9f-a629-4b1c-b8b1-e8c14c1bcc47' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-80315f9f-a629-4b1c-b8b1-e8c14c1bcc47' class='xr-section-summary' title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])</pre></div></div></li><li class='xr-section-item'><input id='section-de68d490-2721-4b5c-8b82-64e977e11680' class='xr-section-summary-in' type='checkbox' checked><label for='section-de68d490-2721-4b5c-8b82-64e977e11680' class='xr-section-summary' >Coordinates: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>id</span></div><div class='xr-var-dims'>()</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>4</div><input id='attrs-a358bd0d-ddc2-4083-90e7-b7cc9c5b89b4' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-a358bd0d-ddc2-4083-90e7-b7cc9c5b89b4' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-f7240ece-94ad-45ca-aa8a-2a52ac1d28df' class='xr-var-data-in' type='checkbox'><label for='data-f7240ece-94ad-45ca-aa8a-2a52ac1d28df' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array(4)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time (d)</span></div><div class='xr-var-dims'>(time (d))</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>0.0 1.0 2.0 ... 364.0 365.0 366.0</div><input id='attrs-4ce9e8f0-d411-4f33-9c83-325186b030b4' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-4ce9e8f0-d411-4f33-9c83-325186b030b4' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-237faf6f-1ced-40e0-a252-a3faef7506b9' class='xr-var-data-in' type='checkbox'><label for='data-237faf6f-1ced-40e0-a252-a3faef7506b9' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([ 0., 1., 2., ..., 364., 365., 366.])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-7f04a3a0-ecf9-4ee8-9c70-d473a2749c7f' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-7f04a3a0-ecf9-4ee8-9c70-d473a2749c7f' class='xr-section-summary' title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
],
"text/plain": [
"<xarray.DataArray 'px' (time (d): 367)>\n",
Expand Down
6 changes: 3 additions & 3 deletions src/obl/obl.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ module subroutine obl_acc_body(self, cb)
integer(I4B) :: i
real(DP) :: r2, irh, rinv2, t0, t1, t2, t3, fac1, fac2

associate(n => self%nbody)
associate(n => self%nbody, xh => self%xh, vh => self%vh, ah => self%ah)
do i = 1, n
r2 = dot_product(self%xh(:, i), self%xh(:, i))
irh = 1.0_DP / sqrt(r2)
Expand All @@ -26,8 +26,8 @@ module subroutine obl_acc_body(self, cb)
t1 = 1.5_DP * cb%j2rp2
t2 = self%xh(3, i) * self%xh(3, i) * rinv2
t3 = 1.875_DP * cb%j4rp4 * rinv2
fac1 = t0 * (t1 - t3 - (5.0_DP * t1 - (14.0_DP - 21.0_DP * t2) * t3) * t2)
fac2 = 2.0_DP * t0 * (t1 - (2.0_DP - (14.0_DP * t2 / 3.0_DP)) * t3)
fac1 = t0 * (t1 - t3 - (5 * t1 - (14.0_DP - 21.0_DP * t2) * t3) * t2)
fac2 = 2 * t0 * (t1 - (2.0_DP - (14.0_DP * t2 / 3.0_DP)) * t3)
self%aobl(:, i) = fac1 * self%xh(:, i)
self%aobl(3, i) = fac2 * self%xh(3, i) + self%aobl(3, i)
end do
Expand Down
4 changes: 2 additions & 2 deletions src/rmvs/rmvs_getacch.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,15 +31,15 @@ module subroutine rmvs_getacch_tp(self, system, param, t, xhp)
class is (rmvs_pl)
select type (cb => system%cb)
class is (rmvs_cb)
associate(xpc => pl%xh, xpct => self%xh)
associate(xpc => pl%xh, xpct => self%xh, apct => self%ah)
allocate(xh_original, source=tp%xh)
param_planetocen = param
! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter
param_planetocen%loblatecb = .false.
param_planetocen%lextra_force = .false.
param_planetocen%lgr = .false.
! Now compute the planetocentric values of acceleration
call whm_getacch_tp(tp, system, param, t, xhp)
call whm_getacch_tp(tp, system, param_planetocen, t, xhp)

! Now compute any heliocentric values of acceleration
if (tp%lfirst) then
Expand Down

0 comments on commit 27dbf18

Please sign in to comment.