From e75ba5cec5350aeb3c49d7e98a3e838f7c5d7c2e Mon Sep 17 00:00:00 2001
From: Mike Lau <mike.lau@monash.edu>
Date: Wed, 5 Aug 2026 14:37:00 +0200
Subject: [PATCH 1/9] (extern_gwinspiral) change variable name nstar to
 nstar_gw

---
 src/main/extern_gwinspiral.f90  | 48 ++++++++++++++++-----------------
 src/utils/analysis_NSmerger.f90 | 10 +++----
 2 files changed, 29 insertions(+), 29 deletions(-)

diff --git a/src/main/extern_gwinspiral.f90 b/src/main/extern_gwinspiral.f90
index 6f13a2658..73a0a84f7 100644
--- a/src/main/extern_gwinspiral.f90
+++ b/src/main/extern_gwinspiral.f90
@@ -28,7 +28,7 @@ module extern_gwinspiral
  ! local variables
  !
  integer, private :: n_threshhold
- integer, public  :: Nstar(2) = 0  ! give default value in case dump header not read
+ integer, public  :: Nstar_gw(2) = 0  ! give default value in case dump header not read
  real,    private :: fstar1_coef,fstar2_coef
  real,    private :: com(3),comstar1(3),comstar2(3),vcomstar1(3),vcomstar2(3),fstar1(3),fstar2(3)
  logical, private :: isseparate = .true.
@@ -60,12 +60,12 @@ subroutine initialise_gwinspiral(npart,nptmass,ierr)
 
  ierr = 0
  nerr = 0
- if (Nstar(1) > 0) then
-    write(*,"(2(a,i8))") ' Initialising inspiral two stars scenario, Nstar_1 = ',Nstar(1),' Nstar_2 = ', Nstar(2)
- elseif (Nstar(1)==0 .and. Nstar(2)==0 .and. nptmass==2) then
+ if (Nstar_gw(1) > 0) then
+    write(*,"(2(a,i8))") ' Initialising inspiral two stars scenario, Nstar_1 = ',Nstar_gw(1),' Nstar_2 = ', Nstar_gw(2)
+ elseif (Nstar_gw(1)==0 .and. Nstar_gw(2)==0 .and. nptmass==2) then
     write(*,"(2(a,i8))") ' Initialising inspiral on two sink particles'
- elseif (Nstar(2)==0 .and. nptmass==1) then
-    write(*,"(2(a,i8))") ' Initialising inspiral star-particle scenario, Nstar_1 = ',Nstar(1)
+ elseif (Nstar_gw(2)==0 .and. nptmass==1) then
+    write(*,"(2(a,i8))") ' Initialising inspiral star-particle scenario, Nstar_1 = ',Nstar_gw(1)
  else
     ierr = 1
     isseparate = .false.
