calcp.c

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

  result_vector[XI] = vector1[YI]*vector2[ZI] - vector1[ZI]*vector2[YI];
  result_vector[YI] = vector1[ZI]*vector2[XI] - vector1[XI]*vector2[ZI];
  result_vector[ZI] = vector1[XI]*vector2[YI] - vector1[YI]*vector2[XI];
}

/*
Computes the potential at x, y, z due to a unit source on panel
and due to a dipole is returned in pfd.

Note, the code is subtle because there are 5 cases depending on
the placement of the collocation point

    CASE1: evaluation point projection strictly outside of panel (happens most
      of the time).
    CASE2: eval pnt proj. strictly inside panel (happens >= #panels times,
      at least for the self terms)
    CASE3: eval pnt proj. on a panel side, not on a corner (happens
      when paneled faces meet at right angles, also possible other ways)
    CASE4: eval pnt proj. on a panel corner (happens rarely to never)
    CASE5: eval pnt proj. on side extension (happens when paneled 
      faces meet at right angles, also possible other ways).
*/
double calcp(panel, x, y, z, pfd)
charge *panel;
double x, y, z, *pfd;
{
  double r[4], fe[4], xmxv[4], ymyv[4];
  double xc, yc, zc, zsq, xn, yn, zn, znabs, xsq, ysq, rsq, diagsq, dtol;
  double v, arg, st, ct, length, s1, c1, s2, c2, s12, c12, val;
  double *s;
  double rInv, r2Inv, r3Inv, r5Inv, r7Inv, r9Inv, zr2Inv;
  double ss1, ss3, ss5, ss7, ss9;
  double s914, s813, s411, s512, s1215;
  double fs, fd, fdsum;
  int okay, i, next;
  struct edge *edge;
  double *corner;

  /* Put the evaluation point into this panel's coordinates. */
  xc = x - panel->x;
  yc = y - panel->y;
  zc = z - panel->z;

  xn = DotP_Product(panel->X, xc, yc, zc);
  yn = DotP_Product(panel->Y, xc, yc, zc);
  zn = DotP_Product(panel->Z, xc, yc, zc);

  zsq = zn * zn;
  xsq = xn * xn;
  ysq = yn * yn;
  rsq = zsq + xsq + ysq;
  diagsq = panel->max_diag * panel->max_diag;

  if(rsq > (LIMITFOURTH * diagsq)) {
    fs = 0.0; fd = 0.0;
    s = panel->moments;
    /* First, second moments. */
    r2Inv = 1.0 / rsq;
    rInv = sqrt(r2Inv);
    r3Inv = r2Inv * rInv;
    r5Inv = r3Inv * r2Inv;
    zr2Inv = zn * r2Inv;
    ss1 = s[1] * rInv;
    ss3 = -(s[3] + s[10]) * r3Inv;
    ss5 = (xsq * s[10] + (xn * yn * s[7]) + ysq * s[3]) * r5Inv;
    fs = ss1 + ONE3 * ss3 + ss5;
    fdsum = ss1 + ss3 + 5.0 * ss5;
    fd = zr2Inv * fdsum;
    if(rsq < (LIMITSECOND * diagsq)) {
    /* Third and fourth moments added for diagsq/r2 between 40 and 150. */
      s914 = s[9] + s[14];
      s813 = s[8] + s[13];
      s411 = s[4] + s[11];
      s512 = s[5] + s[12];
      s1215 = s[12] + s[15];
      r7Inv = r5Inv * r2Inv;
      r9Inv = r7Inv * r2Inv;
      ss5 = (-xn * s813 - yn * s411 + 0.1 * (s512 + s1215)) * r5Inv;

      ss7 = (FIVE3 *((xn * xsq * s[13] + yn * ysq * s[4]) 
		     + 3.0 * xn * yn * (xn * s[11]  +  yn * s[8]))
		     - xsq * s1215 - ysq * s512 - xn * yn * s914) * r7Inv;

      ss9 = (7.0 * (ONE6 * (xsq * xsq * s[15] + ysq * ysq * s[5])
		    + xsq * ysq * s[12])
	     + SEVEN3 * xn * yn * (xsq * s[14] + ysq * s[9])) * r9Inv;

      fs += ss5 + ss7 + ss9;
      fdsum = 5.0 * ss5 + 7.0 * ss7 + 9.0 * ss9;
      fd += zr2Inv * fdsum;
      num4th++;
    }
    else num2nd++;
  }
  else {
    dtol = EQUIV_TOL * panel->min_diag;
    znabs = fabs(zn);

    /* Always move the evaluation point a little bit off the panel. */
    if(znabs < dtol) { 
      zn = 0.5 * dtol;  /* Half of dtol insures detection for zero dipole. */
      znabs = 0.5 * dtol;
    }

    /* Once per corner computations. */
    for(okay = TRUE, i=0; i < panel->shape; i++) {
      corner = panel->corner[i]; 
      xmxv[i] = xc = xn - corner[XI];
      ymyv[i] = yc = yn - corner[YI];
      zc = zn - corner[ZI];
      fe[i] = xc * xc + zc * zc;
      r[i] = sqrt(yc * yc + fe[i]);
      if(r[i] < (1.005 * znabs)) {  /* If r almost z, on vertex normal. */
	okay = FALSE;
      }
    }

    /* Once per edge computations. */
    fs = 0.0; fd = 0.0;
    for(i=0; i < panel->shape; i++) {
      if(i == (panel->shape - 1)) next = 0;
      else next = i + 1;

      /* Now calculate the edge contributions to a panel. */
      length = panel->length[i];
      ct = (panel->corner[next][XI] - panel->corner[i][XI]) / length;
      st = (panel->corner[next][YI] - panel->corner[i][YI]) / length;

      /* v is projection of eval-i edge onto perpend to next-i edge. */
      /* Exploits the fact that corner points in panel coordinates. */
      v = xmxv[i] * st - ymyv[i] * ct;

      /* arg == zero if eval on next-i edge, but then v = 0. */
      arg = (r[i] + r[next] - length)/(r[i] + r[next] + length);
      if(arg > 0.0) fs -= v * log(arg);

      /* Okay means eval not near a vertex normal, Use Hess-Smith. */
      if(okay) {
	s1 = v * r[i];
	c1 = znabs * (xmxv[i] * ct + ymyv[i] * st);
	s2 = v * r[next];
	c2 = znabs * (xmxv[next] * ct + ymyv[next] * st);
      }
      /* Near a vertex normal, use Newman. */
      else {
	s1 = (fe[i] * st) - (xmxv[i] * ymyv[i] * ct);
	c1 = znabs * r[i] * ct;
	s2 = (fe[next] * st) - (xmxv[next] * ymyv[next] * ct);
	c2 = znabs * r[next] * ct;
      }    

      s12 = (s1 * c2) - (s2 * c1);
      c12 = (c1 * c2) + (s1 * s2);
      val = atan2(s12, c12);
      fd += val;
    }
    /* Adjust the computed values. */

    if(fd < 0.0) fd += TWOPI;
    if(zn < 0.0) fd *= -1.0;
    if(znabs < dtol) fd = 0.0;

    fs -= zn * fd;
    numexact++;
  }

  /* Return values of the source and dipole, normalized by area. */
  fs /= panel->area;
  fd /= panel->area;
  if(pfd != NULL) *pfd = fd;


  if(fs < 0.0) {
    fprintf(stderr, 
	    "\ncalcp: panel potential coeff. less than zero = %g\n", fs);
    fprintf(stderr, "Okay = %d Evaluation Point = %g %g %g\n", okay, x, y, z);
    fprintf(stderr, "Evaluation Point in local coords = %g %g %g\n",xn,yn, zn);
    fprintf(stderr, "Panel Description Follows\n");
    dp(panel);
    /*exit(1);*/
  }


  return (fs);
}


