Fixed incorrect calculation of angular momentum error. The error should be calculated from the difference in the ang. mtm. vector quantities not the magnitudes.

MintonGroup · May 10, 2021 · 4879d19 · 4879d19
1 parent b88dc07
commit 4879d19
Show file tree

Hide file tree

Showing 2 changed files with 4 additions and 4 deletions.
diff --git a/src/io/io_conservation_report.f90 b/src/io/io_conservation_report.f90
@@ -56,7 +56,7 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
          flush(egyiu)
          if (.not.lfirst .and. lterminal) then 
             Lmag_now = norm2(Ltot_now)
-            Lerror = (Lmag_now - Lmag_orig) / Lmag_orig
+            Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig
             Eorbit_error = (Eorbit - Eorbit_orig) / abs(Eorbit_orig)
             Ecoll_error = -Ecollisions / abs(Eorbit_orig)
             Etotal_error = (Eorbit - (Eorbit_orig - Ecollisions)) / abs(Eorbit_orig)

diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90
@@ -104,7 +104,7 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad
                                          (ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), &
                                          (pe_after - pe_before) / abs(Etot_before), &
                                          (Etot_after - Etot_before) / abs(Etot_before), &
-                                         (Lmag_after - Lmag_before) / Lmag_before
+                                         norm2(Ltot_after - Ltot_before) / Lmag_before
       write(*,        "(' ---------------------------------------------------------------------------')")
       write(*,        "('  Second pass to get energy ')")
       write(*,        "(' ---------------------------------------------------------------------------')")
@@ -138,7 +138,7 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad
                                         (ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), &
                                         (pe_after - pe_before) / abs(Etot_before), &
                                         (Etot_after - Etot_before) / abs(Etot_before), &
-                                        (Lmag_after - Lmag_before) / Lmag_before
+                                        norm2(Ltot_after - Ltot_before) / Lmag_before
 
       lmerge = lmerge .or. ((Etot_after - Etot_before) / abs(Etot_before) > 0._DP) 
       if (.not.lmerge) then
@@ -163,7 +163,7 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad
                                        (ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), &
                                        (pe_after - pe_before) / abs(Etot_before), &
                                        (Etot_after - Etot_before) / abs(Etot_before), &
-                                       (Lmag_after - Lmag_before) / Lmag_before
+                                       norm2(Ltot_after - Ltot_before) / Lmag_before
       write(*,        "(' ---------------------------------------------------------------------------')")
       !****************************************************************************************************************