#!/usr/local/bin/perl
#
# shell.pl : draw a seashell
# (c) 2013 jl_morel@bribes.org - http://http://bribes.org/perl
use strict;
use warnings;
use OpenGL qw/ :all /;
use Math::Trig;

#------ Base vectors

my @i = ( 1, 0, 0 );
my @j = ( 0, 1, 0 );
my @k = ( 0, 0, 1 );
my @O = ( 0, 0, 0 );

# ------ Parameters of the seashell

my $tmax = 16 * pi;    # 8 rounds

my $r0 = 0.06;         # initial distance to the axis.
my $kr = 0.05;         # $r = $r0 * exp( $kr * $t )

my $z0 = 0.045;        # initial step
my $kz = 0.075;         # $z = $z0 * exp( $kz * $t )

my @Ap;                # shape of the aperture (a circle here )

my $n  = 32;
my $ra = 0.02;         # initial radius of the aperture
my $ka = 0.08;         # aperture radius = $ra * exp( $ka * $t )

for ( 0 .. 2 * $n - 1 )
{                      # here the aperture is a circle (=polygon with 2*$n sides)
  push @Ap, [ $ra * cos( $_ * pi / $n ), $ra * sin( $_ * pi / $n ), 0 ];
}

my @apex = ( 0, 2, 0 );    # shell apex

##--- Parametric equation of the helico-spiral

sub R {
  my $t = shift;
  return $apex[0] + $r0 * exp( $kr * $t ) * cos($t),
         $apex[1] - $z0 * exp( $kz * $t ),
         $apex[2] + $r0 * exp( $kr * $t ) * sin($t);
}

##--- First derivative

sub Rp {
  my $t = shift;
  return $r0 * exp( $kr * $t ) * ( $kr * cos($t) - sin($t) ),
         -$z0 * $kz * exp( $kz * $t ),
         $r0 * exp( $kr * $t ) * ( cos($t) + $kr * sin($t) );
}

##--- Second derivative

sub Rs {
  my $t = shift;
  return
    $r0 * exp( $kr * $t ) * ( ($kr * $kr - 1) * cos($t) - 2 * $kr * sin($t) ),
    -$z0 * $kz * $kz * exp( $kz * $t ),
    $r0 * exp( $kr * $t ) * ( ($kr * $kr - 1) * sin($t) + 2 * $kr * cos($t)  );
}

my $shell;    # OpenGL display list for the shell

#------ Draw the shell

sub DrawShell {

  my $shell = glGenLists(1);
  glNewList( $shell, GL_COMPILE );
  my @C0;     # previous aperture circle

  for ( my $t = 0 ; $t <= $tmax ; $t += 0.1 ) {
    my @R = R($t);
    my @T = NormalizeVector( Rp($t) );                        # tangent
    my @B = NormalizeVector( CrossProduct( @T, Rs($t) ) );    # binormal
    my @N = CrossProduct( @B, @T );                           # normal

    my @C;    # the current aperture circle
              # in the ( @R; @T, @N, @B ) coordinate system
    foreach (@Ap) {
      push @C,
        [
        AddVectors(
          AddVectors(
            ScaleVector( @N, exp( $ka * $t ) * $_->[0] ),
            ScaleVector( @B, exp( $ka * $t ) * $_->[1] )
          ),
          @R
        )
        ];
    }

    glBegin(GL_TRIANGLE_STRIP);  # draws the band between the circles @C0 and @C
    if ($t) {
      for ( my $i = -1 ; $i < 2 * $n - 1 ; $i++ ) {
        glNormal3d(
          FindUnitNormal( @{ $C0[$i] }, @{ $C0[ $i + 1 ] }, @{ $C[$i] } ) );
        glVertex3f( @{ $C0[$i] } );
        glVertex3f( @{ $C[$i] } );
        glVertex3f( @{ $C0[ $i + 1 ] } );
        glVertex3f( @{ $C[ $i + 1 ] } );
      }
    }
    glEnd();
    @C0 = @C;                    # for the next step
  }
  glEndList;
  return $shell;
}

#------ Draw the scene

my $spin = 0;                    # rotation angle

sub display {
  glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );
  glLoadIdentity();
  gluLookAt( 2, 4, 10, 0, 0, 0, 0, 1, 0 );
  glPushMatrix();
    glScalef( 1.9, 1.9, 1.9 );
    glRotatef( $spin, -sqrt(2) / 2, sqrt(2) / 2, 0 );
    glCallList($shell);            # draw the shell
  glPopMatrix();
  glutSwapBuffers();

  # debug code
#  if ( ( my $e = glGetError() ) != GL_NO_ERROR ) {
#    print "error : ", gluErrorString($e), "\n";
#  }
}

#------ Initialization routine

# pearl
my @mat_ambient  = ( 0.25,     0.20725,  0.20725,  0 );
my @mat_diffuse  = ( 1,        0.829,    0.829,    0 );
my @mat_specular = ( 0.296648, 0.296648, 0.296648, 0 );
my $shininess    = 0.088 * 128;

