From 234ae2006cd9bf7b583c0875a57ccf544bccf5c6 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Tue, 4 Mar 2025 11:38:46 +1100 Subject: [PATCH] (solarsystem) better error handling if ephemeris files not downloaded; fix issue with datetime string for epoch being converted automatically --- src/main/utils_infiles.f90 | 3 ++- src/setup/set_solarsystem.f90 | 29 +++++++++++++++++++---------- src/setup/setup_solarsystem.f90 | 19 ++++++++++++------- src/utils/utils_ephemeris.f90 | 6 ++---- 4 files changed, 35 insertions(+), 22 deletions(-) diff --git a/src/main/utils_infiles.f90 b/src/main/utils_infiles.f90 index 8f197386a..d34c370d5 100644 --- a/src/main/utils_infiles.f90 +++ b/src/main/utils_infiles.f90 @@ -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 diff --git a/src/setup/set_solarsystem.f90 b/src/setup/set_solarsystem.f90 index 7a0575898..f1e5d8a5f 100644 --- a/src/setup/set_solarsystem.f90 +++ b/src/setup/set_solarsystem.f90 @@ -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) = & @@ -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 @@ -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) = & @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/setup/setup_solarsystem.f90 b/src/setup/setup_solarsystem.f90 index 7ca66a9ae..def14d04f 100644 --- a/src/setup/setup_solarsystem.f90 +++ b/src/setup/setup_solarsystem.f90 @@ -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 @@ -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 @@ -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 ! @@ -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 @@ -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 @@ -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 diff --git a/src/utils/utils_ephemeris.f90 b/src/utils/utils_ephemeris.f90 index ffd88e7b1..58ceed89b 100644 --- a/src/utils/utils_ephemeris.f90 +++ b/src/utils/utils_ephemeris.f90 @@ -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 @@ -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 @@ -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)// &