/*
// d2vrml.c
// created by F.Mizutani  since 961018
//  961018 marge d2p and p2vrml
//  961113 cleaning
//  960611 BUG fix: normalize vector
*/

/* unix header file */
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>

#define  DOSCOD  0x1111
#define  UNXCOD  0x2222
#ifdef   DOS
#define  DCODE   DOSCOD
#else
#define  DCODE   UNXCOD
#endif
#define  PCODE   11002
#define  GCODE   10004    /* program code */
#define  PMAX    32       /* parameters   */
#define  DIV     50       /* find neighbor polygon  < step/DIV */
#define	 OBTMAX  64
#define  FALSE   0
#define  TRUE    !FALSE

float    viso; /* iso vale */
float    vvvv;
float    xmax, ymax, zmax, datamax;
float    xmin, ymin, zmin, datamin;

typedef struct {
	float x;
	float y;
	float z;
	float data;
        int   flag;
}IDST;

IDST *ids;

typedef struct {
	float x;
	float y;
	float z;
} XYZ;

typedef struct {
        float x[3];
        float y[3];
        float z[3];
        float vx[3];
        float vy[3];
        float vz[3];
        float vnx, vny, vnz;
} PGD;

typedef struct {
	long    mn;
	float	mw;
	float	x;
	float	y;
	float	z;
} MOL;

PGD     *pgd,  *pgd2;
int     pgdct = 0, pgdct2 = 0;
int     pgdmx;
int     pgdmode;
int     pgs;

long    idsz,   /* Mesh data size. */
        idmx,   /* Mesh size for x direction. */
        idmy,   /* Mesh size for y direction. */
        idmz;   /* Mesh size for z direction. */
float   mbgx, mbgy, mbgz;
float   mstx, msty, mstz;
float   mdvx, mdvy, mdvz;
long    dtype;
long	oct;
long	obt[OBTMAX];
MOL	*mol;
long	molct;
int     sfg;
int     dble = 1;
int     cmode;
FILE    *out;
long    pl[PMAX];
float   pf[PMAX];
char    fileout[1024];

#define         PGMAX           12       /* polygon count */
#define         DLMAX           30       /* display list */
#define         ATOMMAX         100      /* atomic number */
#define         RATOM           0.4      /* Atom radius */
#define         HATOM           0.6      /* H radius      RATOM * HATOM */
#define         WATOM           0.25     /* Atomic wire   RATOM * WATOM */
#define         ADIS            3.0      /* Atomic distance */
#define         ADISH           2.9      /* Atomic distance */
                                         /* 961018 H2O modify 3.0 -> 2.9 */
#define         LATOM           0.5
#define         MESHWIDTH       0.05
#define         AXISWIDTH       0.01
#define         AXISLENGTH      1.05
#define         ARROWWIDTH      0.03
#define         ARROWLENGTH     0.08
#define         DIS             5.0
#define         DIVM            0.1
#define         DIVR            1.0
#define         INIX           -2.0
#define         INIY            5.0
#define         INIZ            1.2
#define         EYEX            0.0
#define         EYEY            0.0
#define         EYEZ            DIS
#define         PI              3.14159
#define		PI2		1.57089

/* added by F.Mizutani */
int           Drag = 0;
int           Scale = 0;
int           xpos = 0;
int           ypos = 0;
int           xsize = 500;
int           ysize = 500;
float         spin_x = 0.0;
float         spin_y = 0.0;
float         scalev = 0.5;

double        eyex = EYEX;
double        eyey = EYEY;
double        eyez = EYEZ;
double        cenx = 0.0;
double        ceny = 0.0;
double        cenz = 0.0;
double        dist = DIS;
double        divm = DIVM;

int             SOLID[PGMAX];
int             WIRE    = 0;
int             COLORID = 1;  /* 0:RED 1:GREEN 2:BLUE 3:YELLOW */
int             SHADE   = 1;
int             AXIS    = 1;
int             ALIAS   = 1;
int             FOG     = 1;
int             MESH    = 0;
int             LMODE   = 0;
int             BOND    = 1;
int             GLASS   = 0;
int             BACK    = 0;
int             ROTATE  = 0;
int             CLIP    = 0;
float           wkx,  wky,  wkz,  wkmax;

int     xbig, ybig, zbig;
int     pgv[PGMAX];

char     *m[105] = { "",
           "H",  "HE", "LI", "BE", "B",  "C",  "N",  "O",  "F",  "NE",
           "NA", "MG", "AL", "SI", "P",  "S",  "CL", "AR", "K",  "CA",
           "SC", "TI", "V",  "CR", "MN", "FE", "CO", "NI", "CU", "ZN",
           "GA", "GE", "AS", "SE", "BR", "KR", "RB", "SR", "Y",  "ZR",
           "NB", "MO", "TC", "RU", "RH", "PD", "AG", "CD", "IN", "SN",
           "SB", "TE", "I",  "XE", "CS", "BA", "LA", "CE", "PR", "ND",
           "PM", "SM", "EU", "GD", "TB", "DY", "HO", "ER", "TM", "YB",
           "LU", "HF", "TA", "W",  "RE", "OS", "IR", "PT", "AU", "HG",
           "TL", "PB", "BI", "PO", "AT", "RN", "FR", "RA", "AC", "TH",
           "PA", "U",  "NP", "PU", "AM", "CM", "BK", "CF", "ES", "FM",
           "MD", "NO", "LR", ""
         };
float    mc[22][10] = {
           {0.4,0.2,0.1, 0.7,0.4,0.2, 1.0,1.0,1.0, 0.95},
           {0.5,0.5,0.5, 0.8,0.8,0.8, 1.0,1.0,1.0, 0.95}, /* H  */
           {0.4,0.2,0.1, 0.7,0.4,0.2, 1.0,1.0,1.0, 0.95}, /* He */
           {0.4,0.2,0.1, 0.7,0.4,0.2, 1.0,1.0,1.0, 0.95}, /* Li */
           {0.4,0.2,0.1, 0.7,0.4,0.2, 1.0,1.0,1.0, 0.95}, /* Be */
           {0.4,0.2,0.1, 0.7,0.4,0.2, 1.0,1.0,1.0, 0.95}, /* B  */
           {0.3,0.3,0.3, 0.6,0.6,0.6, 1.0,1.0,1.0, 0.95}, /* C  */
           {0.1,0.1,0.4, 0.2,0.2,0.7, 1.0,1.0,1.0, 0.95}, /* N  */
           {0.4,0.1,0.1, 0.8,0.1,0.1, 1.0,1.0,1.0, 0.95}, /* O  */
           {0.4,0.2,0.1, 0.7,0.4,0.2, 1.0,1.0,1.0, 0.95}, /* F  */
           {0.4,0.2,0.1, 0.7,0.4,0.2, 1.0,1.0,1.0, 0.95}, /* Ne */
           {0.1,0.4,0.1, 0.2,0.7,0.2, 1.0,1.0,1.0, 0.95}, /* Na */
           {0.1,0.4,0.1, 0.2,0.7,0.2, 1.0,1.0,1.0, 0.95}, /* Mg */
           {0.1,0.4,0.1, 0.2,0.7,0.2, 1.0,1.0,1.0, 0.95}, /* Al */
           {0.4,0.2,0.1, 0.7,0.4,0.2, 1.0,1.0,1.0, 0.95}, /* Si */
           {0.5,0.5,0.1, 0.2,0.7,0.2, 1.0,1.0,1.0, 0.95}, /* P  */
           {0.1,0.4,0.1, 0.2,0.7,0.2, 1.0,1.0,1.0, 0.95}, /* S  */
           {0.4,0.4,0.1, 0.7,0.7,0.2, 1.0,1.0,1.0, 0.95}, /* Cl */
           {0.4,0.2,0.1, 0.7,0.4,0.2, 1.0,1.0,1.0, 0.95}, /* Ar */
           {0.1,0.4,0.1, 0.2,0.7,0.2, 1.0,1.0,1.0, 0.95}, /* K  */
           {0.1,0.4,0.1, 0.2,0.7,0.2, 1.0,1.0,1.0, 0.95}, /* Ca */
           {0.4,0.2,0.1, 0.7,0.4,0.2, 1.0,1.0,1.0, 0.95}  /* Sc */
         };
