' This is a FreeBasic program. Created Nov. 2005 by Gabriel LaFreniere. Updated April 30, 2008. ' Please download the editor/compiler from: http://fbide.freebasic.net ' Sub routine procedures are listed in alphabetical order. granules = 127 '128 granules including granule 0. Dim As Double force(granules+1), slope(granules+1), influence(granules+1) Dim As Double pi = 4 * Atn(1), lambda, phase, Gaussian, speed, ratio, x, y choice$ = "A": lambda = 16: speed = 1: unidirectional = 0: work.page = 1 Screen 20,24,3: Gosub Initialization '********************************************************************* ' MAIN LOOP. '********************************************************************* Do Swap work.page, visible.page Screenset work.page, visible.page Pcopy matrix.page, work.page Screensync 'produces about 85 iterations per second. force(0) = force(1) 'soft reflection (left). 'automatic hard reflection (right). influence(1) = force(0) + force(2) - 2 * force(1) 'influence on granule 1, one iteration in advance. For x = 1 To granules '********************************************************************* ' THIS IS MY TRANSPOSITION OF PHILIPPE DELMOTTE'S VIRTUAL WAVE ' ALGORITHM (JUNE 2005) USING NEWTON'S LAWS AND VERLET'S ALGORITHM. '********************************************************************* influence(x+1) = force(x) + force(x+2) - 2 * force(x+1) 'influence based on force difference. slope(x) = slope(x) + influence(x) * speed 'force difference. Cosine curve for sine waves. force(x) = force(x) + slope(x) 'wave curve. Sine curve for sine waves. '********************************************************************* ' END OF ALGORITHM. - THREE PROGRAM LINES ONLY. '********************************************************************* Circle(x * 8, y.center - slope(x) / ratio), 2.5, green Paint (x * 8, y.center - slope(x) / ratio), green, green Circle(x * 8, y.center - force(x)), 2.5, black Paint (x * 8, y.center - force(x)), black, black Next '********************************************************************* ' KEYBOARD MANAGEMENT. '********************************************************************* key$ = Inkey If Len(key$) Then If Len(key$) = 2 Then key$ = Right(key$, 1) + "+" Else key$ = Ucase(key$) Select Case key$ Case Chr(27), "k+": End 'Escape key or Windows' X quit button. Case "A": choice$ = "A" 'four choices. Case "B": choice$ = "B" Case "C": choice$ = "C" Case "D": choice$ = "D" Case "E": unidirectional = 1 'one-way impulse. Case "F": unidirectional = 0 'bidirectional impulse. Case "G": If speed < 1 Then 'full speed. key$ = "" speed = 1: Screenset matrix.page Color blue, background: Locate 35, 88: Print line35$ Color green.text: Locate 36, 88: Print line36$ End If Case "H": If speed = 1 Then 'slower speed. key$ = "" speed = .9: Screenset matrix.page Color blue, background: Locate 36, 88: Print line36$ Color green.text: Locate 35, 88: Print line35$ End If Case "I": choice$ = "A": speed = 1: unidirectional = 0 'initialization. Case "M": key$ = "": Run "WaveMechanics00.exe" 'main menu. Case "P": Screenset visible.page 'pause. Color red, background: Locate 45, 89 Print "P - Paused. Press any key to resume. " key$ = "": Sleep Color black: Screenset work.page, visible.page Case "R": 'reset via initialization. Case Else: key$ = "" 'avoid initialization. End Select If Len(key$) Then Gosub Initialization Do: Loop While Len(Inkey) 'clear buffer. End If '********************************************************************* ' MOUSE MANAGEMENT. '********************************************************************* Getmouse x.mouse, y.mouse, wheel, click line.number = .5 + y.mouse / 16 If line.number < 46 And x.mouse < 695 Then line.number = 0 Color green.text, white Select Case line.number Case 27: Locate 27, 88 If choice$ <> "A" Then Print line27$ If click Then choice$ = "A": Gosub Initialization End If Case 28: Locate 28, 88 If choice$ <> "B" Then Print line28$ If click Then choice$ = "B": Gosub Initialization End If Case 29: Locate 29, 88 If choice$ <> "C" Then Print line29$ If click Then choice$ = "C": Gosub Initialization End If Case 30: Locate 30, 88 If choice$ <> "D" Then Print line30$ If click Then choice$ = "D": Gosub Initialization End If Case 32: Locate 32, 88 If unidirectional = 0 Then Print line32$ If click Then unidirectional = 1: Gosub Initialization End If Case 33: Locate 33, 88 If unidirectional = 1 Then Print line33$ If click Then unidirectional = 0: Gosub Initialization End If Case 35: Locate 35, 88 If speed < 1 Then Print line35$ If click Then speed = 1: Screenset matrix.page Color blue, background: Locate 35, 88: Print line35$ Color green.text: Locate 36, 88: Print line36$ End If End If Case 36: Locate 36, 88 If speed = 1 Then Print line36$ If click Then speed = .9: Screenset matrix.page Color blue, background: Locate 36, 88: Print line36$ Color green.text: Locate 35, 88: Print line35$ End If End If Case 47: Select Case x.mouse Case Is > 700 Case Is > 576: Locate 47, 73: Print line47c$: Sleep 200'slow. If click Then Sleep 1000 'slower. Case Is > 472: Locate 47, 60: Print line47b$ If click Then Run "WaveMechanics00.exe" 'main menu. Case Is > 318: Locate 47, 41: Print line47a$ If click Then 'initialization. choice$ = "A": speed = 1: unidirectional = 0 Do: Getmouse a,b,c, click: Loop While click Gosub Initialization End If End Select Case 48: Select Case x.mouse Case Is > 700 Case Is > 576: Locate 48, 73: Print line48c$; If click Then Run "WaveMechanics03.exe" 'next program. Case Is > 472: Locate 48, 60: Print line48b$; If click Then End 'quit. Case Is > 318: Locate 48, 41: Print line48a$; If click Then Run "WaveMechanics01.exe" 'previous program. End Select End Select Loop '********************************************************************* ' CHOICE A - GAUSSIAN IMPULSE AND RESULTING RICKER WAVELET. '********************************************************************* Choice.A: If unidirectional Then amplitude = 3 * lambda Else amplitude = 6 * lambda For j = 0 To granules + 1 x = j - granules / 2 Gaussian = 1.008 ^ (-x ^ 2) force(j) = amplitude * Gaussian slope(j) = 0 'no slope for standing waves (bidirectional). ' slope(j) = amplitude * Gaussian * x / 63 'check Jocelyn Marcotte's interesting formula. Next If unidirectional Then 'skip this for standing waves (slope = 0). For j = 1 To granules slope(j) = Sqr(speed) * (force(j) - force(j+1)) 'rightward. ' slope(j) = sqr(speed) * (force(j) - force(j-1)) 'leftward also possible. Next End If Return '********************************************************************* ' CHOICE B - GAUSSIAN-DAMPED SINUSOIDAL IMPULSE. '********************************************************************* Choice.B: If unidirectional Then amplitude = 2 * lambda Else amplitude = 4 * lambda For j = 0 To granules + 1 x = granules / 2 - j phase = 2 * pi * x / lambda Gaussian = 1.002 ^ (-x ^ 2) force(j) = amplitude * Gaussian * Sin(phase) slope(j) = 0 'no slope for standing waves (bidirectional). Next If unidirectional Then For j = 1 To granules slope(j) = Sqr(speed) * (force(j) - force(j+1)) Next End If Return '********************************************************************* ' CHOICE C - SAWTOOTH WAVES. '********************************************************************* Choice.C: For j = 0 To granules + 1 force(j) = 0 Next If unidirectional Then amplitude = lambda Else amplitude = 2 * lambda y.point = amplitude For j = granules / 2 - 2 * lambda To granules / 2 + 2 * lambda If y.point < -amplitude Then y.point = amplitude force(j) = y.point y.point -= amplitude / 8 Next For j = 1 To granules 'sawtooth and square waves need more sophisticated If unidirectional Then 'treatment for standing waves (bidirectional) slope(j) = Sqr(speed) * (force(j) - force(j+1)) 'probably because they contain harmonics. Else slope(j) = Sqr(speed) * (force(j) - force(j+1) / 2 - force(j-1) / 2) End If Next Return '********************************************************************* ' CHOICE D - SQUARE WAVES. '********************************************************************* Choice.D: If unidirectional Then amplitude = lambda Else amplitude = 2 * lambda For j = 0 To granules + 1 force(j) = 0: slope(j) = 0 Next For j = granules / 2 - 2 * lambda To granules / 2 + lambda Step lambda For k = 1 To lambda / 2 force(j+k) = amplitude force(j+k+lambda / 2) = -amplitude Next Next For j = 1 To granules If unidirectional Then slope(j) = Sqr(speed) * (force(j) - force(j+1)) Else slope(j) = Sqr(speed) * (force(j) - force(j+1) / 2 - force(j-1) / 2) End If Next Return '********************************************************************* ' FRIENDLY ADJUSTABLE FRAMES. '********************************************************************* Frame: margin = 10 x.left = 8 * left. - margin - 4 x.right = 8 * left. + 8 * x.text + margin - 6 y.top = top * 16 - margin - 16 y.bottom = y.top + 16 * y.text + 2 * margin - 2 Line (x.left, y.top)-(x.right + 1, y.bottom + 1), gray, B Line (x.left + 1, y.top + 1)-(x.right, y.bottom), black, B Line (x.left + 1, y.bottom)-(x.right, y.bottom), white Line (x.left, y.bottom + 1)-(x.right + 1, y.bottom + 1), white Line (x.right, y.top + 1)-(x.right, y.bottom), white Line (x.right+ 1, y.top)-(x.right + 1, y.bottom + 1), white Return '********************************************************************* ' INITIALIZATION. '********************************************************************* Initialization: Windowtitle "WaveMechanics02 - Waves on a string (one dimension)." matrix.page = 2 ratio = Sqr(speed) * 2 * pi / lambda x.center = 512 y.center = 180 red = Rgb(255,0,0) blue = Rgb(0,0,255) green = Rgb(0,200,0) cyan = Rgb(0,100,100) gold = Rgb(180,150,100) gray = Rgb(125,125,125) buff = Rgb(255,255,200) white = Rgb(255,255,255) background = Rgb(225,225,225) light.green = Rgb(175,255,175) green.text = Rgb(0,125,0) dark.gray = Rgb(75,75,75) Screenset matrix.page, matrix.page Color black, background: Cls Gosub Title Gosub Text Select Case choice$ Case "A": Gosub Choice.A Case "B": Gosub Choice.B Case "C": Gosub Choice.C Case "D": Gosub Choice.D End Select Circle(0, y.center), 2.5, black Paint(0, y.center), black, black Circle(1024, y.center), 2.5, black Paint(1023, y.center), black, black Line(0,y.center)-(1023, y.center), gray Return '********************************************************************* ' TEXT. '********************************************************************* Text: Locate 1, 5 : Print "Left: soft reflection." Locate 1,102: Print "Right: hard reflection." line27$ = " A - Gaussian 1-D Ricker wavelet. " line28$ = " B - Sine waves. " line29$ = " C - Sawtooth waves. " line30$ = " D - Square waves. " line32$ = " E - Unidirectional impulse. " line33$ = " F - Bidirectional impulse. " line35$ = " G - Full speed: perfect waves. " line36$ = " H - Slower speed: quantum effect. " line47a$ = " I - Initialize. " line47b$ = " M - Menu. " line47c$ = " Slow. " line48a$ = " Previous Program. " line48b$ = " Quit (Esc.)." line48c$ = " Next Program. " Locate 22 Locate, 3: Print "Mr. Philippe Delmotte elaborated his new virtual medium on June 2005. His simple" Locate, 3: Print "and truly amazing algorithm is capable of simulating all kind of waves. Although" Locate, 3: Print "many algorithms had already been proposed before 2005, especially for video" Locate, 3: Print "games, none of them could reproduce waves with such a perfection." Print Locate, 3: Print "Because waves on a string are involving only one dimension, each so-called aether" Locate, 3: Print "granule is influenced by its two nearest neighbors according to their force (or" Locate, 3: Print "extension) difference. So the algorithm's first line goes like this:" Print Locate, 3: Print " influence(x) = force(x-1) + force(x+1) - 2 * force(x)" Print Locate, 3: Print "Secondly, because this force or extension difference is actually the curve slope," Locate, 3: Print "influence modifies the slope this way:" Print Locate, 3: Print " slope(x) = slope(x) + influence(x)" Print Locate, 3: Print "Finally, the new force is modified according to the new slope. It is that simple:" Print Locate, 3: Print " force(x) = force(x) + slope(x)" Print Locate, 3: Print "The black curve above indicates the force or extension, according to Hookes' law," Locate, 3: Print "and the green one indicates the slope, which is equivalent to quadrature for sine" Locate, 3: Print "waves. This allows one to easily control Delmotte's algorithm, which was somewhat" Locate, 3: Print "different in its original form." Locate 45, 89: Print "P - Pause. R - Reset." Color green.text x.text = Len(line27$) - 1 'text width (pixels = x * 8). y.text = 10 'number of lines (pixels = y * 16). top = 27 'upper limit. left.= 88 'limit on the left hand side: "Locate top, left". Gosub frame Locate 27,88: Print line27$ Locate 28,88: Print line28$ Locate 29,88: Print line29$ Locate 30,88: Print line30$ Locate 32,88: Print line32$ Locate 33,88: Print line33$ Locate 35,88: Print line35$ Locate 36,88: Print line36$ Locate 47,42: Print "I - Initialize. M - Menu. Slow." Locate 48,42: Print "Previous Program. Quit (Esc.). Next Program."; Color gray Locate 47, 3: Print "Thanks to the creators of FreeBASIC." Locate 48, 3: Print "Gabriel LaFreniere glafreniere.com"; Locate 47,89: Print "April 30, 2008. This program may be" Locate 48,89: Print "freely distributed, copied or modified."; Color blue Select Case choice$ Case "A": Locate 27, 88: Print line27$ Case "B": Locate 28, 88: Print line28$ Case "C": Locate 29, 88: Print line29$ Case "D": Locate 30, 88: Print line30$ End Select If speed = 1 Then Locate 35, 88 Else Locate 36, 88 If speed < 1 Then Print line36$ Else Print line35$ If unidirectional Then Locate 32, 88 Else Locate 33, 88 If unidirectional Then Print line32$ Else Print line33$ Color black Return '********************************************************************* ' ENLARGED TITLE USING MY ORIGINAL DEPIXELATION METHOD. '********************************************************************* Title: title$ = "Waves on a string" Locate 1,1: Print title$ For x1 = 0 To 8 * Len(title$) x2 = x.center + 2 * x1 - 8 * Len(title$) For y1 = 1 To 14 If Point(x1, y1) = black Then y2 = 2 * y1 + 2 Line(x2, y2)-(x2, y2 + 1), gold If Point(x1 + 1, y1 - 1) = 0 Then Line(x2 + 1, y2 - 1)-(x2 + 1, y2 + 0), gold If Point(x1 + 1, y1 + 0) = 0 Then Line(x2 + 1, y2 + 0)-(x2 + 1, y2 + 1), gold If Point(x1 + 1, y1 + 1) = 0 Then Line(x2 + 1, y2 + 1)-(x2 + 1, y2 + 2), gold End If Next If (x1+1) Mod 8 Then Else Line(x2+1, 0)-(x2+1, 34), background 'separate invasive characters such as capital M. Next Line(0, 0)-(8 * Len(title$), 14), background, bf 'matrix title erased. For x1 = x.center - 8 * Len(title$) To x.center + 8 * Len(title$) 'adding light and shades. For y1 = 0 To 34 If Point(x1, y1) = gold And Point(x1-1, y1-1) = background Then Pset(x1-1, y1-1), buff If Point(x1, y1) = gold And Point(x1+1, y1+1) = background Then Pset(x1+1, y1+1), black Next Next For x = 512 - 8 * Len(title$) To 512 + 8 * Len(title$) 'adding luster. For y = 4 To 32 If Point(x, y) = gold Then darken = 9 * Abs(18 - y): Pset(x, y), Rgb(240 - darken, 200 - darken, 120 - darken) Next Next Return