From 1c8d750a27e5ac5bce9ea652f10c29380b5a0aff Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 13 Jul 2021 07:21:39 -0400 Subject: [PATCH] Moved rotation and tide variables to base swiftest_pl and swiftest_cb classes --- src/io/io.f90 | 97 +++++++++-- src/modules/swiftest_classes.f90 | 12 ++ src/modules/symba_classes.f90 | 69 -------- src/symba/symba_io.f90 | 266 ------------------------------- 4 files changed, 95 insertions(+), 349 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index d04e10466..2bb46ae0e 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -120,7 +120,13 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) case ("GR") call io_toupper(param_value) if (param_value == "YES" .or. param_value == 'T') self%lgr = .true. - case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "ROTATION", "TIDES", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters + case ("ROTATION") + call io_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T') self%lrotation = .true. + case ("TIDES") + call io_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T') self%ltides = .true. + case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default write(iomsg,*) "Unknown parameter -> ",param_name iostat = -1 @@ -192,6 +198,11 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) return end if end if + if (self%ltides .and. .not. self%lrotation) then + write(iomsg,*) 'Tides require rotation to be turned on' + iostat = -1 + return + end if write(*,*) "T0 = ",self%t0 write(*,*) "TSTOP = ",self%tstop @@ -217,6 +228,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) write(*,*) "EXTRA_FORCE = ",self%lextra_force write(*,*) "BIG_DISCARD = ",self%lbig_discard write(*,*) "RHILL_PRESENT = ",self%lrhill_present + write(*,*) "ROTATION = ", self%lrotation + write(*,*) "TIDES = ", self%ltides write(*,*) "ENERGY = ",self%lenergy write(*,*) "MU2KG = ",self%MU2KG write(*,*) "TU2S = ",self%TU2S @@ -325,6 +338,8 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(param_name, Afmt) "CHK_CLOSE"; write(param_value, Lfmt) param%lclose; write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "ENERGY"; write(param_value, Lfmt) param%lenergy; write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "GR"; write(param_value, Lfmt) param%lgr; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "ROTATION"; write(param_value, Lfmt) param%lrotation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "TIDES"; write(param_value, Lfmt) param%ltides; write(unit, Afmt) adjustl(param_name), adjustl(param_value) iostat = 0 iomsg = "UDIO not implemented" end associate @@ -610,6 +625,14 @@ module subroutine io_read_body_in(self, param) else self%radius(i) = 0.0_DP end if + if (param%lrotation) then + read(iu, iostat = ierr) self%Ip(:, i) + read(iu, iostat = ierr) self%rot(:, i) + end if + if (param%ltides) then + read(iu, iostat = ierr) self%k2(i) + read(iu, iostat = ierr) self%Q(i) + end if class is (swiftest_tp) read(iu, *, iostat = ierr) self%id(i) end select @@ -811,11 +834,23 @@ module subroutine io_read_frame_body(self, iu, param, form, ierr) read(iu, iostat = ierr) self%vh(2, 1:n) read(iu, iostat = ierr) self%vh(3, 1:n) end select - select type(self) + select type(pl => self) class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body - read(iu, iostat = ierr) self%Gmass(1:n) - self%mass(1:n) = self%Gmass / param%GU - read(iu, iostat = ierr) self%radius(1:n) + read(iu, iostat = ierr) pl%Gmass(1:n) + pl%mass(1:n) = pl%Gmass / param%GU + read(iu, iostat = ierr) pl%radius(1:n) + if (param%lrotation) then + read(iu, iostat = ierr) pl%rot(1, 1:n) + read(iu, iostat = ierr) pl%rot(2, 1:n) + read(iu, iostat = ierr) pl%rot(3, 1:n) + read(iu, iostat = ierr) pl%Ip(1, 1:n) + read(iu, iostat = ierr) pl%Ip(2, 1:n) + read(iu, iostat = ierr) pl%Ip(3, 1:n) + end if + if (param%ltides) then + read(iu, iostat = ierr) pl%k2(1:n) + read(iu, iostat = ierr) pl%Q(1:n) + end if end select end associate @@ -849,6 +884,14 @@ module subroutine io_read_frame_cb(self, iu, param, form, ierr) read(iu, iostat = ierr) self%radius read(iu, iostat = ierr) self%j2rp2 read(iu, iostat = ierr) self%j4rp4 + if (param%lrotation) then + read(iu, iostat = ierr) self%Ip(:) + read(iu, iostat = ierr) self%rot(:) + end if + if (param%ltides) then + read(iu, iostat = ierr) self%k2 + read(iu, iostat = ierr) self%Q + end if if (ierr /=0) then write(*,*) 'Error reading central body data' call util_exit(FAILURE) @@ -1133,10 +1176,22 @@ module subroutine io_write_frame_body(self, iu, param) write(iu) self%vh(2, 1:n) write(iu) self%vh(3, 1:n) end select - select type(self) + select type(pl => self) class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body - write(iu) self%Gmass(1:n) - write(iu) self%radius(1:n) + write(iu) pl%Gmass(1:n) + write(iu) pl%radius(1:n) + if (param%lrotation) then + write(iu) pl%rot(1, 1:n) + write(iu) pl%rot(2, 1:n) + write(iu) pl%rot(3, 1:n) + write(iu) pl%Ip(1, 1:n) + write(iu) pl%Ip(2, 1:n) + write(iu) pl%Ip(3, 1:n) + end if + if (param%ltides) then + write(iu) pl%k2(1:n) + write(iu) pl%Q(1:n) + end if end select end associate @@ -1156,12 +1211,26 @@ module subroutine io_write_frame_cb(self, iu, param) integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - write(iu) self%id - !write(iu) self%name - write(iu) self%Gmass - write(iu) self%radius - write(iu) self%j2rp2 - write(iu) self%j4rp4 + associate(cb => self) + !write(iu) cb%name + write(iu) cb%id + write(iu) cb%Gmass + write(iu) cb%radius + write(iu) cb%j2rp2 + write(iu) cb%j4rp4 + if (param%lrotation) then + write(iu) cb%rot(1) + write(iu) cb%rot(2) + write(iu) cb%rot(3) + write(iu) cb%Ip(1) + write(iu) cb%Ip(2) + write(iu) cb%Ip(3) + end if + if (param%ltides) then + write(iu) cb%k2 + write(iu) cb%Q + end if + end associate return end subroutine io_write_frame_cb diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index a5885255c..8efe47fe9 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -67,6 +67,8 @@ module swiftest_classes logical :: lclose = .false. !! Turn on close encounters logical :: lenergy = .false. !! Track the total energy of the system logical :: loblatecb = .false. !! Calculate acceleration from oblate central body (automatically turns true if nonzero J2 is input) + logical :: lrotation = .false. !! Include rotation states of big bodies + logical :: ltides = .false. !! Include tidal dissipation ! Future features not implemented or in development logical :: lgr = .false. !! Turn on GR @@ -114,6 +116,12 @@ module swiftest_classes real(DP), dimension(NDIM) :: xb = 0.0_DP !! Barycentric position (units DU) real(DP), dimension(NDIM) :: vb = 0.0_DP !! Barycentric velocity (units DU / TU) real(DP), dimension(NDIM) :: agr = 0.0_DP !! Acceleration due to post-Newtonian correction + real(DP), dimension(NDIM) :: Ip = 0.0_DP !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. + real(DP), dimension(NDIM) :: rot = 0.0_DP !! Body rotation vector in inertial coordinate frame (units rad / TU) + real(DP) :: k2 = 0.0_DP !! Tidal Love number + real(DP) :: Q = 0.0_DP !! Tidal quality factor + real(DP), dimension(NDIM) :: L0 = 0.0_DP !! Initial angular momentum of the central body + real(DP), dimension(NDIM) :: dL = 0.0_DP !! Change in angular momentum of the central body contains private procedure, public :: initialize => io_read_cb_in !! I/O routine for reading in central body data @@ -192,6 +200,10 @@ module swiftest_classes real(DP), dimension(:,:), allocatable :: xend !! Position at end of step real(DP), dimension(:,:), allocatable :: vbeg !! Velocity at beginning of step real(DP), dimension(:), allocatable :: density !! Body mass density - calculated internally (units MU / DU**3) + real(DP), dimension(:,:), allocatable :: Ip !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. + real(DP), dimension(:,:), allocatable :: rot !! Body rotation vector in inertial coordinate frame (units rad / TU) + real(DP), dimension(:), allocatable :: k2 !! Tidal Love number + real(DP), dimension(:), allocatable :: Q !! Tidal quality factor !! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the !! component list, such as setup_pl and util_spill_pl contains diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 2dad53567..5b712c9de 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -20,8 +20,6 @@ module symba_classes real(DP) :: MTINY = -1.0_DP !! Smallest mass that is fully gravitating integer(I4B), dimension(:), allocatable :: seed !! Random seeds logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. - logical :: lrotation = .false. !! Include rotation states of big bodies - logical :: ltides = .false. !! Include tidal dissipation contains private procedure, public :: reader => symba_io_param_reader @@ -33,21 +31,12 @@ module symba_classes !******************************************************************************************************************************* !> SyMBA central body particle class type, public, extends(helio_cb) :: symba_cb - real(DP), dimension(NDIM) :: Ip = 0.0_DP !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. - real(DP), dimension(NDIM) :: rot = 0.0_DP !! Body rotation vector in inertial coordinate frame (units rad / TU) - real(DP) :: k2 = 0.0_DP !! Tidal Love number - real(DP) :: Q = 0.0_DP !! Tidal quality factor - real(DP), dimension(NDIM) :: L0 = 0.0_DP !! Initial angular momentum of the central body - real(DP), dimension(NDIM) :: dL = 0.0_DP !! Change in angular momentum of the central body real(DP) :: M0 = 0.0_DP !! Initial mass of the central body real(DP) :: dM = 0.0_DP !! Change in mass of the central body real(DP) :: R0 = 0.0_DP !! Initial radius of the central body real(DP) :: dR = 0.0_DP !! Change in the radius of the central body contains private - procedure, public :: initialize => symba_io_read_cb_in !! I/O routine for reading in particle info data - procedure, public :: read_frame => symba_io_read_frame_cb !! I/O routine for reading out a single frame of time-series data for the central bod - procedure, public :: write_frame => symba_io_write_frame_cb !! I/O routine for writing out a single frame of time-series data for the central body end type symba_cb !******************************************************************************************************************************** @@ -83,11 +72,6 @@ module symba_classes !******************************************************************************************************************************* !> SyMBA massive body class type, public, extends(helio_pl) :: symba_pl - real(DP), dimension(:,:), allocatable :: Ip !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). - !! Principal axis rotation assumed. - real(DP), dimension(:,:), allocatable :: rot !! Body rotation vector in inertial coordinate frame (units rad / TU) - real(DP), dimension(:), allocatable :: k2 !! Tidal Love number - real(DP), dimension(:), allocatable :: Q !! Tidal quality factor logical, dimension(:), allocatable :: lcollision !! flag indicating whether body has merged with another this time step logical, dimension(:), allocatable :: lencounter !! flag indicating whether body is part of an encounter this time step integer(I4B), dimension(:), allocatable :: nplenc !! number of encounters with other planets this time step @@ -103,9 +87,6 @@ module symba_classes private procedure, public :: discard => symba_discard_pl !! Process massive body discards procedure, public :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other - procedure, public :: read_frame => symba_io_read_frame_pl !! I/O routine for reading out a single frame of time-series data for a massive body - procedure, public :: initialize => symba_io_read_pl_in !! I/O routine for reading in a massive body structure from file with SyMBA-specific parameters - procedure, public :: write_frame => symba_io_write_frame_pl !! I/O routine for writing out a single frame of time-series data for a massive body procedure, public :: setup => symba_setup_pl !! Constructor method - Allocates space for number of particle end type symba_pl @@ -234,23 +215,6 @@ module subroutine symba_io_initialize_particle_info(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine symba_io_initialize_particle_info - module subroutine symba_io_read_cb_in(self, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_cb), intent(inout) :: self - class(swiftest_parameters), intent(inout) :: param - end subroutine symba_io_read_cb_in - - module subroutine symba_io_read_frame_cb(self, iu, param, form, ierr) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_cb), intent(inout) :: self !! Swiftest central body object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - integer(I4B), intent(out) :: ierr !! Error code - end subroutine symba_io_read_frame_cb - module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) use swiftest_classes, only : swiftest_parameters implicit none @@ -261,31 +225,6 @@ module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) integer(I4B), intent(out) :: ierr !! Error code end subroutine symba_io_read_frame_info - module subroutine symba_io_read_frame_pl(self, iu, param, form, ierr) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_pl), intent(inout) :: self !! Swiftest particle object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - integer(I4B), intent(out) :: ierr !! Error code - end subroutine symba_io_read_frame_pl - - module subroutine symba_io_read_pl_in(self, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_pl), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine symba_io_read_pl_in - - module subroutine symba_io_write_frame_cb(self, iu, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_cb), intent(in) :: self !! SyMBA massive body object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine symba_io_write_frame_cb - module subroutine symba_io_write_frame_info(self, iu, param) use swiftest_classes, only : swiftest_parameters implicit none @@ -294,14 +233,6 @@ module subroutine symba_io_write_frame_info(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_io_write_frame_info - module subroutine symba_io_write_frame_pl(self, iu, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_pl), intent(in) :: self !! SyMBA massive body object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine symba_io_write_frame_pl - module subroutine symba_setup_pl(self,n) implicit none class(symba_pl), intent(inout) :: self !! SyMBA test particle object diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 67c5f979d..bebb225b5 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -67,12 +67,6 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms case ("FRAGMENTATION") call io_toupper(param_value) if (param_value == "YES" .or. param_value == "T") self%lfragmentation = .true. - case ("ROTATION") - call io_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%lrotation = .true. - case ("TIDES") - call io_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%ltides = .true. case ("MTINY") read(param_value, *) param%mtiny case("SEED") @@ -112,8 +106,6 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms end if write(*,*) "SEED: N,VAL = ",size(param%seed), param%seed(:) end if - write(*,*) "ROTATION = ", param%lrotation - write(*,*) "TIDES = ", param%ltides if (self%mtiny < 0.0_DP) then write(iomsg,*) "MTINY invalid or not set: ", self%mtiny @@ -169,8 +161,6 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms ! For the "SEED" parameter line, the first value will be the size of the seed array and the rest will be the seed array elements write(param_name, Afmt) "PARTICLE_FILE"; write(param_value, Afmt) trim(adjustl(param%particle_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "MTINY"; write(param_value, Rfmt) param%mtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) - write(param_name, Afmt) "ROTATION"; write(param_value, Lfmt) param%lrotation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) - write(param_name, Afmt) "TIDES"; write(param_value, Lfmt) param%ltides; write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "FRAGMENTATION"; write(param_value, Lfmt) param%lfragmentation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) if (param%lfragmentation) then write(param_name, Afmt) "SEED" @@ -197,81 +187,6 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms end subroutine symba_io_param_writer - module subroutine symba_io_read_frame_cb(self, iu, param, form, ierr) - !! author: David A. Minton - !! - !! Reads a frame of output of central body data to the binary output file - !! - !! Adapted from David E. Kaufmann's Swifter routine io_read_frame.f90 - !! Adapted from Hal Levison's Swift routine io_read_frame.F - implicit none - ! Arguments - class(symba_cb), intent(inout) :: self !! Swiftest central body object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - integer(I4B), intent(out) :: ierr !! Error code - - call io_read_frame_cb(self, iu, param, form, ierr) - select type(param) - class is (symba_parameters) - if (param%lrotation) then - read(iu, iostat = ierr) self%Ip(:) - read(iu, iostat = ierr) self%rot(:) - end if - if (param%ltides) then - read(iu, iostat = ierr) self%k2 - read(iu, iostat = ierr) self%Q - end if - end select - if (ierr /=0) then - write(*,*) 'Error reading SyMBA central body data' - call util_exit(FAILURE) - end if - return - end subroutine symba_io_read_frame_cb - - module subroutine symba_io_read_frame_pl(self, iu, param, form, ierr) - !! author: David A. Minton - !! - !! Reads a frame of output of a SyMBA massive body object - !! - !! Adapted from David E. Kaufmann's Swifter routine io_read_frame.f90 - !! Adapted from Hal Levison's Swift routine io_read_frame.F - implicit none - ! Arguments - class(symba_pl), intent(inout) :: self !! Swiftest particle object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - integer(I4B), intent(out) :: ierr !! Error code - - call io_read_frame_body(self, iu, param, form, ierr) - select type(param) - class is (symba_parameters) - associate(pl => self, npl => self%nbody) - if (param%lrotation) then - read(iu, iostat = ierr) pl%rot(1, 1:npl) - read(iu, iostat = ierr) pl%rot(2, 1:npl) - read(iu, iostat = ierr) pl%rot(3, 1:npl) - read(iu, iostat = ierr) pl%Ip(1, 1:npl) - read(iu, iostat = ierr) pl%Ip(2, 1:npl) - read(iu, iostat = ierr) pl%Ip(3, 1:npl) - end if - if (param%ltides) then - read(iu, iostat = ierr) pl%k2(1:npl) - read(iu, iostat = ierr) pl%Q(1:npl) - end if - end associate - end select - - if (ierr /=0) then - write(*,*) 'Error reading SyMBA massive body body data' - call util_exit(FAILURE) - end if - return - end subroutine symba_io_read_frame_pl - module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) !! author: David A. Minton !! @@ -286,159 +201,6 @@ module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) ierr = 0 end subroutine symba_io_read_frame_info - module subroutine symba_io_read_cb_in(self, param) - !! author: David A. Minton - !! - !! Reads in central body data - !! - !! Adapted from David E. Kaufmann's Swifter routine swiftest_init_pl.f90 - !! Adapted from Martin Duncan's Swift routine swiftest_init_pl.f - implicit none - ! Arguments - class(symba_cb), intent(inout) :: self - class(swiftest_parameters), intent(inout) :: param - ! Internals - integer(I4B), parameter :: LUN = 7 !! Unit number of input file - integer(I4B) :: iu = LUN - integer(I4B) :: ierr - logical :: is_ascii - real(DP) :: t - real(QP) :: val - - select type(param) - class is (symba_parameters) - ierr = 0 - is_ascii = (param%in_type == 'ASCII') - if (is_ascii) then - open(unit = iu, file = param%incbfile, status = 'old', form = 'FORMATTED', iostat = ierr) - read(iu, *, iostat = ierr) val - self%Gmass = real(val, kind=DP) - self%mass = real(val / param%GU, kind=DP) - read(iu, *, iostat = ierr) self%radius - read(iu, *, iostat = ierr) self%j2rp2 - read(iu, *, iostat = ierr) self%j4rp4 - if (param%lrotation) then - read(iu, *, iostat = ierr) self%Ip(:) - read(iu, *, iostat = ierr) self%rot(:) - end if - if (param%ltides) then - read(iu, *, iostat = ierr) self%k2 - read(iu, *, iostat = ierr) self%Q - end if - else - open(unit = iu, file = param%incbfile, status = 'old', form = 'UNFORMATTED', iostat = ierr) - call self%read_frame(iu, param, XV, ierr) - end if - close(iu) - if (ierr /= 0) then - write(*,*) 'Error opening massive body initial conditions file ',trim(adjustl(param%incbfile)) - call util_exit(FAILURE) - end if - if (self%j2rp2 /= 0.0_DP) param%loblatecb = .true. - end select - - return - end subroutine symba_io_read_cb_in - - module subroutine symba_io_read_pl_in(self, param) - !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott - !! - !! Read in either test particle or massive body data - !! - !! Adapted from David E. Kaufmann's Swifter routine swiftest_init_pl.f90 and swiftest_init_tp.f90 - !! Adapted from Martin Duncan's Swift routine swiftest_init_pl.f and swiftest_init_tp.f - implicit none - ! Arguments - class(symba_pl), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B), parameter :: LUN = 7 !! Unit number of input file - integer(I4B) :: iu = LUN - integer(I4B) :: i, ierr, nbody - logical :: is_ascii - character(len=:), allocatable :: infile - real(DP) :: t - real(QP) :: val - - select type(param) - class is (symba_parameters) - ierr = 0 - is_ascii = (param%in_type == 'ASCII') - select case(param%in_type) - case(ASCII_TYPE) - open(unit = iu, file = infile, status = 'old', form = 'FORMATTED', iostat = ierr) - read(iu, *, iostat = ierr) nbody - call self%setup(nbody) - if (nbody > 0) then - do i = 1, nbody - read(iu, *, iostat = ierr) self%id(i), val - self%mass(i) = real(val / param%GU, kind=DP) - self%Gmass(i) = real(val, kind=DP) - read(iu, *, iostat = ierr) self%radius(i) - if (param%lrotation) then - read(iu, iostat = ierr) self%Ip(:, i) - read(iu, iostat = ierr) self%rot(:, i) - end if - if (param%ltides) then - read(iu, iostat = ierr) self%k2(i) - read(iu, iostat = ierr) self%Q(i) - end if - if (ierr /= 0 ) exit - read(iu, *, iostat = ierr) self%xh(1, i), self%xh(2, i), self%xh(3, i) - read(iu, *, iostat = ierr) self%vh(1, i), self%vh(2, i), self%vh(3, i) - if (ierr /= 0 ) exit - self%status(i) = ACTIVE - end do - end if - case (REAL4_TYPE, REAL8_TYPE) - open(unit = iu, file = infile, status = 'old', form = 'UNFORMATTED', iostat = ierr) - read(iu, iostat = ierr) nbody - call self%setup(nbody) - if (nbody > 0) then - call self%read_frame(iu, param, XV, ierr) - self%status(:) = ACTIVE - end if - case default - write(*,*) trim(adjustl(param%in_type)) // ' is an unrecognized file type' - ierr = -1 - end select - close(iu) - if (ierr /= 0 ) then - write(*,*) 'Error reading in initial conditions from ',trim(adjustl(infile)) - call util_exit(FAILURE) - end if - end select - - return - end subroutine symba_io_read_pl_in - - module subroutine symba_io_write_frame_cb(self, iu, param) - !! author: David A. Minton - !! - !! Writes a single frame of a SyMBA pl file - implicit none - class(symba_cb), intent(in) :: self !! SyMBA massive body object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - - call io_write_frame_cb(self, iu, param) - select type(param) - class is (symba_parameters) - if (param%lrotation) then - write(iu) self%rot(1) - write(iu) self%rot(2) - write(iu) self%rot(3) - write(iu) self%Ip(1) - write(iu) self%Ip(2) - write(iu) self%Ip(3) - end if - if (param%ltides) then - write(iu) self%k2 - write(iu) self%Q - end if - end select - end subroutine symba_io_write_frame_cb - module subroutine symba_io_write_frame_info(self, iu, param) implicit none class(symba_particle_info), intent(in) :: self !! SyMBA particle info object @@ -446,34 +208,6 @@ module subroutine symba_io_write_frame_info(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_io_write_frame_info - module subroutine symba_io_write_frame_pl(self, iu, param) - !! author: David A. Minton - !! - !! Writes a single frame of a SyMBA pl file - implicit none - class(symba_pl), intent(in) :: self !! SyMBA massive body object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - - call io_write_frame_body(self, iu, param) - select type(param) - class is (symba_parameters) - associate(pl => self, npl => self%nbody) - if (param%lrotation) then - write(iu) pl%rot(1, 1:npl) - write(iu) pl%rot(2, 1:npl) - write(iu) pl%rot(3, 1:npl) - write(iu) pl%Ip(1, 1:npl) - write(iu) pl%Ip(2, 1:npl) - write(iu) pl%Ip(3, 1:npl) - end if - if (param%ltides) then - write(iu) pl%k2(1:npl) - write(iu) pl%Q(1:npl) - end if - end associate - end select - end subroutine symba_io_write_frame_pl end submodule s_symba_io