PureBasic Forums - English

Posted: **Wed May 06, 2009 12:36 pm**

Hi,

I've just written this spring simulation in an attempt to teach myself about differential equations. It seems to work (you can mess around with the globals, damping/stiffness etc and watch the changes) but I just wondered (since I am terrible at maths) if it is mathematically correct? Maybe there is a math geek around who can look over the code? The thing that is particularly bothering me is that the spring's acceleration is a function of both velocity and position, so I wasn't sure whether to pass velocity to the integrator as its x value, or position?

Code: Select all

EnableExplicit

Prototype.d DerivativeFunction(x.d,y.d)

Structure Timer
  Freq.l
  Ticks.l
EndStructure

Declare.l CreateTimer(Freq)
Declare.l WaitTimer(*Timer.Timer)

Declare.l UpdateSpring(t.d)
Declare.l DrawSpring()

Declare.d Integrate(x.d,y.d,h.d,f.DerivativeFunction)
Declare.d SpringVelocity(v.d,y.d)
Declare.d SpringAcceleration(v.d,y.d)

#FPS = 30
#SCREEN_WIDTH = 256
#SCREEN_HEIGHT = 256

Global SPRING_MASS.d = 32.0
Global SPRING_DAMPING.d = 2.0
Global SPRING_STIFFNESS.d = 10.0
Global SPRING_VELOCITY.d = 0.0
Global SPRING_Y_POS.d = 80.0
Global SPRING_REST_LENGTH.d = 80.0

Global SPRING_NUM_COILS.l = 8
Global SPRING_COIL_WIDTH = 16

Define EventID.l

