Added a program to generate the terms required for the FMM calculation at any order.

examples/theory/expandTerms.c

+#include <stdio.h>
+#include <stdlib.h>
+
+/* Naive implementation ... */
+int fact1(int x) {
+  if (x == 0) return 1;
+  else return x*fact1(x-1);
+}
+
+int fact(int a, int b, int c) {
+  return fact1(a)*fact1(b)*fact1(c);
+}
+
+int main(int argc, char* argv[]) {
+
+  int p = atoi(argv[1]);
+  int a,b,c,n,nfact;
+  int aa,bb,cc,k, kfact;
+  int m, mfact;
+
+  printf("== Multipole shifting up to order p=%d ==\n\n", p);
+
+  /* Loop over multipole p */
+  for ( n = 0; n<=p; n++ ) {
+    
+    /* Loop over all possible multipole of that p */
+    for ( a = 0; a<=n; a++) {
+      for ( b = 0; b<=n; b++) {
+	for ( c = 0; c<=n; c++) {
+	 
+	  if (a + b + c == n ) {
+	    	    
+
+	    /* Ok, now implement formula */
+	    
+	    /* Loop over all possible terms on the RHS */
+	    for ( aa = 0; aa<=p; aa++) {
+	      for ( bb = 0; bb<=p; bb++) {
+		for ( cc = 0; cc<=p; cc++) {
+		  
+		  k = aa + bb + cc;
+		  kfact = fact(aa,bb,cc);
+
+		  /* Dehnen & Read eq. 55 */
+		  if ( (k <= n) && (a-aa >=0) && (b-bb >=0) && (c-cc >=0)  && k!=1 && n!=1) {
+		    printf("B->M_%d%d%d += " , a, b, c);
+
+		    
+		    if ( aa > 0 )
+		      printf("dx^%d",aa);
+		    if ( bb > 0 )
+		      printf("dy^%d",bb);
+		    if (cc > 0 )
+		      printf("dz^%d",cc);
+		    if ( kfact > 1 )
+		      printf("/%d.", kfact);
+		    if (  aa > 0 || bb >0 || cc > 0 )
+		      printf (" * ");
+		    
+		    
+		    printf( "A.M_%d%d%d;", a-aa, b-bb, c-cc);
+		    printf("\n");
+
+		    }
+		 
+		}
+	      }
+	    }
+	    
+	  }
+	  
+	}
+      }
+    }
+
+    printf("//----------\n");
+
+  }
+
+  printf("\n");
+  printf("== Field tensors up to order p=%d for x component ==\n\n", p);
+
+  /* Loop over multipole p */
+  for ( n = 0; n<=p; n++ ) {
+    
+    /* Loop over all possible multipole of that p */
+    for ( a = 0; a<=n; a++) {
+      for ( b = 0; b<=n; b++) {
+	for ( c = 0; c<=n; c++) {
+	 
+	  if (a + b + c == n ) {
+	    	    
+
+
+	    /* Ok, now implement formula */
+	    
+	    /* Loop over all possible terms on the RHS */
+	    for ( aa = 0; aa<=p; aa++) {
+	      for ( bb = 0; bb<=p; bb++) {
+		for ( cc = 0; cc<=p; cc++) {
+		  
+		  m = aa + bb + cc;
+		  mfact = fact(aa,bb,cc);
+
+		  /* Dehnen & Read eq. 58 */
+		  if ( (m <= p-n) && m !=1){// && (a-aa >=0) && (b-bb >=0) && (c-cc >=0) ) {
+		    
+		    printf("B->F_%d%d%d +=", a, b, c);		    
+		    printf(" A.M_%d%d%d", aa,bb,cc);
+		    printf("*D_%d%d%d(r_x, r_y, r_z, inv_r);", a+aa+1,b+bb,c+cc);
+		    printf("\n");
+		  }
+		}
+	      }
+	    }
+	  }
+	}
+      }
+    }
+
+    printf("//----------\n");
+
+  }
+
+  printf("\n");
+  printf("== Field tensors up to order p=%d for y component ==\n\n", p);
+
+  /* Loop over multipole p */
+  for ( n = 0; n<=p; n++ ) {
+    
+    /* Loop over all possible multipole of that p */
+    for ( a = 0; a<=n; a++) {
+      for ( b = 0; b<=n; b++) {
+	for ( c = 0; c<=n; c++) {
+	 
+	  if (a + b + c == n ) {
+	    	    
+
+
+	    /* Ok, now implement formula */
+	    
+	    /* Loop over all possible terms on the RHS */
+	    for ( aa = 0; aa<=p; aa++) {
+	      for ( bb = 0; bb<=p; bb++) {
+		for ( cc = 0; cc<=p; cc++) {
+		  
+		  m = aa + bb + cc;
+		  mfact = fact(aa,bb,cc);
+
+		  /* Dehnen & Read eq. 58 */
+		  if ( (m <= p-n) && m !=1){// && (a-aa >=0) && (b-bb >=0) && (c-cc >=0) ) {
+		    
+		    printf("B->F_%d%d%d +=", a, b, c);		    
+		    printf(" A.M_%d%d%d", aa,bb,cc);
+		    printf("*D_%d%d%d(r_x, r_y, r_z, inv_r);", a+aa,b+bb+1,c+cc);
+		    printf("\n");
+		  }
+		}
+	      }
+	    }
+	  }
+	}
+      }
+    }
+
+    printf("//----------\n");
+
+  }
+
+  printf("\n");
+  printf("== Field tensors up to order p=%d for z component ==\n\n", p);
+
+  /* Loop over multipole p */
+  for ( n = 0; n<=p; n++ ) {
+    
+    /* Loop over all possible multipole of that p */
+    for ( a = 0; a<=n; a++) {
+      for ( b = 0; b<=n; b++) {
+	for ( c = 0; c<=n; c++) {
+	 
+	  if (a + b + c == n ) {
+	    	    
+
+
+	    /* Ok, now implement formula */
+	    
+	    /* Loop over all possible terms on the RHS */
+	    for ( aa = 0; aa<=p; aa++) {
+	      for ( bb = 0; bb<=p; bb++) {
+		for ( cc = 0; cc<=p; cc++) {
+		  
+		  m = aa + bb + cc;
+		  mfact = fact(aa,bb,cc);
+
+		  /* Dehnen & Read eq. 58 */
+		  if ( (m <= p-n) && m !=1){// && (a-aa >=0) && (b-bb >=0) && (c-cc >=0) ) {
+		    
+		    printf("B->F_%d%d%d +=", a, b, c);		    
+		    printf(" A.M_%d%d%d", aa,bb,cc);
+		    printf("*D_%d%d%d(r_x, r_y, r_z, inv_r);", a+aa,b+bb,c+cc+1);
+		    printf("\n");
+		  }
+		}
+	      }
+	    }
+	  }
+	}
+      }
+    }
+    printf("//----------\n");
+  }
+
+
+
+  printf("\n");
+  printf("== Accelerations up to order p=%d ==\n\n", p);
+
+  /* Loop over multipole p */
+  for ( n = 0; n<=p; n++ ) {
+    
+    /* Loop over all possible multipole of that p */
+    for ( a = 0; a<=n; a++) {
+      for ( b = 0; b<=n; b++) {
+	for ( c = 0; c<=n; c++) {
+	 
+	  if (a + b + c == n ) {
+	    nfact = fact(a,b,c);
+
+	    /* Ok, now implement formula */
+	    
+
+	    /* Dehnen & Read eq. 57 */
+		    
+	    printf("a += ");
+	      
+	      if ( a > 0 )
+		printf("dx^%d",a);
+	      if ( b > 0 )
+		printf("dy^%d",b);
+	      if (c > 0 )
+		printf("dz^%d",c);
+	      if ( nfact > 1 )
+		printf("/%d.", nfact);
+	      if (  a > 0 || b >0 || c > 0 )
+		printf (" * ");
+
+	      printf("B.F_%d%d%d;", a,b,c);
+	      printf("\n");
+	    
+	  }
+	}
+      }
+    }
+
+    printf("//----------\n");
+
+  }
+
+  return 0;
+
+}