float    mf[55] = {
           1.70,
           1.20/*H */, 1.50/*He*/, 1.70/*Li*/, 1.70/*Be*/, 1.70/*B */,
           1.70/*C */, 1.50/*N */, 1.40/*O */, 1.35/*F */, 1.59/*Ne*/,
           1.70/*Na*/, 1.70/*Mg*/, 1.70/*Al*/, 1.70/*Si*/, 1.90/*P */,
           1.85/*S */, 1.80/*Cl*/, 1.59/*Ar*/, 1.70/*K */, 1.70/*Ca*/,
           1.70/*Sc*/, 1.70/*TI*/, 1.70/*V */, 1.70/*CR*/, 1.70/*MN*/,
           1.70/*FE*/, 1.70/*CO*/, 1.70/*NI*/, 1.70/*CU*/, 1.70/*ZN*/,
           1.70/*GA*/, 1.70/*GE*/, 2.00/*AS*/, 2.00/*SE*/, 1.95/*BR*/,
           2.01/*KR*/, 1.70/*RB*/, 1.70/*SR*/, 1.70/*Y */, 1.70/*ZR*/,
           1.70/*NB*/, 1.70/*MO*/, 1.70/*TC*/, 1.70/*RU*/, 1.70/*RH*/,
           1.70/*PD*/, 1.70/*AG*/, 1.70/*CD*/, 1.70/*IN*/, 1.70/*SN*/,
           2.20/*SB*/, 2.20/*TE*/, 2.15/*I */, 2.20/*XE*/
         };

void swapper(p, n)
char    *p;
int     n;
{
  char d[4];
  int  i, j;

  for (j = 0; j < n; j++){
    for (i = 0; i < 4; i++) d[3 - i] = p[i + 4 * j];
    for (i = 0; i < 4; i++) p[i + 4 * j] = d[i];
  }
}

void m_alloc(size)
int	size;
{
  if ((ids = (IDST *)malloc(size * sizeof(IDST)))==NULL){
    perror("ids malloc error "); exit(1);
  }
  pgdmx = size / 3;
  if ((pgd  = (PGD *)malloc(pgdmx * sizeof(PGD))) == NULL){
    perror("pgd malloc error ");  exit(1);
  }
  if ((pgd2 = (PGD *)malloc(pgdmx * sizeof(PGD))) == NULL){
    perror("pgd malloc error ");  exit(1);
  }
}

/* 
// read mesh data from input file
// idsz :  Mesh data size
// idmx :  Mesh size for x direction
// idmy :  Mesh size for y direction
// idmz :  Mesh size for z direction
// ids  :  Mesh struct data
*/
void InputData(fp)
FILE  *fp;
{
        long           i, j, k, n, dcode, pcode;
        float          x, y, z, data, *dt;
        static float   smax =   (float)(1.0e30);
        static float   smin = - (float)(1.0e30);

        printf("ReadFile ... ");
        fread(&dcode, sizeof(long), 1, fp);
	if (dcode != DOSCOD && dcode != UNXCOD){
	    fclose(fp);
	    printf("This data is not made by isoconv.\n");
	    fclose(fp);
	}
	printf("Done.\n");
	if (dcode == DCODE) sfg = FALSE;
	else                sfg = TRUE;
        fread(&pcode, sizeof(long), 1, fp);
	if (sfg) swapper((char *)&pcode, 1);
	if (pcode != PCODE){
	    fclose(fp);
	    printf("This data is not made by isoconv.\n");
	    exit(1);
	}
        fread(&dtype, sizeof(long), 1, fp);
	fread(&datamax, sizeof(float), 1, fp);
	fread(&datamin, sizeof(float), 1, fp);
	if (sfg) swapper((char *)&dtype, 1);
	if (sfg) swapper((char *)&datamax, 1);
	if (sfg) swapper((char *)&datamin, 1);
        fread(&idmx,  sizeof(long), 1, fp);
        fread(&idmy,  sizeof(long), 1, fp);
        fread(&idmz,  sizeof(long), 1, fp);
	if (sfg) swapper((char *)&idmx, 1);
	if (sfg) swapper((char *)&idmy, 1);
	if (sfg) swapper((char *)&idmz, 1);
        idsz = idmx * idmy * idmz;
        printf("Resolution  X=%d Y=%d Z=%d TOTAL=%d\n",idmx,idmy,idmz,idsz);
	m_alloc(idsz);
        fread(&mbgx, sizeof(float), 1, fp);
        fread(&mbgy, sizeof(float), 1, fp);
        fread(&mbgz, sizeof(float), 1, fp);
        fread(&mstx, sizeof(float), 1, fp);
        fread(&msty, sizeof(float), 1, fp);
        fread(&mstz, sizeof(float), 1, fp);
	if (sfg) swapper((char *)&mbgx, 1);
	if (sfg) swapper((char *)&mbgy, 1);
	if (sfg) swapper((char *)&mbgz, 1);
	if (sfg) swapper((char *)&mstx, 1);
	if (sfg) swapper((char *)&msty, 1);
	if (sfg) swapper((char *)&mstz, 1);
        printf("Begin AXIS  X=%+f Y=%+f Z=%+f\n",mbgx,mbgy,mbgz);
        printf("Step  AXIS  X=%+f Y=%+f Z=%+f\n",mstx,msty,mstz);
        mdvx = mstx / DIV; mdvy = msty / DIV; mdvz = mstz / DIV;

	fread(&molct, sizeof(long), 1, fp);
	if (sfg) swapper((char *)&molct, 1);
	mol = (MOL *)malloc(molct * sizeof(MOL));
	fread(mol, sizeof(MOL), molct, fp);
	if (sfg) swapper((char *)mol, sizeof(MOL)/sizeof(long));
	for (i = 0; i < molct; i++){
	    printf("%3d: No.%3d  W=%+f  X=%+f Y=%+f Z=%+f\n", i+1, 
		mol[i].mn, mol[i].mw, mol[i].x, mol[i].y, mol[i].z);
	}

	if (dtype == 1){
	    printf("Orbital data\n");
	    fread(&oct, sizeof(long), 1,   fp);
	    if (sfg) swapper((char *)&oct, 1);
	    fread(obt,  sizeof(long), oct, fp);
	    if (sfg) swapper((char *)obt, oct);
	}
	else printf("Density data\n");

        dt = (float *)malloc(idmz * sizeof(float));

        xmax    = smin;
        ymax    = smin;
        zmax    = smin;
        /*datamax = smin;*/
        xmin    = smax;
        ymin    = smax;
        zmin    = smax;
        /*datamin = smax;*/

	printf("X loop ");
        for(i = 0; i < idmx; i++){
            printf("[%d]", i + 1);
	    fflush(stdout);
            for (j = 0; j < idmy; j++){
                if (fread(dt, sizeof(float), idmz, fp) != idmz){
                    printf("data read error ");
                    exit(1);
		}
	        if (sfg) swapper((char *)dt, idmz);
                for (k = 0; k < idmz; k++){
                    /*data = dt[k];*/
                    x = mbgx + mstx * i;
                    y = mbgy + msty * j;
                    z = mbgz + mstz * k;
                    n = i + j * idmx + k * idmx * idmy;
                    ids[n].x    = x;
                    ids[n].y    = y;
                    ids[n].z    = z;
                    ids[n].data = dt[k];
/*
   printf("%+f %+f %+f %+f\n",ids[i].x,ids[i].y,ids[i].z,ids[i].data);
*/
                    if (x > xmax) xmax = x; 
                    if (x < xmin) xmin = x; 
                    if (y > ymax) ymax = y; 
                    if (y < ymin) ymin = y; 
                    if (z > zmax) zmax = z; 
                    if (z < zmin) zmin = z; 
                    /*if (data > datamax) datamax = data;*/
                    /*if (data < datamin) datamin = data;*/
                }
            }
	}
        printf("\n");
        printf("x=%+e:%+e\n",xmin,xmax);
        printf("y=%+e:%+e\n",ymin,ymax);
        printf("z=%+e:%+e\n",zmin,zmax);
        printf("data=%+e:%+e\n",datamin,datamax);
}

void FlagSelect()
{
  int     i, j, k, m, n;

  if (viso > 0){
    for (k = 0; k < idmz; k++){
      m = k * idmx * idmy;
      for (j = 0; j < idmy; j++){
        n = j * idmx + m;
        for(i = 0; i < idmx; i++){
          if (ids[n].data > viso) ids[n].flag = 1;
          else                    ids[n].flag = 0;
          n++;
        }
      }
    }
  } 
  else {
    for (k = 0; k < idmz; k++){
      m = k * idmx * idmy; 
      for (j = 0; j < idmy; j++){
        n = j * idmx + m;
        for(i = 0; i < idmx; i++){
          if (ids[n].data < viso) ids[n].flag = 1;
          else                    ids[n].flag = 0;
          n++;
        }
      }
    }
  }
}


