float delta_t=0.05;
  float theta1,theta2,dtheta1,dtheta2;
  float x1,x2,y1,y2;
  
  double_pendul p1,p2;

void setup(){
  theta1 = 2.5;
  theta2 = 2.5;
  theta1 += random(0.64);
  theta2 += random(0.64);
  dtheta1 = 0; dtheta2 = 0;
  size(600,300);
  p1=new double_pendul(theta1,theta2,dtheta1,dtheta2,150.,150.,70.,20.);
  theta2 += random(0.0001);
  p2=new double_pendul(theta1,theta2,dtheta1,dtheta2,450.,150.,70.,20.);
  }
void draw(){

  background(0);
  p1.move();
  p2.move();
  p1.display();
  p2.display();
}

class double_pendul
{
  float theta1, theta2, dtheta1, dtheta2;
  float x1, y1, x2, y2;
  float xc,yc;
  float l,r;
  
 double_pendul(float th1, float th2, float dth1, float dth2, 
  float x, float y, float len, float rr){
  theta1=th1; theta2=th2; dtheta1=dth1; dtheta2=dth2;
  xc=x; yc=y; l=len; r=rr;
 } 
 void move(){ 
 float k1x,k1v,k2x,k2v,k3x,k3v,k4x,k4v;
 float l1x,l1v,l2x,l2v,l3x,l3v,l4x,l4v;

   k1x = delta_t* dtheta1;
  l1x = delta_t* dtheta2;
  k1v = delta_t*fnal1(theta1,theta2,dtheta1,dtheta2);
  l1v = delta_t*fnal2(theta1,theta2,dtheta1,dtheta2);
  k2x = delta_t*(dtheta1+k1v/2);
  l2x = delta_t*(dtheta2+l1v/2);
  k2v = delta_t*fnal1(theta1+k1x/2,theta2+l1x/2,dtheta1+k1v/2,dtheta2+l1v/2);
  l2v = delta_t*fnal2(theta1+k1x/2,theta2+l1x/2,dtheta1+k1v/2,dtheta2+l1v/2);
  k3x = delta_t*(dtheta1+k2v/2);
  l3x = delta_t*(dtheta2+l2v/2);
  k3v = delta_t*fnal1(theta1+k2x/2,theta2+l2x/2,dtheta1+k2v/2,dtheta2+l2v/2);
  l3v = delta_t*fnal2(theta1+k2x/2,theta2+l2x/2,dtheta1+k2v/2,dtheta2+l2v/2);
  k4x = delta_t*(dtheta1+k3v);
  l4x = delta_t*(dtheta2+l3v);
  k4v = delta_t*fnal1(theta1+k3x,theta2+l3x,dtheta1+k3v,dtheta2+l3v);
  l4v = delta_t*fnal2(theta1+k3x,theta2+l3x,dtheta1+k3v,dtheta2+l3v);
  theta1 += k1x/6 + k2x/3 + k3x/3 + k4x/6;
  theta2 += l1x/6 + l2x/3 + l3x/3 + l4x/6;
  dtheta1 += k1v/6 + k2v/3 + k3v/3 + k4v/6;
  dtheta2 += l1v/6 + l2v/3 + l3v/3 + l4v/6;
  x1 = sin(theta1);
  y1 = -cos(theta1);
  x2 = x1+sin(theta2);
  y2 = y1-cos(theta2);
}

 float fnd(float a, float b, float c, float d){
 return a*d-b*c;}
float fnm1(float t1,float t2,float a2)
{ return -0.5*sin(t1-t2)*a2*a2 - sin(t1);}
float fnm2(float t1,float t2,float a1)
{ return sin(t1-t2)*a1*a1-sin(t2);}
float fndenm(float t1,float t2)
{return 1-0.5*cos(t1-t2)*cos(t1-t2);}
float fnal1(float t1,float t2,float a1,float a2)
{ return fnd(fnm1(t1,t2,a2),0.5*cos(t1-t2),fnm2(t1,t2,a1),1)/fndenm(t1,t2);}
float fnal2(float t1,float t2,float a1,float a2)
{ return fnd(1,fnm1(t1,t2,a2),cos(t1-t2),fnm2(t1,t2,a1))/fndenm(t1,t2);}

void display() {
  stroke(255);
  line(xc,yc,l*x1+xc,yc-l*y1);
  line(l*x1+xc,yc-l*y1,l*x2+xc,yc-l*y2);
  
  noStroke();
  ellipse(l*x1+xc,yc-l*y1,r,r);
  ellipse(l*x2+xc,yc-l*y2,r,r);
}
}