/// LSU EE 4702-X (Fall 2024), GPU Programming
//
 /// Homework 2 -- SOLUTION
 //
 //  Assignment: https://www.ece.lsu.edu/gpup/2024/hw02.pdf
 //
 /// Modify this file only.
 //
 /// Put Solution and other Work In Routines:
 //
 //    World::ball_setup_prob_0(), World::ball_setup_prob_1(), and
 //    World::ball_setup_prob_2().
 //
 // Based on the Links Project Code
 //

/// Purpose
//
//   Demonstrate simulation of point masses connected by springs.


/// What Code Does

// Simulates balls connected by springs over a platform. Balls and
// springs can be initialized in different arrangements (called
// scenes). consists of tiles, some are purple-tinted mirrors (showing
// a reflection of the ball), the others show the course syllabus.



///  Keyboard Commands
 //
 /// Object (Eye, Light, Ball) Location or Push
 //   Arrows (←,→,↑,↓) Page Up, Page Down
 //        Move object or push ball, depending on mode.
 //        Shift + KEY: motion is 5x faster.
 //        Ctrl + KEY : motion is 5x slower.
 //   'e': Move eye.
 //   'l': Move light.
 //   'b': Move head (first) ball. (Change position but not velocity.)
 //   'B': Push head ball. (Add velocity.)
 //
 /// Eye Direction and GUI Options
 //
 //   'Home', 'End', 'Delete', 'Insert'
 //         Turn the eye direction. Home should rotate eye direction
 //         up, End should rotate eye down, Delete should rotate eye
 //         left, Insert should rotate eye right. The eye direction
 //         vector is displayed in the upper left.
 //  'C-='  (Ctrl =) Increase green text size.
 //  'C--'  (Ctrl -) Decrease green text size.
 //  'F12'  Write screenshot to PNG file.

 /// Simulation Execution Options
 //  (Also see variables below.)
 //
 //  'p'    Pause simulation. (Press again to resume.)
 //  ' '    (Space bar.) Advance simulation by 1/30 second.
 //  'S- '  (Shift Space bar.) Advance simulation by one time step.

 /// Performance Tuning, Graphics Effects, and Debugging Options
 //
 //  'y'    Toggle value of opt_tryout1. Intended for experiments and debugging.
 //  'Y'    Toggle value of opt_tryout2. Intended for experiments and debugging.
 //  'n'    Toggle use of surface normals for spheres.
 //  'o'    Toggle shadow effect.
 //  'r'    Toggle mirror effect.
 //  '!'    Toggle rendering of spheres (balls).
 //  '@'    Toggle rendering of links.
 //  '#'    Toggle rendering of platform.


 /// Scene and Simulated System Options
 //  (Also see variables below.)
 //
 //  'w'    Twirl balls around axis formed by head and tail. (Prob 2 soln).
 //  '1'    Set up scene 1.
 //  '2'    Set up scene 2.
 //  '3'    Set up scene 3.
 //  '4'    Set up scene 4.
 //  'h'    Freeze position of first (head) ball. (Press again to release.)
 //  't'    Freeze position of last (tail) ball. (Press again to release.)
 //  's'    Stop balls.
 //  'g'    Turn gravity on and off.

 /// Variables
 //   Selected program variables can be modified using the keyboard.
 //   Use "Tab" to cycle through the variable to be modified, the
 //   name of the variable is displayed next to "VAR" on the bottom
 //   line of green text.

 //  'Tab' Cycle to next variable.
 //  '`'   Cycle to previous variable.
 //  '+'   Increase variable value.
 //  '-'   Decrease variable value.
 //
 //  VAR Spring Constant - Set spring constant.
 //  VAR Time Step Duration - Set physics time step.
 //  VAR Air Resistance - Set air resistance.
 //  VAR Light Intensity - The light intensity.
 //  VAR Gravity - Gravitational acceleration. (Turn on/off using 'g'.)

#define MAIN_INCLUDE
#undef VBUFFER_STRICT
//  #define VBUFFER_STRICT_WATCH 0xd10d270000000018
#include <vhelper.h>

#include <vstroke.h>
#include <gp/coord.h>
#include <gp/pstring.h>
#include <gp/misc.h>

#include <vutil-texture.h>
#include <vutil-pipeline.h>
#include "shapes.h"

#include <cuda_runtime.h>

#include <gp/colors.h>

#include <gp/util-containers.h>
#include <gp/coord-containers.h>
#include <gp/cuda-gpuinfo.h>
#include "links.cuh"

#include "links-shdr-common.h"


///
/// Main Data Structures
///
//
// class World: All data about scene.


class World;

enum GPU_Physics_Method { GP_cpu, GP_cuda, GP_ENUM_SIZE };
const char* const gpu_physics_method_str[] = { "CPU", "CUDA" };
const char *sh_names[] = { "PLAIN", "SET 1", "SET 2" };


// Object Holding Ball State
//
class Ball {
public:
  Ball() { init(); }
  Ball(pCoor pos){ init(); position = pos; constants_update(); }
  Ball(pCoor pos, const Ball &bd)
  { *this = bd; position = pos; constants_update(); }

  void init()
    {
      velocity = pVect(0,0,0);
      omega = pVect(0,0,0);
      radius = 1;
      mass = 0;
      mass_inv = 1;
      density = 1.00746;
      fdt_to_do = 0;
      locked = false;
      color = color_lsu_spirit_gold;
      specularity = 0;
      contact = false;
      orientation_set(pVect(0,1,0),0);
      world_time_expire = 3600.0 * 24 * 365 * 1000;
    }

  pCoor position;
  pVect velocity;
  pQuat orientation;
  pMatrix3x3p omatrix;
  pVect omega;                  // Spin rate and axis.
  double world_time_expire;

  int idx; // Position in balls_pos_rad;

  union {
    float4 m_r_fdt;
    struct {
      float mass;
      float radius;
      float fdt_to_do; // Radius / moment of inertia.
      float _;
    }; };

  float mass_min; // Mass below which simulation is unstable.
  float radius_sq, radius_inv;  // Pre-computed based on radius.
  float density;
  float mass_inv;

  bool locked;

  pVect force;
  pVect torque;
  pColor color;
  float specularity;            // Amount of shininess, in [0,1]. 
  bool contact;                 // When true, ball rendered in gray.
  float spring_constant_sum;    // Used to compute minimum mass.

  void orientation_set(pNorm dir, float angle)
    { orientation_set(pQuat(dir,angle)); }
  void orientation_set(pQuat ori)
    {
      orientation = ori;
      omatrix = pMatrix3x3_Rotation(orientation);
    }
  void constants_update();
  pVect point_rot_vel(pNorm tact_dir);
  pVect point_rot_vel(pVect tact_dir);
  void apply_tan_force(pNorm tact_dir, pNorm force_dir, double force_mag);
  void apply_tan_force(pNorm tact_dir, pVect force);
  void push(pVect amt);
  void translate(pVect amt);
  void stop();
  void freeze();
};


void
Ball::constants_update()
{
  assert( radius );
  radius_sq = radius * radius;
  if ( mass >= 0 )
    {
      mass = 4.0 / 3 * M_PI * radius_sq * radius * density;
      mass_inv = 1 / mass;
    }
  radius_inv = 1 / radius;
  // FYI, notice simplifications:
  //   const float mo_inertia = 0.4 * mass * radius_sq;
  //   fdt_to_do = radius / mo_inertia;
  fdt_to_do = radius_inv * 2.5 * mass_inv;
}

pVect
Ball::point_rot_vel(pNorm direction)
{
  /// Return velocity of point on surface of ball.
  //
  return radius * cross( omega, direction );
}

pVect
Ball::point_rot_vel(pVect rel_pos)
{
  /// Return velocity of point relative to center.
  //
  return cross( omega, rel_pos );
}

void
Ball::apply_tan_force(pNorm tact_dir, pVect force)
{
  torque += cross(tact_dir, force);
}

void
Ball::apply_tan_force(pNorm tact_dir, pNorm force_dir, double force_mag)
{
  apply_tan_force(tact_dir, force_mag * force_dir);
}

struct Link_Args {
  pColor natural_color = color_red;
  float distance_relaxed = 0;
  float spring_constant = 0;
  float link_stiffness = 0;
  float radius = 0;
};

class Link {
public:
  Link(Ball *b1, Ball *b2, Link_Args args = {}):ball1(b1),ball2(b2),
     snapped(false)
  { init(args); }
  Link(Link *link, pVect cb1p, pVect cb2p, float drp):
    ball1(link->ball1), ball2(link->ball2), cb1(cb1p), cb2(cb2p),
    is_surface_connection(true), is_renderable(false),
    is_simulatable(true), distance_relaxed(drp), spring_constant(0),
    radius(0),snapped(false){link->is_simulatable = false;}

  Link(Link *link, pVect cb1p, pVect cb2p, Link_Args args = {} ):
    ball1(link->ball1), ball2(link->ball2), cb1(cb1p), cb2(cb2p),
    is_surface_connection(true), is_renderable(false),
    is_simulatable(true),
    distance_relaxed(args.distance_relaxed),
    spring_constant(args.spring_constant),
    radius(0),snapped(false){link->is_simulatable = false;}

  void init(Link_Args& args)
  {
    assert( args.distance_relaxed == 0 );
    natural_color = args.natural_color;
    spring_constant = args.spring_constant;
    radius = args.radius;
    init();
  }
  void init()
    {
      assert( ball1->radius > 0 );
      assert( ball2->radius > 0 );
      is_simulatable = true;
      is_renderable = true;
      pNorm n12(ball1->position,ball2->position);
      const float rad_sum = ball1->radius + ball2->radius;
      is_surface_connection = n12.magnitude >= rad_sum;

      pMatrix3x3 b1rot = ball1->omatrix;
      pMatrix3x3 b2rot = ball2->omatrix;

      pVect b1_y = b1rot * pVect(0,1,0);
      pVect b1_x = b1rot * pVect(1,0,0);
      bool b1_dir_is_y = fabs(dot(b1_y,n12)) < 0.999;
      pVect b1_dir_loc = b1_dir_is_y ? b1_y : b1_x;

      con_x = pNorm(cross(b1_dir_loc,n12));
      con_y = cross(n12,con_x);
      b1_dir = mult_MTV( b1rot, con_y );
      b2_dir = mult_MTV( b2rot, con_y );

      if ( !is_surface_connection )
        {
          distance_relaxed = n12.magnitude;
          cb1 = cb2 = pVect(0,0,0);
          is_renderable = false;
          return;
        }

      distance_relaxed = n12.magnitude - rad_sum;
      cb1 = ball1->radius * mult_MTV( b1rot,  n12 );
      cb2 = ball2->radius * mult_MTV( b2rot, -n12 );
    }
  int serial;
  int idx;  // Position in lis_pos1, etc.
  Ball* const ball1;
  Ball* const ball2;
  union {
    int2 ball_idx;
    struct { int ball1_idx, ball2_idx; };
  };
  pVect cb1, cb2, con_x, con_y, b1_dir, b2_dir;
  bool is_surface_connection;
  bool is_renderable, is_simulatable;
  union {
    float2 dr_sc;
    struct {
      float distance_relaxed;
      float spring_constant;
    };
  };
  float radius;
  bool snapped;
  pColor natural_color;
  pColor color;

  pVect spring_force_12;
  pVect torque1,torque2;

};

class Surface {
public:
  Ball *b0, *b1, *b2;
  pColor color = color_blue;
  pTCoor tc0, tc1, tc2;
  union {
    int4 b_idx;
    struct { int b0_idx, b1_idx, b2_idx, _; };
  };
  Surface* operator ->() { return this; }
};

// Declare containers for Balls and Links. The containers are derived
// from std::vector and use overloaded operator += for push_back. (See
// include/util-containers.h.)
//
typedef pVectorI<Link> Links;
typedef pVectorI<Ball> Balls;

pVect vec_rand() {return pVect(drand48(),drand48(),drand48());}

class HW_Info
{
public:
  HW_Info(){set();};
  int n_balls, n_lines;
  float ball_radius;
  pNorm axis;
  pCoor center_pos,  line_0_bot_pos, line_0_top_pos;

  void set()
  {
    center_pos = pCoor(12.2,4.8f,0) + 2*vec_rand();
    n_balls = 4 + random()%8;
    n_lines = 32 + random()%20;

    pNorm az = pVect( drand48(), 10, drand48() );
    axis = az;
    pNorm ax =
      fabs(az.x) < fabs(az.y) ? pVect(0,az.z,-az.y) : pVect(az.z,0,-az.x);
    pNorm ay = cross(az,ax);

    float r_bot = 10 + 10 * drand48();
    ball_radius = 0.2;
    line_0_bot_pos = center_pos + r_bot * ax;
    line_0_top_pos =
      center_pos + ( 10 + 10 * drand48() ) * az
      + drand48() * r_bot * ax + 3*(drand48()-0.5) * r_bot * ay;
  }
};

enum Shader_Option { SO_Plain, SO_Instances, SO_True, SO_ENUM_SIZE };

class World {
public:
  World(pVulkan_Helper &fb, int argc, char **argv)
    :vh(fb), ff_state(fb.qs), raytrace(fb.qs),
     frame_timer(vh.frame_timer), transform(fb.qs),
     shapes(ff_state), sphere(ff_state),
     opt_tryout1(fb.qs.opt_tryout1),opt_tryout2(fb.qs.opt_tryout2)
  {init(argc,argv);}
  void init(int argc, char **argv);
  void init_graphics();
  void run();
  void frame_callback();
  void render(vk::CommandBuffer& cb);
  void render_objects(vk::CommandBuffer& cb, Render_Option render_option);
  void objects_erase();
  void cb_keyboard();
  void modelview_update();

  pVulkan_Helper& vh;
  VFixed_Function_State_Manager ff_state;
  VRaytrace raytrace;
  pFrame_Timer& frame_timer;
  VTransform transform;
  Shapes shapes;
  Sphere sphere;

  bool opt_want_raytrace;

  pVariable_Control variable_control;
  double world_time;
  double last_frame_wall_time;
  float opt_time_step_duration;
  int time_step_count;
  float opt_gravity_accel;      // Value chosen by user.
  pVect gravity_accel;          // Set to zero when opt_gravity is false;
  bool opt_gravity;
  bool opt_head_lock, opt_tail_lock;
  bool &opt_tryout1, &opt_tryout2;  // For ad-hoc experiments.
  float opt_tryoutf;
  bool opt_sphere_true;
  VBufferV<Shdr_Uni_Common> buf_uni_common;

  // Tiled platform for ball.
  //
  float platform_xmin, platform_xmax, platform_zmin, platform_zmax;
  float platform_pi_xwidth_inv;
  VVertex_Buffer_Set bset_platform_tex, bset_platform_mir;
  VPipeline pipe_platform_mir, pipe_platform_tex;
  RT_Geo_Package rt_geo_platform_tex, rt_geo_platform_mir;

  vk::Sampler sampler;
  VTexture texid_hw;
  VTexture texid_syl;
  VTexture texid_emacs;

  void platform_update();
  bool platform_collision_possible(pCoor pos);

