calcp.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"

#ifndef FALSE
#define FALSE 0
#endif
#ifndef TRUE
#define TRUE 1
#endif
#define XI 0
#define YI 1
#define ZI 2
#define EQUIV_TOL 1.0e-9
/*#define PLANAR_TOL 1.0e-3*/
#define PLANAR_TOL 1.0e-5
#define MAX_REAL 1.0e+20;

/* Obvious Constants. */
#define PI 3.1415927
#define TWOPI 6.2831853
#define HALFPI 1.5707963E+00

/* Constants to save typing. */
#define FIVE3 1.666666666667
#define SEVEN3 2.3333333333333
#define ONE6 0.16666666666667
#define ONE3 0.3333333333333
#define FT3 4.66666667


/* Defines breakpoints panel multipoles. */
#define LIMITFOURTH 9.0
#define LIMITSECOND 36.0

/* Constants used for Hart arctan approximation. */
#define B1 0.24091197
#define B2 3.7851122
#define B3 5.6770721
#define B4 5.6772854
#define B5 5.6770747

#define Dot_Product(V1,V2) (V1[XI]*V2[XI]+V1[YI]*V2[YI]+V1[ZI]*V2[ZI])
#define DotP_Product(V1,R,S,T) ((V1[XI])*(R)+(V1[YI])*(S)+(V1[ZI])*(T))


static int num2nd=0, num4th=0, numexact=0;
static int num2ndsav=0, num4thsav=0, numexactsav=0;

initcalcp(panel_list)
  charge *panel_list;
{
  charge *pq, *npq;
  double vtemp[3];
  double length, maxlength, minlength, length20, length31, sum, sum2, delta;
  double normalize();
  int i, j, next;

#if JACDBG == ON
  FILE *foo;

  foo = fopen("corners.txt", "w");
#endif

  for(i=0, pq = panel_list; pq != NULL; pq = pq->next) i++;
#if JACDBG == ON
  printf("Initial number of panels = %d\n", i);
#endif

  pq = panel_list; 
  while (pq != NULL) {
    /* Calculate edge lengths. */
    maxlength = 0.0;
    minlength = MAX_REAL;
    for(i=0; i < pq->shape; i++) {    
      if(i == (pq->shape -1)) next = 0;
      else next = i + 1;
      for(sum= 0, j = 0; j < 3; j++) {
	delta = pq->corner[next][j] - pq->corner[i][j];
	sum += delta * delta;
      }
      pq->length[i] = length = sqrt(sum);
      maxlength = MAX(maxlength, length);
      minlength = MIN(minlength, length);
    }
    
    /* Get diags and lengths. */
    for(sum= 0, sum2 = 0, i = 0; i < 3; i++) {     
      pq->X[i] = delta = pq->corner[2][i] - pq->corner[0][i];
      sum += delta * delta;
      if(pq->shape == 3) pq->Y[i] = pq->corner[0][i] - pq->corner[1][i];      
      else {
	pq->Y[i] = delta = pq->corner[3][i] - pq->corner[1][i];      
	sum2 += delta * delta;
      }
    }
    length20 = sqrt(sum);
    length31 = sqrt(sum2);

    /* Check on lengths for quad. */
    if(pq->shape == 3) {
      pq->max_diag = maxlength;
      pq->min_diag = minlength;
    }
    else {
      length = MAX(length20, length31);
      pq->max_diag = length;
      pq->min_diag = MIN(length20, length31);
    }

    /* Z-axis is normal to two diags. */
    Cross_Product(pq->X, pq->Y, pq->Z);
/*#if 1 == 0*/
    if(flip_normal(pq)) {
      for(i = 0; i < 3; i++) {
	pq->Z[i] = -(pq->Z[i]);	/* flip the normal */
	pq->X[i] = -(pq->X[i]);	/* flip the x direction */
	/* interchange points 0 and 2 so that corner order will be
	   consistent with X flip (note this is OK for both quads and tris) */
	vtemp[0] = pq->corner[0][i];
	pq->corner[0][i] = pq->corner[2][i];
	pq->corner[2][i] = vtemp[0];
      }
      /* interchange length entries in length vector */
      vtemp[0] = pq->length[0];
      pq->length[0] = pq->length[1];
      pq->length[1] = vtemp[0];
      if(pq->shape == 4) {
	vtemp[0] = pq->length[2];
	pq->length[2] = pq->length[3];
	pq->length[3] = vtemp[0];
      }
    } 
/*#endif*/
    pq->area = 0.5 * normalize(pq->Z);
    normalize(pq->X);

    /* Real Y-axis is normal to X and Z (resulting system is left-handed). */
    Cross_Product(pq->X, pq->Z, pq->Y);

    /* Project the corner points into the plane defined by edge midpoints. */
    if(planarize(pq) == FALSE) {     
      /* Planarization too drastic, crack into two triangles. */
      CALLOC(npq, 1, charge, ON, AMSC);
      npq->next = pq->next;
      pq->next = npq;
      npq->shape = 3;
      pq->shape = 3;
      npq->cond = pq->cond;
      npq->surf = pq->surf;
      VCOPY(npq->corner[0], pq->corner[0]);
      VCOPY(npq->corner[1], pq->corner[2]);
      VCOPY(npq->corner[2], pq->corner[3]);
      continue;
    }

    /* Calculate the centroid. */
    centroid(pq, length20);      

    /* Put corners in the newly defined coord system. */
    for(i=0; i < pq->shape; i++) {
      pq->corner[i][XI] -= pq->x;
      pq->corner[i][YI] -= pq->y;
      pq->corner[i][ZI] -= pq->z;
    }
    for(i=0; i < pq->shape; i++) {
      vtemp[XI] = Dot_Product(pq->corner[i], pq->X);
      vtemp[YI] = Dot_Product(pq->corner[i], pq->Y);
      vtemp[ZI] = Dot_Product(pq->corner[i], pq->Z);
      pq->corner[i][XI] = vtemp[XI];
      pq->corner[i][YI] = vtemp[YI];
      pq->corner[i][ZI] = vtemp[ZI];
      if(fabs(pq->corner[i][ZI]) > (EQUIV_TOL * pq->min_diag)) {
	printf("FATAL PROGRAM ERROR: renormalized z=%g\n", pq->corner[i][ZI]);
	exit(1);
      }
      pq->corner[i][ZI] = 0.0;
    }

#if JACDBG == ON			/* dump corners, center to file */
    fprintf(foo, "%g %g %g: ", pq->x, pq->y, pq->z);
    for(i = 0; i < pq->shape; i++) {
      fprintf(foo, "%g %g %g ", pq->corner[i][XI], pq->corner[i][YI],
	      pq->corner[i][ZI]);
    }
    fprintf(foo, "\n");
#endif

    /* Finally, calculate the moments. */
    ComputeMoments(pq);

    /* Iterate for the next panel. */
    pq = pq->next;

  }

}

