// tr04.c
// bending Warren truss beam
//
// cc -o tr04 tr04.c -lgd -lpng -lm
// 
// The forces and angles shown are illustrative, not at all accurate.
// A better version of this program would get the forces right.
//
// truss by James Warren 1848

#define  NAME   "tr04"

#define  WAVEL     1.0   //
#define  RX       10.0   //  as drawn
#define  RY       40.0   //
#define  R0        1.0   //

#define  PRELOAD -60.0   //  preload
#define  SY     -40.0    //
#define  SM      -1.0    //
#define  SS       0.0    //  

#define  SOFF      20    //

#define  NSTRESS   101   //  Stress scale range
#define  MSTRESS  300.0  //  Stress color multiplier
#define  NPOINTS    12   //  Number of truss points (overlap sides)

#define  RATE        6   //  Flash frame rate ( 1560x )
#define  NPLOT      90   //  Number of plots

#define  XSIZE    1024   //
#define  YSIZE     320   //
#define  YH       160.0  //

#define  PI        3.141592653589793238
#define  FNT       "DejaVuMonoSans"
#define  FSZ        20

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

//
int         i,  j                        ; // general index

//
FILE        *pngout;                     ; // output file
char        rmmkdir[200]                 ; //
char        png2swf[200]                 ; //
char        filename[80]                 ; //
char        plotlabel[80]                ; // plot label

// plot stuff
int         iplot                        ; // plot number
gdImagePtr  im                           ; // working image
int         black, gray, white           ; // color indexes
int         red,  green, blue            ; // color indexes
int         colorstress[NSTRESS+1]       ; // color index for stress

double       tx[200],  ty[200]           ; // top points
double       bx[200],  by[200]           ; // bottom points
double      stx[200], sty[200]           ; // top stress
double      sbx[200], sby[200]           ; // bottom stress

double      sx=0.0                       ; // stress returned by strut
double      sy=0.0                       ; // stress returned by strut

int         sx0, sy0, sx1, sy1           ; // stress vector points

//-----------------------------------------------------------------
void initialize() {
   sprintf( rmmkdir, "rm -rf %sdir ; mkdir %sdir", NAME, NAME );
   system(  rmmkdir )                    ; // make directory
}

//-----------------------------------------------------------------
void openplot() {
   im      = gdImageCreateTrueColor( (int)XSIZE, (int)YSIZE );
   black   = gdImageColorResolve(im,   0,   0,   0);
   gray    = gdImageColorResolve(im, 128, 128, 128);
   white   = gdImageColorResolve(im, 255, 255, 255);
   red     = gdImageColorResolve(im, 255,   0,   0);
   green   = gdImageColorResolve(im,   0, 255,   0);
   blue    = gdImageColorResolve(im,   0,   0, 255);

   // stress colors
   for( i=0 ; i < NSTRESS ; i++ ) {
      double stress= ((double) i ) / ((double)(NSTRESS-1));
      int    ir = (int) (255.0*     stress ) ;
      int    ig = (int) (255.0*(1.0-stress)) ;
      int    ib = ig                         ;
      colorstress[i] = gdImageColorAllocate (im, ir, ig, ib );
   }
   gdImageSetThickness( im, 5 );               // line thickness of 2
}

//-----------------------------------------------------------------
// draw strut

void strut( double x0, double y0, double x1, double y1, double nom ) {

   double length  = sqrt( (x0-x1)*(x0-x1) + (y0-y1)*(y0-y1) );
   double stress  = (length/nom) - 1.0 ;

   // returned through globals
   sx     = (PRELOAD+stress) * (x0-x1)/length ;  
   sy     = (PRELOAD+stress) * (y0-y1)/length ;

   int    istress = (int)( MSTRESS*stress + 0.5*(NSTRESS-1) );
   int    scolor  ;
   if(      istress <  0       ){ scolor = white               ; }
   else if( istress >= NSTRESS ){ scolor = white               ; }
   else                         { scolor = colorstress[istress]; }

   gdImageLine( im, ((int)x0), ((int)y0), ((int)x1), ((int)y1), scolor );
}

//--------------------------------------------------------------------
void closeplot() {
   sprintf( filename, "%sdir/a%04d.png", NAME, iplot );
   pngout = fopen(  filename, "wb");
   gdImagePngEx(im, pngout, 1);
   fclose(pngout);
   gdImageDestroy(im);
}