  pCoor light_location;
  float opt_light_intensity;
  enum { MI_Eye, MI_Light, MI_Ball, MI_Ball_V, MI_COUNT } opt_move_item;
  bool opt_pause;
  bool opt_single_frame;      // Simulate for one frame.
  int opt_single_time_step;  // Simulate for one time step.

  VBufferV<Uni_Lighting> uni_light, uni_light_ambient, uni_light_all;
  VBufferV<Uni_Lighting> uni_light_reflected;
  VBufferV<Uni_Lighting> *puni_light_curr;

  pCoor eye_location;
  pVect eye_direction;
  pMatrix transform_projection;
  pMatrix modelview;

  pVect adj_vector;
  double adj_t_prev;
  double adj_t_stop;
  double adj_duration;

  uint opt_shader;

  bool opt_mirror;
  int opt_mirror_depth;
  bool opt_shadows;
  bool opt_shadow_volumes;      // Make shadow volumes visible.
  bool opt_normal_sphere;
  int opt_hide_stuff;

  int last_setup; // Last scene set up.

  void ball_setup_1(bool big);  // Silly Tree
  void ball_setup_2(); // Hanging thing.
  void ball_setup_3(); // Corona virus particle.
  void ball_setup_4(); // Wheel.
  void ball_setup_5(); // Top.
  void ball_setup_6(); // Wheel with webbing between spokes.
  void ball_setup_7();
  void setup_at_end();
  void time_step_cpu(double);
  void balls_stop();
  void balls_freeze();
  void balls_translate(pVect amt, int idx);
  void balls_translate(pVect amt);
  void balls_push(pVect amt, int idx);
  void balls_push(pVect amt);
  void balls_twirl();
  void lock_update();

  Ball *make_marker(pCoor pos, pColor col, float radius = 0.2 );
  Balls balls_marker; // Written to balls right after setup.

  float opt_spring_constant, opt_spring_damp;
  float opt_friction_dynamic;
  float opt_air_resistance;
  float distance_relaxed;
  int chain_length;
  Balls balls;
  Ball *head_ball, *tail_ball;
  Links links;
  vector<Surface> surfaces;
  int opt_sphere_slices;
  float opt_omega;  // For disk and top.

  Links link_new
  ( Ball *ball1, Ball *ball2,
    float stiffness = 0.2, float spring_constant = 0);
  Links link_new
  ( Ball *ball1, Ball *ball2, Link_Args largs );

  /// CUDA Stuff

  int opt_physics_method;
  GPU_Info gpu_info;
  cudaEvent_t frame_start_ce, frame_stop_ce;
  int opt_block_size;
  bool gpu_const_data_stale;
  bool cpu_phys_data_stale;

  void data_cpu_to_gpu_constants();
  void data_cpu_to_gpu_dynamic();
  void data_gpu_to_cpu_dynamic();
  CPU_GPU_Common c;  // c is for common.  See links.cuh.

  void render_link_start();
  void render_link_gather(Link *link);
  void render_link_render(vk::CommandBuffer& cb, bool shadow_volumes = false);
  void render_link_render_sv(vk::CommandBuffer& cb)
  { render_link_render(cb,true); };

  VPipeline pipe_set_2, pipe_set_2_sv;
  double world_time_link_update;

  VCompute comp_links;
  VBuffer<int> buf_links_as_idx; // Indices for links' acceleration structures.
  VBuffer<vec2> buf_links_tex_coor;
  VBuffer<vec4> buf_links_vertex_g;
  VBuffer<vec4> buf_links_normal_g;
  RT_Geo_Package rt_geo_links;

  int opt_sides, opt_segments;

  vector< VBufferVV<vec4>* > lis_bufs;
  VBufferVV<vec4> lis_pos1, lis_pos2;
  VBufferVV<vec4> lis_v1, lis_v2;
  VBufferVV<vec4> lis_b1_ydir, lis_b2_ydir, lis_tex_rad;

  VVertex_Buffer_Set bset_surfaces;
  VPipeline pipe_surfaces;
  RT_Geo_Package rt_geo_surfaces;
  void render_surfaces(vk::CommandBuffer& cb);

  /// Stuff for 2024 Homework 2

  void ball_setup_hw2( char version );
  void ball_setup_prob_0();
  void ball_setup_prob_1();
  void ball_setup_prob_2();

  HW_Info h;
  bool opt_freeze_config;

};

void
World::ball_setup_hw2( char version )
{
  last_setup = version - '0';
  objects_erase();

  /// Homework 2: Do not modify this routine.
  // Your solution and other work goes in World::ball_setup_prob_0(),
  // World::ball_setup_prob_1(), and World::ball_setup_prob_2().

  if ( !opt_freeze_config ) h.set();

  switch ( version ){
  case '0': ball_setup_prob_0(); break;
  case '1': ball_setup_prob_1(); break;
  case '2': ball_setup_prob_2(); break;
  default: assert( false );
  }
  setup_at_end();
}

enum { OH_Links = 1, OH_Platform = 2, OH_Sphere = 4 };

void
World::init_graphics()
{
  vh.want_raytrace().init();
  vh.opt_record_every_time = true;
  vh.display_cb_set([&](){frame_callback();});
  vh.cbs_cmd_record.push_back( [&](vk::CommandBuffer& cb){ render(cb); });
  sphere.ppuni_light = &puni_light_curr;
  sphere.transform = &transform;

  opt_want_raytrace = true;

  if ( vh.have_raytrace )
    {
      raytrace.init();
      for ( auto b: { &bset_platform_mir, &bset_platform_tex, &bset_surfaces } )
        b->usage_set
          ( vk::BufferUsageFlagBits::eVertexBuffer
            | vk::BufferUsageFlagBits::eStorageBuffer );
      buf_links_as_idx.init
        ( vh.qs, vk::BufferUsageFlagBits::eStorageBuffer, 0,
          vk::MemoryPropertyFlagBits::eDeviceLocal );
      buf_links_tex_coor.init
        ( vh.qs, vk::BufferUsageFlagBits::eStorageBuffer, 0,
          vk::MemoryPropertyFlagBits::eDeviceLocal );
      buf_links_vertex_g.init
        ( vh.qs, vk::BufferUsageFlagBits::eStorageBuffer, 0,
          vk::MemoryPropertyFlagBits::eDeviceLocal );
      buf_links_normal_g.init
        ( vh.qs, vk::BufferUsageFlagBits::eStorageBuffer, 0,
          vk::MemoryPropertyFlagBits::eDeviceLocal );
    }

  buf_uni_common.init(vh.qs,vk::BufferUsageFlagBits::eUniformBuffer);
  sphere.buf_uni_common = buf_uni_common;

  lis_bufs =
    { &lis_pos1, &lis_pos2,
      &lis_v1, &lis_v2,
      &lis_b1_ydir, &lis_b2_ydir,
      &lis_tex_rad };

  for ( auto b: lis_bufs )
    b->init(vh.dev_phys, vh.dev, vk::BufferUsageFlagBits::eStorageBuffer );

  ///
  /// Graphical Model Initialization
  ///

  opt_head_lock = false;
  opt_tail_lock = false;
  opt_tryout1 = opt_tryout2 = false;
  opt_normal_sphere = true;
  opt_sphere_true = true;
  opt_tryoutf = 0;
  variable_control.insert_linear(opt_tryoutf,"Tryout F",0.1);

  eye_location = pCoor(24.2,11.6,-38.7);
  eye_direction = pVect(-0.42,-0.09,0.9);

  platform_xmin = -40; platform_xmax = 40;
  platform_zmin = -40; platform_zmax = 40;

  texid_syl.init
    ( vh.qs, P_Image_Read( vh.gp_root_get() / "vulkan" / "gpup.png", 255 ) );
  texid_emacs.init( vh.qs, P_Image_Read("mult.png",-1) );
  texid_hw.init( vh.qs, P_Image_Read("hw05-assign.png",255) );
  sampler = vh.qs.dev.createSampler
    ( { {},
        vk::Filter::eLinear,  // magFilter
        vk::Filter::eLinear,  // minFilter
        vk::SamplerMipmapMode::eLinear,
        // Also: eRepeat, eMirroredRepeat, eClampToEdge, etc.
        vk::SamplerAddressMode::eRepeat,
        vk::SamplerAddressMode::eRepeat,
        vk::SamplerAddressMode::eRepeat,
        0.0f,  // mipLodBias
        true, // anisotropyEnable,
        16.0f, // maxAnisotropy
        false, vk::CompareOp::eNever,  // compareEnable, compareOp
        0.0f, VK_LOD_CLAMP_NONE,  // min and max LOD
        vk::BorderColor::eFloatOpaqueBlack,
        false // unnormalizedCoordinates
        } );

  sphere.texture_set(sampler,texid_emacs);

  opt_light_intensity = 50;
  light_location = pCoor(platform_xmax,platform_xmax,platform_zmin);
  uni_light.init(vh.qs,vk::BufferUsageFlagBits::eUniformBuffer);
  uni_light_all.init(vh.qs,vk::BufferUsageFlagBits::eUniformBuffer);
  uni_light_reflected.init(vh.qs,vk::BufferUsageFlagBits::eUniformBuffer);
  uni_light_ambient.init(vh.qs,vk::BufferUsageFlagBits::eUniformBuffer);

  variable_control.insert(opt_light_intensity,"Light Intensity");

  opt_move_item = MI_Eye;
  opt_pause = false;
  opt_single_time_step = 0;
  opt_single_frame = false;

  platform_update();
  modelview_update();

  adj_t_stop = 0;
  adj_duration = 0.25;

  opt_sides = 20;
  opt_segments = 16;
  //  variable_control.insert(opt_sides,"Link Sides");
  //  variable_control.insert(opt_segments,"Link Segments",5,1);
  opt_hide_stuff = 0;

  if ( vh.have_raytrace )
    {
      raytrace.buf_light_set(uni_light_all);
    }

  opt_freeze_config = false;
}

void
World::run()
{
  vh.message_loop_spin();

  uni_light.destroy();
  uni_light_ambient.destroy();
  uni_light_all.destroy();
  uni_light_reflected.destroy();
  transform.destroy();
  buf_uni_common.destroy();

  for ( auto b: lis_bufs ) b->destroy();

  pipe_platform_tex.destroy(); pipe_platform_mir.destroy();
  bset_platform_tex.destroy(); bset_platform_mir.destroy();
  pipe_set_2.destroy();
  pipe_set_2_sv.destroy();

  texid_syl.destroy();
  texid_emacs.destroy();
  texid_hw.destroy();
  vh.dev.destroySampler(sampler);
  shapes.destroy();
  sphere.destroy();

  bset_surfaces.destroy();
  pipe_surfaces.destroy();
  rt_geo_surfaces.destroy();

  buf_links_as_idx.destroy();
  buf_links_tex_coor.destroy();
  buf_links_normal_g.destroy();
  buf_links_vertex_g.destroy();
  rt_geo_platform_tex.destroy();
  rt_geo_platform_mir.destroy();
  rt_geo_links.destroy();
  comp_links.destroy();
  raytrace.destroy();

  vh.finish();
}

void
World::render_surfaces(vk::CommandBuffer& cb)
{
  if ( !pipe_surfaces )
    pipe_surfaces
      .init( vh.qs )
      .ds_follow( puni_light_curr )
      .ds_follow( transform )
      .ds_uniform_use( "BIND_UNI_COMMON", buf_uni_common )
      .topology_set( vk::PrimitiveTopology::eTriangleList )
      .create();

  bset_surfaces.reset( pipe_surfaces );

  for ( Surface& s: surfaces )
    {
      Ball &b0 = *s.b0, &b1 = *s.b1, &b2 = *s.b2;
      pNorm snorm = cross(b0.position,b1.position,b2.position);
      bset_surfaces << b0.position << b1.position << b2.position;
      bset_surfaces << s.color << s.color << s.color;
      bset_surfaces << snorm << snorm << snorm;
    }

  bset_surfaces.to_dev();

  if ( vh.have_raytrace && !rt_geo_surfaces )
    rt_geo_surfaces.init( raytrace );

  if ( vh.raytrace_active() )
    {
      rt_geo_surfaces.bind( bset_surfaces, RT_Geometry_Triangles );
      return;
    }

  pipe_surfaces.record_draw( cb, bset_surfaces );
}

void
World::platform_update()
{
  const float tile_count = 19;
  const float ep = 1.00001;
  const float delta_x = ( platform_xmax - platform_xmin ) / tile_count * ep;
  const float zdelta = ( platform_zmax - platform_zmin ) / tile_count * ep;

  const float trmin = 0.05;
  const float trmax = 0.7;
  const float tsmin = 0;
  const float tsmax = 0.4;
  pColor color_platform( color_white * 0.9 );

  if ( !pipe_platform_mir )
    pipe_platform_mir
      .init( ff_state )
      .ds_set( color_lsu_spirit_purple )
      .ds_follow( puni_light_curr )
      .ds_follow( transform )
      .topology_set( vk::PrimitiveTopology::eTriangleList )
      .create();
  bset_platform_mir.reset( pipe_platform_mir ).tolerate_tcoords();

  if ( !pipe_platform_tex )
    pipe_platform_tex
      .init( ff_state )
      .ds_set_material(color_platform)
      .ds_follow( puni_light_curr )
      .ds_follow( transform )
      .topology_set( vk::PrimitiveTopology::eTriangleList )
      .ds_use( sampler, texid_syl )
      .create();
  bset_platform_tex.reset( pipe_platform_tex );

  const pVect platform_normal(0,1,0);
  pTCoor t00(trmin,tsmin), t01(trmin,tsmax), t10(trmax,tsmin), t11(trmax,tsmax);
  bool even = true;

  for ( int i = 0; i < tile_count; i++ )
    {
      const float x0 = platform_xmin + i * delta_x;
      const float x1 = x0 + delta_x;
      const float y = 0;
      for ( float z = platform_zmin; z < platform_zmax; z += zdelta )
        {
          auto& bset = even ? bset_platform_tex : bset_platform_mir;
          const float z1 = z+zdelta;
          pCoor p00(x0,y,z), p01(x0,y,z1), p10(x1,y,z), p11(x1,y,z1);

          bset << t11 << platform_normal << p00;
          bset << t10 << platform_normal << p01;
          bset << t00 << platform_normal << p11;

          bset << t11 << platform_normal << p00;
          bset << t00 << platform_normal << p11;
          bset << t01 << platform_normal << p10;

          even = !even;
        }
    }

  bset_platform_mir.to_dev();
  bset_platform_tex.to_dev();

  if ( vh.have_raytrace && !rt_geo_platform_mir )
    {
      pColor mirror_color(color_lsu_spirit_purple);
      mirror_color.a = -0.5;
      rt_geo_platform_mir
        .init( raytrace )
        .color_front_set( mirror_color )
        .bind( bset_platform_mir, RT_Geometry_Triangles );
      rt_geo_platform_tex
        .init( raytrace )
        .color_front_set( color_platform )
        .bind( VR_BUF_IDX_sampler, sampler, texid_syl )
        .bind( bset_platform_tex, RT_Geometry_Triangles );
    }
}

