' This is a FreeBasic program. Created March 2, 2008 by Gabriel LaFreniere. Updated March 31, 2008.
' Sub routine procedures are listed in alphabetical order.
' Please download the editor/compiler from: http://fbide.freebasic.net

Dim As Ulongint iteration                                             'huge number for long period of time.
algorithm = 1: lambda = 200: work.page = 1: radius = 30
Dim sphere(-radius To radius, -radius To radius)
Dim As Double pi = 4 * Atn(1), omega, step.1, step.2
Dim As Double sine, cosine, motion, energy, potential, memory.1, memory.2
Screen 20,24,3: Gosub Initialization


' ********************************************************************
' MAIN LOOP.
' ********************************************************************

Do
  Swap work.page, visible.page
  Screenset work.page, visible.page
  Pcopy 2, work.page
  
  For j = 0 To skipped.frames

    
' ********************************************************************
' MY OWN SIMPLIFICATION OF PHILIPPE DELMOTTE'S VIRTUAL WAVE ALGORITHM USING NEWTON'S LAWS AND VERLET'S ALGORITHM.
' ********************************************************************

    motion = motion - energy                                          'I was unaware of the link to wavelength in 2005.
    energy = energy +  motion * step.1 ^ 2                            'step squared just once, finally equivalent.
                                          
    iteration += 1                                                    'to be converted in time for accurate sine curve.
  Next
  
  If algorithm = 1 Then
    For j = 0 To skipped.frames


' ********************************************************************
' ANSELME DEWAVRIN'S OCT. 2006 SIMPLIFIED ALGORITHM USING EULER'S METHOD.
' ********************************************************************

      sine  =  sine - cosine * step.1                                 'two simple program lines only. AMAZING!
      cosine = cosine + sine * step.1

    Next
  
  Else
    For j = 0 To skipped.frames


' ********************************************************************
' ANSELME DEWAVRIN'S DEC. 2006 SIMPLIFICATION OF JOCELYN MARCOTTE'S VIRTUAL WAVE ALGORITHM.
' ACCORDING TO DEWAVRIN, THIS ALGORITHM IS RELATED TO IIR (infinite impulse response) NUMERIC FILTER. 
' ********************************************************************

      potential = step.2 * memory.1 - memory.2                        'potential as sine, kinetic energy recoverable.
      memory.2 = memory.1                                             'memorizing potential's two previous states.
      memory.1 = potential

    Next
  End If


' ********************************************************************
' DISPLAYING OSCILLATING SPHERE.
' ********************************************************************

  If algorithm = 1 Then
    y.sphere = y.center - amplitude * sine                            'sphere oscillation according to sine.
    Circle(x.center + lambda / 4,y.center - amplitude * cosine),10,red'red circle for cosine, Euler's method only.
  Else
    y.sphere = y.center - amplitude * potential / abs.lambda          'sphere oscillation according to potential.
  End If

  For x = -radius To radius
    For y = -radius To radius
      If sphere(x,y) Then Pset(x + x.center, y + y.sphere), sphere(x,y)
    Next
  Next

  Circle(x.center, y.sphere), radius, gray                            'sphere contour.
  Circle(x.center, y.sphere),      2, black                           'small black circle to check accuracy.
  Circle(x.center + lambda / 4, y.center - amplitude * motion), 15, 0 'black circle for motion (indicates quadrature).
  Circle(x.center + lambda / 4, y.center - amplitude * motion),  2, 0 'small black circle to check accuracy.


