Ring and bar pendulums

sketch_241002a_ringbar.pde

// Ring and bar double pendulums
//
// Mathematics and Processing code by Waldeck Schützer (@infinitymathart)
import java.lang.Math;

boolean recording = false;
boolean use_border = true;
int FPS = 60;
int duration = 3*60;
float text_t = 20.0*3/2;

int num_rows = 5;
int num_cols = 5;
int border = 50;
float col_width, row_height;
float h_margin, v_margin;
float x_left, y_top;
float t = 0;

// Constants
float g = 9.8;       // Gravitational acceleration
float L1 = 30;      // Length of the bar
float L2 = 8;       // Radius of the ring
float d = 1.0/3.0;   // Distance of the first pivot along the bar (1/3 of the length of the bar)
float rho1 = 0.1; // Density of bar
float rho2 = 0.1; // Density of ring
float m1 = rho1*L1;        // Mass of the bar
float m2 = rho2*TAU*L2;        // Mass of the ring
int oversampling = recording ? 100 : 24;

// Moments of inertia
float I1 = (m1 * L1 * L1) / 9.0;  // Moment of inertia of the bar around the pivot at 1/3 of its length
float I2 = m2 * L2 * L2;          // Moment of inertia of the ring


PFont courier;
PFont times;

// Time step
float dt = 0.05;
float fac = 1.0;
int total_frames = duration * FPS;

class pendulum
{
  // Initial conditions
  double theta1 = 0;  // Initial angle of the bar (with respect to vertical)
  double theta2 = 0;  // Initial angle of the ring (with respect to vertical)
  double theta1_dot = 0;   // Initial angular velocity of the bar
  double theta2_dot = 0;   // Initial angular velocity of the ring
  float scl = 10*fac;  // Drawing scale
  float x0 = 0; // Pendulum position
  float y0 = 0;

void step(double dt)
{
    // Compute the current accelerations using the current state
    double deltaTheta = theta1 - theta2;

    // Compute accelerations for theta1 and theta2 based on current angles and velocities
    double theta1_ddot = (-g * (2 * m1 + m2) * Math.sin(theta1)
                         - m2 * g * Math.sin(theta1 - 2 * theta2)
                         - 2 * Math.sin(deltaTheta) * m2 * (theta2_dot * theta2_dot * L2 + theta1_dot * theta1_dot * L1 * Math.cos(deltaTheta)))
                        / I1;

    double theta2_ddot = (2 * Math.sin(deltaTheta) * (theta1_dot * theta1_dot * L1 * (m1 + m2)
                         + g * (m1 + m2) * Math.cos(theta1)
                         + theta2_dot * theta2_dot * L2 * m2 * Math.cos(deltaTheta)))
                        / I2;

    // Predictor step: Use explicit velocity to estimate the new positions
    double theta1_pred = theta1 + theta1_dot * dt + 0.5 * theta1_ddot * dt * dt;
    double theta2_pred = theta2 + theta2_dot * dt + 0.5 * theta2_ddot * dt * dt;

    // Corrector step: Iteratively solve for theta1_next and theta2_next
    double theta1_next = theta1_pred;
    double theta2_next = theta2_pred;

    // Use an iterative solver (e.g., Newton's method or fixed-point iteration)
    for (int i = 0; i < 10; i++) {
        // Compute new accelerations at the predicted future positions
        double theta1_ddot_next = (-g * (2 * m1 + m2) * Math.sin(theta1_next)
                                  - m2 * g * Math.sin(theta1_next - 2 * theta2_next)
                                  - 2 * Math.sin(theta1_next - theta2_next) * m2 * (theta2_dot * theta2_dot * L2 + theta1_dot * theta1_dot * L1 * Math.cos(theta1_next - theta2_next)))
                                 / I1;

        double theta2_ddot_next = (2 * Math.sin(theta1_next - theta2_next) * (theta1_dot * theta1_dot * L1 * (m1 + m2)
                                  + g * (m1 + m2) * Math.cos(theta1_next)
                                  + theta2_dot * theta2_dot * L2 * m2 * Math.cos(theta1_next - theta2_next)))
                                 / I2;

        // Update the next positions using the implicit formula
        theta1_next = theta1 + theta1_dot * dt + 0.5 * theta1_ddot_next * dt * dt;
        theta2_next = theta2 + theta2_dot * dt + 0.5 * theta2_ddot_next * dt * dt;
    }

    // After convergence, update velocities using the new accelerations
    double theta1_ddot_new = (-g * (2 * m1 + m2) * Math.sin(theta1_next)
                             - m2 * g * Math.sin(theta1_next - 2 * theta2_next)
                             - 2 * Math.sin(theta1_next - theta2_next) * m2 * (theta2_dot * theta2_dot * L2 + theta1_dot * theta1_dot * L1 * Math.cos(theta1_next - theta2_next)))
                            / I1;

    double theta2_ddot_new = (2 * Math.sin(theta1_next - theta2_next) * (theta1_dot * theta1_dot * L1 * (m1 + m2)
                             + g * (m1 + m2) * Math.cos(theta1_next)
                             + theta2_dot * theta2_dot * L2 * m2 * Math.cos(theta1_next - theta2_next)))
                            / I2;

    theta1_dot += 0.5 * (theta1_ddot + theta1_ddot_new) * dt;
    theta2_dot += 0.5 * (theta2_ddot + theta2_ddot_new) * dt;

    // Finally, update positions
    theta1 = theta1_next;
    theta2 = theta2_next;
}
 
