/* Program to generate wireframe model of 3D objects
   Niels Lund, September 1998
*/

/* 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 written to disk (balle.pnt)
*/


/* Special version for the Ballerina Small Satellite
   Niels Lund, October 1998
*/


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

int calpoint(int, int, int);
int calrect(double, double, int);
double inerti(int, int, int);

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

void main(int argc, char* argv)
{
  double al, dx, dy, mkg, msum, msum_x, msum_y, msum_z, z_cg, z_mount, r_mount;
  double xmnt, ymnt, zmnt, rmnt, wfac, Ixx, Iyy, Izz;
  int i, j, k, l, m, n, w, nr, nmax, imax;
  char s[256];
  char *ss;

  al = 109.0 * pi / 180.0 ;
  /* al is the characteristic angle between
     the normals to the tetrahedron faces */

  pointnum = 1;

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

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


  nmax = 15;
  if (argc>1) nmax = argv[1];
  if (nmax > 90) nmax = 90;

  out1 = fopen("balle.pnt","wt") ;	/* Output file */
  if ( out1 == NULL )
  {
	printf("File opening error!!!\n") ;
	exit(0) ;
  }

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

  fprintf(out1, "# Wireframe model of the Ballerina satellite. Niels Lund, dec 98");

  imax = 8;	/* generator line for a WATCH detector */
  h[0] = 0.0L;
  h[1] = 11.6L;
  h[2] = 14.6L;
  h[3] = 19.0L;
  h[4] = 19.0L;
  h[5] = 22.5L;
  h[6] = 22.5L;
  h[7] = 28.0L;
  ra[0] = 4.0L;
  ra[1] = 4.0L;
  ra[2] = 7.0L;
  ra[3] = 7.0L;
  ra[4] = 8.5L;
  ra[5] = 8.5L;
  ra[6] = 12.0L;
  ra[7] = 13.0L;

/*  WATCH 0 */

  mkg = weights[elem] = 4.0;		/* mass of unit (kg) */
  z_cg = 20.5;		/* CG distance from reference point */
  z_mount = 19.0;       /* Mounting surface distance from ref. pnt. */
  r_mount = 8.0;	/* Radius of mounting circle */

  displ[0] = 4.0L;
  displ[1] = 0.0L;
  displ[2] = 50.0L;

	BA_TO_WA[0][0] = cos(pi/2.0 - 0.5*al);
	BA_TO_WA[1][0] = 0.0;
	BA_TO_WA[2][0] = -sin(pi/2.0 - 0.5*al);

	BA_TO_WA[0][1] = 0.0L;
	BA_TO_WA[1][1] = 1.0L;
	BA_TO_WA[2][1] = 0.0L;

	BA_TO_WA[0][2] = sin(pi/2.0 - 0.5*al);
	BA_TO_WA[1][2] = 0.0;
	BA_TO_WA[2][2] = cos(pi/2.0 - 0.5*al);

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = 10.0;
  geo_size[elem][1] = 10.0;
  geo_size[elem][2] = 10.0;
  elem++;

  pointnum = calpoint(nmax, imax, pointnum);

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

  linenum += nmax * imax;

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

/*  WATCH 1 */

  weights[elem] = 4.0;		/* mass of unit (kg) */
  displ[0] = 7.0L;
  displ[1] = 0.0L;
  displ[2] = 33.0L;

	BA_TO_WA[0][0] = cos(pi/2.0 + 0.5*al) ;   /* WATCH 1 x-axis in BALLERINA */
	BA_TO_WA[1][0] = 0.0 ;
	BA_TO_WA[2][0] = -sin(pi/2.0 + 0.5*al) ;

	BA_TO_WA[0][1] = 0.0 ;   /* WATCH 1 y-axis in BALLERINA */
	BA_TO_WA[1][1] = 1.0 ;
	BA_TO_WA[2][1] = 0.0 ;

	BA_TO_WA[0][2] = sin(pi/2.0 + 0.5*al) ;   /* WATCH 1 z-axis in BALLERINA */
	BA_TO_WA[1][2] = 0.0 ;
	BA_TO_WA[2][2] = cos(pi/2.0 + 0.5*al) ;

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = 10.0;
  geo_size[elem][1] = 10.0;
  geo_size[elem][2] = 10.0;
  elem++;

  pointnum = calpoint(nmax, imax, pointnum);

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

  linenum += nmax * imax;

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

/*  WATCH 2 */

  weights[elem] = 4.0;		/* mass of unit (kg) */
  displ[0] = -2.0L;
  displ[1] = -2.5L;
  displ[2] = 19.0L;

     /* WATCH 2 x-axis in BALLERINA */
	BA_TO_WA[0][0] = cos(-pi/2.0) * cos(-0.5*al);
	BA_TO_WA[1][0] = -cos(-pi/2.0) * sin(-0.5*al);
	BA_TO_WA[2][0] = sin(-pi/2.0);

     /* WATCH 2 y-axis in BALLERINA */
	BA_TO_WA[0][1] = sin(-0.5*al);
	BA_TO_WA[1][1] = cos(-0.5*al);
	BA_TO_WA[2][1] = 0.0;

     /* WATCH 2 z-axis in BALLERINA */
	BA_TO_WA[0][2] = sin(-pi/2.0) * cos(-0.5*al);
	BA_TO_WA[1][2] = -sin(-pi/2.0) * sin(-0.5*al);
	BA_TO_WA[2][2] = cos(-pi/2.0);

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = 10.0;
  geo_size[elem][1] = 10.0;
  geo_size[elem][2] = 10.0;
  elem++;

  pointnum = calpoint(nmax, imax, pointnum);

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

  linenum += nmax * imax;

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

/*  WATCH 3 */

  weights[elem] = 4.0;		/* mass of unit (kg) */
  displ[0] = -2.0L;
  displ[1] = 2.5L;
  displ[2] = 19.0L;

     /* WATCH 3 x-axis in BALLERINA */
	BA_TO_WA[0][0] = -cos(-pi/2.0) * cos(-0.5*al);
	BA_TO_WA[1][0] = -cos(-pi/2.0) * sin(-0.5*al);
	BA_TO_WA[2][0] = sin(-pi/2.0);

     /* WATCH 3 y-axis in BALLERINA */
	BA_TO_WA[0][1] = sin(-0.5*al);
	BA_TO_WA[1][1] = -cos(-0.5*al);
	BA_TO_WA[2][1] = 0.0;

     /* WATCH 3 z-axis in BALLERINA */
	BA_TO_WA[0][2] = sin(-pi/2.0) * cos(-0.5*al);
	BA_TO_WA[1][2] = sin(-pi/2.0) * sin(-0.5*al);
	BA_TO_WA[2][2] = cos(-pi/2.0);

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = 10.0;
  geo_size[elem][1] = 10.0;
  geo_size[elem][2] = 10.0;
  elem++;

  pointnum = calpoint(nmax, imax, pointnum);

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

  linenum += nmax * imax;

  xmnt = displ[0] + BA_TO_WA[0][2] * z_mount;
  ymnt = displ[1] + BA_TO_WA[1][2] * z_mount;
  zmnt = displ[2] + BA_TO_WA[2][2] * z_mount;
  printf("\n WATCH 3.   Mount center (x,y,z cm): %4.1lf, %4.1lf, %4.1lf;  Radius: %3.1lf",
	xmnt, ymnt, zmnt, r_mount);
  printf("\n WATCH 3.   Mass %4.1lf kg, CG-Mount sep: %4.1lf cm,  Vector: %4.3lf, %4.3lf, %4.3lf\n",
	mkg, z_cg-z_mount, BA_TO_WA[0][2], BA_TO_WA[1][2], BA_TO_WA[2][2]);

/*  Separation Ring */

  mkg = weights[elem] = 2.0;        /* mass of unit (kg) */
  z_cg = 1.5;		  /* CG distance from reference point */
  z_mount = 0.0;       /* Mounting surface distance from ref. pnt. */
  r_mount = 11.5;	/* Radius of mounting circle */

  imax = 2;
  nmax = 15;

  h[0] = 0.0L;
  h[1] = 3.0L;
  ra[0] = 11.5;
  ra[1] = 11.5;

  displ[0] = 0.0L;
  displ[1] = 0.0L;
  displ[2] = 0.0L;

	BA_TO_WA[0][0] = 1.0 ;   /* Sep. Ring x-axis in BALLERINA */
	BA_TO_WA[1][0] = 0.0 ;
	BA_TO_WA[2][0] = 0.0 ;

	BA_TO_WA[0][1] = 0.0 ;     /* y-axis in BALLERINA */
	BA_TO_WA[1][1] = 1.0 ;
	BA_TO_WA[2][1] = 0.0 ;

	BA_TO_WA[0][2] = 0.0 ;     /* 0 z-axis in BALLERINA */
	BA_TO_WA[1][2] = 0.0 ;
	BA_TO_WA[2][2] = 1.0 ;

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = 11.5;
  geo_size[elem][1] = 11.5;
  geo_size[elem][2] = 1.5;
  elem++;

  pointnum = calpoint(nmax, imax, pointnum);

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

  linenum += nmax * imax;

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

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

/* Telescope */

  mkg = weights[elem] = 7.0;		/* mass of unit (kg) */
  z_cg = 64.0;		/* CG distance from reference point */
  z_mount = 65.0;       /* Mounting surface distance from ref. pnt. */
  r_mount = 16.0;	/* Radius of mounting circle */

  imax = 3;
  h[0] = 0.0L;		/* Focal plane end of telescope */
  h[1] = 20.0L;		/* Intermediate level */
  h[2] = 75.0L;		/* Telescope aperture end */
  ra[0] = 5.0L;		/* radius at focal plane end */
  ra[1] = 5.0L;		/* radius at intermediate level */
  ra[2] = 16.0L;	/* radius at aperture end */

  displ[0] = -13.0L;
  displ[1] = -0.0L;
  displ[2] = 5.0L;

	BA_TO_WA[0][0] = 1.0L;
	BA_TO_WA[1][0] = 0.0L;
	BA_TO_WA[2][0] = 0.0L;

	BA_TO_WA[0][1] = 0.0L;
	BA_TO_WA[1][1] = 1.0L;
	BA_TO_WA[2][1] = 0.0L;

	BA_TO_WA[0][2] = 0.0L;
	BA_TO_WA[1][2] = 0.0L;
	BA_TO_WA[2][2] = 1.0L;

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = 12.5;
  geo_size[elem][1] = 12.5;
  geo_size[elem][2] = 12.5;
  elem++;

  pointnum = calpoint(nmax, imax, pointnum);

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

  linenum += nmax * imax;

  xmnt = displ[0] + BA_TO_WA[0][2] * z_mount;
  ymnt = displ[1] + BA_TO_WA[1][2] * z_mount;
  zmnt = displ[2] + BA_TO_WA[2][2] * z_mount;
  printf("\n Telescope  Mount center (x,y,z cm): %4.1lf, %4.1lf, %4.1lf;  Radius: %3.1lf",
	xmnt, ymnt, zmnt, r_mount);
  printf("\n Telescope  Mass %4.1lf kg, CG-Mount sep: %4.1lf cm,  Vector: %4.3lf, %4.3lf, %4.3lf\n",
	mkg, z_cg-z_mount, BA_TO_WA[0][2], BA_TO_WA[1][2], BA_TO_WA[2][2]);

/*  Focal Plane Electronics */

  mkg = weights[elem] = 4.0;		/* mass of unit (kg) */
  z_cg = 5.0;		/* CG distance from reference point */
  z_mount = 0.0;       /* Mounting surface distance from ref. pnt. */

  imax = 2;
  nmax = 4;

  h[0] = 0.0L;
  h[1] = 10.0L;
  dx = 9.0L;
  dy = 8.0L;
  ra[0] = sqrt(dx*dx + dy*dy);
  ra[1] = sqrt(dx*dx + dy*dy);

  displ[0] = -21.0L;
  displ[1] = 0.0L;
  displ[2] = 5.0L;

	BA_TO_WA[0][0] = 1.0L;
	BA_TO_WA[1][0] = 0.0L;
	BA_TO_WA[2][0] = 0.0L;

	BA_TO_WA[0][1] = 0.0L;
	BA_TO_WA[1][1] = 1.0L;
	BA_TO_WA[2][1] = 0.0L;

	BA_TO_WA[0][2] = 0.0L;
	BA_TO_WA[1][2] = 0.0L;
	BA_TO_WA[2][2] = 1.0L;

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = dx;
  geo_size[elem][1] = dy;
  geo_size[elem][2] = 5.0;
  elem++;

  pointnum = calrect(dx, dy, pointnum);

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

  linenum += nmax * imax;

  xmnt = displ[0] + BA_TO_WA[0][2] * z_mount;
  ymnt = displ[1] + BA_TO_WA[1][2] * z_mount;
  zmnt = displ[2] + BA_TO_WA[2][2] * z_mount;
  printf("\n Foc.Pl.El. Mount center (x,y,z cm): %4.1lf, %4.1lf, %4.1lf;  Rect: %3.0lfx%3.0lf",
	xmnt, ymnt, zmnt, dx+dx, dy+dy);
  printf("\n Foc.Pl.El. Mass %4.1lf kg, CG-Mount sep: %4.1lf cm,  Vector: %4.3lf, %4.3lf, %4.3lf\n",
	mkg, z_cg-z_mount, BA_TO_WA[0][2], BA_TO_WA[1][2], BA_TO_WA[2][2]);

/*  Telescope support bracket */

  mkg = 0.2;		/* mass of unit (kg) */
  z_cg = 0.5;		/* CG distance from reference point */
  z_mount = 0.0;       /* Mounting surface distance from ref. pnt. */

  imax = 2;
  nmax = 4;

  h[0] = 0.0L;
  h[1] = 1.0L;
  dx = 2.0L;
  dy = 2.0L;
  ra[0] = sqrt(dx*dx + dy*dy);
  ra[1] = sqrt(dx*dx + dy*dy);

  displ[0] = -13.0L;
  displ[1] = 15.0L;
  displ[2] = 65.0L;

	BA_TO_WA[0][0] = cos(pi/2.0);	/* Bracket x-axis in BALLERINA */
	BA_TO_WA[1][0] = 0.0;
	BA_TO_WA[2][0] = -sin(pi/2.0);

	BA_TO_WA[0][1] = 0.0L;
	BA_TO_WA[1][1] = 1.0L;
	BA_TO_WA[2][1] = 0.0L;

	BA_TO_WA[0][2] = sin(pi/2.0);
	BA_TO_WA[1][2] = 0.0;
	BA_TO_WA[2][2] = cos(pi/2.0);

  pointnum = calrect(dx, dy, pointnum);

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

  linenum += nmax * imax;

  xmnt = displ[0] + BA_TO_WA[0][2] * z_mount;
  ymnt = displ[1] + BA_TO_WA[1][2] * z_mount;
  zmnt = displ[2] + BA_TO_WA[2][2] * z_mount;
  printf("\n Tel. Mount1   Mount center (x,y,z cm): %4.1lf, %4.1lf, %4.1lf;  Rect: %3.0lfx%3.0lf",
	xmnt, ymnt, zmnt, dx+dx, dy+dy);
  printf("\n Tel. Mount1  Mass %4.1lf kg", mkg);


/*  Telescope support bracket */

  mkg = 0.2;		/* mass of unit (kg) */
  z_cg = 0.5;		/* CG distance from reference point */
  z_mount = 0.0;       /* Mounting surface distance from ref. pnt. */

  imax = 2;
  nmax = 4;

  h[0] = 0.0L;
  h[1] = 1.0L;
  dx = 2.0L;
  dy = 2.0L;
  ra[0] = sqrt(dx*dx + dy*dy);
  ra[1] = sqrt(dx*dx + dy*dy);

  displ[0] = -13.0L;
  displ[1] = -15.0L;
  displ[2] = 65.0L;

  pointnum = calrect(dx, dy, pointnum);

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

  linenum += nmax * imax;

  xmnt = displ[0] + BA_TO_WA[0][2] * z_mount;
  ymnt = displ[1] + BA_TO_WA[1][2] * z_mount;
  zmnt = displ[2] + BA_TO_WA[2][2] * z_mount;
  printf("\n Tel. Mount2   Mount center (x,y,z cm): %4.1lf, %4.1lf, %4.1lf;  Rect: %3.0lfx%3.0lf",
	xmnt, ymnt, zmnt, dx+dx, dy+dy);
  printf("\n Tel. Mount2  Mass %4.1lf kg", mkg);

/*  Focal Plane Electronics */

  mkg = weights[elem] = 4.0;		/* mass of unit (kg) */
  z_cg = 5.0;		/* CG distance from reference point */
  z_mount = 0.0;       /* Mounting surface distance from ref. pnt. */

  imax = 2;
  nmax = 4;

  h[0] = 0.0L;
  h[1] = 10.0L;
  dx = 9.0L;
  dy = 8.0L;
  ra[0] = sqrt(dx*dx + dy*dy);
  ra[1] = sqrt(dx*dx + dy*dy);

  displ[0] = -21.0L;
  displ[1] = 0.0L;
  displ[2] = 5.0L;

	BA_TO_WA[0][0] = 1.0L;
	BA_TO_WA[1][0] = 0.0L;
	BA_TO_WA[2][0] = 0.0L;

	BA_TO_WA[0][1] = 0.0L;
	BA_TO_WA[1][1] = 1.0L;
	BA_TO_WA[2][1] = 0.0L;

	BA_TO_WA[0][2] = 0.0L;
	BA_TO_WA[1][2] = 0.0L;
	BA_TO_WA[2][2] = 1.0L;

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = dx;
  geo_size[elem][1] = dy;
  geo_size[elem][2] = 5.0;
  elem++;

  pointnum = calrect(dx, dy, pointnum);

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

  linenum += nmax * imax;

  xmnt = displ[0] + BA_TO_WA[0][2] * z_mount;
  ymnt = displ[1] + BA_TO_WA[1][2] * z_mount;
  zmnt = displ[2] + BA_TO_WA[2][2] * z_mount;
  printf("\n Foc.Pl.El. Mount center (x,y,z cm): %4.1lf, %4.1lf, %4.1lf;  Rect: %3.0lfx%3.0lf",
	xmnt, ymnt, zmnt, dx+dx, dy+dy);
  printf("\n Foc.Pl.El. Mass %4.1lf kg, CG-Mount sep: %4.1lf cm,  Vector: %4.3lf, %4.3lf, %4.3lf\n",
	mkg, z_cg-z_mount, BA_TO_WA[0][2], BA_TO_WA[1][2], BA_TO_WA[2][2]);

/*  Satellite Body */

  mkg = 18.8;		/* mass of unit (kg) */
  z_cg = 30.0;		/* CG distance from reference point */
  z_mount = 0.0;       /* Mounting surface distance from ref. pnt. */

  imax = 2;
  nmax = 4;

  h[0] = 0.0L;
  h[1] = 77.0L;
  dx = 30.0L;
  dy = 30.0L;
  ra[0] = sqrt(dx*dx + dy*dy);
  ra[1] = sqrt(dx*dx + dy*dy);

  displ[0] = 0.0L;
  displ[1] = 0.0L;
  displ[2] = 3.0L;

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = dx;
  geo_size[elem][1] = dy;
  geo_size[elem][2] = 35.0;
  elem++;

  pointnum = calrect(dx, dy, pointnum);

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

  linenum += nmax * imax;

  xmnt = displ[0] + BA_TO_WA[0][2] * z_mount;
  ymnt = displ[1] + BA_TO_WA[1][2] * z_mount;
  zmnt = displ[2] + BA_TO_WA[2][2] * z_mount;
  printf("\n Sat. Body  Mount center (x,y,z cm): %4.1lf, %4.1lf, %4.1lf;  Rect: %3.0lfx%3.0lf",
	xmnt, ymnt, zmnt, dx+dx, dy+dy);
  printf("\n Sat. Body  Mass %4.1lf kg, CG-Mount sep: %4.1lf cm,  Vector: %4.3lf, %4.3lf, %4.3lf\n",
	mkg, z_cg-z_mount, BA_TO_WA[0][2], BA_TO_WA[1][2], BA_TO_WA[2][2]);

/*  Battery */

  mkg = weights[elem] = 8.0;		/* mass of unit (kg) */
  z_cg = 4.3;		/* CG distance from reference point */
  z_mount = 0.0;       /* Mounting surface distance from ref. pnt. */

  imax = 2;
  nmax = 4;

  h[0] = 0.0L;
  h[1] = 8.5L;
  dx = 5.5L;
  dy = 12.5L;
  ra[0] = sqrt(dx*dx + dy*dy);
  ra[1] = sqrt(dx*dx + dy*dy);

  displ[0] = 4.0L;
  displ[1] = 0.0L;
  displ[2] = 5.0L;

	BA_TO_WA[0][0] = 1.0L;
	BA_TO_WA[1][0] = 0.0L;
	BA_TO_WA[2][0] = 0.0L;

	BA_TO_WA[0][1] = 0.0L;
	BA_TO_WA[1][1] = 1.0L;
	BA_TO_WA[2][1] = 0.0L;

	BA_TO_WA[0][2] = 0.0L;
	BA_TO_WA[1][2] = 0.0L;
	BA_TO_WA[2][2] = 1.0L;

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = dx;
  geo_size[elem][1] = dy;
  geo_size[elem][2] = 4.25;
  elem++;

  pointnum = calrect(dx, dy, pointnum);

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

  linenum += nmax * imax;

  xmnt = displ[0] + BA_TO_WA[0][2] * z_mount;
  ymnt = displ[1] + BA_TO_WA[1][2] * z_mount;
  zmnt = displ[2] + BA_TO_WA[2][2] * z_mount;
  printf("\n Battery    Mount center (x,y,z cm): %4.1lf, %4.1lf, %4.1lf;  Rect: %3.0lfx%3.0lf",
	xmnt, ymnt, zmnt, dx+dx, dy+dy);
  printf("\n Battery    Mass %4.1lf kg, CG-Mount sep: %4.1lf cm,  Vector: %4.3lf, %4.3lf, %4.3lf\n",
	mkg, z_cg-z_mount, BA_TO_WA[0][2], BA_TO_WA[1][2], BA_TO_WA[2][2]);

/*  Electronics Box 1 */

  mkg = weights[elem] = 11.5;		/* mass of unit (kg) */
  z_cg = 15.0;		/* CG distance from reference point */
  z_mount = 0.0;       /* Mounting surface distance from ref. pnt. */

  imax = 2;
  nmax = 4;

  h[0] = 0.0L;
  h[1] = 30.0L;
  dx = 9.5L;
  dy = 12.0L;
  ra[0] = sqrt(dx*dx + dy*dy);
  ra[1] = sqrt(dx*dx + dy*dy);

  displ[0] = 1.0L;
  displ[1] = -17.0L;
  displ[2] = 23.0L;

	BA_TO_WA[0][0] = cos(0.25);	/* Elbox 1 x-axis in BALLERINA */
	BA_TO_WA[1][0] = 0.0;
	BA_TO_WA[2][0] = -sin(0.25);

	BA_TO_WA[0][1] = 0.0L;
	BA_TO_WA[1][1] = 1.0L;
	BA_TO_WA[2][1] = 0.0L;

	BA_TO_WA[0][2] = sin(0.25);
	BA_TO_WA[1][2] = 0.0;
	BA_TO_WA[2][2] = cos(0.25);

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = dx;
  geo_size[elem][1] = dy;
  geo_size[elem][2] = 15.0;
  elem++;

  pointnum = calrect(dx, dy, pointnum);

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

  linenum += nmax * imax;

  xmnt = displ[0] + BA_TO_WA[0][2] * z_mount;
  ymnt = displ[1] + BA_TO_WA[1][2] * z_mount;
  zmnt = displ[2] + BA_TO_WA[2][2] * z_mount;
  printf("\n El. Box 1  Mount center (x,y,z cm): %4.1lf, %4.1lf, %4.1lf;  Rect: %3.0lfx%3.0lf",
	xmnt, ymnt, zmnt, dx+dx, dy+dy);
  printf("\n El. Box 1  Mass %4.1lf kg, CG-Mount sep: %4.1lf cm,  Vector: %4.3lf, %4.3lf, %4.3lf\n",
	mkg, z_cg-z_mount, BA_TO_WA[0][2], BA_TO_WA[1][2], BA_TO_WA[2][2]);

/*  Electronics Box 2 */

  mkg = weights[elem] = 11.5;		/* mass of unit (kg) */
  z_cg = 15.0;		/* CG distance from reference point */
  z_mount = 0.0;       /* Mounting surface distance from ref. pnt. */

  imax = 2;
  nmax = 4;

  h[0] = 0.0L;
  h[1] = 30.0L;
  dx = 9.5L;
  dy = 12.0L;
  ra[0] = sqrt(dx*dx + dy*dy);
  ra[1] = sqrt(dx*dx + dy*dy);

  displ[0] = 1.0L;
  displ[1] = 17.0L;
  displ[2] = 23.0L;

  geo_pos[elem][0] = displ[0] + z_cg * BA_TO_WA[0][2];
  geo_pos[elem][1] = displ[1] + z_cg * BA_TO_WA[1][2];
  geo_pos[elem][2] = displ[2] + z_cg * BA_TO_WA[2][2];
  geo_size[elem][0] = dx;
  geo_size[elem][1] = dy;
  geo_size[elem][2] = 15.0;
  elem++;

  pointnum = calrect(dx, dy, pointnum);

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

  linenum += nmax * imax;

  xmnt = displ[0] + BA_TO_WA[0][2] * z_mount;
  ymnt = displ[1] + BA_TO_WA[1][2] * z_mount;
  zmnt = displ[2] + BA_TO_WA[2][2] * z_mount;
  printf("\n El. Box 2  Mount center (x,y,z cm): %4.1lf, %4.1lf, %4.1lf;  Rect: %3.0lfx%3.0lf",
	xmnt, ymnt, zmnt, dx+dx, dy+dy);
  printf("\n El. Box 2  Mass %4.1lf kg, CG-Mount sep: %4.1lf cm,  Vector: %4.3lf, %4.3lf, %4.3lf\n",
	mkg, z_cg-z_mount, BA_TO_WA[0][2], BA_TO_WA[1][2], BA_TO_WA[2][2]);

  /* *************************** */
  /* Summary print */

  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",
	     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();
}

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

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] += BA_TO_WA[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, delta, 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] += BA_TO_WA[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);
}