#include "shape.h"

int
trifac(VERTEX *ifp, int nrows)
{
	// Procedure to form the vertice triplets that make up a triangulation facet
	int i;
	int j, j0, j1, j2, j3;
	int nnode;
	int ntriplets;
	int node[2 * nrows][4 * nrows + 1];

	// Number each vertice
	for (i = 0; i < 2 * nrows; i++) {
		for (j = 0; j < 4 * nrows + 1; j++)
			node[i][j] = 0;
	}
	nnode = 0;
	node[0][0] = nnode;
	for (i = 1; i <= nrows; i++) {
		for (j = 0; j < 4 * i; j++) {
			nnode++;
			node[i][j] = nnode;
			if (j == 0)
				node[i][4 * i] = nnode;
		}
	}
	for (i = nrows - 1; i > 0; i--) {
		for (j = 0; j < 4 * i; j++) {
			nnode++;
			node[2 * nrows - i][j] = nnode;
			if (j == 0)
				node[2 * nrows - i][4 * i] = nnode;
		}
	}
	node[2 * nrows][0] = nnode + 1;

    //if (debug) {
    //    for (i = 0; i < 2 * nrows; i++) {
    //        printf("#");
    //        for (j = 0; j < 4 * nrows + 1; j++)
    //            printf(" %3d ", node[i][j]);
    //        printf("\n");
    //    }
    //}

	// Associate vertices to form facets
	j2 = 0;
	ntriplets = 0;
	for (j1 = 1; j1 <= nrows; j1++) {
		for (j3 = 1; j3 <= 4; j3++) {
			j0 = (j3 - 1) * j1;
			ifp[ntriplets].v[0] = node[j1 - 1][j0 - (j3 - 1)];
			ifp[ntriplets].v[1] = node[j1][j0];
			ifp[ntriplets].v[2] = node[j1][j0 + 1];

            if (debug) {
                printf("#a-- %3d %3d %3d %3d -- %3d -- %3d %3d %3d \n", j0, j1, j2, j3,
                             ntriplets, ifp[ntriplets].v[0], ifp[ntriplets].v[1], ifp[ntriplets].v[2]);
            }

			ntriplets++;
			for (j2 = (j0 + 1); j2 <= (j0 + j1 - 1); j2++) {
				ifp[ntriplets].v[0] = node[j1][j2];
				ifp[ntriplets].v[1] = node[j1 - 1][j2 - (j3 - 1)];
				ifp[ntriplets].v[2] = node[j1 - 1][(j2 - 1) - (j3 - 1)];

                if (debug) {
                    printf("#b-- %3d %3d %3d %3d -- %3d -- %3d %3d %3d \n", j0, j1, j2,
                                 j3, ntriplets, ifp[ntriplets].v[0], ifp[ntriplets].v[1], ifp[ntriplets].v[2]);
                }

				ntriplets++;
				ifp[ntriplets].v[0] = node[j1 - 1][j2 - (j3 - 1)];
				ifp[ntriplets].v[1] = node[j1][j2];
				ifp[ntriplets].v[2] = node[j1][j2 + 1];

                if (debug) {
                    printf("#c-- %3d %3d %3d %3d -- %3d -- %3d %3d %3d \n", j0, j1, j2,
                                 j3, ntriplets, ifp[ntriplets].v[0], ifp[ntriplets].v[1], ifp[ntriplets].v[2]);
                }

				ntriplets++;
			}
		}
	}
	j2 = 0;
	// Create the facets for the lower hemisphere
	for (j1 = nrows + 1; j1 <= 2 * nrows; j1++) {
		for (j3 = 1; j3 <= 4; j3++) {
			j0 = (j3 - 1) * (2 * nrows - j1);
			ifp[ntriplets].v[0] = node[j1][j0];
			ifp[ntriplets].v[1] = node[j1 - 1][(j0 + 1) + (j3 - 1)];
			ifp[ntriplets].v[2] = node[j1 - 1][j0 + (j3 - 1)];

            if (debug) {
                printf("#d-- %3d %3d %3d %3d -- %3d -- %3d %3d %3d \n", j0, j1, j2, j3,
                             ntriplets, ifp[ntriplets].v[0], ifp[ntriplets].v[1], ifp[ntriplets].v[2]);
            }

			ntriplets++;
			for (j2 = (j0 + 1); j2 <= (j0 + 2 * nrows - j1); j2++) {
				ifp[ntriplets].v[0] = node[j1][j2];
				ifp[ntriplets].v[1] = node[j1 - 1][j2 + (j3 - 1)];
				ifp[ntriplets].v[2] = node[j1][j2 - 1];

                if (debug) {
                    printf("#e-- %3d %3d %3d %3d -- %3d -- %3d %3d %3d \n", j0, j1, j2,
                                 j3, ntriplets, ifp[ntriplets].v[0], ifp[ntriplets].v[1], ifp[ntriplets].v[2]);
                }

				ntriplets++;
				ifp[ntriplets].v[0] = node[j1][j2];
				ifp[ntriplets].v[1] = node[j1 - 1][(j2 + 1) + (j3 - 1)];
				ifp[ntriplets].v[2] = node[j1 - 1][j2 + (j3 - 1)];

                if (debug) {
                    printf("#f-- %3d %3d %3d %3d -- %3d -- %3d %3d %3d \n", j0, j1, j2,
                                 j3, ntriplets, ifp[ntriplets].v[0], ifp[ntriplets].v[1], ifp[ntriplets].v[2]);
                }

				ntriplets++;

			}
		}
	}
	return (ntriplets);
}
