Merge branch 'master' of https://bitbucket.org/danielprice/phantom

AUTHORS

@@ -43,6 +43,7 @@ Rebecca Nealon <rebecca.nealon@leicester.ac.uk>
 Zachary Pellow <zpel1@student.monash.edu>
 Alison Young <ayoung@astro.ex.ac.uk>
 Cox, Samuel <sc676@leicester.ac.uk>
+Daniel Price <dprice@cantab.net>
 Hauke WOrpel <hworpel@aip.de>
 Jorge Cuadra <jcuadra@astro.puc.cl>
 Lionel Siess <Lionel.Siess@ulb.ac.be>

src/utils/analysis_disc.f90

@@ -127,7 +127,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyz,pmass,npart,time,iunit)
     print*,'Resetting assume_Ltot_is_same_as_zaxis=.true. in analysis'
  endif
 
- call disc_analysis(xyzh,vxyz,npart,pmass,time,nr,rmin,rmax,H_R,G,M_star,q_index,&
+ call disc_analysis(xyzh,vxyz,npart,pmass,time,nr,rmin,rmax,G,M_star,&
                      tilt,tilt_acc,twist,twistprev,psi,H,rad,h_smooth,sigma,unitlx,unitly,unitlz,&
                      Lx,Ly,Lz,ecc,ninbin,assume_Ltot_is_same_as_zaxis,xyzmh_ptmass,vxyz_ptmass,nptmass)
 

src/utils/analysis_disc_planet.f90

@@ -28,7 +28,7 @@ module analysis
  character(len=20), parameter, public :: analysistype = 'disc_planet'
  public :: do_analysis
 
- integer, parameter :: nr = 50
+ integer, parameter :: nr = 100
 
  private
 
@@ -39,7 +39,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyz,pmass,npart,time,iunit)
  use physcon, only:pi
  use part,    only:xyzmh_ptmass,vxyz_ptmass,nptmass
  use options, only:iexternalforce
- use centreofmass, only:get_total_angular_momentum
+ use centreofmass, only:get_total_angular_momentum,reset_centreofmass
  use vectorutils, only:cross_product3D,rotatevec
  character(len=*), intent(in) :: dumpfile
  real,             intent(inout) :: xyzh(:,:),vxyz(:,:)
@@ -55,17 +55,18 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyz,pmass,npart,time,iunit)
  real :: unitlx(nr),unitly(nr),unitlz(nr),ecc(nr)
  real :: psi(nr),tilt_acc(nr),Lx(nr),Ly(nr),Lz(nr),twistprev(nr)
  real :: L_tot(3),L_p(3),L_inner_mag,L_outer_mag
- real :: L_p_mag,L_ratio_inner,L_ratio_outer,e_planet,ecc_planet(3)
- real :: rad_planet,twist_inner,twist_outer,tilt_inner,tilt_outer,minclin,minclout,mannulus
- real :: m_red,mu,rotate_about_y,rotate_about_z,planet_mass,pos_planet(3),vel_planet(3)
- real :: temp(3),temp_mag,term(3),tilt_planet,twist_planet,L_tot_mag
- real :: unitl_in(3),unitl_out(3),eff_tilt
+ real :: L_p_mag,L_ratio_inner,L_ratio_outer,e_sinks,ecc_sinks(3)
+ real :: rad_sinks,twist_inner,twist_outer,tilt_inner,tilt_outer,minclin,minclout,mannulus
+ real :: m_red,mu,rotate_about_y,rotate_about_z,sinks_mass,pos_sinks(3),vel_sinks(3)
+ real :: temp(3),temp_mag,term(3),tilt_sinks,twist_sinks,L_tot_mag
+ real :: unitl_in(3),unitl_out(3),eff_tilt,eff_outertilt
+ real :: break_radius,break_disc_at,eff_innertilt,dlhs,drhs
  integer :: ninbin(nr),n_count_inner,n_count_outer,nptmassinit
  logical :: assume_Ltot_is_same_as_zaxis
 
  integer, parameter :: iparams = 10
  integer, parameter :: iprec   = 24
- integer, parameter :: iplanet = 23
+ integer, parameter :: isinks = 23
  logical :: do_precession,ifile,iexist
 
  do_precession = .false.
@@ -116,7 +117,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyz,pmass,npart,time,iunit)
 
  assume_Ltot_is_same_as_zaxis = .false.
 