' ********************************************************************
' DISPLAYING ACCURATE PROGRESSIVE SINE CURVE.
' ********************************************************************

  For x = 0 To 1023                                                   'horizontal scan (screen res. 1024 x 768 pixels).
    distance = (x.center - x)
    omega =         2 * pi * distance  /     lambda                   'pulsation in radians according to distance.
    omega = omega + 2 * pi * iteration / abs.lambda                   'adding the time in radians according to iteration.
    y_previous = y: y = y.center + amplitude * Sin(omega)             'vertical scan according to amplitude.
    Line(x - 1, y_previous)-(x, y), black                             'compare Euler's method to sine curve.
  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                                         'chr(27) Escape key or Windows' X quit button.
      Case "A":  algorithm = 1
      Case "B":  algorithm = 2
      Case "C":  If correction Then correction = 0 Else correction = 1
      Case "I":  lambda = 200: multiple = 0: algorithm = 1            'initialization.
                 correction = 0: skipped.frames = 0
      Case "M":  key$ = "": Run "WaveMechanics00.exe"                 'main menu.
      Case "P":  Screenset visible.page 
                 Color red: Locate 42, 89: Print "P - Paused. Press any key to resume.   "
                 key$ = "": Sleep                                     'pause.
                 Color black: Screenset work.page, visible.page
      Case "+":  multiple = multiple * 2                              'increase skipped frames.
                 If multiple = 0 Then multiple = 1
                 If multiple > 4096 Then multiple = 4096
                 skipped.frames = multiple * abs.lambda - 1
      Case "-":  multiple = multiple / 2                              'reduce skipped frames.
                 skipped.frames = multiple * abs.lambda - 1
                 If multiple < 1 Then multiple = 0: skipped.frames = 0
      Case "=":  If skipped.frames = 0 Then                           'add or remove skipped frames.
                   skipped.frames = 1: multiple = 4096
                   skipped.frames = multiple * abs.lambda - 1
                 Else
                   multiple = 0: skipped.frames = 0
                 End If
      Case "M+": skipped.frames += 1                                  'fine tuning (righ/left arrow keys).
                 key$ = "": Screenset 2, 2: Color green.text
                 Locate 43, 89: Print "Press +/-/= to skip frames:"; skipped.frames
                 Screenset work.page, visible.page
      Case "K+": skipped.frames -= 1
                 If skipped.frames < 0 Then skipped.frames = 0
                 key$ = "": Screenset 2, 2: Color green.text
                 Locate 43, 89: Print "Press +/-/= to skip frames:"; skipped.frames; "      "
                 Screenset work.page, visible.page
      Case "0":                                                       'reset.
      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 37: Locate 37, 88                                            'Delmotte's algorithm.
             If algorithm = 2 Then
               Print line37$
               If click Then algorithm = 1: Gosub Initialization
             End If               
    Case 38: Locate 38, 88                                            'Delmotte's algorithm.
             If algorithm = 1 Then
               Print line38$
               If click Then algorithm = 2: Gosub Initialization
             End If               
    Case 39: Locate 39, 88                                            'add correction.
             If correction Then
               Print line39b$
             Else
               Print line39a$
             End If
             If click Then
               If correction Then correction = 0 Else correction = 1
               Do: Getmouse a,b,c, click: Loop While click            'avoid unnecessary initialization. 
               Gosub Initialization
             End If
    Case 43: Locate 43, 88                                            'add or remove skipped frames.
             Print line43$;
             Locate 43,116: Print skipped.frames
             If click Then
               If skipped.frames = 0 Then
                 skipped.frames = 1: multiple = 4096
                 skipped.frames = multiple * abs.lambda - 1
               Else
                 multiple = 0: skipped.frames = 0
               End If
               Do: Getmouse a,b,c, click: Loop While click
               Gosub Initialization
             End If
    Case 45: Locate 45, 88                                            'reset.
             Print line45$
             If click Then
               Do: Getmouse a,b,c, click: Loop While click
               Gosub Initialization
             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
                                correction = 0: skipped.frames = 0    'initialization.
                                lambda = 200: algorithm = 1: 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 "WaveMechanics02.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 "WaveMechanics00.exe" 'previous program.
             End Select
  End Select