void
World::modelview_update()
{
  pMatrix_Translate center_eye(-eye_location);
  pMatrix_Rotation rotate_eye(eye_direction,pVect(0,0,-1));
  modelview = rotate_eye * center_eye;
}

void
World::render_objects(vk::CommandBuffer& cb, Render_Option option)
{
  //  const float shininess_ball = 5;
  pColor spec_color(1,1,1);
  //  pColor spec_light = color_white * opt_light_intensity;

  const bool hide_links = opt_hide_stuff & OH_Links;
  const bool hide_platform = opt_hide_stuff & OH_Platform;
  const bool hide_sphere = opt_hide_stuff & OH_Sphere;

  sphere.init(opt_sphere_slices);

  if ( option == RO_Shadow_Volumes )
    {
      if ( !hide_sphere )
        {
          sphere.light_pos = light_location;
          sphere.render_bunch_gather_sv(world_time);
          for ( Ball *ball: balls )
            sphere.render_shadow_volume(cb,ball->radius,ball->position);

          sphere.render_bunch_render_sv(cb);
        }

      if ( !hide_links )
        {
          render_link_start();
          for ( Link *link: links )
            if ( link->is_renderable ) render_link_gather(link);
          render_link_render_sv(cb);
        }
    }
  else
    {
      sphere.rt_geo_sphere.rt_wanted_set(!hide_sphere);
      if ( !hide_sphere )
        {
          sphere.render_bunch_gather(world_time);
          for ( Ball *ball: balls )
            {
              if ( ball->contact )
                sphere.color = color_gray;
              else if ( ball->mass > 0 && ball->mass < ball->mass_min )
                sphere.color = color_red;
              else if ( ball->mass > 0 && ball->locked )
                sphere.color = color_pale_green;
              else
                sphere.color = ball->color;

              sphere.color.a = ball->specularity ? -ball->specularity : 1;
              sphere.render
                (cb, ball->radius,ball->position,
                 pMatrix_Rotation(ball->orientation));
            }

          sphere.render_bunch_render(cb, opt_sphere_true);
        }

      rt_geo_links.rt_wanted_set(!hide_links);
      if ( !hide_links )
        {
          render_link_start();
          for ( Link *link: links )
            if ( link->is_renderable ) render_link_gather(link);
          render_link_render(cb);
        }
    }

  rt_geo_platform_tex.rt_wanted_set(!hide_platform);
  rt_geo_platform_mir.rt_wanted_set(!hide_platform);

  if ( option != RO_Shadow_Volumes )
    render_surfaces(cb);

  if ( vh.raytrace_active() ) return;
  if ( hide_platform ) return;

  if ( option == RO_Shadow_Volumes ) return;

  //
  // Render Platform
  //
  pipe_platform_tex.record_draw(cb,bset_platform_tex);
}

void
World::render(vk::CommandBuffer& cb)
{
  const bool do_raytrace = vh.raytrace_active();

  /// Emit a Graphical Representation of Simulation State
  //

  // Understanding of the code below not required for introductory
  // lectures.

  // That said, much of the complexity of the code is to show
  // the ball shadow and reflection.


  const int win_width = vh.get_width();
  const int win_height = vh.get_height();
  const float aspect = float(win_width) / win_height;
  const float n_dist = 0.01;
  const float xr = .8 * n_dist;

  transform.eye_from_global_set( modelview );
  transform.clip_from_eye_set
    ( pMatrix_Frustum(-xr,xr,-xr/aspect,xr/aspect,n_dist,5000) );
  // pMatrix_Frustum: left, right, bottom, top, near, far

  vh.clearValues[0].color =
    vk::ClearColorValue( array<float, 4>( { { 0, 0, 0, 0.5 } } ) );

  /// Lighting
  //
  auto& light_all = uni_light_all->cgl_LightSource[0];

  light_all.position = transform.eye_from_global * light_location;

  light_all.constantAttenuation = 0.5;
  light_all.linearAttenuation = 1.0;
  light_all.quadraticAttenuation = 0;

  pColor ambient_color(0x555555);

  uni_light_all->cgl_LightModel.ambient = ambient_color;

  light_all.diffuse = opt_light_intensity * color_white;
  light_all.ambient = color_black;
  light_all.specular = opt_light_intensity * color_white;

  uni_light_reflected = uni_light_all.get_val();
  uni_light_reflected->cgl_LightSource[0].position =
    transform.eye_from_global * pMatrix_Scale(1,-1,1) * light_location;

  uni_light_ambient = uni_light_all.get_val();
  auto& light_amb = uni_light_ambient->cgl_LightSource[0];
  light_amb.diffuse = color_black;
  light_amb.specular = color_black;

  uni_light = uni_light_all.get_val();
  uni_light->cgl_LightModel.ambient = color_black;

  uni_light.to_dev();
  uni_light_ambient.to_dev();
  uni_light_all.to_dev();
  uni_light_reflected.to_dev();

  puni_light_curr = &uni_light_all;

  buf_uni_common->tryout =
    ivec4(opt_tryout1,opt_tryout2,opt_normal_sphere,false);
  buf_uni_common->tryoutf = vec4(opt_tryoutf, 0, 0, 0);
  buf_uni_common->light_pos_g = light_location;
  if ( do_raytrace )
    {
      raytrace.buf_uni_common->opt_tryout = buf_uni_common->tryout;
      raytrace.buf_uni_common->opt_tryoutf = buf_uni_common->tryoutf;
      raytrace.buf_uni_common->opt_shadows = opt_shadows;
      raytrace.buf_uni_common->opt_mirror = opt_mirror;
      raytrace.buf_light_set(uni_light_all);
    }

  const double time_now = time_wall_fp();
  const bool blink_visible = int64_t(time_now*3) & 1;
# define BLINK(txt,pad) ( blink_visible ? txt : pad )

  vh.fbprintf("%s\n",frame_timer.frame_rate_text_get());

  vh.fbprintf
    ("Raytrace: %s ('R')   Compiled: %s\n",
     !vh.have_raytrace ? "NA" :
     opt_want_raytrace ? ( raytrace.rt_ready ? "ON " : "???" ) : "OFF",
#ifdef __OPTIMIZE__
     "WITH OPTIMIZATION"
#else
     BLINK("WITHOUT OPTIMIZATION","")
#endif
     );

  vh.fbprintf
    ("Time Step: %8d  World Time: %11.6f  %s\n",
     time_step_count, world_time,
     opt_pause ? BLINK("PAUSED, 'p' to unpause, SPC or S-SPC to step.","") :
     "Press 'p' to pause."
     );

  vh.fbprintf
    ("Eye location: [%5.1f, %5.1f, %5.1f]  "
     "Eye direction: [%+.2f, %+.2f, %+.2f] ",
     eye_location.x, eye_location.y, eye_location.z,
     eye_direction.x, eye_direction.y, eye_direction.z);

  if ( balls )
    {
      Ball& ball = *balls[0];

      vh.fbprintf
        ("Head Ball Pos  [%5.1f,%5.1f,%5.1f] Vel [%+5.1f,%+5.1f,%+5.1f]\n",
         ball.position.x,ball.position.y,ball.position.z,
         ball.velocity.x,ball.velocity.y,ball.velocity.z );
    }

  vh.fbprintf
    ("Scene: %d ('0-7')  Config: %s ('c') "
     "Tryout 1: %s  ('y')  Tryout 2: %s  ('Y')\n",
     last_setup, opt_freeze_config ? BLINK("FROZEN","      ") : "RANDOM",
     opt_tryout1 ? BLINK("ON ","   ") : "OFF",
     opt_tryout2 ? BLINK("ON ","   ") : "OFF" );

#if 0
  vh.fbprintf
    ("Hide: %c%c%c ('!@#')  Effect: %c%c ('or')  Sphere: %s  ('z')  "
     "Tryout 1: %s  ('y')  Tryout 2: %s  ('Y')  Normals: %s ('n')\n",
     opt_hide_stuff & OH_Sphere ? 'S' : '_',
     opt_hide_stuff & OH_Links ? 'L' : '_',
     opt_hide_stuff & OH_Platform ? 'P' : '_',
     opt_shadows ? 'S' : '_',
     opt_mirror ? 'M' : '_',
     opt_sphere_true ? "TRUE" : "TRI" ,
     opt_tryout1 ? BLINK("ON ","   ") : "OFF",
     opt_tryout2 ? BLINK("ON ","   ") : "OFF",
     opt_normal_sphere ? "SPHERE" : "TRI");
#endif

  pVariable_Control_Elt* const cvar = variable_control.current;
  vh.fbprintf("VAR %s = %.5f  (TAB or SH-TAB to change, +/- to adjust)\n",
                      cvar->name,cvar->get_val());

  // Render Marker for Light Source
  //
  shapes.record_tetrahedron(cb,transform,light_location,0.15);

  //
  // Render ball reflection.  (Will be blended with dark tiles.)
  //

  const bool hide_platform = opt_hide_stuff & OH_Platform;

  if ( do_raytrace )
    {
      raytrace.tform_update(transform);
      render_objects(cb,RO_Normally);
      raytrace.buf_uni_common.to_dev();

      return;
    }

  buf_uni_common.to_dev();

  vk::ColorComponentFlags clr_all
    ( vk::ColorComponentFlagBits::eR | vk::ColorComponentFlagBits::eG |
      vk::ColorComponentFlagBits::eB | vk::ColorComponentFlagBits::eA );

  if ( !hide_platform )
    {
      // Save current state of fixed-functionality state.
      //
      ff_state.push();

      if ( opt_mirror )
        {
          /// Mirror Effect
          //
          //  - Step 1:  Write location of mirrors into stencil buffer.
          //
          //  - Step 2a: Prepare a matrix that transforms to mirrored location.
          //  - Step 2b: Set up blending to blend mirrored objects with mirrors.
          //
          //  - Step 3:  Render all objects except mirrors.

          ///  Step 1: Write location of mirrors into stencil buffer.
          //
          //   Write stencil value 2 at location of dark (mirrored) tiles.
          //

          vk::StencilOpState s_op_apply
            ( vk::StencilOp::eReplace, vk::StencilOp::eKeep,
              vk::StencilOp::eKeep, vk::CompareOp::eNever, 2, 2, 2 );
          // sfail, spass, dfail, compare;  cmask, wmask, ref

          ff_state.pipelineDepthStencilStateCreateInfo_get()
            .setStencilTestEnable(true)
            .setFront(s_op_apply).setBack(s_op_apply);

          // Render mirrored tiles (but only to write stencil buffer).
          //
          pipe_platform_mir.record_draw(cb,bset_platform_mir);

          cb.nextSubpass(vk::SubpassContents::eInline); vh.qs.subpass++;

          // Prepare to write only stenciled locations.
          //
          vk::StencilOpState s_op_use
            ( vk::StencilOp::eKeep, vk::StencilOp::eKeep,
              vk::StencilOp::eKeep, vk::CompareOp::eEqual, 2, 2, 2 );
          // sfail, spass, dfail, compare;  cmask, wmask, ref
          //
          // Stencil test passes if stencil value = 2.
          // Don't change stencil buffer.

          ff_state.pipelineDepthStencilStateCreateInfo_get()
            .setFront(s_op_use).setBack(s_op_use);

          ///  Step 2a: Prepare a matrix that transforms to mirrored location.
          //   Use a transform that reflects objects to other side of platform.
          //
          const pMatrix prev = transform.eye_from_global;
          const pMatrix_Scale reflect(1,-1,1);
          transform.eye_from_global = prev * reflect;

          // Reflected front face should still be treated as the front face.
          //
          ff_state.pipelineRasterizationStateCreateInfo_get()
            .setFrontFace(vk::FrontFace::eClockwise);

          puni_light_curr = &uni_light_reflected;

          /// Step 3:  Render all objects.

          render_objects(cb,RO_Mirrored);

          transform.eye_from_global = prev;

          ff_state.pipelineRasterizationStateCreateInfo_get()
            .setFrontFace(vk::FrontFace::eCounterClockwise);
        }

      // blend operation, blend factors, blend constants
      ff_state.pipelineColorBlendStateCreateInfo_get()
        .setBlendConstants( array<float,4>{.5,.5,.5,.5} );

      ff_state.pipelineColorBlendAttachmentState_get() =
        vk::PipelineColorBlendAttachmentState
        ( true,                   // blendEnable
          vk::BlendFactor::eConstantAlpha,  // srcColorBlendFactor
          vk::BlendFactor::eConstantAlpha,  // dstColorBlendFactor
          vk::BlendOp::eAdd,       // colorBlendOp
          vk::BlendFactor::eOne,   // srcAlphaBlendFactor
          vk::BlendFactor::eZero,  // dstAlphaBlendFactor
          vk::BlendOp::eAdd,       // alphaBlendOp
          clr_all      // colorWriteMask
          );

      puni_light_curr = &uni_light_all;;
      pipe_platform_mir.record_draw(cb,bset_platform_mir);

      ff_state.pop();
    }

  if ( !opt_shadows )
    {
      // Render.
      //
      render_objects(cb,RO_Normally);
    }
  else
    {
      //
      // First pass, render using only ambient light.
      //
      puni_light_curr = &uni_light_ambient;

      // Send balls, tiles, and platform to opengl.
      //
      render_objects(cb,RO_Normally);

      //
      // Second pass, find un-shadowed areas and add on light0.
      //


      /// Write Shadow Volumes
      //
      vk::ClearAttachment clr_att
        ( vk::ImageAspectFlagBits::eStencil, 0,
          vk::ClearDepthStencilValue(0,0) );
      vk::ClearRect clr_rect( vk::Rect2D({}, vh.s_extent), 0, 1 );
      cb.clearAttachments( clr_att, clr_rect );

      ff_state.push();

      vk::StencilOpState s_op_apply_f
        ( vk::StencilOp::eKeep, vk::StencilOp::eIncrementAndWrap,
          vk::StencilOp::eKeep, vk::CompareOp::eAlways, -1, -1, 0 );
      vk::StencilOpState s_op_apply_b
        ( vk::StencilOp::eKeep, vk::StencilOp::eDecrementAndWrap,
          vk::StencilOp::eKeep, vk::CompareOp::eAlways, -1, -1, 0 );
      // sfail, spass, dfail, compare; c-mask, w-mask, ref
      // ref is -1 (all 1's) because we want to operate on all bits in
      // stencil value.

      ff_state.pipelineDepthStencilStateCreateInfo_get()
        .setStencilTestEnable(true)
        .setDepthWriteEnable(false)
        .setFront(s_op_apply_f).setBack(s_op_apply_b);

      if ( !opt_shadow_volumes )
        {
          ff_state.pipelineColorBlendAttachmentState_get()
            .setColorWriteMask( {} );
        }
      else
        {
          ff_state.pipelineColorBlendAttachmentState_get() =
            vk::PipelineColorBlendAttachmentState
            ( true,                   // blendEnable
              vk::BlendFactor::eConstantAlpha,  // srcColorBlendFactor
              vk::BlendFactor::eConstantAlpha,  // dstColorBlendFactor
              vk::BlendOp::eAdd,       // colorBlendOp
              vk::BlendFactor::eOne,   // srcAlphaBlendFactor
              vk::BlendFactor::eZero,  // dstAlphaBlendFactor
              vk::BlendOp::eAdd,       // alphaBlendOp
              clr_all      // colorWriteMask
              );

          ff_state.pipelineColorBlendStateCreateInfo_get()
            .setBlendConstants( array<float,4>{.5,.5,.5,.5} );

          ff_state.pipelineDepthStencilStateCreateInfo_get()
            .setDepthWriteEnable(true);
        }

      // Write stencil with shadow locations based on shadow volumes
      // cast by light0 (light_location).  Shadowed locations will
      // have a positive stencil value.
      //
      render_objects(cb,RO_Shadow_Volumes);

      cb.nextSubpass(vk::SubpassContents::eInline); vh.qs.subpass++;

      // Use stencil test to prevent writes to shadowed areas.
      //
      vk::StencilOpState s_op_apply_fb
        ( vk::StencilOp::eKeep, vk::StencilOp::eKeep,
          vk::StencilOp::eKeep, vk::CompareOp::eEqual, -1, -1, 0 );
      // sfail, spass, dfail, compare;  c-mask, w-mask, ref
      //
      // Stencil test passes if stencil value == 0.

      ff_state.pipelineDepthStencilStateCreateInfo_get()
        .setFront(s_op_apply_fb).setBack(s_op_apply_fb)
        .setDepthCompareOp(vk::CompareOp::eEqual);

      // Allow pixels to be re-written.
      //
      vk::ColorComponentFlags cmask_all
        ( vk::ColorComponentFlagBits::eR | vk::ColorComponentFlagBits::eG |
          vk::ColorComponentFlagBits::eB | vk::ColorComponentFlagBits::eA );

      ff_state.pipelineColorBlendAttachmentState_get() =
        vk::PipelineColorBlendAttachmentState
        ( true,                   // blendEnable
          vk::BlendFactor::eOne,  // srcColorBlendFactor
          vk::BlendFactor::eOne,  // dstColorBlendFactor
          vk::BlendOp::eAdd,       // colorBlendOp
          vk::BlendFactor::eOne,   // srcAlphaBlendFactor
          vk::BlendFactor::eOne,  // dstAlphaBlendFactor
          vk::BlendOp::eAdd,       // alphaBlendOp
          cmask_all      // colorWriteMask
          );

      puni_light_curr = &uni_light;

      // Render.
      //
      render_objects(cb,RO_Normally);

      ff_state.pop();
    }

  puni_light_curr = &uni_light_all;
}