/*
  determine if normal needs to be flipped to get dielectric bdry cond right
*/
oldflip_normal(panel)
charge *panel;
{
  int i;
  double x, y, z;
  double ctr_minus_n[3], ctr_plus_n[3], norm_minus, norm_plus, norm, norm_sq;
  surface *surf = panel->surf;
  int ref_inside = surf->ref_inside, flip_normal;
  double *ref = surf->ref, *normal;
  char *surf_name = surf->name;

  if(surf->type != DIELEC && surf->type != BOTH) return(FALSE);

  /* get panel corner (relative to reference point) and normal */
  x = panel->corner[0][0] - ref[0]; 
  y = panel->corner[0][1] - ref[1]; 
  z = panel->corner[0][2] - ref[2];
  norm_sq = x*x + y*y + z*z;
  norm = sqrt(norm_sq);
  normal = panel->Z;

  /* add the (scaled) normal and negative normal to the panel corner */
  /* negative normal result should be closer to ref point if(ref_inside) */
  ctr_minus_n[0] = x - 0.1*norm*normal[0];
  ctr_minus_n[1] = y - 0.1*norm*normal[1];
  ctr_minus_n[2] = z - 0.1*norm*normal[2];
  ctr_plus_n[0] = x + 0.1*norm*normal[0];
  ctr_plus_n[1] = y + 0.1*norm*normal[1];
  ctr_plus_n[2] = z + 0.1*norm*normal[2];

  /* get norms of test points, one inside (minus) other out (plus) */
  norm_minus = ctr_minus_n[0]*ctr_minus_n[0];
  norm_plus = ctr_plus_n[0]*ctr_plus_n[0];
  for(i = 1; i < 3; i++) {
    norm_minus += ctr_minus_n[i]*ctr_minus_n[i];
    norm_plus += ctr_plus_n[i]*ctr_plus_n[i];
  }

  flip_normal = FALSE;
  if(norm_minus > norm_sq) {
    if(norm_plus > norm_sq) {
      fprintf(stderr, 
	      "flip_normal: both test points on non-reference side\n");
      fprintf(stderr, "  Surface: %s\n", hack_path(surf_name));
      fprintf(stderr, "  Translation: (%g %g %g)\n", surf->trans[0],
	      surf->trans[1], surf->trans[2]);
      fprintf(stderr, "  Reference point: (%g %g %g)\n",
	      ref[0], ref[1], ref[2]);
      fprintf(stderr, "  Panel corner: (%g %g %g)\n",
	      panel->corner[0][0], panel->corner[0][1], panel->corner[0][2]);
      fprintf(stderr, "  Normal: (%g %g %g)\n",
	      normal[0], normal[1], normal[2]);
      exit(1);
    }
    if(ref_inside) flip_normal = TRUE;
    }
  else if(norm_plus < norm_sq) {
    if(norm_minus < norm_sq) {
      fprintf(stderr, 
	      "flip_normal: both test points on reference point side\n");
      fprintf(stderr, "  Surface: %s\n", hack_path(surf_name));
      fprintf(stderr, "  Translation: (%g %g %g)\n", surf->trans[0],
	      surf->trans[1], surf->trans[2]);
      fprintf(stderr, "  Reference point: (%g %g %g)\n",
	      ref[0], ref[1], ref[2]);
      fprintf(stderr, "  Panel corner: (%g %g %g)\n",
	      panel->corner[0][0], panel->corner[0][1], panel->corner[0][2]);
      fprintf(stderr, "  Normal: (%g %g %g)\n",
	      normal[0], normal[1], normal[2]);
      exit(1);
    }	
    if(!ref_inside) flip_normal = TRUE;
  }

  return(flip_normal);
}

