Merge pull request #596 from danieljprice/set-binary-fix

(set-binary) update npartoftype array when particle types are set according to a specific star in a binary

src/setup/set_star.f90

@@ -315,6 +315,8 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,&
     do i=npart_old+1,npart
        call set_particle_type(i,itype+istar_offset)
     enddo
+    npartoftype(itype+istar_offset) = npartoftype(itype+istar_offset) + npart - npart_old
+    npartoftype(igas) = npartoftype(igas) - (npart - npart_old)
  endif
  !
  ! Print summary to screen
@@ -399,7 +401,7 @@ end subroutine set_stars
 !+
 !-----------------------------------------------------------------------
 subroutine shift_star(npart,xyz,vxyz,x0,v0,itype,corotate)
- use part,        only:get_particle_type,set_particle_type,igas
+ use part,        only:get_particle_type,set_particle_type,igas,npartoftype
  use vectorutils, only:cross_product3D
  integer, intent(in) :: npart
  real, intent(inout) :: xyz(:,:),vxyz(:,:)
@@ -432,6 +434,8 @@ subroutine shift_star(npart,xyz,vxyz,x0,v0,itype,corotate)
        if (mytype /= itype+istar_offset) cycle over_parts
        ! reset type back to gas
        call set_particle_type(i,igas)
+       npartoftype(itype+istar_offset) = npartoftype(itype+istar_offset) - 1
+       npartoftype(igas) = npartoftype(igas) + 1
     endif
     xyz(1:3,i) = xyz(1:3,i) + x0(:)
     vxyz(1:3,i) = vxyz(1:3,i) + v0(:)
@@ -809,24 +813,18 @@ subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label)
 
  select case(star%iprofile)
  case(imesa)
-    ! core softening options
-    call read_inopt(star%isinkcore,'isinkcore'//trim(c),db,errcount=nerr)
-
-    if (star%isinkcore) then
-       call read_inopt(lcore_lsun,'lcore'//trim(c),db,errcount=nerr,min=0.,err=ierr)
-       if (ierr==0) star%lcore = lcore_lsun*real(solarl/unit_luminosity)
-    endif
-
     call read_inopt(star%isoftcore,'isoftcore'//trim(c),db,errcount=nerr,min=0)
 
     if (star%isoftcore <= 0) then ! sink particle core without softening
+       call read_inopt(star%isinkcore,'isinkcore'//trim(c),db,errcount=nerr)
        if (star%isinkcore) then
           call read_inopt(mcore_msun,'mcore'//trim(c),db,errcount=nerr,min=0.,err=ierr)
           if (ierr==0) star%mcore = mcore_msun*real(solarm/umass)
           call read_inopt(hsoft_rsun,'hsoft'//trim(c),db,errcount=nerr,min=0.,err=ierr)
           if (ierr==0) star%hsoft = hsoft_rsun*real(solarr/udist)
        endif
     else
+       star%isinkcore = .true.
        call read_inopt(star%outputfilename,'outputfilename'//trim(c),db,errcount=nerr)
        if (star%isoftcore==2) then
           star%isofteningopt=3
@@ -844,6 +842,11 @@ subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label)
           if (ierr==0) star%mcore = mcore_msun*real(solarm/umass)
        endif
     endif
+
+    if (star%isinkcore) then
+       call read_inopt(lcore_lsun,'lcore'//trim(c),db,errcount=nerr,min=0.,err=ierr)
+       if (ierr==0) star%lcore = lcore_lsun*real(solarl/unit_luminosity)
+    endif
  case(ievrard)
     call read_inopt(star%ui_coef,'ui_coef'//trim(c),db,errcount=nerr,min=0.)
  case(:0)

src/setup/set_star_utils.f90

@@ -289,7 +289,7 @@ subroutine set_star_density(lattice,id,master,rmin,Rstar,Mstar,hfact,&
  !
  ! set particle type as gas particles
  !
- npartoftype(igas) = npart   ! npart is number on this thread only
+ npartoftype(igas) = npartoftype(igas) + npart - npart_old   ! npart is number on this thread only
  do i=npart_old+1,npart_old+npart
     call set_particle_type(i,igas)
  enddo