diff --git a/build/Makefile_setups b/build/Makefile_setups
index 28352dda4..1c4b919a5 100644
--- a/build/Makefile_setups
+++ b/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
diff --git a/scripts/bots.sh b/scripts/bots.sh
index 131288648..46fe18b8f 100755
--- a/scripts/bots.sh
+++ b/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' \
diff --git a/src/main/utils_datafiles.f90 b/src/main/utils_datafiles.f90
index b138d0b16..94cb3c3f6 100644
--- a/src/main/utils_datafiles.f90
+++ b/src/main/utils_datafiles.f90
@@ -18,7 +18,7 @@ module datautils
 ! :Dependencies: None
 !
  implicit none
- public :: find_datafile
+ public :: find_datafile,download_datafile
 
  private
 
diff --git a/src/main/utils_filenames.f90 b/src/main/utils_filenames.f90
index f0acb635c..c2a6308eb 100644
--- a/src/main/utils_filenames.f90
+++ b/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
diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90
index 83cb9aea8..7a617d3c3 100644
--- a/src/setup/set_star.f90
+++ b/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
 !-----------------------------------------------------------------------
diff --git a/src/setup/setup_solarsystem.f90 b/src/setup/setup_solarsystem.f90
index 72015427a..333519a44 100644
--- a/src/setup/setup_solarsystem.f90
+++ b/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
diff --git a/src/utils/moddump_rad_to_LTE.f90 b/src/utils/moddump_rad_to_LTE.f90
index 865ddf1f9..c3c309830 100644
--- a/src/utils/moddump_rad_to_LTE.f90
+++ b/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
-    
-    
\ No newline at end of file
+
diff --git a/src/utils/utils_ephemeris.f90 b/src/utils/utils_ephemeris.f90
new file mode 100644
index 000000000..66664e39c
--- /dev/null
+++ b/src/utils/utils_ephemeris.f90
@@ -0,0 +1,215 @@
+!--------------------------------------------------------------------------!
+! 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 ephemeris
+!
+! Download and read data files from the JPL horizons ephemeris service
+!
+! :References:
+! https://ssd-api.jpl.nasa.gov/doc/horizons.html
+!
+! :Owner: Daniel Price
+!
+! :Runtime parameters: None
+!
+! :Dependencies: datautils, fileutils
+!
+ implicit none
+ integer, parameter :: nelem = 7
+
+ public :: get_ephemeris,nelem
+
+ private
+
+contains
+
+!-----------------------------------------------------------------------
+!+
+!  read data files from the JPL Horizons ephemeris server
+!+
+!-----------------------------------------------------------------------
+function get_ephemeris(object,ierr) result(elems)
+ use datautils, only:download_datafile
+ character(len=*), intent(in)  :: object  ! name of the solar system object
+ integer,          intent(out) :: ierr
+ real :: elems(nelem)
+ character(len=512) :: url
+ character(len=30)  :: localfile,filepath
+ logical :: iexist
+ integer :: ierr2
+
+ ierr = 0
+ localfile = trim(object)//'.txt'
+ inquire(file=localfile,exist=iexist)
+ if (.not.iexist) then
+    call construct_horizons_api_url(trim(object),url,ierr)
+    if (ierr /= 0) then
+       print*,' ERROR: could not query horizons database for object='//trim(object)
+       return
+    endif
+    call download_datafile(url=url,dir=localfile,filename='',filepath=filepath,ierr=ierr)
+ else
+    filepath = localfile
+ endif
+
+ call read_ephemeris_file(filepath,elems,ierr2)
+ if (ierr2 /= 0) ierr = ierr2
+
+end function get_ephemeris
+
+!-----------------------------------------------------------------------
+!+
+!  construct API call to the JPL Horizons ephemeris server
+!+
+!-----------------------------------------------------------------------
+subroutine construct_horizons_api_url(object,url,ierr)
+ character(len=*), intent(in)  :: object  ! name of the solar system object
+ character(len=*), intent(out) :: url     ! url for query
+ integer,          intent(out) :: ierr
+ character(len=3)  :: cmd
+ character(len=10) :: start_epoch,end_epoch
+ integer(kind=8)   :: values(8),year,month,day
+
+ ierr = 0
+ select case(trim(adjustl(object)))
+ case('pluto')
+    cmd = '999' ! pluto barycentre
+ case('neptune')
+    cmd = '899' ! neptune barycentre
+ case('uranus')
+    cmd = '799' ! uranus barycentre
+ case('saturn')
+    cmd = '699' ! saturn barycentre
+ case('jupiter')
+    cmd = '599' ! jupiter barycentre
+ case('mars')
+    cmd = '499' ! mars barycentre
+ case('earth')
+    cmd = '399' ! earth-moon barycentre
+ case('venus')
+    cmd = '299' ! venus barycentre
+ case('mercury')
+    cmd = '199' ! mercury barycentre
+ case default
+    ierr = 1
+ end select
+
+ !if (present(epoch)) then
+!    start_epoch = epoch
+ !else
+ call date_and_time(values=values)
+ year = values(1); month = values(2); day = values(3)
+ write(start_epoch,"(i4.4,'-',i2.2,'-',i2.2)") year,month,day
+ write(end_epoch,"(i4.4,'-',i2.2,'-',i2.2)") year,month,day+1
+ !endif
+
+ url = "'https://ssd.jpl.nasa.gov/api/horizons.api?format=text&COMMAND='"//trim(cmd)// &
+       "'&OBJ_DATA='YES'&MAKE_EPHEM='YES'&EPHEM_TYPE='ELEMENTS'&CENTER='500@10'&START_TIME='"&
+        //trim(start_epoch)//"'&STOP_TIME='"//trim(end_epoch)// &
+       "'&STEP_SIZE='1DAYS'&REF_SYSTEM='ICRF'&REF_PLANE='ECLIPTIC''"
+
+end subroutine construct_horizons_api_url
+
+!-----------------------------------------------------------------------
+!+
+!  read required information from the ephemeris data file
+!+
+!-----------------------------------------------------------------------
+subroutine read_ephemeris_file(file,elems,ierr)
+ character(len=*), intent(in)  :: file
+ real,             intent(out) :: elems(nelem)
+ integer, intent(out) :: ierr
+ integer :: iu,j
+ character(len=80)  :: line
+ character(len=*), parameter :: tag(nelem) = &
+    (/'GM km^3/s^2',&
+      'A          ',&
+      'EC         ',&
+      'IN         ',&
+      'OM         ',&
+      'W          ',&
+      'TA         '/)
+ logical :: got_elem(nelem)
+
+ ! give default parameters
+ got_elem(:) = .false.
+ elems(:) = 0.
+
+ open(newunit=iu,file=file,status='old',iostat=ierr)
+ if (ierr /= 0) then
+    print "(a)",' ERROR opening '//trim(file)
+    return
+ endif
+ print "(/,a)",' > reading ephemeris from '//trim(file)
+ do while(ierr == 0)
+    read(iu,"(a)",iostat=ierr) line
+    do j=1,nelem
+       if (.not.got_elem(j)) then
+          call read_value(line,tag(j),elems(j),got_elem(j))
+       endif
+    enddo
+ enddo
+ ierr = 0
+ print "(a,/)"
+ do j=1,nelem
+    if (.not.got_elem(j)) then
+       print*,' ERROR: could not find '//trim(tag(j))//' in '//trim(file)
+       ierr = ierr + 1
+    endif
+ enddo
+ close(iu)
+
+end subroutine read_ephemeris_file
+
+!-----------------------------------------------------------------------
+!+
+!  utility routine to read a var = blah string from the ephemeris file
+!+
+!-----------------------------------------------------------------------
+subroutine read_value(line,tag,val,got_val)
+ use fileutils, only:string_delete,string_replace
+ character(len=*), intent(in)    :: line, tag
+ real,             intent(out)   :: val
+ logical,          intent(inout) :: got_val
+ character(len=len(line)) :: string
+ integer :: ieq,j,ierr
+
+ val = 0.
+ if (.not.got_val) then
+    ! make a copy of the line string
+    string = line
+
+    call string_delete(string,'(planet)') ! pluto has GM(planet) km^3/s^2
+
+    ! delete double spaces
+    call string_replace(string,'  ',' ')
+    call string_delete(string,'(')
+    call string_delete(string,')')
+    call string_delete(string,',')
+
+    ! add space at start of string if necessary
+    if (string(1:1) /= ' ') string = ' '//trim(string)
+
+    ! search for ' TAG = value' in the line, including spaces
+    ! to avoid confusion of A = val with TA= val
+    j = max(index(string,' '//trim(tag)//' ='),index(string,' '//trim(tag)//'='))
+
+    ! if we get a match read the value from what is after the equals sign
+    if (j > 0) then
+       string = string(j:)  ! start
+       ieq = index(string,'=')
+       read(string(ieq+1:),*,iostat=ierr) val
+       if (ierr == 0) then
+          ! mark this quantity as having already been read
+          write(*,"(1x,a,g0)",advance='no') trim(tag)//' = ',val
+          got_val = .true.
+       endif
+    endif
+ endif
+
+end subroutine read_value
+
+end module ephemeris
diff --git a/src/utils/utils_mpc.f90 b/src/utils/utils_mpc.f90
index 9a14e5f30..269a5c296 100644
--- a/src/utils/utils_mpc.f90
+++ b/src/utils/utils_mpc.f90
@@ -28,7 +28,7 @@ module mpc
 ! :References:
 ! https://minorplanetcenter.net/Extended_Files/Extended_MPCORB_Data_Format_Manual.pdf
 !
-! :Owner: Not Committed Yet
+! :Owner: Daniel Price
 !
 ! :Runtime parameters: None
 !