@@ -100,21 +100,21 @@ subroutine gw_still_inspiralling(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptma
  fstar2_coef = 0.
  stopped_now = .false.
  if ( isseparate ) then
-    if (Nstar(1) == 0 .and. nptmass==2) then
+    if (Nstar_gw(1) == 0 .and. nptmass==2) then
        comstar1  = xyzmh_ptmass(1:3,2)
        mstar1    = xyzmh_ptmass(4,2)
        vcomstar1 = vxyz_ptmass(1:3,2)
     else
-       call get_centreofmass(comstar1,vcomstar1,Nstar(1),xyzh(:,1:Nstar(1)),&
-                          vxyzu(:,1:Nstar(1)),mass=mstar1)
+       call get_centreofmass(comstar1,vcomstar1,Nstar_gw(1),xyzh(:,1:Nstar_gw(1)),&
+                          vxyzu(:,1:Nstar_gw(1)),mass=mstar1)
     endif
-    if (Nstar(2) == 0 .and. nptmass>=1) then
+    if (Nstar_gw(2) == 0 .and. nptmass>=1) then
        comstar2  = xyzmh_ptmass(1:3,1)
        mstar2    = xyzmh_ptmass(4,1)
        vcomstar2 = vxyz_ptmass(1:3,1)
     else
-       call get_centreofmass(comstar2,vcomstar2,Nstar(2),xyzh(:,Nstar(1)+1:npart),&
-                          vxyzu(:,Nstar(1)+1:npart),mass=mstar2)
+       call get_centreofmass(comstar2,vcomstar2,Nstar_gw(2),xyzh(:,Nstar_gw(1)+1:npart),&
+                          vxyzu(:,Nstar_gw(1)+1:npart),mass=mstar2)
     endif
     com  = (comstar1*mstar1 + comstar2*mstar2)/(mstar1+mstar2)
     dirstar1 = comstar1 - com ! The directional vector to star 1
@@ -124,11 +124,11 @@ subroutine gw_still_inspiralling(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptma
     !
     k1 = 0
     !$omp parallel default(none) &
-    !$omp shared(nstar,xyzh,com,dirstar1) &
+    !$omp shared(Nstar_gw,xyzh,com,dirstar1) &
     !$omp private(i,dx,dy,dz,dir) &
     !$omp reduction(+:k1)
     !$omp do
-    do i=1,nstar(1)
+    do i=1,Nstar_gw(1)
        dx  = xyzh(1,i) - com(1)
        dy  = xyzh(2,i) - com(2)
        dz  = xyzh(3,i) - com(3)
@@ -143,11 +143,11 @@ subroutine gw_still_inspiralling(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptma
     !
     k2 = 0
     !$omp parallel default(none) &
-    !$omp shared(nstar,npart,xyzh,com,dirstar1) &
+    !$omp shared(Nstar_gw,npart,xyzh,com,dirstar1) &
     !$omp private(i,dx,dy,dz,dir) &
     !$omp reduction(+:k2)
     !$omp do
-    do i=nstar(1)+1,nstar(1)+nstar(2)
+    do i=Nstar_gw(1)+1,Nstar_gw(1)+Nstar_gw(2)
        dx  = xyzh(1,i) - com(1)
        dy  = xyzh(2,i) - com(2)
        dz  = xyzh(3,i) - com(3)
@@ -220,20 +220,20 @@ subroutine get_gw_force_i(i,fextxi,fextyi,fextzi,phi)
  real,    intent(inout) :: fextxi,fextyi,fextzi,phi
 
  if (i > 0 .and. isseparate ) then
-    if (i <= nstar(1)) then
+    if (i <= Nstar_gw(1)) then
        fextxi = fstar1(1)
        fextyi = fstar1(2)
        fextzi = fstar1(3)
-    elseif (nstar(2) > 0) then
+    elseif (Nstar_gw(2) > 0) then
        fextxi = fstar2(1)
        fextyi = fstar2(2)
        fextzi = fstar2(3)
     endif
- elseif (i == -1 .and. nstar(2)==0 .and. isseparate) then
+ elseif (i == -1 .and. Nstar_gw(2)==0 .and. isseparate) then
     fextxi = fstar2(1)  ! acceleration applied to sink particle 1 (star 2)
     fextyi = fstar2(2)
     fextzi = fstar2(3)
- elseif (i == -2 .and. nstar(1)==0 .and. isseparate) then
+ elseif (i == -2 .and. Nstar_gw(1)==0 .and. isseparate) then
     fextxi = fstar1(1) ! acceleration applied to sink particle 2 (star 1)
     fextyi = fstar1(2)
     fextzi = fstar1(3)
@@ -283,8 +283,8 @@ subroutine write_headeropts_gwinspiral(hdr,ierr)
  integer,      intent(out)   :: ierr
 
  ierr = 0
- call add_to_header(Nstar(1),'Nstar_1',hdr,ierr)
- call add_to_header(Nstar(2),'Nstar_2',hdr,ierr)
+ call add_to_header(Nstar_gw(1),'Nstar_1',hdr,ierr)
+ call add_to_header(Nstar_gw(2),'Nstar_2',hdr,ierr)
 
 end subroutine write_headeropts_gwinspiral
 
@@ -301,8 +301,8 @@ subroutine read_headeropts_gwinspiral(hdr,nptmass,ierr)
  integer :: ierr1,ierr2
 
  ierr  = 0
- call extract('Nstar_1',Nstar(1),hdr,ierr1)
- call extract('Nstar_2',Nstar(2),hdr,ierr2)
+ call extract('Nstar_1',Nstar_gw(1),hdr,ierr1)
+ call extract('Nstar_2',Nstar_gw(2),hdr,ierr2)
 
  if (ierr1 /= 0 .or. ierr2 /= 0) then
     ! if there are two sink particles and Nstar_1
diff --git a/src/utils/analysis_NSmerger.f90 b/src/utils/analysis_NSmerger.f90
index c20ea975b..cd81b2d8d 100644
--- a/src/utils/analysis_NSmerger.f90
+++ b/src/utils/analysis_NSmerger.f90
@@ -25,7 +25,7 @@ module analysis
  use physcon,         only: pi
  use centreofmass,    only: get_centreofmass
  use readwrite_dumps, only:opened_full_dump
- use extern_gwinspiral, only:Nstar
+ use extern_gwinspiral, only:Nstar_gw
  implicit none
  character(len=20), parameter, public :: analysistype = 'NSmerger'
  !
@@ -75,8 +75,8 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
     !
     !--Get the index range for each star
     ! (Note from DJP: This should now be automatically read from the dump header)
-    if (Nstar(1) <= 0) call fatal('analysis_NSmerger','Require Nstar(1) > 0 in header of dump file')
-    if (Nstar(2) <= 0) call fatal('analysis_NSmerger','Require Nstar(2) > 0 in header of dump file')
+    if (Nstar_gw(1) <= 0) call fatal('analysis_NSmerger','Require Nstar_gw(1) > 0 in header of dump file')
+    if (Nstar_gw(2) <= 0) call fatal('analysis_NSmerger','Require Nstar_gw(2) > 0 in header of dump file')
     !
     !--Prompt for the density_cut off
     if (choice > 1) then
@@ -148,8 +148,8 @@ subroutine trace_com(dumpfile,xyzh,vxyzu,time,npart,iunit)
  endif
  !
  !--Get centre of masses of the stars
- call get_centreofmass(com1,vcom1,nstar(1),xyzh(:,1:nstar(1)),vxyzu(:,1:nstar(1)))
- call get_centreofmass(com2,vcom2,nstar(2),xyzh(:,nstar(1)+1:npart),vxyzu(:,nstar(1)+1:npart))
+ call get_centreofmass(com1,vcom1,Nstar_gw(1),xyzh(:,1:Nstar_gw(1)),vxyzu(:,1:Nstar_gw(1)))
+ call get_centreofmass(com2,vcom2,Nstar_gw(2),xyzh(:,Nstar_gw(1)+1:npart),vxyzu(:,Nstar_gw(1)+1:npart))
  !
  rad = sqrt( (com1(1)-com2(1))**2 + (com1(2)-com2(2))**2 + (com1(3)-com2(3))**2 )
  !

From 6b6e41d3775dfd797a37e6ea061b85f7a04d7f6e Mon Sep 17 00:00:00 2001
From: Mike Lau <mike.lau@monash.edu>
Date: Wed, 5 Aug 2026 14:38:09 +0200
Subject: [PATCH 2/9] (extern_gwinspiral) add boost factor fow GW inspiral

---
 src/main/extern_gwinspiral.f90 | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/main/extern_gwinspiral.f90 b/src/main/extern_gwinspiral.f90
index 73a0a84f7..287ea7195 100644
--- a/src/main/extern_gwinspiral.f90
+++ b/src/main/extern_gwinspiral.f90
@@ -24,6 +24,7 @@ module extern_gwinspiral
  ! Runtime parameters
  !
  real,    public  :: stopratio = 0.005
+ real,    public  :: gw_boostfac = 1.  ! multiplier to force applied to stars
  !
  ! local variables
  !
@@ -203,8 +204,8 @@ subroutine get_gw_force()
     !
     ! Compute the drag force vectors for each star
     !
-    fstar1 = fstar1_coef * vcomstar1 / (vstar1sq*separation**5)
-    fstar2 = fstar2_coef * vcomstar2 / (vstar2sq*separation**5)
+    fstar1 = gw_boostfac * fstar1_coef * vcomstar1 / (vstar1sq*separation**5)
+    fstar2 = gw_boostfac * fstar2_coef * vcomstar2 / (vstar2sq*separation**5)
  endif
 
 end subroutine get_gw_force
@@ -255,6 +256,7 @@ subroutine write_options_gwinspiral(iunit)
  integer, intent(in) :: iunit
 
  call write_inopt(stopratio,'stop_ratio','ratio of particles crossing CoM to indicate a merger',iunit)
+ call write_inopt(gw_boostfac,'gw_boostfac','multiplicative factor to boost GW force',iunit)
 
 end subroutine write_options_gwinspiral
 
@@ -269,6 +271,7 @@ subroutine read_options_gwinspiral(db,nerr)
  integer,      intent(inout) :: nerr
 
  call read_inopt(stopratio,'stop_ratio',db,errcount=nerr,min=0.,max=1.)
+ call read_inopt(gw_boostfac,'gw_boostfac',db,errcount=nerr,min=0.)
 
 end subroutine read_options_gwinspiral
 

From 8d4ab262bbd7a31b02caf507c6c81c66e69d0405 Mon Sep 17 00:00:00 2001
From: Mike Lau <mike.lau@monash.edu>
Date: Wed, 5 Aug 2026 14:58:03 +0200
Subject: [PATCH 3/9] (extern_gwinspiral) require G=1 in code units

---
 src/main/externalforces.f90 | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/main/externalforces.f90 b/src/main/externalforces.f90
index 4cf0f32af..659fab2f6 100644
--- a/src/main/externalforces.f90
+++ b/src/main/externalforces.f90
@@ -818,7 +818,8 @@ subroutine initialise_externalforces(iexternalforce,ierr)
  end select
 
  select case(iexternalforce)
- case(iext_star,iext_binary,iext_corot_binary,iext_prdrag,iext_spiral,iext_lensethirring,iext_einsteinprec,iext_gnewton)
+ case(iext_star,iext_binary,iext_corot_binary,iext_prdrag,iext_spiral,iext_lensethirring,&
+      iext_einsteinprec,iext_gnewton,iext_gwinspiral)
     !
     !--check that G=1 in code units
     !

From 024559f1f9088d3e516a120f39ef51c5b14cc5a1 Mon Sep 17 00:00:00 2001
From: Mike Lau <mike.lau@monash.edu>
Date: Wed, 5 Aug 2026 15:00:12 +0200
Subject: [PATCH 4/9] (extern_gwinspiral) add option to binary setup

---
 src/setup/setup_binary.f90 | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/setup/setup_binary.f90 b/src/setup/setup_binary.f90
index ff63b9392..5b3650833 100644
--- a/src/setup/setup_binary.f90
+++ b/src/setup/setup_binary.f90
@@ -30,7 +30,7 @@ module setup
  implicit none
  public :: setpart
 
- logical :: relax,write_rho_to_file,corotate
+ logical :: relax,write_rho_to_file,corotate,gwinspiral
  integer, parameter :: max_stars = max_hier_levels
  type(star_t)  :: star(max_stars)
  type(orbit_t) :: orbit
@@ -53,7 +53,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,&
  use setorbit,       only:set_defaults_orbit,set_orbit
  use sethierarchical,only:set_hierarchical,set_hierarchical_default_options,hs,generate_hierarchy_string
  use options,        only:iexternalforce,alphau
- use externalforces, only:iext_corotate,iext_geopot,iext_star,omega_corotate,mass1,accradius1
+ use externalforces, only:iext_corotate,iext_geopot,iext_star,iext_gwinspiral,omega_corotate,mass1,accradius1
+ use extern_gwinspiral, only:Nstar_gw
  use io,             only:master,fatal
  use setstar,        only:set_defaults_stars,set_stars,shift_stars
  use eos,            only:X_in,Z_in,use_var_comp
@@ -109,6 +110,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,&
  call set_hierarchical_default_options()
  relax = .true.
  corotate = .false.
+ gwinspiral = .false.
  use_var_comp = .false.
  write_rho_to_file = .true.
 
@@ -181,6 +183,10 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,&
     accradius1 = xyzmh_ptmass(ihacc,nptmass+1)
     xyzmh_ptmass(:,nptmass) = xyzmh_ptmass(:,nptmass+1)
     vxyz_ptmass(:,nptmass) = vxyz_ptmass(:,nptmass+1)
+ elseif (nstar == 2 .and. (gwinspiral .or. iexternalforce==iext_gwinspiral)) then
+    Nstar_gw(1) = star(1)%np
+    Nstar_gw(2) = npart - star(1)%np
+    iexternalforce = iext_gwinspiral
  elseif (nstar >= 1 .and. set_oblateness .and. star(1)%iprofile == 0) then
     ! set J2 for sink particle 1 to be equal to oblateness of Saturn
     xyzmh_ptmass(iJ2,1) = 0.01629
@@ -222,6 +228,7 @@ subroutine write_setupfile(filename)
     call write_hierarchical_setupfile(iunit,nstar)
  elseif (nstar == 2) then
     call write_inopt(corotate,'corotate','set stars in corotation',iunit)
+    call write_inopt(gwinspiral,'gwinspiral','set gravitational wave inspiral',iunit)
     call write_options_orbit(orbit,iunit)
 
     write(iunit,"(/,a)") '# timestepping'
@@ -261,6 +268,7 @@ subroutine read_setupfile(filename,ierr)
     call read_hierarchical_setupfile(db,nerr,nstar)
  elseif (nstar == 2) then
     call read_inopt(corotate,'corotate',db,errcount=nerr)
+    call read_inopt(gwinspiral,'gwinspiral',db,errcount=nerr)
     m1 = in_code_units(star(1)%m,ierr,unit_type='mass')
     m2 = in_code_units(star(2)%m,ierr1,unit_type='mass')
     if (ierr /= 0 .or. ierr1 /= 0) then

From e03a05a3e19428ca09963c610910d7f4ec322fd4 Mon Sep 17 00:00:00 2001
From: Mike Lau <mike.lau@monash.edu>
Date: Wed, 5 Aug 2026 15:39:42 +0200
Subject: [PATCH 5/9] (extern_gwinspiral) set gwinspiral in moddump_binary

---
 src/utils/moddump_binary.f90     | 16 ++++++++++++++--
 src/utils/moddump_sinkbinary.f90 |  4 ++--
 2 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/src/utils/moddump_binary.f90 b/src/utils/moddump_binary.f90
index 0d5c5fbca..07b868b5c 100644
--- a/src/utils/moddump_binary.f90
+++ b/src/utils/moddump_binary.f90
@@ -35,9 +35,10 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  use centreofmass,      only:reset_centreofmass,get_centreofmass
  use prompting,         only:prompt
  use options,           only:iexternalforce
- use externalforces,    only:omega_corotate,iext_corotate
+ use externalforces,    only:omega_corotate,iext_corotate,iext_gwinspiral
  use extern_corotate,   only:icompanion_grav,companion_xpos,companion_mass,primarycore_xpos,&
                              primarycore_mass,primarycore_hsoft,hsoft
+ use extern_gwinspiral, only:Nstar_gw
  use infile_utils,      only:open_db_from_file,inopts,read_inopt,close_db
  use table_utils,       only:yinterp
  use readwrite_mesa,    only:read_mesa
@@ -62,7 +63,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  real                      :: time2,hfact2
  real                      :: xyzmh1_stash(nsinkproperties),xyzmh2_stash(nsinkproperties),vxyz1_stash(3),vxyz2_stash(3)
  real, allocatable         :: r(:),den(:),pres(:),temp(:),enitab(:),Xfrac(:),Yfrac(:),mu(:),m(:)
- logical                   :: use_corotating_frame,iprimary_grav_ans
+ logical                   :: use_corotating_frame,iprimary_grav_ans,gwinspiral
  character(len=20)         :: filename = 'binary.in'
  character(len=100)        :: densityfile,dumpname
  type(inopts), allocatable :: db(:)
@@ -178,6 +179,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
        hacc = 0.
        hsoft = 0.
        use_corotating_frame = .false.
+       gwinspiral = .false.
 
        call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)
        call delete_dead_or_accreted_particles(npart,npartoftype)  !removes the dead or accreted particles for a correct total mass computation
@@ -206,6 +208,9 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
           call prompt('Enter softening length for companion', hsoft, 0.)
        endif
        call prompt('Do you want to transform to a corotating frame and simulate corotating binary?', use_corotating_frame)
+       call prompt('Do you want to add gravitational radiation reaction?', gwinspiral)
+
+       if (gwinspiral) iexternalforce = iext_gwinspiral
 
        ! set the binary
        if (use_corotating_frame) then
@@ -302,6 +307,11 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
              vxyz_ptmass(1:3,nptmass1+nptmass2) = vxyz2_stash(1:3)
           endif
 
+          if (gwinspiral) then
+             Nstar_gw(1) = nstar1
+             Nstar_gw(2) = nstar2
+          endif
+
        else
           nptmass = nptmass1 + 1
           xyzmh_ptmass(1:3,nptmass) = xyzmh2_stash(1:3)
@@ -309,6 +319,8 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
           xyzmh_ptmass(4,nptmass) = m2
           xyzmh_ptmass(ihacc,nptmass) = hacc
           xyzmh_ptmass(ihsoft,nptmass) = hsoft
+
+          if (gwinspiral) Nstar_gw(1) = nstar1
        endif
 
        if (nptmass1 == 1) then
diff --git a/src/utils/moddump_sinkbinary.f90 b/src/utils/moddump_sinkbinary.f90
index 4fc04bbc8..e20db4177 100644
--- a/src/utils/moddump_sinkbinary.f90
+++ b/src/utils/moddump_sinkbinary.f90
@@ -32,7 +32,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  use timestep, only:dtmax,tmax
  use options,   only:iexternalforce
  use externalforces, only:iext_gwinspiral
- use extern_gwinspiral, only:Nstar
+ use extern_gwinspiral, only:Nstar_gw
  integer, intent(inout) :: npart
  integer, intent(inout) :: npartoftype(:)
  real,    intent(inout) :: massoftype(:)
@@ -77,7 +77,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)
 
  if (iexternalforce==iext_gwinspiral) then
-    Nstar(1) = npart
+    Nstar_gw(1) = npart
  endif
 
  period = 2.*pi*sqrt(a**3/(m1 + m2))

From a63bdad067f5e31375a7298d755fcfd6b44c7978 Mon Sep 17 00:00:00 2001
From: Mike Lau <mike.lau@monash.edu>
Date: Wed, 5 Aug 2026 17:28:21 +0200
Subject: [PATCH 6/9] (extern_gwinspiral) fixes

---
 src/utils/moddump_binary.f90     | 3 +--
 src/utils/moddump_binarystar.f90 | 6 +++---
 2 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/src/utils/moddump_binary.f90 b/src/utils/moddump_binary.f90
index 07b868b5c..8aacd96c5 100644
--- a/src/utils/moddump_binary.f90
+++ b/src/utils/moddump_binary.f90
@@ -208,8 +208,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
           call prompt('Enter softening length for companion', hsoft, 0.)
        endif
        call prompt('Do you want to transform to a corotating frame and simulate corotating binary?', use_corotating_frame)
-       call prompt('Do you want to add gravitational radiation reaction?', gwinspiral)
-
+       if (.not. use_corotating_frame) call prompt('Do you want to add gravitational radiation reaction?', gwinspiral)
        if (gwinspiral) iexternalforce = iext_gwinspiral
 
        ! set the binary
diff --git a/src/utils/moddump_binarystar.f90 b/src/utils/moddump_binarystar.f90
index f6c49fe13..f4cfff682 100644
--- a/src/utils/moddump_binarystar.f90
+++ b/src/utils/moddump_binarystar.f90
@@ -745,11 +745,11 @@ end subroutine add_vrotational
 !  Save nstar so it can be properly written to the header
 !
 subroutine save_nstar(Nstar1,Nstar2)
- use extern_gwinspiral, only:Nstar
+ use extern_gwinspiral, only:Nstar_gw
  integer, intent(in) :: Nstar1,Nstar2
 
- Nstar(1) = Nstar1
- Nstar(2) = Nstar2
+ Nstar_gw(1) = Nstar1
+ Nstar_gw(2) = Nstar2
 
 end subroutine save_nstar
 !-----------------------------------------------------------------------

From f1f6a616f634e08ed1829949bc7dd7de4357874f Mon Sep 17 00:00:00 2001
From: Mike Lau <mike.lau@monash.edu>
Date: Fri, 7 Aug 2026 10:06:04 +0200
Subject: [PATCH 7/9] remove redundant moddump_sinkbinary

---
 src/utils/moddump_sinkbinary.f90 | 90 --------------------------------
 1 file changed, 90 deletions(-)
 delete mode 100644 src/utils/moddump_sinkbinary.f90

diff --git a/src/utils/moddump_sinkbinary.f90 b/src/utils/moddump_sinkbinary.f90
deleted file mode 100644
index e20db4177..000000000
--- a/src/utils/moddump_sinkbinary.f90
+++ /dev/null
@@ -1,90 +0,0 @@
-!--------------------------------------------------------------------------!
-! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. !
-! Copyright (c) 2007-2026 The Authors (see AUTHORS)                        !
-! See LICENCE file for usage and distribution conditions                   !
-! http://phantomsph.github.io/                                             !
-!--------------------------------------------------------------------------!
-module moddump
-!
-! Add a sink particle binary to the dump
-!
-! :References: None
-!
-! :Owner: Mike Lau
-!
-! :Runtime parameters: None
-!
-! :Dependencies: centreofmass, extern_gwinspiral, externalforces, options,
-!   part, physcon, prompting, setbinary, timestep, units
-!
- implicit none
- character(len=*), parameter, public :: moddump_flags = ''
-
-contains
-
-subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
- use part,      only:xyzmh_ptmass,vxyz_ptmass,nptmass,igas
- use setbinary, only:set_binary
- use units,     only:umass,udist
- use physcon,   only:solarm,solarr,pi
- use prompting, only:prompt
- use centreofmass, only:reset_centreofmass
- use timestep, only:dtmax,tmax
- use options,   only:iexternalforce
- use externalforces, only:iext_gwinspiral
- use extern_gwinspiral, only:Nstar_gw
- integer, intent(inout) :: npart
- integer, intent(inout) :: npartoftype(:)
- real,    intent(inout) :: massoftype(:)
- real,    intent(inout) :: xyzh(:,:),vxyzu(:,:)
- real :: m1,m2,a,e,racc,period
- integer :: i,ierr
-
- ! find current mass from existing particles
- m2 = npart*massoftype(igas)
- print*,' Mass of star from existing file in Msun = ',m2*umass/solarm
- call reset_centreofmass(npart,xyzh,vxyzu)
- !
- ! set up a sink particle binary
- !
- m1 = 1.4
- call prompt('enter mass of point mass in Msun ',m1,0.)
- print*,' code unit of distance in Rsun = ',udist/solarr
-
- a = 2500.
- call prompt('enter semi-major axis in code units',a,0.)
-
- racc = a/200.
- call prompt('enter accretion radius of point mass',racc,0.01)
- print*,' accretion radius in solar radii = ',racc*udist/solarr
-
- e = 0.
- call prompt('enter eccentricity ',e,0.)
-
- !
- ! add sink particle binary
- !
- nptmass = 0
- call set_binary(m1,m2,a,e,racc,racc,xyzmh_ptmass,vxyz_ptmass,nptmass,ierr)
- !
- ! delete one of the sink particles and replace it with our polytrope
- !
- nptmass = nptmass - 1
- do i=1,npart
-    xyzh(1:3,i) = xyzh(1:3,i) + xyzmh_ptmass(1:3,2)
-    vxyzu(1:3,i) = vxyzu(1:3,i) + vxyz_ptmass(1:3,2)
- enddo
- call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)
-
- if (iexternalforce==iext_gwinspiral) then
-    Nstar_gw(1) = npart
- endif
-
- period = 2.*pi*sqrt(a**3/(m1 + m2))
- print*,' orbital period = ',period
- tmax = 1000.*period
- dtmax = 0.1*period
-
-end subroutine modify_dump
-
-end module moddump

