Merge branch 'master' of https://github.com/danieljprice/phantom

danieljprice · Jul 13, 2020 · 331cb33 · 331cb33
2 parents e1678df + 68e841f
commit 331cb33
Show file tree

Hide file tree

Showing 4 changed files with 41 additions and 39 deletions.
diff --git a/AUTHORS b/AUTHORS
@@ -15,8 +15,8 @@ Mark Hutchison <markahutch@gmail.com>
 David Liptai <dliptai@hotmail.com>
 Terrence Tricco <ttricco@cita.utoronto.ca>
 Mike Lau <mike.lau@monash.edu>
-Christophe Pinte <christophe.pinte@univ-grenoble-alpes.fr>
 Bec Nealon <rebecca.nealon@leicester.ac.uk>
+Christophe Pinte <christophe.pinte@univ-grenoble-alpes.fr>
 Sergei Biriukov <svbiriukov@gmail.com>
 Giovanni Dipierro <giovanni.dipierro@leicester.ac.uk>
 Hauke Worpel <hworpel@aip.de>

diff --git a/src/main/inject_asteroidwind.f90 b/src/main/inject_asteroidwind.f90
@@ -11,17 +11,18 @@
 !
 !  REFERENCES: None
 !
-!  OWNER: David Liptai
+!  OWNER: Bec Nealon
 !
 !  $Id$
 !
 !  RUNTIME PARAMETERS:
+!    dndt_type     -- injection rate (0=const, 1=cos(t), 2=r^(-2))
 !    mdot          -- mass injection rate in grams/second
 !    npartperorbit -- particle injection rate in particles/binary orbit
 !    vlag          -- percentage lag in velocity of wind
-!    dndt_type     -- injection rate (0=const, 1=cos(t), 2=r^(-2))
 !
-!  DEPENDENCIES: infile_utils, io, part, partinject, physcon, random, units
+!  DEPENDENCIES: externalforces, infile_utils, io, options, part,
+!    partinject, physcon, random, units
 !+
 !--------------------------------------------------------------------------
 module inject
@@ -221,7 +222,7 @@ function dndt_func(t,r,ra,rp,ecc,scaling) result(dndt)
  end select
 
  dndt = dndt*scaling
- 
+
 end function dndt_func
 
 !-----------------------------------------------------------------------
@@ -362,13 +363,13 @@ real function get_E(period,ecc,deltat)
  M_guess = M_ref - 2.*tol
 
  do while (abs(M_ref - M_guess) > tol)
-   M_guess = E_guess - ecc*sin(E_guess)
-   if (M_guess > M_ref) then
-      E_right = E_guess
-   else
-      E_left = E_guess
-   endif
-   E_guess = 0.5*(E_left + E_right)
+    M_guess = E_guess - ecc*sin(E_guess)
+    if (M_guess > M_ref) then
+       E_right = E_guess
+    else
+       E_left = E_guess
+    endif
+    E_guess = 0.5*(E_left + E_right)
  enddo
 
  get_E = E_guess

diff --git a/src/setup/setup_asteroidwind.f90 b/src/setup/setup_asteroidwind.f90
@@ -20,14 +20,15 @@
 !    eccentricity  -- eccentricity
 !    gastemp       -- gas temperature in K
 !    hacc1         -- white dwarf (sink) accretion radius (solar radii)
+!    ipot          -- wd modelled by sink (0) or externalforce(1)
 !    m1            -- mass of white dwarf (solar mass)
 !    m2            -- mass of asteroid (ceres mass)
 !    norbits       -- number of orbits
 !    rasteroid     -- radius of asteroid (km)
 !    semia         -- semi-major axis (solar radii)
 !
-!  DEPENDENCIES: eos, infile_utils, inject, io, part, physcon, setbinary,
-!    spherical, timestep, units
+!  DEPENDENCIES: eos, externalforces, infile_utils, inject, io, options,
+!    part, physcon, setbinary, spherical, timestep, units
 !+
 !--------------------------------------------------------------------------
 module setup
@@ -146,33 +147,33 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
 ! point mass for the asteroid
 
  if (ipot == 0) then
-   !
-   !--Set a binary orbit given the desired orbital parameters
-   !
-   call set_binary(m1,m2,semia,ecc,hacc1,hacc2,xyzmh_ptmass,vxyz_ptmass,nptmass,ierr)
-   xyzmh_ptmass(ihsoft,2) = rasteroid ! Asteroid radius softening
+    !
+    !--Set a binary orbit given the desired orbital parameters
+    !
+    call set_binary(m1,m2,semia,ecc,hacc1,hacc2,xyzmh_ptmass,vxyz_ptmass,nptmass,ierr)
+    xyzmh_ptmass(ihsoft,2) = rasteroid ! Asteroid radius softening
 
  else
 
-   !
-   !--Set the asteroid on orbit around the fixed potential
-   !
-   mass1                = m1
-   accradius1           = hacc1
-   iexternalforce       = 11
-   call update_externalforce(iexternalforce,time,0.)
-
-   ! Orbit and position
-   nptmass = 1
-   xyzmh_ptmass(1:3,1) = (/semia*(1. + ecc),0.,0./)
-   vxyz_ptmass(1:3,1)  = (/0.,sqrt(semia*(1.-ecc**2)*(m1+m2))/xyzmh_ptmass(1,1),0./)
-
-   xyzmh_ptmass(4,1)      = m2
-   xyzmh_ptmass(ihacc,1)  = hacc2        ! Asteroid should not accrete
-   xyzmh_ptmass(ihsoft,1) = rasteroid ! Asteroid radius softening
- endif   
-
- ! both	of these are reset in the first	call to	inject_particles
+    !
+    !--Set the asteroid on orbit around the fixed potential
+    !
+    mass1                = m1
+    accradius1           = hacc1
+    iexternalforce       = 11
+    call update_externalforce(iexternalforce,time,0.)
+
+    ! Orbit and position
+    nptmass = 1
+    xyzmh_ptmass(1:3,1) = (/semia*(1. + ecc),0.,0./)
+    vxyz_ptmass(1:3,1)  = (/0.,sqrt(semia*(1.-ecc**2)*(m1+m2))/xyzmh_ptmass(1,1),0./)
+
+    xyzmh_ptmass(4,1)      = m2
+    xyzmh_ptmass(ihacc,1)  = hacc2        ! Asteroid should not accrete
+    xyzmh_ptmass(ihsoft,1) = rasteroid ! Asteroid radius softening
+ endif
+
+ ! both        of these are reset in the first        call to        inject_particles
  massoftype(igas) = 1.0
  hfact = 1.2
  call inject_particles(time,0.,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,npartoftype,dtinj)

diff --git a/src/setup/setup_star.f90 b/src/setup/setup_star.f90
@@ -48,7 +48,7 @@
 !    ui_coef            -- specific internal energy (units of GM/R)
 !    unsoftened_profile -- Path to MESA profile for softening
 !    use_exactN         -- find closest particle number to np
-!    write_rho_to_file  -- write density profile to
+!    write_rho_to_file  -- write density profile to file
 !
 !  DEPENDENCIES: centreofmass, dim, eos, eos_idealplusrad,
 !    extern_densprofile, externalforces, infile_utils, io, kernel, options,