mulmulti.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 *Irn, *Mphi;		/* (1/r)^n+1, m*phi vect's */
double *Ir, *phi;		/* 1/r and phi arrays, used to update above */
double *Rho, *Rhon;		/* rho and rho^n array */
double *Beta, *Betam;		/* beta and beta*m array */
double *tleg;		/* Temporary Legendre storage. */
double **factFac;		/* factorial factor array: (n-m+1)...(n+m) */

/* 
   Used various places.  Returns number of coefficients in the multipole 
   expansion. 
*/
int multerms(order)
int order;
{
  return(costerms(order) + sinterms(order));
}

/*
  returns number of cos(m*phi)-weighted terms in the real (not cmpx) multi exp
*/
int costerms(order)
int order;
{
  return(((order + 1) * (order + 2)) / 2);
}

/*
  returns number of sin(m*phi)-weighted terms in the real (not cmpx) multi exp
*/
int sinterms(order)
int order;
{
  return((((order + 1) * (order + 2)) / 2) - (order+1));
}


/*
  takes two sets of cartesian absolute coordinates; finds rel. spherical coor.
*/
void xyz2sphere(x, y, z, x0, y0, z0, rho, cosA, beta)
double x, y, z, x0, y0, z0, *rho, *cosA, *beta;
{
  /* get relative coordinates */
  x -= x0;			/* "0" coordinates play the role of origin */
  y -= y0;
  z -= z0;
  /* get spherical coordinates */
  *rho = sqrt(x*x + y*y + z*z);

  if(*rho == 0.0) *cosA = 1.0;
  else *cosA = z/(*rho);

  if(x == 0.0 && y == 0.0) *beta = 0.0;
  else *beta = atan2(y, x);

}

/*
  gives the linear index into vector from n and m used by all routines dealing
  with moments (cosine parts) and Leg. function evals, e.g. Mn^m and Pn^m(cosA)
  used for all cases except for the sine part of arrays (use sindex() instead)
  assumed entry order: (n,m) = (0,0) (1,0) (1,1) (2,0) (2,1) (2,2) (3,0)...
*/
/* REPLACED BY MACRO CINDEX(N, M) 24July91 */

#if TRUE == FALSE

int index(n, m)
int n, m;
{
  if(m > n) {
    fprintf(stderr, "index: m = %d > n = %d\n", m, n);
    exit(1);
  }
  if(n < 0 || m < 0) {
    fprintf(stderr, "index: n = %d or m = %d negative\n", n, m);
    exit(1);
  }
  return(m + (n*(n+1))/2);
}

/*
  gives the linear index into vector from n and m used by all routines dealing
  with moments (sine parts), e.g. Mn^m 
  assumes an array with all m = 0 (Mn^0) entries omitted to save space
  assumed entry order: (n,m) = (1,1) (2,1) (2,2) (3,1) (3,2) (3,3) (4,1)...
*/
/* REPLACED BY MACRO SINDEX(N, M, CTERMS) 24July91 */
int sindex(n, m, cterms)
int n, m, cterms;		/* cterms is costerms(order) */
{
  if(m > n) {
    fprintf(stderr, "sindex: m = %d > n = %d\n", m, n);
    exit(1);
  }
  if(n < 0 || m < 0) {
    fprintf(stderr, "sindex: n = %d or m = %d negative\n", n, m);
    exit(1);
  }
  if(m == 0) {
    fprintf(stderr, "sindex: attempt to index M%d^0\n", n);
    exit(1);
  }
  return(cterms + m + (n*(n+1))/2 - (n+1));
}

#endif				/* end of index(), sindex() comment out */

/*
  returns i = sqrt(-1) to the power of the argument
*/
double iPwr(e)
int e;				/* exponent, computes i^e */
{
  if(e == 0) return(1.0);
  if(e % 2 != 0) {
    fprintf(stderr, "iPwr: odd exponent %d\n", e);
    exit(1);
  }
  else {
    e = e/2;			/* get power of negative 1 */
    if(e % 2 == 0) return(1.0);
    else return(-1.0);
  }
}

/*
  returns factorial of the argument (x!)
*/
double fact(x)
int x;
{
  double ret = 1.0;
  if(x == 0 || x == 1) return(1.0);
  else if(x < 0) {
    fprintf(stderr, "fact: attempt to take factorial of neg no. %d\n", x);
    exit(1);
  }
  else {
    while(x > 1) {
      ret *= x;
      x--;
    }
    return(ret);
  }
}

/*
  produces factorial factor array for mulMulti2P
*/
void evalFactFac(array, order)
int order;
double **array;
{
  int n, m;			/* array[n][m] = (m+n)!/(n-m)! */

  /* do first column of lower triangular part - always 1's */
  for(n = 0; n < order+1; n++) array[n][0] = 1.0;

  /* do remaining columns of lower triangular part */
  /*  use (n-m)!/(n+m)! = 1/(n-m+1)...(n+m) since m \leq n */
  /*  (array entry is divided into number to effect mul by (n-m)!/(n+m)!) */
  for(n = 1; n <= order; n++) {
    for(m = 1; m <= n; m++) {
      array[n][m] = (n-(m-1)) * array[n][m-1] * (n+m);
    }
  }

#if DISFAF == ON
  fprintf(stdout, "FACTORIAL FACTOR ARRAY:\n");
  dumpMat(array, order+1, order+1);
#endif

}

