diff --git a/src/multipole.h b/src/multipole.h index 950ec991dad9d94a06fcd136100c7245f9750e97..05e9916c0f5d864e8b9c4cd9400348bc1d78d04f 100644 --- a/src/multipole.h +++ b/src/multipole.h @@ -1295,249 +1295,284 @@ __attribute__((nonnull)) INLINE static void gravity_M2M( #endif #if SELF_GRAVITY_MULTIPOLE_ORDER > 1 - /* Shift 2nd order term (Simplify 1st order terms that are all 0) */ - m_a->M_200 = m_b->M_200 + /*X_100(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_000; - m_a->M_020 = m_b->M_020 + /*X_010(dx) * m_b->M_010 +*/ X_020(dx) * m_b->M_000; - m_a->M_002 = m_b->M_002 + /*X_001(dx) * m_b->M_001 +*/ X_002(dx) * m_b->M_000; - m_a->M_110 = m_b->M_110 + - /*X_100(dx) * m_b->M_010 + X_010(dx) * m_b->M_100 +*/ X_110(dx) * - m_b->M_000; - m_a->M_101 = m_b->M_101 + - /*X_100(dx) * m_b->M_001 + X_001(dx) * m_b->M_100 +*/ X_101(dx) * - m_b->M_000; - m_a->M_011 = m_b->M_011 + - /*X_010(dx) * m_b->M_001 + X_001(dx) * m_b->M_010 +*/ X_011(dx) * - m_b->M_000; + /* Shift 2nd order terms (1st order mpole (all 0) commented out) */ + m_a->M_002 = + m_b->M_002 /* + X_001(dx) * m_b->M_001 */ + X_002(dx) * m_b->M_000; + m_a->M_011 = + m_b->M_011 /* + X_001(dx) * m_b->M_010 */ /* + X_010(dx) * m_b->M_001 */ + + X_011(dx) * m_b->M_000; + m_a->M_020 = + m_b->M_020 /* + X_010(dx) * m_b->M_010 */ + X_020(dx) * m_b->M_000; + m_a->M_101 = + m_b->M_101 /* + X_001(dx) * m_b->M_100 */ /* + X_100(dx) * m_b->M_001 */ + + X_101(dx) * m_b->M_000; + m_a->M_110 = + m_b->M_110 /* + X_010(dx) * m_b->M_100 */ /* + X_100(dx) * m_b->M_010 */ + + X_110(dx) * m_b->M_000; + m_a->M_200 = + m_b->M_200 /* + X_100(dx) * m_b->M_100 */ + X_200(dx) * m_b->M_000; #endif + #if SELF_GRAVITY_MULTIPOLE_ORDER > 2 - /* Shift 3rd order term (Simplify 1st order terms that are all 0) */ - m_a->M_300 = m_b->M_300 + X_100(dx) * m_b->M_200 + - /*X_200(dx) * m_b->M_100 +*/ X_300(dx) * m_b->M_000; - m_a->M_030 = m_b->M_030 + X_010(dx) * m_b->M_020 + - /*X_020(dx) * m_b->M_010 +*/ X_030(dx) * m_b->M_000; - m_a->M_003 = m_b->M_003 + X_001(dx) * m_b->M_002 + - /*X_002(dx) * m_b->M_001 +*/ X_003(dx) * m_b->M_000; - m_a->M_210 = - m_b->M_210 + X_100(dx) * m_b->M_110 + X_010(dx) * m_b->M_200 + - /*X_200(dx) * m_b->M_010 + X_110(dx) * m_b->M_100 + */ X_210(dx) * - m_b->M_000; - m_a->M_201 = m_b->M_201 + X_100(dx) * m_b->M_101 + X_001(dx) * m_b->M_200 + - /*X_200(dx) * m_b->M_001 + X_101(dx) * m_b->M_100 +*/ X_201(dx) * - m_b->M_000; - m_a->M_120 = m_b->M_120 + X_010(dx) * m_b->M_110 + X_100(dx) * m_b->M_020 + - /*X_020(dx) * m_b->M_100 + X_110(dx) * m_b->M_010 +*/ X_120(dx) * - m_b->M_000; - m_a->M_021 = m_b->M_021 + X_010(dx) * m_b->M_011 + X_001(dx) * m_b->M_020 + - /*X_020(dx) * m_b->M_001 + X_011(dx) * m_b->M_010 +*/ X_021(dx) * - m_b->M_000; - m_a->M_102 = m_b->M_102 + X_001(dx) * m_b->M_101 + X_100(dx) * m_b->M_002 + - /*X_002(dx) * m_b->M_100 + X_101(dx) * m_b->M_001 +*/ X_102(dx) * - m_b->M_000; - m_a->M_012 = m_b->M_012 + X_001(dx) * m_b->M_011 + X_010(dx) * m_b->M_002 + - /*X_002(dx) * m_b->M_010 + X_011(dx) * m_b->M_001 +*/ X_012(dx) * - m_b->M_000; - m_a->M_111 = m_b->M_111 + X_100(dx) * m_b->M_011 + X_010(dx) * m_b->M_101 + - X_001(dx) * m_b->M_110 + /*X_110(dx) * m_b->M_001 + */ - /*X_101(dx) * m_b->M_010 + X_011(dx) * m_b->M_100 +*/ - X_111(dx) * m_b->M_000; + /* Shift 3rd order terms (1st order mpole (all 0) commented out) */ + m_a->M_003 = m_b->M_003 + + X_001(dx) * m_b->M_002 /* + X_002(dx) * m_b->M_001 */ + + X_003(dx) * m_b->M_000; + m_a->M_012 = m_b->M_012 + + X_001(dx) * m_b->M_011 /* + X_002(dx) * m_b->M_010 */ + + X_010(dx) * m_b->M_002 /* + X_011(dx) * m_b->M_001 */ + + X_012(dx) * m_b->M_000; + m_a->M_021 = m_b->M_021 + X_001(dx) * m_b->M_020 + + X_010(dx) * m_b->M_011 /* + X_011(dx) * m_b->M_010 */ + /* + X_020(dx) * m_b->M_001 */ + + X_021(dx) * m_b->M_000; + m_a->M_030 = m_b->M_030 + + X_010(dx) * m_b->M_020 /* + X_020(dx) * m_b->M_010 */ + + X_030(dx) * m_b->M_000; + m_a->M_102 = m_b->M_102 + + X_001(dx) * m_b->M_101 /* + X_002(dx) * m_b->M_100 */ + + X_100(dx) * m_b->M_002 /* + X_101(dx) * m_b->M_001 */ + + X_102(dx) * m_b->M_000; + m_a->M_111 = m_b->M_111 + X_001(dx) * m_b->M_110 + + X_010(dx) * m_b->M_101 /* + X_011(dx) * m_b->M_100 */ + + X_100(dx) * m_b->M_011 /* + X_101(dx) * m_b->M_010 */ + /* + X_110(dx) * m_b->M_001 */ + + X_111(dx) * m_b->M_000; + m_a->M_120 = m_b->M_120 + + X_010(dx) * m_b->M_110 /* + X_020(dx) * m_b->M_100 */ + + X_100(dx) * m_b->M_020 /* + X_110(dx) * m_b->M_010 */ + + X_120(dx) * m_b->M_000; + m_a->M_201 = m_b->M_201 + X_001(dx) * m_b->M_200 + + X_100(dx) * m_b->M_101 /* + X_101(dx) * m_b->M_100 */ + /* + X_200(dx) * m_b->M_001 */ + + X_201(dx) * m_b->M_000; + m_a->M_210 = m_b->M_210 + X_010(dx) * m_b->M_200 + + X_100(dx) * m_b->M_110 /* + X_110(dx) * m_b->M_100 */ + /* + X_200(dx) * m_b->M_010 */ + + X_210(dx) * m_b->M_000; + m_a->M_300 = m_b->M_300 + + X_100(dx) * m_b->M_200 /* + X_200(dx) * m_b->M_100 */ + + X_300(dx) * m_b->M_000; #endif #if SELF_GRAVITY_MULTIPOLE_ORDER > 3 - /* Shift 4th order terms (Simplify 1st order terms that are all 0) */ - m_a->M_004 = m_b->M_004 + X_001(dx) * m_b->M_003 + X_002(dx) * m_b->M_002 + - /*X_003(dx) * m_b->M_001 +*/ X_004(dx) * m_b->M_000; - m_a->M_013 = m_b->M_013 + X_001(dx) * m_b->M_012 + X_002(dx) * m_b->M_011 + - /*X_003(dx) * m_b->M_010 +*/ X_010(dx) * m_b->M_003 + - X_011(dx) * m_b->M_002 + /*X_012(dx) * m_b->M_001 +*/ + + /* Shift 4th order terms (1st order mpole (all 0) commented out) */ + m_a->M_004 = m_b->M_004 + X_001(dx) * m_b->M_003 + + X_002(dx) * m_b->M_002 /* + X_003(dx) * m_b->M_001 */ + + X_004(dx) * m_b->M_000; + m_a->M_013 = m_b->M_013 + X_001(dx) * m_b->M_012 + + X_002(dx) * m_b->M_011 /* + X_003(dx) * m_b->M_010 */ + + X_010(dx) * m_b->M_003 + + X_011(dx) * m_b->M_002 /* + X_012(dx) * m_b->M_001 */ + X_013(dx) * m_b->M_000; m_a->M_022 = m_b->M_022 + X_001(dx) * m_b->M_021 + X_002(dx) * m_b->M_020 + - X_010(dx) * m_b->M_012 + X_011(dx) * m_b->M_011 + - /*X_012(dx) * m_b->M_010*/ +X_020(dx) * m_b->M_002 + - /*X_021(dx) * m_b->M_001*/ +X_022(dx) * m_b->M_000; + X_010(dx) * m_b->M_012 + + X_011(dx) * m_b->M_011 /* + X_012(dx) * m_b->M_010 */ + + X_020(dx) * m_b->M_002 /* + X_021(dx) * m_b->M_001 */ + + X_022(dx) * m_b->M_000; m_a->M_031 = m_b->M_031 + X_001(dx) * m_b->M_030 + X_010(dx) * m_b->M_021 + - X_011(dx) * m_b->M_020 + X_020(dx) * m_b->M_011 + - /*X_021(dx) * m_b->M_010 + X_030(dx) * m_b->M_001 + */ - X_031(dx) * m_b->M_000; - m_a->M_040 = m_b->M_040 + X_010(dx) * m_b->M_030 + X_020(dx) * m_b->M_020 + - /*X_030(dx) * m_b->M_010 +*/ X_040(dx) * m_b->M_000; - m_a->M_103 = m_b->M_103 + X_001(dx) * m_b->M_102 + X_002(dx) * m_b->M_101 + - /*X_003(dx) * m_b->M_100 +*/ X_100(dx) * m_b->M_003 + - X_101(dx) * m_b->M_002 + /*X_102(dx) * m_b->M_001 +*/ + X_011(dx) * m_b->M_020 + + X_020(dx) * m_b->M_011 /* + X_021(dx) * m_b->M_010 */ + /* + X_030(dx) * m_b->M_001 */ + + X_031(dx) * m_b->M_000; + m_a->M_040 = m_b->M_040 + X_010(dx) * m_b->M_030 + + X_020(dx) * m_b->M_020 /* + X_030(dx) * m_b->M_010 */ + + X_040(dx) * m_b->M_000; + m_a->M_103 = m_b->M_103 + X_001(dx) * m_b->M_102 + + X_002(dx) * m_b->M_101 /* + X_003(dx) * m_b->M_100 */ + + X_100(dx) * m_b->M_003 + + X_101(dx) * m_b->M_002 /* + X_102(dx) * m_b->M_001 */ + X_103(dx) * m_b->M_000; m_a->M_112 = m_b->M_112 + X_001(dx) * m_b->M_111 + X_002(dx) * m_b->M_110 + X_010(dx) * m_b->M_102 + - X_011(dx) * m_b->M_101 + /*X_012(dx) * m_b->M_100 + */ + X_011(dx) * m_b->M_101 /* + X_012(dx) * m_b->M_100 */ + X_100(dx) * m_b->M_012 + - X_101(dx) * m_b->M_011 + /*X_102(dx) * m_b->M_010 +*/ - X_110(dx) * m_b->M_002 + - /*X_111(dx) * m_b->M_001 +*/ X_112(dx) * m_b->M_000; + X_101(dx) * m_b->M_011 /* + X_102(dx) * m_b->M_010 */ + + X_110(dx) * m_b->M_002 /* + X_111(dx) * m_b->M_001 */ + + X_112(dx) * m_b->M_000; m_a->M_121 = m_b->M_121 + X_001(dx) * m_b->M_120 + X_010(dx) * m_b->M_111 + X_011(dx) * m_b->M_110 + - X_020(dx) * m_b->M_101 + /*X_021(dx) * m_b->M_100 +*/ + X_020(dx) * m_b->M_101 /* + X_021(dx) * m_b->M_100 */ + X_100(dx) * m_b->M_021 + X_101(dx) * m_b->M_020 + - X_110(dx) * m_b->M_011 + - /*X_111(dx) * m_b->M_010 + X_120(dx) * m_b->M_001 +*/ X_121(dx) * - m_b->M_000; - m_a->M_130 = m_b->M_130 + X_010(dx) * m_b->M_120 + X_020(dx) * m_b->M_110 + - /*X_030(dx) * m_b->M_100 +*/ X_100(dx) * m_b->M_030 + - X_110(dx) * m_b->M_020 + /*X_120(dx) * m_b->M_010 +*/ + X_110(dx) * m_b->M_011 /* + X_111(dx) * m_b->M_010 */ + /* + X_120(dx) * m_b->M_001 */ + + X_121(dx) * m_b->M_000; + m_a->M_130 = m_b->M_130 + X_010(dx) * m_b->M_120 + + X_020(dx) * m_b->M_110 /* + X_030(dx) * m_b->M_100 */ + + X_100(dx) * m_b->M_030 + + X_110(dx) * m_b->M_020 /* + X_120(dx) * m_b->M_010 */ + X_130(dx) * m_b->M_000; m_a->M_202 = m_b->M_202 + X_001(dx) * m_b->M_201 + X_002(dx) * m_b->M_200 + - X_100(dx) * m_b->M_102 + X_101(dx) * m_b->M_101 + - /*X_102(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_002 + - /*X_201(dx) * m_b->M_001 +*/ X_202(dx) * m_b->M_000; + X_100(dx) * m_b->M_102 + + X_101(dx) * m_b->M_101 /* + X_102(dx) * m_b->M_100 */ + + X_200(dx) * m_b->M_002 /* + X_201(dx) * m_b->M_001 */ + + X_202(dx) * m_b->M_000; m_a->M_211 = m_b->M_211 + X_001(dx) * m_b->M_210 + X_010(dx) * m_b->M_201 + X_011(dx) * m_b->M_200 + X_100(dx) * m_b->M_111 + - X_101(dx) * m_b->M_110 + X_110(dx) * m_b->M_101 + - /*X_111(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_011 + - /*X_201(dx) * m_b->M_010 + X_210(dx) * m_b->M_001 +*/ X_211(dx) * - m_b->M_000; + X_101(dx) * m_b->M_110 + + X_110(dx) * m_b->M_101 /* + X_111(dx) * m_b->M_100 */ + + X_200(dx) * m_b->M_011 /* + X_201(dx) * m_b->M_010 */ + /* + X_210(dx) * m_b->M_001 */ + + X_211(dx) * m_b->M_000; m_a->M_220 = m_b->M_220 + X_010(dx) * m_b->M_210 + X_020(dx) * m_b->M_200 + - X_100(dx) * m_b->M_120 + X_110(dx) * m_b->M_110 + - /*X_120(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_020 + - /*X_210(dx) * m_b->M_010 +*/ X_220(dx) * m_b->M_000; + X_100(dx) * m_b->M_120 + + X_110(dx) * m_b->M_110 /* + X_120(dx) * m_b->M_100 */ + + X_200(dx) * m_b->M_020 /* + X_210(dx) * m_b->M_010 */ + + X_220(dx) * m_b->M_000; m_a->M_301 = m_b->M_301 + X_001(dx) * m_b->M_300 + X_100(dx) * m_b->M_201 + - X_101(dx) * m_b->M_200 + X_200(dx) * m_b->M_101 + - /*X_201(dx) * m_b->M_100 + X_300(dx) * m_b->M_001 + */ - X_301(dx) * m_b->M_000; + X_101(dx) * m_b->M_200 + + X_200(dx) * m_b->M_101 /* + X_201(dx) * m_b->M_100 */ + /* + X_300(dx) * m_b->M_001 */ + + X_301(dx) * m_b->M_000; m_a->M_310 = m_b->M_310 + X_010(dx) * m_b->M_300 + X_100(dx) * m_b->M_210 + - X_110(dx) * m_b->M_200 + X_200(dx) * m_b->M_110 + - /*X_210(dx) * m_b->M_100 + X_300(dx) * m_b->M_010 +*/ - X_310(dx) * m_b->M_000; - m_a->M_400 = m_b->M_400 + X_100(dx) * m_b->M_300 + X_200(dx) * m_b->M_200 + - /*X_300(dx) * m_b->M_100 +*/ X_400(dx) * m_b->M_000; + X_110(dx) * m_b->M_200 + + X_200(dx) * m_b->M_110 /* + X_210(dx) * m_b->M_100 */ + /* + X_300(dx) * m_b->M_010 */ + + X_310(dx) * m_b->M_000; + m_a->M_400 = m_b->M_400 + X_100(dx) * m_b->M_300 + + X_200(dx) * m_b->M_200 /* + X_300(dx) * m_b->M_100 */ + + X_400(dx) * m_b->M_000; #endif #if SELF_GRAVITY_MULTIPOLE_ORDER > 4 - /* Shift 5th order terms */ + + /* Shift 5th order terms (1st order mpole (all 0) commented out) */ m_a->M_005 = m_b->M_005 + X_001(dx) * m_b->M_004 + X_002(dx) * m_b->M_003 + - X_003(dx) * m_b->M_002 + /*X_004(dx) * m_b->M_001 +*/ + X_003(dx) * m_b->M_002 /* + X_004(dx) * m_b->M_001 */ + X_005(dx) * m_b->M_000; m_a->M_014 = m_b->M_014 + X_001(dx) * m_b->M_013 + X_002(dx) * m_b->M_012 + - X_003(dx) * m_b->M_011 + /*X_004(dx) * m_b->M_010 +*/ + X_003(dx) * m_b->M_011 /* + X_004(dx) * m_b->M_010 */ + X_010(dx) * m_b->M_004 + X_011(dx) * m_b->M_003 + - X_012(dx) * m_b->M_002 + /*X_013(dx) * m_b->M_001 +*/ + X_012(dx) * m_b->M_002 /* + X_013(dx) * m_b->M_001 */ + X_014(dx) * m_b->M_000; m_a->M_023 = m_b->M_023 + X_001(dx) * m_b->M_022 + X_002(dx) * m_b->M_021 + X_003(dx) * m_b->M_020 + X_010(dx) * m_b->M_013 + - X_011(dx) * m_b->M_012 + X_012(dx) * m_b->M_011 + - /*X_013(dx) * m_b->M_010 +*/ X_020(dx) * m_b->M_003 + - X_021(dx) * m_b->M_002 + - /*X_022(dx) * m_b->M_001 +*/ X_023(dx) * m_b->M_000; + X_011(dx) * m_b->M_012 + + X_012(dx) * m_b->M_011 /* + X_013(dx) * m_b->M_010 */ + + X_020(dx) * m_b->M_003 + + X_021(dx) * m_b->M_002 /* + X_022(dx) * m_b->M_001 */ + + X_023(dx) * m_b->M_000; m_a->M_032 = m_b->M_032 + X_001(dx) * m_b->M_031 + X_002(dx) * m_b->M_030 + X_010(dx) * m_b->M_022 + X_011(dx) * m_b->M_021 + X_012(dx) * m_b->M_020 + X_020(dx) * m_b->M_012 + - X_021(dx) * m_b->M_011 + /*X_022(dx) * m_b->M_010 +*/ - X_030(dx) * m_b->M_002 + - /*X_031(dx) * m_b->M_001 +*/ X_032(dx) * m_b->M_000; + X_021(dx) * m_b->M_011 /* + X_022(dx) * m_b->M_010 */ + + X_030(dx) * m_b->M_002 /* + X_031(dx) * m_b->M_001 */ + + X_032(dx) * m_b->M_000; m_a->M_041 = m_b->M_041 + X_001(dx) * m_b->M_040 + X_010(dx) * m_b->M_031 + X_011(dx) * m_b->M_030 + X_020(dx) * m_b->M_021 + - X_021(dx) * m_b->M_020 + X_030(dx) * m_b->M_011 + - /*X_031(dx) * m_b->M_010 + X_040(dx) * m_b->M_001 +*/ - X_041(dx) * m_b->M_000; + X_021(dx) * m_b->M_020 + + X_030(dx) * m_b->M_011 /* + X_031(dx) * m_b->M_010 */ + /* + X_040(dx) * m_b->M_001 */ + + X_041(dx) * m_b->M_000; m_a->M_050 = m_b->M_050 + X_010(dx) * m_b->M_040 + X_020(dx) * m_b->M_030 + - X_030(dx) * m_b->M_020 + /* X_040(dx) * m_b->M_010 +*/ + X_030(dx) * m_b->M_020 /* + X_040(dx) * m_b->M_010 */ + X_050(dx) * m_b->M_000; m_a->M_104 = m_b->M_104 + X_001(dx) * m_b->M_103 + X_002(dx) * m_b->M_102 + - X_003(dx) * m_b->M_101 + /*X_004(dx) * m_b->M_100 +*/ + X_003(dx) * m_b->M_101 /* + X_004(dx) * m_b->M_100 */ + X_100(dx) * m_b->M_004 + X_101(dx) * m_b->M_003 + - X_102(dx) * m_b->M_002 + /*X_103(dx) * m_b->M_001 +*/ + X_102(dx) * m_b->M_002 /* + X_103(dx) * m_b->M_001 */ + X_104(dx) * m_b->M_000; m_a->M_113 = m_b->M_113 + X_001(dx) * m_b->M_112 + X_002(dx) * m_b->M_111 + X_003(dx) * m_b->M_110 + X_010(dx) * m_b->M_103 + - X_011(dx) * m_b->M_102 + X_012(dx) * m_b->M_101 + - /*X_013(dx) * m_b->M_100 +*/ X_100(dx) * m_b->M_013 + - X_101(dx) * m_b->M_012 + - X_102(dx) * m_b->M_011 + /*X_103(dx) * m_b->M_010 +*/ + X_011(dx) * m_b->M_102 + + X_012(dx) * m_b->M_101 /* + X_013(dx) * m_b->M_100 */ + + X_100(dx) * m_b->M_013 + X_101(dx) * m_b->M_012 + + X_102(dx) * m_b->M_011 /* + X_103(dx) * m_b->M_010 */ + X_110(dx) * m_b->M_003 + - X_111(dx) * m_b->M_002 + /* X_112(dx) * m_b->M_001 +*/ + X_111(dx) * m_b->M_002 /* + X_112(dx) * m_b->M_001 */ + X_113(dx) * m_b->M_000; m_a->M_122 = m_b->M_122 + X_001(dx) * m_b->M_121 + X_002(dx) * m_b->M_120 + X_010(dx) * m_b->M_112 + X_011(dx) * m_b->M_111 + X_012(dx) * m_b->M_110 + X_020(dx) * m_b->M_102 + - X_021(dx) * m_b->M_101 + /*X_022(dx) * m_b->M_100 +*/ + X_021(dx) * m_b->M_101 /* + X_022(dx) * m_b->M_100 */ + X_100(dx) * m_b->M_022 + X_101(dx) * m_b->M_021 + X_102(dx) * m_b->M_020 + X_110(dx) * m_b->M_012 + - X_111(dx) * m_b->M_011 + /*X_112(dx) * m_b->M_010 +*/ - X_120(dx) * m_b->M_002 + - /*X_121(dx) * m_b->M_001+*/ X_122(dx) * m_b->M_000; - m_a->M_131 = - m_b->M_131 + X_001(dx) * m_b->M_130 + X_010(dx) * m_b->M_121 + - X_011(dx) * m_b->M_120 + X_020(dx) * m_b->M_111 + X_021(dx) * m_b->M_110 + - X_030(dx) * m_b->M_101 + - /*X_031(dx) * m_b->M_100 +*/ X_100(dx) * m_b->M_031 + - X_101(dx) * m_b->M_030 + X_110(dx) * m_b->M_021 + X_111(dx) * m_b->M_020 + - X_120(dx) * - m_b->M_011 + /*X_121(dx) * m_b->M_010 + X_130(dx) * m_b->M_001 +*/ - X_131(dx) * m_b->M_000; + X_111(dx) * m_b->M_011 /* + X_112(dx) * m_b->M_010 */ + + X_120(dx) * m_b->M_002 /* + X_121(dx) * m_b->M_001 */ + + X_122(dx) * m_b->M_000; + m_a->M_131 = m_b->M_131 + X_001(dx) * m_b->M_130 + X_010(dx) * m_b->M_121 + + X_011(dx) * m_b->M_120 + X_020(dx) * m_b->M_111 + + X_021(dx) * m_b->M_110 + + X_030(dx) * m_b->M_101 /* + X_031(dx) * m_b->M_100 */ + + X_100(dx) * m_b->M_031 + X_101(dx) * m_b->M_030 + + X_110(dx) * m_b->M_021 + X_111(dx) * m_b->M_020 + + X_120(dx) * m_b->M_011 /* + X_121(dx) * m_b->M_010 */ + /* + X_130(dx) * m_b->M_001 */ + + X_131(dx) * m_b->M_000; m_a->M_140 = m_b->M_140 + X_010(dx) * m_b->M_130 + X_020(dx) * m_b->M_120 + - X_030(dx) * m_b->M_110 + /*X_040(dx) * m_b->M_100 +*/ + X_030(dx) * m_b->M_110 /* + X_040(dx) * m_b->M_100 */ + X_100(dx) * m_b->M_040 + X_110(dx) * m_b->M_030 + - X_120(dx) * m_b->M_020 + /*X_130(dx) * m_b->M_010 +*/ + X_120(dx) * m_b->M_020 /* + X_130(dx) * m_b->M_010 */ + X_140(dx) * m_b->M_000; m_a->M_203 = m_b->M_203 + X_001(dx) * m_b->M_202 + X_002(dx) * m_b->M_201 + X_003(dx) * m_b->M_200 + X_100(dx) * m_b->M_103 + - X_101(dx) * m_b->M_102 + X_102(dx) * m_b->M_101 + - /*X_103(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_003 + - X_201(dx) * m_b->M_002 + - /* X_202(dx) * m_b->M_001 +*/ X_203(dx) * m_b->M_000; + X_101(dx) * m_b->M_102 + + X_102(dx) * m_b->M_101 /* + X_103(dx) * m_b->M_100 */ + + X_200(dx) * m_b->M_003 + + X_201(dx) * m_b->M_002 /* + X_202(dx) * m_b->M_001 */ + + X_203(dx) * m_b->M_000; m_a->M_212 = m_b->M_212 + X_001(dx) * m_b->M_211 + X_002(dx) * m_b->M_210 + X_010(dx) * m_b->M_202 + X_011(dx) * m_b->M_201 + X_012(dx) * m_b->M_200 + X_100(dx) * m_b->M_112 + X_101(dx) * m_b->M_111 + X_102(dx) * m_b->M_110 + X_110(dx) * m_b->M_102 + - X_111(dx) * m_b->M_101 + /*X_112(dx) * m_b->M_100 +*/ + X_111(dx) * m_b->M_101 /* + X_112(dx) * m_b->M_100 */ + X_200(dx) * m_b->M_012 + - X_201(dx) * m_b->M_011 + /*X_202(dx) * m_b->M_010 +*/ - X_210(dx) * m_b->M_002 + - /*X_211(dx) * m_b->M_001 +*/ X_212(dx) * m_b->M_000; - m_a->M_221 = - m_b->M_221 + X_001(dx) * m_b->M_220 + X_010(dx) * m_b->M_211 + - X_011(dx) * m_b->M_210 + X_020(dx) * m_b->M_201 + X_021(dx) * m_b->M_200 + - X_100(dx) * m_b->M_121 + X_101(dx) * m_b->M_120 + X_110(dx) * m_b->M_111 + - X_111(dx) * m_b->M_110 + - X_120(dx) * m_b->M_101 + /*X_121(dx) * m_b->M_100 +*/ - X_200(dx) * m_b->M_021 + X_201(dx) * m_b->M_020 + X_210(dx) * m_b->M_011 + - /*X_211(dx) * m_b->M_010 + X_220(dx) * m_b->M_001 +*/ X_221(dx) * - m_b->M_000; + X_201(dx) * m_b->M_011 /* + X_202(dx) * m_b->M_010 */ + + X_210(dx) * m_b->M_002 /* + X_211(dx) * m_b->M_001 */ + + X_212(dx) * m_b->M_000; + m_a->M_221 = m_b->M_221 + X_001(dx) * m_b->M_220 + X_010(dx) * m_b->M_211 + + X_011(dx) * m_b->M_210 + X_020(dx) * m_b->M_201 + + X_021(dx) * m_b->M_200 + X_100(dx) * m_b->M_121 + + X_101(dx) * m_b->M_120 + X_110(dx) * m_b->M_111 + + X_111(dx) * m_b->M_110 + + X_120(dx) * m_b->M_101 /* + X_121(dx) * m_b->M_100 */ + + X_200(dx) * m_b->M_021 + X_201(dx) * m_b->M_020 + + X_210(dx) * m_b->M_011 /* + X_211(dx) * m_b->M_010 */ + /* + X_220(dx) * m_b->M_001 */ + + X_221(dx) * m_b->M_000; m_a->M_230 = m_b->M_230 + X_010(dx) * m_b->M_220 + X_020(dx) * m_b->M_210 + X_030(dx) * m_b->M_200 + X_100(dx) * m_b->M_130 + - X_110(dx) * m_b->M_120 + X_120(dx) * m_b->M_110 + - /*X_130(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_030 + - X_210(dx) * m_b->M_020 + - /*X_220(dx) * m_b->M_010 +*/ X_230(dx) * m_b->M_000; + X_110(dx) * m_b->M_120 + + X_120(dx) * m_b->M_110 /* + X_130(dx) * m_b->M_100 */ + + X_200(dx) * m_b->M_030 + + X_210(dx) * m_b->M_020 /* + X_220(dx) * m_b->M_010 */ + + X_230(dx) * m_b->M_000; m_a->M_302 = m_b->M_302 + X_001(dx) * m_b->M_301 + X_002(dx) * m_b->M_300 + X_100(dx) * m_b->M_202 + X_101(dx) * m_b->M_201 + X_102(dx) * m_b->M_200 + X_200(dx) * m_b->M_102 + - X_201(dx) * m_b->M_101 + /*X_202(dx) * m_b->M_100 +*/ - X_300(dx) * m_b->M_002 + - /*X_301(dx) * m_b->M_001 +*/ X_302(dx) * m_b->M_000; - m_a->M_311 = - m_b->M_311 + X_001(dx) * m_b->M_310 + X_010(dx) * m_b->M_301 + - X_011(dx) * m_b->M_300 + X_100(dx) * m_b->M_211 + X_101(dx) * m_b->M_210 + - X_110(dx) * m_b->M_201 + X_111(dx) * m_b->M_200 + X_200(dx) * m_b->M_111 + - X_201(dx) * m_b->M_110 + - X_210(dx) * m_b->M_101 + /*X_211(dx) * m_b->M_100 +*/ - X_300(dx) * - m_b->M_011 + /*X_301(dx) * m_b->M_010 + X_310(dx) * m_b->M_001 +*/ - X_311(dx) * m_b->M_000; + X_201(dx) * m_b->M_101 /* + X_202(dx) * m_b->M_100 */ + + X_300(dx) * m_b->M_002 /* + X_301(dx) * m_b->M_001 */ + + X_302(dx) * m_b->M_000; + m_a->M_311 = m_b->M_311 + X_001(dx) * m_b->M_310 + X_010(dx) * m_b->M_301 + + X_011(dx) * m_b->M_300 + X_100(dx) * m_b->M_211 + + X_101(dx) * m_b->M_210 + X_110(dx) * m_b->M_201 + + X_111(dx) * m_b->M_200 + X_200(dx) * m_b->M_111 + + X_201(dx) * m_b->M_110 + + X_210(dx) * m_b->M_101 /* + X_211(dx) * m_b->M_100 */ + + X_300(dx) * m_b->M_011 /* + X_301(dx) * m_b->M_010 */ + /* + X_310(dx) * m_b->M_001 */ + + X_311(dx) * m_b->M_000; m_a->M_320 = m_b->M_320 + X_010(dx) * m_b->M_310 + X_020(dx) * m_b->M_300 + X_100(dx) * m_b->M_220 + X_110(dx) * m_b->M_210 + X_120(dx) * m_b->M_200 + X_200(dx) * m_b->M_120 + - X_210(dx) * m_b->M_110 + /*X_220(dx) * m_b->M_100 +*/ - X_300(dx) * m_b->M_020 + - /*X_310(dx) * m_b->M_010 +*/ X_320(dx) * m_b->M_000; + X_210(dx) * m_b->M_110 /* + X_220(dx) * m_b->M_100 */ + + X_300(dx) * m_b->M_020 /* + X_310(dx) * m_b->M_010 */ + + X_320(dx) * m_b->M_000; m_a->M_401 = m_b->M_401 + X_001(dx) * m_b->M_400 + X_100(dx) * m_b->M_301 + X_101(dx) * m_b->M_300 + X_200(dx) * m_b->M_201 + - X_201(dx) * m_b->M_200 + X_300(dx) * m_b->M_101 + - /*X_301(dx) * m_b->M_100 + X_400(dx) * m_b->M_001 +*/ - X_401(dx) * m_b->M_000; + X_201(dx) * m_b->M_200 + + X_300(dx) * m_b->M_101 /* + X_301(dx) * m_b->M_100 */ + /* + X_400(dx) * m_b->M_001 */ + + X_401(dx) * m_b->M_000; m_a->M_410 = m_b->M_410 + X_010(dx) * m_b->M_400 + X_100(dx) * m_b->M_310 + X_110(dx) * m_b->M_300 + X_200(dx) * m_b->M_210 + - X_210(dx) * m_b->M_200 + X_300(dx) * m_b->M_110 + - /*X_310(dx) * m_b->M_100 + X_400(dx) * m_b->M_010 +*/ - X_410(dx) * m_b->M_000; + X_210(dx) * m_b->M_200 + + X_300(dx) * m_b->M_110 /* + X_310(dx) * m_b->M_100 */ + /* + X_400(dx) * m_b->M_010 */ + + X_410(dx) * m_b->M_000; m_a->M_500 = m_b->M_500 + X_100(dx) * m_b->M_400 + X_200(dx) * m_b->M_300 + - X_300(dx) * m_b->M_200 + /*X_400(dx) * m_b->M_100 +*/ + X_300(dx) * m_b->M_200 /* + X_400(dx) * m_b->M_100 */ + X_500(dx) * m_b->M_000; #endif #if SELF_GRAVITY_MULTIPOLE_ORDER > 5 diff --git a/theory/Multipoles/generate_multipoles/multipoles.py b/theory/Multipoles/generate_multipoles/multipoles.py index 40834ac8592e98fb5ff73e7ea43717e03cef1da8..9c9f5513a728ca4fad5ad8771cf6fde1603f59f9 100644 --- a/theory/Multipoles/generate_multipoles/multipoles.py +++ b/theory/Multipoles/generate_multipoles/multipoles.py @@ -240,7 +240,7 @@ print("-------------------------------------------------\n") if order > 0: print("#if SELF_GRAVITY_MULTIPOLE_ORDER > %d" % (order - 1)) -print("/* Shift %s order terms */" % ordinal(order)) +print("/* Shift %s order terms (1st order mpole (all 0) commented out) */" % ordinal(order)) # Create all the terms relevent for this order for i in range(order + 1): @@ -262,11 +262,18 @@ for i in range(order + 1): and jj + jjj == j and kk + kkk == k ): - print( - "+ X_%d%d%d(dx) * m_b->M_%d%d%d" - % (ii, jj, kk, iii, jjj, kkk), - end=" ", - ) + if iii + jjj + kkk == 1: + print( + "/* + X_%d%d%d(dx) * m_b->M_%d%d%d */" + % (ii, jj, kk, iii, jjj, kkk), + end=" ", + ) + else: + print( + "+ X_%d%d%d(dx) * m_b->M_%d%d%d" + % (ii, jj, kk, iii, jjj, kkk), + end=" ", + ) print(";")