  // Velocity Verlet Integrator (symplectic)
  void step_verlet(float dt)
  {
    // Compute the current accelerations using the equations of motion
    double deltaTheta = theta1 - theta2;

    // Bar acceleration
    double theta1_ddot = (-g * (2 * m1 + m2) * Math.sin(theta1)
                         - m2 * g * Math.sin(theta1 - 2 * theta2)
                         - 2 * Math.sin(deltaTheta) * m2 * (theta2_dot * theta2_dot * L2 + theta1_dot * theta1_dot * L1 * Math.cos(deltaTheta)))
                        / I1;

    // Ring acceleration
    double theta2_ddot = (2 * Math.sin(deltaTheta) * (theta1_dot * theta1_dot * L1 * (m1 + m2)
                         + g * (m1 + m2) * Math.cos(theta1)
                         + theta2_dot * theta2_dot * L2 * m2 * Math.cos(deltaTheta)))
                        / I2;

    // Update positions (angles) using the current velocities and accelerations
    theta1 += theta1_dot * dt + 0.5 * theta1_ddot * dt * dt;
    theta2 += theta2_dot * dt + 0.5 * theta2_ddot * dt * dt;

    // Compute the new accelerations based on the updated positions
    deltaTheta = theta1 - theta2;
    
    double new_theta1_ddot = (-g * (2 * m1 + m2) * Math.sin(theta1)
                             - m2 * g * Math.sin(theta1 - 2 * theta2)
                             - 2 * Math.sin(deltaTheta) * m2 * (theta2_dot * theta2_dot * L2 + theta1_dot * theta1_dot * L1 * Math.cos(deltaTheta)))
                            / I1;

    double new_theta2_ddot = (2 * Math.sin(deltaTheta) * (theta1_dot * theta1_dot * L1 * (m1 + m2)
                             + g * (m1 + m2) * Math.cos(theta1)
                             + theta2_dot * theta2_dot * L2 * m2 * Math.cos(deltaTheta)))
                            / I2;

    // Update velocities using the average of the old and new accelerations
    theta1_dot += 0.5 * (theta1_ddot + new_theta1_ddot) * dt;
    theta2_dot += 0.5 * (theta2_ddot + new_theta2_ddot) * dt;
  }

  // Symplectic Euler Integrator
  void step_euler(float dt)
  {
    // Compute the accelerations using the moments of inertia and equations of motion
    double deltaTheta = theta1 - theta2;
  
    // Bar acceleration (considering the inertia of the bar and gravity)
    double theta1_ddot = (-g * (2 * m1 + m2) * Math.sin(theta1)
                         - m2 * g * Math.sin(theta1 - 2 * theta2)
                         - 2 * Math.sin(deltaTheta) * m2 * (theta2_dot * theta2_dot * L2 + theta1_dot * theta1_dot * L1 * Math.cos(deltaTheta)))
                        / I1;
  
    // Ring acceleration (considering the inertia of the ring and the swinging bar)
    double theta2_ddot = (2 * Math.sin(deltaTheta) * (theta1_dot * theta1_dot * L1 * (m1 + m2)
                         + g * (m1 + m2) * Math.cos(theta1)
                         + theta2_dot * theta2_dot * L2 * m2 * Math.cos(deltaTheta)))
                        / I2;
  
    // Symplectic integration: update velocities first using the accelerations
    theta1_dot += theta1_ddot * dt;
    theta2_dot += theta2_ddot * dt;

    // Then, update the positions (angles) using the updated velocities
    theta1 += theta1_dot * dt;
    theta2 += theta2_dot * dt;

  }
  
  void draw()
  {
    pushMatrix();
    translate(x0, y0);
    // Compute positions of the bar and ring
    float x1 = scl*(2 * L1 / 3) * sin((float)theta1); // Position of the endpoint of the bar where the ring is attached
    float y1 = -scl*(2 * L1 / 3) * cos((float)theta1);
    float x2 = -scl*(L1 / 3) * sin((float)theta1);    // Position of the pivot connecting the ring
    float y2 = scl*(L1 / 3) * cos((float)theta1);
    float x3 = x2 - scl*L2 * sin((float)theta2);      // Position of the ring
    float y3 = y2 + scl*L2 * cos((float)theta2);
  
    // Draw the bar
    stroke(255);
    strokeWeight(0.5*fac);
    line(x1, y1, x2, y2); // Line representing the bar from its pivot point to the ring pivot
  
    // Draw the ring as a circle
    noFill();
    stroke(255);
    strokeWeight(1*fac);
    ellipse(x3, y3, scl*L2 * 2, scl*L2 * 2); // Ring
  
    // Draw the pivot points
    fill(255);
    noStroke();
    ellipse(0, 0, scl*2, scl*2); // Fixed pivot point
    ellipse(x2, y2, scl*2, scl*2); // Pivot where the ring connects
    popMatrix();
  }

}