- call disc_analysis(xyzh,vxyz,npart,pmass,time,nr,rmin,rmax,H_R,G,M_star,q_index,&
+ call disc_analysis(xyzh,vxyz,npart,pmass,time,nr,rmin,rmax,G,M_star,&
                      tilt,tilt_acc,twist,twistprev,psi,H,rad,h_smooth,sigma,unitlx,unitly,unitlz,&
                      Lx,Ly,Lz,ecc,ninbin,assume_Ltot_is_same_as_zaxis,xyzmh_ptmass,vxyz_ptmass,nptmass)
 
@@ -172,6 +173,11 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyz,pmass,npart,time,iunit)
     nptmassinit = 1
  endif
 
+ ! Prepare for calculations below
+ if (nptmass > 0) then
+    call reset_centreofmass(npart,xyzh,vxyz,nptmass,xyzmh_ptmass,vxyz_ptmass)
+ endif
+
  ! Prepare for rotations later on
  call get_total_angular_momentum(xyzh,vxyz,npart,L_tot,xyzmh_ptmass,vxyz_ptmass,nptmass)
 
@@ -187,49 +193,54 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyz,pmass,npart,time,iunit)
 ! Set up the radius array
  dr = (rmax-rmin)/real(nr-1)
 
-! Calculating and printing information for the planet
+! Calculating and printing information for the sinks
  if (nptmass>nptmassinit) then
     do i=nptmassinit+1,nptmass
-       write(filename,"(a,i3.3)")"planet_",i-1
+       write(filename,"(a,i3.3)")"sinks_",i-1
        inquire(file=filename,exist=iexist)
        if (.not.iexist .or. numfile == 0) then
-          open(iplanet,file=filename,status="replace")
-          write(iplanet,"('#',27(1x,'[',i2.2,1x,a11,']',2x))") &
+          open(isinks,file=filename,status="replace")
+          write(isinks,"('#',13(1x,'[',i2.2,1x,a11,']',2x))") &
                1,'time', &
                2,'rad(sph)', &
-               3,'x', &
-               4,'y', &
-               5,'z', &
-               6,'lx', &
-               7,'ly', &
-               8,'lz', &
-               9,'vx', &
-               10,'vy', &
-               11,'vz', &
-               12,'Tilt p', &
-               13,'Twist p', &
-               14,'Ecc p', &
-               15,'Tilt in', &
-               16,'Tilt out', &
-               17,'Twist in', &
-               18,'Twist out', &
-               19,'Lin/Lp', &
-               20,'Lout/Lp',&
-               21,'lx in', &
-               22,'ly in', &
-               23,'lz in', &
-               24,'lx out', &
-               25,'ly out', &
-               26,'lz out', &
-               27,'eff tilt'
+               3,'Tilt sink', &
+               4,'Twist sink', &
+               5,'Tilt ID', &
+               6,'Tilt OD', &
+               7,'Twist ID', &
+               8,'Twist OD', &
+               9,'Lin/Lp', &
+               10,'Lout/Lp',&
+               11,'Rel. misalignment', &
+               12,'Sink-ID', &
+               13,'Sink-OD'
+
        else
-          open(iplanet,file=filename,status="old",position="append")
+          open(isinks,file=filename,status="old",position="append")
        endif
 
-       planet_mass = xyzmh_ptmass(4,i)
-       pos_planet = xyzmh_ptmass(1:3,i) - xyzmh_ptmass(1:3,1)
-       rad_planet = sqrt(dot_product(pos_planet,pos_planet))
-       vel_planet = vxyz_ptmass(1:3,i) - vxyz_ptmass(1:3,1)
+
+       ! Properties of sink
+       sinks_mass = xyzmh_ptmass(4,i)
+       pos_sinks = xyzmh_ptmass(1:3,i) - xyzmh_ptmass(1:3,1)
+       rad_sinks = sqrt(dot_product(pos_sinks,pos_sinks))
+       vel_sinks = vxyz_ptmass(1:3,i) - vxyz_ptmass(1:3,1)
+
+       ! Identify if there is a break in the disc based on sigma
+       break_radius = -1
+       do ii=2,nr-2
+          dlhs = (sigma(ii) - sigma(ii-1))/(rad(ii) - rad(ii-1))
+          drhs = (sigma(ii+1) - sigma(ii))/(rad(ii+1) - rad(ii))
+          if (dlhs < 0. .and. drhs > 0.) then
+             dlhs = (sigma(ii-1) - sigma(ii-2))/(rad(ii-1) - rad(ii-2))
+             drhs = (sigma(ii+2) - sigma(ii+1))/(rad(ii+2) - rad(ii+1))
+             if (dlhs < 0. .and. drhs > 0.) break_radius = rad(ii)
+          endif
+       enddo
+
+       ! If there is a break in the disc use this to decide the inner and outer disc
+       break_disc_at = rad_sinks
+       if (break_radius > rad(1)) break_disc_at = break_radius
 
        tilt_inner = 0.
        tilt_outer = 0.