' ********************************************************************
' WAVELENGTH SELECTION.
' ********************************************************************

  If y.mouse < y.center + 18 And y.mouse > y.center - 18 Then
    Line(0, y.center - 16)-(1023, y.center + 16), black, b
    Line(0, y.center - 15)-(1023, y.center + 15), white, bf
    Line(x.center, y.center - 15)-(x.center + lambda / 4, y.center + 15), gray, bf
    Line(x.center, y.center - 16)-(x.center + lambda, y.center + 16), black, b
    Color black, white
    If lambda > 0 Then Locate 18, 44 Else Locate 18, 67
    Print "Lambda: ";: Print lambda; " pixels."
    new = (x.mouse - x.center) * 4: If Abs(new) < 24 Then new = 24
    If click Then                                                     'on click.
      lambda = new
      multiple = 0: skipped.frames = 0
      Gosub Initialization
    End If
    If wheel = -1 Then previous.wheel = wheel
    If wheel - previous.wheel Then                                    'on mouse wheel.
      lambda = lambda + wheel - previous.wheel
      If lambda > 0 And lambda < 24 Then lambda = 24
      If lambda < 0 And lambda >-24 Then lambda =-24
      multiple = 0: skipped.frames = 0
      Gosub Initialization
    End If
  End If

  Color black, background
  Locate 42, 113: Print iteration                                     'iteration number up to billions.
  if algorithm = 1 then
    locate 32, 118: print using "#.######"; sqr(sine^2+cosine^2)
  end if
  If wheel = -1 Then Else previous.wheel = wheel
Loop


'*********************************************************************
' FRIENDLY ADJUSTABLE FRAMES.
' ********************************************************************

Frame:
margin = 14                                                           'internal margin.
x.left  = 8 * left. - 8 - margin                                      'final dimensions. Lateral margins are larger.
x.right = 8 * left. - 8 + 8 * x.text + margin
y.top = top * 16 - 16 - .6 * margin
y.bottom = y.top + 16 * y.text + 2 * .5 * margin

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


'*********************************************************************
' GRAPHICS.
' ********************************************************************

Graphics:
For x = -radius To radius                                             'memorizing a matrix sphere.
  x1.squared = x ^ 2
  x2.squared = (x - 1000 + radius) ^ 2
  For y = -radius To radius
    potential.squared = y ^ 2
    memory.2.squared = (y - 1000 + radius) ^ 2
    If Sqr(x1.squared + potential.squared) <= radius Then
      diagonal = Sqr(x2.squared + memory.2.squared)
      shade = 4.5 * (diagonal - 1414 + 2 * radius + 10)               'gray shades for the sphere.
      If shade < 0 Then shade = 0 Else If shade > 255 Then shade = 255
      sphere(x,y) = Rgb(shade,shade,shade)
      If sphere(x,y) = 0 Then sphere(x,y) = 1                         '1 as printable black.
    End If
  Next
Next

For x = 0 To 512 Step abs.lambda / 4
  Line(x.center + x, y.center - 10)-(x.center + x, y.center), gray    'lambda / 4 marks.
  Line(x.center - x, y.center - 10)-(x.center - x, y.center), gray
Next

For x = 0 To 512 Step abs.lambda
  Line(x.center + x, y.center - 16)-(x.center + x, y.center), black   'lambda marks.
  Line(x.center - x, y.center - 16)-(x.center - x, y.center), black
Next

Line(0, y.center)-(1023, y.center), gray
Line(x.center + lambda / 4, y.center - amplitude)-(x.center + lambda / 4, y.center + amplitude), gray
Line(-1, y.center - amplitude)-(1024, y.center + amplitude), white, b
Line(x.center, y.center - amplitude)-(x.center + lambda, y.center + amplitude), gray, b
Return


'*********************************************************************
' INITIALIZATION.
' ********************************************************************