void NormalizeVector()
{
    int     i, j, k, l;
    int     ct;
    int     nb[100], xyz[100];
    float   length, x, y, z;
    PGD     *ppp;
    int     pppct;

    return;

    if (pgdmode == 0){ ppp = pgd;  pppct = pgdct;  }
    if (pgdmode == 1){ ppp = pgd2; pppct = pgdct2; }
    for (i = 0; i < pppct; i++){
        if (i % 100 == 0){ printf("[%d]", i); fflush(stdout); }
        if (i % 1000 == 0){ printf("\n"); fflush(stdout); }
        for (j = 0; j < 3; j++){
            if (ppp[i].vx[j] != 0.0 || ppp[i].vy[j] != 0.0 || ppp[i].vz[j] == 0.0) continue;
            x = ppp[i].vnx; y = ppp[i].vny; z = ppp[i].vnz;
            ct = 1;
            for (k = i + 1; k < pppct; k++){
                if (ppp[i].z[0] + mstz * 2 < ppp[k].z[0]) break;
                for (l = 0; l < 3; l++){
                    if (fabs(ppp[i].x[j] - ppp[k].x[l]) < mdvx &&
                        fabs(ppp[i].y[j] - ppp[k].y[l]) < mdvy &&
                        fabs(ppp[i].z[j] - ppp[k].z[l]) < mdvz &&
                        ct < 100){
                        nb[ct] = k; xyz[ct] = l; ct++;
                        x += ppp[k].vnx; y += ppp[k].vny; z += ppp[k].vnz;
			break;
                    }
                }
            }
	    length = sqrt(x*x + y*y + z*z);
            ppp[i].vx[j] = (float)(x / length);
            ppp[i].vy[j] = (float)(y / length);
            ppp[i].vz[j] = (float)(z / length);
            for (k = 1; k < ct; k++){
                ppp[nb[k]].vx[xyz[k]] = ppp[i].vx[j];
                ppp[nb[k]].vy[xyz[k]] = ppp[i].vy[j];
                ppp[nb[k]].vz[xyz[k]] = ppp[i].vz[j]; 
            }
        }
    }
    printf("\n");
    if (pgdmode == 0) pgdct  = pppct;
    if (pgdmode == 1) pgdct2 = pppct;
}

XYZ CalcNormal(v)
XYZ v[3];
{
        int     i;
        XYZ     e[2], n;
        float   length;
        float   x, y, z;

        for (i = 0; i < 2; i++){
                e[i].x = v[i + 1].x - v[0].x;
                e[i].y = v[i + 1].y - v[0].y;
                e[i].z = v[i + 1].z - v[0].z;
	}
        x = e[0].y * e[1].z - e[0].z * e[1].y;
        y = e[0].z * e[1].x - e[0].x * e[1].z;
        z = e[0].x * e[1].y - e[0].y * e[1].x;
        length = (float)sqrt(x*x + y*y + z*z);
        n.x = x / length;
        n.y = y / length;
        n.z = z / length;
	return n;
}

void SetTriangle(v)
XYZ   v[3];
{
    XYZ   n;
    int   i;
    float ck;
    PGD   *ppp;
    int   pppct;

    if (pgdmode == 0){ ppp = pgd;  pppct = pgdct;  }
    if (pgdmode == 1){ ppp = pgd2; pppct = pgdct2; }
    for (i = 0; i < 3; i++){
      ppp[pppct].x[i] = v[i].x;
      ppp[pppct].y[i] = v[i].y;
      ppp[pppct].z[i] = v[i].z;
    }
    pppct++;
    if (pgdmode == 0) pgdct  = pppct;
    if (pgdmode == 1) pgdct2 = pppct;
    return;

    n = CalcNormal(v);
    ck = n.x * n.x + n.y * n.y + n.z * n.z;
    if (0.95 < ck && ck < 1.05){
        if (pppct < pgdmx){
            ppp[pppct].vnx = n.x;
            ppp[pppct].vny = n.y;
            ppp[pppct].vnz = n.z;
            for (i = 0; i < 3; i++){
                ppp[pppct].x[i] = v[i].x;
                ppp[pppct].y[i] = v[i].y;
                ppp[pppct].z[i] = v[i].z;
                ppp[pppct].vx[i] = 0.0;
                ppp[pppct].vy[i] = 0.0;
                ppp[pppct].vz[i] = 0.0;
            }
            pppct++;
        }
    }
/*    else printf("E:%+f %+f %+f\n",n.x,n.y,n.z); */
    else printf("!");
    if (pgdmode == 0) pgdct  = pppct;
    if (pgdmode == 1) pgdct2 = pppct;
}

XYZ CalcVertex(n1, n2)
int	n1;
int	n2;
{
        float   a;
        XYZ     v;

        a = (float)(viso - ids[n1].data)/(ids[n2].data - ids[n1].data);
	v.x = ids[n1].x + a*(ids[n2].x - ids[n1].x);
	v.y = ids[n1].y + a*(ids[n2].y - ids[n1].y);
	v.z = ids[n1].z + a*(ids[n2].z - ids[n1].z);
	return v;
}

void DisplayTriangle(wnode)
int	*wnode;
{
        XYZ   v[3];

        v[0] = CalcVertex(wnode[0],wnode[1]);
        v[2] = CalcVertex(wnode[2],wnode[3]);
        v[1] = CalcVertex(wnode[4],wnode[5]);
        SetTriangle(v);
}

void DisplayCorner(wnode)
int wnode[4];
{
        XYZ   v[3];

        v[0] = CalcVertex(wnode[0],wnode[1]);
        v[2] = CalcVertex(wnode[0],wnode[2]);
        v[1] = CalcVertex(wnode[0],wnode[3]);
        SetTriangle(v);
}

void Display2Tri(wnode)
int	*wnode;
{
        XYZ   v[3];

        v[0] = CalcVertex(wnode[0],wnode[1]);
        v[2] = CalcVertex(wnode[0],wnode[2]);
        v[1] = CalcVertex(wnode[3],wnode[5]);
        SetTriangle(v);

        v[0] = CalcVertex(wnode[0],wnode[2]);
        v[2] = CalcVertex(wnode[3],wnode[4]);
        v[1] = CalcVertex(wnode[3],wnode[5]);
        SetTriangle(v);
}

void SetCornerOn(cnode, i, j)
int cnode[8][8];
int i,j;
{
        int     wnode[4];

        if (j < 4){
		wnode[0] = cnode[i][j];
		wnode[1] = cnode[i][(j+1)%4];
		wnode[2] = cnode[i][j+4];
		wnode[3] = cnode[i][(j+3)%4];
	}
        else {
		wnode[0] = cnode[i][j];
		wnode[1] = cnode[i][4+(j+1)%4];
		wnode[2] = cnode[i][4+(j+3)%4];
		wnode[3] = cnode[i][j-4];
	}
        DisplayCorner(wnode);
}

void SetCornerOff(cnode, i, j)
int cnode[8][8];
int i,j;
{
        int     wnode[4];

        if (j < 4){
		wnode[0] = cnode[i][j];
		wnode[1] = cnode[i][(j+1)%4];
		wnode[2] = cnode[i][(j+3)%4];
		wnode[3] = cnode[i][j+4];
	}
        else {
		wnode[0] = cnode[i][j];
		wnode[1] = cnode[i][4+(j+1)%4];
		wnode[2] = cnode[i][j-4];
		wnode[3] = cnode[i][4+(j+3)%4];
	}
        DisplayCorner(wnode);
}

/*
// １のパターン
*/
void Flag1(cnode)
int cnode[8][8];
{
        int     i;

        for (i = 0; i < 8; i++){
	    if(ids[cnode[i][0]].flag == 1){
                SetCornerOn(cnode,i,0);
		break;
	    }
	}
}

/*
// １のパターンの反転パターン
*/
void Flag7(cnode)
int cnode[8][8];
{
        int     i;

        for (i = 0; i < 8; i++){
	    if(ids[cnode[i][0]].flag == 0){
                SetCornerOff(cnode,i,0);
		break;
	    }
	}
}

/*
// 選ばれた２点が、格子の底面にある場合
*/
void F2_Patern2_1(cnode,cyc,p0,p1)
int cnode[8][8];
int cyc,p0,p1;
{
        int     wnode[6];

	wnode[0] = cnode[cyc][p0];
	wnode[1] = cnode[cyc][p0+4];
	wnode[2] = cnode[cyc][(p0+3)%4];
	wnode[3] = cnode[cyc][p1];
	wnode[4] = cnode[cyc][(p1+1)%4];
	wnode[5] = cnode[cyc][p1+4];
        Display2Tri(wnode);
}

/*
// 選ばれた２点が、格子の上面にある場合
*/
void F2_Patern2_2(cnode, cyc, p0, p1)
int cnode[8][8];
int cyc,p0,p1;
{
        int     wnode[6];

	wnode[0] = cnode[cyc][p0];
	wnode[1] = cnode[cyc][4+(p0+3)%4];
	wnode[2] = cnode[cyc][p0-4];
	wnode[3] = cnode[cyc][p1];
	wnode[4] = cnode[cyc][p1-4];
	wnode[5] = cnode[cyc][4+(p1+1)%4];
        Display2Tri(wnode);
}

/*
// 選ばれた２点が、格子の側面にある場合
*/
void F2_Patern2_3(cnode, cyc, p0, p1)
int cnode[8][8];
int cyc,p0,p1;
{
        int     wnode[6];

	wnode[0] = cnode[cyc][p0];
	wnode[1] = cnode[cyc][(p0+3)%4];
	wnode[2] = cnode[cyc][(p0+1)%4];
	wnode[3] = cnode[cyc][p1];
	wnode[4] = cnode[cyc][4+(p1+1)%4];
	wnode[5] = cnode[cyc][4+(p1+3)%4];
        Display2Tri(wnode);
}