void
World::cb_keyboard()
{
  const int key = vh.keyboard_key_get();
  if ( !key ) return;
  pVect adjustment(0,0,0);
  pVect user_rot_axis(0,0,0);
  const bool kb_mod_s = vh.keyboard_shift;
  const bool kb_mod_c = vh.keyboard_control;
  const float move_amt = kb_mod_s ? 2.0 : kb_mod_c ? 0.08 : 0.4;

  switch ( key ) {
  case FB_KEY_LEFT: adjustment.x = -move_amt; break;
  case FB_KEY_RIGHT: adjustment.x = move_amt; break;
  case FB_KEY_PAGE_UP: adjustment.y = move_amt; break;
  case FB_KEY_PAGE_DOWN: adjustment.y = -move_amt; break;
  case FB_KEY_DOWN: adjustment.z = move_amt; break;
  case FB_KEY_UP: adjustment.z = -move_amt; break;
  case FB_KEY_DELETE: user_rot_axis.y = 1; break;
  case FB_KEY_INSERT: user_rot_axis.y =  -1; break;
  case FB_KEY_HOME:
    if ( opt_move_item == MI_Light ) adjustment.x = 1; else user_rot_axis.x = 1;
    break;
  case FB_KEY_END: user_rot_axis.x = -1; break;
  case '0': case '1': case '2': ball_setup_hw2( key ); break;
  case '3': ball_setup_3(); break;
  case '4': ball_setup_4(); break;
  case '5': ball_setup_5(); break;
  case '6': ball_setup_6(); break;
  case '7': ball_setup_7(); break;
  case 'a':
    opt_physics_method++;
    if ( opt_physics_method == GP_ENUM_SIZE ) opt_physics_method = GP_cpu;
    break;
  case 'b': opt_move_item = MI_Ball; break;
  case 'B': opt_move_item = MI_Ball_V; break;

  case 'c': opt_freeze_config = !opt_freeze_config; break;

  case 'e': case 'E': opt_move_item = MI_Eye; break;
  case 'g': case 'G': opt_gravity = !opt_gravity; break;
  case 'h': case 'H': opt_head_lock = !opt_head_lock; break;
#if 0
  case '!': opt_hide_stuff ^= OH_Sphere; break;
  case '@': opt_hide_stuff ^= OH_Links; break;
  case '#': opt_hide_stuff ^= OH_Platform; break;
#endif
  case 't': case 'T': opt_tail_lock = !opt_tail_lock; break;
  case 'l': case 'L': opt_move_item = MI_Light; break;
  //  case 'r': opt_mirror = !opt_mirror; break;
  case 'R':
    opt_want_raytrace = !opt_want_raytrace;
    world_time_link_update = 0;
    break;
#if 0
  case 'o': opt_shadows = !opt_shadows; break;
  case 'O':
    opt_shadow_volumes = !opt_shadow_volumes;
    if ( opt_shadow_volumes ) opt_shadows = true;
    break;
  case 'n': case 'N': opt_normal_sphere = !opt_normal_sphere; break;
#endif
  case 'p': case 'P': opt_pause = !opt_pause; break;
  case 's': case 'S': balls_stop(); break;
#if 0
  case 'v': case 'V': opt_shader++;
    if ( opt_shader >= sizeof(sh_names)/sizeof(sh_names[0]) ) opt_shader = 0;
    break;
  case 'w': case 'W': balls_twirl(); break;
#endif
  case 'y': opt_tryout1 = !opt_tryout1; break;
  case 'Y': opt_tryout2 = !opt_tryout2; break;
  //  case 'z': opt_sphere_true = !opt_sphere_true; break;
  case ' ':
    if ( kb_mod_s && kb_mod_c ) opt_single_time_step = 10;
    else if ( kb_mod_s ) opt_single_time_step = 100;
    else if ( kb_mod_c ) opt_single_time_step = 1; else opt_single_frame = true;
    opt_pause = true;
    break;
  case FB_KEY_TAB:
    if ( !kb_mod_s ) { variable_control.switch_var_right(); break; }
  case 96: variable_control.switch_var_left(); break;
  case '-':case '_': variable_control.adjust_lower(); break;
  case '+':case '=': variable_control.adjust_higher(); break;
  default: printf("Unknown key, %d\n",key); break;
  }

  gravity_accel.y = opt_gravity ? -opt_gravity_accel : 0;

  opt_mirror_depth = raytrace.recursion_max_set( opt_mirror_depth + 1 ) - 1;
  lock_update();
  sphere.bunch_invalidate();

  // Update eye_direction based on keyboard command.
  //
  if ( user_rot_axis.x || user_rot_axis.y )
    {
      pMatrix_Rotation rotall(eye_direction,pVect(0,0,-1));
      user_rot_axis *= invert(rotall);
      eye_direction *= pMatrix_Rotation(user_rot_axis, M_PI * 0.03);
      modelview_update();
    }

  // Update eye_location based on keyboard command.
  //
  if ( adjustment.x || adjustment.y || adjustment.z )
    {
      const double angle =
        fabs(eye_direction.y) > 0.99
        ? 0 : atan2(eye_direction.x,-eye_direction.z);
      pMatrix_Rotation rotall(pVect(0,1,0),-angle);
      adjustment *= rotall;

      switch ( opt_move_item ){
      case MI_Ball:
        adj_vector = adjustment;
        if ( adj_t_stop == 0 )
          adj_t_prev = adj_t_stop = world_time;
        adj_t_stop += adj_duration;
        break;
      case MI_Ball_V: balls_push(adjustment,0); break;
      case MI_Light:
        if ( key == FB_KEY_HOME )
          {
            pVect l_to_b(light_location,balls[0]->position);
            adjustment = 0.25 * l_to_b;
          }
        light_location += adjustment;
        break;
      case MI_Eye: eye_location += adjustment; break;
      default: break;
      }
      modelview_update();
    }

}




void
World::init(int argc, char **argv)
{
  chain_length = 14;

  const bool all_vars = false;

  opt_time_step_duration = 0.00003;
  if ( all_vars )
    variable_control.insert(opt_time_step_duration,"Time Step Duration");

  distance_relaxed = 15.0 / chain_length;
  opt_spring_constant = 500;
  variable_control.insert(opt_spring_constant,"Spring Constant");

  opt_spring_damp = 0.7;
  if ( all_vars )
    variable_control.insert(opt_spring_damp,"Spring Damping").set_max(1);

  opt_friction_dynamic = 0.2;
  if ( all_vars )
    variable_control.insert(opt_friction_dynamic,"Dynamic Friction");

  opt_omega = 10;
  if ( all_vars )
    variable_control.insert(opt_omega,"Initial Rotation Rate");

  opt_gravity_accel = 9.8;
  opt_gravity = true;
  gravity_accel = pVect(0,-opt_gravity_accel,0);
  variable_control.insert(opt_gravity_accel,"Gravity");

  opt_air_resistance = 0.01;
  variable_control.insert(opt_air_resistance,"Air Resistance");

  world_time = 0;
  time_step_count = 0;
  last_frame_wall_time = time_wall_fp();
  frame_timer.work_unit_set("Steps / s");

  world_time_link_update = 0;

  opt_mirror = true;
  opt_mirror_depth = 3;
  if ( all_vars )
    variable_control.insert(opt_mirror_depth,"Mirror Depth",1,0);
  opt_shadows = true;
  opt_shadow_volumes = false;
  opt_shader = SO_Plain;

  /// Init CUDA

  // Get information about GPU and its ability to run CUDA.
  //
  int device_count;
  cudaGetDeviceCount(&device_count); // Get number of GPUs.
  if ( device_count == 0 )
    {
      fprintf(stderr,"No GPU found, exiting.\n");
      exit(1);
    }

  // Print information about the available GPUs.
  //
  gpu_info_print();

  // Use course library routine to choose a good GPU for CUDA.
  //
  int dev = gpu_choose_index();
  CE(cudaSetDevice(dev));
  printf("Using GPU %d\n",dev);

  gpu_info.get_gpu_info(dev);
  cuda_setup(&gpu_info);

  // Print information about time_step routine.
  //
  printf("\nCUDA Routine Resource Usage:\n");

  int block_size_max = 4096;

  for ( int i=0; i<gpu_info.num_kernels; i++ )
    {
      printf("For %s:\n", gpu_info.ki[i].name);
      printf("  %6zd shared, %zd const, %zd loc, %d regs; "
             "%d max threads per block.\n",
             gpu_info.ki[i].cfa.sharedSizeBytes,
             gpu_info.ki[i].cfa.constSizeBytes,
             gpu_info.ki[i].cfa.localSizeBytes,
             gpu_info.ki[i].cfa.numRegs,
             gpu_info.ki[i].cfa.maxThreadsPerBlock);
      set_min(block_size_max,gpu_info.ki[i].cfa.maxThreadsPerBlock);
    }

  opt_physics_method = GP_cuda;
  gpu_const_data_stale = true;
  cpu_phys_data_stale = false;

  CE( cudaEventCreate(&frame_start_ce) );
  CE( cudaEventCreate(&frame_stop_ce) );

  opt_block_size = 32;

  if ( all_vars )
    variable_control.insert_power_of_2
      (opt_block_size, "Block Size", 32, block_size_max);

  c.alloc_n_balls = 0;
  c.alloc_n_links = 0;

  init_graphics();

  opt_sphere_slices = 35;
  sphere.opt_texture = true;
  sphere.init(opt_sphere_slices);
  if ( all_vars )
    variable_control.insert(opt_sphere_slices,"Sphere Slices");

  const char* const links_shader_code_path = "links-shdr-links.cc";

  const auto bvec_vs = vk::ShaderStageFlagBits::eVertex;
  const auto bvec_vgs =
    vk::ShaderStageFlagBits::eVertex | vk::ShaderStageFlagBits::eGeometry;

  pipe_set_2
    .init( ff_state )
    .ds_follow( transform )
    .ds_follow( puni_light_curr )
    .ds_uniform_use( "BIND_UNI_COMMON", buf_uni_common )
    .ds_set_material( color_dim_gray, color_salmon )
    .ds_use( sampler, texid_hw )
    .ds_storage_follow( "BIND_LINKS_POS1", lis_pos1, bvec_vs )
    .ds_storage_follow( "BIND_LINKS_POS2", lis_pos2, bvec_vs )
    .ds_storage_follow( "BIND_LINKS_V1", lis_v1, bvec_vs )
    .ds_storage_follow( "BIND_LINKS_V2", lis_v2, bvec_vs )
    .ds_storage_follow( "BIND_LINKS_B1_YDIR", lis_b1_ydir, bvec_vs )
    .ds_storage_follow( "BIND_LINKS_B2_YDIR", lis_b2_ydir, bvec_vs )
    .ds_storage_follow( "BIND_LINKS_TEX_RAD", lis_tex_rad, bvec_vgs )
    .shader_inputs_info_set( )
    .shader_code_set
    (links_shader_code_path,
     "vs_main_2(); ",     // Name of vertex shader main routine.
     "gs_main_2();",      // Name of geometry shader main routine.
     "fs_main();"         // Name of fragment shader main routine.
     )
    .topology_set( vk::PrimitiveTopology::eLineStrip )
    .create();

  if ( vh.have_raytrace )
    comp_links
      .init( vh.qs )
      .ds_uniform_use( "BIND_UNI_COMMON", buf_uni_common )
      .ds_uniform_use( "BIND_TRANSFORM", raytrace.buf_uni_transform )
      .ds_storage_follow( "BIND_LINKS_POS1", lis_pos1 )
      .ds_storage_follow( "BIND_LINKS_POS2", lis_pos2 )
      .ds_storage_follow( "BIND_LINKS_V1", lis_v1 )
      .ds_storage_follow( "BIND_LINKS_V2", lis_v2 )
      .ds_storage_follow( "BIND_LINKS_B1_YDIR", lis_b1_ydir )
      .ds_storage_follow( "BIND_LINKS_B2_YDIR", lis_b2_ydir )
      .ds_storage_follow( "BIND_LINKS_TEX_RAD", lis_tex_rad )
      .ds_storage_follow( "BIND_LINKS_AS_IDX", buf_links_as_idx )
      .ds_storage_follow( "BIND_LINKS_TEX_COOR", buf_links_tex_coor )
      .ds_storage_follow( "BIND_LINKS_VERTEX_G", buf_links_vertex_g )
      .ds_storage_follow( "BIND_LINKS_NORMAL_G", buf_links_normal_g )
      .shader_code_set(links_shader_code_path, "comp_main(); ")
      .create();

  pipe_set_2_sv
    .init( ff_state ).ds_use( uni_light )
    .ds_follow( transform )
    .ds_uniform_use( "BIND_UNI_COMMON", buf_uni_common )
    .ds_storage_follow( "BIND_LINKS_POS1", lis_pos1 )
    .ds_storage_follow( "BIND_LINKS_POS2", lis_pos2 )
    .ds_storage_follow( "BIND_LINKS_V1", lis_v1 )
    .ds_storage_follow( "BIND_LINKS_V2", lis_v2 )
    .ds_storage_follow( "BIND_LINKS_B1_YDIR", lis_b1_ydir )
    .ds_storage_follow( "BIND_LINKS_B2_YDIR", lis_b2_ydir )
    .ds_storage_follow( "BIND_LINKS_TEX_RAD", lis_tex_rad )
    .shader_inputs_info_set( )
    .shader_code_set
    (links_shader_code_path,
     "vs_main_2_sv(); ",     // Name of vertex shader main routine.
     "gs_main_2_sv();",      // Name of geometry shader main routine.
     "fs_main_sv();"         // Name of fragment shader main routine.
     )
    .topology_set( vk::PrimitiveTopology::eLineStrip )
    .create();

  ball_setup_hw2('1');
  return;
  // For benchmarking.
  opt_pause = true;
  opt_hide_stuff = OH_Links | OH_Platform;
  opt_shadows = false;
  opt_mirror = false;
}

