/* Program to generate wireframe model of 3D objects
   Niels Lund, September 1998
   Modified for parameter file input 'net_work.txt'. August 2004
*/

/* Each object in the model is defined by a number of generator lines
   which are rotated around an axis to form the object.
   The generator lines can consist of up to 9 straight line segments
   defined by <imax> points in the arrays h[10] and ra[10].
   Cylinders and cones are produced by using many generators (defined
   by <nmax>), the subroutine call used is:
              pointnum = calpoint(nmax, imax, pointnum);
   (pointnum is a running index for the points in the wireframe model)

   A simple box can be described by a generator consisting of a single
   straight line defined by two points (imax=2), repeated four times
   (nmax=4). To make boxes with rectangular cross section use the
   subroutine:
              pointnum = calrect(dx, dy, pointnum);


   The resulting wireframe model is output to the file: 'wirefram.dat'
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

int calpoint(int, int, int);
int calrect(double, double, int);
double inerti(int, int, int);
int datain(int *type, int *imax, double h[], double r[], double *z_mount, double *r_mount,
   double *wgt, double *z_cg, int *nmax_i);

FILE *out1, *out2, *mhfil;
double dir_matrix[3][3], displ[3], xw[3], xb[2000][3], h[12], ra[12];
double geo_pos[20][3], geo_size[20][3], weights[20], cg[3];
double pi = 3.141593;
int pointnum, linenum, facenum, elem;

main(int argc, char *argv[])
{
  double dx, dy, mkg, msum, msum_x, msum_y, msum_z, z_cg, z_mount, r_mount;
  double xmnt, ymnt, zmnt, Ixx, Iyy, Izz, ri[12];
  double mean_r;
  int i, nmax, imax, nmax_c, status, type, nmax_i;
  char s[256], alfil[200];

  pointnum = 1;

  linenum = 0;
  facenum = 0;
  elem = 0;

  msum = msum_x = msum_y = msum_z = 0.0L;
  Ixx = Iyy = Izz = 0.0L;


  nmax_c = 15;
  if (argc>2) nmax_c = atoi(argv[2]);
  if (nmax_c > 90) nmax_c = 90;

  out1 = fopen("wirefram.dat","wt") ;	/* Output file */
  if ( out1 == NULL )
  {
     printf("File opening error on output file 'wireframe.dat'!\n") ;
     exit(0) ;
  }
  fprintf(out1, "# Wireframe model for 'ThreeD.class' JAVA program\n");
  fprintf(out1, "# Niels Lund. August 2004\n");

  out2 = fopen("temp.lin","wt") ; /* temporary file */
  if ( out1 == NULL )
  {
     printf("File opening error on temporary file\n") ;
     exit(0) ;
  }

  if (argc>1) {
  	mhfil = fopen(argv[1], "rt");
   strcpy(alfil, argv[1]);
  } else {
      sprintf(alfil, "net_gen.txt");
  		mhfil = fopen("net_gen.txt", "rt");
  }
  if (mhfil == NULL) {
     printf("No input file found! (%s)\n", alfil);
     exit(0);
  }