From e296f2e4d99f3b487609a504c64600b12d0cc149 Mon Sep 17 00:00:00 2001
From: Mike Lau <mike.lau@monash.edu>
Date: Fri, 7 Aug 2026 10:16:28 +0200
Subject: [PATCH 8/9] (gwinspiral) fix build failure

---
 src/main/extern_gwinspiral.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/main/extern_gwinspiral.f90 b/src/main/extern_gwinspiral.f90
index 287ea7195..fd882ed11 100644
--- a/src/main/extern_gwinspiral.f90
+++ b/src/main/extern_gwinspiral.f90
@@ -271,7 +271,7 @@ subroutine read_options_gwinspiral(db,nerr)
  integer,      intent(inout) :: nerr
 
  call read_inopt(stopratio,'stop_ratio',db,errcount=nerr,min=0.,max=1.)
- call read_inopt(gw_boostfac,'gw_boostfac',db,errcount=nerr,min=0.)
+ call read_inopt(gw_boostfac,'gw_boostfac',db,errcount=nerr,min=0.,default=1.)
 
 end subroutine read_options_gwinspiral
 

From 0ce46a2bf2aca49b4e8cb49ad4822fdca83292b0 Mon Sep 17 00:00:00 2001
From: Mike Lau <mike.lau@monash.edu>
Date: Fri, 7 Aug 2026 16:07:37 +0200
Subject: [PATCH 9/9] (gwinspiral) fix test failure

