'Grid Cell Fields - Gary Gaulin, 2011
'From "An oscillatory interference model of grid cell firing", N Burgess, C Barry, and J O'Keefe
'   http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2678278/
'Navigational memory addressing of mammal and probably insect brains
' use neurons as voltage controlled oscillators.
'Speed and direction information combine with ~8-10 Hz Theta Waves
'  to produce hexagonal grid interference patterns of varying sizes
'  that sum to address a single "Place" in time.
'Knowledge of AM/FM radio electronics and AC phase shifting is very helpful.
'No "Neural Network" or conventional AI experience is here required.
'Video primer, lecture starts 7 1/2 minutes in:
' http://www.youtube.com/watch?v=pATzOJLptGo
'Other good open access resources:
' http://www.icn.ucl.ac.uk/nburgess/GridCellModelPoster05.pdf
' http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2408670/
' http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2653064/

Option Explicit         'Tells program compiler to report when error is made coding variables.
Const PI = 3.14159265358979
Const PI2 = PI / 2
Const Radian = PI * 2
Const Degrees = Radian / 360
Const ThetaRhythm = 10 * 2 * PI

Dim EarlyExit As Boolean
Dim Speed As Double
Dim Sec As Double
Dim T As Double
Dim Tstep As Double
Dim H As Double     'Heaviside function  (H(x)=x if x>0; H(x)=0 otherwise)
Dim B As Double
Dim Vsoma As Double
Dim Vdend As Double
Dim Vdend1 As Double
Dim Vdend2 As Double
Dim Vdend3 As Double
Dim Vdend4 As Double
Dim Vdend5 As Double
Dim Vdend6 As Double
Dim Vsum As Double
Dim Vsum1 As Double
Dim Vsum2 As Double
Dim Vsum3 As Double
Dim Vsum4 As Double
Dim Vsum5 As Double
Dim Vsum6 As Double
Dim Wsoma As Double
Dim Wdend As Double
Dim Wdend1 As Double
Dim Wdend2 As Double
Dim Wdend3 As Double
Dim Wdend4 As Double
Dim Wdend5 As Double
Dim Wdend6 As Double
Dim Adend As Double
Dim Asoma As Double
Dim Pdend As Double
Dim HeadingDirection As Double
Dim Ysoma As Double
Dim Ydend As Double
Dim Ysum As Double
Dim WaveHeight As Long

Dim Pic
Dim Xpixels As Long
Dim Ypixels As Long
Dim Xctr As Long
Dim Yctr As Long
Dim X As Long
Dim Y As Long
Dim Xp As Long
Dim Yp As Long
Dim C As Long
Dim Rd As Long
Dim Gr As Long
Dim Bl As Long
Dim A As Double

Private Sub Form_Load()
   Call PhasePrecessionChart
End Sub

Private Sub Form_Unload(Cancel As Integer)
  EarlyExit = True      'End program without anything left running or in RAM.
  DoEvents
End Sub

Private Sub GridCellFieldsCommand_Click()
   Call GridCellFields
End Sub

Private Sub SaveBitmapsCommand_Click()
   SavePicture FieldsPicture.Image, App.Path & "\GridCellFields.BMP"
   SavePicture PhasePicture.Image, App.Path & "\PhasePrecessionChart.BMP"
End Sub

Private Function PrintPixel(XS As Long, YS As Long, Num As Double)
If Num = 0 Then
    Rd = 0
    Gr = 0
    Bl = 120
  GoTo PrintIt
End If
If Num < 0.25 Then
A = Num * 2 * PI
    Rd = 0
    Gr = Sin(A) * 255
    Bl = 129 + (Sin(A) * 256)
  GoTo PrintIt
End If
If Num <= 0.5 Then
A = (Num - 0.25) * 2 * PI
    Rd = 0
    Gr = 255
    Bl = (Cos(A)) * 255
  GoTo PrintIt
End If
If Num <= 0.7 Then
A = (Num - 0.5) * 2.5 * PI
    Bl = 0
    Gr = 255
    Rd = Sin(A) * 255
  GoTo PrintIt
