Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

(solarsystem) use km as length unit for solar system setup #629

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

(solarsystem) use km as length unit for solar system setup #629

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/main/utils_infiles.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -696,8 +696,9 @@ subroutine read_inopt_from_line(line,name,valstring,ierr,comment)
!
!--for time strings, assume they are of the form hh:mm:ss
! convert to a number of seconds as a real
! but be careful to ignore datetime strings like 2020-10-04 12:00:00
!
if (index(valstring,':') /= 0) then
if (index(valstring,':') /= 0 .and. index(valstring,'-')==0) then
nsec = 0
read(valstring,"(i3.3,1x,i2.2,1x,i2.2)",iostat=ierr) nhr,nmin,nsec
if (ierr/=0) then
Expand Down
29 changes: 19 additions & 10 deletions src/setup/set_solarsystem.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,10 +105,11 @@ end subroutine set_minor_planets
! from the JPL server
!+
!----------------------------------------------------------------
subroutine add_sun_and_planets(nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,epoch)
subroutine add_sun_and_planets(nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,nerr,epoch)
integer, intent(inout) :: nptmass
real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:)
real, intent(in) :: mtot
integer, intent(out) :: nerr
character(len=*), intent(in), optional :: epoch
integer, parameter :: nbodies = 10
character(len=*), parameter :: planet_name(nbodies) = &
Expand All @@ -122,14 +123,16 @@ subroutine add_sun_and_planets(nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,epoch)
'saturn ', &
'uranus ', &
'neptune'/)
integer :: i
integer :: i,ierr

nerr = 0
do i=1,nbodies
if (present(epoch)) then
call add_body(planet_name(i),nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,epoch=epoch)
call add_body(planet_name(i),nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,ierr,epoch=epoch)
else
call add_body(planet_name(i),nptmass,xyzmh_ptmass,vxyz_ptmass,mtot)
call add_body(planet_name(i),nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,ierr)
endif
nerr = nerr + ierr
enddo

end subroutine add_sun_and_planets
Expand All @@ -140,10 +143,11 @@ end subroutine add_sun_and_planets
! from the JPL server
!+
!----------------------------------------------------------------
subroutine add_dwarf_planets(nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,epoch)
subroutine add_dwarf_planets(nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,nerr,epoch)
integer, intent(inout) :: nptmass
real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:)
real, intent(in) :: mtot
integer, intent(out) :: nerr
character(len=*), intent(in), optional :: epoch
integer, parameter :: nbodies = 5
character(len=*), parameter :: planet_name(nbodies) = &
Expand All @@ -152,14 +156,16 @@ subroutine add_dwarf_planets(nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,epoch)
'eris ', &
'makemake', &
'haumea '/)
integer :: i
integer :: i,ierr

nerr = 0
do i=1,nbodies
if (present(epoch)) then
call add_body(planet_name(i),nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,epoch=epoch)
call add_body(planet_name(i),nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,ierr,epoch=epoch)
else
call add_body(planet_name(i),nptmass,xyzmh_ptmass,vxyz_ptmass,mtot)
call add_body(planet_name(i),nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,ierr)
endif
nerr = nerr + ierr
enddo

end subroutine add_dwarf_planets
Expand All @@ -170,7 +176,7 @@ end subroutine add_dwarf_planets
! from the JPL server
!+
!----------------------------------------------------------------
subroutine add_body(body_name,nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,epoch)
subroutine add_body(body_name,nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,ierr,epoch)
use ephemeris, only:get_ephemeris,nelem
use units, only:umass,udist
use physcon, only:gg,km,solarm,solarr,earthm,au
Expand All @@ -180,13 +186,15 @@ subroutine add_body(body_name,nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,epoch)
integer, intent(inout) :: nptmass
real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:)
real, intent(in) :: mtot
integer, intent(out) :: ierr
character(len=*), intent(in), optional :: epoch
real :: elems(nelem)
real :: xyz_tmp(size(xyzmh_ptmass(:,1)),2),vxyz_tmp(3,2),gm_cgs
real :: mbody,rbody,a,e,inc,O,w,f
integer :: ierr,ntmp
integer :: ntmp
logical :: got_elem(nelem)

