mulmulti.c

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

      else {
	cosA = vector[CINDEX(1,0)]/vector[CINDEX(1,1)]; /* cosA= -cot(theta) */
	sinMA *= vector[CINDEX(1,1)]; /*(-sin(alpha))^(m-1)->(-sin(alpha))^m*/
	vector[CINDEX(m, m)] = fact * sinMA;
	if(m != order) {		/* do if not on last row */
	  vector[CINDEX(m+1, m)] = vector[CINDEX(1, 0)]*(2*m+1)
	      *vector[CINDEX(m, m)];
	}
      }
    }
    for(x = 2; x < order-m+1; x++) { /* generate row of evals recursively */
      vector[CINDEX(x+m, m)] = 
	  ((2*(x+m)-1)*vector[CINDEX(1, 0)]*vector[CINDEX(x+m-1, m)]
	   - (x + 2*m - 1)*vector[CINDEX(x+m-2, m)])/x;
    }
  }
}
    
  
/* 
  Returns a matrix which gives a cube's multipole expansion when *'d by chg vec
  OPTIMIZATIONS USING is_dummy HAVE NOT BEEN COMPLETELY IMPLEMENTED
*/
double **mulQ2Multi(chgs, is_dummy, numchgs, x, y, z, order)
double x, y, z; 
charge **chgs;
int numchgs, order, *is_dummy;
{
  double **mat;
  double cosA;			/* cosine of elevation coordinate */
  int i, j, k, kold, n, m, start;
  int cterms = costerms(order), terms = multerms(order);

  /* Allocate the matrix. */
  CALLOC(mat, terms, double*, ON, AQ2M);
  for(i=0; i < terms; i++) 
      CALLOC(mat[i], numchgs, double, ON, AQ2M);

  /* get Legendre function evaluations, one set for each charge */
  /*  also get charge coordinates, set up for subsequent evals */
  for(j = 0; j < numchgs; j++) { /* for each charge */

    /* get cosA for eval; save rho, beta in rho^n and cos/sin(m*beta) arrays */
    xyz2sphere(chgs[j]->x, chgs[j]->y, chgs[j]->z,
	       x, y, z, &(Rho[j]), &cosA, &(Beta[j]));
    Rhon[j] = Rho[j]; /* init powers of rho_i's */
    Betam[j] = Beta[j];		/* init multiples of beta */
    evalLegendre(cosA, tleg, order);	/* write moments to temporary array */

    /* write a column of the matrix with each set of legendre evaluations */
    for(i = 0; i < cterms; i++) mat[i][j] = tleg[i]; /* copy for cos terms */
  }

#if DALQ2M == ON
  fprintf(stdout,
	  "\nQ2M MATRIX BUILD:\n    AFTER LEGENDRE FUNCTION EVALUATON\n");
  dumpMat(mat, terms, numchgs);
#endif

  /* some of this can be avoided using is_dummy to skip unused columns */
  /* add the rho^n factors to the cos matrix entries. */
  for(i = 1, k = kold = 2; i < cterms; i++) { /* loop on rows of matrix */
    for(j = 0; j < numchgs; j++) mat[i][j] *= Rhon[j]; /* mul in factor */
    k -= 1;
    if(k == 0) {		/* so that effective n varys appropriately */
      kold = k = kold + 1;
      for(j = 0; j < numchgs; j++) Rhon[j] *= Rho[j];	/* r^n-1 -> r^n */
    }
  }

#if DALQ2M == ON
  fprintf(stdout,"    AFTER ADDITION OF RHO^N FACTORS\n");
  dumpMat(mat, terms, numchgs);
#endif

  /* copy result to lower (sine) part of matrix */
  for(n = 1; n <= order; n++) { /* loop on rows of matrix */
    for(m = 1; m <= n; m++) {
      for(j = 0; j < numchgs; j++) { /* copy a row */
	mat[SINDEX(n, m, cterms)][j] = mat[CINDEX(n, m)][j];
      }
    }
  }

#if DALQ2M == ON
  fprintf(stdout,"    AFTER COPYING SINE (LOWER) HALF\n");
  dumpMat(mat, terms, numchgs);
#endif

  /* add factors of cos(m*beta) and sin(m*beta) to matrix entries */
  for(m = 1; m <= order; m++) {	/* lp on m in Mn^m (no m=0 since cos(0)=1) */
    for(n = m; n <= order; n++) { /* loop over rows with same m */
      for(j = 0; j < numchgs; j++) { /* add factors to a row */
	mat[CINDEX(n, m)][j] *= (2.0*cos(Betam[j]));   /* note factors of 2 */
	mat[SINDEX(n, m, cterms)][j] *= (2.0*sin(Betam[j]));
      }
    }
    for(j = 0; j < numchgs; j++) Betam[j] += Beta[j]; /* (m-1)*beta->m*beta */
  }

  /* THIS IS NOT VERY GOOD: zero out columns corresponding to dummy panels */
  for(j = 0; j < numchgs; j++) {
    if(is_dummy[j]) {
      for(i = 0; i < terms; i++) {
	mat[i][j] = 0.0;
      }
    }
  }

#if DISQ2M == ON
  dispQ2M(mat, chgs, numchgs, x, y, z, order);
#endif

  return(mat);
}