Ball*
World::make_marker(pCoor position, pColor color, float radius )
{
  Ball* const ball = new Ball;
  ball->position = position;
  ball->locked = true;
  ball->velocity = pVect(0,0,0);
  ball->radius = radius;
  ball->mass = -1;
  ball->contact = false;
  ball->color = color;
  balls_marker.push_back(ball);
  return ball;
}

void
World::lock_update()
{
  // This routine called when options like opt_head_lock might have
  // changed.

  // Update locked status.
  //
  if ( head_ball ) head_ball->locked = opt_head_lock;
  if ( tail_ball ) tail_ball->locked = opt_tail_lock;

  // Re-compute minimum mass needed for stability.
  //
  for ( Ball *ball: balls ) ball->spring_constant_sum = 0;
  const double dtis = pow( opt_time_step_duration, 2 );
  for ( Link *link: links )
    {
      Ball* const b1 = link->ball1;
      Ball* const b2 = link->ball2;
      b1->spring_constant_sum += opt_spring_constant;
      b2->spring_constant_sum += opt_spring_constant;
    }
  for ( Ball *ball: balls )
    {
      ball->mass_min = ball->spring_constant_sum * dtis;
      ball->constants_update();
    }
  gpu_const_data_stale = true;
}

void
World::balls_twirl()
{
  if ( !head_ball || !tail_ball ) return;

  pNorm axis(head_ball->position, tail_ball->position);
  const float omega = 0.1 * opt_omega;

  for ( Ball *ball: balls )
    {
      pVect b_to_top(ball->position,head_ball->position);
      const float dist_along_axis = dot(b_to_top,axis);
      const float lsq = b_to_top.mag_sq() - dist_along_axis * dist_along_axis;
      ball->omega += omega * axis;
      if ( lsq <= 1e-5 ) continue;
      const float dist_to_axis = sqrt(lsq);
      pNorm rot_dir = cross(b_to_top,axis);
      ball->velocity += omega * dist_to_axis * rot_dir;
    }
}

void
World::objects_erase()
{
  balls.erase();
  links.erase();
  surfaces.clear();
}

Links
World::link_new
(Ball *ball1, Ball *ball2, float stiffness, float spring_constant_p)
{
  Link_Args largs
    { .spring_constant = spring_constant_p, .link_stiffness = stiffness };
  return link_new(ball1,ball2,largs);
}

Links
World::link_new
(Ball *ball1, Ball *ball2, Link_Args largs )
{
  const float spring_constant_p = largs.spring_constant;
  const float stiffness = largs.link_stiffness;
  Links links_rv;
  assert( ball1->radius > 0 );
  assert( ball2->radius > 0 );
  const float spring_constant = spring_constant_p ?: opt_spring_constant;

  Link* const rlink = new Link(ball1,ball2);
  rlink->radius = largs.radius;
  links_rv += rlink;

  pNorm n12(ball1->position,ball2->position);
  const float rad_sum = ball1->radius + ball2->radius;
  pMatrix3x3 b1rot = ball1->omatrix;
  pMatrix3x3 b2rot = ball2->omatrix;
  pCoor ctr = ball1->position
    + ( ball1->radius + 0.5 * ( n12.magnitude - rad_sum ) ) * n12;

  pVect con_x = rlink->con_x;
  pVect con_y = rlink->con_y;

  const float lrad = stiffness * ball1->radius;

  const int n_springs = 4;
  const double delta_theta = 2 * M_PI / n_springs;

  for ( int i=0; i<n_springs; i++ )
    {
      const double theta = i * delta_theta;
      pVect convec = lrad * ( cos(theta) * con_x + sin(theta) * con_y );
      pCoor con = ctr + convec;
      pVect cb1 = mult_MTV( b1rot, con - ball1->position );
      pVect cb2 = mult_MTV( b2rot, con - ball2->position );
      links_rv +=
        new Link(rlink, cb1, cb2,
                 { .distance_relaxed=0, .spring_constant=spring_constant } );
    }
  return links_rv;
}

///
/// CUDA Code
///

float4 tof4(pCoor c)
{
  float4 f4;
  f4.x = c.x;
  f4.y = c.y;
  f4.z = c.z;
  f4.w = c.w;
  return f4;
}

float4 tof4(pVect c)
{
  float4 f4;
  f4.x = c.x;
  f4.y = c.y;
  f4.z = c.z;
  f4.w = 0;
  return f4;
}

pCoor topc(float4 c)
{
  pCoor f4;
  f4.x = c.x;
  f4.y = c.y;
  f4.z = c.z;
  f4.w = c.w;
  return f4;
}

void
World::data_cpu_to_gpu_constants()
{
  c.platform_xmin = platform_xmin;
  c.platform_xmax = platform_xmax;
  c.platform_zmax = platform_zmax;
  c.platform_zmin = platform_zmin;

  c.n_balls = balls.size();
  c.n_links = links.size();
  c.n_surfaces = surfaces.size();
  c.opt_tryout1 = opt_tryout1;
  c.opt_tryout2 = opt_tryout2;
  c.opt_tryoutf = opt_tryoutf;
  c.opt_spring_constant = opt_spring_constant;
  c.opt_spring_damp = opt_spring_damp;
  c.opt_air_resistance = opt_air_resistance;
  c.opt_friction_dynamic = opt_friction_dynamic;
  c.opt_head_lock = opt_head_lock;
  c.opt_tail_lock = opt_tail_lock;
  c.gravity_accel = gravity_accel;

  data_cpu_to_gpu_common(&c);
}


