PureBasic Forum
http://forums.purebasic.com/english/

Double Pendulum
http://forums.purebasic.com/english/viewtopic.php?f=16&t=70957
Page 1 of 1

Author:  Fig [ Sun Jul 01, 2018 8:05 pm ]
Post subject:  Double Pendulum

https://en.wikipedia.org/wiki/Double_pendulum
Feel free to change arms length, pendulum's mass, departure angles...
Code:
#ArmLength1=100:#ArmLength2=120
#Mass1=3:#Mass2=2
#MaxHistoricNumber=200
;departure angles
Angle1.f=-#PI/2:Angle2.f=-#PI/2
#x=800:#y=600
If InitSprite() = 0 Or InitKeyboard() = 0 Or InitMouse() = 0 Or OpenWindow(0, 0, 0, #X, #Y, "Double Pendulum", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)=0 Or OpenWindowedScreen(WindowID(0),0,0,#X,#Y,0,0,0,#PB_Screen_WaitSynchronization)=0
    MessageRequester("Error", "Can't open the sprite system", 0)
    End
EndIf
Structure xy
    x.i
    y.i
EndStructure
NewList trace.xy()
x0.i=#x/2:y0.i=#y/2:x1.i=0:y1.i=0:x2.i=0:y2.i=0
acceleration1.f=0.0:acceleration2.f=0.0
vitesse1.f=0:vitesse2.f=0
g.f=1
Repeat
    FlipBuffers()
    ClearScreen(#Black)
    ExamineKeyboard()
   
    num1.f=-g*(2*#Mass1+#Mass2)*Sin(Angle1)
    num2.f=-g*#Mass2*Sin(Angle1-2*Angle2)
    num3.f=-2*#Mass2*Sin(Angle1-Angle2)
    num4.f=vitesse2*vitesse2*#ArmLength2+vitesse1*vitesse1*#ArmLength1*Cos(Angle1-Angle2)
    num.f=num1+num2+num3*num4
    den.f=2*#Mass1+#Mass2-#Mass2*Cos(2*Angle1-2*Angle2)
    den.f=den*#ArmLength2
    acceleration1=num/den
   
    num1.f=2*Sin(Angle1-Angle2)
    num2.f=vitesse1*vitesse1*#ArmLength1*(#Mass1+#Mass2)
    num3.f= g*(#Mass1+#Mass2)*Cos(Angle1)
    num4.f=vitesse2*vitesse2*#ArmLength2*#Mass2*Cos(Angle1-Angle2)
    num.f=(num2+num3+num4)
    num=num1*num
    den.f  =#ArmLength2*(2*#Mass1+#Mass2-(#Mass2*Cos(2*Angle1-2*Angle2)))
    acceleration2=num/den
   
    vitesse1+acceleration1
    vitesse2+acceleration2
    Angle1+vitesse1
    Angle2+vitesse2
    x1=x0+#ArmLength1*Sin(Angle1)
    y1=y0+#ArmLength1*Cos(Angle1)
    x2=x1+#ArmLength2*Sin(Angle2)
    y2=y1+#ArmLength2*Cos(Angle2)
    If FirstElement(trace())
        xold.i=trace()\x:yold.i=trace()\y
        If ListSize(trace())>#MaxHistoricNumber:DeleteElement(trace()):EndIf
        LastElement(trace())
    EndIf
    AddElement(Trace())
    Trace()\x=x2
    Trace()\y=y2
   
    StartDrawing(ScreenOutput())
   
    ForEach trace()
        colorR.i=ListIndex(trace())*255/ListSize(trace())
        LineXY(trace()\x,trace()\y,xold,yold,RGB(colorR,0,0))
        xold=trace()\x:yold=trace()\y
    Next
    LineXY(x0,y0,x1,y1,#White)
    LineXY(x1,y1,x2,y2,#White)
    Circle(x1,y1,#Mass1)
    Circle(x2,y2,#Mass2)
    StopDrawing()
    While WindowEvent():Wend
Until KeyboardPushed(#PB_Key_Escape)

Author:  luis [ Sun Jul 01, 2018 11:36 pm ]
Post subject:  Re: Double Pendulum

Very nice, I like it.
Nice how the trace slowly fades away.
Interesting demo :)

Author:  kvitaliy [ Mon Jul 02, 2018 6:00 am ]
Post subject:  Re: Double Pendulum

Very nice!

Author:  IceSoft [ Mon Jul 02, 2018 8:33 am ]
Post subject:  Re: Double Pendulum

little bit off topic:
I will write a Chipmunk4PB example too.
Thanks for the idea.

Author:  RSBasic [ Mon Jul 02, 2018 10:10 am ]
Post subject:  Re: Double Pendulum

Great Image

Author:  idle [ Mon Jul 02, 2018 11:02 am ]
Post subject:  Re: Double Pendulum

That's really great.

Author:  davido [ Mon Jul 02, 2018 7:21 pm ]
Post subject:  Re: Double Pendulum

@Fig,
Another fine demonstration.
Thank you.

Author:  IceSoft [ Thu Jul 05, 2018 1:45 pm ]
Post subject:  Re: Double Pendulum

Here the source part (not really optimized) of the "Double Pendulum" for the Chipmunk4PB wrapper framework.
Hope I can upload the next Chipmunk4PB wrapper framework version in few days.

Code:
Procedure initFuncDoublePendulum()
  #ArmLength1 = 100:
  #ArmLength2 = 120
 
  initFunc()
  cpSpaceSetGravity(*space, cpv(A,0, 100));
  cpSpaceSetIterations(*space, 120) 
 
  ; First stick
  *pos1.cpVect = cpv(pos1.cpVect, wx/2,  200);     
  *sp1.cpSprite =  CreateSegmentElement(*space, *pos1\x,  *pos1\y, *pos1\x + #ArmLength1, *pos1\y, 0, #NormalBody)
  cpShapeSetMass(*sp1\shape, 10)
  *constraint1 = cpSpaceAddConstraint(*space, cpPinJointNew(*sp1\body, cpSpaceGetStaticBody(*space), *pos1, *pos1));
 
  ; Second stick
  *pos2.cpVect = cpv(pos2.cpVect, *pos1\x + #ArmLength1,  *pos1\y);
  *sp2.cpSprite = CreateSegmentElement(*space, *pos2\x,  *pos2\y, *pos2\x + #ArmLength2, *pos2\y, 0, #NormalBody)
  cpShapeSetMass(*sp2\shape, 20)
  *constraint2 = cpSpaceAddConstraint(*space, cpPinJointNew(*sp2\body, *sp1\body, *pos2, *pos2));
EndProcedure

Procedure updateFuncDoublePendulum()
  updateFunc()
EndProcedure

Procedure destroyFuncDoublePendulum()
  destroyFunc()
EndProcedure

AddElement(ChipmunkDemo())
ChipmunkDemo()\title = "Double Pendulum"
ChipmunkDemo()\timestep = 1.0/180
ChipmunkDemo()\initFunc = @initFuncDoublePendulum()
ChipmunkDemo()\updateFunc = @updateFuncDoublePendulum()
ChipmunkDemo()\destroyFunc = @destroyFuncDoublePendulum()

Author:  Michael Vogel [ Thu Jul 05, 2018 11:52 pm ]
Post subject:  Re: Double Pendulum

Wow :shock:

Well done, even with some quirks (e.g. #ArmLength2=35)...
If anyone is able to predict the behaviour of a triple pendulum I will get nervous...

Author:  IceSoft [ Fri Jul 06, 2018 6:40 am ]
Post subject:  Re: Double Pendulum

Hint:
A n Pendulum with many small sections "simulates" a rope ;-)
I added an example to the Chipmunk4PB framework.

Author:  Michael Vogel [ Sat Jul 07, 2018 10:11 am ]
Post subject:  Re: Double Pendulum

I've checked if vector drawing would be fast enough for Fig's nice code (replaced also the sprite commands by default event handling)...

Press 'Space' to change the drawing mode:

Code:
Global ActivePendulum
Global DrawingMode
Global DrawingModeChange

Procedure LineXY_(x1,y1,x2,y2,color)

   If DrawingMode
      MovePathCursor(x1,y1)
      AddPathLine(x2,y2)
      VectorSourceColor($FF000000|color)
      StrokePath(2)

   Else
      LineXY(x1,y1,x2,y2,color)
   EndIf

EndProcedure
Procedure Circle_(x,y,radius)

   If DrawingMode
      AddPathCircle(x,y,radius)
      ;VectorSourceColor($FF000000|color)
      FillPath()

   Else
      Circle(x,y,radius)

   EndIf

EndProcedure

Procedure Pendulum(nil)

   #ArmLength1=100:#ArmLength2=120
   #Mass1=7:#Mass2=4
   #MaxHistoricNumber=200

   ;departure angles
   Angle1.f=-#PI/2:Angle2.f=-#PI/2
   #x=800:#y=600

   Structure xy
      x.i
      y.i
   EndStructure
   NewList trace.xy()

   x0.i=#x/2:y0.i=#y/2:x1.i=0:y1.i=0:x2.i=0:y2.i=0
   acceleration1.f=0.0:acceleration2.f=0.0
   vitesse1.f=0:vitesse2.f=0
   g.f=1

   topleft.xy
   topright.xy

   ActivePendulum=#True

   While ActivePendulum

      num1.f=-g*(2*#Mass1+#Mass2)*Sin(Angle1)
      num2.f=-g*#Mass2*Sin(Angle1-2*Angle2)
      num3.f=-2*#Mass2*Sin(Angle1-Angle2)
      num4.f=vitesse2*vitesse2*#ArmLength2+vitesse1*vitesse1*#ArmLength1*Cos(Angle1-Angle2)
      num.f=num1+num2+num3*num4
      den.f=2*#Mass1+#Mass2-#Mass2*Cos(2*Angle1-2*Angle2)
      den.f=den*#ArmLength2
      acceleration1=num/den

      num1.f=2*Sin(Angle1-Angle2)
      num2.f=vitesse1*vitesse1*#ArmLength1*(#Mass1+#Mass2)
      num3.f= g*(#Mass1+#Mass2)*Cos(Angle1)
      num4.f=vitesse2*vitesse2*#ArmLength2*#Mass2*Cos(Angle1-Angle2)
      num.f=(num2+num3+num4)
      num=num1*num
      den.f  =#ArmLength2*(2*#Mass1+#Mass2-(#Mass2*Cos(2*Angle1-2*Angle2)))
      acceleration2=num/den

      vitesse1+acceleration1
      vitesse2+acceleration2
      Angle1+vitesse1
      Angle2+vitesse2
      x1=x0+#ArmLength1*Sin(Angle1)
      y1=y0+#ArmLength1*Cos(Angle1)
      x2=x1+#ArmLength2*Sin(Angle2)
      y2=y1+#ArmLength2*Cos(Angle2)
      If FirstElement(trace())
         xold.i=trace()\x:yold.i=trace()\y
         If ListSize(trace())>#MaxHistoricNumber:DeleteElement(trace()):EndIf
         LastElement(trace())
      EndIf
      AddElement(Trace())
      Trace()\x=x2
      Trace()\y=y2

      If DrawingMode
         StartVectorDrawing(CanvasVectorOutput(0))
         FillVectorOutput()   
      Else
         StartDrawing(CanvasOutput(0))
         Box(0,0,#X,#Y,#Black)
      EndIf

      #Marker=$CBA762
      #Taboo=#ArmLength1/2
      topleft\y=9999
      topright\y=9999

      ForEach trace()
         colorR.i=ListIndex(trace())*255/ListSize(trace())
         LineXY_(trace()\x,trace()\y,xold,yold,RGB(colorR,0,0))
         xold=trace()\x
         yold=trace()\y
         If xold<x0-#Taboo
            If yold<topleft\y
               topleft\y=yold
               topleft\x=xold
            EndIf
         ElseIf xold>x0+#Taboo
            If yold<topright\y
               topright\y=yold
               topright\x=xold
            EndIf
         EndIf
      Next

      LineXY_(topleft\x-6,topleft\y,topleft\x+6,topleft\y,#Marker)
      LineXY_(topright\x-6,topright\y,topright\x+6,topright\y,#Marker)

      LineXY_(x0,y0,x1,y1,#White)
      LineXY_(x1,y1,x2,y2,#White)

      Circle_(x1,y1,#Mass1)
      Circle_(x2,y2,#Mass2)

      If DrawingMode
         StopVectorDrawing()
      Else
         StopDrawing()
      EndIf

      DrawingMode=DrawingModeChange
   Wend

   ActivePendulum=666

EndProcedure
Procedure Exit()
   ActivePendulum=#Null
   While ActivePendulum<>666
      Delay(1)
   Wend
   Delay(5)
   End
EndProcedure

OpenWindow(0,0,0,#X,#Y,"Double Pendulum",#PB_Window_SystemMenu|#PB_Window_ScreenCentered)
CanvasGadget(0,0,0,#X,#Y)
AddKeyboardShortcut(0,#PB_Shortcut_Escape,666)
AddKeyboardShortcut(0,#PB_Shortcut_Space,1000)

CreateThread(@Pendulum(),#Null)

Repeat

   Select WaitWindowEvent()
   Case #PB_Event_CloseWindow
      Exit()
   Case #PB_Event_Menu
      If EventMenu()=666
         Exit()
      Else
         DrawingModeChange!1
      EndIf
   EndSelect

ForEver


Author:  IceSoft [ Mon Jul 09, 2018 8:44 am ]
Post subject:  Re: Double Pendulum

Here are a Chipmunk2D Example (exe only! for the moment)
You can set 2-4 Sticks with the [-/+] Keys

Double Pendulum(screenshot):
Image

4 Pendulum(screenshot):]
Image

Zip (Exe only)
https://drive.google.com/file/d/1Moarsl ... sp=sharing

Author:  zefiro_flashparty [ Fri Sep 07, 2018 8:51 pm ]
Post subject:  Re: Double Pendulum

:lol:
very interesting that,

Author:  Psychophanta [ Sun Oct 07, 2018 12:01 pm ]
Post subject:  Re: Double Pendulum

Great one. :!:

Page 1 of 1 All times are UTC + 1 hour
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group
http://www.phpbb.com/