zbufinout.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"
#include "zbufGlobal.h"

double black, white;		/* densities corresponding to shades */

/*
  loads axes' lines
*/
void setupLine(axi, index, x1, y1, z1, x2, y2, z2)
double x1, y1, z1, x2, y2, z2;
double ***axi;
int index;
{
  axi[index][0][0] = x1;
  axi[index][0][1] = y1;
  axi[index][0][2] = z1;
  axi[index][1][0] = x2;
  axi[index][1][1] = y2;
  axi[index][1][2] = z2;
}

#if 1 == 0
/*
  replaces Song's make_charges_all_patches when only patches are
  to be dispalyed
  - requires minor mods to main in patran.c
*/
charge *make_charges_from_patches()
{
  int i, j, k;
  PATCH *patch_pntr;
  extern PATCH *start_patch;
  charge *start_chg, *charge_pntr;
  GRID *current_grid;
  extern GRID *start_grid;

  start_chg = NULL;

  /* for each patch, */
  patch_pntr = start_patch;
  while(patch_pntr != NULL) {

    /* create a charge struct */
    if(start_chg == NULL) {
      CALLOC(charge_pntr, 1, charge *, ON, AMSC);
      start_chg = charge_pntr;
    }
    else {
      CALLOC(charge_pntr->next, 1, charge *, ON, AMSC);
      charge_pntr = charge_pntr->next;
    }

    /* for each patch grid, find and load the coordinates 
       (assumes 4 sides, just like Song does throughout) */
    for(i = 0; i < 4; i++) {
      current_grid = start_grid;
      while(current_grid != NULL) {
	if(patch_pntr->corner[i] == current_grid->ID) {
	  if(i == 0) {
	    CALLOC(charge_pntr->corner0, 3, double *, ON, AMSC)
	    for(k = 0; k < 3; k++) 
		charge_pntr->corner0[k] = current_grid->coord[k];
	  }
	  else if(i == 1) {
	    CALLOC(charge_pntr->corner1, 3, double *, ON, AMSC)
	    for(k = 0; k < 3; k++) 
		charge_pntr->corner1[k] = current_grid->coord[k];
	  }
	  else if(i == 2) {
	    CALLOC(charge_pntr->corner2, 3, double *, ON, AMSC)
	    for(k = 0; k < 3; k++) 
		charge_pntr->corner2[k] = current_grid->coord[k];
	  }
	  else if(i == 3) {
	    CALLOC(charge_pntr->corner3, 3, double *, ON, AMSC)
	    for(k = 0; k < 3; k++) 
		charge_pntr->corner3[k] = current_grid->coord[k];
	  }
	  else {
	    fprintf(stderr, "make_charges_from_patches: panel has more than 4 sides\n");
	    exit(1);
	  }
	  break;
	}
	current_grid = current_grid->next;
      }
    }

    /* put in the shape - always rectangular */
    charge_pntr->shape = 4;
    patch_pntr = patch_pntr->next;
  }
  return(start_chg);
}
#endif

/*
  set the grey levels taking into account the rd_, rc_, and rb_ options
  - levels are scaled to match the range of densities present
  - sets the values of the extremal densities (`black' and `white')
*/
void figure_grey_levels(face_list, chgs, chglist, use_density)
double *chgs;
charge *chglist;
face **face_list;
int use_density;
{
  int first, i;
  double dif;
  charge *panel;
  extern int rc_, rb_, rd_;
  extern double black, white;
  extern ITER *kq_num_list;

  /* find the minimum and maximum charge density values */
  first = TRUE;
  for(panel = chglist; panel != NULL; panel = panel->next) {
    if(panel->dummy) continue;

    /* skip if panel is on conductor in q picture kill list */
    if(panel->surf->type == CONDTR || panel->surf->type == BOTH) {
      if(want_this_iter(kq_num_list, panel->cond)) continue;
    }

    /* skip if removing DIELEC's and panel is on dielectric i/f
       - type both interfaces are considered conductors for removal purposes */
    if(panel->surf->type == DIELEC && rd_ == TRUE) continue;

    if(first == TRUE) {
      first = FALSE;
      if(use_density) black = white = chgs[panel->index]/panel->area;
      else black = white = chgs[panel->index];
    }
    else {
      if(chgs[panel->index] == 0.0) {
	dif = 0.0;
      }
      if(use_density) {
	black = MAX(black, chgs[panel->index]/panel->area);
	white = MIN(white, chgs[panel->index]/panel->area);
      }
      else {
	black = MAX(black, chgs[panel->index]);
	white = MIN(white, chgs[panel->index]);
      }	  
    }

  }

  /* assign the grey levels - 0.0 = white, 1.0 = black */
  dif = black - white;
  for(panel = chglist, i = 0; panel != NULL; panel = panel->next) {
    if(panel->dummy) continue;

    /* skip if panel is on conductor in q picture kill list */
    if(panel->surf->type == CONDTR || panel->surf->type == BOTH) {
      if(want_this_iter(kq_num_list, panel->cond)) continue;
    }

    /* skip if removing DIELEC's and panel is on dielectric i/f
       - type both interfaces are considered conductors for removal purposes */
    if(panel->surf->type == DIELEC && rd_ == TRUE) continue;

    if(use_density)
	face_list[i]->greylev = (chgs[panel->index]/panel->area - white)/dif;
    else face_list[i]->greylev = (chgs[panel->index] - white)/dif;
    i++;
  }
}