void
World::data_cpu_to_gpu_dynamic()
{
  const int n_balls = balls.size();
  const int n_links = links.size();
  const int n_surfaces = surfaces.size();
  const bool need_alloc =
    n_balls > c.alloc_n_balls || n_links > c.alloc_n_links
    || n_surfaces > c.alloc_n_surfaces;
  const bool need_init = c.alloc_n_balls == 0;

#define ALLOC_MOVE_ZERO(n_balls,balls,memb,do_move,do_zero)                   \
  {                                                                           \
    const int size_elt_bytes = sizeof(c.balls.memb[0]);                       \
    const int size_bytes = n_balls * size_elt_bytes;                          \
    if ( need_alloc )                                                         \
      {                                                                       \
        if ( !need_init )                                                     \
          {                                                                   \
            CE( cudaFree( c.balls.memb ) );                                   \
            free(c.h_##balls.memb);                                           \
          }                                                                   \
        CE( cudaMalloc( &c.balls.memb, size_bytes ) );                        \
        c.h_##balls.memb = (decltype(c.balls.memb)) malloc( size_bytes );     \
      }                                                                       \
    if ( do_move )                                                            \
      {                                                                       \
        for ( int i=0; i<n_balls; i++ )                                       \
          c.h_##balls.memb[i] = balls[i]->memb;                               \
        CE( cudaMemcpy                                                        \
            ( c.balls.memb, c.h_##balls.memb,                                 \
              size_bytes, cudaMemcpyDefault ) );                              \
      }                                                                       \
    else if ( do_zero )                                                       \
      {                                                                       \
        CE( cudaMemset( c.balls.memb, 0, size_bytes ) );                      \
      }                                                                       \
  }

#define MOVEc(num,ptr,memb) \
  ALLOC_MOVE_ZERO(num,ptr,memb,gpu_const_data_stale,0)
#define MOVE(num,ptr,memb) ALLOC_MOVE_ZERO(num,ptr,memb,1,0)
#define ALLOC(num,ptr,memb) ALLOC_MOVE_ZERO(num,ptr,memb,0,0)
#define ZERO(num,ptr,memb) ALLOC_MOVE_ZERO(num,ptr,memb,0,1)
  MOVEc(n_balls,balls,position);
  MOVEc(n_balls,balls,velocity);
  ZERO(n_balls,balls,force);

  MOVEc(n_balls,balls,orientation);
  MOVEc(n_balls,balls,omatrix);
  MOVEc(n_balls,balls,omega);
  ZERO(n_balls,balls,torque);
  MOVEc(n_balls,balls,m_r_fdt);
  MOVEc(n_balls,balls,locked);

  MOVEc(n_links,links,ball_idx);
  MOVEc(n_links,links,cb1);
  MOVEc(n_links,links,cb2);
  MOVEc(n_links,links,dr_sc);
  MOVEc(n_links,links,is_simulatable);
  MOVEc(n_links,links,is_surface_connection);

  MOVEc(n_surfaces,surfaces,b_idx);

#undef MOVE
#undef MOVEc
#undef ALLOC
#undef ZERO
#undef ALLOC_MOVE
  if ( need_alloc )
    {
      c.alloc_n_balls = n_balls;
      c.alloc_n_links = n_links;
      c.alloc_n_surfaces = n_surfaces;
      data_cpu_to_gpu_common(&c);
    }
  gpu_const_data_stale = false;
}

void
World::data_gpu_to_cpu_dynamic()
{
  const int n_balls = balls.size();

#define MOVE(n_balls,balls,memb)                                              \
  {                                                                           \
    const int size_elt_fr_bytes = sizeof(c.balls.memb[0]);                    \
    const int size_bytes = n_balls * size_elt_fr_bytes;                       \
    CE( cudaMemcpy                                                            \
        ( c.h_##balls.memb, c.balls.memb, size_bytes,                         \
          cudaMemcpyDeviceToHost ) );                                         \
    for ( int i=0; i<n_balls; i++ )                                           \
      balls[i]->memb = c.h_##balls.memb[i];                                   \
  }

  MOVE(n_balls,balls,position);
  MOVE(n_balls,balls,velocity);
  MOVE(n_balls,balls,orientation);
  MOVE(n_balls,balls,omatrix);
  MOVE(n_balls,balls,omega);

#undef MOVE

}


///
/// Physical Simulation Code
///

 /// Initialize Simulation
//

void
World::ball_setup_prob_0()
{
  /// Homework 2
  //
  //  Use this routine for experimentation.
  //  Put the solution to Problem 1 in World::ball_setup_prob_1() and
  //  put the solution to Problem 2 in World::ball_setup_prob_2().

  // The axis used to construct the hyperbolid.
  //
  const pNorm axis{h.axis};

  // A point on the axis used to construct the hyperboloid.
  //
  const pCoor center_pos{h.center_pos};
  //
  // This will also be the center of the circle that forms the base
  // of the hyperboloid.

  // The endpoints of the first line used to construct the hyperboloid.
  //
  const pCoor line_0_bot_pos{h.line_0_bot_pos};
  const pCoor line_0_top_pos{h.line_0_top_pos};
  //
  // The balls along this line are yellow (in the unmodified code).

  // The number of balls in one "line".
  //
  const int n_balls = h.n_balls;
  //
  // The yellow balls form one line.

  // The number of lines.
  //
  const int n_lines = h.n_lines;

  const float ball_radius{h.ball_radius};


  // Drop a marker ball at the center position.
  //
  make_marker( center_pos, color_green, ball_radius*5 );
  //
  // Feel free to insert markers at additional points to help debug.

  // Settings for the link.  No need to modify this.
  Link_Args largs{
    .spring_constant = opt_spring_constant / 10,
    .link_stiffness = 50,
    .radius = ball_radius * 0.7f
  };

  // Container to hold pointers to balls. Organized to facilitate putting
  // links between balls.
  //
  vector<vector<Ball*>> balls_lines;

  // Compute a vector that will advance from one vector to another.
  //
  pVect grid_pitch( 5 * ball_radius, 0, 0 );

  for ( int j=0; j<n_lines; j++ )
    {
      // Construct a container to hold the current line.
      //
      vector<Ball*>& balls_line = balls_lines.emplace_back();

      // Compute the endpoints of the current line to draw.
      // Needs to be modified for Problem 1 and maybe Problem 2.
      //
      pCoor line_bot_pos = line_0_bot_pos + j * grid_pitch;
      pCoor line_top_pos = line_0_top_pos + j * grid_pitch;
      //
      // In the first iteration these should match line_0_top_pos,
      // line_0_bot_pos.

      // Note: For Problem 1, there is no need to modify the code below,
      // in particular in the "i" loop and the code that follows.

      // Compute a vector that advances along a line, from bottom to top.
      //
      pVect line_pitch = (1.0/(n_balls-1)) * ( line_top_pos - line_bot_pos );

      // Construct a line consisting of n_balls balls.
      //
      for ( int i=0; i<n_balls; i++ )
        {
          // Compute the position along the line between line_bot_pos
          // and line_top_pos.
          //
          pCoor pos = line_bot_pos + i * line_pitch;

          // Create a new ball, and add it to "our" list, balls_line,
          // and the simulators list, balls.
          //
          Ball* const ball = new Ball(pos);
          balls_line.push_back( ball );   // Our list of balls, for linking.
          balls.push_back( ball );        // List of simulated balls.
          //
          // Since ball is a pointer it is still possible to set
          // members of ball, such as ball->color, after ball was
          // added to the lists.

          ball->velocity = pVect(0,0,0);
          ball->radius = ball_radius;
          ball->color =
            j == 0 ? color_gold :
            i == 0 ? color_pale_green :
            i == n_balls - 1 ? color_blue :
            color_light_slate_gray;
        }
    }

  // Put links between balls in adjacent lines.
  //
  for ( int j=0; j<balls_lines.size(); j++ )
    {
      auto &bl0 = balls_lines[j], &bl1 = balls_lines[(j+1)%n_lines];
      for ( int i = 0; i<bl0.size(); i++ )
        {
          // Add a link along line j. (From bottom to top.)
          //
          if ( i ) links += link_new( bl0[i-1], bl0[i], largs );

          // Add a link from line j to line j+1.
          //
          links += link_new( bl0[i], bl1[i], largs );
        }
    }

  // Choose a head and tail balls. These can be locked in place.
  //
  tail_ball = balls_lines[0].front();
  head_ball = balls_lines[0].back();
  opt_head_lock = false;  opt_tail_lock = false;
}

void
World::ball_setup_prob_1()
{
  /// Homework 2
  //
  //  Put solution to Problem 1 in this routine.

  // The axis used to construct the hyperbolid.
  //
  const pNorm axis{h.axis};

  // A point on the axis used to construct the hyperboloid.
  //
  const pCoor center_pos{h.center_pos};
  //
  // This will also be the center of the circle that forms the base
  // of the hyperboloid.

  // The endpoints of the first line used to construct the hyperboloid.
  //
  const pCoor line_0_bot_pos{h.line_0_bot_pos};
  const pCoor line_0_top_pos{h.line_0_top_pos};
  //
  // The balls along this line are yellow (in the unmodified code).

  // The number of balls in one "line".
  //
  const int n_balls = h.n_balls;
  //
  // The yellow balls form one line.

  // The number of lines.
  //
  const int n_lines = h.n_lines;

  const float ball_radius{h.ball_radius};


  // Drop a marker ball at the center position.
  //
  make_marker( center_pos, color_green, ball_radius*5 );
  //
  // Feel free to insert markers at additional points to help debug.

  // Settings for the link.  No need to modify this.
  Link_Args largs{
    .spring_constant = opt_spring_constant / 10,
    .link_stiffness = 50,
    .radius = ball_radius * 0.7f
  };

  // Container to hold pointers to balls. Organized to facilitate putting
  // links between balls.
  //
  vector<vector<Ball*>> balls_lines;

  /// SOLUTION -- Problem 1
  //
  //  Solution Plan:
  //
  //    - Find axes of bottom circle.
  //    - Find center of top circle.
  //    - Find axes of top circle. (Instead could find twist angle.)
  //
  /// Find Axes for the Bottom Circle
  //
  //  The x axis is directly constructed using center_pos and line_0_top_pos.
  pVect vx_bot( center_pos, line_0_bot_pos );
  //
  //  The y axis is found using the just-computed vx_bot and the given axis.
  pVect vy_bot = cross(axis,vx_bot);
  //
  /// Find Center of Top Circle
  //
  //  Compute the coordinates of the center of the top circle.
  //  The top circle is in the plane defined by point line_0_top_pos
  //  and the normal in variable axis.
  //
  //  Line from bottom plane (and circle) to top plane (and circle).
  pVect v0tb( line_0_bot_pos, line_0_top_pos );
  //
  //  The dot product gives the distance from one plane to the other.
  pVect vht = dot(v0tb,axis) * axis;
  //
  //  Use the distance (vht) to find the center of the top circle.
  pCoor ctr_top_pos = center_pos + vht;
  //
  /// Find Axes of Top Circle
  //
  //  The bottom circle and the top circle are in parallel planes and
  //  so vx_bot and vy_bot could be used for both. If vx_bot,vy_bot
  //  were used for the upper circle we'd need to compute a "twist"
  //  angle. The problem with using a twist angle is that we'd need to
  //  compute a sine and cosine twice per iteration, once for theta
  //  and once for theta+twist_angle. Rather than do that a second set
  //  of axes will be computed so that the same angle could be used
  //  for both.
  //
  pVect vx_top( ctr_top_pos, line_0_top_pos );
  pVect vy_top = cross(axis,vx_top);
  //
  //  Place a marker at the center of the top circle.
  make_marker( ctr_top_pos, color_blue, ball_radius * 5 );
  // 
  const float delta_theta = 2 * M_PI / n_lines;

  // Compute a vector that will advance from one vector to another.
  //
  pVect grid_pitch( 5 * ball_radius, 0, 0 );

  for ( int j=0; j<n_lines; j++ )
    {
      // Construct a container to hold the current line.
      //
      vector<Ball*>& balls_line = balls_lines.emplace_back();

      // Compute the endpoints of the current line to draw.
      // Needs to be modified for Problem 1 and maybe Problem 2.
      //
#if 0
      // Original code (before solution)
      pCoor line_bot_pos = line_0_bot_pos + j * grid_pitch;
      pCoor line_top_pos = line_0_top_pos + j * grid_pitch;
#endif
      //
      // In the first iteration these should match line_0_top_pos,
      // line_0_bot_pos.

      /// SOLUTION -- Problem 1
      //
      //  Use the upper and lower circle axes to find the top and bottom
      //  endpoints of the line rotated by theta.
      //
      const float theta = j * delta_theta;
      pCoor line_bot_pos =
        center_pos + vx_bot * cosf(theta) + vy_bot * sinf(theta);
      pCoor line_top_pos =
        ctr_top_pos + vx_top * cosf(theta) + vy_top * sinf(theta);

      // Note: For Problem 1, there is no need to modify the code below,
      // in particular in the "i" loop and the code that follows.

      // Compute a vector that advances along a line, from bottom to top.
      //
      pVect line_pitch = (1.0/(n_balls-1)) * ( line_top_pos - line_bot_pos );

      // Construct a line consisting of n_balls balls.
      //
      for ( int i=0; i<n_balls; i++ )
        {
          // Compute the position along the line between line_bot_pos
          // and line_top_pos.
          //
          pCoor pos = line_bot_pos + i * line_pitch;

          // Create a new ball, and add it to "our" list, balls_line,
          // and the simulators list, balls.
          //
          Ball* const ball = new Ball(pos);
          balls_line.push_back( ball );   // Our list of balls, for linking.
          balls.push_back( ball );        // List of simulated balls.
          //
          // Since ball is a pointer it is still possible to set
          // members of ball, such as ball->color, after ball was
          // added to the lists.

          ball->velocity = pVect(0,0,0);
          ball->radius = ball_radius;
          ball->color =
            j == 0 ? color_gold :
            i == 0 ? color_pale_green :
            i == n_balls - 1 ? color_blue :
            color_light_slate_gray;
        }
    }

  // Put links between balls in adjacent lines.
  //
  for ( int j=0; j<balls_lines.size(); j++ )
    {
      auto &bl0 = balls_lines[j], &bl1 = balls_lines[(j+1)%n_lines];
      for ( int i = 0; i<bl0.size(); i++ )
        {
          // Add a link along line j. (From bottom to top.)
          //
          if ( i ) links += link_new( bl0[i-1], bl0[i], largs );

          // Add a link from line j to line j+1.
          //
          links += link_new( bl0[i], bl1[i], largs );
        }
    }

  // Choose a head and tail balls. These can be locked in place.
  //
  tail_ball = balls_lines[0].front();
  head_ball = balls_lines[0].back();
  opt_head_lock = false;  opt_tail_lock = false;
}

void
World::ball_setup_prob_2()
{
  /// Homework 2
  //
  //  Put solution to Problem 2 in this routine.

  // The axis used to construct the hyperbolid.
  //
  const pNorm axis{h.axis};

  // A point on the axis used to construct the hyperboloid.
  //
  const pCoor center_pos{h.center_pos};
  //
  // This will also be the center of the circle that forms the base
  // of the hyperboloid.

  // The endpoints of the first line used to construct the hyperboloid.
  //
  const pCoor line_0_bot_pos{h.line_0_bot_pos};
  const pCoor line_0_top_pos{h.line_0_top_pos};
  //
  // The balls along this line are yellow (in the unmodified code).

  // The number of balls in one "line".
  //
  const int n_balls = h.n_balls;
  //
  // The yellow balls form one line.

  // The number of lines.
  //
  const int n_lines = h.n_lines;

  const float ball_radius{h.ball_radius};


  // Drop a marker ball at the center position.
  //
  make_marker( center_pos, color_green, ball_radius*5 );
  //
  // Feel free to insert markers at additional points to help debug.

  // Settings for the link.  No need to modify this.
  Link_Args largs{
    .spring_constant = opt_spring_constant / 10,
    .link_stiffness = 50,
    .radius = ball_radius * 0.7f
  };

  // Container to hold pointers to balls. Organized to facilitate putting
  // links between balls.
  //
  vector<vector<Ball*>> balls_lines;

  // Compute a vector that will advance from one vector to another.
  //
  pVect grid_pitch( 5 * ball_radius, 0, 0 );

  // Compute a vector that advances along a line, from bottom to top.
  //
  pVect line_pitch = (1.0/(n_balls-1)) * ( line_0_top_pos - line_0_bot_pos );

  /// SOLUTION -- Problem 2
  //
  //  Define out the sample transformations.
  #if 0
  // Translate so that coordinate P is moved to P + grid_pitch.
  pMatrix_Translate transl( grid_pitch );

  // The following are examples of transformation matrices that may or
  // may not be useful.
  //
  // Rotate pi/2 radians around the x axis.
  pMatrix_Rotation rot( pVect(1,0,0), M_PI / 2 );
  //
  pNorm ax(1,2,3), ay(4,5,6), az(7,8,9);
  // Rotate so that direction ax becomes pVect(1,0,0) (x-axis), etc.
  pMatrix_Rows rot_to_local( ax, ay, az );
  //
  // Rotate so that direction pVect(1,0,0) becomes ax.
  pMatrix_Cols rot_from_local( ax, ay, az );
  //
  // Examples of pMatrix constructors shown above. 
  #endif

  /// SOLUTION -- Problem 2
  //
  // Translate coordinate space so center_pos is origin.
  //
  pMatrix_Translate to_ctr( -center_pos );
  //
  // Rotate around axis.
  //
  const float delta_theta = 2 * M_PI / n_lines;
  pMatrix_Rotation rot( axis, delta_theta );
  //
  // Translate back so that center_pos is where it was.
  //
  pMatrix_Translate from_ctr( center_pos );
  //
  // Combine the three to get one transformation matrix.
  //
  pMatrix transform = from_ctr * rot * to_ctr;

  for ( int j=0; j<n_lines; j++ )
    {
      // Construct a container to hold the current line.
      //
      vector<Ball*>& balls_line = balls_lines.emplace_back();

      // Construct a line consisting of n_balls balls.
      //
      for ( int i=0; i<n_balls; i++ )
        {
          pCoor pos;

          // Compute the position along the line between line_bot_pos
          // and line_top_pos.
          //
          if ( j == 0 )
            {
              // Use the line_0 end point to find the position on the
              // first line.
              //
              pos = line_0_bot_pos + i * line_pitch;
            }
          else
            {
              // Use find the position on ball i on this line using a
              // transformation of the position of ball i one the
              // previous line.
              //
              pCoor pos_prev = balls_lines[j-1][i]->position;
              pos = transform * pos_prev;
            }

          // Create a new ball, and add it to "our" list, balls_line,
          // and the simulators list, balls.
          //
          Ball* const ball = new Ball(pos);
          balls_line.push_back( ball );   // Our list of balls, for linking.
          balls.push_back( ball );        // List of simulated balls.
          //
          // Since ball is a pointer it is still possible to set
          // members of ball, such as ball->color, after ball was
          // added to the lists.

          ball->velocity = pVect(0,0,0);
          ball->radius = ball_radius;
          ball->color =
            j == 0 ? color_gold :
            i == 0 ? color_pale_green :
            i == n_balls - 1 ? color_blue :
            color_light_slate_gray;
        }
    }

  // Put links between balls in adjacent lines.
  //
  for ( int j=0; j<balls_lines.size(); j++ )
    {
      auto &bl0 = balls_lines[j], &bl1 = balls_lines[(j+1)%n_lines];
      for ( int i = 0; i<bl0.size(); i++ )
        {
          // Add a link along line j. (From bottom to top.)
          //
          if ( i ) links += link_new( bl0[i-1], bl0[i], largs );

          // Add a link from line j to line j+1.
          //
          links += link_new( bl0[i], bl1[i], largs );
        }
    }

  // Choose a head and tail balls. These can be locked in place.
  //
  tail_ball = balls_lines[0].front();
  head_ball = balls_lines[0].back();
  opt_head_lock = false;  opt_tail_lock = false;
}



void
World::ball_setup_1(bool big)
{
  // Arrange balls to form sort of tree.

  last_setup = 1;

  double ball_rad = 1;  // Warning: Re-assigned, changes behavior of nb.

  pCoor first_pos(17.2,ball_rad,-20.2);

  pVect dir_up(0.1,3,0.1);
  pVect dir_z(3*distance_relaxed,0,0);
  pVect dir_x(0,0,3*distance_relaxed);

  // Remove objects from the simulated objects lists, balls and links.
  // The delete operator is used on objects in the lists.
  //
  objects_erase();

  auto nb =
    [&] (pCoor pos) {
      Ball* const ball = new Ball;
      ball->position = pos;
      ball->locked = false;
      ball->velocity = pVect(0,0,0);
      ball->radius = ball_rad;
      ball->contact = false;
      balls += ball;
      return ball;
      };

  tail_ball = nb(first_pos);
  tail_ball->specularity = 0.5;

  pCoor last_pos = tail_ball->position;

  const int n_tiers = big ? 32 : 8;
  const double shrink_final = 0.17;
  const double shrink = pow(shrink_final,1.0/n_tiers);

  Balls tail_balls;
  for ( int i=1; i<n_tiers; i++ )
    {
      ball_rad *= shrink;
      Ball* const ball = nb(last_pos + (i==1?1.0:ball_rad) * dir_up);
      ball->specularity = 0.5;
      tail_balls += ball;
      links += link_new(balls[balls-2],ball,0.8);
      last_pos = ball->position;
    }

  head_ball = balls[balls-1];

  pNorm dir_zn(dir_z);
  pNorm dir_xn(dir_x);
  for ( Ball *ball: tail_balls )
    {
      if ( ball == head_ball ) break;
      pNorm to_top(ball->position,head_ball->position);
      const float r = to_top.magnitude;
      int num_hairs = max(4,int(7*r));
      ball_rad = ball->radius / 4;
      const float delta_theta = 2 * M_PI / num_hairs;
      for ( int i=0; i<num_hairs; i++ )
        {
          const float theta = i * delta_theta;
          pCoor pos = ball->position +
            r * cos(theta) * dir_xn + r * sin(theta) * dir_zn;
          Ball* const ballr = nb(pos);
          links += link_new(ballr,ball);
        }
    }

  opt_head_lock = false;
  opt_tail_lock = true;
  setup_at_end();
}

void
World::ball_setup_2()
{
  // Arrange balls to form a suspended thing.

  last_setup = 2;

  pCoor first_pos(17.2,14.8f,-20.2);

  pVect dir_dn(0.5,first_pos.y/8,0.5);
  pVect dir_z(3*distance_relaxed,0,0);
  pVect dir_x(0,0,3*distance_relaxed);

  // Remove objects from the simulated objects lists, balls and links.
  // The delete operator is used on objects in the lists.
  //
  objects_erase();

  float ball_rad = 0.3;  // Warning: Re-assigned, changes behavior of nb.

  auto nb =
    [&] (pCoor pos) {
      Ball* const ball = new Ball;
      ball->position = pos;
      ball->locked = false;
      ball->velocity = pVect(0,0,0);
      ball->radius = ball_rad;
      ball->contact = false;
      ball->specularity = 0.5;
      balls += ball;
      return ball;
      };

  head_ball = nb(first_pos);
  for ( int j: {-1,1} )
    for ( int i: {-1,0,1} )
      {
        Ball* const ball = nb(first_pos - dir_dn + i * dir_x + j * dir_z );
        if ( i != 0 ) links += link_new(head_ball,ball);
        if ( i != -1 ) links += link_new(balls[balls-2],ball);
        if ( j == 1 && i != 0 )  links += link_new(balls[balls-4],ball);
      }

  Ball* const tail_start = nb(first_pos - dir_dn);
  links += link_new(balls[balls-6],tail_start);
  links += link_new(balls[balls-3],tail_start);
  pCoor last_pos = tail_start->position;
  Balls tail_balls;
  for ( int i=1; i<5; i++ )
    {
      Ball* const ball = nb(last_pos - 0.5 * dir_dn);
      tail_balls += ball;
      links += link_new(balls[balls-2],ball);
      last_pos = ball->position;
    }

  tail_ball = balls[balls-1];

  pNorm dir_zn(dir_z);
  pNorm dir_xn(dir_x);
  int num_hairs = 20;
  ball_rad = 0.1;
  for ( Ball *ball: tail_balls )
    {
      const float delta_theta = 2 * M_PI / num_hairs;
      for ( int i=0; i<num_hairs; i++ )
        {
          const float theta = i * delta_theta;
          const float r = ball->radius * 10;
          pCoor pos = ball->position +
            r * cos(theta) * dir_xn + r * sin(theta) * dir_zn;
          Ball* const ballr = nb(pos);
          links += link_new(ballr,ball);
        }
    }

  opt_head_lock = true;    // Head ball will be frozen in space.
  opt_tail_lock = false;   // Tail ball can move freely.
  setup_at_end();
}

void
World::ball_setup_3()
{
  // Arrange links to form what I would eventually learn to be a
  // corona virus.

  const bool sarcov = true;

  last_setup = 3;
  objects_erase();

  pCoor first_pos(17.2,14.8f,-10.2);

  float ball_radius = 1;

  auto nb =
    [&] (pCoor pos) {
      Ball* const ball = new Ball;
      ball->position = pos;
      ball->locked = false;
      ball->velocity = pVect(0,0,0);
      ball->radius = ball_radius;
      ball->contact = false;
      balls += ball;
      return ball;
      };

  head_ball = nb(first_pos);
  head_ball->specularity = 0.5;
  ball_radius = ball_radius * ( sarcov ? 0.08 : .1 );
  const int n_spikes = sarcov ? 100 : 500;

  for ( int i=0; i<n_spikes; i++ )
    {
      pNorm dir(0.5-drand48(),0.5-drand48(),0.5-drand48());
      const float len = head_ball->radius * ( sarcov ? 1.5 : 3.0 );
      Ball* const ball = nb(first_pos + len * dir);
      links += link_new(ball,head_ball,2);
    }

  tail_ball = balls[1];

  opt_head_lock = false;    // Head ball will be frozen in space.
  opt_tail_lock = false;   // Tail ball can move freely.

  setup_at_end();
}

void
World::ball_setup_4()
{
  // Render a wheel (disk of balls with one in the center).

  last_setup = 4;
  objects_erase();

  pCoor center_pos = pCoor(17.2,14.8f,-10.2) + 2*vec_rand();

  const int n_balls = 4 * ( 4 + random()%4 );
  float ring_outer_radius = 5 + 2 * drand48();
  float ball_radius = 0.2 * 2 * M_PI * ring_outer_radius / n_balls;

  pNorm c_to_eye(center_pos,eye_location);
  pNorm az = c_to_eye + 0.5 * vec_rand() + pVect(0,0.5,0);

  pNorm ax = fabs(az.x) < fabs(az.y)
    ? pVect(0,az.z,-az.y) : pVect(az.z,0,-az.x);
  pVect ay(az,ax);

  // Default settings for balls.
  //
  Ball ball_default;
  ball_default.specularity = 0.5;
  ball_default.velocity = pVect(0,0,0);
  ball_default.radius = ball_radius;
  ball_default.color = color_light_slate_gray;

  // Scale axes by radius so that they can be used to find points on circle.
  //
  pVect vx = ring_outer_radius * ax;
  pVect vy = ring_outer_radius * ay;

  const float d_theta = 2 * M_PI / n_balls;

  // Place balls in a circle.
  //
  for ( int i=0; i<n_balls; i++ )
    {
      const float theta = d_theta * i;
      pCoor pos = center_pos + vx * cosf(theta) + vy * sinf(theta);

      // Construct new ball and add it to the ball list.
      //
      Ball* ball = balls.emplace_back( new Ball(pos,ball_default) );

      // Set color of first ball to blue, others to white.
      ball->color = i == 0 ? color_blue : color_white;
    }

  const float link_stiffness = 10;

  // Connect each ball to its neighbor with a link.
  //
  for ( int i=0; i<n_balls; i++ )
    links += link_new( balls[i], balls[(i+1)%n_balls], link_stiffness );

  // Place a ball in the center and connect it to all the other balls.
  if ( true )
    {
      Ball* const center_ball =
        balls.emplace_back( new Ball(center_pos,ball_default) );

      for ( int i=0; i<n_balls; i++ )
        links += link_new( center_ball, balls[i], 2 );
    }

  pNorm spin_axis = az;
  pCoor c_of_mass = center_pos;

  for ( Ball* ball: balls )
    {
      pVect cm_to_ball( c_of_mass, ball->position );

      ball->velocity = opt_omega * cross( spin_axis, cm_to_ball );
      ball->omega = opt_omega * spin_axis;
    }

  tail_ball = balls.back();
  head_ball = balls.front();

  opt_head_lock = false; // If true, head ball locked in place.
  opt_tail_lock = false;

  setup_at_end();
}

void
World::ball_setup_5()
{
  // Construct a top. (Top, as in a child's toy.)

  last_setup = 5;
  objects_erase();

  pCoor center_pos(17.2,14.8f,-10.2);

  const int n_balls = 20;
  const float crad = 5;

  const float ball_radius = 0.3 * 2 * M_PI * crad / n_balls;

  // Default settings for balls.
  //
  Ball ball_default;
  ball_default.specularity = 0.5;
  ball_default.velocity = pVect(0,0,0);
  ball_default.radius = ball_radius;
  ball_default.color = pColor(0xffbbff); // Plum 1. That's the color's name.

  const float d_theta = 2 * M_PI / n_balls;

  pVect vz = pVect(1,6,1) - 2*vec_rand();
  pNorm ax( vz.y, -vz.x, 0 );
  pNorm ay = cross(vz,ax);
  pVect vx = crad * ax;
  pVect vy = crad * ay;

  pNorm spin_axis = vz;
  float omega = opt_omega;
  opt_spring_damp = .4;

  // Construct Top's Main Ring
  //
  for ( int i=0; i<n_balls; i++ )
    {
      const float theta = d_theta * i;
      pCoor pos = center_pos + vx * cosf(theta) + vy * sinf(theta);
      Ball* ball = balls.emplace_back( new Ball(pos,ball_default) );

      // Find the coordinate on the spin axis closest to the new ball.
      //
      pCoor axis_cp =
        center_pos + spin_axis * dot( spin_axis, pVect(center_pos,pos) );
      //
      // Use this coordinate to compute a vector giving the ball's
      // motion when the top is spinning around axis spin_axis.
      //
      pVect cp_to_p(axis_cp,pos);
      ball->velocity = omega * cross(spin_axis,cp_to_p);
    }

  // Make one of the balls on the top's main ring a different color.
  balls[0]->color = color_blue;

  // Connect the balls on the main ring with links.
  //
  for ( int i=0; i<n_balls; i++ )
    links += link_new( balls[i], balls[(i+1)%n_balls], 2 );

  // Complete the top with balls above and below the ring, and one in
  // the center.
  //
  pNorm az(vz);
  pCoor top_pos = center_pos + az * crad * 1.5;
  pCoor bot_pos = center_pos - az * crad * 1.5;

  head_ball = balls.emplace_back( new Ball(top_pos,ball_default) );
  tail_ball = balls.emplace_back( new Ball(bot_pos,ball_default) );
  Ball* ctr_ball = balls.emplace_back( new Ball(center_pos,ball_default) );

  for ( int i=0; i<n_balls; i += 2 )
    {
      links += link_new( head_ball, balls[i], 2 );
      links += link_new( tail_ball, balls[i], 2 );
    }

  for ( int i=0; i<n_balls; i++ )
    links += link_new( ctr_ball, balls[i], 2 );

  for ( auto b: balls ) b->omega = omega * spin_axis;

  opt_head_lock = false;    // Head ball will be frozen in space.
  opt_tail_lock = false;   // Tail ball can move freely.

  setup_at_end();
}

void
World::ball_setup_6()
{
  // Render what will become a parachute, probably.

  last_setup = 6;
  objects_erase();

  pCoor center_pos = pCoor(17.2,14.8f,-10.2) + 2*vec_rand();

  const int n_balls = 4 * ( 4 + random()%4 );
  float ring_outer_radius = 5 + 2 * drand48();
  float ring_inner_radius = ring_outer_radius * ( 0.1 + drand48() * 0.4 );
  float ball_radius = 0.2 * 2 * M_PI * ring_outer_radius / n_balls;

  pNorm c_to_eye(center_pos,eye_location);
  pNorm az = c_to_eye + 0.5 * vec_rand() + pVect(0,0.5,0);

  pNorm ax = fabs(az.x) < fabs(az.y)
    ? pVect(0,az.z,-az.y) : pVect(az.z,0,-az.x);
  pVect ay(az,ax);

  // Default settings for balls.
  //
  Ball ballo_default;
  ballo_default.specularity = 0.5;
  ballo_default.velocity = pVect(0,0,0);
  ballo_default.radius = ball_radius;
  ballo_default.color = color_light_slate_gray;
  Ball balli_default = ballo_default;
  balli_default.radius = ballo_default.radius * 0.6;

  const float d_theta = 2 * M_PI / n_balls;

  const float link_stiffness = 10;

  vector<Ball*> balls_outer, balls_inner;

  // Place balls in a circle.
  //
  for ( int i=0; i<n_balls; i++ )
    {
      const float theta = d_theta * i;
      pVect pos_n = cosf(theta) * ax + sinf(theta) * ay;
      pCoor pos_o = center_pos + ring_outer_radius * pos_n;
      pCoor pos_i = center_pos + ring_inner_radius * pos_n;

      // Construct new ball and add it to the ball list.
      //
      Ball* ball_o = balls_outer.emplace_back( new Ball(pos_o,ballo_default) );
      Ball* ball_i = balls_inner.emplace_back( new Ball(pos_i,balli_default) );
      Link *link =
        new Link( ball_o, ball_i,
                  { .spring_constant = opt_spring_constant/200 } );
      links.push_back(link);

      // Set color of first ball to blue, others to white.
      ball_i->color = ball_o->color =
        i == 0 ? color_blue :
        theta <= M_PI/2 ? color_lsu_spirit_gold :
        color_white;
    }

  // Connect each ball to its neighbor with a link.
  //
  for ( auto bc: { &balls_outer, &balls_inner } )
    for ( int i=0; i<n_balls; i++ )
      {
        balls.push_back( bc[0][i] );
        links += link_new( bc[0][i], bc[0][(i+1)%n_balls], link_stiffness );
      }

  for ( int i=0; i<n_balls; i++ )
    {
      const int j = ( i + 1 ) % n_balls;
      Ball* b00 = balls_inner[i];
      Ball* b01 = balls_outer[i];
      Ball* b10 = balls_inner[j];
      Ball* b11 = balls_outer[j];
      surfaces.emplace_back(b00,b01,b11,color_light_gray);
      surfaces.emplace_back(b10,b00,b11,color_light_gray);
    }

  pNorm spin_axis = az;
  pCoor c_of_mass = center_pos;
  float omega_init = opt_omega * 0.1;

  if ( 0 )
  for ( Ball* ball: balls )
    {
      pVect cm_to_ball( c_of_mass, ball->position );

      ball->velocity = omega_init * cross( spin_axis, cm_to_ball );
      ball->omega = omega_init * spin_axis;
    }

  tail_ball = balls.back();
  head_ball = balls.front();

  opt_head_lock = false; // If true, head ball locked in place.
  opt_tail_lock = false;

  setup_at_end();
}

void
World::ball_setup_7()
{
  // Render what will become a parachute, probably.

  last_setup = 7;
  objects_erase();
  setup_at_end();
}

void
World::setup_at_end()
{
  balls.take(balls_marker);
  lock_update();
  time_step_count = 0;
  world_time = 0;
  world_time_link_update = -1;
  for ( uint i=0; i<balls.size(); i++ ) balls[i]->idx = i;
  for ( Link *l : links )
    { l->ball1_idx = l->ball1->idx; l->ball2_idx = l->ball2->idx; }
  for ( auto& s: surfaces )
    { s.b0_idx = s.b0->idx;  s.b1_idx = s.b1->idx;  s.b2_idx = s.b2->idx; }
}


 /// Advance Simulation State by delta_t Seconds
//
void
World::time_step_cpu(double delta_t)
{
  const int n_balls = balls;
  const int n_links = links;

  // Smoothly move ball in response to user input.
  //
  if ( adj_t_stop )
    {
      const double dt = min(world_time,adj_t_stop) - adj_t_prev;
      pVect adj = dt/adj_duration * adj_vector;
      balls_translate(adj,0);
      adj_t_prev = world_time;
      if ( world_time >= adj_t_stop ) adj_t_stop = 0;
    }

#pragma omp parallel for
  for ( int bi=0; bi<n_balls; bi++ )
    {
      Ball* const ball = balls[bi];
      assert( ball->fdt_to_do );
      ball->force = ball->mass * gravity_accel;
      ball->torque = pVect(0,0,0);
    }

  const int n_surfaces = surfaces.size();
  for ( int i=0; i<n_surfaces; i++ )
    {
      Surface& s = surfaces[i];
      Ball& b0 = *s.b0;
      Ball& b1 = *s.b1;
      Ball& b2 = *s.b2;
      pNorm snorm = cross(b0.position,b1.position,b2.position);
      pVect vel = (1.0/3) * ( b0.velocity + b1.velocity + b2.velocity );
      float speed = dot(snorm,vel);
      float fmag =
        -speed * fabs(speed) * snorm.magnitude * 0.5f * opt_air_resistance;
      pVect f = fmag * snorm;
      b0.force += f;
      b1.force += f;
      b2.force += f;
    }

#pragma omp parallel for
  for ( int i=0; i<n_links; i++ )
    {
      Link* const link = links[i];
      if ( !link->is_simulatable ) continue;

      //
      // Spring Force from Neighbor Balls
      //

      Ball* const ball1 = link->ball1;
      Ball* const ball2 = link->ball2;

      /// Notes:
      //
      // link->cb1
      //   Coordinate where link touches ball1, in ball1's local space.
      //
      // ball1->omatrix
      //   Ball's orientation transformation ..
      //   .. which transforms vectors in ball1's local space to object space.
      //

      // Find position and velocity of the point where the link touches
      // the surface of ball 1 ...
      //
      pVect dir1 = ball1->omatrix * link->cb1;
      pCoor pos1 = ball1->position + dir1;
      pVect vel1 = ball1->velocity + ball1->point_rot_vel(dir1);

      // ... and ball 2.
      //
      pVect dir2 = ball2->omatrix * link->cb2;
      pCoor pos2 = ball2->position + dir2;
      pVect vel2 = ball2->velocity + ball2->point_rot_vel(dir2);

      // Construct a normalized (Unit) Vector from ball to neighbor
      // based on link connection points and ball centers.
      //
      pNorm link_dir(pos1,pos2);
      pNorm c_to_c(ball1->position,ball2->position);

      const float link_length = link_dir.magnitude;

      // Compute the speed of ball's end of link towards neighbor's end of link.
      //
      pVect delta_v = vel2 - vel1;
      float delta_s = dot( delta_v, link_dir );

      // Compute by how much the spring is stretched (positive value)
      // or compressed (negative value).
      //
      const float spring_stretch = link_length - link->distance_relaxed;

      // Determine whether spring is gaining energy (whether its length
      // is getting further from its relaxed length).
      //
      const bool gaining_e = ( delta_s > 0.0 ) == ( spring_stretch > 0 );

      const float blend_length_factor = 0.1; // Just a guess.
      const float blend_s = blend_length_factor * c_to_c.magnitude;
      const float blend = min( 1.0f, fabs(delta_s) / blend_s );

      // Use a smaller spring constant when spring is loosing energy,
      // a quick and dirty way of simulating energy loss due to spring
      // friction.
      //
      const float sc_base = link->spring_constant ?: opt_spring_constant;
      const float spring_constant =
        gaining_e ? sc_base :
        sc_base * ( 1.0f - blend + blend * opt_spring_damp );

      const float force_mag = spring_constant * spring_stretch;
      pVect spring_force_12 = force_mag * link_dir;

      // Apply forces affecting linear momentum.
      //
      link->spring_force_12 = spring_force_12;

      if ( ! link->is_surface_connection ) continue;

      // Apply torque.
      //
      link->torque1 = cross(pNorm(dir1), spring_force_12);
      link->torque2 = cross(pNorm(dir2), spring_force_12);
    }

  // Note: Because two links can reference the same ball this should
  // not be done in parallel.
  for ( Link *link: links )
    {
      if ( !link->is_simulatable ) continue;
      Ball* const ball1 = link->ball1;
      Ball* const ball2 = link->ball2;
      ball1->force += link->spring_force_12;
      ball2->force -= link->spring_force_12;
      ball1->torque += link->torque1;
      ball2->torque -= link->torque2;
    }

  ///
  /// Update Position of Each Ball
  ///

#pragma omp parallel for
  for ( int bi=0; bi<n_balls; bi++ )
    {
      Ball* const ball = balls[bi];
      if ( ball->locked )
        {
          ball->velocity = pVect(0,0,0);
          ball->omega = pVect(0,0,0);
          continue;
        }

      // Update Velocity
      //
      // This code assumes that force on ball is constant over time
      // step. This is clearly wrong when balls are moving with
      // respect to each other because the springs are changing
      // length. This inaccuracy will make the simulation unstable
      // when spring constant is large for the time step.
      //
      const float mass = max( ball->mass, ball->mass_min );

      pVect delta_v = ( ball->force / mass ) * delta_t;
      const float r = ball->radius;

      if ( platform_collision_possible(ball->position)
           && fabs(ball->position.y) < r )
        {
          const bool above = ball->position.y >= 0;
          const float spring_constant_plat =
            ball->velocity.y < 0 == above ? 100000 : 50000;
          const float fric_coefficient = opt_friction_dynamic;
          const float force_away = ( above ? 1 : -1 )
            * ( r - fabs(ball->position.y) ) * spring_constant_plat;
          const float delta_v_away = force_away / mass * delta_t;
          const float fric_force_mag = fric_coefficient * force_away;
          pNorm surface_v(ball->velocity.x,0,ball->velocity.z);
          const float delta_v_surf = fric_force_mag / mass * delta_t;

          if ( delta_v_surf > surface_v.magnitude )
            {
              // Ignoring other forces?
              delta_v =
                pVect(-ball->velocity.x,delta_v.y,-ball->velocity.z);
            }
          else
            {
              delta_v -= delta_v_surf * surface_v;
            }
          delta_v.y += delta_v_away;
        }

      ball->velocity += delta_v;

      // Air Resistance
      //
      const pNorm ball_dir(ball->velocity);
      const float delta_s_air_raw =
        opt_air_resistance * ball_dir.mag_sq * delta_t;
      const float delta_s_air = min( ball_dir.magnitude, delta_s_air_raw );
      ball->velocity -= delta_s_air * ball_dir;

      // Update Position
      //
      // Assume that velocity is constant.
      //
      ball->position += ball->velocity * delta_t;

      ball->omega += delta_t * ball->fdt_to_do * ball->torque;

      pNorm axis(ball->omega);

      // Update Orientation
      //
      // If ball isn't spinning fast skip expensive rotation.
      //
      if ( axis.mag_sq < 0.000001 ) continue;

      // Update ball orientation.
      //
      ball->orientation_set
        ( pQuat(axis,delta_t * axis.magnitude) * ball->orientation );
    }
}

bool
World::platform_collision_possible(pCoor pos)
{
  // Assuming no motion in x or z axes.
  //
  return pos.x >= platform_xmin && pos.x <= platform_xmax
    && pos.z >= platform_zmin && pos.z <= platform_zmax;
}

 /// External Modifications to State
//
//   These allow the user to play with state while simulation
//   running.

// Move the ball.
//
void Ball::translate(pVect amt) {position += amt;}

// Add velocity to the ball.
//
void Ball::push(pVect amt) {velocity += amt;}

// Set the velocity to zero.
//
void Ball::stop() {velocity = pVect(0,0,0); }

// Set the velocity and rotation (not yet supported) to zero.
//
void Ball::freeze() {velocity = pVect(0,0,0); }



void World::balls_translate(pVect amt,int b){head_ball->translate(amt);}
void World::balls_push(pVect amt,int b){head_ball->push(amt);}
void World::balls_translate(pVect amt)
{ for ( Ball *ball: balls ) ball->translate(amt); }
void World::balls_push(pVect amt)
{ for ( Ball *ball: balls ) ball->push(amt); }
void World::balls_stop()
{ for ( Ball *ball: balls ) ball->stop(); }
void World::balls_freeze(){balls_stop();}


void
World::render_link_start()
{
  if ( world_time_link_update == world_time ) return;

  // Empty the containers in preparation for a new pass.
  for ( auto b: lis_bufs ) b->clear();
}

void
World::render_link_gather(Link *link)
{
  // Place data needed to render link (the argument to this function)
  // in arrays. Just before rendering those arrays will be placed in
  // buffer objects.


  // If we've already collected links, just return.  This happens
  // when doing multiple passes to render shadows and reflections.
  //
  if ( world_time_link_update == world_time ) return;

  Ball *const ball1 = link->ball1;
  Ball *const ball2 = link->ball2;

  // Direction of link relative to center off ball.
  //
  pVect dir1 = ball1->omatrix * link->cb1;

  // Place on surface of ball where link connects.
  //
  pCoor pos1 = ball1->position + dir1;

  // Direction and connection location for ball 2.
  //
  pVect dir2 = ball2->omatrix * link->cb2;
  pCoor pos2 = ball2->position + dir2;

  pVect p1p2(pos1,pos2);
  pNorm p1p2n(p1p2);

  // Radius of link.
  //
  const float rad = link->radius ?: ball1->radius * 0.3;

  // Compute scale factors for texture.
  const float tex_aspect_ratio = 8.5 / 11;  // Yes, inches. Can't remember mm.
  const float page_width_o = 2.5;  // Width of texture in object-space coords.
  const float tex_t_scale = link->distance_relaxed / page_width_o;
  const float tex_angle_scale = rad * tex_aspect_ratio / page_width_o;

  pNorm dirn1(dir1);
  pNorm dirn2(dir2);

  // Vectors used to describe the cubic Hermite curve.
  //
  pVect4 v1 = p1p2n.magnitude * dirn1;
  pVect4 v2 = p1p2n.magnitude * dirn2;

  // Convert link's local x and y axes to global coordinates.
  //
  pVect4 b1_ydir = ball1->omatrix * link->b1_dir;
  pVect4 b2_ydir = ball2->omatrix * link->b2_dir;


  // Place data needed to render links in lists. Note that
  // the += operator appends the right-hand side to the container
  // on the left-hand side.
  //
  link->idx = lis_pos1.size();
  lis_pos1 << pos1;
  lis_pos2 << pos2;
  lis_v1 << v1;
  lis_v2 << v2;
  lis_b1_ydir << b1_ydir;
  lis_b2_ydir << b2_ydir;
  lis_tex_rad << vec4( tex_t_scale, tex_angle_scale, rad, 0);
}



void
World::render_link_render(vk::CommandBuffer& cb, bool shadow_volumes)
{
  // Perform a rendering pass to render all the links collected
  // by the routine render_link_2.


  // Check whether this is the first time we are rendering this data.
  //
  const bool first_render =
    world_time_link_update != world_time
    || opt_segments != buf_uni_common->links_n_segs
    || opt_sides != buf_uni_common->sides;
  world_time_link_update = world_time;
  const size_t n_instances = lis_pos1.size();
  const bool do_rt = vh.raytrace_active();
  const int segments = shadow_volumes ? max(2,opt_segments/4) : opt_segments;
  const int segments_sv = max(2,opt_segments/4);

  if ( !n_instances ) return;

  if ( first_render )
    {
      for ( auto b: lis_bufs ) b->to_dev();
      buf_uni_common->sides = opt_sides;
      buf_uni_common->links_n_segs_sv = segments_sv;
      buf_uni_common->delta_tee_sv = 1.0f / segments_sv;
      buf_uni_common->links_n_segs = opt_segments;
      buf_uni_common->delta_tee = 1.0f / opt_segments;
    }

  if ( do_rt && first_render )
    {
      const int n_tri_p_segment = 2 * opt_sides;
      const int n_tri_p_link = n_tri_p_segment * opt_segments;
      const int n_tri_p_pass = n_tri_p_link * n_instances;

      const int n_vtx_p_segment = 2 * ( opt_sides + 1 );
      const int n_vtx_p_link = n_vtx_p_segment * opt_segments;
      const int n_vtx_p_pass = n_vtx_p_link * n_instances;

      if ( !rt_geo_links )
        rt_geo_links
          .init( raytrace )
          .color_set( color_dim_gray, color_salmon )
          .bind( VR_BUF_IDX_sampler, sampler, texid_hw )
          .grouping_set( RT_Geometry_Indexed );

      buf_links_as_idx.assure(3*n_tri_p_pass);
      buf_links_tex_coor.assure(n_vtx_p_pass);
      buf_links_vertex_g.assure(n_vtx_p_pass);
      buf_links_normal_g.assure(n_vtx_p_pass);

      rt_geo_links
        .bind( VR_BUF_IDX_indices, buf_links_as_idx, 3*n_tri_p_pass )
        .bind( VR_BUF_IDX_pos, buf_links_vertex_g, n_vtx_p_pass )
        .bind( VR_BUF_IDX_normal, buf_links_normal_g )
        .bind( VR_BUF_IDX_tcoor, buf_links_tex_coor );
    }

  if ( do_rt )
    {
      QSet& qs = vh.qs;
      buf_uni_common.to_dev();
      qs.cb.begin( { vk::CommandBufferUsageFlagBits::eOneTimeSubmit } );
      comp_links.record_dispatch(qs.cb,opt_segments,n_instances);
      qs.cb.end();
      qs.q.submit( vk::SubmitInfo{ 0, nullptr, nullptr, 1, &qs.cb }, nullptr );
      qs.q.waitIdle();
      return;
    }


  VPipeline& pipe = shadow_volumes ? pipe_set_2_sv : pipe_set_2;

  pipe.record_draw_instanced(cb, segments+1, n_instances);
}


void
World::frame_callback()
{
  // This routine called whenever window needs to be updated, but
  // before commands are recorded.

  // Get any waiting keyboard commands.
  //
  cb_keyboard();

  if ( vh.have_raytrace ) vh.opt_want_raytrace_now = opt_want_raytrace;

  // This is a hack for 2023 HW 2.
  while ( !balls.empty() && balls.back()->world_time_expire <= world_time )
    balls.pop_back();
  //
  // Balls must be last in list and not part of a link.

  const double time_now = time_wall_fp();
  const int time_step_count_start = time_step_count;
  const int nmp = gpu_info.cuda_prop.multiProcessorCount;
  int bl_p_sm_max = 1024;
  for ( int i=0; i<gpu_info.num_kernels; i++ )
    set_min
      (bl_p_sm_max,
       gpu_info.get_max_active_blocks_per_mp(i, opt_block_size, 0) );
  const int grid_size = bl_p_sm_max * nmp;

  const bool opt_cuda = opt_physics_method != GP_cpu;
  if ( opt_cuda )
    {
      CE( cudaEventRecord(frame_start_ce,0) );
      data_cpu_to_gpu_constants();
      data_cpu_to_gpu_dynamic();
    }

  if ( !opt_pause || opt_single_frame || opt_single_time_step )
    {
      /// Advance simulation state.

      // Amount of time since the user saw the last frame.
      //
      const double wall_delta_t = time_now - last_frame_wall_time;

      // Compute amount by which to advance simulation state for this frame.
      //
      const double duration =
        vh.video_frame_duration() ?:
        opt_single_time_step ? opt_single_time_step * opt_time_step_duration :
        opt_single_frame ? 1/60.0 :
        wall_delta_t;

      const double world_time_target = world_time + duration;
      const double wall_time_limit = time_now + 0.05;

      while ( world_time < world_time_target )
        {
          if ( opt_cuda )
            {
              launch_time_step(opt_time_step_duration,grid_size,opt_block_size);
            }
          else
            {
              time_step_cpu(opt_time_step_duration);
            }
          time_step_count++;
          world_time += opt_time_step_duration;
          const double time_right_now = time_wall_fp();
          if ( time_right_now > wall_time_limit ) break;
        }

      // Reset these, just in case they were set.
      //
      opt_single_frame = false;
      opt_single_time_step = 0;
    }

  if ( opt_cuda )
    {
      data_gpu_to_cpu_dynamic();
      CE( cudaEventRecord(frame_stop_ce,0) );
    }

  frame_timer.phys_end();
  frame_timer.work_amt_set(time_step_count-time_step_count_start);

  float cuda_time = 0;
  if ( opt_cuda )
    {
      CE(cudaEventSynchronize(frame_stop_ce));
      CE(cudaEventElapsedTime(&cuda_time,frame_start_ce,frame_stop_ce));
    }
  frame_timer.cuda_frame_time_set(cuda_time);

  last_frame_wall_time = time_now;

}

int
main(int argv, char **argc)
{
  pVulkan_Helper pvulkan_helper(argv,argc);
  World world(pvulkan_helper,argv,argc);

  world.run();
}