VERSION 5.00
Begin VB.Form Form1 
   BorderStyle     =   1  'Fixed Single
   Caption         =   "Optics Simulation"
   ClientHeight    =   6630
   ClientLeft      =   45
   ClientTop       =   330
   ClientWidth     =   10395
   LinkTopic       =   "Form1"
   MaxButton       =   0   'False
   ScaleHeight     =   6630
   ScaleWidth      =   10395
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Reset"
      Height          =   255
      Left            =   4560
      TabIndex        =   14
      Top             =   5400
      Width           =   855
   End
   Begin VB.CheckBox Check3 
      Caption         =   "Pause Simulation"
      Height          =   255
      Left            =   5760
      TabIndex        =   13
      Top             =   5880
      Value           =   1  'Checked
      Width           =   1935
   End
   Begin VB.CheckBox Check2 
      Caption         =   "Show B field"
      Height          =   255
      Left            =   5760
      TabIndex        =   12
      Top             =   5520
      Width           =   1455
   End
   Begin VB.PictureBox Pic1 
      Appearance      =   0  'Flat
      BackColor       =   &H80000005&
      ForeColor       =   &H80000008&
      Height          =   4530
      Left            =   120
      ScaleHeight     =   4500
      ScaleWidth      =   4500
      TabIndex        =   9
      Top             =   120
      Width           =   4530
   End
   Begin VB.ComboBox DeltaTime 
      Height          =   315
      Left            =   2520
      Style           =   2  'Dropdown List
      TabIndex        =   7
      Top             =   6120
      Width           =   1575
   End
   Begin VB.ComboBox DeltaSpace 
      Height          =   315
      Left            =   2520
      Style           =   2  'Dropdown List
      TabIndex        =   5
      Top             =   5760
      Width           =   1815
   End
   Begin VB.CheckBox Check1 
      Caption         =   "Grid Lines"
      Height          =   255
      Left            =   5760
      TabIndex        =   4
      Top             =   5160
      Value           =   1  'Checked
      Width           =   1455
   End
   Begin VB.ComboBox WaveType 
      Height          =   315
      Left            =   1680
      Style           =   2  'Dropdown List
      TabIndex        =   3
      Top             =   5400
      Width           =   2775
   End
   Begin VB.PictureBox Pic2 
      Appearance      =   0  'Flat
      BackColor       =   &H80000005&
      ForeColor       =   &H80000008&
      Height          =   4530
      Left            =   4800
      ScaleHeight     =   4500
      ScaleWidth      =   4500
      TabIndex        =   0
      Top             =   120
      Width           =   4530
   End
   Begin VB.Timer Timer1 
      Interval        =   100
      Left            =   7920
      Top             =   5280
   End
   Begin VB.Label Label8 
      Alignment       =   2  'Center
      BackColor       =   &H000000FF&
      Caption         =   "E-Field -Y"
      ForeColor       =   &H00FFFFFF&
      Height          =   255
      Left            =   3360
      TabIndex        =   17
      Top             =   4920
      Width           =   1095
   End
   Begin VB.Label Label7 
      Alignment       =   2  'Center
      BackColor       =   &H00FF0000&
      Caption         =   "E-Field +Y"
      ForeColor       =   &H00FFFFFF&
      Height          =   255
      Left            =   3360
      TabIndex        =   16
      Top             =   4680
      Width           =   1095
   End
   Begin VB.Label Label6 
      Caption         =   "Initial Configuration"
      Height          =   255
      Left            =   120
      TabIndex        =   15
      Top             =   5400
      Width           =   1455
   End
   Begin VB.Label Label5 
      Caption         =   "DT"
      Height          =   255
      Left            =   2040
      TabIndex        =   11
      Top             =   6120
      Width           =   375
   End
   Begin VB.Label Label4 
      Caption         =   "DS"
      Height          =   255
      Left            =   2040
      TabIndex        =   10
      Top             =   5760
      Width           =   375
   End
   Begin VB.Label Label3 
      Caption         =   "1800 nm"
      Height          =   255
      Left            =   9600
      TabIndex        =   8
      Top             =   2160
      Width           =   735
   End
   Begin VB.Line Line1 
      X1              =   9480
      X2              =   9480
      Y1              =   120
      Y2              =   4680
   End
   Begin VB.Label TimeLabel 
      Caption         =   "Iteration:"
      Height          =   255
      Left            =   5760
      TabIndex        =   6
      Top             =   6240
      Width           =   2175
   End
   Begin VB.Label Label2 
      Caption         =   "S(X,Y,Z) : Magnitude of Poynting vector"
      Height          =   255
      Left            =   5280
      TabIndex        =   2
      Top             =   4800
      Width           =   3015
   End
   Begin VB.Label Label1 
      Caption         =   "EY(X,Y) : E-Field along the Y-axis"
      BeginProperty Font 
         Name            =   "Microsoft Sans Serif"
         Size            =   8.25
         Charset         =   0
         Weight          =   400
         Underline       =   0   'False
         Italic          =   0   'False
         Strikethrough   =   0   'False
      EndProperty
      Height          =   255
      Left            =   360
      TabIndex        =   1
      Top             =   4800
      Width           =   2655
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
'*** Optics Simulation - Dan Piehl ***
'
'This program is designed to simulation some optical phenomena caused by
'the propagation of electromagnetic waves.  A 3-dimensional array is used
'for the electric and magnetic components of each cell in the simulation
'space.
'
Option Explicit
Option Base 1
Option Compare Text

