Merge pull request #424 from danieljprice/solar_system

Solar system

build/Makefile_setups

@@ -78,6 +78,14 @@ ifeq ($(SETUP), asteroidwind)
     KNOWN_SETUP=yes
 endif
 
+ifeq ($(SETUP), solarsystem)
+#   orbits of minor planets
+    ISOTHERMAL=yes
+    SETUPFILE=utils_mpc.f90 setup_solarsystem.f90
+    KNOWN_SETUP=yes
+    DUST=yes
+endif
+
 ifeq ($(SETUP), galdisc)
 #   galactic disc simulations
     IND_TIMESTEPS=yes

src/main/utils_binary.f90

@@ -7,6 +7,8 @@
 module binaryutils
 !
 ! Utility routines for binary setup, also used in asteroid injection
+! Main thing is routines to compute eccentric anomaly from mean anomaly
+! by solving Kepler's equation
 !
 ! :References: None
 !
@@ -19,25 +21,40 @@ module binaryutils
  implicit none
  real, parameter :: pi = 4.*atan(1.)
 
+ public :: get_E,get_E_from_mean_anomaly
+ public :: get_orbit_bits
+
 contains
 
 !---------------------------------------------------------------
 !+
-! Get eccentric anomaly (this function uses bisection
-! to guarantee convergence, which is not guaranteed for
-! small M or E)
+! Get eccentric anomaly given time since pericentre
 !+
 !---------------------------------------------------------------
 subroutine get_E(period,ecc,deltat,E)
  real, intent(in)  :: period,ecc,deltat
  real, intent(out) :: E
- real :: mu,M_ref,M_guess
- real :: E_left,E_right,E_guess
- real, parameter :: tol = 1.e-10
+ real :: mu,M_ref
 
  mu = 2.*pi/period
  M_ref = mu*deltat ! mean anomaly
 
+ E = get_E_from_mean_anomaly(M_ref,ecc)
+
+end subroutine get_E
+
+!---------------------------------------------------------------
+!+
+! Get eccentric anomaly from mean anomaly (this function uses
+! bisection to guarantee convergence, which is not guaranteed for
+! small M or E)
+!+
+!---------------------------------------------------------------
+real function get_E_from_mean_anomaly(M_ref,ecc) result(E)
+ real, intent(in) :: M_ref,ecc
+ real :: E_left,E_right,E_guess,M_guess
+ real, parameter :: tol = 1.e-10
+
  ! first guess
  E_left = 0.
  E_right = 2.*pi
@@ -56,7 +73,7 @@ subroutine get_E(period,ecc,deltat,E)
 
  E = E_guess
 
-end subroutine get_E
+end function get_E_from_mean_anomaly
 
 !-----------------------------------------------------------------------
 !+

src/main/utils_datafiles.f90

@@ -52,20 +52,25 @@ function find_datafile(filename,dir,env_var,url,verbose) result(filepath)
  inquire(file=trim(filename),exist=iexist)
  if (iexist) then
     filepath = trim(filename)
+    if (isverbose) print "(a)",' reading '//trim(filepath)
  else
     ierr = 0
     mydir  = ' '
     if (present(env_var)) then
        my_env_var = env_var
+       call get_environment_variable(my_env_var,env_dir)
+    elseif (present(dir)) then
+       env_dir = dir
     else
-       my_env_var = 'DATA_DIR'
+       env_dir = './'
     endif
-    call get_environment_variable(my_env_var,env_dir)
     if (len_trim(env_dir) > 0) then
        mydir = trim(env_dir)
        if (present(dir)) mydir = trim(mydir)//'/'//trim(dir)//'/'
-       if (isverbose) print "(a)",' Reading '//trim(filename)//' in '//trim(mydir)//&
+       if (isverbose .and. present(env_var)) then
+          print "(a)",' Reading '//trim(filename)//' in '//trim(mydir)//&
                                   ' (from '//trim(my_env_var)//' setting)'
+       endif
        filepath = trim(mydir)//trim(filename)
        inquire(file=trim(filepath),exist=iexist)
        if (.not.iexist) then

src/setup/set_binary.f90

