mullocal.c

很经典的电磁计算（电容计算）软件 MIT90年代开发

/*
Copyright (c) 1990 Massachusetts Institute of Technology, Cambridge, MA.
All rights reserved.

This Agreement gives you, the LICENSEE, certain rights and obligations.
By using the software, you indicate that you have read, understood, and
will comply with the terms.

Permission to use, copy and modify for internal, noncommercial purposes
is hereby granted.  Any distribution of this program or any part thereof
is strictly prohibited without prior written consent of M.I.T.

Title to copyright to this software and to any associated documentation
shall at all times remain with M.I.T. and LICENSEE agrees to preserve
same.  LICENSEE agrees not to make any copies except for LICENSEE'S
internal noncommercial use, or to use separately any portion of this
software without prior written consent of M.I.T.  LICENSEE agrees to
place the appropriate copyright notice on any such copies.

Nothing in this Agreement shall be construed as conferring rights to use
in advertising, publicity or otherwise any trademark or the name of
"Massachusetts Institute of Technology" or "M.I.T."

M.I.T. MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.  By
way of example, but not limitation, M.I.T. MAKES NO REPRESENTATIONS OR
WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR
THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS OR DOCUMENTATION WILL
NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS.
M.I.T. shall not be held liable for any liability nor for any direct,
indirect or consequential damages with respect to any claim by LICENSEE
or any third party on account of or arising from this Agreement or use
of this software.
*/

#include "mulGlobal.h"

/*
  globals used for temporary storage
*/
double **facFrA;		/* array of factorial fractions */
double *cosmkB;			/* array used to look up cos[(m+-k)beta] */
double *sinmkB;			/* array used to look up sin[(m+-k)beta] */

/*
  initializes the factorial fraction array used in M2L, L2L matrix calculation
*/
void evalFacFra(array, order)
int order;			/* array is 2*order+1 x 2*order+1 */
double **array;			/* array[num][den] = num!/den! */
{
  int d, i;
  for(i = 0; i <= 2*order; i++) {
    array[i][i] = 1.0; /* do main diagonal */
    if(i > 0 && i < 2*order) array[i+1][i] = i+1; /* do first sub diagonal */
  }
  for(d = 3; d <= 2*order; d++) { /* loop on lower triangular rows */
    for(i = 1; i < d-1; i++) {	/* loop on columns */
      array[d][i] = array[d-1][i]*array[d][d-1];
    }
  }
  /* invert lower part entries and copy to top */
  for(d = 2; d <= 2*order; d++) {
    for(i = 1; i <= d-1; i++) {
      array[i][d] = 1/array[d][i];
    }
  }
  /* copy 1st row and column from computed values */
  for(d = 1; d <= 2*order; d++) {
    array[0][d] = array[1][d];
    array[d][0] = array[d][1];
  }

#if DISFAC == ON
  fprintf(stdout, "FACTORIAL FRACTION ARRAY:\n");
  dumpMat(array, 2*order+1, 2*order+1);
#endif

}

/*
  initializes sqrt((m+n)!/(n-m)!) lookup table (for L2L)
*/
void evalSqrtFac(arrayout, arrayin, order)
int order;
double **arrayout, **arrayin;
{
  int n, m;			/* arrayout[n][m] = sqrt((m+n)!/(n-m)!) */

  /* set up first column, always ones */
  for(n = 0; n < order+1; n++) arrayout[n][0] = 1.0;

  /* set up lower triangular (n+m)!/(n-m)! */
  for(n = 1; n <= order; n++) {
    for(m = 1; m <= n; m++) {
      arrayout[n][m] = sqrt(arrayin[n][m]);
    }
  }

#if DISSFA == ON
  fprintf(stdout, "SQUARE ROOT FACTORIAL ARRAY:\n");
  dumpMat(arrayout, order+1, order+1);
#endif

}


/*
  initializes cos[(m+-k)beta] and sin[(m+-k)beta] lookup tables (M2L and L2L)
*/
void evalSinCos(beta, order)
int order;
double beta;
{
  int i;
  double temp = beta;

  for(i = 1; i <= 2*order; beta += temp, i++) {
    sinmkB[i] = sin(beta);
    cosmkB[i] = cos(beta);
  }
}

/*
  looks up sin[(m+-k)beta]
*/
double sinB(sum)
int sum;
{
  if(sum < 0) return(-sinmkB[abs(sum)]);
  else return(sinmkB[sum]);
}

/*
  looks up cos[(m+-k)beta]
*/
double cosB(sum)
int sum;
{
  return(cosmkB[abs(sum)]);
}