Const Max = 50          'Largest array subscript value
Dim LimitX As Integer   'Upper bound on X subscript
Dim LimitY As Integer   'Upper bound on Y subscript
Dim LimitZ As Integer   'Upper bound on Z subscript
Dim Center As Integer   'Center of simulation space
Const PI = 3.14159265358979 'Approximate Pi value

Dim TimeCounter As Long 'Iteration counter
Dim TimeClock As Double 'Elapsed time accumulator

Dim EX(0 To Max, 0 To Max, 0 To Max) As Double 'E vector field
Dim EY(0 To Max, 0 To Max, 0 To Max) As Double
Dim EZ(0 To Max, 0 To Max, 0 To Max) As Double

Dim BX(0 To Max, 0 To Max, 0 To Max) As Double 'B vector field
Dim BY(0 To Max, 0 To Max, 0 To Max) As Double
Dim BZ(0 To Max, 0 To Max, 0 To Max) As Double

Dim DeltaEX(0 To Max, 0 To Max, 0 To Max) As Double 'Changes to E vector field
Dim DeltaEY(0 To Max, 0 To Max, 0 To Max) As Double
Dim DeltaEZ(0 To Max, 0 To Max, 0 To Max) As Double

Dim DeltaBX(0 To Max, 0 To Max, 0 To Max) As Double 'Changes to B vector field
Dim DeltaBY(0 To Max, 0 To Max, 0 To Max) As Double
Dim DeltaBZ(0 To Max, 0 To Max, 0 To Max) As Double