@@ -247,15 +258,15 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyz,pmass,npart,time,iunit)
        do ii=1,nr
           if (ninbin(ii) > 0) then
              mannulus = 2.*pi*rad(ii)*sigma(ii)*dr
-             if (rad(ii) < rad_planet) then
-                n_count_inner = n_count_inner + 1
+             if (rad(ii) < break_disc_at) then
+                n_count_inner = n_count_inner + ninbin(ii)
                 minclin = minclin + mannulus
                 tilt_inner = tilt_inner + tilt(ii)*mannulus
                 twist_inner = twist_inner + twist(ii)*mannulus
                 L_inner_mag = L_inner_mag + sqrt(Lx(ii)**2 + Ly(ii)**2 + Lz(ii)**2)
                 unitl_in = unitl_in + (/unitlx(ii),unitly(ii),unitlz(ii)/)*mannulus
              else
-                n_count_outer = n_count_outer + 1
+                n_count_outer = n_count_outer + ninbin(ii)
                 minclout = minclout + mannulus
                 tilt_outer = tilt_outer + tilt(ii)*mannulus
                 twist_outer = twist_outer + twist(ii)*mannulus
@@ -265,53 +276,81 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyz,pmass,npart,time,iunit)
           endif
        enddo
 
-       ! Average the tilt and twist inside and outside the planet orbit
-       tilt_inner = tilt_inner/minclin!real(n_count_inner)
-       tilt_outer = tilt_outer/minclout!real(n_count_outer)
-       twist_inner = twist_inner/minclin!real(n_count_inner)
-       twist_outer = twist_outer/minclout!real(n_count_outer)
-       unitl_in = unitl_in/minclin!real(n_count_inner)
-       unitl_out = unitl_out/minclout!real(n_count_outer)
-
-       ! Rotate planet vector such that Ltot is parallel to z-axis
-       call cross_product3D(xyzmh_ptmass(1:3,i),vxyz_ptmass(1:3,i),L_p)
-       L_p = L_p*planet_mass
+       ! Average the tilt and twist of inner and outer disc
+       if (n_count_inner > 0) then
+          tilt_inner = tilt_inner/minclin
+          twist_inner = twist_inner/minclin
+          unitl_in = unitl_in/minclin
+       else
+          tilt_inner = 0.
+          twist_inner = 0.
+          unitl_in = 0.
+       endif
+       if (n_count_outer > 0) then
+          tilt_outer = tilt_outer/minclout
+          twist_outer = twist_outer/minclout
+          unitl_out = unitl_out/minclout
+       else
+          tilt_outer = 0.
+          twist_outer = 0.
+          unitl_out = 0.
+       endif
+
+       call cross_product3D(pos_sinks,vel_sinks,L_p)
+       L_p = L_p*sinks_mass
        L_p_mag = sqrt(dot_product(L_p,L_p))
-       call rotatevec(L_p,(/0.,0.,1.0/),-rotate_about_z)
-       call rotatevec(L_p,(/0.,1.0,0./),rotate_about_y)
-       tilt_planet = acos(L_p(3)/L_p_mag)
+
+       if (.not.assume_Ltot_is_same_as_zaxis) then
+          ! Rotate sinks vector such that Ltot is parallel to z-axis
+          call rotatevec(L_p,(/0.,0.,1.0/),-rotate_about_z)
+          call rotatevec(L_p,(/0.,1.0,0./),rotate_about_y)
+       endif
+       tilt_sinks = acos(L_p(3)/L_p_mag)
 
        ! Angular momentum ratios
-       L_p_mag = sqrt(dot_product(L_p,L_p))
        L_ratio_inner = L_inner_mag/L_p_mag
        L_ratio_outer = L_outer_mag/L_p_mag
 
        ! Calculate the effective tilt between the inner and outer disc
        eff_tilt = acos(dot_product(unitl_in,unitl_out))
+       ! Calculate the effective tilt between discs and binary
+       ! NB: here we assume that each sink in the binary has the same
+       ! direction of AM
+       eff_innertilt = acos(dot_product(unitl_in,L_p/L_p_mag))
+       eff_outertilt = acos(dot_product(unitl_out,L_p/L_p_mag))
 
-       ! For now, measure twist of planet as in disc
-       twist_planet = atan2(L_p(2),L_p(1))
-       if (twist_planet < 0.) twist_planet = twist_planet + 2.*pi
+       ! For now, measure twist of sinks as in disc
+       twist_sinks = atan2(L_p(2),L_p(1))
+       if (twist_sinks < 0.) twist_sinks = twist_sinks + 2.*pi
 
        ! Calculate the eccentricity of the point mass
        ! These are calculated using the relative position and velocities,
        ! rotations above do not affect these