/* 
  Used for all but no local downward pass. 
*/
double **mulMulti2Local(x, y, z, xp, yp, zp, order)
int order;
double x, y, z, xp, yp, zp;	/* multipole and local cube centers */
{
  int i, j, k, n, m;
  int terms = multerms(order);	/* the number of non-zero moments */
  int ct = costerms(order);	/* the number of non-zero cos (bar) moments */
  double **mat;			/* the transformation matrix */
  double rho, cosA, beta;	/* spher. position of multi rel to local */
  double rhoJ, rhoN;		/* rho^j and (-1)^n*rho^(n+1) in main loop */
  double rhoFac;		/* = rhoJ*rhoN intermediate storage */
  double temp1, temp2, temp3;
  double iPwr(), sinB(), cosB();
  extern double *tleg, *Ir, *Irn, *phi, *Mphi; /* external temporary storage */

  /* allocate the multi to local transformation matrix */
  CALLOC(mat, terms, double*, ON, AM2L);
  for(i = 0; i < terms; i++)
      CALLOC(mat[i], terms, double, ON, AM2L);

  /* find relative spherical coordinates */
  xyz2sphere(x, y, z, xp, yp, zp, &rho, &cosA, &beta);

  /* generate legendre function evaluations */
  evalLegendre(cosA, tleg, 2*order); /* multi->loc needs 2x legendres */

  /* generate sin[(m+-k)beta] and cos[(m+-k)beta] look up arrays */
  /*  other lookup arrays generated in mulMultiAlloc() */
  evalSinCos(beta, order);

  /* generate multi to local transformation matrix; uses NB12 pg30 */
  /*  rhoFac factor divides could be reduced to once per loop */
  for(j = 0, rhoJ = 1.0; j <= order; rhoJ *= rho, j++) {
    for(k = 0; k <= j; k++) {	/* loop on Nj^k's, local exp moments */
      for(n = 0, rhoN = rho; n <= order; rhoN *= (-rho), n++) {
	for(m = 0; m <= n; m++) { /* loop on On^m's, multipole moments */

	  /* generate a bar(N)j^k and dblbar(N)j^k entry */
	  rhoFac = rhoJ*rhoN;	/* divisor to give (-1)^n/rho^(j+n+1) factor */
	  if(k == 0) {		/* use abbreviated formulae in this case */

	    /* generate only bar(N)j^0 entry (dblbar(N)j^0 = 0 always) */
	    if(m != 0) {
	      temp1 = tleg[CINDEX(j+n, m)]*facFrA[j+n-m][n+m];
	      mat[CINDEX(j, 0)][CINDEX(n, m)] += temp1*cosB(m)/rhoFac;
	      mat[CINDEX(j, 0)][SINDEX(n, m, ct)] += temp1*sinB(m)/rhoFac;
	    }
	    else mat[CINDEX(j, 0)][CINDEX(n, 0)] 
		+= tleg[CINDEX(j+n, 0)]*facFrA[j+n][n]/rhoFac;
	  }
	  else {
	    temp1 = tleg[CINDEX(j+n, abs(m-k))]
		*facFrA[j+n-abs(m-k)][n+m]*iPwr(abs(k-m)-k-m);
	    temp2 = tleg[CINDEX(j+n, m+k)]*facFrA[j+n-m-k][n+m];
	    temp3 = tleg[CINDEX(j+n, k)]*facFrA[j+n-k][n]*2;

	    /* generate bar(N)j^k entry */
	    if(m != 0) {
	      mat[CINDEX(j, k)][CINDEX(n, m)] 
		  += (temp1*cosB(m-k)+temp2*cosB(m+k))/rhoFac;
	      mat[CINDEX(j, k)][SINDEX(n, m, ct)] 
		  += (temp1*sinB(m-k)+temp2*sinB(m+k))/rhoFac;
	    }
	    else mat[CINDEX(j, k)][CINDEX(n, 0)] += temp3*cosB(k)/rhoFac;

	    /* generate dblbar(N)j^k entry */
	    if(m != 0) {
	      mat[SINDEX(j, k, ct)][CINDEX(n, m)] 
		  += (-temp1*sinB(m-k)+temp2*sinB(m+k))/rhoFac;
	      mat[SINDEX(j, k, ct)][SINDEX(n, m, ct)] 
		  += (temp1*cosB(m-k)-temp2*cosB(m+k))/rhoFac;
	    }
	    else mat[SINDEX(j, k, ct)][CINDEX(n, 0)] += temp3*sinB(k)/rhoFac;
	  }
	}
      }
    }
  }

#if DISM2L == ON
  dispM2L(mat, x, y, z, xp, yp, zp, order);
#endif

  return(mat);
}

/* 
  Used only for true (Greengard) downward pass - similar to Multi2Local
*/
double **mulLocal2Local(x, y, z, xc, yc, zc, order)
int order;
double x, y, z, xc, yc, zc;	/* parent and child cube centers */
{
  int i, j, k, n, m;
  int terms = multerms(order);	/* the number of non-zero moments */
  int ct = costerms(order);	/* the number of non-zero cos (bar) moments */
  double **mat;			/* the transformation matrix */
  double rho, cosA, beta;	/* spher. position of multi rel to local */
  double rhoJ, rhoN;		/* rho^j and (-1)^n*rho^(n+1) in main loop */
  double rhoFac;		/* = rhoJ*rhoN intermediate storage */
  double temp1, temp2, temp3;
  double iPwr(), sinB(), cosB();
  extern double *tleg, *Ir, *Irn, *phi, *Mphi; /* external temporary storage */

  /* allocate the local to local transformation matrix */
  CALLOC(mat, terms, double*, ON, AL2L);
  for(i = 0; i < terms; i++)
      CALLOC(mat[i], terms, double, ON, AL2L);

  /* find relative spherical coordinates */
  xyz2sphere(x, y, z, xc, yc, zc, &rho, &cosA, &beta);

  /* generate legendre function evaluations */
  evalLegendre(cosA, tleg, 2*order); /* local->local needs 2x legendres */

  /* generate sin[(m+-k)beta] and cos[(m+-k)beta] look up arrays */
  /*  other lookup arrays generated in mulMultiAlloc() */
  evalSinCos(beta, order);

  /* generate local to local transformation matrix; uses NB12 pg36Y */
  /*  rhoFac factor divides could be reduced to once per loop */
  for(j = 0, rhoJ = 1.0; j <= order; rhoJ *= (-rho), j++) {