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