/*
// ２、３、４のパターン
*/
void Flag2(cnode)
int cnode[8][8];
{
        int     i, j, cyc;
        int     pnode[2];

        for (i = 0; i < 8; i++){
	    if (ids[cnode[i][0]].flag == 1){
		cyc = i;
		pnode[0] = 0;
		for(j=1;j<8;j++){
	    	    if (ids[cnode[i][j]].flag == 1){
			pnode[1] = j;
			break;
		    }
		}
		break;
	    }
	}
/*
	printf("Flag2:cyc=%d,0-%d \n",cyc,pnode[1]);
*/
	switch(pnode[1]) {
	case 1:
                F2_Patern2_1(cnode,cyc,pnode[0],pnode[1]);
		break;
	case 3:
                F2_Patern2_1(cnode,cyc,pnode[1],pnode[0]);
		break;
	case 4:
                F2_Patern2_3(cnode,cyc,pnode[0],pnode[1]);
		break;
	default:
                SetCornerOn(cnode,cyc,pnode[0]);
                SetCornerOn(cnode,cyc,pnode[1]);
		printf("[WARNING PAT2]");
	}
}
/*
// 選ばれた２点が、格子の底面にある場合
*/
void F6_Patern2_1(cnode, cyc, p0, p1)
int cnode[8][8];
int cyc,p0,p1;
{
        int     wnode[6];

	wnode[0] = cnode[cyc][p0];
	wnode[2] = cnode[cyc][p0+4];
	wnode[1] = cnode[cyc][(p0+3)%4];
	wnode[3] = cnode[cyc][p1];
	wnode[5] = cnode[cyc][(p1+1)%4];
	wnode[4] = cnode[cyc][p1+4];
        Display2Tri(wnode);
}

/*
// 選ばれた２点が、格子の上面にある場合
*/
void F6_Patern2_2(cnode, cyc, p0, p1)
int cnode[8][8];
int cyc,p0,p1;
{
        int     wnode[6];

	wnode[0] = cnode[cyc][p0];
	wnode[2] = cnode[cyc][4+(p0+3)%4];
	wnode[1] = cnode[cyc][p0-4];
	wnode[3] = cnode[cyc][p1];
	wnode[5] = cnode[cyc][p1-4];
	wnode[4] = cnode[cyc][4+(p1+1)%4];
        Display2Tri(wnode);
}

/*
// 選ばれた２点が、格子の側面にある場合
*/
void F6_Patern2_3(cnode, cyc, p0, p1)
int cnode[8][8];
int cyc,p0,p1;
{
        int     wnode[6];

	wnode[0] = cnode[cyc][p0];
	wnode[2] = cnode[cyc][(p0+3)%4];
	wnode[1] = cnode[cyc][(p0+1)%4];
	wnode[3] = cnode[cyc][p1];
	wnode[5] = cnode[cyc][4+(p1+1)%4];
	wnode[4] = cnode[cyc][4+(p1+3)%4];
        Display2Tri(wnode);
}
/*
// ２、３、４のパターンの反転パターン
*/
void Flag6(cnode)
int cnode[8][8];
{
        int     i, j, cyc;
        int     pnode[2];

        for (i = 0; i < 8; i++){
	    if (ids[cnode[i][0]].flag == 0){
		cyc = i;
		pnode[0] = 0;
		for(j=1;j<8;j++){
		    if (ids[cnode[i][j]].flag == 0){
			pnode[1] = j;
			break;
		    }
		}
		break;
	    }
	}
/*
	printf("Flag6:cyc=%d,0-%d \n",cyc,pnode[1]);
*/
	switch(pnode[1]) {
	case 1:
                F6_Patern2_1(cnode,cyc,pnode[0],pnode[1]);
		break;
	case 3:
                F6_Patern2_1(cnode,cyc,pnode[1],pnode[0]);
		break;
	case 4:
                F6_Patern2_3(cnode,cyc,pnode[0],pnode[1]);
		break;
	default:
                SetCornerOff(cnode,cyc,pnode[0]);
                SetCornerOff(cnode,cyc,pnode[1]);
		printf("[WARNING PAT6]");
	}
}
/*
// 選ばれた6点で三角形を描く
*/
void F3_Patern5(cnode, cyc, p1, p2, p3, p4, p5, p6)
int cnode[8][8];
int cyc,p1,p2,p3,p4,p5,p6;
{
        int     wnode[6];

	wnode[0] = cnode[cyc][p1];
	wnode[1] = cnode[cyc][p2];
	wnode[2] = cnode[cyc][p3];
	wnode[3] = cnode[cyc][p4];
	wnode[4] = cnode[cyc][p5];
	wnode[5] = cnode[cyc][p6];
        DisplayTriangle(wnode);
}
/*
// Patern 7
*/
void F3_Patern7(cnode, cyc, pnode)
int cnode[8][8];
int pnode[3];
{
        SetCornerOn(cnode,cyc,pnode[0]);
        SetCornerOn(cnode,cyc,pnode[1]);
        SetCornerOn(cnode,cyc,pnode[2]);
}

/*
//  For Case flag=3
//  Patern 5,6,7
*/
void Flag3(cnode)
int cnode[8][8];
{
        int     i, j, k, cyc;
        int     pnode[3];

	k = 1;
        for (i = 0; i < 8; i++){
	    if (ids[cnode[i][0]].flag == 1){
		cyc = i;
		pnode[0] = 0;
                for(j = 1; j < 8; j++){
		    if (ids[cnode[i][j]].flag == 1){
			pnode[k++] = j;
		    }
		}
		break;
	    }
	}
/*
	printf("Flag3:cyc=%d,0-%d-%d \n",cyc,pnode[1],pnode[2]);
*/
	switch(pnode[1]) {
	case 1:
		switch(pnode[2]) {
		case 2:  /* Patern 5:0-1-2 */
                        F3_Patern5(cnode,cyc,0,4,2,6,1,5);
                        F3_Patern5(cnode,cyc,0,4,0,3,2,6);
                        F3_Patern5(cnode,cyc,2,6,0,3,2,3);
			break;
		case 3:  /* Patern 5:0-1-3 */
                        F3_Patern5(cnode,cyc,1,5,0,4,3,7);
                        F3_Patern5(cnode,cyc,1,5,3,7,1,2);
                        F3_Patern5(cnode,cyc,1,2,3,7,2,3);
			break;
		case 4:  /* Patern 5:0-1-4 */
                        F3_Patern5(cnode,cyc,4,7,0,3,1,2);
                        F3_Patern5(cnode,cyc,4,7,1,2,4,5);
                        F3_Patern5(cnode,cyc,4,5,1,2,1,5);
			break;
		case 5:  /* Patern 5:0-1-5 */
                        F3_Patern5(cnode,cyc,0,3,1,2,5,6);
                        F3_Patern5(cnode,cyc,0,3,5,6,0,4);
                        F3_Patern5(cnode,cyc,0,4,5,6,4,5);
			break;
		case 6:  /* Patern 6:0-1-6 */
		case 7:  /* Patern 6:0-1-7 */
                        F2_Patern2_1(cnode,cyc,pnode[0],pnode[1]);
                        SetCornerOn(cnode,cyc,pnode[2]);
			break;
		}
		break;
	case 2:
		switch(pnode[2]) {
		case 3:  /* Patern 5:0-2-3 */
                        F3_Patern5(cnode,cyc,0,4,3,7,2,6);
                        F3_Patern5(cnode,cyc,0,4,2,6,0,1);
                        F3_Patern5(cnode,cyc,0,1,2,6,1,2);
			break;
		case 4:  /* Patern 6:0-2-4 */
                        F2_Patern2_3(cnode,cyc,pnode[0],pnode[2]);
                        SetCornerOn(cnode,cyc,pnode[1]);
			break;
		case 5:  /* Patern 7:0-2-5 */
                        F3_Patern7(cnode,cyc,pnode);
			break;
		case 6:  /* Patern 6:0-2-6 */
                        F2_Patern2_3(cnode,cyc,pnode[1],pnode[2]);
                        SetCornerOn(cnode,cyc,pnode[0]);
			break;
		case 7:  /* Patern 7:0-2-7 */
                        F3_Patern7(cnode,cyc,pnode);
			break;
		}
		break;
	case 3:
		switch(pnode[2]) {
		case 4:  /* Patern 5:0-3-4 */
                        F3_Patern5(cnode,cyc,2,3,0,1,4,5);
                        F3_Patern5(cnode,cyc,2,3,4,5,3,7);
                        F3_Patern5(cnode,cyc,3,7,4,5,4,7);
			break;
		case 5:  /* Patern 6:0-3-5 */
		case 6:  /* Patern 6:0-3-6 */
                        F2_Patern2_1(cnode,cyc,pnode[1],pnode[0]);
                        SetCornerOn(cnode,cyc,pnode[2]);
			break;
		case 7:  /* Patern 5:0-3-7 */
                        F3_Patern5(cnode,cyc,6,7,2,3,0,1);
                        F3_Patern5(cnode,cyc,6,7,0,1,4,7);
                        F3_Patern5(cnode,cyc,4,7,0,1,0,4);
			break;
		}
		break;
	case 4:
		switch(pnode[2]) {
		case 5:  /* Patern 5:0-4-5 */
                        F3_Patern5(cnode,cyc,5,6,4,7,0,3);
                        F3_Patern5(cnode,cyc,5,6,0,3,1,5);
                        F3_Patern5(cnode,cyc,1,5,0,3,0,1);
			break;
		case 6:  /* Patern 6:0-4-6 */
                        F2_Patern2_3(cnode,cyc,pnode[0],pnode[1]);
                        SetCornerOn(cnode,cyc,pnode[2]);
			break;
		case 7:  /* Patern 5:0-4-7 */
                        F3_Patern5(cnode,cyc,0,1,4,5,6,7);
                        F3_Patern5(cnode,cyc,0,1,6,7,0,3);
                        F3_Patern5(cnode,cyc,0,3,6,7,3,7);
			break;
		}
		break;
	case 5:
		switch(pnode[2]) {
		case 6:  /* Patern 6:0-5-6 */
                        F2_Patern2_2(cnode,cyc,pnode[1],pnode[2]);
                        SetCornerOn(cnode,cyc,pnode[0]);
			break;
		case 7:  /* Patern 7:0-5-7 */
                        F3_Patern7(cnode,cyc,pnode);
			break;
		}
		break;
	case 6:
		switch(pnode[2]) {
		case 7:  /* Patern 6:0-6-7 */
                        F2_Patern2_2(cnode,cyc,pnode[1],pnode[2]);
                        SetCornerOn(cnode,cyc,pnode[0]);
			break;
		}
	} /* end of switch(pnode[1]) */
}