@@ -55,19 +55,19 @@ module setbinary
 subroutine set_binary(m1,m2,semimajoraxis,eccentricity, &
                       accretion_radius1,accretion_radius2, &
                       xyzmh_ptmass,vxyz_ptmass,nptmass,ierr,omega_corotate,&
-                      posang_ascnode,arg_peri,incl,f,verbose)
- use binaryutils, only:get_E
+                      posang_ascnode,arg_peri,incl,f,mean_anomaly,verbose)
+ use binaryutils, only:get_E,get_E_from_mean_anomaly
  real,    intent(in)    :: m1,m2
  real,    intent(in)    :: semimajoraxis,eccentricity
  real,    intent(in)    :: accretion_radius1,accretion_radius2
  real,    intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:)
  integer, intent(inout) :: nptmass
  integer, intent(out)   :: ierr
- real,    intent(in),  optional :: posang_ascnode,arg_peri,incl,f
+ real,    intent(in),  optional :: posang_ascnode,arg_peri,incl,f,mean_anomaly
  real,    intent(out), optional :: omega_corotate
  logical, intent(in),  optional :: verbose
  integer :: i1,i2,i
- real    :: mtot,dx(3),dv(3),Rochelobe1,Rochelobe2,period
+ real    :: mtot,dx(3),dv(3),Rochelobe1,Rochelobe2,period,bigM
  real    :: x1(3),x2(3),v1(3),v2(3),omega0,cosi,sini,xangle,reducedmass,angmbin
  real    :: a,E,E_dot,P(3),Q(3),omega,big_omega,inc,ecc,tperi,term1,term2,theta
  logical :: do_verbose
