diff --git a/Makefile.Defines b/Makefile.Defines index 7bb6ae32e..6b35a612e 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -61,17 +61,17 @@ OPTREPORT = -qopt-report=5 #gfortran flags GDEBUG = -g -Og -fbacktrace -fbounds-check GPRODUCTION = -O3 -GPAR = -fopenmp -ftree-parallelize-loops=4 +GPAR = -fopenmp #-ftree-parallelize-loops=4 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporaries #FFLAGS = $(IDEBUG) $(HEAPARR) $(SIMDVEC) $(PAR) -FFLAGS = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(PAR) $(SIMDVEC) $(HEAPARR) -FORTRAN = ifort +#FFLAGS = -init=snan,arrays -no-wrap-margin -O0 -g -traceback $(STRICTREAL) $(PAR) $(SIMDVEC) $(HEAPARR) +#FORTRAN = ifort #AR = xiar -#FORTRAN = gfortran -#FFLAGS = -ffree-line-length-none $(GPRODUCTION) $(GPAR) #$(GDEBUG) $(GMEM) +FORTRAN = gfortran +FFLAGS = -ffree-line-length-none $(GPAR) $(GPRODUCTION) #$(GMEM) AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/examples/in_form_el/cb.in b/examples/in_form_el/cb.in new file mode 100644 index 000000000..673a79459 --- /dev/null +++ b/examples/in_form_el/cb.in @@ -0,0 +1,5 @@ +0 +0.0002959122081920778 +0.004650467260962157 +4.7535806948127355e-12 +-2.2473967953572827e-18 diff --git a/examples/in_form_el/param.in b/examples/in_form_el/param.in new file mode 100644 index 000000000..409cd2d2b --- /dev/null +++ b/examples/in_form_el/param.in @@ -0,0 +1,35 @@ +! VERSION Swiftest parameter input from file param.in +T0 0.0 +TSTOP 365250000000.0 +DT 5.0 +ISTEP_OUT 7305000 +ISTEP_DUMP 7305000 +OUT_FORM EL +OUT_TYPE REAL8 +OUT_STAT REPLACE +IN_TYPE ASCII +PL_IN pl.in +TP_IN tp.in +CB_IN cb.in +BIN_OUT bin.dat +CHK_QMIN 0.004650467260962157 +CHK_RMIN 0.004650467260962157 +CHK_RMAX 1000.0 +CHK_EJECT 1000.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 0.004650467260962157 1000.0 +MU2KG 1.988409870698051e+30 +TU2S 86400.0 +DU2M 149597870700.0 +IN_FORM EL +ENC_OUT enc.dat +EXTRA_FORCE NO +DISCARD_OUT discard.out +BIG_DISCARD NO +CHK_CLOSE YES +RHILL_PRESENT YES +FRAGMENTATION NO +ROTATION NO +TIDES NO +ENERGY NO +GR YES diff --git a/examples/in_form_el/pl.in b/examples/in_form_el/pl.in new file mode 100644 index 000000000..55748b408 --- /dev/null +++ b/examples/in_form_el/pl.in @@ -0,0 +1,33 @@ +8 +1 4.9125474498983625e-11 0.0014751256791533402 +1.6306381826061646e-05 +0.38710383731080894 0.24071190040066506 8.640070276952633 +30.683987164281632 57.949871930463694 151.08868778597147 +2 7.243452483873647e-10 0.0067591154147628296 +4.0453784346544176e-05 +0.7233315774753523 0.02498272794314777 2.4772002857938187 +125.76075677093282 6.742571718202682 350.2972657114627 +3 8.997011382166019e-10 0.010044745408888887 +4.25875607065041e-05 +0.9999904186618421 0.05232751342856539 2.0722848768764717 +113.84624589363827 312.12453079628943 148.9860363695549 +4 9.549535102761465e-11 0.007246758104365299 +2.2657408050928896e-05 +1.5236293140305446 0.005382832120155795 4.468216946001687 +325.6098790099478 152.40711385323957 110.04936074450407 +5 2.825345908631355e-07 0.3552714410085414 +0.0004673261703049093 +5.203095539549931 0.03063821583296288 1.3661365597246853 +113.82182400586191 303.52091573107697 96.03226696213453 +6 8.459715183006416e-08 0.43765112940564116 +0.00038925687730393614 +9.579329287028964 0.08182512713593486 2.54624964218005 +100.32672664745415 306.8523855342038 64.14978871000847 +7 1.2920249163736674e-08 0.4695173932519177 +0.00016953449859497232 +19.185526866180975 0.02725895554029442 0.6468816262831063 +75.84341829565291 203.1828812930829 337.91763726094723 +8 1.5243589003230834e-08 0.7812638143548571 +0.00016458790412449367 +30.033406636357004 0.014808837377897958 1.0172902231424474 +112.31047122167676 61.6982073585126 219.66063286932612 diff --git a/examples/in_form_el/tp.in b/examples/in_form_el/tp.in new file mode 100644 index 000000000..48e97fb82 --- /dev/null +++ b/examples/in_form_el/tp.in @@ -0,0 +1,76 @@ +25 +1030101 +1.9298376430729873 0.012639757841142886 24.991702977377773 +341.74785371597466 248.1210122966632 133.62225322652563 +1030102 +1.9576788467130908 0.027970945719088838 25.55968189226752 +168.9126943096874 174.5978360268346 251.14179333855816 +1030103 +1.967626111006107 0.09226955476091864 1.5279774530000956 +343.6507149897438 15.943105536008861 203.03043240705532 +1030104 +1.9481905063928624 0.0011928078627614026 13.712706682327513 +65.70102045298484 272.5349689212673 248.3878548054673 +1030105 +1.9770133412534745 0.07760502570359305 21.113956523745326 +352.7420976555242 265.31463697173825 90.20915048805514 +1030106 +1.9218154288668319 0.024539743705313734 6.3122997131518 +96.26526838131817 245.40058043017729 203.1458881543752 +1030107 +1.8804043534140047 0.051960914508786316 5.641080857843637 +135.37074231372475 310.57579433118855 92.01889766237802 +1030108 +1.905574620628057 0.0094526649542939 27.795518144975087 +257.055868633706 262.18164024081847 174.33454380356292 +1030109 +1.9610371451706512 0.0642394856575368 11.238178270217231 +155.06613544539204 179.7982091248266 227.65351567825778 +1030110 +1.8442099887042422 0.06774466132861966 7.7755827872476635 +53.89722157548661 231.27772356054047 224.74925586092095 +1030111 +1.918614204096997 0.09311188233380624 22.656238093717935 +8.91787503082778 121.09392635435141 61.064615587866086 +1030112 +1.846083421031052 0.058817512099736215 23.17733139808657 +204.2210767808533 325.0096873835248 80.8104389106319 +1030113 +1.877211316979287 0.09545803603787674 22.75125883423326 +338.55088656488704 260.97440743659024 358.27993618149003 +1030114 +1.847699999560908 0.04422349734135026 24.400107433270254 +175.41066398268853 271.9213024922804 40.10860995923193 +1030115 +1.9693855972875431 0.01087613268562867 18.3688586644925 +288.07225817085316 351.2861239674453 138.17614312177193 +1030116 +1.8499142747340205 0.05785154439430367 18.680424125046468 +170.34777680935684 121.2501486818594 116.63207933400314 +1030117 +1.8484833703421926 0.027361686628703055 12.35340287511839 +77.72361263316543 170.96254103643864 153.80655599509745 +1030118 +1.8844066519945433 0.012253504795668403 28.11812027401767 +5.816089316051807 234.40531004184038 284.0432510837761 +1030119 +1.8321110668963587 0.05520820027904241 4.246848772229706 +195.98962278787414 196.26877791689589 296.01468603344676 +1030120 +1.9298252957677076 0.04051603959563857 11.38835718089586 +14.020062127093711 350.7131191785245 61.91415736777295 +1030121 +1.8296584872076656 0.02265081514507832 19.242848763922552 +170.21169721403766 255.88594273947075 131.4644380844818 +1030122 +1.8699958208259342 0.08248881864547655 6.775274873263883 +131.47697608382563 36.70815817760613 166.48187957273242 +1030123 +1.8427387543178881 0.07305783679736073 5.034586865380743 +311.1450141807796 214.40263993828856 332.6408661254218 +1030124 +1.8567511528637752 0.014548178017731906 8.946207161705699 +27.813385444122652 111.5675183218652 201.4816583454089 +1030125 +1.8846946839402372 0.06619567127431686 13.271555732235106 +271.1533808977903 174.59924852630596 257.0822687250529 diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index f9a7378c0..17eb90d87 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -174,17 +174,17 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date, ds): p6.append(pldata[key].vectors()['vz'][0] * VCONV) p7.append(pldata[key].elements()['a'][0] * DCONV) p8.append(pldata[key].elements()['e'][0]) - p9.append(pldata[key].elements()['incl'][0] * np.pi / 180.0) - p10.append(pldata[key].elements()['Omega'][0] * np.pi / 180.0) - p11.append(pldata[key].elements()['w'][0] * np.pi / 180.0) - p12.append(pldata[key].elements()['M'][0] * np.pi / 180.0) + p9.append(pldata[key].elements()['incl'][0]) + p10.append(pldata[key].elements()['Omega'][0]) + p11.append(pldata[key].elements()['w'][0]) + p12.append(pldata[key].elements()['M'][0]) elif param['OUT_FORM'] == 'EL': p1.append(pldata[key].elements()['a'][0] * DCONV) p2.append(pldata[key].elements()['e'][0]) - p3.append(pldata[key].elements()['incl'][0] * np.pi / 180.0) - p4.append(pldata[key].elements()['Omega'][0] * np.pi / 180.0) - p5.append(pldata[key].elements()['w'][0] * np.pi / 180.0) - p6.append(pldata[key].elements()['M'][0] * np.pi / 180.0) + p3.append(pldata[key].elements()['incl'][0]) + p4.append(pldata[key].elements()['Omega'][0]) + p5.append(pldata[key].elements()['w'][0]) + p6.append(pldata[key].elements()['M'][0]) p7.append(pldata[key].vectors()['x'][0] * DCONV) p8.append(pldata[key].vectors()['y'][0] * DCONV) p9.append(pldata[key].vectors()['z'][0] * DCONV) diff --git a/python/swiftest/swiftest/tool.py b/python/swiftest/swiftest/tool.py index a96610bc2..37e1c6ba9 100644 --- a/python/swiftest/swiftest/tool.py +++ b/python/swiftest/swiftest/tool.py @@ -8,10 +8,10 @@ """ def wrap_angle(angle): - while np.any(angle >= 2 * np.pi): - angle[angle >= 2 * np.pi] -= 2 * np.pi + while np.any(angle >= 360.0 ): + angle[angle >= 360.0] -= 360.0 while np.any(angle < 0.0): - angle[angle < 0.0] += 2 * np.pi + angle[angle < 0.0] += 360.0 return angle def follow_swift(ds, ifol=None, nskp=None): @@ -50,15 +50,9 @@ def follow_swift(ds, ifol=None, nskp=None): intxt = input('Input the print frequency\n') nskp = int(intxt) - dr = 180.0 / np.pi fol['obar'] = fol['capom'] + fol['omega'] fol['obar'] = fol['obar'].fillna(0) fol['obar'] = wrap_angle(fol['obar']) - fol['obar'] = fol['obar'] * dr - fol['inc'] = fol['inc'] * dr - fol['capom'] = fol['capom'] * dr - fol['omega'] = fol['omega'] * dr - fol['capm'] = fol['capm'] * dr fol['peri'] = fol['a'] * (1.0 - fol['e']) fol['apo'] = fol['a'] * (1.0 + fol['e']) diff --git a/src/io/io.f90 b/src/io/io.f90 index 16af0ca10..67d4665ab 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -1101,6 +1101,13 @@ module function io_read_frame_body(self, iu, param) result(ierr) end do end select + + if (param%in_form == EL) then + self%inc(1:n) = self%inc(1:n) * DEG2RAD + self%capom(1:n) = self%capom(1:n) * DEG2RAD + self%omega(1:n) = self%omega(1:n) * DEG2RAD + self%capm(1:n) = self%capm(1:n) * DEG2RAD + end if end associate ierr = 0 @@ -1205,11 +1212,6 @@ module function io_read_frame_system(self, iu, param) result(ierr) goto 667 end if - if (param%in_form == EL) then - call self%pl%el2xv(self%cb) - call self%tp%el2xv(self%cb) - end if - return 667 continue @@ -1482,10 +1484,10 @@ module subroutine io_write_frame_body(self, iu, param) case (EL) write(iu, err = 667, iomsg = errmsg) self%a(1:n) write(iu, err = 667, iomsg = errmsg) self%e(1:n) - write(iu, err = 667, iomsg = errmsg) self%inc(1:n) - write(iu, err = 667, iomsg = errmsg) self%capom(1:n) - write(iu, err = 667, iomsg = errmsg) self%omega(1:n) - write(iu, err = 667, iomsg = errmsg) self%capm(1:n) + write(iu, err = 667, iomsg = errmsg) self%inc(1:n) * RAD2DEG + write(iu, err = 667, iomsg = errmsg) self%capom(1:n) * RAD2DEG + write(iu, err = 667, iomsg = errmsg) self%omega(1:n) * RAD2DEG + write(iu, err = 667, iomsg = errmsg) self%capm(1:n) * RAD2DEG case (XV) write(iu, err = 667, iomsg = errmsg) self%xh(1, 1:n) write(iu, err = 667, iomsg = errmsg) self%xh(2, 1:n) diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index 17c5d9b19..ad886e4f7 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -22,7 +22,8 @@ module swiftest_globals real(DP), parameter :: PI3BY2 = 4.712388980384689857693965074919254326296_DP !! Definition of /(3 \pi / 2\) real(DP), parameter :: TWOPI = 6.283185307179586476925286766559005768394_DP !! Definition of 2 \pi real(DP), parameter :: THIRD = 0.333333333333333333333333333333333333333_DP !! Definition of 1 / 3 - real(DP), parameter :: DEGRAD = 180.0_DP/PI !! Definition of conversion factor from degrees to radians + real(DP), parameter :: DEG2RAD = PI / 180.0_DP !! Definition of conversion factor from degrees to radians + real(DP), parameter :: RAD2DEG = 180.0_DP / PI !! Definition of conversion factor from degrees to radians integer(I4B), parameter :: LOWERCASE_BEGIN = iachar('a') !! ASCII character set parameter for lower to upper conversion - start of lowercase integer(I4B), parameter :: LOWERCASE_END = iachar('z') !! ASCII character set parameter for lower to upper conversion - end of lowercase diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 0a2664feb..049a0acf3 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -141,6 +141,10 @@ module subroutine setup_initialize_system(self, param) call self%set_msys() call self%pl%set_mu(self%cb) call self%tp%set_mu(self%cb) + if (param%in_form == EL) then + call self%pl%el2xv(self%cb) + call self%tp%el2xv(self%cb) + end if call self%pl%eucl_index() if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb) self%pl%lfirst = param%lfirstkick diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 08b5728e8..710936001 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -417,22 +417,24 @@ module subroutine symba_util_rearray_pl(self, system, param) class(symba_parameters), intent(in) :: param !! Current run configuration parameters ! Internals class(symba_pl), allocatable :: tmp !! The discarded body list. - integer(I4B) :: i, j, k + integer(I4B) :: i, j, k, npl logical, dimension(:), allocatable :: lmask class(symba_plplenc), allocatable :: plplenc_old logical :: lencounter - associate(pl => self, pl_adds => system%pl_adds) + associate(pl => self, pl_adds => system%pl_adds, nadd => system%pl_adds%nbody) + npl = pl%nbody ! Deallocate any temporary variables if (allocated(pl%xbeg)) deallocate(pl%xbeg) if (allocated(pl%xend)) deallocate(pl%xend) ! Remove the discards and destroy the list, as the system already tracks pl_discards elsewhere - allocate(lmask, source=pl%ldiscard(:)) - lmask(:) = lmask(:) .or. pl%status(:) == INACTIVE + allocate(lmask(npl)) + lmask(1:npl) = pl%ldiscard(1:npl) .or. pl%status(1:npl) == INACTIVE allocate(tmp, mold=self) call pl%spill(tmp, lspill_list=lmask, ldestructive=.true.) + npl = pl%nbody call tmp%setup(0,param) deallocate(tmp) deallocate(lmask) @@ -442,39 +444,40 @@ module subroutine symba_util_rearray_pl(self, system, param) call plplenc_old%copy(system%plplenc_list) ! Add in any new bodies - if (pl_adds%nbody > 0) then + if (nadd > 0) then ! Append the adds to the main pl object - call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)]) + call pl%append(pl_adds, lsource_mask=[(.true., i=1, nadd)]) + npl = pl%nbody - allocate(lmask(pl%nbody)) - lmask(:) = pl%status(1:pl%nbody) == NEW_PARTICLE + allocate(lmask(npl)) + lmask(1:npl) = pl%status(1:npl) == NEW_PARTICLE - call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, pl%nbody)], lmask)) + call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, npl)], lmask)) deallocate(lmask) end if ! Reset all of the status flags for this body - where(pl%status(1:pl%nbody) /= INACTIVE) - pl%status(1:pl%nbody) = ACTIVE - pl%ldiscard(1:pl%nbody) = .false. - pl%lcollision(1:pl%nbody) = .false. - pl%lmtiny(1:pl%nbody) = pl%Gmass(1:pl%nbody) > param%GMTINY - pl%lmask(1:pl%nbody) = .true. + where(pl%status(1:npl) /= INACTIVE) + pl%status(1:npl) = ACTIVE + pl%ldiscard(1:npl) = .false. + pl%lcollision(1:npl) = .false. + pl%lmtiny(1:npl) = pl%Gmass(1:npl) > param%GMTINY + pl%lmask(1:npl) = .true. elsewhere - pl%ldiscard(1:pl%nbody) = .true. - pl%lmask(1:pl%nbody) = .false. + pl%ldiscard(1:npl) = .true. + pl%lmask(1:npl) = .false. end where - pl%nplm = count(pl%lmtiny(1:pl%nbody) .and. pl%lmask(1:pl%nbody)) + pl%nplm = count(pl%lmtiny(1:npl) .and. pl%lmask(1:npl)) ! Reindex the new list of bodies call pl%sort("mass", ascending=.false.) call pl%eucl_index() ! Reset the kinship trackers - pl%kin(1:pl%nbody)%nchild = 0 - pl%kin(1:pl%nbody)%parent = [(i, i=1, pl%nbody)] + pl%kin(1:npl)%nchild = 0 + pl%kin(1:npl)%parent = [(i, i=1, npl)] ! Re-build the encounter list lencounter = pl%encounter_check(system, param%dt, 0)