'*** Initialize the Simulation ***
'
Private Sub InitializeSimulation(Structure As Integer)
  Dim X As Integer    'Index along X-axis
  Dim Y As Integer    'Index along Y-axis
  Dim Z As Integer    'Index along Z-axis
  Dim EPX As Double   'Electric field at point P (X component)
  Dim EPY As Double   'Electric field at point P (Y component)
  Dim EPZ As Double   'Electric field at point P (Z component)
  Dim BPX As Double   'Magnetic field at point P (X component)
  Dim BPY As Double   'Magnetic field at point P (Y component)
  Dim BPZ As Double   'Magnetic field at point P (Z component)
  Dim D As Double     'Cylindrical Distance from center of simulation
  Dim M As Double     'Magnitude variable
  Dim R As Double     'Radius variable
  Dim A As Double
  Dim XT As Double
  '--- Reset time and field arrays ---
  TimeCounter = 0                 'Reset time to zero
  TimeClock = 0
  For X = 0 To LimitX              'Loop through every cell
    For Y = 0 To LimitY
      For Z = 0 To LimitZ
        EX(X, Y, Z) = 0: EY(X, Y, Z) = 0: EZ(X, Y, Z) = 0  'Reset E field
        BX(X, Y, Z) = 0: BY(X, Y, Z) = 0: BZ(X, Y, Z) = 0  'Reset B field
      Next
    Next
  Next
  If Structure = 0 Then Exit Sub  'For empty space, take no further action
  
  
  '--- Create a plane wave (3 wavelengths across the simulation) ---
  If Structure = 1 Then
    For X = 0 To LimitX
      EPY = 10000 * Sin(6 * PI * CDbl(X) / (LimitX + 1))
      For Y = 0 To LimitY
        For Z = 0 To LimitZ
          EY(X, Y, Z) = EPY   'Set electric field
          BZ(X, Y, Z) = EPY   'Set the magnetic field
        Next
      Next
    Next
    Exit Sub
  End If

  
  '--- Wave Pulses ---
  If Structure = 2 Or Structure = 3 Then
    A = -PI / 4           'Set the angle of propagation
    For X = 0 To LimitX
      For Y = 0 To LimitY
        M = 3 * (Sin(A) * X + Cos(A) * Y) / (LimitX + 1)
        If Structure = 2 Then   'For a simple wave pulse, truncate the wave
          R = 0
          If Abs(X + Y - 2 * Center) < Center / 3 Then
            If Abs(X - Y + Center / 2) < Center / 2 Then R = 1
          End If
        Else                    'Otherwise, apply Gaussian distribution
          R = Exp(-30 * ((X - 0.7 * Center) ^ 2 + (Y - 1.3 * Center) ^ 2) / (LimitX + 1) / (LimitY + 1))
        End If
        EPY = R * 10000 * Sin(A) * Sin(4 * PI * M)
        EPX = R * -10000 * Cos(A) * Sin(4 * PI * M)
        For Z = 0 To LimitZ
          EX(X, Y, Z) = EPX
          EY(X, Y, Z) = EPY
          BZ(X, Y, Z) = -R * 10000 * Sin(4 * PI * M)
        Next
      Next
    Next
    Exit Sub
  End If
  
  '--- Create single-point disturbance ---
  If Structure = 4 Then
    EY(Center, Center, LimitZ \ 2) = 10000
    Exit Sub
  End If
  
  '--- Superposition Demo ---
  If Structure = 5 Then
    R = LimitX \ 8
    For X = 0 To LimitX
      For Y = 0 To LimitY
        For Z = 0 To LimitZ
          If X >= Center Then
            XT = X - Center - R
          Else
            XT = X - Center + R
          End If
            
          
          D = Sqr((XT + 0.5) ^ 2 + (Y - Center + 0.5) ^ 2)
          M = 7000 * Exp(-100 * D ^ 2 / (LimitX ^ 2))
          
          EPX = M * -(Y - Center + 0.5) / D
          EPY = M * (XT + 0.5) / D
          EPZ = 0
          BPX = 0
          BPY = 0
          BPZ = M
          EX(X, Y, Z) = EPX
          EY(X, Y, Z) = EPY
          EZ(X, Y, Z) = EPZ
          BX(X, Y, Z) = BPX
          BY(X, Y, Z) = BPY
          BZ(X, Y, Z) = BPZ
         
        Next
      Next
    Next
    Exit Sub
  End If
    

          
  
  
  '--- Create Cylindrical Disturbance ---
  R = CDbl(LimitX + 1) / 2 - 2
  For X = -R - 1 To R
    For Y = -R - 1 To R
      For Z = 0 To LimitZ
        D = Sqr((X + 0.5) ^ 2 + (Y + 0.5) ^ 2)
        If Structure = 6 Then
          M = 10000 * Exp(-0.3 * D ^ 2)
          EPX = 0
          EPY = M
          EPZ = 0
          BPX = 0
          BPY = 0
          BPZ = 0
        Else
          If D < Center / 2 Then
            M = 10000 / Sqr(D) * Sin(PI * D / (LimitX + 1) * 8)
          Else
            M = 0
          End If
          EPX = M * -(Y + 0.5) / D
          EPY = M * (X + 0.5) / D
          EPZ = 0
          BPX = 0
          BPY = 0
          BPZ = M
        End If
        EX(Center + X, Center + Y, Z) = EPX
        EY(Center + X, Center + Y, Z) = EPY
        EZ(Center + X, Center + Y, Z) = EPZ
        BX(Center + X, Center + Y, Z) = BPX
        BY(Center + X, Center + Y, Z) = BPY
        BZ(Center + X, Center + Y, Z) = BPZ
      Next
    Next
  Next
  
  
End Sub

