gravity_binary.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: gravity_binary.f90 #
5!# #
6!# Copyright (C) 2010-2024 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Björn Sperling <sperling@astrophysik.uni-kiel.de> #
9!# Anna Feiler <afeiler@astrophysik.uni-kiel.de> #
10!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
11!# #
12!# This program is free software; you can redistribute it and/or modify #
13!# it under the terms of the GNU General Public License as published by #
14!# the Free Software Foundation; either version 2 of the License, or (at #
15!# your option) any later version. #
16!# #
17!# This program is distributed in the hope that it will be useful, but #
18!# WITHOUT ANY WARRANTY; without even the implied warranty of #
19!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
20!# NON INFRINGEMENT. See the GNU General Public License for more #
21!# details. #
22!# #
23!# You should have received a copy of the GNU General Public License #
24!# along with this program; if not, write to the Free Software #
25!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
26!# #
27!#############################################################################
40!----------------------------------------------------------------------------!
52!----------------------------------------------------------------------------!
62 USE common_dict
63#ifdef PARALLEL
64#ifdef HAVE_MPI_MOD
65 USE mpi
66#endif
67#endif
68 IMPLICIT NONE
69#ifdef PARALLEL
70#ifdef HAVE_MPIF_H
71 include 'mpif.h'
72#endif
73#endif
74 !--------------------------------------------------------------------------!
75 PRIVATE
77 CLASS(marray_base), ALLOCATABLE :: pot_prim,pot_sec
78 REAL, DIMENSION(:,:,:), POINTER :: r_sec
79 REAL, DIMENSION(:,:,:,:), POINTER :: posvec_sec
80 REAL, DIMENSION(:,:,:,:), POINTER :: posvec_sec_tmp
81 REAL, DIMENSION(:,:,:,:), POINTER :: fr_sec
82 REAL, DIMENSION(:,:,:,:,:), POINTER :: fposvec_sec
83 REAL, POINTER :: mass2
84 REAL :: excent
85 REAL :: semaaxis
86 REAL :: period
87 REAL :: eps1,eps2
88 REAL :: switchon2
89 REAL :: omega_rot
90 CONTAINS
91 PROCEDURE :: initgravity_binary
93 PROCEDURE :: updatepositions
94 PROCEDURE :: infogravity
95 PROCEDURE :: setoutput
96 PROCEDURE :: getmass_primary
97 PROCEDURE :: getmass_secondary
99 final :: finalize
100 END TYPE
101
102 PUBLIC :: &
103 ! types
105 ! module procedures
106 kepler
107 !--------------------------------------------------------------------------!
108
109CONTAINS
115 SUBROUTINE initgravity_binary(this,Mesh,Physics,config,IO)
118 IMPLICIT NONE
119 !------------------------------------------------------------------------!
120 CLASS(gravity_binary), INTENT(INOUT) :: this
121 CLASS(mesh_base), INTENT(IN) :: Mesh
122 CLASS(physics_base), INTENT(IN) :: Physics
123 TYPE(dict_typ),POINTER :: config
125 TYPE(dict_typ),POINTER :: IO
126 !------------------------------------------------------------------------!
127 INTEGER :: mesh_rot
128 INTEGER :: err
129 !------------------------------------------------------------------------!
130 ! init base class
131 CALL this%InitGravity(mesh,physics,"binary point masses",config,io)
132
133 SELECT TYPE(physics)
134 TYPE IS(physics_euler)
135 ! do nothing
136 TYPE IS(physics_eulerisotherm)
137 ! do nothing
138 CLASS DEFAULT
139 CALL this%Error("InitGravity_binary", &
140 "Physics not supported for binary gravity term")
141 END SELECT
142
143 ALLOCATE(&
144 this%mass,this%massloss,this%accrate,this%mass2,this%pos(3,2),this%r0(3), &
145 this%omega2, this%omega, &
146 this%r_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX),&
147 this%r_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX),&
148 this%posvec_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
149 this%posvec_prim_tmp(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
150 this%posvec_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
151 this%posvec_sec_tmp(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
152 this%pot,this%pot_prim,this%pot_sec, &
153 this%fr_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3), &
154 this%fr_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3), &
155 this%fposvec_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3,3),&
156 this%fposvec_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3,3),&
157 ! last two entries (EAST,NORTH)x(dim1,dim2)
158 stat = err)
159 IF (err.NE.0) CALL this%Error("gravity_pointmass::InitGravity_pointmass", "Unable allocate memory!")
160
161 ! initialize potentials
162 this%pot = marray_base(4)
163 this%pot_prim = marray_base(4)
164 this%pot_sec = marray_base(4)
165 this%omega = marray_base()
166 this%omega2 = marray_base(2)
167 this%pot%data1d(:) = 0.0
168 this%pot%data1d(:) = 0.0
169 this%pot_prim%data1d(:) = 0.0
170 this%pot_sec%data1d(:) = 0.0
171 this%omega2%data1d(:) = 0.0
172
173 ! initialize omega with zero except for ghost cells to avoid zero division
174 WHERE (mesh%without_ghost_zones%mask1d(:))
175 this%omega%data1d(:) = 0.0
176 ELSEWHERE
177 this%omega%data1d(:) = 1.0
178 END WHERE
179
180
181 ! mass of primary component
182 CALL getattr(config, "mass1", this%mass, 1.0)
183 CALL setattr(io, "mass1", ref(this%mass))
184
185 ! mass of secondary component
186 CALL getattr(config, "mass2", this%mass2, 1.0)
187 CALL setattr(io, "mass2", ref(this%mass2))
188
189 ! excentricity
190 CALL getattr(config, "excentricity", this%excent, 0.0)
191
192 ! semi major axis
193 CALL getattr(config, "semimajoraxis", this%semaaxis, -1.0)
194
195 ! Softening parameter
196 CALL getattr(config, "softening1", this%eps1, 0.0)
197
198 CALL getattr(config, "softening2", this%eps2, 0.0)
199
200 ! soft switch on for primary component
201 CALL getattr(config, "switchon1", this%switchon, -1.0)
202
203 ! soft switch on for secondary component
204 CALL getattr(config, "switchon2", this%switchon2, -1.0)
205
206 ! rotating reference frame
207 CALL getattr(config, "omega_rot", this%omega_rot, 0.0)
208
209 ! check whether rotating mesh is activated
210 CALL getattr(config, "mesh_rot", mesh_rot, 0)
211 IF (mesh_rot.EQ.1) THEN
212 this%omega_rot = mesh%OMEGA
213 this%r0(:) = mesh%rotcent(:)
214 IF (mesh%omega.EQ.0.0) &
215 CALL this%Warning("InitGravity_binary","Rotating mesh enabled but OMEGA set to zero in mesh")
216 ELSE
217 this%r0(:) = 0.0
218 END IF
219
220 ! reset mass flux
221 this%mdot = 0.0
222 this%massloss = 0.0
223 this%accrate = 0.0
224
225 ! set period
226 this%period = 2.*pi*sqrt(this%semaaxis/(this%mass+this%mass2)/physics%constants%GN)*this%semaaxis
227
228 !initialise output
229 CALL this%SetOutput(mesh,physics,config,io)
230 CALL this%InfoGravity()
231
232 END SUBROUTINE initgravity_binary
233
234
236 SUBROUTINE infogravity(this)
237 IMPLICIT NONE
238 !------------------------------------------------------------------------!
239 CLASS(gravity_binary), INTENT(IN) :: this
240 !------------------------------------------------------------------------!
241 CHARACTER(LEN=32) :: mass1_str,mass2_str,excent_str,sema_str,omega_str
242 !------------------------------------------------------------------------!
243 WRITE (mass1_str,'(ES10.2)') this%mass
244 WRITE (mass2_str,'(ES10.2)') this%mass2
245 WRITE (excent_str,'(ES10.2)') this%excent
246 WRITE (sema_str,'(ES10.2)') this%semaaxis
247 CALL this%Info(" primary mass: " // trim(mass1_str))
248 CALL this%Info(" secondary mass: " // trim(mass2_str))
249 CALL this%Info(" excentricity: " // trim(excent_str))
250 CALL this%Info(" semi major axis: " // trim(sema_str))
251 IF(this%period.NE.0.0) THEN
252 WRITE(omega_str,'(ES10.2)') 2.0*pi/this%period
253 CALL this%Info(" pattern speed: " // trim(omega_str))
254 END IF
255 ! some sanity checks
256 IF (this%excent.LT.0.OR.this%excent.GE.1) &
257 CALL this%Error("InitGravity_binary","excentricity out of range, should be 0 <= e < 1 => check setup")
258 IF (this%semaaxis.LE.0) &
259 CALL this%Error("InitGravity_binary","semi major axis out of range, should be 0 < a => check setup")
260 END SUBROUTINE infogravity
261
262
263 SUBROUTINE setoutput(this,Mesh,Physics,config,IO)
264 IMPLICIT NONE
265 !------------------------------------------------------------------------!
266 CLASS(gravity_binary), INTENT(INOUT) :: this
267 CLASS(mesh_base), INTENT(IN) :: Mesh
268 CLASS(physics_base), INTENT(IN) :: Physics
269 TYPE(dict_typ), POINTER :: config,IO
270 !------------------------------------------------------------------------!
271 INTEGER :: valwrite
272 !------------------------------------------------------------------------!
273
274 CALL this%gravity_pointmass%SetOutput(mesh,physics,config,io)
275
276 ! output cartesian positions of both components of the binary system
277 CALL getattr(config, "output/binpos", valwrite, 0)
278 IF (valwrite .EQ. 1) &
279 CALL setattr(io, "binpos", this%pos)
280
281 END SUBROUTINE setoutput
282
283
293 SUBROUTINE updategravity_single(this,Mesh,Physics,Fluxes,pvar,time,dt)
294 IMPLICIT NONE
295 !------------------------------------------------------------------------!
296 CLASS(gravity_binary), INTENT(INOUT) :: this
297 CLASS(mesh_base), INTENT(IN) :: Mesh
298 CLASS(physics_base), INTENT(IN) :: Physics
299 CLASS(fluxes_base), INTENT(IN) :: Fluxes
300 CLASS(marray_compound),INTENT(INOUT) :: pvar
301 REAL, INTENT(IN) :: time,dt
302 !------------------------------------------------------------------------!
303 REAL :: GM1, GM2
304 INTEGER :: i,j,k
305 !------------------------------------------------------------------------!
306 ! update the positions of the binary system
307 CALL this%UpdatePositions(time)
308
309 ! compute curvilinear components of the position vectors
310 ! 1. primary component
311 this%posvec_prim_tmp(:,:,:,1) = this%pos(1,1)
312 this%posvec_prim_tmp(:,:,:,2) = this%pos(2,1)
313 this%posvec_prim_tmp(:,:,:,3) = this%pos(3,1)
314 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,east,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,1,:))
315 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,north,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,2,:))
316 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,top,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,3,:))
317 CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,this%posvec_prim_tmp,this%posvec_prim)
318 ! 2. secondary component
319 this%posvec_sec_tmp(:,:,:,1) = this%pos(1,2)
320 this%posvec_sec_tmp(:,:,:,2) = this%pos(2,2)
321 this%posvec_sec_tmp(:,:,:,3) = this%pos(3,2)
322 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,east,:),this%posvec_sec_tmp(:,:,:,:),this%fposvec_sec(:,:,:,1,:))
323 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,north,:),this%posvec_sec_tmp(:,:,:,:),this%fposvec_sec(:,:,:,2,:))
324 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,top,:),this%posvec_sec_tmp(:,:,:,:),this%fposvec_sec(:,:,:,3,:))
325 CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,this%posvec_sec_tmp,this%posvec_sec)
326
327 ! subtract the result from the position vectors of all cell bary centers
328 ! this gives you the curvilinear components of all vectors pointing
329 ! from the point mass of each component of the binary system to the bary
330 ! centers of all cells on the mesh
331 this%posvec_prim(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - this%posvec_prim(:,:,:,:)
332 this%fposvec_prim(:,:,:,1,:) = mesh%posvec%faces(:,:,:,east,:) - this%fposvec_prim(:,:,:,1,:) ! EAST
333 this%fposvec_prim(:,:,:,2,:) = mesh%posvec%faces(:,:,:,north,:) - this%fposvec_prim(:,:,:,2,:) ! NORTH
334 this%fposvec_prim(:,:,:,3,:) = mesh%posvec%faces(:,:,:,top,:) - this%fposvec_prim(:,:,:,3,:) ! TOP
335
336 this%posvec_sec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - this%posvec_sec(:,:,:,:)
337 this%fposvec_sec(:,:,:,1,:) = mesh%posvec%faces(:,:,:,east,:) - this%fposvec_sec(:,:,:,1,:) ! EAST
338 this%fposvec_sec(:,:,:,2,:) = mesh%posvec%faces(:,:,:,north,:) - this%fposvec_sec(:,:,:,2,:) ! NORTH
339 this%fposvec_sec(:,:,:,3,:) = mesh%posvec%faces(:,:,:,top,:) - this%fposvec_sec(:,:,:,3,:) ! TOP
340
341 ! compute the distances between each component of the binary system
342 ! and all cell bary centers
343 this%r_prim(:,:,:) = sqrt(this%posvec_prim(:,:,:,1)**2 &
344 + this%posvec_prim(:,:,:,2)**2 + this%posvec_prim(:,:,:,3)**2)
345 this%fr_prim(:,:,:,1) = sqrt(this%fposvec_prim(:,:,:,1,1)**2 &
346 + this%fposvec_prim(:,:,:,1,2)**2 + this%fposvec_prim(:,:,:,1,3)**2) ! shifted EAST-faces
347 this%fr_prim(:,:,:,2) = sqrt(this%fposvec_prim(:,:,:,2,1)**2 &
348 + this%fposvec_prim(:,:,:,2,2)**2 + this%fposvec_prim(:,:,:,2,3)**2) ! shifted NORTH-faces
349 this%fr_prim(:,:,:,3) = sqrt(this%fposvec_prim(:,:,:,3,1)**2 &
350 + this%fposvec_prim(:,:,:,3,2)**2 + this%fposvec_prim(:,:,:,3,3)**2) ! shifted TOP-faces
351 this%r_sec(:,:,:) = sqrt(this%posvec_sec(:,:,:,1)**2 &
352 + this%posvec_sec(:,:,:,2)**2 + this%posvec_sec(:,:,:,3)**2)
353 this%fr_sec(:,:,:,1) = sqrt(this%fposvec_sec(:,:,:,1,1)**2 &
354 + this%fposvec_sec(:,:,:,1,2)**2 + this%fposvec_sec(:,:,:,1,3)**2) ! shifted EAST-faces
355 this%fr_sec(:,:,:,2) = sqrt(this%fposvec_sec(:,:,:,2,1)**2 &
356 + this%fposvec_sec(:,:,:,2,2)**2 + this%fposvec_sec(:,:,:,2,3)**2) ! shifted NORTH-faces
357 this%fr_sec(:,:,:,3) = sqrt(this%fposvec_sec(:,:,:,3,1)**2 &
358 + this%fposvec_sec(:,:,:,3,2)**2 + this%fposvec_sec(:,:,:,3,3)**2) ! shifted NORTH-faces
359
360 ! compute square of Keplerian velocities for updated time value
361 ! with respect to PRIMARY star
362 gm1 = physics%Constants%GN*this%GetMass_primary(time)
363 this%omega2%data4d(:,:,:,1) = gm1 / (this%r_prim(:,:,:)**3 + this%eps1**3)
364
365 ! and SECONDARY star
366 gm2 = physics%Constants%GN*this%GetMass_secondary(time)
367 this%omega2%data4d(:,:,:,2) = gm2 / (this%r_sec(:,:,:)**3 + this%eps2**3)
368
369 ! potential calculation for binary and use pointmassroutine for it
370 CALL this%gravity_pointmass%CalcPotential(mesh,physics,this%mass,this%r_prim,this%fr_prim,this%pot_prim%data4d)
371 CALL this%gravity_pointmass%CalcPotential(mesh,physics,this%mass2,this%r_sec,this%fr_sec,this%pot_sec%data4d)
372
373 this%pot%data1d(:) = this%pot_prim%data1d(:) + this%pot_sec%data1d(:)
374
375 ! set curvilinear components of the gravitational acceleration
376!NEC$ collapse
377 DO k=mesh%KGMIN,mesh%KGMAX
378 DO j=mesh%JGMIN,mesh%JGMAX
379!NEC$ ivdep
380 DO i=mesh%IGMIN,mesh%IGMAX
381 this%accel%data4d(i,j,k,1:physics%VDIM) = -this%omega2%data4d(i,j,k,1) * this%posvec_prim(i,j,k,1:physics%VDIM)&
382 -this%omega2%data4d(i,j,k,2) * this%posvec_sec(i,j,k,1:physics%VDIM)
383 END DO
384 END DO
385 END DO
386
387 END SUBROUTINE updategravity_single
388
389
400 SUBROUTINE updatepositions(this,time)
401 IMPLICIT NONE
402 !------------------------------------------------------------------------!
403 CLASS(gravity_binary), INTENT(INOUT) :: this
404 REAL, INTENT(IN) :: time
405 !------------------------------------------------------------------------!
406 REAL :: r,phi,r1
407 INTEGER :: err
408 !------------------------------------------------------------------------!
409 CALL kepler(time,this%period,this%excent,this%semaaxis,r,phi)
410 IF(r.LT.0) CALL this%Error("UpdatePositions_binary","Solving Kepler-Equation failed!")
411
412 !transform to rotating reference frame if necessary
413 phi = phi - this%omega_rot * time
414
415 ! cartesian coordinates of SECONDARY component (right at t=0)
416 ! \todo Generalize this implementation to arbitrary orbit orientations;
417 ! it works currently only if the angular momentum axis points along the z-axis.
418 r1 = this%mass/(this%mass+this%mass2)*r
419 this%pos(1,2) = r1*cos(phi)
420 this%pos(2,2) = r1*sin(phi)
421 this%pos(3,2) = 0.0
422 ! and PRIMARY COMPONENT
423 this%pos(:,1) = -this%pos(:,2) * this%mass2/this%mass
424 ! shift with respect to center of rotation
425 this%pos(:,1) = this%pos(:,1) + this%r0(:)
426 this%pos(:,2) = this%pos(:,2) + this%r0(:)
427 END SUBROUTINE updatepositions
428
440 PURE SUBROUTINE calcdiskheight_single(this,Mesh,Physics,pvar,bccsound,h_ext,height)
441 IMPLICIT NONE
442 !------------------------------------------------------------------------!
443 CLASS(gravity_binary), INTENT(INOUT) :: this
444 CLASS(mesh_base), INTENT(IN) :: mesh
445 CLASS(physics_base), INTENT(IN) :: physics
446 CLASS(marray_compound), INTENT(INOUT) :: pvar
447 TYPE(marray_base), INTENT(INOUT) :: bccsound,h_ext,height
448 !------------------------------------------------------------------------!
449 ! compute the effective omega
450 this%omega%data1d(:) = sqrt(this%omega2%data2d(:,1) + this%omega2%data2d(:,2))
451 ! compute the disk height
452 height%data1d(:) = getdiskheight(bccsound%data1d(:),this%omega%data1d(:))
453 END SUBROUTINE calcdiskheight_single
454
455
457 FUNCTION getmass_primary(this,time) RESULT(mass)
458 IMPLICIT NONE
459 !------------------------------------------------------------------------!
460 CLASS(gravity_binary), INTENT(IN) :: this
461 REAL, INTENT(IN) :: time
462 !------------------------------------------------------------------------!
463 REAL :: mass
464 !------------------------------------------------------------------------!
465 IF(time.LE.this%switchon) THEN
466 mass = this%mass * sin(0.5*pi*time/this%switchon)**2
467 ELSE
468 mass = this%mass
469 END IF
470 END FUNCTION getmass_primary
471
473 FUNCTION getmass_secondary(this,time) RESULT(mass)
474 IMPLICIT NONE
475 !------------------------------------------------------------------------!
476 CLASS(gravity_binary), INTENT(IN) :: this
477 REAL, INTENT(IN) :: time
478 !------------------------------------------------------------------------!
479 REAL :: mass
480 !------------------------------------------------------------------------!
481 IF(time.LE.this%switchon2) THEN
482 mass = this%mass2 * sin(0.5*pi*time/this%switchon2)**2
483 ELSE
484 mass = this%mass2
485 END IF
486 END FUNCTION getmass_secondary
487
489 SUBROUTINE finalize(this)
490 IMPLICIT NONE
491 !------------------------------------------------------------------------!
492 TYPE(gravity_binary), INTENT(INOUT) :: this
493 !------------------------------------------------------------------------!
494 DEALLOCATE(this%r0,this%r_sec,this%posvec_sec,this%posvec_sec_tmp,&
495 this%pot_prim,this%pot_sec,this%fr_sec,this%fposvec_sec,this%mass2)
496
497 CALL this%Finalize_base()
498 END SUBROUTINE finalize
499
556 PURE SUBROUTINE kepler(time,period,excent,semax,r,phi)
557 USE roots
558 IMPLICIT NONE
559 !------------------------------------------------------------------------!
560 REAL, INTENT(IN) :: time
561 REAL, INTENT(IN) :: period
562 REAL, INTENT(IN) :: excent
563 REAL, INTENT(IN) :: semax
564 REAL, INTENT(OUT):: r
565 REAL, INTENT(OUT):: phi
566 !------------------------------------------------------------------------!
567 REAL, DIMENSION(2):: plist
568 REAL :: tau,e,cose
569 INTEGER :: err
570 !------------------------------------------------------------------------!
571 ! mean anomaly
572 tau = 2.*pi*modulo(time,period)/period
573 plist(1)=excent
574 plist(2)=tau
575 IF ((tau.GE.0.0) .AND. (tau.LE.pi)) THEN
576 ! intersection lies ABOVE abscissa
577 CALL getroot(func,max(0.0,tau),min(pi,excent+tau),e,err,plist)
578 ELSE
579 ! intersection lies BELOW abscissa
580 CALL getroot(func,max(pi,tau-excent),min(2.*pi,tau),e,err,plist)
581 END IF
582 ! return a negative number for distance r
583 IF (err.NE.0) THEN
584 r = -1
585 phi = 0
586 RETURN
587 END IF
588 ! compute new distance and position angle
589 cose = cos(e)
590 ! distance to the origin of primary star position
591 r = semax * (1.0 - excent*cose)
592 ! position angle of primary star
593 ! phi1 = 2.0*ATAN(SQRT((1.0+excent)/(1.0-excent))*TAN(0.5*E))
594 ! A little bit more complicated but this avoids one evaluation of TAN;
595 ! uses TAN(x/2) = +/- SQRT((1-cos(x))/(1+cos(x)))
596 IF (cose.NE.-1.0) THEN
597 phi = 2.0*atan(sign(sqrt(((1.+excent)*(1.-cose)) &
598 /((1.-excent)*(1.+cose))),pi-e))
599 ELSE
600 ! since lim_{x->0} atan(q/x) = pi/2 for q>0
601 phi = pi
602 END IF
603 END SUBROUTINE kepler
604
606 PURE SUBROUTINE func(x,fx,plist)
607 IMPLICIT NONE
608 !------------------------------------------------------------------------!
609 REAL, INTENT(IN) :: x
610 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
611 !plist[1]=excent, plist[2]=tau
612 REAL, INTENT(OUT) :: fx
613 !------------------------------------------------------------------------!
614 fx = x - plist(1)*sin(x) - plist(2)
615 END SUBROUTINE func
616
617END MODULE gravity_binary_mod
Dictionary for generic data types.
Definition: common_dict.f90:61
type(logging_base), save this
base module for numerical flux functions
Definition: fluxes_base.f90:39
generic gravity terms module providing functionaly common to all gravity terms
source terms module for gravitational acceleration due to two pointmasses
subroutine updategravity_single(this, Mesh, Physics, Fluxes, pvar, time, dt)
calculates the resultant gravitational accerleration of a binary system
real function getmass_primary(this, time)
get mass of the primary component, eventually scaled by switch on factor
subroutine infogravity(this)
write binary parameters to screen
real function getmass_secondary(this, time)
get mass of the primary component, eventually scaled by switch on factor
subroutine updatepositions(this, time)
compute new positions for both components of the binary system
pure subroutine func(x, fx, plist)
find root of the Kepler equation
pure subroutine calcdiskheight_single(this, Mesh, Physics, pvar, bccsound, h_ext, height)
compute the pressure scale height of a disk for the binary potential
subroutine initgravity_binary(this, Mesh, Physics, config, IO)
Constructor of gravity binary module.
pure subroutine, public kepler(time, period, excent, semax, r, phi)
solve Kepler-Equation and return new distance r and position angle phi
source terms module for gravitational acceleration due to a point mass at the center of the coordinat...
elemental real function, public getdiskheight(cs, omega)
return pressure scale height of a geometrically thin Keplerian disk
Basic fosite module.
subroutine finalize(this)
Destructor of logging_base class.
base class for mesh arrays
Definition: marray_base.f90:36
derived class for compound of mesh arrays
basic mesh module
Definition: mesh_base.f90:72
integer, parameter east
named constant for eastern boundary
Definition: mesh_base.f90:101
subroutine setoutput(this, config, IO)
Setup mesh fields for i/o.
Definition: mesh_base.f90:828
integer, parameter top
named constant for top boundary
Definition: mesh_base.f90:101
integer, parameter north
named constant for northern boundary
Definition: mesh_base.f90:101
Basic physics module.
physics module for 1D,2D and 3D non-isothermal Euler equations
physics module for 1D,2D and 3D isothermal Euler equations
root finding subroutines
Definition: roots.f90:45
basic mesh array class
Definition: marray_base.f90:69
mesh data structure
Definition: mesh_base.f90:122