-       m_red = planet_mass*xyzmh_ptmass(4,1)/(planet_mass + xyzmh_ptmass(4,1))
-       mu = 1.0*(planet_mass + xyzmh_ptmass(4,1))  !mu = GM, G=1
-       call cross_product3D(vel_planet,L_p/planet_mass,term)
-       ecc_planet = term/mu - pos_planet/rad_planet
+       m_red = sinks_mass*xyzmh_ptmass(4,1)/(sinks_mass + xyzmh_ptmass(4,1))
+       mu = 1.0*(sinks_mass + xyzmh_ptmass(4,1))  !mu = GM, G=1
+       call cross_product3D(vel_sinks,L_p/sinks_mass,term)
+       ecc_sinks = term/mu - pos_sinks/rad_sinks
 
-       if (ecc_planet(1) < -1.0) then
+       if (ecc_sinks(1) < -1.0) then
           print*,'Warning: eccentricity is undefined, returning e=1.'
           term = 0.
        endif
-       e_planet = sqrt(dot_product(ecc_planet,ecc_planet))
+       e_sinks = sqrt(dot_product(ecc_sinks,ecc_sinks))
+
+       write(isinks,'(13(es18.10,1X))') time,break_disc_at, &
+            tilt_sinks,twist_sinks,tilt_inner,tilt_outer, &
+            twist_inner,twist_outer,L_ratio_inner,L_ratio_outer, &
+            eff_tilt,eff_innertilt,eff_outertilt
+
+       if (numfile == 0) then
+          print*,'For sink',i,'located at',rad_sinks
+          print*,'and for a reference radius of',break_disc_at
+          print*,'AM RATIOS: Lin/Lp=',L_ratio_inner,' and Lout/Lp=',L_ratio_outer
+          print*,'Mass inside',real(n_count_inner)*pmass
+          print*,'Mass outside',real(n_count_outer)*pmass
+       endif
 
-       write(iplanet,'(27(es18.10,1X))') time, rad_planet,xyzmh_ptmass(1:3,i), &
-            L_p(:),vel_planet(1:3),tilt_planet,twist_planet,e_planet, &
-            tilt_inner,tilt_outer,twist_inner,twist_outer,L_ratio_inner,L_ratio_outer, &
-            unitl_in(:),unitl_out(:),eff_tilt
-       close(iplanet)
+       close(isinks)
     enddo
  endif
 

src/utils/grid2pdf.f90

@@ -18,7 +18,7 @@
 !
 !  USAGE: grid2pdf format filename(s)
 !
-!  DEPENDENCIES: dim, io, io_grid, pdfs, rhomach
+!  DEPENDENCIES: io, io_grid, pdfs, rhomach
 !+
 !--------------------------------------------------------------------------
 program grid2pdf

src/utils/io_grid.f90

@@ -18,7 +18,7 @@
 !
 !  RUNTIME PARAMETERS: None
 !
-!  DEPENDENCIES: fileutils, hdf5utils, io
+!  DEPENDENCIES: fileutils, hdf5utils, io, readwrite_griddata
 !+
 !--------------------------------------------------------------------------
 module io_grid

src/utils/utils_disc.f90

@@ -40,7 +40,7 @@ module discanalysisutils
 !+
 !----------------------------------------------------------------
 
-subroutine disc_analysis(xyzh,vxyz,npart,pmass,time,nbin,rmin,rmax,H_R,G,M_star,q_index,&
+subroutine disc_analysis(xyzh,vxyz,npart,pmass,time,nbin,rmin,rmax,G,M_star,&
                      tilt,tilt_acc,twist,twistprev,psi,H,bin,h_smooth,sigma,unitlx,unitly,unitlz,Lx,Ly,Lz,&
                      ecc,ninbin,assume_Ltot_is_same_as_zaxis,xyzmh_ptmass,vxyz_ptmass,nptmass)
  use physcon,        only:pi
@@ -51,7 +51,7 @@ subroutine disc_analysis(xyzh,vxyz,npart,pmass,time,nbin,rmin,rmax,H_R,G,M_star,
  use prompting,      only:prompt
  real, intent(inout)              :: xyzh(:,:),vxyz(:,:),pmass,time
  integer, intent(in)              :: nbin,npart
- real, intent(in)                 :: rmin,rmax,H_R,G,M_star,q_index
+ real, intent(in)                 :: rmin,rmax,G,M_star
  logical, intent(in)              :: assume_Ltot_is_same_as_zaxis
  real, optional, intent(inout)    :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:)
  integer, optional, intent(inout) :: nptmass