Merge pull request #430 from danieljprice/solar_system

Solar system

build/Makefile_setups

@@ -81,7 +81,7 @@ endif
 ifeq ($(SETUP), solarsystem)
 #   orbits of minor planets
     ISOTHERMAL=yes
-    SETUPFILE=utils_mpc.f90 setup_solarsystem.f90
+    SETUPFILE=utils_ephemeris.f90 utils_mpc.f90 setup_solarsystem.f90
     KNOWN_SETUP=yes
     DUST=yes
 endif

scripts/bots.sh

@@ -207,7 +207,7 @@ for edittype in $bots_to_run; do
                'shout' )
                  sed -e 's/SQRT(/sqrt(/g' \
                      -e 's/NINT(/nint(/g' \
-                     -e 's/STOP/stop/g' \
+                     -e 's/ STOP/ stop/g' \
                      -e 's/ATAN/atan/g' \
                      -e 's/ACOS(/acos(/g' \
                      -e 's/ASIN(/asin(/g' \

src/main/utils_datafiles.f90

@@ -18,7 +18,7 @@ module datautils
 ! :Dependencies: None
 !
  implicit none
- public :: find_datafile
+ public :: find_datafile,download_datafile
 
  private
 

src/main/utils_filenames.f90

@@ -21,7 +21,7 @@ module fileutils
  implicit none
  public :: getnextfilename,numfromfile,basename,get_ncolumns,skip_header,get_column_labels,split
  public :: strip_extension,is_digit,files_are_sequential
- public :: ucase,lcase,make_tags_unique,get_nlines,string_delete
+ public :: ucase,lcase,make_tags_unique,get_nlines,string_delete,string_replace,nospaces
 
  private
 
@@ -465,6 +465,31 @@ pure subroutine string_delete(string,skey)
  enddo
 end subroutine string_delete
 
+!---------------------------------------------------------------------------
+!
+! subroutine to replace a matching section of a string with another
+! string, possibly of differing length
+!
+!---------------------------------------------------------------------------
+subroutine string_replace(string,skey,sreplacewith)
+ character(len=*), intent(inout) :: string
+ character(len=*), intent(in)    :: skey,sreplacewith
+ character(len=len(string)) :: remstring
+ integer :: ipos,imax,lensub,i
+
+ ipos = index(trim(string),skey)
+ lensub = len(skey)
+ imax   = len(string)
+ i = 0
+ do while (ipos > 0 .and. i <= imax)
+    i = i + 1  !  only allow as many replacements as characters
+    remstring = string(ipos+lensub:len_trim(string))
+    string = string(1:ipos-1)//sreplacewith//remstring
+    ipos = index(trim(string),skey)
+ enddo
+
+end subroutine string_replace
+
 !---------------------------------------------------------------------------
 !
 ! Split a string into substrings based on a delimiter

src/setup/set_star.f90

@@ -374,7 +374,7 @@ subroutine write_dist(item_in,dist_in,udist)
     write(*,'(1x,2(a,es12.5),a)') item_in, dist_in*udist,' cm     = ',dist_in,' km'
  else
     write(*,'(1x,2(a,es12.5),a)') item_in, dist_in*udist,' cm     = ',dist_in*udist/solarr,' R_sun'
-  endif
+ endif
 
 end subroutine write_dist
 !-----------------------------------------------------------------------

src/setup/setup_solarsystem.f90

@@ -10,14 +10,14 @@ module setup
 !
 ! :References: https://minorplanetcenter.net/data
 !
-! :Owner: Not Committed Yet
+! :Owner: Daniel Price
 !
 ! :Runtime parameters:
 !   - dumpsperorbit : *number of dumps per orbit*
 !   - norbits       : *number of orbits*
 !
-! :Dependencies: infile_utils, io, mpc, part, physcon, setbinary, timestep,
-!   units
+! :Dependencies: centreofmass, datautils, ephemeris, infile_utils, io, mpc,
+!   part, physcon, setbinary, timestep, units
 !
  implicit none
  public :: setpart
