diff --git a/Makefile b/Makefile index c3cce41c2..5f3459f1a 100644 --- a/Makefile +++ b/Makefile @@ -170,12 +170,21 @@ fast: ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make fastdir - cd $(SWIFTEST_HOME)/src/util; \ rm -f Makefile.Defines Makefile; \ ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make fastdir + cd $(SWIFTEST_HOME)/src/orbel; \ + rm -f Makefile.Defines Makefile; \ + ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ + ln -s $(SWIFTEST_HOME)/Makefile .; \ + make fastdir + cd $(SWIFTEST_HOME)/src/drift; \ + rm -f Makefile.Defines Makefile; \ + ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ + ln -s $(SWIFTEST_HOME)/Makefile .; \ + make fastdir cd $(SWIFTEST_HOME)/src/orbel; \ rm -f Makefile.Defines Makefile; \ @@ -207,6 +216,12 @@ fast: $(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \ rm -f *.o *.smod + cd $(SWIFTEST_HOME)/src/helio; \ + $(FORTRAN) $(FFASTFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c helio_drift.f90; \ + $(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \ + $(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \ + rm -f *.o *.smod + libdir: $(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c *.f90; \ $(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \ diff --git a/Makefile.Defines b/Makefile.Defines index 16ce3afc3..ea2772b04 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -57,7 +57,7 @@ SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdta PAR = -qopenmp -parallel HEAPARR = -heap-arrays 4194304 OPTREPORT = -qopt-report=5 -IPRODUCTION = -no-wrap-margin -O3 -ipo -qopt-prefetch=0 -sox $(PAR) $(SIMDVEC) #$(HEAPARR) +IPRODUCTION = -no-wrap-margin -O3 -qopt-prefetch=0 -sox $(PAR) $(SIMDVEC) #$(HEAPARR) #gfortran flags GDEBUG = -g -Og -fbacktrace -fbounds-check -ffree-line-length-none diff --git a/examples/symba_mars_disk/param.in b/examples/symba_mars_disk/param.in index 023f31647..fd0cc9134 100644 --- a/examples/symba_mars_disk/param.in +++ b/examples/symba_mars_disk/param.in @@ -1,6 +1,6 @@ !Parameter file for the SyMBA-RINGMOONS test T0 0.0 -TSTOP 1200.0 +TSTOP 12000.0 DT 600.0 CB_IN cb.in PL_IN mars.in diff --git a/src/orbel/orbel.f90 b/src/orbel/orbel.f90 index a815b92c7..5de45e296 100644 --- a/src/orbel/orbel.f90 +++ b/src/orbel/orbel.f90 @@ -254,40 +254,24 @@ real(DP) pure function orbel_flon(e,icapn) biga = (-0.5_DP * b + sq)**(1.0_DP / 3.0_DP) bigb = -(+0.5_DP * b + sq)**(1.0_DP / 3.0_DP) x = biga + bigb - ! write(6,*) 'cubic = ',x**3 +a*x +b orbel_flon = x ! If capn is VSMALL (or zero) no need to go further than cubic even for ! e =1. - if( capn < VSMALL) go to 100 - - do i = 1,IMAX - x2 = x * x - f = a0 + x * (a1 + x2 * (a3 + x2 * (a5 + x2 * (a7 + x2 * (a9 + x2 * (a11 + x2)))))) - fp = b1 + x2 * (b3 + x2 * (b5 + x2 * (b7 + x2 * (b9 + x2 * (b11 + 13 * x2))))) - dx = -f / fp - ! write(6,*) 'i,dx,x,f : ' - ! write(6,432) i,dx,x,f - 432 format(1x,i3,3(2x,1p1e22.15)) - orbel_flon = x + dx - ! if we have converged here there's no point in going on - if(abs(dx) <= VSMALL) go to 100 + if( capn >= VSMALL) then + do i = 1,IMAX + x2 = x**2 + f = a0 + x * (a1 + x2 * (a3 + x2 * (a5 + x2 * (a7 + x2 * (a9 + x2 * (a11 + x2)))))) + fp = b1 + x2 * (b3 + x2 * (b5 + x2 * (b7 + x2 * (b9 + x2 * (b11 + 13 * x2))))) + dx = -f / fp + orbel_flon = x + dx + ! if we have converged here there's no point in going on + if(abs(dx) <= VSMALL) exit x = orbel_flon - end do - - ! abnormal return here - we've gone thru the loop - ! imax times without convergence - if(iflag == 1) then - orbel_flon = -orbel_flon - capn = -capn + end do end if - !write(*,*) 'flon : returning without complete convergence' - diff = e * sinh(orbel_flon) - orbel_flon - capn - !write(*,*) 'n, f, ecc*sinh(f) - f - n : ' - !write(*,*) capn,orbel_flon,diff - return ! normal return here, but check if capn was originally negative - 100 if(iflag == 1) then + if(iflag == 1) then orbel_flon = -orbel_flon capn = -capn end if