@@ -146,6 +146,10 @@ subroutine set_binary(m1,m2,semimajoraxis,eccentricity, &
        ! (https://en.wikipedia.org/wiki/Eccentric_anomaly#From_the_true_anomaly)
        theta = f*pi/180.
        E = atan2(sqrt(1. - ecc**2)*sin(theta),(ecc + cos(theta)))
+    elseif (present(mean_anomaly)) then
+       ! get eccentric anomaly from mean anomaly by solving Kepler equation
+       bigM = mean_anomaly*pi/180.
+       E = get_E_from_mean_anomaly(bigM,ecc)
     else
        ! set binary at apastron
        tperi = 0.5*period ! time since periastron: half period = apastron
@@ -173,8 +177,10 @@ subroutine set_binary(m1,m2,semimajoraxis,eccentricity, &
              'inclination     (i, deg):',incl, &
              'angle asc. node (O, deg):',posang_ascnode, &
              'arg. pericentre (w, deg):',arg_peri
-       if (present(f)) print "(2x,a,g12.4)", &
+       if (present(f)) print "(2x,a,1pg14.6)", &
              'true anomaly    (f, deg):',f
+       if (present(mean_anomaly)) print "(2x,a,1pg14.6)", &
+             'mean anomaly    (M, deg):',mean_anomaly
     endif
 
     ! Rotating everything

src/setup/setup_sedov.f90

@@ -18,7 +18,7 @@ module setup
 !   - shuffle_parts : *relax particles by shuffling*
 !
 ! :Dependencies: boundary, infile_utils, io, kernel, mpidomain, mpiutils,
-!   options, part, physcon, setup_params, timestep, unifdis,
+!   options, part, physcon, prompting, setup_params, timestep, unifdis,
 !   utils_shuffleparticles
 !
  implicit none

src/setup/setup_solarsystem.f90

@@ -0,0 +1,202 @@
+!--------------------------------------------------------------------------!
+! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. !
+! Copyright (c) 2007-2023 The Authors (see AUTHORS)                        !
+! See LICENCE file for usage and distribution conditions                   !
+! http://phantomsph.bitbucket.io/                                          !
+!--------------------------------------------------------------------------!
+module setup
+!
+! Setup asteroid orbits using data from the IAU Minor Planet Center
+!
+! :References: https://minorplanetcenter.net/data
+!
+! :Owner: Not Committed Yet
+!
+! :Runtime parameters:
+!   - dumpsperorbit : *number of dumps per orbit*
+!   - norbits       : *number of orbits*
+!
+! :Dependencies: infile_utils, io, mpc, part, physcon, setbinary, timestep,
+!   units
+!
+ implicit none
+ public :: setpart
+
+ real :: norbits
+ integer :: dumpsperorbit
+
+ private
+
+contains
+
+!----------------------------------------------------------------
+!+
+!  setup for solar system orbits
+!+
+!----------------------------------------------------------------
+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
+ integer,           intent(in)    :: id
+ integer,           intent(inout) :: npart
+ integer,           intent(out)   :: npartoftype(:)
+ real,              intent(out)   :: xyzh(:,:)
+ real,              intent(inout) :: massoftype(:)
+ real,              intent(inout) :: polyk,gamma,hfact
+ real,              intent(inout) :: time
+ character(len=20), intent(in)    :: fileprefix
+ real,              intent(out)   :: vxyzu(:,:)
+ character(len=120) :: filename
+ integer :: ierr,nbodies,i
+ logical :: iexist
+ real    :: period,semia,mtot,hpart
+ integer, parameter :: max_bodies = 2000000
+ type(mpc_entry), allocatable :: dat(:)
+!
+! default runtime parameters
+!
+ norbits       = 1000.
+ dumpsperorbit = 1
+!
+! read runtime parameters from setup file
+!
+ if (id==master) print "(/,65('-'),1(/,a),/,65('-'),/)",' solar system'
+ filename = trim(fileprefix)//'.setup'
+ inquire(file=filename,exist=iexist)
+ if (iexist) call read_setupfile(filename,ierr)
+ if (.not. iexist .or. ierr /= 0) then
+    if (id==master) then
+       call write_setupfile(filename)
+       print*,' Edit '//trim(filename)//' and rerun phantomsetup'
+    endif
+    stop
+ endif
+!
+! set units
+!
+ call set_units(mass=solarm,dist=au,G=1.d0)
+!
+! general parameters
+!
+ time  = 0.
+ polyk = 0.
+ gamma = 1.
+!
+!--space available for injected gas particles
+!
+ npart = 0
+ npartoftype(:) = 0
+ xyzh(:,:)  = 0.
+ vxyzu(:,:) = 0.
+ nptmass = 0
+
+ semia  = 5.2  !  Jupiter
+ 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'
+
+ do i=1,nbodies
+    !
+    ! 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,&
+                    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.)
+    !
+    ! now delete the point masses but set a dust particle as the secondary
+    !
+    nptmass = 0
+    xyzh(1:3,i)  = xyzmh_ptmass(1:3,2)
+    xyzh(4,i)    = hpart  ! give a random length scale as the smoothing length
+    vxyzu(1:3,i) = vxyz_ptmass(1:3,2)
+    call set_particle_type(i,idust)
+ enddo
+ !
+ ! restore the Sun
+ !
+ nptmass = 1
+ !
+ ! set mass of all the minor bodies equal
+ !
+ npart = nbodies
+ ndustlarge = 1
+ ndusttypes = 1
+ npartoftype(idust) = nbodies
+ massoftype(idust) = 1.e-20
+ grainsize(1:ndustlarge) = km/udist         ! assume km-sized bodies
+ graindens(1:ndustlarge) = 2./unit_density  ! 2 g/cm^3
+
+ hfact = 1.2
+
+ if (ierr /= 0) call fatal('setup','ERROR during setup')
+
+end subroutine setpart
+
+!----------------------------------------------------------------
+!+
+!  write setup parameters to file
+!+
+!----------------------------------------------------------------
+subroutine write_setupfile(filename)
+ use infile_utils, only:write_inopt
+ character(len=*), intent(in) :: filename
+ integer, parameter :: iunit = 20
+
+ print "(a)",' writing setup options file '//trim(filename)
+ open(unit=iunit,file=filename,status='replace',form='formatted')
+
+ write(iunit,"(a)") '# input file for solar system setup routines'
+ call write_inopt(norbits,'norbits','number of orbits',iunit)
+ call write_inopt(dumpsperorbit,'dumpsperorbit','number of dumps per orbit',iunit)
+
+ close(iunit)
+
+end subroutine write_setupfile
+
+!----------------------------------------------------------------
+!+
+!  read setup parameters from file
+!+
+!----------------------------------------------------------------
+subroutine read_setupfile(filename,ierr)
+ use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db
+ use io,           only:error
+ character(len=*), intent(in)  :: filename
+ integer,          intent(out) :: ierr
+ integer, parameter :: iunit = 21
+ integer :: nerr
+ type(inopts), allocatable :: db(:)
+
+ nerr = 0
+ ierr = 0
+ call open_db_from_file(db,filename,iunit,ierr)
+ call read_inopt(norbits,      'norbits',      db,min=0.,errcount=nerr)
+ call read_inopt(dumpsperorbit,'dumpsperorbit',db,min=0 ,errcount=nerr)
+ call close_db(db)
+
+ if (nerr > 0) then
+    print "(1x,i2,a)",nerr,' error(s) during read of setup file: re-writing...'
+    ierr = nerr
+ endif
+
+end subroutine read_setupfile
+
+end module setup