@@ -35,15 +35,16 @@ module setup
 !+
 !----------------------------------------------------------------
 subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,time,fileprefix)
- use part,      only:nptmass,xyzmh_ptmass,vxyz_ptmass,idust,ndustlarge,set_particle_type,&
-                     grainsize,graindens,ndusttypes
- use setbinary, only:set_binary
- use units,     only:set_units,umass,udist,unit_density
- use physcon,   only:solarm,au,pi,km
- use io,        only:master,fatal
- use timestep,  only:tmax,dtmax
- use mpc,       only:read_mpc,mpc_entry
- use datautils, only:find_datafile
+ use part,         only:nptmass,xyzmh_ptmass,vxyz_ptmass,idust,set_particle_type,&
+                        grainsize,graindens,ndustlarge,ndusttypes
+ use setbinary,    only:set_binary
+ use units,        only:set_units,umass,udist,unit_density
+ use physcon,      only:solarm,au,pi,km
+ use io,           only:master,fatal
+ use timestep,     only:tmax,dtmax
+ use mpc,          only:read_mpc,mpc_entry
+ use datautils,    only:find_datafile
+ use centreofmass, only:reset_centreofmass
  integer,           intent(in)    :: id
  integer,           intent(inout) :: npart
  integer,           intent(out)   :: npartoftype(:)
@@ -97,16 +98,14 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  vxyzu(:,:) = 0.
  nptmass = 0
 
- semia  = 5.2  !  Jupiter
+ semia  = 1.  !  Earth
  mtot   = solarm/umass
  hpart  = 100.*au/udist
 
  period = 2.*pi*sqrt(semia**3/mtot)
  tmax   = norbits*period
  dtmax  = period/dumpsperorbit
-!
-! read the orbital data from the data file
-!
+
  filename = find_datafile('Distant.txt',url='https://www.minorplanetcenter.net/iau/MPCORB/')
  call read_mpc(filename,nbodies,dat=dat)
  print "(a,i0,a)",' read orbital data for ',nbodies,' minor planets'
@@ -116,7 +115,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
     ! for each solar system object get the xyz positions from the orbital parameters
     !
     !print*,i,'aeiOwM=',dat(i)%a,dat(i)%ecc,dat(i)%inc,dat(i)%O,dat(i)%w,dat(i)%M
-    call set_binary(mtot,epsilon(0.),dat(i)%a,dat(i)%ecc,1.,1.e-15,&
+    call set_binary(mtot,epsilon(0.),dat(i)%a,dat(i)%ecc,0.02,1.e-15,&
                     xyzmh_ptmass,vxyz_ptmass,nptmass,ierr,incl=dat(i)%inc,&
                     arg_peri=dat(i)%w,posang_ascnode=dat(i)%O,&
                     mean_anomaly=dat(i)%M,verbose=.false.)
@@ -144,12 +143,75 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  grainsize(1:ndustlarge) = km/udist         ! assume km-sized bodies
  graindens(1:ndustlarge) = 2./unit_density  ! 2 g/cm^3
 
+ !
+ ! add the planets
+ !
+ call set_solarsystem_planets(nptmass,xyzmh_ptmass,vxyz_ptmass)
+
+ call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)
  hfact = 1.2
 
  if (ierr /= 0) call fatal('setup','ERROR during setup')
 
 end subroutine setpart
 