'*** Plot the Lattice on the VB Form ***
'
Private Sub ShowLattice()
  Dim V As Double
  Dim X As Integer
  Dim Y As Integer
  Dim SX As Single
  Dim SY As Single
  Dim PX As Single
  Dim PY As Single
  Dim C As Long
  Dim GridLines As Boolean
  Dim ShowBField As Boolean
  SX = Pic1.ScaleWidth / (LimitX + 1)   'Compute X & Y scale factor
  SY = Pic1.ScaleWidth / (LimitY + 1)
  GridLines = (Check1.Value = 1)        'Get gridline setting
  ShowBField = (Check2.Value = 1)       'Get B-field setting
  
  For X = 0 To LimitX
    For Y = 0 To LimitY
      PX = X * SX             'Determine X,Y screen location
      PY = (LimitY - Y) * SY
      
      '--- Draw E-field cells ---
      V = 240 * EY(X, Y, LimitZ \ 2) / 8000
      If V < 0 Then
        If V < -255 Then V = -255
        C = RGB(-V, 0, 0)           'Show negative E in red
      Else
        If V > 255 Then V = 255
        C = RGB(0, 0, V)            'Show positive E in blue
      End If
      Pic1.Line (PX, PY)-(PX + SX, PY + SY), C, BF
      If GridLines Then
        Pic1.Line (PX, PY)-(PX + SX, PY + SY), RGB(64, 64, 64), B
      End If
      If ShowBField Then
        V = BZ(X, Y, LimitZ \ 2)
        If V > 2000 Then          'Show positive B-field
          V = 256 * V / 10000
          If V > 255 Then V = 255
          Pic1.Line (PX + SX \ 2, PY + SY \ 2)-(PX + SX \ 2 + 15, PY + SY \ 2 + 15), RGB(0, V, 0), BF
        ElseIf V < -2000 Then     'Show negative B-field
          V = 256 * -V / 10000
          Pic1.Line (PX + SX \ 2 - 30, PY + SY \ 2 - 30)-(PX + SX \ 2 + 45, PY + SY \ 2 + 45), RGB(0, V, 0)
          Pic1.Line (PX + SX \ 2 + 30, PY + SY \ 2 - 30)-(PX + SX \ 2 - 45, PY + SY \ 2 + 45), RGB(0, V, 0)
        End If
      End If
      
      '--- Draw S cells ---
      V = EX(X, Y, LimitZ \ 2) * EX(X, Y, LimitZ \ 2) _
        + EY(X, Y, LimitZ \ 2) * EY(X, Y, LimitZ \ 2) _
        + EZ(X, Y, LimitZ \ 2) * EZ(X, Y, LimitZ \ 2)
      V = 255 * V / (50000000)
      If V < 0 Then V = 0
      If V > 255 Then V = 255
      C = RGB(V, V, V)
      Pic2.Line (PX, PY)-(PX + SX, PY + SY), C, BF
      If GridLines Then
        Pic2.Line (PX, PY)-(PX + SX, PY + SY), RGB(64, 64, 64), B
      End If
    
    Next
  Next
End Sub



'*** Handle transition from one state to the next ***