ierr = 0
if (trim(adjustl(lcase(body_name)))=='sun') then
!
! add the Sun, ideally would add here its motion
Expand All @@ -207,6 +215,7 @@ subroutine add_body(body_name,nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,epoch)
endif
if (ierr > 0 .or. .not.any(got_elem)) then
print "(a)",' ERROR: could not read ephemeris data for '//body_name
ierr = 1
return ! skip if error reading ephemeris file
endif
a = elems(1)*km/udist
Expand Down
19 changes: 12 additions & 7 deletions src/setup/setup_solarsystem.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,8 +43,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
grainsize,graindens,ndustlarge,ndusttypes,ndustsmall,ihacc,igas
use setbinary, only:set_binary
use units, only:set_units,umass,udist,unit_density,unit_velocity,utime,in_code_units
use physcon, only:solarm,au,pi,km,solarr,ceresm,earthm,earthr,days
use io, only:master,fatal
use physcon, only:solarm,pi,au,km,solarr,ceresm,earthm,earthr,days
use io, only:master,fatal,warning
use timestep, only:tmax,dtmax
use centreofmass, only:reset_centreofmass
use setsolarsystem,only:set_minor_planets,add_sun_and_planets,add_body
Expand All @@ -62,7 +62,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
real, intent(inout) :: time
character(len=20), intent(in) :: fileprefix
real, intent(out) :: vxyzu(:,:)
integer :: ierr,i
integer :: ierr,i,nerr
!integer :: values(8),year,month,day
logical :: iexist
real :: period,semia,mtot,dx
Expand Down Expand Up @@ -98,7 +98,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
!
! set units
!
call set_units(mass=solarm,dist=au,G=1.d0)
call set_units(mass=solarm,dist=km,G=1.d0)
!
! general parameters
!
Expand Down Expand Up @@ -141,12 +141,15 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
!
! add the planets
!
call add_sun_and_planets(nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,epoch)
ierr = 0
call add_sun_and_planets(nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,nerr,epoch)
if (nerr > 0) ierr = ierr + nerr
!
! add the bringer of death
!
if (np_apophis > 0) then
call add_body('apophis',nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,epoch)
call add_body('apophis',nptmass,xyzmh_ptmass,vxyz_ptmass,mtot,nerr,epoch)
if (nerr > 0) call warning('apophis','missing some information')

r_apophis = xyzmh_ptmass(5,nptmass)
m_apophis = 4./3.*pi*(rho_0/unit_density)*r_apophis**3
Expand Down Expand Up @@ -188,7 +191,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
!
call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)

if (ierr /= 0) call fatal('setup','ERROR during setup')
if (ierr /= 0) call fatal('setup','ERRORS during setup')

end subroutine setpart

Expand Down Expand Up @@ -236,7 +239,9 @@ subroutine read_setupfile(filename,ierr)
call read_inopt(dtmax_in,'dtmax_in',db,errcount=nerr)
call read_inopt(asteroids,'asteroids',db,errcount=nerr)
call read_inopt(np_apophis,'np_apophis',db,errcount=nerr)
print*,' prev epoch = ',trim(epoch)
call read_inopt(epoch,'epoch',db,errcount=nerr)
print*,' GOT EPOCH=',epoch
call close_db(db)

if (nerr > 0) then
Expand Down
6 changes: 2 additions & 4 deletions src/utils/utils_ephemeris.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ subroutine construct_horizons_api_url(object,url,ierr,epoch)
integer, intent(out) :: ierr
character(len=*), intent(in), optional :: epoch
character(len=8) :: cmd
character(len=10) :: start_epoch,end_epoch
character(len=20) :: start_epoch,end_epoch
integer :: values(8),year,month,day

ierr = 0
Expand Down Expand Up @@ -124,8 +124,7 @@ subroutine construct_horizons_api_url(object,url,ierr,epoch)

! take the input epoch but only if it parses into YYYY-MM-DD correctly
if (present(epoch)) then
!print "(a)", ' parsing EPOCH='//trim(epoch)
read(epoch,"(i4.4,1x,i2.2,1x,i2.2)") year,month,day
read(epoch(1:10),"(i4.4,1x,i2.2,1x,i2.2)",iostat=ierr) year,month,day
if (ierr == 0) then
start_epoch = epoch
else
Expand All @@ -136,7 +135,6 @@ subroutine construct_horizons_api_url(object,url,ierr,epoch)

! end one day later
write(end_epoch,"(i4.4,'-',i2.2,'-',i2.2)") year,month,day+1

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)// &
Expand Down
Loading