+!----------------------------------------------------------------
+!+
+!  setup the solar system planets by querying their ephemeris
+!  from the JPL server
+!+
+!----------------------------------------------------------------
+subroutine set_solarsystem_planets(nptmass,xyzmh_ptmass,vxyz_ptmass)
+ use ephemeris, only:get_ephemeris,nelem
+ use units,     only:umass,udist
+ use physcon,   only:gg,km,solarm,earthm,au
+ use setbinary, only:set_binary
+ integer, intent(inout) :: nptmass
+ real,    intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:)
+ integer,          parameter :: nplanets = 9
+ character(len=*), parameter :: planet_name(nplanets) = &
+     (/'mercury', &
+       'venus  ', &
+       'earth  ', &
+       'mars   ', &
+       'jupiter', &
+       'saturn ', &
+       'uranus ', &
+       'neptune', &
+       'pluto  '/)  ! for nostalgia's sake
+ real    :: elems(nelem),xyz_tmp(size(xyzmh_ptmass(:,1)),2),vxyz_tmp(3,2),gm_cgs
+ real    :: msun,mplanet,a,e,inc,O,w,f
+ integer :: i,ierr,ntmp
+
+ msun = solarm/umass
+ do i=1,nplanets
+    elems = get_ephemeris(planet_name(i),ierr)
+    if (ierr /= 0) then
+       print "(a)",' ERROR: could not read ephemeris data for '//planet_name(i)
+       cycle  ! skip if error reading ephemeris file
+    endif
+    gm_cgs  = elems(1)*km**3
+    mplanet = (gm_cgs/gg)/umass
+    a   = elems(2)*km/udist
+    e   = elems(3)
+    inc = elems(4)
+    O   = elems(5)
+    w   = elems(6)
+    f   = elems(7)
+    print*,' mplanet/mearth = ',mplanet*umass/earthm,' a = ',a*udist/au,' au'
+    ntmp = 0
+    call set_binary(msun,mplanet,a,e,0.01,0.01,&
+                    xyz_tmp,vxyz_tmp,ntmp,ierr,incl=inc,&
+                    arg_peri=w,posang_ascnode=O,f=f,verbose=.false.)
+    nptmass = nptmass + 1
+    xyzmh_ptmass(:,nptmass) = xyz_tmp(:,2)
+    vxyz_ptmass(:,nptmass)  = vxyz_tmp(:,2)
+ enddo
+
+ print*,' nptmass = ',nptmass
+
+end subroutine set_solarsystem_planets
+
 !----------------------------------------------------------------
 !+
 !  write setup parameters to file

src/utils/moddump_rad_to_LTE.f90

@@ -5,30 +5,21 @@
 ! http://phantomsph.bitbucket.io/                                          !
 !--------------------------------------------------------------------------!
 module moddump
- !
- ! Convert a radiative dump into a non-radiative dump, assuming LTE (i.e., ieos=12)
- ! INSTRUCTIONS: 1) Set gmw (mean molecular weight) in the code below to the
- !                  same value assumed in the radiative dump
- !               1) Make this moddump with a radiative setup (e.g. SETUP=radstar),
- !                  otherwise the rad array is not read.
- !               2) Run this moddump
- !               3) Write/modify the infile to remove radiation-related options,
- !                  make sure ieos=12 and mu is correct
- !               4) Recompile Phantom with the desired non-radiative setup (e.g.
- !                  SETUP=star).
- !
- ! :References: None
- !
- ! :Owner: Mike Lau
- !
- ! :Runtime parameters: None
- !
- ! :Dependencies: 
- !
+!
+! moddump
+!
+! :References: None
+!
+! :Owner: Mike Lau
+!
+! :Runtime parameters: None
+!
+! :Dependencies: dim, eos, io, part
+!
  implicit none
-    
- contains
-    
+
+contains
+
 subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  use dim, only:do_radiation
  use io,  only:fatal
@@ -39,7 +30,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  real,    intent(inout) :: massoftype(:)
  real,    intent(inout) :: xyzh(:,:),vxyzu(:,:)
  integer                :: i
- 
+
  if (.not. do_radiation) call fatal("moddump_rad_to_LTE","Not compiled with radiation")
  do i=1,npart
     if (isnan(rad(iradxi,i))) call fatal("moddump_rad_to_LTE","rad array contains NaNs")
@@ -51,7 +42,6 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  print*,'mu has been changed to',gmw  ! mu should not change from what was assumed with radiation
 
 end subroutine modify_dump
-    
+
 end module moddump
-
-
+