/*
// Patern5 for flag=5
*/
void F5_Patern5(cnode, cyc, p1, p2, p3, p4, p5, p6)
int cnode[8][8];
int cyc,p1,p2,p3,p4,p5,p6;
{
        int     wnode[6];

	wnode[0] = cnode[cyc][p1];
	wnode[1] = cnode[cyc][p2];
	wnode[4] = cnode[cyc][p3];
	wnode[5] = cnode[cyc][p4];
	wnode[2] = cnode[cyc][p5];
	wnode[3] = cnode[cyc][p6];
        DisplayTriangle(wnode);
}
/*
// Patern 7 for flag=5
*/
void F5_Patern7(cnode, cyc, pnode)
int cnode[8][8];
int pnode[3];
{
        SetCornerOff(cnode,cyc,pnode[0]);
        SetCornerOff(cnode,cyc,pnode[1]);
        SetCornerOff(cnode,cyc,pnode[2]);
}

/*
//  For Case flag=5
//  flag=3 の反転パターン
//  Patern 5,6,7
*/
void Flag5(cnode)
int cnode[8][8];
{
        int     i, j, k, cyc;
        int     pnode[3];

	k = 1;
        for (i = 0; i < 8; i++){
	    if (ids[cnode[i][0]].flag == 0){
		cyc = i;
		pnode[0] = 0;
                for(j = 1; j < 8; j++){
		    if (ids[cnode[i][j]].flag == 0){
			pnode[k++] = j;
		    }
		}
		break;
	    }
	}
/*
	printf("Flag5:cyc=%d,0-%d-%d \n",cyc,pnode[1],pnode[2]);
*/
	switch(pnode[1]) {
	case 1:
		switch(pnode[2]) {
		case 2:  /* Patern 5:0-1-2 */
                        F5_Patern5(cnode,cyc,0,4,2,6,1,5);
                        F5_Patern5(cnode,cyc,0,4,0,3,2,6);
                        F5_Patern5(cnode,cyc,2,6,0,3,2,3);
			break;
		case 3:  /* Patern 5:0-1-3 */
                        F5_Patern5(cnode,cyc,1,5,0,4,3,7);
                        F5_Patern5(cnode,cyc,1,5,3,7,1,2);
                        F5_Patern5(cnode,cyc,1,2,3,7,2,3);
			break;
		case 4:  /* Patern 5:0-1-4 */
                        F5_Patern5(cnode,cyc,4,7,0,3,1,2);
                        F5_Patern5(cnode,cyc,4,7,1,2,4,5);
                        F5_Patern5(cnode,cyc,4,5,1,2,1,5);
			break;
		case 5:  /* Patern 5:0-1-5 */
                        F5_Patern5(cnode,cyc,0,3,1,2,5,6);
                        F5_Patern5(cnode,cyc,0,3,5,6,0,4);
                        F5_Patern5(cnode,cyc,0,4,5,6,4,5);
			break;
		case 6:  /* Patern 6:0-1-6 */
		case 7:  /* Patern 6:0-1-7 */
                        F6_Patern2_1(cnode,cyc,pnode[0],pnode[1]);
                        SetCornerOff(cnode,cyc,pnode[2]);
			break;
		}
		break;
	case 2:
		switch(pnode[2]) {
		case 3:  /* Patern 5:0-2-3 */
                        F5_Patern5(cnode,cyc,0,4,3,7,2,6);
                        F5_Patern5(cnode,cyc,0,4,2,6,0,1);
                        F5_Patern5(cnode,cyc,0,1,2,6,1,2);
			break;
		case 4:  /* Patern 6:0-2-4 */
                        F6_Patern2_3(cnode,cyc,pnode[0],pnode[2]);
                        SetCornerOff(cnode,cyc,pnode[1]);
			break;
		case 5:  /* Patern 7:0-2-5 */
                        F5_Patern7(cnode,cyc,pnode);
			break;
		case 6:  /* Patern 6:0-2-6 */
                        F6_Patern2_3(cnode,cyc,pnode[1],pnode[2]);
                        SetCornerOff(cnode,cyc,pnode[0]);
			break;
		case 7:  /* Patern 7:0-2-7 */
                        F5_Patern7(cnode,cyc,pnode);
			break;
		}
		break;
	case 3:
		switch(pnode[2]) {
		case 4:  /* Patern 5:0-3-4 */
                        F5_Patern5(cnode,cyc,2,3,0,1,4,5);
                        F5_Patern5(cnode,cyc,2,3,4,5,3,7);
                        F5_Patern5(cnode,cyc,3,7,4,5,4,7);
			break;
		case 5:  /* Patern 6:0-3-5 */
		case 6:  /* Patern 6:0-3-6 */
                        F6_Patern2_1(cnode,cyc,pnode[1],pnode[0]);
                        SetCornerOff(cnode,cyc,pnode[2]);
			break;
		case 7:  /* Patern 5:0-3-7 */
                        F5_Patern5(cnode,cyc,6,7,2,3,0,1);
                        F5_Patern5(cnode,cyc,6,7,0,1,4,7);
                        F5_Patern5(cnode,cyc,4,7,0,1,0,4);
			break;
		}
		break;
	case 4:
		switch(pnode[2]) {
		case 5:  /* Patern 5:0-4-5 */
                        F5_Patern5(cnode,cyc,5,6,4,7,0,3);
                        F5_Patern5(cnode,cyc,5,6,0,3,1,5);
                        F5_Patern5(cnode,cyc,1,5,0,3,0,1);
			break;
		case 6:  /* Patern 6:0-4-6 */
                        F6_Patern2_3(cnode,cyc,pnode[0],pnode[1]);
                        SetCornerOff(cnode,cyc,pnode[2]);
			break;
		case 7:  /* Patern 5:0-4-7 */
                        F5_Patern5(cnode,cyc,0,1,4,5,6,7);
                        F5_Patern5(cnode,cyc,0,1,6,7,0,3);
                        F5_Patern5(cnode,cyc,0,3,6,7,3,7);
			break;
		}
		break;
	case 5:
		switch(pnode[2]) {
		case 6:  /* Patern 6:0-5-6 */
                        F6_Patern2_2(cnode,cyc,pnode[1],pnode[2]);
                        SetCornerOff(cnode,cyc,pnode[0]);
			break;
		case 7:  /* Patern 7:0-5-7 */
                        F5_Patern7(cnode,cyc,pnode);
			break;
		}
		break;
	case 6:
		switch(pnode[2]) {
		case 7:  /* Patern 6:0-6-7 */
                        F6_Patern2_2(cnode,cyc,pnode[1],pnode[2]);
                        SetCornerOff(cnode,cyc,pnode[0]);
			break;
		}
	} /* end of switch(pnode[1]) */
}

/*
// Patern 9
// F4:0-1-3-4 のcaseを基準
*/
void F4_Patern9(cnode, cyc, p1, p2, p3, p4, p5, p7)
int cnode[8][8];
int cyc,p1,p2,p3,p4,p5,p7;
{
        F3_Patern5(cnode,cyc,p1,p5,p4,p5,p4,p7);
        F3_Patern5(cnode,cyc,p1,p5,p4,p7,p1,p2);
        F3_Patern5(cnode,cyc,p1,p2,p4,p7,p3,p7);
        F3_Patern5(cnode,cyc,p1,p2,p3,p7,p2,p3);
}