/*
  Allocates space for temporary vectors.
*/
void mulMultiAlloc(maxchgs, order, depth)
int maxchgs, order, depth;
{
  int x;

#if DISSYN == ON
  extern int *multicnt, *localcnt, *evalcnt;
#endif
#if DMTCNT == ON
  extern int **Q2Mcnt, **Q2Lcnt, **Q2Pcnt, **L2Lcnt;
  extern int **M2Mcnt, **M2Lcnt, **M2Pcnt, **L2Pcnt, **Q2PDcnt;
#endif
  extern double *sinmkB, *cosmkB, **facFrA;

  if(maxchgs > 0) {
    CALLOC(Rho, maxchgs, double, ON, AMSC); /* rho array */
    CALLOC(Rhon, maxchgs, double, ON, AMSC); /* rho^n array */
    CALLOC(Beta, maxchgs, double, ON, AMSC); /* beta array */
    CALLOC(Betam, maxchgs, double, ON, AMSC); /* beta*m array */
    CALLOC(Irn, maxchgs, double, ON, AMSC);	/* (1/r)^n+1 vector */
    CALLOC(Ir, maxchgs, double, ON, AMSC); /* 1/r vector */
    CALLOC(Mphi, maxchgs, double, ON, AMSC); /* m*phi vector */
    CALLOC(phi, maxchgs, double, ON, AMSC); /* phi vector */
  }
  CALLOC(tleg, costerms(2*order), double, ON, AMSC);
	       	/* temp legendre storage (2*order needed for local exp) */
  CALLOC(factFac, order+1, double*, ON, AMSC);
  for(x = 0; x < order+1; x++) {
    CALLOC(factFac[x], order+1, double, ON, AMSC);
  }
  evalFactFac(factFac, order);	/* get factorial factors for mulMulti2P */

#if DISSYN == ON
  /* for counts of local/multipole expansions and eval mat builds by level */
  CALLOC(localcnt, depth+1, int, ON, AMSC);
  CALLOC(multicnt, depth+1, int, ON, AMSC);
  CALLOC(evalcnt, depth+1, int, ON, AMSC);
#endif

#if DMTCNT == ON
  /* for counts of transformation matrices by level  */
  CALLOC(Q2Mcnt, depth+1, int*, ON, AMSC);
  CALLOC(Q2Lcnt, depth+1, int*, ON, AMSC);
  CALLOC(Q2Pcnt, depth+1, int*, ON, AMSC);
  CALLOC(L2Lcnt, depth+1, int*, ON, AMSC);
  CALLOC(M2Mcnt, depth+1, int*, ON, AMSC); 
  CALLOC(M2Lcnt, depth+1, int*, ON, AMSC); 
  CALLOC(M2Pcnt, depth+1, int*, ON, AMSC); 
  CALLOC(L2Pcnt, depth+1, int*, ON, AMSC);
  CALLOC(Q2PDcnt, depth+1, int*, ON, AMSC);
  for(x = 0; x < depth+1; x++) {
    CALLOC(Q2Mcnt[x], depth+1, int, ON, AMSC);
    CALLOC(Q2Lcnt[x], depth+1, int, ON, AMSC);
    CALLOC(Q2Pcnt[x], depth+1, int, ON, AMSC);
    CALLOC(L2Lcnt[x], depth+1, int, ON, AMSC);
    CALLOC(M2Mcnt[x], depth+1, int, ON, AMSC);
    CALLOC(M2Lcnt[x], depth+1, int, ON, AMSC);
    CALLOC(M2Pcnt[x], depth+1, int, ON, AMSC);
    CALLOC(L2Pcnt[x], depth+1, int, ON, AMSC);
    CALLOC(Q2PDcnt[x], depth+1, int, ON, AMSC);
  }
#endif

  /* from here down could be switched out when the fake dwnwd pass is used */
  CALLOC(facFrA, 2*order+1, double*, ON, AMSC);
  for(x = 0; x < 2*order+1; x++) {
    CALLOC(facFrA[x], 2*order+1, double, ON, AMSC);
  }
  /* generate table of factorial fraction evaluations (for M2L and L2L) */
  evalFacFra(facFrA, order);
  CALLOC(sinmkB, 2*order+1, double, ON, AMSC); /* sin[(m+-k)beta] */
  CALLOC(cosmkB, 2*order+1, double, ON, AMSC); /* cos[(m+-k)beta] */
  cosmkB[0] = 1.0;		/* look up arrays used for local exp */
  /* generate array of sqrt((n+m)!(n-m)!)'s for L2L
  evalSqrtFac(sqrtFac, factFac, order); */
}

/*
  returns an vector of Legendre function evaluations of form Pn^m(cosA)
  n and m have maximum value order
  vector entries correspond to (n,m) = (0,0) (1,0) (1,1) (2,0) (2,1)...
*/
void evalLegendre(cosA, vector, order)
double cosA, *vector;
int order;
{
  int x;
  int n, m;			/* as in Pn^m, both <= order */
  double sinMA;			/* becomes sin^m(alpha) in higher order P's */
  double fact;			/* factorial factor */

  /* do evaluations of first four functions separately w/o recursions */
  vector[CINDEX(0, 0)] = 1.0;	/* P0^0 */
  if(order > 0) {
    vector[CINDEX(1, 0)] = cosA;	/* P1^0 */
    vector[CINDEX(1, 1)] = sinMA = -sqrt(1-cosA*cosA); /* P1^1 = -sin(alpha) */
  }
  if(order > 1) vector[CINDEX(2, 1)] = 3*sinMA*cosA; /* P2^1 = -3sin()cos() */

  /* generate remaining evaluations by recursion on lower triangular array */
  fact = 1.0;
  for(m = 0; m < order+1; m++) {
    if(m != 0 && m != 1) {	/* set up first two evaluations in row */
      fact *= (2*m - 1); /* (2(m-1)-1)!! -> (2m-1)!! */
      /* use recursion on m */
      if(vector[CINDEX(1, 1)] == 0.0) {
	vector[CINDEX(m, m)] = 0.0;
	if(m != order) vector[CINDEX(m+1, m)] = 0.0;	/* if not last row */
      }