Initialization:
Windowtitle "WaveMechanics01. Oscillations. Anselme Dewavrin's algorithms. Philippe Delmotte and Jocelyn Marcotte's virtual wave medium."
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)
white = Rgb(255,255,255)
background =  Rgb(225,225,225)
light.green = Rgb(175,255,175)
green.text =  Rgb(0,100,0)
dark.gray =   Rgb(75,75,75)
abs.lambda = Abs(lambda)
If correction Then
  step.1 = Sqr(2 - (Sin(4 * pi / lambda) / Sin(2 * pi / lambda)))     'my adaptation of Dewarin's step (march 2008).
  step.2 = Sin(4 * pi / lambda) / Sin(2 * pi / lambda)                'Dewavrin's original step for IIR algorithm.
Else
  step.1 = 2 * pi / abs.lambda                                        'basic step in radians for Delmotte's algorithm.
  step.2 = 2 - (2 * pi / lambda) ^ 2                                  'basic step for Jocelyn Marcotte's algorithm.
End If
memory.1 = 0
memory.2 = 2 * pi
amplitude = 120
y.center = 279
x.center = 512
iteration = 0
sine   = 0
cosine = 1
energy = 0
motion = 1
energy = energy + step.1 ^ 2 * motion                                 'needs initial synchronization.
Screenset 2, visible.page: Color black, background: Cls
Gosub Title
Gosub Graphics
Gosub Text
Return


'*********************************************************************
' TEXT SECTION.
' ********************************************************************

Text:
x.text = Len("energy = energy + motion / step ^ 2 ")                  'text width (pixels = x * 8).
y.text = 2                                                            'number of lines (pixels = y * 16).
top = 29                                                              'upper limit.
left.= 89                                                             'limit on the left hand side: "Locate top, left".
Locate top, left.: ? "motion = motion - energy"
Locate    , left.: ? "energy = energy + motion * step ^ 2"
Gosub Frame

If algorithm = 1 Then
  top = 34
  Locate top, left.: ? "sine  =  sine - cosine * step"                'Euler's method (Delmotte).
  Locate    , left.: ? "cosine = cosine + sine * step"
  Locate 40,114: Print Using "###.#########"; step.1
  a$ = "step = sqr(2 - (sin(4 * pi / lambda) / sin(2 * pi / lambda)))"
  locate 32, 89: print "sqr(sine ^ 2 + cosine ^ 2) ="

  If correction Then
    Locate 40, 89: Print "Step using formula below: ";
  Else
    Locate 40, 89: Print "step =  2 * pi / lambda = ";
  End If
Else
  y.text = 3: top = 33
      potential = step.2 * memory.1 - memory.2                        'potential energy as sine, so kinetic recoverable.
      memory.2 = memory.1                                             'memorizing potential's two previous states.
      memory.1 = potential
  Locate top, left.: ? "potential = memory1 * step - memory2"         'infinite impulse filter (Marcotte).
  Locate    , left.: ? "memory2 = memory1"
  Locate    , left.: ? "memory1 = potential"
  Locate 40,114: Print Using "###.#########"; step.2
  a$ = "step = Sin(4 * pi / lambda) / Sin(2 * pi / lambda)"
  If correction Then
    Locate 40, 89: Print "Step using formula below: ";
  Else
    Locate 40, 89: Print "step = 2-(2*pi/lambda)^2 =";
  End If
End If  
Gosub Frame
x.text = Len(a$): top = 45: left.= 13: y.text = 1
Locate top, left.: ? a$
Gosub Frame