/*
// Patern 11
// F4:0-1-3-7 のcaseを基準
*/
void F4_Patern11(cnode, cyc, p0, p1, p2, p3, p4, p5, p6, p7, p8)
int cnode[8][8];
int cyc,p0,p1,p2,p3,p4,p5,p6,p7,p8;
{
        F3_Patern5(cnode,cyc,p1,p5,p0,p4,p1,p2);
        F3_Patern5(cnode,cyc,p0,p4,p6,p7,p1,p2);
        F3_Patern5(cnode,cyc,p0,p4,p4,p7,p6,p7);
        F3_Patern5(cnode,cyc,p1,p2,p6,p7,p2,p3);
}

/*
// Patern 13
*/
void F4_Patern13(cnode, cyc, pnode)
int cnode[8][8];
int pnode[4];
{
        SetCornerOn(cnode,cyc,pnode[0]);
        SetCornerOn(cnode,cyc,pnode[1]);
        SetCornerOn(cnode,cyc,pnode[2]);
        SetCornerOn(cnode,cyc,pnode[3]);
}

/*
// Patern 14
// F4:0-2-3-4 のcaseを基準
*/
void F4_Patern14(cnode, cyc, p0, p1, p2, p3, p4, p5, p6, p7, p8)
int cnode[8][8];
int cyc,p0,p1,p2,p3,p4,p5,p6,p7,p8;
{
        F3_Patern5(cnode,cyc,p0,p1,p4,p5,p1,p2);
        F3_Patern5(cnode,cyc,p4,p5,p3,p7,p1,p2);
        F3_Patern5(cnode,cyc,p1,p2,p3,p7,p2,p6);
        F3_Patern5(cnode,cyc,p4,p5,p4,p7,p3,p7);
}


/*
//  For Case flag=4
//  Patern 8,9,10,11,12,13,14
*/
void Flag4(cnode)
int cnode[8][8];
{
        int     i, j, k, cyc;
        int     pnode[4];

	k = 1;
        for (i = 0; i < 8; i++){
		if (ids[cnode[i][0]].flag == 1){
			cyc = i;
			pnode[0] = 0;
                        for(j = 1; j < 8; j++){
				if (ids[cnode[i][j]].flag == 1){
					pnode[k++] = j;
				}
			}
			break;
		}
	}
/*
	printf("Flag4:cyc=%d,0-%d-%d-%d \n",cyc,pnode[1],pnode[2],pnode[3]);
*/
	switch(pnode[1]) {
	  case 1:
		switch(pnode[2]) {
		  case 2:
			switch(pnode[3]) {
			  case 3:  /* patern 8:0-1-2-3 */
                                F3_Patern5(cnode,cyc,0,4,3,7,1,5);
                                F3_Patern5(cnode,cyc,1,5,3,7,2,6);
				break;
			  case 4:  /* Patern 11:0-1-2-4 */
                                F4_Patern11(cnode,cyc,1,2,3,0,5,6,7,4);
				break;
			  case 5:  /* patern 9:0-1-2-5 */
                                F4_Patern9(cnode,cyc,2,3,0,5,6,4);
				break;
			  case 6:  /* Patern 14:0-1-2-6 */
                                F4_Patern14(cnode,cyc,2,3,0,1,6,7,4,5);
				break;
			  case 7:  /* Patern 12:0-1-2-7 */
                                F3_Patern5(cnode,cyc,0,4,2,6,1,5);
                                F3_Patern5(cnode,cyc,0,4,0,3,2,6);
                                F3_Patern5(cnode,cyc,2,6,0,3,2,3);
                                SetCornerOn(cnode,cyc,7);
				break;
			}
			break;
		  case 3:
			switch(pnode[3]) {
			  case 4:  /* patern 9:0-1-3-4 */
                                F4_Patern9(cnode,cyc,1,2,3,4,5,7);
				break;
			  case 5:  /* Patern 14:0-1-3-5 */
                                F4_Patern14(cnode,cyc,1,2,3,0,5,6,7,4);
				break;
			  case 6:  /* Patern 12:0-1-3-6 */
                                F3_Patern5(cnode,cyc,1,5,0,4,3,7);
                                F3_Patern5(cnode,cyc,1,5,3,7,1,2);
                                F3_Patern5(cnode,cyc,1,2,3,7,2,3);
                                SetCornerOn(cnode,cyc,6);
				break;
			  case 7:  /* Patern 11:0-1-3-7 */
                                F4_Patern11(cnode,cyc,0,1,2,3,4,5,6,7);
				break;
			}
			break;
		  case 4:
			switch(pnode[3]) {
			  case 5:  /* patern 8 0-1-4-5 */
                                F3_Patern5(cnode,cyc,4,7,0,3,5,6);
                                F3_Patern5(cnode,cyc,5,6,0,3,1,2);
				break;
			  case 6:  /* Patern 12:0-1-4-6 */
                                F3_Patern5(cnode,cyc,4,7,0,3,1,2);
                                F3_Patern5(cnode,cyc,4,7,1,2,4,5);
                                F3_Patern5(cnode,cyc,4,5,1,2,1,5);
                                SetCornerOn(cnode,cyc,6);
				break;
			  case 7:  /* Patern 14:0-1-4-7 */
                                F4_Patern14(cnode,cyc,4,5,1,0,7,6,2,3);
				break;
			}
			break;
		  case 5:
			switch(pnode[3]) {
			  case 6:  /* Patern 11:0-1-5-6 */
                                F4_Patern11(cnode,cyc,1,0,4,5,2,3,7,6);
				break;
			  case 7:  /* Patern 12:0-1-5-7 */
                                F3_Patern5(cnode,cyc,0,3,1,2,5,6);
                                F3_Patern5(cnode,cyc,0,3,5,6,0,4);
                                F3_Patern5(cnode,cyc,0,4,5,6,4,5);
                                SetCornerOn(cnode,cyc,7);
				break;
			}
			break;
		  case 6:
			switch(pnode[3]) {
			  case 7:  /* patern 10:0-1-6-7 */
                                F2_Patern2_1(cnode,cyc,0,1);
                                F2_Patern2_2(cnode,cyc,6,7);
				break;
			}
			break;
		  case 7:
			break;
		}/* end of switch(pnode[2]) */
		break;
	  case 2:
		switch(pnode[2]) {
		  case 3:
			switch(pnode[3]) {
			  case 4:  /* Patern 14:0-2-3-4 */
                                F4_Patern14(cnode,cyc,0,1,2,3,4,5,6,7,8);
				break;
			  case 5:  /* Patern 12:0-2-3-5 */
                                F3_Patern5(cnode,cyc,0,4,3,7,2,6);
                                F3_Patern5(cnode,cyc,0,4,2,6,0,1);
                                F3_Patern5(cnode,cyc,0,1,2,6,1,2);
                                SetCornerOn(cnode,cyc,5);
				break;
			  case 6:  /* Patern 11:0-2-3-6 */
                                F4_Patern11(cnode,cyc,3,0,1,2,7,4,5,6);
				break;
			  case 7:  /* patern 9:0-2-3-7 */
                                F4_Patern9(cnode,cyc,0,1,2,7,4,6);
				break;
			}
			break;
		  case 4:
			switch(pnode[3]) {
			  case 5:  /* Patern 12:0-2-4-5 */
                                F3_Patern5(cnode,cyc,5,6,4,7,0,3);
                                F3_Patern5(cnode,cyc,5,6,0,3,1,5);
                                F3_Patern5(cnode,cyc,1,5,0,3,0,1);
                                SetCornerOn(cnode,cyc,2);
				break;
			  case 6:  /* patern 10:0-2-4-6 */
                                F2_Patern2_3(cnode,cyc,0,4);
                                F2_Patern2_3(cnode,cyc,2,6);
				break;
			  case 7:  /* Patern 12:0-2-4-7 */
                                F3_Patern5(cnode,cyc,0,1,4,5,6,7);
                                F3_Patern5(cnode,cyc,0,1,6,7,0,3);
                                F3_Patern5(cnode,cyc,0,3,6,7,3,7);
                                SetCornerOn(cnode,cyc,2);
				break;
			}
			break;
		  case 5:
			switch(pnode[3]) {
			  case 7:  /* patern 13:0-2-5-7 */
                                F4_Patern13(cnode,cyc,pnode);
				break;
			}
			break;
		}/* end of switch(pnode[2]) */
		break;
	  case 3:
		switch(pnode[2]) {
		  case 4:
			switch(pnode[3]) {
			  case 5:  /* Patern 11:0-3-4-5 */
                                F4_Patern11(cnode,cyc,0,3,7,4,1,2,6,5);
				break;
			  case 6:  /* Patern 12:0-3-4-6 */
                                F3_Patern5(cnode,cyc,2,3,0,1,4,5);
                                F3_Patern5(cnode,cyc,2,3,4,5,3,7);
                                F3_Patern5(cnode,cyc,3,7,4,5,4,7);
                                SetCornerOn(cnode,cyc,6);
				break;
			  case 7:  /* patern 8:0-3-4-7 */
                                F3_Patern5(cnode,cyc,4,5,6,7,0,1);
                                F3_Patern5(cnode,cyc,0,1,6,7,2,3);
				break;
			}
			break;
		  case 5:
			switch(pnode[3]) {
			  case 6:  /* patern 10:0-3-5-6 */
                                F2_Patern2_1(cnode,cyc,3,0);
                                F2_Patern2_2(cnode,cyc,5,6);
				break;
			  case 7:  /* Patern 12:0-3-5-7 */
                                F3_Patern5(cnode,cyc,6,7,2,3,0,1);
                                F3_Patern5(cnode,cyc,6,7,0,1,4,7);
                                F3_Patern5(cnode,cyc,4,7,0,1,0,4);
                                SetCornerOn(cnode,cyc,5);
				break;
			}
			break;
		  case 6:
			switch(pnode[3]) {
			  case 7:  /* Patern 14:0-3-6-7 */
                                F4_Patern14(cnode,cyc,3,2,6,7,0,1,5,4);
				break;
			}
			break;
		}/* end of switch(pnode[2]) */
		break;
	  case 4:
		switch(pnode[2]) {
		  case 5:
			switch(pnode[3]) {
			  case 6:  /* Patern 14:0-4-5-6 */
                                F4_Patern14(cnode,cyc,4,7,6,5,0,3,2,1);
				break;
			  case 7:  /* patern 9:0-4-5-7 */
                                F4_Patern9(cnode,cyc,7,6,5,0,3,1);
				break;
			}
			break;
		  case 6:
			switch(pnode[3]) {
			  case 7:  /* Patern 11:0-4-6-7 */
                                F4_Patern11(cnode,cyc,7,6,5,4,3,2,1,0);
				break;
			}
			break;
		}/* end of switch(pnode[2]) */
		break;
	  case 5:
		switch(pnode[2]) {
		  case 6:
			break;
		}/* end of switch(pnode[2]) */
		break;
	} /* end of switch(pnode[1]) */
}