---
 src/setup/setup_binary.f90   | 2 +-
 src/utils/moddump_binary.f90 | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/setup/setup_binary.f90 b/src/setup/setup_binary.f90
index 5b3650833..365a1e4e7 100644
--- a/src/setup/setup_binary.f90
+++ b/src/setup/setup_binary.f90
@@ -268,7 +268,7 @@ subroutine read_setupfile(filename,ierr)
     call read_hierarchical_setupfile(db,nerr,nstar)
  elseif (nstar == 2) then
     call read_inopt(corotate,'corotate',db,errcount=nerr)
-    call read_inopt(gwinspiral,'gwinspiral',db,errcount=nerr)
+    call read_inopt(gwinspiral,'gwinspiral',db,errcount=nerr,default=.false.)
     m1 = in_code_units(star(1)%m,ierr,unit_type='mass')
     m2 = in_code_units(star(2)%m,ierr1,unit_type='mass')
     if (ierr /= 0 .or. ierr1 /= 0) then
diff --git a/src/utils/moddump_binary.f90 b/src/utils/moddump_binary.f90
index 8aacd96c5..a933a7cf5 100644
--- a/src/utils/moddump_binary.f90
+++ b/src/utils/moddump_binary.f90
@@ -307,8 +307,8 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
           endif
 
           if (gwinspiral) then
-             Nstar_gw(1) = nstar1
-             Nstar_gw(2) = nstar2
+             Nstar_gw(1) = nstar2
+             Nstar_gw(2) = nstar1
           endif
 
        else
@@ -319,7 +319,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
           xyzmh_ptmass(ihacc,nptmass) = hacc
           xyzmh_ptmass(ihsoft,nptmass) = hsoft
 
-          if (gwinspiral) Nstar_gw(1) = nstar1
+          if (gwinspiral) Nstar_gw(1) = npart
        endif
 
        if (nptmass1 == 1) then