dumpnums(flag, size)
int flag, size;
{
  double total;

  if(flag == ON) {		/* if first call */
    num2ndsav = num2nd;
    num4thsav = num4th;
    numexactsav = numexact;
  }
  else {
    total = num2ndsav + num4thsav + numexactsav;
#if MULDAT == ON
    fprintf(stdout, "Potential coefficient counts\n multipole only:\n");
    fprintf(stdout, 
	    "  2nd order: %d %.3g%%; 4th: %d %.3g%%; Integral: %d %.3g%%\n",
	    num2nd, 100*(num2ndsav/total), num4th, 100*(num4thsav/total), 
	    numexact, 100*(numexactsav/total));
#endif
    total = num2nd + num4th + numexact;
#if MULDAT == ON
    fprintf(stdout, " multipole plus adaptive:\n");
    fprintf(stdout, 
	    "  2nd order: %d %.3g%%; 4th: %d %.3g%%; Integral: %d %.3g%%\n",
	    num2nd, 100*(num2nd/total), num4th, 100*(num4th/total), 
	    numexact, 100*(numexact/total));
#endif
    fprintf(stdout, "Percentage of multiplies done by multipole: %.3g%%\n",
	    100*(size*size - total)/(size*size));
    if(size*size == total) 
	fprintf(stdout, "Warning: no multipole acceleration\n");
  }
}