sub init {
  glClearColor( 0, 0, 0, 0 );    # Black background
  glShadeModel(GL_SMOOTH);       # Smooth shading
  glEnable(GL_MULTISAMPLE);      # Enable multisample antialiasing
  glEnable(GL_DEPTH_TEST);       # Enable hidden surface removal
  glEnable(GL_LIGHTING);

  my @light_diffuse  = ( 0.8, 0.8, 0.8, 1.0 );
  my @light_ambient  = ( 0.8, 0.8, 0.8, 1.0 );
  my @light_specular = ( 0.5, 0.5, 0.5, 1.0 );

  # set light 0
  glEnable(GL_LIGHT0);
  my @light0_position = ( -20.0, 20.0, 0.0, 0 );
  glLightfv_p( GL_LIGHT0, GL_POSITION, @light0_position );
  glLightfv_p( GL_LIGHT0, GL_DIFFUSE,  @light_diffuse );
  glLightfv_p( GL_LIGHT0, GL_AMBIENT,  @light_ambient );
  glLightfv_p( GL_LIGHT0, GL_SPECULAR, @light_specular );

  # set light 1
  glEnable(GL_LIGHT1);
  my @light1_position = ( 10.0, 50.0, 0.0, 0 );
  glLightfv_p( GL_LIGHT1, GL_POSITION, @light1_position );
  glLightfv_p( GL_LIGHT1, GL_DIFFUSE,  @light_diffuse );
  glLightfv_p( GL_LIGHT1, GL_AMBIENT,  @light_ambient );
  glLightfv_p( GL_LIGHT1, GL_SPECULAR, @light_specular );

  # set material
  glMaterialfv_p( GL_FRONT, GL_AMBIENT,  @mat_ambient );
  glMaterialfv_p( GL_FRONT, GL_DIFFUSE,  @mat_diffuse );
  glMaterialfv_p( GL_FRONT, GL_SPECULAR, @light_specular );
  glMaterialf( GL_FRONT, GL_SHININESS, $shininess );

  $shell = DrawShell();
}

#------ GLUT Callback called when the window is resized

sub reshape {
  my ( $w, $h ) = @_;
  glViewport( 0, 0, $w, $h );
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective( 45, $h ? $w / $h : 0, 1, 20 );
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
}

#------ Routine for rotating the scene

my $WaitUntil = 0;

sub spinDisplay {
  my $TimeNow = glutGet(GLUT_ELAPSED_TIME);
  if ( $TimeNow >= $WaitUntil ) {
    $spin += 1;
    $spin = $spin - 360 if ( $spin > 360 );
    glutPostRedisplay();
    $WaitUntil = $TimeNow + 1000 / 25;    # 25 frames/s
  }
}

#------ GLUT callback for the mouse

sub mouse {
  my ( $button, $state, $x, $y ) = @_;
  if ( $button == GLUT_LEFT_BUTTON ) {
    glutIdleFunc( \&spinDisplay ) if ( $state == GLUT_DOWN );
  }
  elsif ( $button == GLUT_RIGHT_BUTTON ) {
    glutIdleFunc(undef) if ( $state == GLUT_DOWN );
  }
}

#------ Main program

glutInit();
glutInitDisplayMode(
      GLUT_DOUBLE         # Double buffering
    | GLUT_RGB            # RGB color mode
    | GLUT_DEPTH          # Hidden surface removal
    | GLUT_MULTISAMPLE    # Multisample antialiasing
);
glutInitWindowSize( 300, 300 );
glutCreateWindow("Seashell");
init();
glutDisplayFunc( \&display );
glutReshapeFunc( \&reshape );
glutMouseFunc( \&mouse );
glutIdleFunc( \&spinDisplay );
glutMainLoop();

# ============ Some 3D math routines

#------ Returns the cross product of 2 vectors

sub CrossProduct {
  return (
    $_[1] * $_[5] - $_[2] * $_[4],
    $_[3] * $_[2] - $_[0] * $_[5],
    $_[0] * $_[4] - $_[1] * $_[3]
  );
}

#------ Returns the length of a vector

sub GetVectorLength {
  return sqrt( $_[0] * $_[0] + $_[1] * $_[1] + $_[2] * $_[2] );
}

#------ Returns the sum of two vectors

sub AddVectors {
  return ( $_[0] + $_[3], $_[1] + $_[4], $_[2] + $_[5] );
}

#------ Returns the vector scaled by the last parameter

sub ScaleVector {
  return ( $_[0] * $_[3], $_[1] * $_[3], $_[2] * $_[3] );
}

#------ Returns a normalized vector (length = 1)

sub NormalizeVector {
  return ( 0, 0, 0 ) if ( my $norm = GetVectorLength(@_) ) == 0;
  return ScaleVector( @_, 1 / $norm );
}

#------ Returns the unit normal vector of a triangle specified by the three
# points P1, P2, and P3.

sub FindUnitNormal {
  return NormalizeVector(
    CrossProduct(
      $_[0] - $_[3], $_[1] - $_[4], $_[2] - $_[5],    # P1-P2
      $_[3] - $_[6], $_[4] - $_[7], $_[5] - $_[8]     # P2-P3
    )
  );
}