/* dual.c - sdt - 20-OCT-1985 - create a transformation between the
 *                              purple octohedron and the blue cube,
 *                              to illustrate that they are dual
 *                              Platonic solids.
 *
 * usage:  dual [filename]
 * if no file is given, the formatted data will go to stdout.
 * if the named file exists, program will abort.
 * example:  dual x.x
 *           crystal -r x.x
 *
 * mods 03-NOV-1985 to take advantage of the auto reverse mode (-r)
 * added to crystal (the program that displays 3d movies).
 */

#include "stdio.h"

   static int nlines = 9 ;
/*    primitives based on octohedron; (a,b,c) -> (x,y,z) later.
 */
   static float a[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 0.5, 0.0 } ;
   static float b[] = { 0.0, 0.0, 0.5, 1.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0 } ;
   static float c[] = { 1.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0 } ;
   static int upo[] = {   2,   1,   1,   2,   1,   1,   2,   1,   1,   0 } ;

/*    primitives based on cube; (d,e,f) -> (x,y,z) later.
 */
   static float d[] = { 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 } ;
   static float e[] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0 } ;
   static float f[] = { 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0 } ;
   static int upc[] = {   2,   1,   1,   2,   1,   1,   2,   1,   1,   0 } ;

/*    the format string to use when writing x,y,z,up.
 */
   static char *fmt = { "%.3f %.3f %.3f %1d\n" } ;

/*    scaled to inscribe the unit sphere
 */
   float x[10], y[10], z[10] ;
   double pi ;

main( argc, argv )
int argc;
char *argv[] ;
{
   FILE *fp ;
   extern FILE *fopen() ;
   extern int fclose() ;
   extern double atan() ;

   pi = atan( 1.0 ) * 4.0 ;

   if( argc == 1 ) {  /* no args */
      dual( stdout ) ;
   }
   else {
      if( (fp = fopen( *++argv, "r" )) != NULL )
         abort( "dual: file %s exists - will not overwrite!", *argv ) ;
      if( (fp = fopen( *argv, "w" )) == NULL )
         abort( "dual: can't open %s for output", *argv ) ;
      dual( fp ) ;
      fclose( fp ) ;
   }
}

dual( fp )
FILE *fp ;
{
   int i, j ; /* loop index */
   int jmax = 40 ;
   extern int scale() ;
   extern double cos() ;
   float sigma ;

   for( j = 0 ; j < jmax ; j++ ) {
      fprintf( stderr, "j=%d\n", j ) ;
      sigma = (float)j / (float)jmax ;
      sigma = (1.0 - cos( sigma * pi )) * 0.5 ;
      octo( sigma );
      scale() ;
      spit( fp, upo ) ;
   }

   for( j = 0 ; j < jmax ; j++ ) {
      fprintf( stderr, "j=%d\n", j ) ;
      sigma = 1.0 - (float)j / (float)jmax ;
      sigma = (1.0 - cos( sigma * pi )) * 0.5 ;
      cube( sigma ) ;
      scale() ;
      spit( fp, upc ) ;
   }
}

int scale() /* return 0 if it worked, else 1 */
{
   int i ;
   extern double sqrt() ;
   double sqrad ;  /* square of radius */
   double maxsq ;  /* max of the squares */
   double maxradi ; /* one over maximum squared radius */

   for( maxsq = 0.0, i = 0 ; i <= nlines ; i++ )
      if( maxsq < (sqrad = x[i] * x[i] + y[i] * y[i] + z[i] * z[i]) )
         maxsq = sqrad ;
   if( maxsq > 0.0 ) {
      maxradi = 1.0 / sqrt( maxsq ) ;
      for( i = 0 ; i <= nlines ; i++ ) {
         x[i] *= maxradi ;
         y[i] *= maxradi ;
         z[i] *= maxradi ;
      }
      return( 0 ) ;
   }
   else
      return( 1 ) ;
}

spit( fp, up )
FILE *fp ;
int *up ;
{
   int i ;

   fprintf( fp, "%d\n", (nlines + 1) * 8 - 1) ;
   for( i = 0 ; i <= nlines ; i++ )
      fprintf( fp, fmt,  x[i],  y[i],  z[i], up[i] ) ;
   for( i = 0 ; i <= nlines ; i++ )
      fprintf( fp, fmt, -y[i],  x[i],  z[i], up[i] ) ;
   for( i = 0 ; i <= nlines ; i++ )
      fprintf( fp, fmt, -x[i], -y[i],  z[i], up[i] ) ;
   for( i = 0 ; i <= nlines ; i++ )
      fprintf( fp, fmt,  y[i], -x[i],  z[i], up[i] ) ;
   for( i = 0 ; i <= nlines ; i++ )
      fprintf( fp, fmt,  x[i],  y[i], -z[i], up[i] ) ;
   for( i = 0 ; i <= nlines ; i++ )
      fprintf( fp, fmt, -y[i],  x[i], -z[i], up[i] ) ;
   for( i = 0 ; i <= nlines ; i++ )
      fprintf( fp, fmt, -x[i], -y[i], -z[i], up[i] ) ;
   for( i = 0 ; i <= nlines ; i++ )
      fprintf( fp, fmt,  y[i], -x[i], -z[i], up[i] ) ;
}

octo( sigma )
float sigma ; /* 0.0 to 1.0 */
{
   int i ;
   float s, t ;

   for( i = 0 ; i <= nlines ; i++ ) {
      x[i] = a[i] ;
      y[i] = b[i] ;
      z[i] = c[i] ;
   }
   s = sigma * 0.5 ;       /* 0.0 to 0.5 */
   t = 1.0 - sigma * 0.5 ; /* 1.0 to 0.5 */

   x[0] = x[9] = s ;
   z[0] = z[9] = t ;

   y[1] = s ;
   z[1] = t ;

   y[3] = t ;
   z[3] = s ;

   x[4] = s ;
   y[4] = t ;

   x[6] = t ;
   y[6] = s ;

   x[7] = t ;
   z[7] = s ;
}

cube( sigma )
float sigma ; /* 0.0 to 1.0 */
{
   int i ;
   float s ;

   for( i = 0 ; i <= nlines ; i++ ) {
      x[i] = d[i] ;
      y[i] = e[i] ;
      z[i] = f[i] ;
   }
   s = 1.0 - sigma ; /* 1.0 to 0.0 */

   y[0] = y[7] = y[9] = s ;

   x[1] = x[3] = s ;

   z[4] = z[6] = s ;
}