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

#------ Base vectors (and colors)

my @i = ( 1, 0, 0 );    # red
my @j = ( 0, 1, 0 );    # green
my @k = ( 0, 0, 1 );    # blue
my @O = ( 0, 0, 0 );    # black

#------ Klein bottle parametric equations

sub r {
  my $u = shift;
  return ( 1 - cos($u) / 2 );
}

sub klein {
  my ( $u, $v ) = @_;
  return
    $u <= pi
    ? 3 / 2 * cos($u) * ( 1 + sin($u) ) + r($u) * cos($u) * cos($v)
    : 3 / 2 * cos($u) * ( 1 + sin($u) ) + r($u) * cos( $v + pi ),
    $u <= pi
    ? 4 * sin($u) + r($u) * sin($u) * cos($v)
    : 4 * sin($u),
    r($u) * sin($v);
}

my $bottle;           # OpenGL display list for the bottle
my $inc = pi / 90;    # parameters increment for $u and $v

#------ Draws the bottle

sub DrawBottle {

  my $bottle = glGenLists(1);
  glNewList( $bottle, GL_COMPILE );

  glColor3f(@i);      # red
  my @C0;             # init for $u = 0;
  for ( my $v = 0 ; $v <= 2 * pi ; $v += $inc ) {
    push @C0, [ klein( 0, $v ) ];
  }

  for ( my $u = 0 ; $u <= 2 * pi ; $u += $inc ) {
    my @C1;           # current value
    for ( my $v = 0 ; $v <= 2 * pi ; $v += $inc ) {
      push @C1, [ klein( $u, $v ) ];
    }

    glBegin(GL_TRIANGLE_STRIP);

    # draws the band between @C0 and @C1
    for ( my $i = -1 ; $i < $#C0 ; $i++ ) {
      if ( $u < pi * 1.75 ) {    # This damn surface is non-orientable!
        glNormal3d(
          FindUnitNormal( @{ $C0[$i] }, @{ $C0[ $i + 1 ] }, @{ $C1[$i] } ) );
      }
      else {
        glNormal3d(
          FindUnitNormal( @{ $C0[ $i + 1 ] }, @{ $C0[$i] }, @{ $C1[$i] } ) );
      }
      glVertex3f( @{ $C0[$i] } );
      glVertex3f( @{ $C1[$i] } );
      glVertex3f( @{ $C0[ $i + 1 ] } );
      glVertex3f( @{ $C1[ $i + 1 ] } );
    }
    glEnd();
    @C0 = @C1;    # for the next step
  }

  # draws the wireframe
  my $linewidth = 0.015;

  glColor3f(@j);                       # green
  for ( my $v = 0 ; $v <= 2 * pi ; $v += 18 * $inc ) {
    my @P1 = klein( 0, $v );
    for ( my $u = $inc ; $u <= 2 * pi ; $u += $inc ) {
      my @P2 = klein( $u, $v );
      DrawSegment( @P1, @P2, $linewidth );
      DrawPoint( @P1, $linewidth );    # junction between two segments
      @P1 = @P2;
    }
  }

  glColor3f(@k);                       # blue
  for ( my $u = 0 ; $u <= 2 * pi ; $u += 6 * $inc ) {
    my @P1 = klein( $u, 0 );
    for ( my $v = $inc ; $v <= 2 * pi ; $v += $inc ) {
      my @P2 = klein( $u, $v );
      DrawSegment( @P1, @P2, $linewidth );
      DrawPoint( @P1, $linewidth );    # junction between two segments
      @P1 = @P2;
    }
  }

  glEndList();
  return $bottle;
}

#------ Draws the scene

my $spin = 0;

sub display {
  glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );
  glLoadIdentity();
  gluLookAt( 2.0, 4.0, 10.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 );

  my $k = 0.85;
  glScalef( $k, $k, $k );
  glRotatef( $spin, 1, 0, 1 );
  glRotatef( 90, @k );
  glCallList($bottle);    # draw the bottle
  glutSwapBuffers();

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

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

sub reshape {
  my ( $w, $h ) = @_;
  glViewport( 0, 0, $w, $h );
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();    #  define the projection
  gluPerspective( 45.0, $h ? $w / $h : 0, 1.0, 20.0 );
  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.0;
    $spin = $spin - 360.0 if ( $spin > 360.0 );
    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 );
  }
}

#------ Initialization routine

my @light0_position    = ( 2.0, -8.0, -10.0, 0 );
my @mat_amb_diff_color = ( 0.8, 0.8,  0.8,   1 );
my @light_diffuse      = ( 0.3, 0.3,  0.3,   1 );
my @light_ambient      = ( 0.2, 0.2,  0.2,   1 );

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

  # Light and material
  glLightfv_p( GL_LIGHT0, GL_POSITION, @light0_position );
  glLightfv_p( GL_LIGHT0, GL_DIFFUSE,  @light_diffuse );
  glLightfv_p( GL_LIGHT0, GL_AMBIENT,  @light_ambient );
  glMaterialfv_p( GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, @mat_amb_diff_color );
  glEnable(GL_COLOR_MATERIAL);    # Material track the current color

  $bottle = DrawBottle();
}

#------ Main

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

# ====== Utility Functions

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

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

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

sub DotProduct {
  return $_[0] * $_[3] + $_[1] * $_[4] + $_[2] * $_[5];
}

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

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

#------ 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
    )
  );
}

#------ Draw a point in space
# We draw a sphere with a small radius
# Usage: DrawPoint( $x, $y, $z [, $r]);

sub DrawPoint {
  my ( $x, $y, $z, $r ) = @_;
  $r ||= 0.025;
  glPushMatrix();
  glTranslatef( $x, $y, $z );
  glutSolidSphere( $r, 10, 10 );    # The point!
  glPopMatrix();
}

#------ Draw a segment in space
# We draw a cylinder with a small radius
# Usage: DrawSegment( @A, @B [, $r]);

sub DrawSegment {
  my ( $Ax, $Ay, $Az, $Bx, $By, $Bz, $r ) = @_;
  $r ||= 0.025;
  my @u = ( $Bx - $Ax, $By - $Ay, $Bz - $Az );
  my $len = GetVectorLength(@u);
  @u = NormalizeVector(@u);
  my @v = CrossProduct( @k, @u );
  my $alpha = rad2deg acos DotProduct( @k, @u );

  # Setup the quadric object
  my $Qobj = gluNewQuadric();
  gluQuadricDrawStyle( $Qobj, GLU_FILL );
  gluQuadricNormals( $Qobj, GLU_SMOOTH );
  gluQuadricOrientation( $Qobj, GLU_OUTSIDE );
  gluQuadricTexture( $Qobj, GL_FALSE );

  glPushMatrix();
  glTranslatef( $Ax, $Ay, $Az );
  glRotatef( $alpha, @v );
  gluCylinder( $Qobj, $r, $r, $len, 10, 1 );
  glPopMatrix();
}