void SelectCell(stpoint)
int stpoint;
{
  int     node[8], cnode[8][8];
  int     idz,i,j;
  int     flag;

  node[0] = stpoint;
  node[1] = node[0]+1;
  node[2] = node[0]+idmx+1;
  node[3] = node[0]+idmx;
  idz = idmx*idmy;
  node[4] = node[0]+idz;
  node[5] = node[1]+idz;
  node[6] = node[2]+idz;
  node[7] = node[3]+idz;

  flag = 0;
  for (i=0;i<8;i++){
    flag += ids[node[i]].flag;
  }
  if (flag == 0) return;

  /*
  // nodeのindexをcyclicにcnodeにセット
  */
  for (i=0;i<4;i++){
    for (j=0;j<4;j++){
      cnode[i][j] = node[(i+j)%4];
    }
    for (j=4;j<8;j++){
      cnode[i][j] = node[4+(i+j)%4];
    }
  }
  for (i=4;i<8;i++){
    for (j=0;j<4;j++){
      cnode[i][j] = node[4+(i+8-j)%4];
    }
    for (j=4;j<8;j++){
      cnode[i][j] = node[(i+8-j)%4];
    }
  }

  switch(flag) {
  case 1: Flag1(cnode); break;
  case 2: Flag2(cnode); break;
  case 3: Flag3(cnode); break;
  case 4: Flag4(cnode); break;
  case 5: Flag5(cnode); break;
  case 6: Flag6(cnode); break;
  case 7: Flag7(cnode); break;
  }
}

void MakeObject(sw)
int  sw;
{
  int     i, j, ct;
  PGD     *p;

  if (sw == 0){ p = pgd;  ct = pgdct;  }
  if (sw == 1){ p = pgd2; ct = pgdct2; }

  fprintf(out, "DEF Polygon_%d Separator {\n", sw + 1);
  fprintf(out, " Material {\n");
  if (sw == 0) fprintf(out, "  diffuseColor 1.0 0.3 0.3\n");
  else         fprintf(out, "  diffuseColor 0.3 1.0 0.3\n");
  fprintf(out, "  specularColor 0.2 0.2 0.2\n");
  fprintf(out, "  shininess 0.2\n");
  fprintf(out, "  transparency 0.5\n");
  fprintf(out, " }\n");
  fprintf(out, " Coordinate3 {\n");
  fprintf(out, "  point [\n");
  printf("Output VRML polygon ");
  for (i = 0; i < ct; i++){
    if (i % 100 == 0){ printf("[%d]", i); fflush(stdout); }
    for (j = 0; j < 3; j++){
      fprintf(out, "   %+5.3f %+5.3f %+5.3f",p[i].x[j],p[i].y[j],p[i].z[j]);
      if (i == ct - 1 && j == 2) fprintf(out, "\n");
      else                       fprintf(out, ",\n");
    }
  }
  fprintf(out, "  ]\n");
  fprintf(out, " }\n");
  printf("\n");

  fprintf(out, " IndexedFaceSet {\n");
  fprintf(out, "  coordIndex [\n");
  for (i = 0; i < ct; i++){
    fprintf(out, "   %d, %d, %d, -1", i*3, i*3+1, i*3+2);
    if (i == ct - 1) fprintf(out, "\n");
    else             fprintf(out, ",\n");
  }
  fprintf(out, "  ]\n");
  fprintf(out, " }\n");
  fprintf(out, "}\n");
}

void MolColor(sw)
int  sw;
{
  int    x;

  if (sw < 22) x = sw;
  else         x =  0;
  fprintf(out, "  Material {\n");
  fprintf(out, "   diffuseColor %4.2f %4.2f %4.2f\n",
               mc[x][3],mc[x][4],mc[x][5]);
  fprintf(out, "  }\n");
}

void MakeBond()
{
  int    i, j, err;
  float  d, d2, d2h, rw, rd, r;
  float  xd, yd, zd, dd, lg, lg2;
  float  xa, ya, za;
  int    x, y, c, fg;
  char   buf[256];

  d   = ADIS;
  d2  = d * d;
  d2h = ADISH * ADISH;
  rd  = 180.0 / PI;
  rw  = RATOM * WATOM;
  lg2 = LATOM;
  printf("Make BOND\n");
  fprintf(out, "DEF ATOM Separator {\n");
  for (i = 0; i < molct; i++){
    fprintf(out, " DEF ATOM_No%d_Mol%d Separator {\n", i+1, mol[i].mn);
    fprintf(out, "  Translation { translation %+5.3f %+5.3f %+5.3f }\n",
                 mol[i].x, mol[i].y, mol[i].z);
    MolColor(mol[i].mn);
    if (mol[i].mn == 1) r = RATOM * HATOM;
    else                r = RATOM;
    fprintf(out, "  Sphere { radius %+5.3f }\n", r);
    fprintf(out, " }\n");
  }

  for (i = 0; i < molct; i++){
    fprintf(out, " DEF ATOM_LINE_No%d Separator {\n", i+1);
    MolColor(1);
    for (j = i + 1; j < molct; j++){
      xd = mol[j].x - mol[i].x;
      if (fabs(xd) > d) continue;
      yd = mol[j].y - mol[i].y;
      if (fabs(yd) > d) continue;
      zd = mol[j].z - mol[i].z;
      if (fabs(zd) > d) continue;
      dd = xd * xd + yd * yd + zd * zd;
      if (dd > d2) continue;
      if (mol[i].mn == 1 && mol[j].mn == 1 && dd > d2h) continue;
      lg = sqrt(dd);
      printf("Line %d-%d  Length=%+f\n", i, j, lg);
      xa = 0.0; ya = 0.0; za = 0.0;
      if (xd == 0.0){
        if (yd == 0.0){
          if (zd > 0.0) xa =  PI2;
          else          xa = -PI2;
        }
        else if (zd == 0.0){
          if (yd > 0.0) ;
          else xa = PI;
        }
        else xa = (float)atan2((double)zd, (double)yd);
      }
      else if (yd == 0.0){                                 
        if (zd == 0.0){
          if (xd > 0.0) za = -PI2;
          else          za =  PI2;
        }
        else {
          xa =  PI2;
          za = (float)atan2((double)zd, (double)xd) - PI2;
        }
      }
      else if (zd == 0.0){
        za = (float)atan2((double)yd, (double)xd) - PI2;
      }
      else {
        xa = (float)atan2((double)zd, (double)yd);
        za = (float)acos((double)(xd / lg)) - PI2;
      }
      fprintf(out, "  DEF ATOM_LINE_No%d_No%d Separator {\n", i+1, j+1);
      fprintf(out, "   Translation { translation %+5.3f %+5.3f %+5.3f }\n",
                   mol[i].x,mol[i].y,mol[i].z);
      if (xa != 0.0) fprintf(out, "   Rotation { rotation 1 0 0 %+5.3f }\n", xa);
      if (za != 0.0) fprintf(out, "   Rotation { rotation 0 0 1 %+5.3f }\n", za);
      fprintf(out, "   Translation { translation 0 %+5.3f 0 }\n", lg / 2.0);
      fprintf(out, "   Cylinder {\n");
      fprintf(out, "    parts  ALL\n");
      fprintf(out, "    radius %+5.3f\n", rw);
      fprintf(out, "    height %+5.3f\n", lg);
      fprintf(out, "   }\n");
      fprintf(out, "  }\n");
    }
    fprintf(out, " }\n");
  }  
  fprintf(out, "}\n");
}                      

