gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Implement Keplerian orbit for the satellite in binaries problem.

Browse Source 
- now the satellite changes positions and velocity accordind to Keplers
   laws; correct Keplerian orbit has been implemented in
   shape_binaries();

 - correct initialization of the satellite in init_binaries();
This commit is contained in:
Grzegorz Kowal 2011-05-17 20:23:44 -03:00
parent 68a39a28e4
commit e68a189916
1 changed files with 119 additions and 61 deletions

114
src/problem.F90

@ -709,7 +709,7 @@ module problem
use config , only : im, jm, km
use config , only : dens, pres, csnd2, gamma
use config , only : dnfac, dnrat
use config , only : rstar, vstar, rsat, dsat, vsat, tsat
use config , only : rstar, vstar, rsat, dsat, vsat, tsat, esat
use constants, only : dpi
use coords , only : ax, ay, az, adr
use scheme , only : prim2cons
@ -732,13 +732,13 @@ module problem
! local variables
!
integer :: i, j, k
real :: dr
real :: dnamb, pramb
real :: dnstar, prstar, vrstar, rc
real :: dnsat , prsat , vrsat , rs, xs, ys
real :: dr, rc, rs, xs, ys
real :: om = 0.0d0, sn = 0.0d0, cs = 1.0d0
real :: asat, bsat, xsat
real :: xsh, ysh, xvl, yvl
real, save :: dnamb, pramb
real, save :: dnstar, prstar, vrstar
real, save :: dnsat , prsat , vrsat
real, save :: asat, bsat, xsat = - 1.0d0
real, save :: xsh, ysh, xvl, yvl
! local arrays
!
@ -749,6 +749,10 @@ module problem
!
!-------------------------------------------------------------------------------
!
! calculate the interior and orbit parameters only once
!
if (xsat .le. 0.0d0) then
! calculate pressure from sound speed
!
#ifdef ISO
@ -758,7 +762,7 @@ module problem
pres = csnd2 * dens / gamma
#endif /* ADI */
! calculate parameters
! calculate star interior parameters
!
dnamb = dens
pramb = pres
@ -768,13 +772,27 @@ module problem
prsat = pramb * dnfac / dnrat
vrstar = vstar / rstar
vrsat = vsat / rsat
if (tsat .gt. 0.0d0) then
! calculate orbit parameters
!
asat = dsat / (1.0d0 - esat)
xsat = asat * esat
bsat = sqrt(asat * asat - xsat * xsat)
! angular speed and triginometric functions
!
om = dpi / tsat
sn = 0.0d0
cs = 1.0d0
! calculate satellite initial position and speed
!
xsh = asat * cs - xsat
ysh = bsat * sn
xvl = - asat * sn * om
yvl = bsat * cs * om
end if
xvl = - asat * om * sn
yvl = bsat * om * cs
! obtain the diagonal size of a cell at the current level
!
@ -1386,19 +1404,22 @@ module problem
! local variables
!
integer :: i, j, k
real :: dr
real :: dnamb, pramb
real :: dnstar, vrstar, prstar, rc
real :: dnsat , vrsat , prsat , rs, xs, ys
real :: dr, rs, xs, ys, rc
real :: om = 0.0d0, sn = 0.0d0, cs = 1.0d0
real :: rhs, res
real, save :: dnamb, pramb
real, save :: dnstar, vrstar, prstar
real, save :: dnsat , vrsat , prsat
#ifdef ADI
real :: ekin, ekstar, eksat
real, save :: ekstar, eksat
real :: ekin
#ifdef MHD
real :: emag
#endif /* MHD */
#endif /* ADI */
real :: om = 0.0d0, ang = 0.0d0, sn = 0.0d0, cs = 1.0d0
real :: asat, bsat, xsat
real :: xsh, ysh, xvl, yvl
real, save :: asat, bsat, xsat = -1.0d0
real, save :: xsh, ysh, xvl, yvl
real, save :: tpre = -1.0d0, ang = 0.0d0
! local arrays
!
@ -1408,7 +1429,11 @@ module problem
!
!-------------------------------------------------------------------------------
!
! calculate parameters
! calculate the star and satellite orbit parameters only once
!
if (xsat .lt. 0.0d0) then
! calculate parameters for the star interiors
!
#ifdef ISO
pres = csnd2 * dens
@ -1428,19 +1453,52 @@ module problem
ekstar = 0.5d0 * dnstar * vstar * vstar
eksat = 0.5d0 * dnsat * vsat * vsat
#endif /* ADI */
if (tsat .gt. 0.0d0) then
om = dpi / tsat
ang = om * t
sn = sin(ang)
cs = cos(ang)
end if
! calculate orbit parameters
!
asat = dsat / (1.0d0 - esat)
xsat = asat * esat
bsat = sqrt(asat * asat - xsat * xsat)
end if
! recalculate the angles and coordinates
!
if (t .gt. tpre) then
! calculate the new value of eccentric anomaly
!
rhs = dpi * t / tsat
res = 1.0d0
i = 1
do while(res .ge. 1.0e-14 .and. i .le. 1000)
sn = esat * sin(ang)
ang = rhs + sn
res = abs((ang - sn) - rhs)
i = i + 1
end do
! calculate trigonometric functions
!
sn = sin(ang)
cs = cos(ang)
! calculate the angular velocity
!
om = dpi / tsat / (1.0d0 - esat * cs)
! calculate position coordinates and velocity components
!
xsh = asat * cs - xsat
ysh = bsat * sn
xvl = - asat * om * sn
yvl = bsat * om * cs
xvl = - asat * sn * om
yvl = bsat * cs * om
! substitute the previous time
!
tpre = t
end if
! obtain the diagonal size of a cell at the current level
!