pendulum[] pendulums = new pendulum[num_rows*num_cols];

float ease(float p, float g) {
  if (p < 0.5) 
    return 0.5 * pow(2*p, g);
  else
    return 1 - 0.5 * pow(2*(1 - p), g);
}

float thereandback(float p)
{
  return constrain(p*(1-p)*4,0,1);
}

void settings()
{
  if (recording)
     size(2160, 2160);
  else
     size(800, 800, P2D);
  smooth(8);
}

void setup() 
{
  randomSeed(1337);
  fac = width/800.0;
  border = int(border * fac);
  
  col_width = width/(num_cols+1);
  row_height = height/(num_rows+1);
  h_margin = width - col_width * num_cols/2;
  v_margin = height - row_height * num_rows/2;
  x_left = (h_margin - width);
  y_top = (v_margin - height);
  courier = createFont("Courier New",48*fac,true);
  times = createFont("Times New Roman",24*fac,false);
  
  background(0);
  float rho = 0.00125*TAU/(num_rows*num_cols);
  for(int i=0; i<num_rows*num_cols; i++)
  {
    pendulums[i] = new pendulum();
    pendulums[i].x0 = (i % num_rows)*col_width + col_width/2 + x_left;
    pendulums[i].y0 = (i / num_rows)*row_height + row_height/2 + y_top;
    pendulums[i].scl = 2*fac;
    pendulums[i].theta1 = PI/2 + i*rho;
    pendulums[i].theta2 = PI/2 + i*rho;
  }
}

void draw() 
{
  t = (float)(frameCount-1)/total_frames;
  blendMode(BLEND);
  fill(0,50);
  noStroke();
  rect(0,0,width,height);
  translate(width / 2, height / 2); // Translate to center of screen
  //lights();

  // step the model with oversampling
  for(pendulum p: pendulums)
  {
    for(int i=0; i<oversampling; i++)
       p.step(dt/oversampling);
    p.draw();
  }
  // texts
  fill(255);stroke(255);
  textAlign(CENTER, CENTER);
  textFont(courier);
  textSize(15*fac);
  text("@infinitymathart • 2024",0.0*width,0.40*height);

  textFont(times);
  float t1 = thereandback( (text_t*t));
  if (t1 >0)
  {
    fill(255,255*t1);
    text("These are ring and bar double pendulums.", 0.0*width, -0.4*height);
  }
  t1 = thereandback( (text_t*t-1) );
  if (t1>0)
  {
    fill(255,255*t1);
    text("Although deterministic, they exhibit unpredictable behavior.", 0.0*width, -0.4*height);
  }
  t1 = thereandback( (text_t*t-2) );
  if (t1>0)
  {
    fill(255,255*t1);
    text("These started with only a minute difference in position,", 0.0*width, -0.4*height);
  }
  t1 = thereandback( (text_t*t-3) );
  if (t1>0)
  {
    fill(255,255*t1);
    text("and initially appear to be moving in perfect sync,", 0.0*width, -0.4*height);
  }
  t1 = thereandback( (text_t*t-4) );
  if (t1>0)
  {
    fill(255,255*t1);
    text("but will soon diverge into mesmerizing chaos.", 0.0*width, -0.4*height);
  }
  t1 = thereandback( (text_t*t-5) );
  if (t1>0)
  {
    fill(255,255*t1);
    text("Keep watching and prepare to be amazed!", 0.0*width, -0.4*height);
  }

  if(use_border)
  {
    fill(0);
    noStroke();
    pushMatrix();
    translate(-width/2,-height/2);
    rect(0, 0, width, border);
    rect(0, 0, border, height);
    rect(0, height-border, width, border);
    rect(width-border, 0, border, height);
    noFill();
    stroke(255);
    strokeWeight(1);
    rect(border,border,width-2*border,height-2*border);
    popMatrix();
  }

  if (recording)
  {
     println(frameCount,"/",total_frames);
     saveFrame("/tmp/p/frame_####.png");
     if (frameCount == total_frames)
     {
       stop();
       exit();
     }
  }

}

/* License:
 *
 * Copyright (c) 2024 Waldeck Schützer
 *
 * All rights reserved.
 *
 * This code after the template and the related animations are the property of the
 * copyright holder. Any reproduction, distribution, or use of this material,
 * in whole or in part, without the express written permission of the copyright
 * holder is strictly prohibited.
 */