/*
  attempt to read charge density information to set panel grey levels
  file format:
     (....lines of garbage....)
     Panel charges, iteration <iter>
     <chg->index> q/A = <chargedensity> (...garbage...)
     <chg->index> q/A = <chargedensity> (...garbage...)
     (...)
     <chg->index> q/A = <chargedensity> (...garbage...)
     End panel charges
     (....lines of garbage....)
  CHARGE DENSITIES ARE NOW READ DIRECTLY FROM VECTORS IN SYSTEM STRUCTS
  THIS FUNCTION IS NOT USED
*/
void get_charge_densities(q, file, iter)
double *q;
char *file;
int iter;
{
  int index, linecnt, header_found, type;
  char str1[BUFSIZ], str2[BUFSIZ], str3[BUFSIZ], linein[BUFSIZ];
  double density;
  FILE *fp, *fopen();

  if((fp = fopen(file, "r")) == NULL) {
    fprintf(stderr, 
	    "get_charge_densities: can't open charge file\n  `%s'\nto read\n", 
	    file);
    exit(1);
  }

  linecnt = 0;
  header_found = FALSE;
  while(fgets(linein, sizeof(linein), fp) != NULL) {
    linecnt++;

    if(!header_found) {
      /* look for "Panel charges, iteration <iter>" */
      if(sscanf(linein, "%s %s %s %d", str1, str2, str3, &index) == 4) {
	if(!strcmp(str1, "Panel") && !strcmp(str2, "charges,")
	      && !strcmp(str3, "iteration") && index == iter) {
	  header_found = TRUE;
	}
      }
      continue;
    }

    /* check for q/A line */
    if(sscanf(linein, "%d %s %s %lf", &index, str1, str2, &density) == 4) {
      if(!strcmp(str1, "q/A") && !strcmp(str2, "=")) {
	q[index] = density;
      }
    }
    /* check for end of list line */
    else if(sscanf(linein, "%s %s %s", str1, str2, str3) == 3) {
      if(!strcmp(str1, "End") && !strcmp(str2, "panel") 
	 && !strcmp(str3, "charges")) {
	break;
      }
    }
    else {
      fprintf(stderr, 
	      "get_charge_densities: bad charge file format, line %d:\n%s\n",
	      linecnt, linein);
      exit(1);
    }
  }

  if(!header_found) {
    fprintf(stderr, 
	    "get_charge_densities: can't find iteration %d data in\n `%s'\n",
	    iter, file);
    exit(1);
  }

}

/*
  figures the corner coordinates in absolute coordinates
*/
void getAbsCoord(vec, panel, num)
double *vec;
charge *panel;
int num;
{
  double *cor = panel->corner[num];
  double *x = panel->X, *y = panel->Y, *z = panel->Z;

  vec[0] = panel->x + cor[0]*x[0] + cor[1]*y[0] + cor[2]*z[0];
  vec[1] = panel->y + cor[0]*x[1] + cor[1]*y[1] + cor[2]*z[1];
  vec[2] = panel->z + cor[0]*x[2] + cor[1]*y[2] + cor[2]*z[2];

}

  
/*
  transfer fastcap panel info to face structs
*/
face **fastcap2faces(numfaces, chglist, q, use_density)
int *numfaces, use_density;	/* use_density = TRUE => use q/A not q */
charge *chglist;
double *q;
{
  int i, j, dummy;
  int autmom, autlev, numMom, numLev;
  char infile[BUFSIZ];
  double dot(), getPlane(), relperm;
  charge *chgp;
  face *head, *tail, **faces;
  extern int x_, k_, q_iter, q_, rc_, rd_, rb_;
  extern double ***axes;
  extern double axeslen;
  extern char *q_file;
  extern double black, white, linewd;
  surface *surf_list, *input_surfaces();
  ssystem *sys;
  int qindex=1, cindex=1;
  double tavg[3], lavg[3];
  extern ITER *kq_num_list;

  /* transfer info to face structs (a waste but saves wrtting new fnt end) */
  for(chgp = chglist, head = NULL, *numfaces = 0; chgp != NULL; 
      chgp = chgp->next) {
    if(chgp->dummy) continue;

#if RMWEDGE == ON
    /* remove all dielectric panels in the first quadrant */
    tavg[0] = tavg[1] = tavg[2] = 0.0;
    for(i = 0; i < chgp->shape; i++) {
      getAbsCoord(lavg, chgp, i);
      for(j = 0; j < 3; j++) tavg[j] += lavg[j];
    }
    if(chgp->surf->type == DIELEC && tavg[0] > 1.3 && tavg[1] > 1.3 
       && tavg[2] > 1.3)
	continue;
#endif

    /* skip if panel is on conductor in q picture kill list */
    if(chgp->surf->type == CONDTR || chgp->surf->type == BOTH) {
      if(want_this_iter(kq_num_list, chgp->cond)) continue;
    }

    /* skip if removing DIELEC's and panel is on dielectric i/f
       - type both interfaces are considered conductors for removal purposes */
    if(chgp->surf->type == DIELEC && rd_ == TRUE) continue;

    /* create and link in a new face */
    if(head == NULL) {
      CALLOC(head, 1, face, ON, AMSC);
      tail = head;
    }
    else {
      CALLOC(tail->next, 1, face, ON, AMSC);
      tail = tail->next;