Locate 4, 3
?"Robert Hooke wrote in 1678: ® Ut tensio sic vis ¯ (As the extension, so the force). Any mechanical oscillating system works ":Locate,3
?"like a pendulum. Firstly, motion as kinetic energy is stored as energy inside a field of force. Then this stored energy,    ":Locate,3
?"which exerts a force like a spring, is converted again into motion. Because energy cannot be created or destroyed, this     ":Locate,3
?"process may go on endlessly on condition that force and extension are proportional according to Hooke's law. Otherwise,     ":Locate,3
?"energy is progressively transferred into heat or waves such as sound. The important point is that energy transmission in    ":Locate,3
?"a step-by-step process introduces quantum properties."
Locate 28, 3
?"Mr. Philippe Delmotte from Lille, France, invented a computerized virtual medium   ": Locate, 3
?"for waves in June 2005. I simplified Mr. Delmotte's algorithm in November 2005 in  ": Locate, 3
?"order to obtain pure oscillations without waves. The recent algorithm is:"
?: Locate, 3
If algorithm = 1 Then
  ?"Mr. Delmotte deduced his algorithm from Verlet's one. He also took Newton's laws ": Locate, 3
  ?"into account. However, Mr. Anselme Dewavrin also from Lille, France, informed me ": Locate, 3
  ?"in October 2006 that my iterative calculus was clearly similar to Euler's method ": Locate, 3
  ?"using sine and cosine. As far as I know, he truly invented this simplification:  "
Else
  ?"Mr. Jocelyn Marcotte invented his own algorithm for computerized virtual waves in": Locate, 3
  ?"jan 2006. In Dec 2006, Mr. Dewavrin also discovered that Mr. Marcotte's algorithm": Locate, 3
  ?"was similar to the IIR (infinite impulse response) electronic filter. Once again,": Locate, 3
  ?"he succeeded in simplifying the algorithm for oscillations and he obtained this: "
End If
Locate 37, 3
?"The sphere above follows Mr. Delmotte's sine or Mr. Marcotte's potential. Amazingly,":Locate, 3
?"all three algorithms introduce exactly the same quantum anomaly. It becomes clearly": Locate, 3
?"visible after millions of iterations, especially for shorter wavelengths. In order ": Locate, 3
?"to accelerate the process, you may activate the oscilloscope by pressing the + key.": Locate, 3
?"In addition, the sine and cosine (the red circle above) are not perfectly exact. It": Locate, 3
?"turns out that most mechanical waves truly behave this way because the medium must ": Locate, 3
?"transmit energy step by step. The accurate step (press C to activate) is given by:"

line37$  = " A - Philippe Delmotte's algorithm.     "
line38$  = " B - Jocelyn Marcotte's algorithm.      "
line39a$ = " C - Add step correction.               "
line39b$ = " C - Remove step correction.            "
line43$  = " Press +/-/= to skip frames:            "
line45$  = " 0 - Reset.                             "
line47a$ = " I - Initialize    "
line47b$ = " M - Menu    "
line47c$ = " Slow          "
line48a$ = " Previous Program  "
line48b$ = " Quit (Esc.) "
line48c$ = " Next Program  "

Locate 19, 82: Print "Click here to select wavelength: ";: Print lambda; " pixels."
Locate 42, 89: Print "P - Pause.   Iteration:"

Locate 44, 89: Print "Fine tuning: use left/right arrow keys."
Locate 45, 88: Print "Press 0 to reset."

Color green.text
Locate 37, 88: Print line37$
Locate 38, 88: Print line38$
Locate 39, 88: If correction Then Print line39b$ Else Print line39a$
Locate 43, 88: Print line43$;
Locate 43,116: Print skipped.frames
Locate 45, 88: Print line45$
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 "March 31, 2008. This program may be"
Locate 48, 89: Print "freely distributed, copied or modified.";
Color blue
If algorithm = 1 Then
  Locate 37, 88
  Print line37$
Else
  Locate 38, 88
  Print line38$
End If
Color black
Screenset work.page, visible.page
Return


'*********************************************************************
' ENLARGED TITLE USING MY ORIGINAL DEPIXELATION METHOD. 
'*********************************************************************

Title:
title$ = "Oscillations"
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 x = x.center - 8 * Len(title$) To x.center + 8 * Len(title$)      'adding light and shades.
  For y = 0 To 34
    If Point(x, y) = gold And Point(x-1, y-1) = background Then Pset(x-1, y-1), white
    If Point(x, y) = gold And Point(x+1, y+1) = background Then Pset(x+1, y+1), black
  Next
Next
Return