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

Commit

Permalink
Added central body to output file
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Jun 22, 2021
1 parent 7c8e877 commit 117f8ce
Showing 1 changed file with 52 additions and 54 deletions.
106 changes: 52 additions & 54 deletions src/io/io_write_frame.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,10 @@
integer(I4B) :: i, j, ierr
integer(I4B), save :: iu = lun, iout_form = xv
real(DP) :: a, e, inc, capom, omega, capm, mu
real(DP), dimension(2:swiftest_plA%nbody) :: a_pl, e_pl, inc_pl, capom_pl, omega_pl, capm_pl
real(DP), dimension(1:swiftest_tpA%nbody) :: a_tp, e_tp, inc_tp, capom_tp, omega_tp, capm_tp
real(DP), dimension(swiftest_plA%nbody) :: a_pl, e_pl, inc_pl, capom_pl, omega_pl, capm_pl
real(DP), dimension(swiftest_tpA%nbody) :: a_tp, e_tp, inc_tp, capom_tp, omega_tp, capm_tp

associate(out_stat => param%out_stat, out_form => param%out_form, out_type => param%out_type, outfile => param%outfile)
associate(out_stat => param%out_stat, out_form => param%out_form, out_type => param%out_type, outfile => param%outfile, npl => swiftest_plA%nbody, ntp => swiftest_tpA%nbody)
if (lfirst) then
select case(out_stat)
case('APPEND')
Expand Down Expand Up @@ -55,10 +55,10 @@
end if
end if