double tilelength(nq)
charge *nq;
{
  return nq->max_diag;
}



/*
Evaluation of moments of quadrilateral surface relative to
local system, array S(15).  First initialize array
Note that S(2)=S(6)=0 due to transfer above
*/

ComputeMoments(pp)
charge *pp;
{
  int order=MAXORDER;
  int i, j, nside,  N, M, N1, M1, M2, MN1, MN2;
  double dx, dy, dxdy, dydx, x, y, z, SI, *xp, *yp, *xpn, *ypn;
  double ypp[3];
  static double *XP[4], *YP[4], **I;
  static int maxorder = 0;
  static double CS[16] = { 0.0, 1.0, 1.0, 1.5, 1.5, 3.75, 1.0, 3.0, 
			   1.5, 7.5, 1.5, 1.5, 3.75, 1.5, 7.5, 3.75 };
  double *multi, sumc, sums, sign, **createBinom(), **createPmn();
  int m, n, r, halfn, flrm, ceilm, numterms, rterms;
  
  /* Allocate temporary storage and initialize arrays. */
  if(order > maxorder) {
    for(i = 0; i < 4; i++) {
      CALLOC(XP[i], (order+3), double, ON, AQ2P);
      CALLOC(YP[i], (order+3), double, ON, AQ2P);
    }
    /* Allocate the euclidean moments matrix, Imn. */
    CALLOC(I, (order+1), double*, ON, AQ2P);
    for(i = 0; i <= order; i++) CALLOC(I[i], (order+1), double, ON, AQ2P);

    maxorder = order;
  }

  /* First zero out the Moments matrix. */
  for(i = 0; i <= order; i++) {
    for(j = 0; j <= order; j++) {
      I[i][j] = 0.0;
    }
  }
    
  /* Compute powers of x and y at corner pts. */
  for(i = 0; i < pp->shape; i++) {
    xp = XP[i];
    yp = YP[i];
    xp[1] = pp->corner[i][XI];
    yp[1] = pp->corner[i][YI];
    for(j = 2; j <= order+2; j++) {
      xp[j] = xp[j-1] * xp[1];
      yp[j] = yp[j-1] * yp[1];
    }
  }

  /* First moment, easy, just the panel area. */
  I[0][0] = pp->area;

  /* By using centroid, (1,0) and (0,1) are zero, so begin with (2,0). */
  for(nside = 0; nside < pp->shape; nside++) {
    xp = XP[nside];
    yp = YP[nside];
    if(nside == (pp->shape - 1)) {
      xpn = XP[0];
      ypn = YP[0];  
    }
    else {
      xpn = XP[nside + 1];
      ypn = YP[nside + 1];
    }

    dx = xpn[1] - xp[1];
    dy = ypn[1] - yp[1];

    if(fabs(dx) >= fabs(dy)) {
      dydx = dy/dx;
      for(M = 2; M <= order; M++) {
	M1 = M + 1;
	M2 = M + 2;

	SI = ((xpn[M1] * ypn[1]) - (xp[M1] * yp[1])) / M1
	     + dydx * (xp[M2] - xpn[M2]) / (M1 * M2);
	I[M][0] += SI;

	for(N = 1; N <= M; N++) {
	  N1 = N + 1;
	  MN1 = M - N + 1;
	  SI = (xpn[MN1] * ypn[N1] - xp[MN1] * yp[N1]) / (MN1 * N1)
		     - (dydx * N * SI) / MN1;
	  I[M-N][N] += SI;
	}
      }
    }
    else {
      dxdy = dx/dy;
      for(M = 2; M <= order; M++) {
	M1 = M + 1;
	M2 = M + 2;
	SI = (dxdy / (M1 * M2)) * (ypn[M2] - yp[M2]);
	I[0][M] += SI;
	for(N = 1; N <= M; N++) {
	  MN1 = M - N + 1;
	  MN2 = MN1 + 1;
          SI = dxdy * ((xpn[N] * ypn[MN2] - xp[N] * yp[MN2]) / (MN1 * MN2) 
			- (N * SI / MN1));
	  I[N][M-N] += SI;
	}
      }
    }
  }

  /* Now Create the S vector for calcp. */
  for(i = 0, M = 0; M <= 4; M++) {
    for(N = 0; N <= (4 - M); N++) {
      i++;
      pp->moments[i] = I[M][N] * CS[i];
    }
  }

}