void MakePolygon()
{
  if (pgv[0] == 1)            MakeObject(0);
  if (pgs > 1 && pgv[1] == 1) MakeObject(1);
}

void MakeParts()
{
/*    MakeAxis();*/
  MakeBond();
  MakePolygon();
}

MakeVRML()
{
    int     winx=800, winy=800;
    int     i, j, k;
    long    bg, ed;
    char    filename[256];
    char    buf[256];
    FILE    *fp;

    for (i = 0; i < PGMAX; i++){ pgv[i] = 1; SOLID[i] = 1; }
/*
    if (argv[i][0] == '-' || argv[i][0] == '/'){
            switch(toupper(argv[i][1])){
            case 'W': for (i = 0; i < PGMAX; i++) SOLID[i] = 0;
                      break;
            case 'S': SHADE = 0; break;
            case 'F': FOG   = 0; break;
            case 'A': AXIS  = 0; break;
            case 'U': ALIAS = 0; break;
            case 'C':
                if ('0' <= argv[i][2] && argv[i][2] <= '9') COLORID = argv[i][2] - '0';
                break;
            case 'Q':
                if ('0' <= argv[i][2] && argv[i][2] <= '9') pgv[argv[i][2] - '1'] = 0;
                break;
            }
        }
*/
    printf("-----------------------\n");
    printf("Make VRML\n");
    printf("Moleculer count = %d\n", molct);
    for (i = 0; i < molct; i++)
        printf("%d: %d  X=%+e  Y=%+e  Z=%+e\n",i+1,mol[i].mn,mol[i].x,mol[i].y,mol[i].z);
                                 
    idmx    = pl[2];
    idmy    = pl[3];
    idmz    = pl[4];
    xbig    = idmx / 2;
    ybig    = idmz / 2;
    zbig    = idmz / 2;

    mbgx    = pf[0];
    mbgy    = pf[1];
    mbgz    = pf[2];
    mstx    = pf[3];
    msty    = pf[4];
    mstz    = pf[5];
    xmax    = pf[6];
    xmin    = pf[7];
    ymax    = pf[8];
    ymin    = pf[9];
    zmax    = pf[10];
    zmin    = pf[11];
    datamax = pf[12];
    datamin = pf[13];

    wkx = xmax - xmin;
    wky = ymax - ymin;
    wkz = zmax - zmin;
    wkmax = wkx;
    if(wky > wkmax) wkmax = wky;
    if(wkz > wkmax) wkmax = wkz;
    wkmax = (float)0.5 * (float)sqrt(2*wkmax*wkmax);
/*    printf("wkmax=%+f \n",wkmax);*/

    dist = wkmax * 1.2;
    eyex = dist / INIX;
    eyey = dist / INIY;
    eyez = dist / INIZ;
    divm = dist / 50.0;

    out = fopen(fileout, "wt");
    fprintf(out, "#VRML V1.0 ascii\n");
    fprintf(out, "DEF Title Info {\n");
    fprintf(out, " string \"MO display\"\n");
    fprintf(out, "}\n");
    fprintf(out, "DEF ViewerSpeed Info {\n");
    fprintf(out, " string \"10\"\n");
    fprintf(out, "}\n");
/*    fprintf(out, "DirectionalLight {\n");
    fprintf(out, " on TRUE\n");
    fprintf(out, " intensity 1\n");
    fprintf(out, " color 1 1 1\n");
    fprintf(out, " direction -1 -1 -1\n");
    fprintf(out, " location 100.0 100.0 100.0\n");
    fprintf(out, "}\n");
    fprintf(out, "DirectionalLight {\n");
    fprintf(out, " on TRUE\n");
    fprintf(out, " intensity 0.7\n");
    fprintf(out, " color 1 1 1\n");
    fprintf(out, " direction 1 1 1\n");
    fprintf(out, " location -100.0 -100.0 -100.0\n");
    fprintf(out, "}\n");*/
    MakeParts();
}   

/*
// input argument
//  *argv[0] : program name
//  *argv[1] : input  file name
//  *argv[2] : iso value
*/

void main(argc, argv)
int	argc;
char	**argv;
{
    int     winx=800, winy=800;
    int     i, j, k;
    int     db, ac;
    char    filein[1024];
    FILE    *fi;
    float   dt[18];

    printf("\n");
    printf("D2VRML  Copyright(C) 1996  Computer Center of IMS\n\n");
    db = 0;
    if (argc != 3 && argc != 4) {
        printf("USAGE: D2VRML infile iso_value\n");
        exit(1);
    }
    strcpy(filein,  argv[1]);
    sscanf(argv[2],"%f",&vvvv);
    viso = vvvv;
    ac = 3;
    while (1){
        if (ac >= argc) break;
        if (argv[ac][0] == '-' || argv[ac][0] == '/'){
            if (toupper(argv[ac][1]) == 'D') dble = 2;
	    if (toupper(argv[ac][1]) == 'T') dble = 3;
        }
        ac++;
    }
    strcpy(fileout, filein);
    for (i = strlen(fileout) - 1; i >= 0; i--){
	if (fileout[i] == '.') break;
    }
    if (i > 0){
	fileout[i] = 0x00;
	strcat(fileout, ".wrl");
    }
    printf("argc= %d infile = %s outfile = %s iso_value=%+f\n",
            argc, filein, fileout, viso);

    if((fi = fopen(filein, "rb"))==NULL){
            perror("fopen error ");
            exit(1);
    }
    InputData(fi);
    fclose(fi);

    printf("-----------------------\n");
    printf("Make polygon data\n");
    pgdmode = 0;
    FlagSelect();
    printf("Z loop ");
    for (k = 0; k < idmz - 1; k++){
        printf("[%d]", k + 1);
	fflush(stdout);
        for (j = 0; j < idmy - 1; j++){
            for (i = 0; i < idmx - 1; i++){
                SelectCell(i + idmx * (j + idmy * k));
            }
        }
    }
    printf("\nPolygon total = %d\n", pgdct);
    printf("Make normal vector\n");
    NormalizeVector();

    pl[0]  = (long)GCODE;  /* program code  */
    pl[1]  = dtype;  /* data type     */
    pl[2]  = idmx;   /* x resolution  */
    pl[3]  = idmy;   /* y resolution  */
    pl[4]  = idmz;   /* z resolution  */
    pl[5]  = 1;      /* object count  */
    pl[6]  = oct;    /* orbital count */
    pl[7]  = molct;  /* moleculer count */
    pl[14] = obt[0];
    pl[15] = obt[1];
    pl[16] = obt[2];
    pl[17] = obt[3];
    pl[18] = obt[4];
    pl[19] = obt[5];
    pl[20] = pgdct;  /* polygon count 1 */

    pf[0]  = mbgx;
    pf[1]  = mbgy;
    pf[2]  = mbgz;
    pf[3]  = mstx;
    pf[4]  = msty;
    pf[5]  = mstz;
    pf[6]  = xmax;
    pf[7]  = xmin;
    pf[8]  = ymax;
    pf[9]  = ymin;
    pf[10] = zmax;
    pf[11] = zmin;
    pf[12] = datamax;
    pf[13] = datamin;
    pf[20] = viso;     /* iso value 1 */

    if (DCODE == DOSCOD) cmode = FALSE;
    else                 cmode = TRUE;
    if (db == 0 && dtype == 1) db = 1;
    if (db == 1) printf("Double data mode\n");
    if (cmode)   printf("Swap byte mode\n");
    else         printf("Not swap byte mode\n");

    if (db == 1){
	pgdmode = 1;
        viso = -vvvv;
        printf("-----------------------\n");
        printf("Make polygon data No.2\n");
	printf("isovalue = %+f\n", viso);
        FlagSelect();
	printf("Z loop ");
        for (k = 0; k < idmz - 1; k++){
            printf("[%d]", k + 1);
	    fflush(stdout);
            for (j = 0; j < idmy - 1; j++){
                for (i = 0; i < idmx - 1; i++){
                    SelectCell(i + idmx * (j + idmy * k));
                }
            }
        }
        printf("\nNo.2 polygon total = %d\n", pgdct2);
        printf("Make normal vector No.2\n");
        NormalizeVector();
        pl[5]++;         /* object count     */
        pl[21] = pgdct;  /* polygon count 1  */
        pf[21] = viso;   /* iso value 1      */
    }
    pgs = pl[5];
    MakeVRML();
}