Private Sub Transition()
  Dim X As Integer    'X,Y,Z subscript counter
  Dim Y As Integer
  Dim Z As Integer
  Dim X1 As Integer   'X,Y,Z neighbor subscripts
  Dim X2 As Integer
  Dim Y1 As Integer
  Dim Y2 As Integer
  Dim Z1 As Integer
  Dim Z2 As Integer
  Dim D As Double     'Change in numerical value
  Dim TC As Double    'Time-based constant (related to delta-t)
  On Error Resume Next  'Skip errors when floating-point limit is reached
    
  '--- Perform time-based calculations for each iteration ---
  D = DeltaTime.ListIndex + 1
  TC = (1 + LimitX) / (2 * D * 100)
  TimeCounter = TimeCounter + 1
  TimeClock = TimeClock + 60 / D
  TimeLabel = "Iteration: " & TimeCounter & " (" & TimeClock & " as)"
  
  '--- Determine changes needed based on neighboring cells ---
  For X = 0 To LimitX
    X1 = X - 1: If X1 < 0 Then X1 = LimitX
    X2 = X + 1: If X2 > LimitX Then X2 = 0
    For Y = 0 To LimitY
      Y1 = Y - 1: If Y1 < 0 Then Y1 = LimitY
      Y2 = Y + 1: If Y2 > LimitY Then Y2 = 0
         
      For Z = 0 To LimitZ
        Z1 = Z - 1: If Z1 < 0 Then Z1 = LimitZ
        Z2 = Z + 1: If Z2 > LimitZ Then Z2 = 0
        
        'dBZ/dy - dBY/dz = dEX / dt
        'dBX/dz - dBZ/dx = dEY / dt
        'dBY/dx - dBX/dy = dEZ / dt
        D = BZ(X, Y2, Z) - BZ(X, Y1, Z) - (BY(X, Y, Z2) - BY(X, Y, Z1))
        DeltaEX(X, Y, Z) = D * TC
        
        D = BX(X, Y, Z2) - BX(X, Y, Z1) - (BZ(X2, Y, Z) - BZ(X1, Y, Z))
        DeltaEY(X, Y, Z) = D * TC
        
        D = BY(X2, Y, Z) - BY(X1, Y, Z) - (BX(X, Y2, Z) - BX(X, Y1, Z))
        DeltaEZ(X, Y, Z) = D * TC
        
        'dEZ/dy - dEY/dz = - dBX / dt
        'dEX/dz - dEZ/dx = - dBY / dt
        'dEY/dx - dEX/dy = - dBZ / dt
        D = EY(X, Y, Z2) - EY(X, Y, Z1) - (EZ(X, Y2, Z) - EZ(X, Y1, Z))
        DeltaBX(X, Y, Z) = D * TC
        D = EZ(X2, Y, Z) - EZ(X1, Y, Z) - (EX(X, Y, Z2) - EX(X, Y, Z1))
        DeltaBY(X, Y, Z) = D * TC
        D = EX(X, Y2, Z) - EX(X, Y1, Z) - (EY(X2, Y, Z) - EY(X1, Y, Z))
        DeltaBZ(X, Y, Z) = D * TC
      Next
    Next
  Next

  '--- Apply changes to E & B vectors ---
  For X = 0 To LimitX
    For Y = 0 To LimitY
      For Z = 0 To LimitZ
        EX(X, Y, Z) = (EX(X, Y, Z) + DeltaEX(X, Y, Z))
        EY(X, Y, Z) = (EY(X, Y, Z) + DeltaEY(X, Y, Z))
        EZ(X, Y, Z) = (EZ(X, Y, Z) + DeltaEZ(X, Y, Z))
        BX(X, Y, Z) = (BX(X, Y, Z) + DeltaBX(X, Y, Z))
        BY(X, Y, Z) = (BY(X, Y, Z) + DeltaBY(X, Y, Z))
        BZ(X, Y, Z) = (BZ(X, Y, Z) + DeltaBZ(X, Y, Z))
      Next
    Next
  Next


End Sub

'*** Event Procedure to Reset Simulation ***
'
Private Sub Command1_Click()
  WaveType_Click
End Sub

'*** List item selection ***
'
Private Sub WaveType_Click()
  Dim I As Integer
  I = WaveType.ListIndex    'Determine simulation structure
  If I < 0 Then Exit Sub    'If none selected, don't continue
  InitializeSimulation I    'Initialize the simulation
  
End Sub

'*** Change to spatial granularity ***
'
Private Sub DeltaSpace_Click()
  Dim I As Integer
  I = 10 * (DeltaSpace.ListIndex + 1) 'Determine matrix size (10,20,30,40,50)
  LimitX = I - 1                      'Set limits on X,Y,Z subscripts
  LimitY = I - 1
  LimitZ = 0
  Center = I \ 2                      'Calculate simulation center
  WaveType_Click                      'Reset the simulation
End Sub

'*** Start-up Procedure ***
'
Private Sub Form_Load()
  
  '--- Fill the drop-down lists with appropriate values ---
  With DeltaSpace
    .AddItem "180 nm (10 X 10)"
    .AddItem "90 nm  (20 X 20)"
    .AddItem "60 nm  (30 X 30)"
    .AddItem "45 nm  (40 X 40)"
    .AddItem "36 nm  (50 X 50)"
  End With
  With DeltaTime
    .AddItem "60 as"
    .AddItem "30 as"
    .AddItem "20 as"
    .AddItem "15 as"
    .AddItem "12 as"
    .AddItem "10 as"
  End With
  With WaveType
    .AddItem "Empty Space"
    .AddItem "Plane Wave (600 nm)"
    .AddItem "Simple Wave Pulse (200 nm)"
    .AddItem "Gaussian Wave Pulse (200 nm)"
    .AddItem "Point Disturbance"
    .AddItem "Superposition"
    .AddItem "Gaussian Disturbance"
    .AddItem "Cylindrical Wave (450 nm)"
  End With
  
  '--- Set the initial drop-down list values ---
  DeltaTime.ListIndex = 5
  DeltaSpace.ListIndex = 1
  WaveType.ListIndex = 0
End Sub


'*** Run Timed Events ***
'
Private Sub Timer1_Timer()
  If Check3.Value = 0 Then    'If simulation is not paused...
    Transition                'perform transition to next point in time
  End If
  ShowLattice                 'Refresh the display
End Sub