/*
  determine if normal needs to be flipped to get dielectric bdry cond right
  - this function uses 0.0 as a breakpoint when really machine precision
    weighted checks should be done (really not an issue if ref point far)
*/
flip_normal(panel)
charge *panel;
{
  int i;
  double x, y, z;
  double ctr_minus_n[3], ctr_plus_n[3], norm_minus, norm_plus, norm, norm_sq;
  surface *surf = panel->surf;
  int ref_inside = surf->ref_inside, flip_normal;
  double *ref = surf->ref, *normal, angle, norm_n;
  char *surf_name = surf->name;

  if(surf->type != DIELEC && surf->type != BOTH) return(FALSE);

  /* get panel corner (relative to reference point) and normal */
  x = panel->corner[0][0] - ref[0]; 
  y = panel->corner[0][1] - ref[1]; 
  z = panel->corner[0][2] - ref[2];
  norm_sq = x*x + y*y + z*z;
  norm = sqrt(norm_sq);
  normal = panel->Z;
  norm_n = 
      sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);

  /* figure the dot product between the normal and the corner-ref point
     - ref_inside and angle <= 90 => flip
     - ref_inside and angle  > 90 => no flip
     - !ref_inside and angle <= 90 => no flip 
     - !ref_inside and angle > 90 => flip */
  /*    figure angle */
  angle = (x*normal[0] + y*normal[1] + z*normal[2])/(norm*norm_n);
  if(ref_inside && angle <= 0.0) flip_normal = TRUE;
  else if(ref_inside && angle > 0.0) flip_normal = FALSE;
  else if(!ref_inside && angle <= 0.0) flip_normal = FALSE;
  else if(!ref_inside && angle > 0.0) flip_normal = TRUE;
  else {
    fprintf(stderr, 
	    "flip_normal: inconclusive test for normal flipping\n");
    fprintf(stderr, "  Surface: %s\n", hack_path(surf_name));
    fprintf(stderr, "  Translation: (%g %g %g)\n", surf->trans[0],
	    surf->trans[1], surf->trans[2]);
    fprintf(stderr, "  Reference point: (%g %g %g)\n",
	    ref[0], ref[1], ref[2]);
    fprintf(stderr, "  Panel corner: (%g %g %g)\n",
	    panel->corner[0][0], panel->corner[0][1], panel->corner[0][2]);
    fprintf(stderr, "  Normal: (%g %g %g)\n",
	    normal[0], normal[1], normal[2]);
    exit(1);
  }
    
  return(flip_normal);
}
  