Define *Timer = CreateTimer(1000 / #FPS)
Define Pings.l
Define i.l

If InitSprite() = 0 ;Or InitKeyboard() = 0
  MessageRequester("OOPS", "Can't initialise sprite stuff.", 0)
  End
EndIf


If Not OpenWindow(0, 100, 200, 256, 256, "Spring Simulator")
  MessageRequester("OOPS", "Can't open main window.", 0)
  End
EndIf

If Not OpenWindowedScreen(WindowID(0),0,0,#SCREEN_WIDTH,#SCREEN_HEIGHT,0,0,0)
  MessageRequester("OOPS", "Can't open screen.", 0)
  End 
EndIf


Repeat
   
  Pings = WaitTimer(*Timer)
  For i = 1 To Pings
    UpdateSpring(0.5)
  Next
   

  EventID.l = WindowEvent()


  Select EventID
    Case #PB_Event_CloseWindow
      End
  EndSelect

	ClearScreen(RGB(0,0,0))

  If StartDrawing(ScreenOutput())
    DrawingMode(#PB_2DDrawing_Transparent)
    DrawSpring()
    StopDrawing()
    FlipBuffers()
  EndIf


ForEver

End

Procedure.l UpdateSpring(t.d)
  SPRING_VELOCITY = SPRING_VELOCITY + Integrate(SPRING_VELOCITY,0.0,t,@SpringAcceleration())
  SPRING_Y_POS = SPRING_Y_POS + Integrate(SPRING_VELOCITY,0.0,t,@SpringVelocity())
EndProcedure

Procedure.l DrawSpring()

    Protected CoilHeight.d = (SPRING_Y_POS + SPRING_REST_LENGTH) / SPRING_NUM_COILS
    Protected CoilY.d = 16.0 + 0.5 * CoilHeight
    Protected CoilX.d = #SCREEN_WIDTH / 2.0 - 0.5 * SPRING_COIL_WIDTH
    Protected i.l
    
    LineXY(#SCREEN_WIDTH / 2.0,16.0,#SCREEN_WIDTH / 2.0 - 0.5 * SPRING_COIL_WIDTH,16.0 + 0.5 * CoilHeight,RGB(0,255,0))
    
    For i = 0 To SPRING_NUM_COILS - 2
      LineXY(CoilX,CoilY,CoilX + SPRING_COIL_WIDTH,CoilY + 0.5 * CoilHeight,RGB(0,255,0))
      LineXY(CoilX + SPRING_COIL_WIDTH,CoilY + 0.5 * CoilHeight,CoilX,CoilY + CoilHeight,RGB(0,255,0))
      CoilY = CoilY + CoilHeight
    Next
    
    LineXY(CoilX,CoilY,CoilX + 0.5 * SPRING_COIL_WIDTH,CoilY + 0.5 * CoilHeight,RGB(0,255,0))

    Box(#SCREEN_WIDTH / 2 - 16,0,32,16,RGB(0,0,255))
    Box(#SCREEN_WIDTH / 2 - 0.5 * SPRING_MASS,16 + SPRING_Y_POS + SPRING_REST_LENGTH,SPRING_MASS,SPRING_MASS,RGB(255,0,0))
    ;Box(#SCREEN_WIDTH / 2 - 1,16,2,SPRING_Y_POS + SPRING_REST_LENGTH,RGB(0,255,0))
    ;Define YPosText$ = StrD(SPRING_Y_POS + SPRING_REST_LENGTH)
    ;DrawText(#SCREEN_WIDTH / 2 - 0.5 * SPRING_MASS - TextWidth(YPosText$),SPRING_Y_POS + SPRING_REST_LENGTH,YPosText$,RGB(255,255,255))
EndProcedure

Procedure.l CreateTimer(Freq)
  Protected *Timer.Timer = AllocateMemory(SizeOf(Timer))
  
  If *Timer
    *Timer\Ticks = ElapsedMilliseconds()
    *Timer\Freq = Freq
  EndIf
  
  ProcedureReturn *Timer
  
EndProcedure

Procedure.l WaitTimer(*Timer.Timer)
  Protected Elapsed
  Protected Pings

  If *Timer
    Repeat
      Elapsed = ElapsedMilliseconds()
      If Elapsed - *Timer\Ticks > *Timer\Freq
        Break
      EndIf
      Delay(1)
    ForEver
  
    Pings = (Elapsed - *Timer\Ticks) / *Timer\Freq  
    *Timer\Ticks = *Timer\Ticks + *Timer\Freq * Pings
  
  EndIf

  ProcedureReturn Pings

EndProcedure

Procedure.l DestroyTimer(*Timer.Timer)
  If *Timer
    FreeMemory(*Timer)
    ProcedureReturn 1
  EndIf
  ProcedureReturn 0
EndProcedure


; the slope functions

Procedure.d SpringVelocity(v.d,y.d)
  ProcedureReturn v
EndProcedure

Procedure.d SpringAcceleration(v.d,y.d)
  ProcedureReturn (-(SPRING_DAMPING * v) - (SPRING_STIFFNESS * SPRING_Y_POS))  / SPRING_MASS
EndProcedure

; Runge-Kutta integrator

Procedure.d Integrate(x.d,y.d,h.d,f.DerivativeFunction)

  Protected u1.d
  Protected u2.d
  Protected u3.d
  Protected u4.d
  
  Protected y1.d
     
  u1 = h * f(x,y)
  u2 = h * f(x + h / 2.0, y + u1 / 2.0)
  u3 = h * f(x + h / 2.0, y + u2 / 2.0)
  u4 = h * f(x + h, y + u3)
  
  y1 = y + (1.0 / 6.0) * (u1 + (2.0 * u2) + (2.0 * u3) + u4)

  ProcedureReturn y1

EndProcedure

As an aside: well done and thanks to the PB team, I just installed Eeebuntu and PureBASIC on my Eee pc, sweet setup!

Posted: **Wed May 06, 2009 1:57 pm**

Looks good to me, nice work!

Posted: **Wed May 06, 2009 8:10 pm**

Heya,

Thanks!

I guess I just needed some reassurance as I've written some truly *awful* math algorithms in the past and when this thing appeared to work after the second attempt I was like "it can't be true!"

Posted: **Mon May 11, 2009 6:19 am**

Great post ProphetOfDoom. Eventhough you already posted the exact string of codes, it seems that I still do not understand on how it works.

Anyways, thanks a lot.

Simulation pret immobilier

Posted: **Fri May 15, 2009 2:12 pm**

Hey friends,

Thanks for the nice comments. If you want to know how it works I would suggest reading Calculus For Dummies (not that I'm calling you dummies lol - I just think it's a great book!) and then check out this website:

http://www.myphysicslab.com/

I'm still a beginner myself tho so don't ask me any tricky questions lol.

Posted: **Tue May 26, 2009 4:04 am**

I don't know if i'll ever have a use for it, but it's still neat to look at.

FYI Works in linux as well.

PureBasic Forums - English

Spring Simulation

Spring Simulation