End If
A = (Num - 0.7) * 3.3334 * PI2
    Bl = 0
    Gr = Cos(A) * 255
    Rd = 255
PrintIt:
  Pic.PSet (XS, YS), RGB(Rd, Gr, Bl)
End Function

Private Sub PhasePrecessionChart()
  Set Pic = PhasePicture
    Xpixels = Pic.Width
    Ypixels = Pic.Height
    Ysum = Pic.Height * 0.25
    Ydend = Pic.Height * 0.615
    Ysoma = Pic.Height * 0.87
    Asoma = 1
    Adend = 1
    Pdend = 0
    Wsoma = ThetaRhythm
    Wdend = 11.5 * 2 * PI
    WaveHeight = (Pic.Height / 8.2)
For T = 0 To 2 Step 0.0002
    Xp = T * Pic.Width / 2
  Vsoma = Asoma * Sin(Wsoma * T)
    Yp = -Vsoma * WaveHeight
    Pic.PSet (Xp, Yp + Ysoma), vbRed
  Vdend = Adend * Sin(Wdend * (T + Pdend))
    Yp = -Vdend * WaveHeight
    Pic.PSet (Xp, Yp + Ydend), vbBlue
  Vsum = Vsoma + Vdend
    Yp = -Vsum * WaveHeight
    Pic.PSet (Xp, Yp + Ysum), vbBlack
Next T
End Sub

Private Sub GridCellFields()
  Set Pic = FieldsPicture
      Pic.Cls
      Xpixels = Pic.Width
      Xctr = Xpixels / 2
      Ypixels = Pic.Height
      Yctr = Ypixels / 2
                                        Pdend = 0
                                      Sec = 78 / 30
                              Tstep = Sec / Xctr
                           Speed = 30
                       B = 0.05 * 2 * PI
              Wsoma = 10 * 2 * PI
              Asoma = 0.5
              Adend = 0.5
For HeadingDirection = 0 To Radian Step 0.0008
     Wdend1 = Wsoma + (B * Speed * Cos(HeadingDirection - 0))
     Wdend2 = Wsoma + (B * Speed * Cos(HeadingDirection - (60 * Degrees)))
     Wdend3 = Wsoma + (B * Speed * Cos(HeadingDirection - (120 * Degrees)))
     Wdend4 = Wsoma + (B * Speed * Cos(HeadingDirection - (180 * Degrees)))
     Wdend5 = Wsoma + (B * Speed * Cos(HeadingDirection - (240 * Degrees)))
     Wdend6 = Wsoma + (B * Speed * Cos(HeadingDirection - (300 * Degrees)))
  For T = 0 To Sec + 1.1 Step Tstep
      Vsoma = Asoma * Cos(Wsoma * T)
     Vdend1 = Adend * Cos(Wdend1 * (T + Pdend))
     Vdend2 = Adend * Cos(Wdend2 * (T + Pdend))
     Vdend3 = Adend * Cos(Wdend3 * (T + Pdend))
     Vdend4 = Adend * Cos(Wdend4 * (T + Pdend))
     Vdend5 = Adend * Cos(Wdend5 * (T + Pdend))
     Vdend6 = Adend * Cos(Wdend6 * (T + Pdend))
           Vsum1 = Vsoma + Vdend1
           Vsum2 = Vsoma + Vdend2
           Vsum3 = Vsoma + Vdend3
           Vsum4 = Vsoma + Vdend4
           Vsum5 = Vsoma + Vdend5
           Vsum6 = Vsoma + Vdend6
       H = Vsum1 * Vsum2 * Vsum3 * Vsum4 * Vsum5 * Vsum6
    If H < 0 Then H = 0
                     Xp = Cos(HeadingDirection) * T * Xctr / Sec
                     Yp = Sin(HeadingDirection) * T * Xctr / Sec
     Call PrintPixel(Xp + Xctr, Yp + Yctr, H)
  Next T
      DoEvents
   If EarlyExit = True Then Exit Sub
Next HeadingDirection
End Sub