call io_write_hdr(iu, t, swiftest_plA%nbody, swiftest_tpA%nbody, iout_form, out_type)
call io_write_hdr(iu, t, npl, ntp, iout_form, out_type)
select case (iout_form)
case (EL)
do i = 2, swiftest_plA%nbody
do i = 2, npl
mu = swiftest_plA%mass(1) + swiftest_plA%mass(i)
j = swiftest_plA%id(i)
call orbel_xv2el(swiftest_plA%xh(:,i), swiftest_plA%vh(:,i), mu, a, e, inc, capom, omega, capm)
Expand All @@ -70,7 +70,7 @@
capm_pl(i) = capm
end do
mu = swiftest_plA%mass(1)
do i = 1, swiftest_tpA%nbody
do i = 1, ntp
j = swiftest_tpA%id(i)
call orbel_xv2el(swiftest_tpA%xh(:,i), swiftest_tpA%vh(:,i), mu, a, e, inc, capom, omega, capm)
a_tp(i) = a
Expand All @@ -80,62 +80,60 @@
omega_tp(i) = omega
capm_tp(i) = capm
end do
write(LUN) swiftest_plA%id(2:swiftest_plA%nbody)
write(LUN) a_pl(2:swiftest_plA%nbody)
write(LUN) e_pl(2:swiftest_plA%nbody)
write(LUN) inc_pl(2:swiftest_plA%nbody)
write(LUN) capom_pl(2:swiftest_plA%nbody)
write(LUN) omega_pl(2:swiftest_plA%nbody)
write(LUN) capm_pl(2:swiftest_plA%nbody)
write(LUN) swiftest_plA%mass(2:swiftest_plA%nbody)
write(LUN) swiftest_plA%radius(2:swiftest_plA%nbody)
write(LUN) swiftest_plA%id(1:npl)
write(LUN) a_pl(1:npl)
write(LUN) e_pl(1:npl)
write(LUN) inc_pl(1:npl)
write(LUN) capom_pl(1:npl)
write(LUN) omega_pl(1:npl)
write(LUN) capm_pl(1:npl)
write(LUN) swiftest_plA%mass(1:npl)
write(LUN) swiftest_plA%radius(1:npl)
if (param%lrotation) then
write(LUN) swiftest_plA%rot(1,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%rot(2,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%rot(3,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%Ip(1,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%Ip(2,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%Ip(3,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%rot(1,1:npl)
write(LUN) swiftest_plA%rot(2,1:npl)
write(LUN) swiftest_plA%rot(3,1:npl)
write(LUN) swiftest_plA%Ip(1,1:npl)
write(LUN) swiftest_plA%Ip(2,1:npl)
write(LUN) swiftest_plA%Ip(3,1:npl)
end if
if (swiftest_tpA%nbody > 0) then
write(LUN) swiftest_tpA%id(1:swiftest_tpA%nbody)
write(LUN) a_tp(1:swiftest_tpA%nbody)
write(LUN) e_tp(1:swiftest_tpA%nbody)
write(LUN) inc_tp(1:swiftest_tpA%nbody)
write(LUN) capom_tp(1:swiftest_tpA%nbody)
write(LUN) omega_tp(1:swiftest_tpA%nbody)
write(LUN) capm_tp(1:swiftest_tpA%nbody)
if (ntp > 0) then
write(LUN) swiftest_tpA%id(1:ntp)
write(LUN) a_tp(1:ntp)
write(LUN) e_tp(1:ntp)
write(LUN) inc_tp(1:ntp)
write(LUN) capom_tp(1:ntp)
write(LUN) omega_tp(1:ntp)
write(LUN) capm_tp(1:ntp)
end if

case (XV)
write(LUN) swiftest_plA%id(2:swiftest_plA%nbody)
write(LUN) swiftest_plA%xh(1,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%xh(2,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%xh(3,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%vh(1,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%vh(2,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%vh(3,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%mass(2:swiftest_plA%nbody)
write(LUN) swiftest_plA%radius(2:swiftest_plA%nbody)
write(LUN) swiftest_plA%id(1:npl)
write(LUN) swiftest_plA%xh(1,1:npl)
write(LUN) swiftest_plA%xh(2,1:npl)
write(LUN) swiftest_plA%xh(3,1:npl)
write(LUN) swiftest_plA%vh(1,1:npl)
write(LUN) swiftest_plA%vh(2,1:npl)
write(LUN) swiftest_plA%vh(3,1:npl)
write(LUN) swiftest_plA%mass(1:npl)
write(LUN) swiftest_plA%radius(1:npl)
if (param%lrotation) then
write(LUN) swiftest_plA%rot(1,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%rot(2,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%rot(3,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%Ip(1,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%Ip(2,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%Ip(3,2:swiftest_plA%nbody)
write(LUN) swiftest_plA%rot(1,1:npl)
write(LUN) swiftest_plA%rot(2,1:npl)
write(LUN) swiftest_plA%rot(3,1:npl)
write(LUN) swiftest_plA%Ip(1,1:npl)
write(LUN) swiftest_plA%Ip(2,1:npl)
write(LUN) swiftest_plA%Ip(3,1:npl)
end if

if (swiftest_tpA%nbody > 0) then
write(LUN) swiftest_tpA%id(:)
write(LUN) swiftest_tpA%xh(1,:)
write(LUN) swiftest_tpA%xh(2,:)
write(LUN) swiftest_tpA%xh(3,:)
write(LUN) swiftest_tpA%vh(1,:)
write(LUN) swiftest_tpA%vh(2,:)
write(LUN) swiftest_tpA%vh(3,:)
if (ntp > 0) then
write(LUN) swiftest_tpA%id(1:ntp)
write(LUN) swiftest_tpA%xh(1,1:ntp)
write(LUN) swiftest_tpA%xh(2,1:ntp)
write(LUN) swiftest_tpA%xh(3,1:ntp)
write(LUN) swiftest_tpA%vh(1,1:ntp)
write(LUN) swiftest_tpA%vh(2,1:ntp)
write(LUN) swiftest_tpA%vh(3,1:ntp)
end if

end select

close(unit = iu, iostat = ierr)
Expand Down

0 comments on commit 117f8ce

Please sign in to comment.