/*
Changes the corner points so that they lie in the plane defined by the
panel diagonals and any midpoint of an edge.
*/
planarize(pq)
charge *pq;
{
  double origin[3], corner[3], delta[4][3], px, py, dx, dy, dz;
  int i, j, numcorners = pq->shape;
  double tolsq = PLANAR_TOL * pq->min_diag; 

  tolsq *= tolsq;

  /* Triangular panels are just fine already. */
  if(numcorners != 4) return(TRUE);

  /* Pick edge midpoint as origin. */
  for(i=0; i < 3; i++) origin[i] = 0.5 * (pq->corner[1][i] + pq->corner[0][i]);

  for(i=0; i < numcorners; i++) {
    for(j=0; j < 3; j++) corner[j] = pq->corner[i][j] - origin[j];
    px = Dot_Product(corner, pq->X);
    py = Dot_Product(corner, pq->Y);

    dx = px * pq->X[XI] + py * pq->Y[XI] + origin[XI] - pq->corner[i][XI];
    dy = px * pq->X[YI] + py * pq->Y[YI] + origin[YI] - pq->corner[i][YI];
    dz = px * pq->X[ZI] + py * pq->Y[ZI] + origin[ZI] - pq->corner[i][ZI];

    delta[i][XI] = dx;
    delta[i][YI] = dy;
    delta[i][ZI] = dz;
    
    /* Moved too much, crack the panel. */
    if((dx * dx + dy * dy + dz * dz) > tolsq) return(FALSE);
  }

  /* Still Here? Must be planarizing. */
  for(i=0; i < numcorners; i++) {
    for(j=0; j < 3; j++) {
      pq->corner[i][j] += delta[i][j];
    }
  }
  return(TRUE);
}


/* 
Determines centroid of a panel (defined as the point which make the
first moments vanish.  Calculation begins by projection into the
coordinate system defined by the panel normal as the z-axis and
edge02 as the x-axis.
*/
centroid(pp, x2)
charge *pp;
double x2;
{
  double vertex1[3], vertex3[3];
  double sum, dl, x1, y1, x3, y3, xc, yc;
  int i;

  /* Use vertex 0 as the origin. */
  for(i=0; i< 3; i++) {
    vertex1[i] = pp->corner[1][i] - pp->corner[0][i];
    if(pp->shape == 4) vertex3[i] = pp->corner[3][i] - pp->corner[0][i];
    else vertex3[i] = pp->corner[2][i] - pp->corner[0][i];
  }

  /* Project into the panel axes. */
  y1 = Dot_Product(vertex1, pp->Y);
  y3 = Dot_Product(vertex3, pp->Y);
  x1 = Dot_Product(vertex1, pp->X);
  x3 = Dot_Product(vertex3, pp->X);

  yc = ONE3 * (y1 + y3);
  xc = ONE3 * (x2 + ((x1 * y1 - x3 * y3)/(y1 - y3)));

  pp->x = pp->corner[0][XI] + xc * pp->X[XI] + yc * pp->Y[XI];
  pp->y = pp->corner[0][YI] + xc * pp->X[YI] + yc * pp->Y[YI];
  pp->z = pp->corner[0][ZI] + xc * pp->X[ZI] + yc * pp->Y[ZI];

}

double normalize(vector)
  double vector[3];
{
  double length;
  int i;

  length = sqrt( vector[0]*vector[0] 
		+ vector[1]*vector[1] 
		+ vector[2]*vector[2]);
    
  for (i=0; i<3; i++) vector[i] = vector[i] / length;

  return length;
}

/* Assumes the vectors are normalized. */
int If_Equal(vector1, vector2)
  double vector1[3], vector2[3];
{
  int i;

  for (i=0; i<3; i++)
    if (fabs(vector1[i] - vector2[i]) > EQUIV_TOL) return FALSE;   
  return TRUE;
}

/* Calculates result_vector = vector1 X vector2. */
Cross_Product(vector1, vector2, result_vector)
  double vector1[], vector2[], result_vector[];
{