Merge branch 'master' into boundaries

src/evolution.F90

@@ -249,6 +249,8 @@ module evolution
 
       end do ! meta blocks
 
+    end if ! maxlev > 1
+
 ! find the smallest spacial step
 !
 #if NDIMS == 2
@@ -258,8 +260,6 @@ module evolution
     dxmin = min(adx(lev), ady(lev), adz(lev))
 #endif /* NDIMS == 3 */
 
-    end if ! maxlev > 1
-
 ! iterate over all data blocks in order to find the maximum speed among them
 !
     pdata => list_data

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       :: dr, rc, rs, xs, ys
-    real    :: dnamb, pramb
-    real    :: dnstar, prstar, vrstar, rc
-    real    :: dnsat , prsat , vrsat , rs, xs, ys
     real       :: om = 0.0d0, sn = 0.0d0, cs = 1.0d0
-    real    :: asat, bsat, xsat
+    real, save :: dnamb, pramb
-    real    :: xsh, ysh, xvl, yvl
+    real, save :: dnstar, prstar, vrstar
+    real, save :: dnsat , prsat , vrsat
+    real, save :: asat, bsat = - 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 (bsat .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,26 @@ module problem
       prsat  = pramb * dnfac / dnrat
       vrstar = vstar / rstar
       vrsat  = vsat  / rsat
-    if (tsat .gt. 0.0d0) then
+
-      om      = dpi / tsat
+! calculate orbit parameters
+!
+      asat   = dsat / (1.0d0 - esat)
+      bsat   = asat * sqrt(1.0d0 - esat * esat)
+
+! angular speed and triginometric functions
+!
       sn     = 0.0d0
       cs     = 1.0d0
+      om     = dpi / tsat / (1.0d0 - esat * cs)
+
+! calculate satellite initial position and speed
+!
+      xsh    =   asat * (cs - esat)
+      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 +1403,22 @@ module problem
 ! local variables
 !
     integer    :: i, j, k
-    real    :: dr
+    real       :: dr, rs, xs, ys, rc
-    real    :: dnamb, pramb
+    real       :: om = 0.0d0, sn = 0.0d0, cs = 1.0d0
-    real    :: dnstar, vrstar, prstar, rc
+    real       :: rhs, res
-    real    :: dnsat , vrsat , prsat , rs, xs, ys
+    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, save :: asat, bsat = -1.0d0
-    real    :: asat, bsat, xsat
+    real, save :: xsh, ysh, xvl, yvl
-    real    :: xsh, ysh, xvl, yvl
+    real, save :: tpre = -1.0d0, ang = 0.0d0
 
 ! local arrays
 !
@@ -1408,7 +1428,11 @@ module problem
 !
 !-------------------------------------------------------------------------------
 !
-! calculate parameters
+! calculate the star and satellite orbit parameters only once
+!
+    if (bsat .lt. 0.0d0) then
+
+! calculate parameters for the star interiors
 !
 #ifdef ISO
       pres   = csnd2 * dens
@@ -1428,19 +1452,51 @@ 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
+! calculate orbit parameters
-      ang     = om * t
+!
+      asat   = dsat / (1.0d0 - esat)
+      bsat   = asat * sqrt(1.0d0 - esat * esat)
+
+    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-15 .and. i .le. 100)
+        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)
-    end if
+
-    asat    = dsat / (1.0d0 - esat)
+! calculate the angular velocity
-    xsat    = asat * esat
+!
-    bsat    = sqrt(asat * asat - xsat * xsat)
+      om = dpi / tsat / (1.0d0 - esat * cs)
-    xsh     = asat * cs - xsat
+
+! calculate position coordinates and velocity components
+!
+      xsh  =   asat * (cs - esat)
       ysh  =   bsat * sn
-    xvl     = - asat * om * sn
+      xvl  = - asat * sn * om
-    yvl     =   bsat * om * cs
+      yvl  =   bsat * cs * om
+
+! substitute the previous time
+!
+      tpre = t
+
+    end if
 
 ! obtain the diagonal size of a cell at the current level
 !