//---------------------------------------------------------------------
int main() {


   // -----------------------------------------------------------------

   initialize();

   // define truss points ( two off screen on each side );
   double  deltax = ((double)XSIZE)/((double)(NPOINTS));
   double  deltah = 0.5*deltax ;
   double  deltay = deltax ;
   double  startx = -1.6*deltax ;
   double  starty = YH - 0.5*deltay ;
   double  diag   = sqrt( deltah*deltah + deltay*deltay );  
   double  p                            ;
   double  r0    = R0                   ;
   double  wavel = WAVEL                ;

   for( iplot=0 ;  iplot < NPLOT ; iplot++ ) {
      double ipn = 6*((double)iplot)/((double)NPLOT) ; // 0 to almost 6

      if(      ipn < 1.0 ) {
         r0    = 0.0   ;
         wavel = WAVEL ;
      }
      else if( ipn < 2.0 ) {
         r0    = R0*(ipn-1.0) ;
         wavel = WAVEL ;
      }
      else if( ipn < 3.0 ) {
         r0    = R0    ;
         wavel = WAVEL ;
      }
      else if( ipn < 4.0 ) {
         r0    = R0    ;
         wavel = 0.5*WAVEL*(5.0-ipn) ;
      }
      else if( ipn < 5.0 ) {
         r0    = R0    ;
         wavel = 0.5*WAVEL ;
      }
      else {
         r0 = R0*(6.0-ipn);
         wavel = 0.5 * WAVEL ;
      }
      double  wavem = 2.0*PI/(wavel*XSIZE) ;

      openplot();

      //    t0 ----- t1
      //      \     /  
      //       \   /
      //        \ /
      //         b0 ----- b1

   
      for( i=0 ; i<NPOINTS+5 ; i++ ) {        // strut points loop
         double d0 = deltax*((double)i) ;
         double dh = d0 + deltah        ;

         tx[i]   = startx + d0     ;
         ty[i]   = starty          ;
         p       = wavem*tx[i]     ;
         tx[i]  -= r0 * RX * sin( p )   ;
         ty[i]  += r0 * RY * cos( p )   ;
   
         bx[i]   = startx + dh     ;
         by[i]   = starty + deltay ;
         p       = wavem*bx[i]     ;
         bx[i]  += r0 * RX * sin( p )   ;
         by[i]  += r0 * RY * cos( p )   ;
      }
   
      for( i=0 ; i<NPOINTS+4 ; i++ ) {
         stx[i] = 0.0 ;
         sty[i] = -SY ;
         sbx[i] = 0.0 ;
         sby[i] = +SY ;
      }
   
      // diagonals
      gdImageSetThickness( im,  6 );                 // line thickness
      for( i=0 ; i<NPOINTS+4 ; i++ ) {
         strut( tx[i], ty[i], bx[i],   by[i],  diag ); // t0 to b0
         stx[i]   += SS*sx ;
         sty[i]   += SS*sy ;
         sbx[i]   -= SS*sx ;
         sby[i]   -= SS*sy ;
   
         strut( bx[i], by[i], tx[i+1], ty[i+1], diag ); // b0 to t1
         sbx[i]   += SS*sx ;
         sby[i]   += SS*sy ;
         stx[i+1] -= SS*sx ;
         sty[i+1] -= SS*sy ;
      }
   
      // mains
      gdImageSetThickness( im, 10 );                 // line thickness
      for( i=0 ; i<NPOINTS+4 ; i++ ) {
         strut( tx[i], ty[i], tx[i+1], ty[i+1], deltax ); // t0 to t1
         stx[i]   += SM*sx ;
         sty[i]   -= SM*sy ;
         stx[i+1] -= SM*sx ;
         sty[i+1] += SM*sy ;
   
         strut( bx[i], by[i], bx[i+1], by[i+1], deltax ); // b0 to b1
         sbx[i]   += SM*sx ;
         sby[i]   -= SM*sy ;
         sbx[i+1] -= SM*sx ;
         sby[i+1] += SM*sy ;
      }
   
      // stress arrows
      gdImageSetThickness( im,  4 );                // line thickness
   
      for( i=1 ; i<NPOINTS+3 ; i++ ) {
   
         sx0 = ((int)tx[i]) ;
         sy0 = ((int)ty[i]) - SOFF  ;
         // sx1 = sx0 - (int)stx[i] ;
         sx1 = sx0 ;
         sy1 = sy0 - (int)sty[i] ;
          
         gdImageLine( im, sx0, sy0, sx1, sy1, white );
         
         sx0 = ((int)bx[i]) ;
         sy0 = ((int)by[i]) + SOFF ;
         // sx1 = sx0 - (int)sbx[i] ;
         sx1 = sx0 ;
         sy1 = sy0 - (int)sby[i] ;
          
         gdImageLine( im, sx0, sy0, sx1, sy1, white );
      }  
      closeplot();
   }

   // make flash video
   sprintf( png2swf, "png2swf -o %s.swf -r%3d %sdir/*.png", NAME, RATE, NAME );
   system( png2swf );

   return(0);
}