Skip to content

Commit

Permalink
(CE-analysis) divv output choice of 0 gives no output
Browse files Browse the repository at this point in the history
  • Loading branch information
themikelau committed Aug 18, 2023
1 parent 0b80203 commit 2ee6dcd
Showing 1 changed file with 12 additions and 8 deletions.
20 changes: 12 additions & 8 deletions src/utils/analysis_common_envelope.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1411,10 +1411,10 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
'13) JstarS' !option to calculate JstarS

quantities_to_calculate = (/1,2,4,5/)
call prompt('Choose first quantity to compute ',quantities_to_calculate(1),1,Nquantities)
call prompt('Choose second quantity to compute ',quantities_to_calculate(2),1,Nquantities)
call prompt('Choose third quantity to compute ',quantities_to_calculate(3),1,Nquantities)
call prompt('Choose fourth quantity to compute ',quantities_to_calculate(4),1,Nquantities)
call prompt('Choose first quantity to compute ',quantities_to_calculate(1),0,Nquantities)
call prompt('Choose second quantity to compute ',quantities_to_calculate(2),0,Nquantities)
call prompt('Choose third quantity to compute ',quantities_to_calculate(3),0,Nquantities)
call prompt('Choose fourth quantity to compute ',quantities_to_calculate(4),0,Nquantities)
endif

! Calculations performed outside loop over particles
Expand All @@ -1424,7 +1424,7 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
com_vxyz = 0.
do k=1,4
select case (quantities_to_calculate(k))
case(1,2,3,6,8,9,13) ! Nothing to do
case(0,1,2,3,6,8,9,13) ! Nothing to do
case(4,5,11,12) ! Fractional difference between gas and orbital omega
if (quantities_to_calculate(k) == 4 .or. quantities_to_calculate(k) == 5) then
com_xyz = (xyzmh_ptmass(1:3,1)*xyzmh_ptmass(4,1) + xyzmh_ptmass(1:3,2)*xyzmh_ptmass(4,2)) &
Expand Down Expand Up @@ -1501,6 +1501,9 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
endif
quant(k,i) = JstarS

case(0) ! Skip
quant(k,i) = 0.

case(1,9) ! Total energy (kin + pot + therm)
rhopart = rhoh(xyzh(4,i), particlemass)
call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i))
Expand Down Expand Up @@ -2508,7 +2511,7 @@ subroutine planet_profile(num,dumpfile,particlemass,xyzh,vxyzu)
real, dimension(3) :: planet_com,planet_vcom,vnorm,ri,Rvec
real, allocatable :: R(:),z(:),rho(:)

call get_planetIDs(nplanet,planetIDs)
if (dump_number ==0 ) call get_planetIDs(nplanet,planetIDs)
allocate(R(nplanet),z(nplanet),rho(nplanet))

! Find highest density in planet
Expand All @@ -2526,7 +2529,7 @@ subroutine planet_profile(num,dumpfile,particlemass,xyzh,vxyzu)
vnorm = planet_vcom / sqrt(dot_product(planet_vcom,planet_vcom))

! Write to file
file_name = trim(dumpfile)//".dist"
file_name = trim(dumpfile)//".planetpart"
open(newunit=iu, file=file_name, status='replace')

! Record R and z cylindrical coordinates w.r.t. planet_com
Expand All @@ -2535,7 +2538,8 @@ subroutine planet_profile(num,dumpfile,particlemass,xyzh,vxyzu)
z(i) = dot_product(ri, vnorm)
Rvec = ri - z(i)*vnorm
R(i) = sqrt(dot_product(Rvec,Rvec))
write(iu,"(es13.6,2x,es13.6,2x,es13.6)") R(i),z(i),rho(i)
! write(iu,"(es13.6,2x,es13.6,2x,es13.6)") R(i),z(i),rho(i)
write(iu,"(es13.6,2x,es13.6,2x,es13.6,2x,es13.6,2x,es13.6)") xyzh(1,i),xyzh(2,i),xyzh(3,i),rho(i),vxyzu(4,i)
enddo

close(unit=iu)
Expand Down

0 comments on commit 2ee6dcd

Please sign in to comment.