double **mulMulti2Multi(x, y, z, xp, yp, zp, order)
double x, y, z, xp, yp, zp;	/* cube center, parent cube center */
int order;
{
  double **mat, rho, rhoPwr, cosA, beta, mBeta, temp1, temp2; 
  double iPwr(), fact();
  int r, j, k, m, n, c;
  int cterms = costerms(order), sterms = sinterms(order);
  int terms = cterms + sterms;

  /* Allocate the matrix (terms x terms ) */
  CALLOC(mat, terms, double*, ON, AM2M);
  for(r=0; r < terms; r++) CALLOC(mat[r], terms, double, ON, AM2M);
  for(r = 0; r < terms; r++)
      for(c = 0; c < terms; c++) mat[r][c] = 0.0;

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

  /* get the requisite Legendre function evaluations */
  evalLegendre(cosA, tleg, order);

  /* for each new moment (Nj^k) stuff the appropriate matrix entries */
  /* done completely brute force, one term at a time; uses exp in nb 12, p29 */
  for(j = 0; j <= order; j++) {
    for(k = 0; k <= j; k++) {
      for(n = 0, rhoPwr = 1.0; n <= j; n++, rhoPwr *= rho) {
	for(m = 0, mBeta = 0.0; m <= n; m++, mBeta += beta) {

	  if(k == 0) {		/* figure terms for Nj^0, ie k = 0 */
	    if(m <= j-n) {	/* if O moments are nonzero */
	      temp1 = fact(j)*rhoPwr*iPwr(2*m)*tleg[CINDEX(n, m)];
	      temp1 /= (fact(j-n+m)*fact(n+m));
	      mat[CINDEX(j, k)][CINDEX(j-n, m)] += temp1*cos(mBeta);
	      if(m != 0) {		/* if sin term is non-zero */
		mat[CINDEX(j, k)][SINDEX(j-n, m, cterms)] += temp1*sin(mBeta);
	      }
	    }
	  }
	  else {		/* figure terms for Nj^k, k != 0 */
	    temp1 = fact(j+k)*rhoPwr*tleg[CINDEX(n, m)]/fact(n+m);
	    temp2 = temp1*iPwr(2*m)/fact(j-n+k+m);
	    temp1 = temp1*iPwr(k-m-abs(k-m))/fact(j-n+abs(k-m));

	    /* write the cos(kPhi) coeff, bar(N)j^k */
	    if(m != 0) {
	      if(k-m < 0 && abs(k-m) <= j-n) {	/* use conjugates here */
		mat[CINDEX(j, k)][CINDEX(j-n, m-k)] += temp1*cos(mBeta);
		mat[CINDEX(j, k)][SINDEX(j-n, m-k,cterms)] += temp1*sin(mBeta);
	      }
	      else if(k-m == 0) {	/* double to compensate for 2Re sub. */
		mat[CINDEX(j, k)][CINDEX(j-n, k-m)] += 2*temp1*cos(mBeta);
		/* sin term is always zero */
	      }
	      else if(k-m > 0 && k-m <= j-n) {
		mat[CINDEX(j, k)][CINDEX(j-n, k-m)] += temp1*cos(mBeta);
		mat[CINDEX(j, k)][SINDEX(j-n, k-m,cterms)] -= temp1*sin(mBeta);
	      }
	      if(k+m <= j-n) {
		mat[CINDEX(j, k)][CINDEX(j-n, k+m)] += temp2*cos(mBeta);
		mat[CINDEX(j, k)][SINDEX(j-n, k+m,cterms)] += temp2*sin(mBeta);
	      }
	    }			/* do if m = 0 and O moments not zero */
	    else if(k <= j-n) mat[CINDEX(j, k)][CINDEX(j-n, k)] += temp2;

	    /* write the sin(kPhi) coeff, dblbar(N)j^k, if it is non-zero */
	    if(m != 0) {
	      if(k-m < 0 && abs(k-m) <= j-n) {	/* use conjugates here */
		mat[SINDEX(j, k,cterms)][CINDEX(j-n, m-k)] += temp1*sin(mBeta);
		mat[SINDEX(j, k, cterms)][SINDEX(j-n, m-k, cterms)] 
		    -= temp1*cos(mBeta);
	      }
	      else if(k-m == 0) {/* double to compensate for 2Re sub */
		mat[SINDEX(j, k, cterms)][CINDEX(j-n, k-m)] 
		    += 2*temp1*sin(mBeta);
		/* sine term is always zero */
	      }
	      else if(k-m > 0 && k-m <= j-n) {
		mat[SINDEX(j, k,cterms)][CINDEX(j-n, k-m)] += temp1*sin(mBeta);
		mat[SINDEX(j, k, cterms)][SINDEX(j-n, k-m, cterms)] 
		    += temp1*cos(mBeta);
	      }
	      if(k+m <= j-n) {
		mat[SINDEX(j, k,cterms)][CINDEX(j-n, k+m)] -= temp2*sin(mBeta);
		mat[SINDEX(j, k, cterms)][SINDEX(j-n, k+m, cterms)] 
		    += temp2*cos(mBeta);
	      }
	    }			/* do if m = 0 and moments not zero */
	    else if(k <= j-n) mat[SINDEX(j, k,cterms)][SINDEX(j-n, k, cterms)] 
		+= temp2;
	  }
	}
      }
    }
  }
#if DISM2M == ON
  dispM2M(mat, x, y, z, xp, yp, zp, order);
#endif
  return(mat);
}