next_elem:
  nmax_i = 0;
  status = datain(&type, &imax, h, ri, &z_mount, &r_mount, &weights[elem], &z_cg, &nmax_i);

  printf("status: %2d\n", status);

  if (status < 0) return(0);
  if (status > 1) goto complete;

  printf("%2d  type: %1d, imax: %2d, z_mount: %5.1lf, r_mount: %5.1lf, weight: %5.1lf, z_cg: %5.1lf\n",
      elem, type, imax, z_mount, r_mount, weights[elem], z_cg);
  printf("%2d  displacement:   %7.4lf  %7.4lf  %7.4lf\n", elem, displ[0], displ[1], displ[2]);
  printf("%2d  X dir. cosines:   %7.4lf  %7.4lf  %7.4lf\n",
     elem, dir_matrix[0][0], dir_matrix[0][1],	dir_matrix[0][2]);
  printf("%2d  Y dir. cosines:   %7.4lf  %7.4lf  %7.4lf\n",
     elem, dir_matrix[1][0], dir_matrix[1][1], dir_matrix[1][2]);
  printf("%2d  Z dir. cosines:   %7.4lf  %7.4lf  %7.4lf\n",
     elem, dir_matrix[2][0], dir_matrix[2][1], dir_matrix[2][2]);

  geo_pos[elem][0] = displ[0] + z_cg * dir_matrix[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * dir_matrix[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * dir_matrix[2][2];

  mkg = weights[elem];

  if (type == 1) {	/* cylindrical elements */
     if (nmax_i > 0) nmax = nmax_i; else nmax = nmax_c;
     mean_r = 0.0;
     geo_size[elem][0] = abs(h[0] - h[imax-1]) / 2.0;
     for (i=1; i<imax; i++) mean_r += (ri[i-1] + ri[i]) * abs(h[i-1] - h[i]);
     geo_size[elem][1] = mean_r * 0.5 / abs(h[0] - h[imax-1]);
     geo_size[elem][2] = geo_size[elem][1];

     for (i=0; i<12; i++) ra[i] = ri[i];

     elem++;

     pointnum = calpoint(nmax, imax, pointnum);

  } else {		/* box elements */

     nmax = 4;
     dx = ri[0] / 2.0;
     dy = ri[1] / 2.0;

     ra[0] = sqrt(dx*dx + dy*dy);
     ra[1] = sqrt(dx*dx + dy*dy);

     geo_size[elem][0] = dx;
     geo_size[elem][1] = dy;
     geo_size[elem][2] = abs(h[0] - h[1]) * 0.5;

     elem++;

     pointnum = calrect(dx, dy, pointnum);
  }

  msum += mkg;
  msum_x += mkg * (displ[0] + dir_matrix[0][2] * z_cg);
  msum_y += mkg * (displ[1] + dir_matrix[1][2] * z_cg);
  msum_z += mkg * (displ[2] + dir_matrix[2][2] * z_cg);


  linenum += nmax * imax;

  xmnt = displ[0] + dir_matrix[0][2] * z_mount;
  ymnt = displ[1] + dir_matrix[1][2] * z_mount;
  zmnt = displ[2] + dir_matrix[2][2] * z_mount;

  printf("\nElement %2d.   Mount center (x,y,z cm): %4.1lf, %4.1lf, %4.1lf;  Radius: %3.1lf",
	elem, xmnt, ymnt, zmnt, r_mount);
  printf("\nElement %2d.   Mass %4.1lf kg, CG-Mount sep: %4.1lf cm,  Vector: %4.3lf, %4.3lf, %4.3lf\n",
	elem, mkg, z_cg-z_mount, dir_matrix[0][2], dir_matrix[1][2], dir_matrix[2][2]);

  goto next_elem;

  /* *************************** */
  /* Summary print */
complete:
  printf("\n points: %4d, line segments: %4d, faces: %3d\n",
             pointnum, linenum, facenum);

  cg[0] = msum_x/msum;
  cg[1] = msum_y/msum;
  cg[2] = msum_z/msum;

  printf("\n Total mass: %4.1lf kg, CG-position (x,y,z cm): %4.1lf, %4.1lf, %4.1lf",
	     msum, cg[0], cg[1], cg[2]);

  for (i=0; i<elem; i++) {
	Ixx = Ixx + inerti(1, 2, i);
	Iyy = Iyy + inerti(0, 2, i);
	Izz = Izz + inerti(0, 1, i);
  }

  printf("\n Moments of inertia: Ixx: %4.1lf, Iyy: %4.1lf, Izz: %4.1lf\n",
	     Ixx, Iyy, Izz);

  fclose(out2);

  out2 = fopen("temp.lin","rt") ;
  if ( out2 == NULL )
  {
	printf("File opening error!!!\n") ;
	exit(0) ;
  }

  while (fgets(s, imax, out2) != NULL) fputs(s, out1);

  fcloseall();
  return(0);
}

/* **************************************************************** */

int calpoint(int nmax, int imax, int pointnum)
{
  double angle, delta, cang, sang;
  int n, i, m, k, mmax;
  m = 0;
  mmax = nmax * imax;

  delta = 2.0 * pi / nmax;
  angle = delta / 2.0L;

  for (n=0; n<nmax; n++) {
    sang = sin(angle);
    cang = cos(angle);
    for (i=0; i<imax; i++) {
      xb[m][0] = ra[i] * cang;
      xb[m][1] = ra[i] * sang;
      xb[m][2] = h[i];
/*      printf("\nv %8.3lf %8.3lf %8.3lf", xb[m][0], xb[m][1], xb[m][2]); */
      m++;
    }
    angle += delta;
  }
/* transform individual points to new positions corresponding to each WATCH */

  for (m=0; m<mmax; m++) {
    for ( i = 0 ; i < 3 ; i++ ) {
     	xw[i] = 0.0 ;
      for ( k = 0 ; k < 3 ; k++ ) xw[i] += dir_matrix[i][k]*xb[m][k];
    }
    fprintf(out1, "\nv %8.3lf %8.3lf %8.3lf",
            xw[0]+displ[0], xw[1]+displ[1], xw[2]+displ[2]);
  }

/* generate list of line segments and faces */

  for (n=0; n<nmax; n++) {  /* generate line elements along "frembringere" */
    for (i=1; i<imax; i++) {
      m = pointnum + n * imax + i - 1;
      fprintf(out2, "\nl %4d %4d", m, m+1);
      linenum++;
    }
  }

  for (i=0; i<imax; i++) {  /* generate "face" polygons */
    m = pointnum + i;
    fprintf(out2, "\nf");
    for (n=0; n<nmax; n++) fprintf(out2, " %4d", m + n*imax);
    facenum++;
  }

  pointnum += nmax * imax;

  return(pointnum);
}

/* *********************************************************** */

int calrect(double dx, double dy, int pointnum)
{
  double angle, cang, sang;
  int n, i, m, k, mmax, nmax, imax;
  m = 0;
  nmax = 4;
  imax = 2;
  mmax = nmax * imax;

  angle = atan2(dy, dx);

  sang = sin(angle);
  cang = cos(angle);

  for (i=0; i<imax; i++) {
    xb[m][0] = ra[i] * cang;
    xb[m][1] = ra[i] * sang;
    xb[m][2] = h[i];
    m++;
  }

  for (i=0; i<imax; i++) {
    xb[m][0] = ra[i] * cang;
    xb[m][1] = -ra[i] * sang;
    xb[m][2] = h[i];
    m++;
  }

  for (i=0; i<imax; i++) {
    xb[m][0] = -ra[i] * cang;
    xb[m][1] = -ra[i] * sang;
    xb[m][2] = h[i];
    m++;
  }

  for (i=0; i<imax; i++) {
    xb[m][0] = -ra[i] * cang;
    xb[m][1] = ra[i] * sang;
    xb[m][2] = h[i];
    m++;
  }

/* transform individual points to new positions corresponding to each WATCH */

  for (m=0; m<mmax; m++) {
    for ( i = 0 ; i < 3 ; i++ ) {
     	xw[i] = 0.0 ;
      for ( k = 0 ; k < 3 ; k++ ) xw[i] += dir_matrix[i][k]*xb[m][k];
    }
    fprintf(out1, "\nv %8.3lf %8.3lf %8.3lf",
            xw[0]+displ[0], xw[1]+displ[1], xw[2]+displ[2]);
  }

/* generate list of line segments and faces */

  for (n=0; n<nmax; n++) {  /* generate line elements along "frembringere" */
    for (i=1; i<imax; i++) {
      m = pointnum + n * imax + i - 1;
      fprintf(out2, "\nl %4d %4d", m, m+1);
      linenum++;
    }
  }

  for (i=0; i<imax; i++) {  /* generate "face" polygons */
    m = pointnum + i;
    fprintf(out2, "\nf");
    for (n=0; n<nmax; n++) fprintf(out2, " %4d", m + n*imax);
    facenum++;
  }

  pointnum += nmax * imax;

  return(pointnum);
}

/* *********************************************************** */

double inerti(int k, int m, int elem)
{
	double iner;

	iner = weights[elem] *
		(geo_size[elem][k]*geo_size[elem][k]
		 + geo_size[elem][m]*geo_size[elem][m]) / 2.0;

	iner += weights[elem] *
		((geo_pos[elem][k]-cg[k])*(geo_pos[elem][k]-cg[k])
		 + (geo_pos[elem][m]-cg[m])*(geo_pos[elem][m]-cg[m]));

	return(iner);
}

/* ***************************************************************************************** */

int datain(int *type, int *imax, double h[], double r[], double *z_mount, double *r_mount,
   double *wgt, double *z_cg, int *nmax_i)
{
      char mhline[410];
      int i, n;
      double dum0, dum1, dum2;

      for (i=0; i<200; i++) {
            fgets(mhline, 200, mhfil);
            if (feof(mhfil)) goto q_end;
            sscanf(mhline, "%d", &n);
            printf("n=%d, inline: %s", n, mhline);
            if (n==1) {
               sscanf(mhline, "%d %d %d", &n, type, imax);
               if ((*type == 2) && (*imax != 2)) {
                  printf("Bad input data: (type = 2 but no_of_segment endpoints != 2)  %s\n", mhline);
                  return(-1);
               }
               if ((*type < 1) || (*type > 2)) {
                  printf("Bad input data: (type not 1[cylinder] or 2[box])  %s\n", mhline); 
                  return(-1);
               }
               if (*imax > 12) {
                  printf("Bad input data: (imax > 12)  %s\n", mhline); 
                  return(-1);
               }
            }
            if (n==2) {sscanf(mhline, "%d %lf %lf", &n, z_mount, r_mount); }
            if (n==3) {
               sscanf(mhline, "%d %lf %lf %lf",&n, &dum0, &dum1, &dum2);
               displ[0] = dum0;
               displ[1] = dum1;
               displ[2] = dum2;
            }
            if (n==4) {
               sscanf(mhline, "%d %lf %lf %lf",&n, &dum0, &dum1, &dum2); 
               dir_matrix[0][0] = dum0;
               dir_matrix[0][1] = dum1;
               dir_matrix[0][2] = dum2;
            }
            if (n==5) {
               sscanf(mhline, "%d %lf %lf %lf",&n, &dum0, &dum1, &dum2); 
               dir_matrix[1][0] = dum0;
               dir_matrix[1][1] = dum1;
               dir_matrix[1][2] = dum2;
            }
            if (n==6) {
               sscanf(mhline, "%d %lf %lf %lf",&n, &dum0, &dum1, &dum2); 
               dir_matrix[2][0] = dum0;
               dir_matrix[2][1] = dum1;
               dir_matrix[2][2] = dum2;
            }
            if (n==7) sscanf(mhline, "%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
               &n, &h[0], &h[1], &h[2], &h[3], &h[4], &h[5], &h[6], &h[7], &h[8], &h[9], &h[10], &h[11]);
            if (n==8) sscanf(mhline, "%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
               &n, &r[0], &r[1], &r[2], &r[3], &r[4], &r[5], &r[6], &r[7], &r[8], &r[9], &r[10], &r[11]);
            if (n==9) sscanf(mhline, "%d %lf %lf", &n, wgt, z_cg);
            if (n==10) sscanf(mhline, "%d %d", &n, nmax_i);
            if (n == 99) return(0);

            if ((n > 10) || (n < 0)) {
               printf("Bad input data: (line_identifier outside range 0 to 10). %s\n", mhline);
               return(-1);
            }
      }
q_end: return(2);
}