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!----------------------------------------------------------------------------!
53MODULE gravity_binary_mod
54 USE gravity_pointmass_mod
55 USE gravity_base_mod
56 USE fluxes_base_mod
57 USE physics_base_mod
58 USE mesh_base_mod
59 USE logging_base_mod
60 USE marray_base_mod
61 USE marray_compound_mod
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
76 TYPE, EXTENDS(gravity_pointmass) :: gravity_binary
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
92 PROCEDURE :: updategravity_single
93 PROCEDURE :: updatepositions
94 PROCEDURE :: infogravity
95 PROCEDURE :: setoutput
96 PROCEDURE :: getmass_primary
97 PROCEDURE :: getmass_secondary
98 PROCEDURE :: calcdiskheight_single
99 final :: finalize
100 END TYPE
101
102 PUBLIC :: &
103 ! types
104 gravity_binary, &
105 ! module procedures
106 kepler
107 !--------------------------------------------------------------------------!
108
109CONTAINS
115 SUBROUTINE initgravity_binary(this,Mesh,Physics,config,IO)
116 USE physics_euler_mod, ONLY: physics_euler
117 USE physics_eulerisotherm_mod, ONLY: physics_eulerisotherm
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
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
common_dict::this
type(logging_base), save this
Definition: common_dict.f90:149
fluxes_base_mod
base module for numerical flux functions
Definition: fluxes_base.f90:39
gravity_base_mod
generic gravity terms module providing functionaly common to all gravity terms
Definition: gravity_base.f90:46
gravity_binary_mod
source terms module for gravitational acceleration due to two pointmasses
Definition: gravity_binary.f90:53
gravity_binary_mod::updategravity_single
subroutine updategravity_single(this, Mesh, Physics, Fluxes, pvar, time, dt)
calculates the resultant gravitational accerleration of a binary system
Definition: gravity_binary.f90:294
gravity_binary_mod::getmass_primary
real function getmass_primary(this, time)
get mass of the primary component, eventually scaled by switch on factor
Definition: gravity_binary.f90:458
gravity_binary_mod::infogravity
subroutine infogravity(this)
write binary parameters to screen
Definition: gravity_binary.f90:237
gravity_binary_mod::getmass_secondary
real function getmass_secondary(this, time)
get mass of the primary component, eventually scaled by switch on factor
Definition: gravity_binary.f90:474
gravity_binary_mod::updatepositions
subroutine updatepositions(this, time)
compute new positions for both components of the binary system
Definition: gravity_binary.f90:401
gravity_binary_mod::func
pure subroutine func(x, fx, plist)
find root of the Kepler equation
Definition: gravity_binary.f90:607
gravity_binary_mod::calcdiskheight_single
pure subroutine calcdiskheight_single(this, Mesh, Physics, pvar, bccsound, h_ext, height)
compute the pressure scale height of a disk for the binary potential
Definition: gravity_binary.f90:441
gravity_binary_mod::initgravity_binary
subroutine initgravity_binary(this, Mesh, Physics, config, IO)
Constructor of gravity binary module.
Definition: gravity_binary.f90:116
gravity_binary_mod::kepler
pure subroutine, public kepler(time, period, excent, semax, r, phi)
solve Kepler-Equation and return new distance r and position angle phi
Definition: gravity_binary.f90:557
gravity_pointmass_mod
source terms module for gravitational acceleration due to a point mass at the center of the coordinat...
Definition: gravity_pointmass.f90:48
gravity_pointmass_mod::getdiskheight
elemental real function, public getdiskheight(cs, omega)
return pressure scale height of a geometrically thin Keplerian disk
Definition: gravity_pointmass.f90:505
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
logging_base_mod::finalize
subroutine finalize(this)
Destructor of logging_base class.
Definition: logging_base.f90:160
marray_base_mod
base class for mesh arrays
Definition: marray_base.f90:36
marray_compound_mod
derived class for compound of mesh arrays
Definition: marray_compound.f90:35
mesh_base_mod
basic mesh module
Definition: mesh_base.f90:72
mesh_base_mod::east
integer, parameter east
named constant for eastern boundary
Definition: mesh_base.f90:101
mesh_base_mod::setoutput
subroutine setoutput(this, config, IO)
Setup mesh fields for i/o.
Definition: mesh_base.f90:828
mesh_base_mod::top
integer, parameter top
named constant for top boundary
Definition: mesh_base.f90:101
mesh_base_mod::north
integer, parameter north
named constant for northern boundary
Definition: mesh_base.f90:101
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
physics_euler_mod
physics module for 1D,2D and 3D non-isothermal Euler equations
Definition: physics_euler.f90:41
physics_eulerisotherm_mod
physics module for 1D,2D and 3D isothermal Euler equations
Definition: physics_eulerisotherm.f90:44
roots
root finding subroutines
Definition: roots.f90:45
common_dict::dict_typ
Definition: common_dict.f90:123
fluxes_base_mod::fluxes_base
Definition: fluxes_base.f90:61
gravity_binary_mod::gravity_binary
Definition: gravity_binary.f90:76
gravity_pointmass_mod::gravity_pointmass
Definition: gravity_pointmass.f90:77
marray_base_mod::marray_base
basic mesh array class
Definition: marray_base.f90:69
marray_compound_mod::marray_compound
data types and methods
Definition: marray_compound.f90:46
mesh_base_mod::mesh_base
mesh data structure
Definition: mesh_base.f90:122
physics_base_mod::physics_base
Definition: physics_base.f90:69
physics_euler_mod::physics_euler
Definition: physics_euler.f90:54
physics_eulerisotherm_mod::physics_eulerisotherm
Definition: physics_eulerisotherm.f90:56
Generated by   doxygen 1.9.4