/* 
  builds multipole evaluation matrix; used only for fake downward pass 
*/
double **mulMulti2P(x, y, z, chgs, numchgs, order)
double x, y, z;			/* multipole expansion origin */
charge **chgs;
int numchgs, order;
{
  double **mat;
  double cosTh;			/* cosine of elevation coordinate */
  double factorial;		/* 1/factorial = (n-m)!/(n+m)! */
  int i, j, k, m, n, kold, start;
  int cterms = costerms(order), sterms = sinterms(order);
  int terms = cterms + sterms;

  CALLOC(mat, numchgs, double*, ON, AM2P);
  for(i=0; i < numchgs; i++) 
      CALLOC(mat[i], terms, double, ON, AM2P);

  /* get Legendre function evaluations, one set for each charge */
  /*   also get charge coordinates to set up rest of matrix */
  for(i = 0; i < numchgs; i++) { /* for each charge, do a legendre eval set */
    xyz2sphere(chgs[i]->x, chgs[i]->y, chgs[i]->z,
	       x, y, z, &(Ir[i]), &cosTh, &(phi[i]));

    Irn[i] = Ir[i]; /* initialize (1/r)^n+1 vec. */
    Mphi[i] = phi[i];		/* initialize m*phi vector */

    evalLegendre(cosTh, mat[i], order);	/* wr moms to 1st (cos) half of row */

  }

#if DALM2P == ON
  fprintf(stdout,
	  "\nM2P MATRIX BUILD:\n    AFTER LEGENDRE FUNCTION EVALUATON\n");
  dumpMat(mat, numchgs, terms);
#endif

  /* add the (1/r)^n+1 factors to the left (cos(m*phi)) half of the matrix */
  for(j = 0, k = kold = 1; j < cterms; j++) { /* loop on columns of matrix */
    for(i = 0; i < numchgs; i++) mat[i][j] /= Irn[i]; /* divide by r^n+1 */
    k -= 1;
    if(k == 0) {		/* so that n changes as appropriate */
      kold = k = kold + 1;
      for(i = 0; i < numchgs; i++) Irn[i] *= Ir[i]; /* r^n -> r^n+1 */
    }
  }

#if DALM2P == ON
  fprintf(stdout,
	  "    AFTER ADDITION OF (1/R)^N+1 FACTORS\n");
  dumpMat(mat, numchgs, terms);
#endif

  /* add the factorial fraction factors to the left (cos(m*phi)) part of mat */
  /*  note that (n-m)!/(n+m)! = 1/(n-m+1)...(n+m) since m \leq n */
  for(n = 1; n <= order; n++) {
    for(m = 1; m <= n; m++) {
      for(i = 0; i < numchgs; i++) mat[i][CINDEX(n, m)] /= factFac[n][m];
    }
  }

#if DALM2P == ON
  fprintf(stdout,
	  "    AFTER ADDITION OF FACTORIAL FRACTION FACTORS\n");
  dumpMat(mat, numchgs, terms);
#endif

  /* copy left half of matrix to right half for sin(m*phi) terms */
  for(i = 0; i < numchgs; i++) { /* loop on rows of matrix */
    for(n = 1; n <= order; n++) { 
      for(m = 1; m <= n; m++) {	/* copy a row */
	mat[i][SINDEX(n, m, cterms)] = mat[i][CINDEX(n, m)];
      }
    }
  }

#if DALM2P == ON
  fprintf(stdout,
	  "    AFTER COPYING SINE (RIGHT) HALF\n");
  dumpMat(mat, numchgs, terms);
#endif

  /* add factors of cos(m*phi) and sin(m*phi) to left and right halves resp. */
  for(m = 1; m <= order; m++) {	/* lp on m in Mn^m (no m=0 since cos(0)=1) */
    for(n = m; n <= order; n++) { /* loop over cols with same m */
      for(i = 0; i < numchgs; i++) { /* add factors to a column */
	mat[i][CINDEX(n, m)] *= cos(Mphi[i]);
	mat[i][SINDEX(n, m, cterms)] *= sin(Mphi[i]);
      }
    }
    for(i = 0; i < numchgs; i++) Mphi[i] += phi[i]; /* (m-1)*phi->m*phi */
  }

#if DISM2P == ON
  dispM2P(mat, x, y, z, chgs, numchgs, order);
#endif

  return(mat);
}