/* Debugging Print Routines follow. */

dp(panel)
charge *panel;
{
  int i;
  double c[4][3];

  printf("shape=%d maxdiag=%g mindiag=%g area=%g\n", 
	 panel->shape, 
	 panel->max_diag, panel->min_diag, panel->area);

  fprintf(stdout, "surface: `%s'\n", panel->surf->name);

  printf("x=%g y=%g z=%g\n", panel->x, panel->y, panel->z);
  printf("X= %g %g %g\n", panel->X[0], panel->X[1], panel->X[2]);
  printf("Y= %g %g %g\n", panel->Y[0], panel->Y[1], panel->Y[2]);
  printf("Z= %g %g %g\n", panel->Z[0], panel->Z[1], panel->Z[2]);

  for(i=0; i < panel->shape; i++)
      printf("corner%d = %g %g %g\n", 
	     i, panel->corner[i][0], panel->corner[i][1], panel->corner[i][2]);

  for(i = 0; i < panel->shape; i++) {
    c[i][0] = panel->x + panel->corner[i][0]*panel->X[0] 
	+ panel->corner[i][0]*panel->X[1] + panel->corner[i][0]*panel->X[2];
    c[i][1] = panel->y + panel->corner[i][1]*panel->Y[0] 
	+ panel->corner[i][1]*panel->Y[1] + panel->corner[i][1]*panel->Y[2];
    c[i][2] = panel->z + panel->corner[i][2]*panel->Z[0] 
	+ panel->corner[i][2]*panel->Z[1] + panel->corner[i][2]*panel->Z[2];
    printf("absolute corner%d = %g %g %g\n", i, c[i][0], c[i][1], c[i][2]);
  }

  for(i=0; i < panel->shape; i++)
      printf("length%d = %g\n", i, panel->length[i]);

  printf("multipole coeffs:  ");
  for(i=0; i < 16; i++) {
    printf("%g  ", panel->moments[i]);
    if( (i % 6) == 0) printf("\n");
  }
  printf("\n");
}



#define DIS 2
#define SCALE 5

testCalcp(pp)
charge *pp;
{

  double offx, offy, offz, x, y, z, mult;
  int i, j, k;

  offx = pp->x;
  offy = pp->y;
  offz = pp->z;

  mult = 0.5 * pp->min_diag;

  printf("\n\nCenter Point %g %g %g\n", offx, offy, offz);
  for(i=0; i < DIS; i++) {
    for(j=0; j < DIS; j++) {
      for(k=0; k < DIS; k++) {
	x = offx + i * mult * SCALE;
	y = offy + j * mult * SCALE;
	z = offz + k * mult * SCALE;
	printf("Eval pt = %g %g %g\n", x, y, z);
        calcp(pp, x, y, z, NULL);
      }
    }
  }
}


fileCorners(pp, f)
FILE *f;
charge *pp;

{
  int i;

  for(i=0; i < pp->shape; i++)
      fprintf(f, "%g %g\n", pp->corner[i][0], pp->corner[i][1]);
}


/* Test the moment code. */
calcpm(multi, x, y, z, origorder, order)
double *multi, x, y, z;
int order;
{
  charge panel, *ppanel;
  double **mat, **mulMulti2P(), potential;
  int i, numterms;

  /* Create a temporary panel which has evaluation point as centroid. */
  panel.x = x;
  panel.y = y;
  panel.z = z;
  ppanel = &panel;

  mat = mulMulti2P(0.0, 0.0, 0.0, &ppanel, 1, origorder);
  numterms = multerms(origorder);

  /* Calculate the potential. */
  for(potential = 0.0, i = 0; i < numterms; i++) {
    potential += mat[0][i] * multi[i];
  }
}



/*

dismat(mat, size)
double **mat;
int size;
{
  int i, j;

  for(i=0; i <= size; i++) {
    printf("\nrow i = ");
    for(j=0; j <= size; j++) {
      printf("%g  ", mat[i][j]);
    }
  }
  printf("\n");
}
*/