'PHOTOELECTRIC EFFECT - CLASSICAL MODEL
'by A.S.Kapoyannis and G.Kalkanis
'University of Athens
'2015
Public Class Classical
    Dim bmpBackBuffer As Bitmap
    Dim directGraphics As Graphics
    Dim w As Integer, h As Integer
    Dim Frequency, Intensity As Double
    Dim randx, randy, xx(400), yy(400), ux(400), uy(400), umax, uxmax, duxmax, uxmaxi(400), umean, tlag1, dtlag1, tlag30, Ic, dIc, Ici(1000) As Double
    Dim umax1i(1000), umax2i(1000), umax3i(1000), umax1ii(1000), umax2ii(1000), umax3ii(1000), dumax1ii(1000), dumax2ii(1000), dumax3ii(1000) As Double
    Dim Count(1000) As Integer
    Dim uoscil(400), toscil(400), phioscil(400) As Double
    Dim Noscil(400) As Integer
    Dim UmaxMethod(3) As Integer
    Dim tlag11i(1000), tlag12i(1000), tlag11ii(1000), tlag12ii(1000), dtlag11ii(1000), dtlag12ii(1000), Nescapetime(1000) As Double
    Dim TlagMethod(3) As Integer
    Dim xxesc(400), yyesc(400), uxesc(400), uyesc(400) As Double
    Dim Nescape(400), Nescapevelocity(400), nescprevious, Ncycle, icycle, nuesccycle(600), nuesccycleed(600) As Integer
    Dim sumu As Double
    Dim mulenergy As Double = 4.138576779 / 100.0
    Dim mulvelocity As Double = 1.205705686
    Dim Frequencymin, Frequencymax, Frequencystep1, Frequencystep2, Intensitymin, Intensitymax, Intensitystep1, Intensitystep2 As Double
    Dim Frequencyi(300), Intensityi(300), Icii(300), dIcii(300) As Double
    Dim yesc(300) As Integer
    Dim Ne, Nesc, Ninside, nescmax, ReplaceEscaped, ShowEnergyEscaped As Integer
    Dim istop, ydraw, startlight, lightreach, imode, NmodeFrequency, NmodeIntensity, Nmode As Integer
    Private Sub Form1_Load(ByVal sender As Object, ByVal e As System.EventArgs) Handles Me.Load
        w = PicMain.ClientRectangle.Width
        h = PicMain.ClientRectangle.Height
        bmpBackBuffer = New Bitmap(w, h)
        directGraphics = PicMain.CreateGraphics
        TextBox14.Text = HScrollBar8.Value
    End Sub
    Private Sub CheckBox4_CheckedChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles CheckBox4.CheckedChanged
        'Constant Frequency
        If CheckBox4.Checked = False And CheckBox5.Checked = False Then CheckBox5.Checked = True
        If CheckBox4.Checked = True Then
            HScrollBar4.Value = HScrollBar1.Value
            TextBox10.Text = HScrollBar4.Value * 10.0
            TextBox1.Text = HScrollBar1.Value * 10.0
            HScrollBar6.Value = 0
            TextBox12.Text = 0
        End If
    End Sub
    Private Sub HScrollBar1_Scroll(ByVal sender As System.Object, ByVal e As System.Windows.Forms.ScrollEventArgs) Handles HScrollBar1.Scroll
        'Minimum Frequency
        If CheckBox4.Checked = True Or HScrollBar4.Value <= HScrollBar1.Value Then
            HScrollBar4.Value = HScrollBar1.Value
            TextBox10.Text = HScrollBar4.Value * 10.0
            HScrollBar6.Value = 0
            TextBox12.Text = 0
        End If
        TextBox1.Text = HScrollBar1.Value * 10.0
    End Sub
    Private Sub HScrollBar4_Scroll(ByVal sender As System.Object, ByVal e As System.Windows.Forms.ScrollEventArgs) Handles HScrollBar4.Scroll
        'Maximum Frequency
        If CheckBox4.Checked = True Or HScrollBar1.Value >= HScrollBar4.Value Then
            HScrollBar1.Value = HScrollBar4.Value
            TextBox1.Text = HScrollBar1.Value * 10.0
            HScrollBar6.Value = 0
            TextBox12.Text = 0
        End If
        TextBox10.Text = HScrollBar4.Value * 10.0
    End Sub
    Private Sub HScrollBar6_Scroll(ByVal sender As System.Object, ByVal e As System.Windows.Forms.ScrollEventArgs) Handles HScrollBar6.Scroll
        'Frequency points between minimum and maximum
        If CheckBox4.Checked = True Or HScrollBar1.Value = HScrollBar4.Value Then
            HScrollBar6.Value = 0
            TextBox12.Text = 0
        End If
        TextBox12.Text = HScrollBar6.Value
    End Sub
    Private Sub CheckBox5_CheckedChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles CheckBox5.CheckedChanged
        'Constant Irradiance
        If CheckBox5.Checked = False And CheckBox4.Checked = False Then CheckBox4.Checked = True
        If CheckBox5.Checked = True Then
            HScrollBar5.Value = HScrollBar2.Value
            TextBox11.Text = HScrollBar5.Value
            TextBox2.Text = HScrollBar2.Value
            HScrollBar7.Value = 0
            TextBox13.Text = 0
        End If
    End Sub
    Private Sub HScrollBar2_Scroll(ByVal sender As System.Object, ByVal e As System.Windows.Forms.ScrollEventArgs) Handles HScrollBar2.Scroll
        'Minimum Irradiance
        If CheckBox5.Checked = True Or HScrollBar5.Value <= HScrollBar2.Value Then
            HScrollBar5.Value = HScrollBar2.Value
            TextBox11.Text = HScrollBar5.Value
        End If
        TextBox2.Text = HScrollBar2.Value
    End Sub
    Private Sub HScrollBar5_Scroll(ByVal sender As System.Object, ByVal e As System.Windows.Forms.ScrollEventArgs) Handles HScrollBar5.Scroll
        'Maximum Irradiance
        If CheckBox5.Checked = True Or HScrollBar2.Value >= HScrollBar5.Value Then
            HScrollBar2.Value = HScrollBar5.Value
            TextBox2.Text = HScrollBar2.Value
        End If
        TextBox11.Text = HScrollBar5.Value
    End Sub
    Private Sub HScrollBar7_Scroll(ByVal sender As System.Object, ByVal e As System.Windows.Forms.ScrollEventArgs) Handles HScrollBar7.Scroll
        'Irradiance points between minimum and maximum
        If CheckBox5.Checked = True Or HScrollBar2.Value = HScrollBar5.Value Then
            HScrollBar7.Value = 0
            TextBox13.Text = 0
        End If
        TextBox13.Text = HScrollBar7.Value
    End Sub
    Private Sub HScrollBar8_Scroll(ByVal sender As System.Object, ByVal e As System.Windows.Forms.ScrollEventArgs) Handles HScrollBar8.Scroll
        '# of Repetitions
        TextBox14.Text = HScrollBar8.Value
    End Sub
    Private Sub CheckBox6_CheckedChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles CheckBox6.CheckedChanged
        'Diagram esc.e-f
        If CheckBox6.Checked = True And CheckBox11.Checked = True Then CheckBox11.Checked = False
    End Sub
    Private Sub CheckBox11_CheckedChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles CheckBox11.CheckedChanged
        'Diagram Ic-Ee
        If CheckBox11.Checked = True And CheckBox6.Checked = True Then CheckBox6.Checked = False
    End Sub
    Private Sub CheckBox7_CheckedChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles CheckBox7.CheckedChanged
        'Diagram Vmax-f
        If CheckBox7.Checked = True And CheckBox9.Checked = True Then CheckBox9.Checked = False
    End Sub
    Private Sub CheckBox9_CheckedChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles CheckBox9.CheckedChanged
        'Diagram Vmax-Ee
        If CheckBox9.Checked = True And CheckBox7.Checked = True Then CheckBox7.Checked = False
    End Sub
    Private Sub CheckBox8_CheckedChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles CheckBox8.CheckedChanged
        'Diagram tlag-f
        If CheckBox8.Checked = True And CheckBox10.Checked = True Then CheckBox10.Checked = False
    End Sub
    Private Sub CheckBox10_CheckedChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles CheckBox10.CheckedChanged
        'Diagram tlag-Ee
        If CheckBox10.Checked = True And CheckBox8.Checked = True Then CheckBox8.Checked = False
    End Sub
    Private Sub Button2_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button2.Click
        Dim i As Integer
        'ShowEnergyEscaped = 0: The Energy Distribution of Escaped Electrons is NOT SHOWN
        'ShowEnergyEscaped = 1: The Energy Distribution of Escaped Electrons is SHOWN
        ShowEnergyEscaped = 0
        'ReplaceEscaped=0: The escaped e- ARE NOT replaced
        'ReplaceEscaped=1: The escaped e- ARE replaced
        ReplaceEscaped = 1
        'UmaxMethod(1)=0: Umax is NOT calculated through maximum uesc
        'UmaxMethod(1)=1: Umax is calculated through maximum uesc
        'UmaxMethod(2)=0: Umax is NOT calculated through linear fit to the distribution of uesc
        'UmaxMethod(2)=1: Umax is calculated through linear fit to the distribution of uesc
        'UmaxMethod(3)=0: Umax is NOT calculated through mean and deviation to the distribution of uesc
        'UmaxMethod(3)=1: Umax is calculated through mean and deviation to the distribution of uesc
        UmaxMethod(1) = 0
        UmaxMethod(2) = 0
        UmaxMethod(3) = 1
        'TlagMethod(1)=0: Tlag is NOT calculated through tlag1 (the time the 1st e- escapes)
        'TlagMethod(1)=1: Tlag is calculated through tlag1 (the time the 1st e- escapes)
        'TlagMethod(2)=0: Tlag is NOT calculated through linear fit to the Nesc-t diagram
        'TlagMethod(2)=1: Tlag is calculated through linear fit to the Nesc-t diagram
        TlagMethod(1) = 0
        TlagMethod(2) = 1
        Frequencymin = HScrollBar1.Value * 10.0
        Frequencymax = HScrollBar4.Value * 10.0
        If CheckBox4.Checked = True Or HScrollBar1.Value = HScrollBar4.Value Then
            NmodeFrequency = 1
        Else
            NmodeFrequency = HScrollBar6.Value + 2
        End If
        Intensitymin = HScrollBar2.Value
        Intensitymax = HScrollBar5.Value
        If CheckBox5.Checked = True Or HScrollBar2.Value = HScrollBar5.Value Then
            NmodeIntensity = 1
        Else
            NmodeIntensity = HScrollBar7.Value + 2
        End If
        HScrollBar1.Enabled = False
        HScrollBar4.Enabled = False
        HScrollBar6.Enabled = False
        HScrollBar2.Enabled = False
        HScrollBar5.Enabled = False
        HScrollBar7.Enabled = False
        HScrollBar8.Enabled = False
        Button4.Enabled = True
        Frequencyi = XPoints(Frequencymin, Frequencymax, NmodeFrequency, 1)
        Intensityi = XPoints(Intensitymin, Intensitymax, NmodeIntensity, 2)
        Ncycle = HScrollBar8.Value
        imode = 1
        icycle = 1
        For i = 1 To Ncycle
            Count(i) = 0
        Next i
        Ne = 100
        InitialConditions()
        Nesc = 0
        For i = 1 To NmodeFrequency
            yesc(i) = 0
        Next i
        nescmax = 50
        Timer1.Enabled = True
    End Sub
    Private Function XPoints(ByVal xmin As Double, ByVal xmax As Double, ByVal N As Integer, ByVal idecide As Integer) As Double()
        'Returns N positive points ordered between xmin and xmax
        Dim i, N1, N2, icurrent, i12 As Integer
        Dim A1(100), A2(100), step1, step2, xi(100), B As Double
        If idecide = 1 Then
            If N = 1 Then xi(1) = xmin
            If N > 1 Then
                step1 = (xmax - xmin) / (N - 1)
                For i = 1 To N
                    xi(i) = xmin + step1 * (i - 1)
                Next i
            End If
        End If
        If idecide = 2 Then
            If N > 1 Then
                N1 = 2 + Math.Ceiling((N - 2) / 2)
                N2 = N - N1
                step1 = (xmax - xmin) / (N1 - 1)
                step2 = (1.0 / xmin - 1.0 / xmax) / (N2 + 1)
            Else
                step1 = 0.0
                step2 = 0.0
            End If
            For i = 1 To N1
                A1(i) = xmin + (i - 1) * step1
            Next i
            For i = 1 To N2
                A2(i) = 1.0 / (1.0 / xmax + i * step2)
            Next i
            xi(1) = xmin
            xi(N) = xmax
            For k = 2 To N - 1
                B = xmax
                For i = 2 To N1
                    If A1(i) < B Then
                        B = A1(i)
                        icurrent = i
                        i12 = 1
                    End If
                Next i
                For i = 1 To N2
                    If A2(i) < B Then
                        B = A2(i)
                        icurrent = i
                        i12 = 2
                    End If
                Next i
                xi(k) = B
                If i12 = 1 Then
                    A1(icurrent) = xmax
                Else
                    A2(icurrent) = xmax
                End If
            Next k
        End If
        Return xi
    End Function
    Private Sub InitialConditions()
        Dim rand As New Random()
        Dim Ni As Integer
        'Set initial positions and velocities
        Ninside = Ne
        Ni = Ne
        Do Until Ni < 1
            randx = 450 * rand.NextDouble() + 25
            randy = 450 * rand.NextDouble() + 25
            xx(Ni) = randx
            yy(Ni) = randy
            randx = 4 * rand.NextDouble() - 2
            randy = 4 * rand.NextDouble() - 2
            ux(Ni) = randx
            uy(Ni) = randy
            Nescape(Ni) = 0
            Nescapevelocity(Ni) = 0
            Noscil(Ni) = 0
            uoscil(Ni) = 0
            toscil(Ni) = 0
            phioscil(Ni) = 0
            Ni = Ni - 1
        Loop
        For i = 1 To Ne
            xxesc(i) = 0
            yyesc(i) = 0
            uxesc(i) = 0
            uyesc(i) = 0
        Next i
        umax = 0.0
        uxmax = 0.0
        umean = 0.0
        tlag1 = 0.0
        tlag30 = 0.0
        startlight = 0
        lightreach = 0
        Button2.Enabled = False
        HScrollBar1.Enabled = False
        HScrollBar2.Enabled = False
    End Sub
    Private Sub CheckBox1_CheckedChanged(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles CheckBox1.CheckedChanged
        If istop = 0 Then
            If CheckBox1.Checked = True Then
                Timer1.Enabled = False
                Timer2.Enabled = True
            End If
            If CheckBox1.Checked = False Then
                Timer2.Enabled = False
                Timer1.Enabled = True
            End If
        Else
            Timer1.Enabled = False
            Timer2.Enabled = False
        End If
    End Sub
    Private Sub Timer1_Tick(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Timer1.Tick
        Timer1.Interval = HScrollBar3.Value
        Simulation()
    End Sub
    Private Sub Timer2_Tick(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Timer2.Tick
        Timer2.Interval = 1
        Nesc = 0
        nescprevious = 0
        Do While Nesc < nescmax And istop = 0
            Simulation()
        Loop
    End Sub
    Private Sub Simulation()
        Dim Eo, c, w, w0, tstep, xmin, xmax, ymin, ymax, uxcut As Double
        Dim N, Ni, Nj, Nx, Ny, ii, jj, i_col_el, Ncount, Nxcen, Nycen, grid, xgrid, ygrid, Nxmin, Nymin, Ncol(1000), tesc, Remain As Integer
        Dim Ngrid(6, 6), Impact, energy, velocitycut, bb, B, D, xx1, yy1, upar, utr, tpr, costh, sinth As Double
        Dim xcen, ycen, dxcen, dycen, Dx, Dy, cosf, p1xmin, p1xmax, p1x, p2x, p1y, p2y, EEE, PPP, ms2ep2 As Double
        Dim ax, ay, bx, by, r1, r2, An(20), Rls(5), Rls1(6), MV(2) As Double
        Dim nu(1000), nuesc(1000), nuescvel(200), ndraw, nndraw, nnndraw, narrow, Nxstart, Nystart, nnx, nnnx As Integer
        Dim nustepin, nustepesc, uesccur, uesccurx, extractionwork, numax As Double
        Dim i, ifrequency, iintensity As Integer
        Dim yem, xem, yi, yemnext As Integer
        Dim P As New Pen(Color.Blue, 1)
        Dim Pp As New Pen(Color.Green, 1)
        Dim Pr As New Pen(Color.Green, 1)
        Dim Pe As New Pen(Color.AntiqueWhite, 10)
        Dim Pbb As New Pen(Color.Black, 2)
        Dim Pred As New Pen(Color.Red, 2)
        Dim Por As New Pen(Color.OrangeRed, 2)
        Dim rand As New Random()
        Dim drawFont As New Font("Arial", 9)
        Dim drawFontscript As New Font("Arial", 7)
        Dim G As Graphics = Graphics.FromImage(bmpBackBuffer)
        Static t As Double = 0
        Static tem As Double = 0
        istop = 0
        If CheckBox4.Checked = True Then
            ifrequency = 1
            iintensity = imode
            Nmode = NmodeIntensity
        End If
        If CheckBox5.Checked = True Then
            ifrequency = imode
            iintensity = 1
            Nmode = NmodeFrequency
        End If
        Frequency = Frequencyi(ifrequency) / 10.0
        Intensity = Intensityi(iintensity)
        TextBox15.Text = "Iteration:" & icycle & " / " & Ncycle
        TextBox9.Text = "Point:" & imode & " / " & Nmode
        TextBox16.Text = "f=" & Frequency * 10.0 & " THz"
        TextBox17.Text = "Ee=" & Intensity & " W/m2"
        'Calcucate Nesc
        Nesc = 0
        For i = 1 To Math.Min(nescmax + 5, Ne)
            If xxesc(i) > 0 Then Nesc = Nesc + 1
        Next i
        If Nesc >= nescmax Then Nesc = 0
        c = 3D
        Eo = Math.Sqrt(2 * Intensity)
        w = 2 * Math.PI * Frequency
        w0 = 2 * Math.PI * 1
        xcen = 250.0
        ycen = 250.0
        Nxcen = xcen
        Nycen = ycen
        grid = 90
        dxcen = 10.0
        dycen = 10.0
        Impact = 10.0
        tesc = 0
        Ni = Ne
        xmin = xcen - 225
        xmax = xcen + 225
        ymin = ycen - 225
        ymax = ycen + 225
        Nxmin = xmin
        Nymin = ymin
        uxcut = 2.2
        velocitycut = Math.Sqrt(30.0)
        extractionwork = velocitycut ^ 2 - 10.0
        Remain = 0
        ii = 0
        If startlight <> 1 Then
            tem = 0
            GoTo 30
        End If
        If c * tem >= xcen Then lightreach = 1
        xem = c * tem
        If xem > xcen Then
            xem = xcen
        End If
30:     N = Ne
        tstep = 0.5
        Ni = N
        Ncount = 0
        i_col_el = 0
        'Set the next step positions
        Do Until Ni < 1
            If Nescape(Ni) = 0 Then
                xx(Ni) = xx(Ni) + ux(Ni) * tstep
                yy(Ni) = yy(Ni) + uy(Ni) * tstep
            End If
            Ni = Ni - 1
        Loop
        For i = 1 To Nesc
            If xxesc(i) < 1000 Then
                xxesc(i) = xxesc(i) + uxesc(i) * tstep
                yyesc(i) = yyesc(i) + uyesc(i) * tstep
            Else
                xxesc(i) = 1000
                yyesc(i) = i * 5
            End If
        Next i
        'Set collision matrix all 0
        Ni = N
        Do Until Ni < 1
            Ncol(Ni) = 0
            Ni = Ni - 1
        Loop
        Ni = N
        Do Until Ni < 1
            '
            'INTERACTION WITH E/M FIELD
            '
            If startlight <> 1 Or lightreach <> 1 Then GoTo 50
            If (xx(Ni) - xcen) ^ 2 + (yy(Ni) - ycen) ^ 2 < 40000 And Nescape(Ni) = 0 Then
                randx = 2 * rand.NextDouble()
                randy = 2 * rand.NextDouble() - 1
                If Noscil(Ni) = 1 Then
                    toscil(Ni) = toscil(Ni) + tstep
                    uy(Ni) = uoscil(Ni) + (Eo / Frequency) * (Math.Cos(phioscil(Ni)) - Math.Cos(2 * Math.PI * Frequency * toscil(Ni) / 500 + phioscil(Ni)))
                Else
                    Noscil(Ni) = 1
                    randx = 2 * Math.PI * rand.NextDouble()
                    phioscil(Ni) = randx
                    uoscil(Ni) = uy(Ni)
                    toscil(Ni) = 0
                    uy(Ni) = uoscil(Ni) + (Eo / Frequency) * (Math.Cos(phioscil(Ni)) - Math.Cos(2 * Math.PI * Frequency * toscil(Ni) / 500 + phioscil(Ni)))
                End If
            Else
                Noscil(Ni) = 0
            End If
            'If the electron does not have enough energy...
50:         If (ux(Ni) ^ 2 + uy(Ni) ^ 2 <= velocitycut ^ 2 Or ux(Ni) <= 0 Or Math.Abs(uy(Ni) / ux(Ni)) > 1 / Math.Sqrt(2)) And Nescape(Ni) = 0 Then
                If xx(Ni) < xmin Then
                    xx(Ni) = xmin + xmin - xx(Ni)
                    ux(Ni) = -ux(Ni)
                    Ncol(Ni) = 1
                    Noscil(Ni) = 0
                End If
                If xx(Ni) > xmax Then
                    xx(Ni) = xmax - (xx(Ni) - xmax)
                    ux(Ni) = -ux(Ni)
                    Ncol(Ni) = 1
                    Noscil(Ni) = 0
                End If
                If yy(Ni) < ymin Then
                    yy(Ni) = ymin + ymin - yy(Ni)
                    uy(Ni) = -uy(Ni)
                    Ncol(Ni) = 1
                    Noscil(Ni) = 0
                End If
                If yy(Ni) > ymax Then
                    yy(Ni) = ymax - (yy(Ni) - ymax)
                    uy(Ni) = -uy(Ni)
                    Ncol(Ni) = 1
                    Noscil(Ni) = 0
                End If
                If (ux(Ni) ^ 2 + uy(Ni) ^ 2 > velocitycut ^ 2 And ux(Ni) > 0 And Math.Abs(uy(Ni) / ux(Ni)) <= 1 / Math.Sqrt(2)) Then
                    Nescape(Ni) = 1
                    Noscil(Ni) = 0
                    GoTo 110
                End If
                'Collision with grid ions
                ii = 1
                Do Until ii > 6
                    jj = 1
                    Do Until jj > 6
                        xgrid = Nxmin + (ii - 1) * grid
                        Dx = -(xgrid - xx(Ni))
                        ygrid = Nymin + (jj - 1) * grid
                        Dy = -(ygrid - yy(Ni))
                        cosf = (Dx * ux(Ni) + Dy * uy(Ni)) / (Math.Sqrt((Dx ^ 2 + Dy ^ 2) * (ux(Ni) ^ 2 + uy(Ni) ^ 2)))
                        If Math.Sqrt(Dx ^ 2 + Dy ^ 2) <= Impact Then
                            'Calculation of the point of impact with the grid ion
                            If ux(Ni) <> 0 Then
                                B = ((yy(Ni) * ux(Ni) - xx(Ni) * uy(Ni)) / ux(Ni)) - ygrid
                                bb = (B * uy(Ni) / ux(Ni)) - xgrid
                                D = bb ^ 2 - (1.0 + (uy(Ni) / ux(Ni)) ^ 2) * (xgrid ^ 2 + B ^ 2 - Impact ^ 2)
                                If ux(Ni) > 0 Then
                                    xx1 = (-bb - Math.Sqrt(D)) / (1.0 + (uy(Ni) / ux(Ni)) ^ 2)
                                End If
                                If ux(Ni) < 0 Then
                                    xx1 = (-bb + Math.Sqrt(D)) / (1.0 + (uy(Ni) / ux(Ni)) ^ 2)
                                End If
                                yy1 = yy(Ni) + (xx1 * uy(Ni) - xx(Ni) * uy(Ni)) / ux(Ni)
                                tpr = (xx(Ni) - xx1) / ux(Ni)
                            Else
                                xx1 = xx(Ni)
                                If uy(Ni) > 0 Then
                                    yy1 = ygrid - Math.Sqrt(Impact ^ 2 - (xx1 - xgrid) ^ 2)
                                End If
                                If uy(Ni) < 0 Then
                                    yy1 = ygrid + Math.Sqrt(Impact ^ 2 - (xx1 - xgrid) ^ 2)
                                End If
                                tpr = (yy(Ni) - yy1) / uy(Ni)
                            End If
                            'Calculation of new velocity
                            ax = (xx1 - xgrid)
                            ay = (yy1 - ygrid)
                            bx = ux(Ni)
                            by = uy(Ni)
                            ux(Ni) = ((-ax ^ 2 + ay ^ 2) * bx - 2 * ax * ay * by) / (ax ^ 2 + ay ^ 2)
                            uy(Ni) = ((ax ^ 2 - ay ^ 2) * by - 2 * ax * ay * bx) / (ax ^ 2 + ay ^ 2)
                            'Calculation of new position after the collision
                            xx(Ni) = xx1 + ux(Ni) * tpr
                            yy(Ni) = yy1 + uy(Ni) * tpr
                            'If it collides with ion no other collision is permitted in this round
                            Ncol(Ni) = 1
                            Noscil(Ni) = 0
                        End If
                        jj = jj + 1
                    Loop
                    ii = ii + 1
                Loop
44:             If (ux(Ni) ^ 2 + uy(Ni) ^ 2 > velocitycut ^ 2 And ux(Ni) > 0 And Math.Abs(uy(Ni) / ux(Ni)) <= 1 / Math.Sqrt(2)) Then
                    Nescape(Ni) = 1
                    Noscil(Ni) = 0
                    GoTo 110
                End If
                'Collision with electrons
                Ncount = Ncount + 1
                ii = 0
                Nj = Ni - 1
                If Ncol(Ni) = 1 Then
                    Noscil(Ni) = 0
                    GoTo 100
                End If
                Do Until Nj < 1
                    If Ncol(Nj) = 1 Or (Nescape(Nj) = 1) Then
                        GoTo 90
                    End If
                    r1 = Math.Sqrt((xx(Ni) - xx(Nj)) ^ 2 + (yy(Ni) - yy(Nj)) ^ 2)
                    If r1 <= Impact / 1.0 Then
                        p1xmin = 0.5 * (ux(Ni) + ux(Nj) - Math.Sqrt((ux(Ni) - ux(Nj)) ^ 2 + (uy(Ni) - uy(Nj)) ^ 2))
                        p1xmax = 0.5 * (ux(Ni) + ux(Nj) + Math.Sqrt((ux(Ni) - ux(Nj)) ^ 2 + (uy(Ni) - uy(Nj)) ^ 2))
                        randx = p1xmin + (p1xmax - p1xmin) * rand.NextDouble()
                        p1x = ux(Ni)
                        p2x = ux(Nj)
                        p1y = uy(Ni)
                        p2y = uy(Nj)
                        ux(Ni) = randx
                        ux(Nj) = p1x + p2x - ux(Ni)
                        EEE = p1x ^ 2 + p1y ^ 2 + p2x ^ 2 + p2y ^ 2 - ux(Ni) ^ 2 - ux(Nj) ^ 2
                        PPP = p1y + p2y
                        randy = rand.NextDouble()
                        ms2ep2 = Math.Sqrt(Math.Abs(2 * EEE - PPP ^ 2))
                        uy(Ni) = 0.5 * (PPP + ms2ep2)
                        uy(Nj) = 0.5 * (PPP - ms2ep2)
                        i_col_el = i_col_el + 1
                        Ncol(Ni) = 1
                        Noscil(Ni) = 0
                        Ncol(Nj) = 1
                        Noscil(Nj) = 0
                        If (ux(Ni) ^ 2 + uy(Ni) ^ 2 > velocitycut ^ 2 And ux(Ni) > 0 And Math.Abs(uy(Ni) / ux(Ni)) <= 1 / Math.Sqrt(2)) Then
                            Nescape(Ni) = 1
                            Noscil(Ni) = 0
                            GoTo 100
                        End If
                        If (ux(Nj) ^ 2 + uy(Nj) ^ 2 > velocitycut ^ 2 And ux(Nj) > 0 And Math.Abs(uy(Nj) / ux(Nj)) <= 1 / Math.Sqrt(2)) Then
                            Nescape(Nj) = 1
                            Noscil(Nj) = 0
                        End If
                    End If
90:                 Nj = Nj - 1
                Loop
            Else
                '
                'THE ELECTRON ESCAPES
                '
                Nescape(Ni) = 1
                Noscil(Ni) = 0
            End If
100:        Ni = Ni - 1
        Loop
110:    Ni = Ne
        Dim sum As Double
        Dim iii As Integer
        Do Until Ni < 1
            If Nescape(Ni) = 1 And Nescapevelocity(Ni) = 0 Then
                ax = ux(Ni)
                ay = uy(Ni)
                r1 = Math.Sqrt(1.0 - velocitycut ^ 2 / (ux(Ni) ^ 2 + uy(Ni) ^ 2))
                'r1 = 1.0
                ux(Ni) = ax * r1
                uy(Ni) = ay * r1
                Nescapevelocity(Ni) = 1
                Nesc = Nesc + 1
                If ReplaceEscaped = 0 Then Ninside = Ninside - 1
                xxesc(Nesc) = xx(Ni)
                yyesc(Nesc) = yy(Ni)
                uxesc(Nesc) = ux(Ni)
                uyesc(Nesc) = uy(Ni)
                If ReplaceEscaped = 1 Then
                    Nescape(Ni) = 0
                    Nescapevelocity(Ni) = 0
                    randx = 450 * rand.NextDouble() + 25
                    randy = 450 * rand.NextDouble() + 25
                    xx(Ni) = randx
                    yy(Ni) = randy
                    randx = 4 * rand.NextDouble() - 2
                    randy = 4 * rand.NextDouble() - 2
                    ux(Ni) = randx
                    uy(Ni) = randy
                End If
            End If
            Ni = Ni - 1
        Loop
        If Nesc >= 1 And nescprevious < 1 Then
            tlag1 = t
        End If
        If Nesc >= 30 And nescprevious < 30 Then
            tlag30 = t
        End If
        'Calculate if electrons have escaped for the current frequency
        If Nesc >= 1 Then yesc(ifrequency) = 1
        '
        'CLEAR PREVIOUS GRAPHICS & DRAW A RESTANGLE
        '
        If CheckBox1.Checked = False Or Nesc >= nescmax Then
            G.Clear(Color.White)
            G.FillRectangle(Brushes.AntiqueWhite, Nxcen - 225 - 5, Nycen - 225 - 5, 460, 460)
        End If
        '
        'DRAW THE GRID IONS
        '
        If CheckBox1.Checked = False Or Nesc >= nescmax Then
            For ii = 1 To 6
                For jj = 1 To 6
                    randx = 2 * Math.PI * rand.NextDouble()
                    randy = 3.0 * rand.NextDouble()
                    xgrid = Nxmin + (ii - 1) * grid + randy * Math.Sin(randx)
                    ygrid = Nymin + (jj - 1) * grid + randy * Math.Cos(randx)
                    G.FillEllipse(Brushes.Gray, xgrid - 3, ygrid - 3, 6, 6)
                Next jj
            Next ii
        End If
        '
        'DRAW THE E/M FIELD
        '
        If CheckBox1.Checked = False Or Nesc >= nescmax Then
            Do Until xem < 0
                yem = Eo * Math.Sin((2 * Math.PI * Frequency * (tem - xem / c)) / 300)
                yemnext = Eo * Math.Sin((2 * Math.PI * Frequency * (tem - (xem - 1) / c)) / 300)
                G.DrawLine(P, xem, 250 - yem, xem - 1, 250 - yemnext)
                xem = xem - 1
            Loop
        End If
        tem = tem + 1
        '
        'DRAW THE ELECTRONS
        '
        If CheckBox1.Checked = False Or Nesc >= nescmax Then
            Ni = Ne
            Do Until Ni < 1
                Nx = xx(Ni)
                Ny = yy(Ni)
                If Nescape(Ni) = 0 Then
                    G.FillEllipse(Brushes.Red, Nx - 2, Ny - 2, 4, 4)
                End If
                Ni = Ni - 1
            Loop
            For i = 1 To Nesc
                Nx = xxesc(i)
                Ny = yyesc(i)
                If Nx <= 700 Then
                    G.FillEllipse(Brushes.Green, Nx - 3, Ny - 3, 6, 6)
                End If
            Next i
        End If
        '
        'DISTRIBUTION OF ENERGY
        '
        If CheckBox1.Checked = False Or Nesc >= nescmax Then
            If CheckBox3.Checked = True Then
                nustepin = 1.5
                nustepesc = 1.5
                ii = 1
                Do Until ii = 200
                    nu(ii) = 0
                    nuesc(ii) = 0
                    nuescvel(ii) = 0
                    ii = ii + 1
                Loop
                Ni = Ne
                numax = 0
                Do Until Ni < 1
                    If Nescape(Ni) = 0 Then
                        ii = (ux(Ni) ^ 2 + uy(Ni) ^ 2) / nustepin
                        nu(ii + 1) = nu(ii + 1) + 1
                        If nu(ii + 1) > numax Then numax = nu(ii + 1)
                    End If
                    Ni = Ni - 1
                Loop
                umax = 0
                uxmax = 0
                For i = 1 To Nesc
                    uesccur = Math.Sqrt(uxesc(i) ^ 2 + uyesc(i) ^ 2)
                    uesccurx = Math.Sqrt(uxesc(i) ^ 2)
                    ii = (uxesc(i) ^ 2 + uyesc(i) ^ 2) / nustepesc
                    nuesc(ii + 1) = nuesc(ii + 1) + 1
                    If (uesccur > umax) Then
                        umax = uesccur
                    End If
                    If (uesccurx > uxmax) Then
                        uxmax = uesccurx
                    End If
                Next i
                ii = 1
                Do Until ii = 200
                    If Ninside > 0 Then
                        If ReplaceEscaped = 0 Then nu(ii) = nu(ii) * Ne / (Ne - Nesc)
                    End If
                    ii = ii + 1
                Loop
                G.FillRectangle(Brushes.Gainsboro, 20, 500, 460, 80)
                Nx = 5
                If Nesc = 0 Then
                    Ny = 3
                Else
                    Ny = 3
                End If
                If ShowEnergyEscaped = 0 Then Nx = Nx * 2
                'Without light, electrons exceeding the workline are not shown 
                If lightreach = 1 Then
                    ndraw = 20 + Nx * velocitycut ^ 2 / nustepesc
                Else
                    ndraw = 20 + Nx * (velocitycut ^ 2 - extractionwork + 1) / nustepesc
                End If
                ii = 1
                Do Until ii = 150
                    If 20 + Nx * ii <= ndraw Then
                        G.FillRectangle(Brushes.Red, 20 + Nx * ii, 580 - Ny * nu(ii), Nx, Ny * nu(ii))
                    End If
                    ii = ii + 1
                Loop
                ii = 1
                If ShowEnergyEscaped = 1 Then
                    Do Until ii = 150
                        If ndraw + Nx * ii + Nx > 480 Then Exit Do
                        G.FillRectangle(Brushes.Green, ndraw + Nx * ii, 580 - Ny * nuesc(ii), Nx, Ny * nuesc(ii))
                        ii = ii + 1
                    Loop
                End If
                'Draw Exstraction Work Line
                ndraw = 20 + Nx + Nx * (velocitycut ^ 2 - extractionwork) / nustepin
                nndraw = 20 + Nx + Nx * (velocitycut ^ 2) / nustepin
                G.DrawLine(Pbb, ndraw, 585, nndraw, 585)
                narrow = 5
                'Left Arrow
                G.DrawLine(Pbb, ndraw, 585, ndraw + narrow, 585 - narrow)
                G.DrawLine(Pbb, ndraw, 585, ndraw + narrow, 585 + narrow)
                'Right Arrow
                G.DrawLine(Pbb, nndraw - narrow, 585 - narrow, nndraw, 585)
                G.DrawLine(Pbb, nndraw - narrow, 585 + narrow, nndraw, 585)
                'Print Exstraction Work Value
                'Dim drawString As String = ("w=" & Format(extractionwork, "###0.0"))
                Dim drawString As String = ("W")
                G.DrawString(drawString, drawFont, Brushes.Black, -10 + (ndraw + nndraw) / 2, 590)
                Dim pbd As New Pen(Color.Black, 1)
                pbd.DashStyle = Drawing2D.DashStyle.DashDot
                G.DrawLine(pbd, ndraw, 580, ndraw, 500)
                G.DrawLine(pbd, nndraw, 580, nndraw, 500)
            End If
        End If
        GoTo 500
        '
        'Draw the fitted curve
        '
        Dim xxdraw As Double
        If Nesc >= 50 Then
            If nustepesc = 0 Then nustepesc = 1
            nnndraw = 20 + Nx * velocitycut ^ 2 / nustepesc
            For xxdraw = 0 To 10 Step 0.1
                ndraw = Ny * f((xxdraw), Rls(1), Rls(2))
                nndraw = Ny * f((xxdraw + 0.1), Rls(1), Rls(2))
                nnx = Nx * (xxdraw) ^ 2 / nustepesc
                nnnx = Nx * (xxdraw + 0.1) ^ 2 / nustepesc
                If (10 + nnndraw + nnnx) > 480 Then Exit For
                G.DrawLine(Por, 5 + nnndraw + nnx, 580 - ndraw, 5 + nnndraw + nnnx, 580 - nndraw)
            Next xxdraw
        End If
        '
        'TIME LAG
        '
500:    i = 1
        Do Until i > 100
            If Nesc >= i And nescprevious < i Then
                Nescapetime(i) = t / 10.0
            End If
            i = i + 1
        Loop
        '
        'FIND ESCAPED ELECTRONS - TIME
        '
        Rls = FindEscapedElectronsTime(Nescapetime, Nesc)
        If Nesc >= 30 And nescprevious < 30 Then
            tlag30 = t / 10.0
        End If
        If Nesc >= nescmax Then
            tesc = tesc + tstep
        End If
        If startlight = 1 And lightreach = 1 Then t = t + tstep
        nescprevious = Nesc
        TextBox3.Text = "t=" & t / 10.0 & " s"
        If Nesc >= nescmax Then
            If Nesc >= nescmax Then
                If UmaxMethod(3) = 1 Then
                    Rls1 = Findumax(uxesc, uyesc, Nesc, 3, G)
                    umax3ii(imode) = Rls1(5) ^ 2 * mulenergy
                    dumax3ii(imode) = 2 * Rls1(5) * Rls1(6) * mulenergy
                End If
                If UmaxMethod(1) = 1 Then
                    Rls1 = Findumax(uxesc, uyesc, Nesc, 1, G)
                    umax1i(icycle) = Rls1(5) ^ 2 * mulenergy
                End If
                If UmaxMethod(2) = 1 Then
                    Rls1 = Findumax(uxesc, uyesc, Nesc, 2, G)
                    umax2ii(imode) = Rls1(5) ^ 2 * mulenergy
                    dumax2ii(imode) = 2 * Rls1(5) * Rls1(6) * mulenergy
                End If
            End If
            tlag11i(icycle) = Nescapetime(1)
            tlag12i(icycle) = -Rls(1) / Rls(3)
            Ici(icycle) = Rls(3)
            If icycle >= 1 Then
                If UmaxMethod(1) = 1 Then
                    MV = MeanValue(umax1i, icycle)
                    umax1ii(imode) = MV(1)
                    dumax1ii(imode) = MV(2)
                End If
                MV = MeanValue(tlag11i, icycle)
                tlag11ii(imode) = MV(1)
                dtlag11ii(imode) = MV(2)
                MV = MeanValue(tlag12i, icycle)
                tlag12ii(imode) = MV(1)
                dtlag12ii(imode) = MV(2)
                MV = MeanValue(Ici, icycle)
                Icii(imode) = MV(1)
                dIcii(imode) = MV(2)
            End If
        End If
        'Compare Methods of umax
        'If (Nesc >= nescmax) And (icycle = Ncycle) Then
        'Dim filename1 As String = "umax.txt"
        'Dim mul As Double = 0.592673898
        'My.Computer.FileSystem.WriteAllText(filename1, Frequency * 10.0 & " " & umax1ii(imode) & " " & dumax1ii(imode) & " " & umax2ii(imode) & " " & dumax2ii(imode) & " " & umax3ii(imode) & " " & dumax3ii(imode) & vbCrLf, True)
        'My.Computer.FileSystem.WriteAllText(filename1, Intensity & " " & umax1ii(imode) & " " & dumax1ii(imode) & " " & umax2ii(imode) & " " & dumax2ii(imode) & " " & umax3ii(imode) & " " & dumax3ii(imode) & vbCrLf, True)
        'End If
        'Compare Methods of tlag
        'If (Nesc >= nescmax) And (icycle = Ncycle) Then
        'Dim filename2 As String = "tlag.txt"
        'My.Computer.FileSystem.WriteAllText(filename2, Frequency * 10.0 & " " & tlag11ii(imode) & " " & dtlag11ii(imode) & " " & tlag12ii(imode) & " " & dtlag12ii(imode) & vbCrLf, True)
        'My.Computer.FileSystem.WriteAllText(filename2, Intensity & " " & 1.0 / Intensity & " " & tlag11ii(imode) & " " & dtlag11ii(imode) & " " & tlag12ii(imode) & " " & dtlag12ii(imode) & vbCrLf, True)
        'End If
        'Record Ic
        'If (Nesc >= nescmax) And (icycle = Ncycle) Then
        'Dim filename3 As String = "Ic.txt"
        'My.Computer.FileSystem.WriteAllText(filename3, Intensity & " " & Icii(imode) & " " & dIcii(imode) & vbCrLf, True)
        'End If
        '
        'DRAW DIAGRAM Escaped e-(Yes-No)-FREQUENCY
        '
        If (CheckBox1.Checked = False Or Nesc >= nescmax) And CheckBox6.Checked = True Then
            If ((Nesc >= 1) Or yesc(ifrequency) = 1) And (ifrequency <= NmodeFrequency) Then
                yesc(ifrequency) = 1
                DrawEscFrequency(Frequencyi, yesc, ifrequency, G)
            Else
                If ifrequency >= 2 Then DrawEscFrequency(Frequencyi, yesc, ifrequency - 1, G)
            End If
        End If
        '
        'DRAW DIAGRAM UMAX-FREQUENCY
        '
        If (CheckBox1.Checked = False Or Nesc >= nescmax) And CheckBox7.Checked = True Then
            If (icycle > 1) Or (icycle = 1 And Nesc >= nescmax) Then DrawUmaxFrequency(Frequencyi, umax1ii, dumax1ii, umax2ii, dumax2ii, umax3ii, dumax3ii, ifrequency, G)
            If (icycle = 1 And Nesc < nescmax) And (ifrequency > 1) Then DrawUmaxFrequency(Frequencyi, umax1ii, dumax1ii, umax2ii, dumax2ii, umax3ii, dumax3ii, ifrequency - 1, G)
        End If
        '
        'DRAW DIAGRAM LAGTIME-FREQUENCY
        '
        If (CheckBox1.Checked = False Or Nesc >= nescmax) And CheckBox8.Checked = True Then
            If (icycle > 1) Or (icycle = 1 And Nesc >= nescmax) Then DrawLagtimeFrequency(Frequencyi, tlag11ii, dtlag11ii, tlag12ii, dtlag12ii, ifrequency, G)
            If (icycle = 1 And Nesc < nescmax) And (ifrequency > 1) Then DrawLagtimeFrequency(Frequencyi, tlag11ii, dtlag11ii, tlag12ii, dtlag12ii, ifrequency - 1, G)
        End If
        '
        'DRAW DIAGRAM ICURRENT-IRRADIANCE
        '
        If (CheckBox1.Checked = False Or Nesc >= nescmax) And CheckBox11.Checked = True Then
            If (icycle > 1) Or (icycle = 1 And Nesc >= nescmax) Then DrawIcIntensity(Intensityi, Icii, dIcii, iintensity, G)
            If (icycle = 1 And Nesc < nescmax) And (iintensity > 1) Then DrawIcIntensity(Intensityi, Icii, dIcii, iintensity - 1, G)
        End If
        '
        'DRAW DIAGRAM UMAX-IRRADIANCE
        '
        If (CheckBox1.Checked = False Or Nesc >= nescmax) And CheckBox9.Checked = True Then
            If (icycle > 1) Or (icycle = 1 And Nesc >= nescmax) Then DrawUmaxIntensity(Intensityi, umax1ii, dumax1ii, umax2ii, dumax2ii, umax3ii, dumax3ii, iintensity, G)
            If (icycle = 1 And Nesc < nescmax) And (iintensity > 1) Then DrawUmaxIntensity(Intensityi, umax1ii, dumax1ii, umax2ii, dumax2ii, umax3ii, dumax3ii, iintensity - 1, G)
        End If
        '
        'DRAW DIAGRAM LAGTIME-IRRADIANCE
        '
        If (CheckBox1.Checked = False Or Nesc >= nescmax) And CheckBox10.Checked = True Then
            If (icycle > 1) Or (icycle = 1 And Nesc >= nescmax) Then DrawLagtimeIntensity(Intensityi, tlag11ii, dtlag11ii, tlag12ii, dtlag12ii, iintensity, G)
            If (icycle = 1 And Nesc < nescmax) And (iintensity > 1) Then DrawLagtimeIntensity(Intensityi, tlag11ii, dtlag11ii, tlag12ii, dtlag12ii, iintensity - 1, G)
        End If
        '
        'DRAW ESCAPED ELECTRONS - TIME DIAGRAM
        '
        If Nesc >= 1 And CheckBox2.Checked = True Then
            If CheckBox1.Checked = False Or Nesc >= nescmax Then DrawEscapedElectronsTime(Nescapetime, Nesc, G)
        End If
        '
        'DRAW THE IMAGE
        '
        If CheckBox1.Checked = False Or Nesc >= nescmax Then
            PicMain.Image = bmpBackBuffer
        End If
        '
        'START UP AGAIN
        '
        If Nesc >= nescmax Then
            If icycle = Ncycle Then
                imode = imode + 1
                For i = 1 To Ncycle
                    Count(i) = 0
                Next i
                '
                'GENERAL STOP
                '
                If imode > Nmode Then
                    'WRITE RESULTS TO FILE
                    WriteToFile()
                    istop = 1
                    Timer1.Enabled = False
                    Timer2.Enabled = False
                    CheckBox1.Enabled = False
                    For i = 1 To 1
                        My.Computer.Audio.PlaySystemSound(Media.SystemSounds.Asterisk)
                    Next i
                End If
                icycle = 0
            End If
            InitialConditions()
            t = 0
            tem = 0
            startlight = 1
            icycle = icycle + 1
        End If
    End Sub
    Private Function Findumax(ByVal ux() As Double, ByVal uy() As Double, _
                           ByVal N As Integer, ByVal Nway As Integer, ByVal G As Graphics) As Double()
        '
        'Fits the ux distribution of escaped electrons and finds uxmax
        'Count is set to 1 when the velocities are added in the current cycle
        '
        Dim i, j, jmax, ii, iimin, iimax, iicut, k, kmax, iround, iroundmax, ndr, ndraw, nndraw As Integer
        Dim uesc(N), nustep, xi(200), yi(200), fi(200), dfdAi(200), dfdBi(200), Nx, Ny As Double
        Dim A, B, S1, S2, S3, S4, S5, D, dA, dB As Double
        Dim Accuracy As Double
        Dim Results(6) As Double
        Accuracy = 0.00000001
        iroundmax = 100
        If icycle = 1 And Count(1) = 0 Then
            sumu = 0
            For ii = 1 To 600
                nuesccycle(ii) = 0
            Next
        End If
        'Find the velocities of escaped electrons
        j = 0
        For i = 1 To N
            j = j + 1
            uesc(j) = Math.Sqrt(ux(i) ^ 2 + uy(i) ^ 2)
            sumu = sumu + uesc(j)
        Next i
        jmax = j
        If Nway = 1 Then
            Results(5) = 0
            For i = 1 To N
                If uesc(i) > Results(5) Then Results(5) = uesc(i)
            Next i
            Results(6) = 0
            Return Results
        End If
        'Form the N-u distribution
        'nustep = 0.1732
        'nustep = 0.5
        'USUAL BIN
        nustep = 1
        'nustep = 0.25
        'nustep = 0.1
        'nustep = 0.05
        If Count(icycle) = 0 Then
            For j = 1 To jmax
                ii = uesc(j) / nustep
                nuesccycle(ii + 1) = nuesccycle(ii + 1) + 1
                'If (ii + 1) > iimax Then iimax = ii + 1
            Next j
            Count(icycle) = 1
        End If
        iimax = 0
        For ii = 1 To 200
            If nuesccycle(ii) > 0 Then iimax = ii
        Next
        'If icycle <> Ncycle Then Return Results
        'The k data points after cutting all the point with ux<=1
        'iicut = (1.0 + nustep) / nustep
        'iicut = (1.0) / nustep
        iicut = (2.0 + nustep) / nustep
        iimin = iicut
        For ii = iimin To iimax
            If nuesccycle(ii + 1) >= nuesccycle(ii) Then
                iimin = ii + 1
            End If
            If nuesccycle(ii + 1) < nuesccycle(ii) Then Exit For
        Next ii
        For ii = iimin To iimax
            k = ii - iimin + 1
            'xi(k) = (ii - 0.5) * nustep
            xi(k) = (ii) * nustep - nustep / 2
            'xi(k) = (ii) * nustep + nustep / 2
            If nuesccycle(ii) = 0 Then Exit For
            yi(k) = nuesccycle(ii)
        Next ii
        kmax = k
        'Find ymax
        S1 = 0
        For k = 1 To kmax
            If yi(k) > S1 Then
                S1 = yi(k)
                ii = k
            End If
        Next k
        S1 = xi(ii)
        S2 = 0
        S3 = 0
        For k = ii To kmax
            S2 = S2 + (xi(k) - xi(ii)) ^ 2 * yi(k)
            S3 = S3 + yi(k)
        Next k
        S4 = Math.Sqrt(S2 / S3)
        Results = LeastSquares(xi, yi, kmax)
        A = Results(1)
        dA = Results(2)
        B = -Results(3)
        dB = Results(4)

        GoTo 100
        G.FillRectangle(Brushes.Gainsboro, 500, 500, 200, 80)
        'Draw the distribution of velocity of escaped electrons
        Ny = yi(1)
        For k = 1 To kmax
            If yi(k) > Ny Then Ny = yi(k)
        Next k
        Nx = 5
        If Ny > 0 Then
            Ny = (icycle * 80 / Ncycle) / Ny
        Else
            Ny = 2
        End If
        For k = 1 To kmax
            nndraw = Ny * yi(k)
            ndraw = Nx * (xi(k) - nustep / 2)
            ndr = Nx * nustep
            G.FillRectangle(Brushes.Green, 500 + ndraw, 580 - nndraw, ndr, nndraw)
        Next k
        'Draw the fitted curve
        Dim nnx, nnnx As Integer
        For xxdraw = 2.0 + nustep To xi(kmax) Step 0.1
            ndraw = Ny * f((xxdraw), Results(1), -Results(3))
            nndraw = Ny * f((xxdraw + 0.1), Results(1), -Results(3))
            nnx = Nx * (xxdraw) / nustep
            nnnx = Nx * (xxdraw + 0.1) / nustep
            G.DrawLine(Pens.Red, 500 + nnx, 580 - ndraw, 500 + nnnx, 580 - nndraw)
        Next xxdraw
100:    If Nway = 2 Then
            Results(5) = -Results(1) / Results(3)
            Results(6) = Results(5) * Math.Sqrt((Results(2) / Results(1)) ^ 2 + (Results(4) / Results(3)) ^ 2)
            Return Results
        End If
        If Nway = 3 Then
            Results(5) = S1 + 2.5 * S4
            Results(6) = S4 / 2
            Return Results
        End If
    End Function
    Private Sub VelocityDistribution(ByVal ux() As Double, ByVal uy() As Double, ByVal N As Integer)
        '
        'Forms and writes to file the E distribution of escaped electrons
        '
        Dim i, j, jmax, ii, iimax As Integer
        Dim uesc(N), nustep, mul As Double
        mul = 0.592673898
        If icycle = 1 Then
            For ii = 1 To 600
                nuesccycleed(ii) = 0
            Next
        End If
        'Find the energies of escaped electrons
        For i = 1 To N
            'uesc(i) = ux(i) ^ 2 + uy(i) ^ 2
            uesc(i) = Math.Sqrt(ux(i) ^ 2 + uy(i) ^ 2)
        Next i
        'Form the N-E distribution
        'nustep = 0.1732
        'nustep = 1.5
        'nustep = 1
        'nustep = 0.25
        nustep = 0.2
        'nustep = 0.05
        For i = 1 To N
            ii = uesc(i) / nustep
            If (ii + 1) <= 600 Then nuesccycleed(ii + 1) = nuesccycleed(ii + 1) + 1
        Next i
        iimax = 0
        For ii = 1 To 600
            If nuesccycleed(ii) > 0 Then iimax = ii
        Next
        'Write nuesccycle to file
        'If icycle = Ncycle Then
        'Dim filename As String
        'filename = "VelocityDistribution.txt"
        'My.Computer.FileSystem.WriteAllText(filename, "f=" & Frequency * 10.0 & "  Ee=" & Intensity & vbCrLf, True)
        'My.Computer.FileSystem.WriteAllText(filename, "ux" & " " & "uy" & vbCrLf, True)
        'My.Computer.FileSystem.WriteAllText(filename, "---------------------------" & vbCrLf, True)
        'For i = 1 To N
        'My.Computer.FileSystem.WriteAllText(filename, ux(i) & " " & uy(i) & vbCrLf, True)
        'Next i
        'My.Computer.FileSystem.WriteAllText(filename, "u" & " " & "N(u)" & vbCrLf, True)
        'My.Computer.FileSystem.WriteAllText(filename, "---------------------------" & vbCrLf, True)
        'For ii = 1 To iimax
        'My.Computer.FileSystem.WriteAllText(filename, Format(mul * (ii * nustep - nustep / 2), "##0.000") & " " & Format(nuesccycleed(ii) / (nustep * mul), "###0.000") & vbCrLf, True)
        'Next ii
        'End If
    End Sub
    Private Function f(ByVal x As Double, ByVal A As Double, ByVal B As Double) As Double
        f = A - B * x
    End Function
    Private Function dfdA(ByVal x As Double, ByVal A As Double, ByVal B As Double) As Double
        dfdA = 1.0
    End Function
    Private Function dfdB(ByVal x As Double, ByVal A As Double, ByVal B As Double) As Double
        dfdB = -x
    End Function
    Private Function dfdx(ByVal x As Double, ByVal A As Double, ByVal B As Double) As Double
        dfdx = (1 - B * x ^ 2) * 2 * A * x * Math.Exp(-B * x ^ 2)
    End Function
    Private Function FindEscapedElectronsTime(ByVal x() As Double, ByVal Nesc As Integer) As Double()
        '
        'FIND ESCAPED ELECTRONS - TIME
        '
        Dim i As Integer
        Dim Rls(5), xls(100), yls(100) As Double
        i = 1
        Do Until i > Nesc
            xls(i) = Nescapetime(i)
            yls(i) = i
            i = i + 1
        Loop
        If Nesc >= 10 Then
            Rls = LeastSquares(xls, yls, Nesc)
        End If
        If Nesc >= 10 Then
            tlag1 = -Rls(1) / Rls(3)
        Else
            tlag1 = Nescapetime(1)
        End If
        Return Rls
    End Function
    Private Sub DrawEscapedElectronsTime(ByVal x() As Double, ByVal Nesc As Integer, ByVal GG As Graphics)
        '
        'DRAW ESCAPED ELECTRONS - TIME DIAGRAM
        '
        Dim i, ndraw, nndraw, Nxstart, Nystart, Narrow, Nxmax, Ny, Noffset, Nxfont, Nyfont As Integer
        Dim Rls(5), xls(100), yls(100), Nx As Double
        Dim Pbb As New Pen(Color.Black, 2)
        Dim Pred As New Pen(Color.Red, 2)
        Dim drawFont As New Font("Arial", 9)
        Dim drawFontscript As New Font("Arial", 7)
        Nxstart = 520
        Nystart = 580
        Nxmax = 130
        Noffset = 5
        Narrow = 5
        Nxfont = 7
        Nyfont = 8
        If Nescapetime(1) > 0 Then
            Nx = 5 * Nescapetime(1) / Nxmax
        Else
            Nx = 1
        End If
        If Nx = 0 Then Nx = 1
        If Nescapetime(Nesc) / Nx > Nxmax Then Nx = Nescapetime(Nesc) / Nxmax
        Ny = 2
        ndraw = Nescapetime(Nesc) / Nx
        GG.FillRectangle(Brushes.Gainsboro, Nxstart, Nystart - Nesc * Ny - Noffset, ndraw + Noffset, Nesc * Ny + Noffset)
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart + ndraw + Noffset, Nystart)
        GG.DrawLine(Pbb, Nxstart + ndraw + Noffset, Nystart, Nxstart + ndraw + Noffset - Narrow, Nystart - Narrow)
        GG.DrawLine(Pbb, Nxstart + ndraw + Noffset, Nystart, Nxstart + ndraw + Noffset - Narrow, Nystart + Narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart, Nystart - Nesc * Ny - Noffset)
        GG.DrawLine(Pbb, Nxstart, Nystart - Nesc * Ny - Noffset, Nxstart - Narrow, Nystart - Nesc * Ny - Noffset + Narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart - Nesc * Ny - Noffset, Nxstart + Narrow, Nystart - Nesc * Ny - Noffset + Narrow)
        GG.DrawLine(Pens.Black, Nxstart, Nystart, Nxstart, Nystart + Narrow)
        GG.DrawLine(Pens.Black, Nxstart - Narrow, Nystart, Nxstart, Nystart)
        GG.DrawString("0", drawFont, Brushes.Black, Nxstart - 2 * Nxfont, Nystart + Narrow)
        GG.DrawString("t(s)", drawFont, Brushes.Black, Nxstart + ndraw + Nxfont, Nystart - Nyfont)
        i = 1
        Do Until i > Nesc
            ndraw = Nescapetime(i) / Nx
            GG.FillEllipse(Brushes.Blue, Nxstart + ndraw - 2, Nystart - Ny * i - 2, 4, 4)
            xls(i) = Nescapetime(i)
            yls(i) = i
            i = i + 1
        Loop
        If Nesc >= 10 Then
            Rls = LeastSquares(xls, yls, Nesc)
            ndraw = (-Rls(1) / Rls(3)) / Nx
            nndraw = ((Nesc - Rls(1)) / Rls(3)) / Nx
            GG.DrawLine(Pred, Nxstart + ndraw, Nystart, Nxstart + nndraw, Nystart - Ny * Nesc)
        End If
        If Nesc >= 10 Then
            tlag1 = -Rls(1) / Rls(3)
        Else
            tlag1 = Nescapetime(1)
        End If
        GG.DrawString("N", drawFont, Brushes.Black, Nxstart - 35, Nystart - Nesc * Ny - Nyfont)
        GG.DrawString("esc", drawFontscript, Brushes.Black, Nxstart - 25, Nystart - Nesc * Ny)
        ndraw = tlag1 / Nx
        GG.DrawString("t", drawFont, Brushes.Black, Nxstart + ndraw - 20, Nystart + Narrow)
        GG.DrawString("lag", drawFontscript, Brushes.Black, Nxstart + ndraw - 15, Nystart + Narrow + 7)
        GG.DrawString("=" & Format(tlag1, "####0") & "s", drawFont, Brushes.Black, Nxstart + ndraw, Nystart + Narrow)
        GG.DrawLine(Pens.Black, Nxstart + ndraw, Nystart, Nxstart + ndraw, Nystart + Narrow)
        ndraw = Nescapetime(Nesc) / Nx
        If Nesc >= 10 Then
            GG.DrawString("Ic=", drawFont, Brushes.Black, Nxstart + ndraw + 10, Nystart - Nesc * Ny - 10)
            GG.DrawString(Format(Rls(3), "##0.000"), drawFont, Brushes.Black, Nxstart + ndraw + 10, Nystart - Nesc * Ny + 5)
            GG.DrawString("μA", drawFont, Brushes.Black, Nxstart + ndraw + 10, Nystart - Nesc * Ny + 20)
        End If
    End Sub
    Private Sub DrawEscFrequency(ByVal x() As Double, ByVal y() As Integer, ByVal N As Integer, ByVal GG As Graphics)
        '
        'DRAW DIAGRAM Esc(Yes-No)-FREQUENCY
        '
        Dim i, ndraw, nndraw, Nxstart, Nystart, Narrow, Nxmax, Nymax, Nyyes, Nyno, Noffset, Nxfont, Nyfont As Integer
        Dim Pbb As New Pen(Color.Black, 2)
        Dim Pbbb As New Pen(Color.Black, 1)
        Dim Pred As New Pen(Color.Red, 2)
        Dim drawFont As New Font("Arial", 10)
        Dim drawFontscript As New Font("Arial", 7)
        Nxstart = 750
        Nystart = 180
        Nxmax = 220
        Nymax = 150
        Nyyes = 120
        Nyno = 20
        Noffset = 20
        Narrow = 5
        Nxfont = 7
        Nyfont = 8
        GG.FillRectangle(Brushes.Gainsboro, Nxstart, Nystart - Nymax, Nxmax, Nymax)
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart + Nxmax, Nystart)
        GG.DrawLine(Pbb, Nxstart + Nxmax, Nystart, Nxstart + Nxmax - Narrow, Nystart - Narrow)
        GG.DrawLine(Pbb, Nxstart + Nxmax, Nystart, Nxstart + Nxmax - Narrow, Nystart + Narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart, Nystart - Nymax)
        GG.DrawLine(Pbb, Nxstart, Nystart - Nymax, Nxstart - Narrow, Nystart - Nymax + Narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart - Nymax, Nxstart + Narrow, Nystart - Nymax + Narrow)
        GG.DrawString("E", drawFont, Brushes.Blue, Nxstart + Nxmax - 70, Nystart - 80)
        GG.DrawString("e", drawFontscript, Brushes.Blue, Nxstart + Nxmax - 59, Nystart - 74)
        GG.DrawString("=" & Format(Intensity, "####0.##"), drawFont, Brushes.Blue, Nxstart + Nxmax - 54, Nystart - 80)
        GG.DrawString(" W/m", drawFont, Brushes.Blue, Nxstart + Nxmax - 13, Nystart - 80)
        GG.DrawString("2", drawFontscript, Brushes.Blue, Nxstart + Nxmax + 27, Nystart - 80)
        GG.DrawString("e- Extraction", drawFont, Brushes.Black, Nxstart - 36, Nystart - Nymax - Noffset)
        GG.DrawString("NO ", drawFont, Brushes.Black, Nxstart - 37, Nystart - Nyno - Nyfont)
        GG.DrawLine(Pens.Black, Nxstart, Nystart - Nyno, Nxstart - Narrow, Nystart - Nyno)
        GG.DrawString("YES", drawFont, Brushes.Black, Nxstart - 39, Nystart - Nyyes - Nyfont)
        GG.DrawLine(Pens.Black, Nxstart, Nystart - Nyyes, Nxstart - Narrow, Nystart - Nyyes)
        GG.DrawString("f(THz)", drawFont, Brushes.Black, Nxstart + Nxmax, Nystart - Nyfont)
        GG.DrawString("0", drawFont, Brushes.Black, Nxstart - Nxfont, Nystart)
        For j = 0 To 1
            nndraw = Nystart + (j - 1) * Nyno - j * Nyyes
            i = 1
            ndraw = 0
            Do Until ndraw >= Nxstart + Nxmax
                ndraw = Nxstart + 2 * i * Narrow
                GG.DrawLine(Pbbb, ndraw - Narrow, nndraw, ndraw, nndraw)
                i = i + 1
            Loop
        Next j
        ndraw = (Nxmax - Noffset)
        GG.DrawString(Format(x(N), "0.####"), drawFont, Brushes.Black, Nxstart + ndraw - Nxfont, Nystart + Narrow)
        GG.DrawLine(Pens.Black, Nxstart + ndraw, Nystart, Nxstart + ndraw, Nystart + Narrow)
        For i = 1 To N
            ndraw = (Nxmax - Noffset) * x(i) / x(N)
            If y(i) = 1 Then nndraw = Nyyes
            If y(i) = 0 Then nndraw = Nyno
            GG.FillEllipse(Brushes.Green, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
        Next i
    End Sub
    Private Sub DrawUmaxFrequency(ByVal x() As Double, ByVal y1() As Double, ByVal dy1() As Double, _
                                  ByVal y2() As Double, ByVal dy2() As Double, ByVal y3() As Double, ByVal dy3() As Double, _
                                  ByVal N As Integer, ByVal GG As Graphics)
        '
        'DRAW DIAGRAM UMAX-FREQUENCY
        '
        Dim i, ndraw, nndraw, nnndraw, Nxstart, Nystart, narrow As Integer
        Dim Nx, Ny, ymax, yimax As Double
        Dim Pbb As New Pen(Color.Black, 2)
        Dim Pred As New Pen(Color.Red, 2)
        Dim drawFont As New Font("Arial", 10)
        Dim drawFontscript As New Font("Arial", 7)
        Nxstart = 750
        Nystart = 400
        ymax = y1(1) + dy1(1)
        yimax = y1(1)
        If UmaxMethod(1) = 1 Then
            For i = 1 To N
                If y1(i) + dy1(i) > ymax Then
                    ymax = y1(i) + dy1(i)
                    yimax = y1(i)
                End If
            Next i
        End If
        If UmaxMethod(2) = 1 Then
            For i = 1 To N
                If y2(i) + dy2(i) > ymax Then
                    ymax = y2(i) + dy2(i)
                    yimax = y2(i)
                End If
            Next i
        End If
        If UmaxMethod(3) = 1 Then
            For i = 1 To N
                If y3(i) + dy3(i) > ymax Then
                    ymax = y3(i) + dy3(i)
                    yimax = y3(i)
                End If
            Next i
        End If
        Nx = 200 / x(N)
        Ny = 180 / ymax
        nndraw = Ny * ymax
        ndraw = Nx * x(N)
        GG.FillRectangle(Brushes.Gainsboro, Nxstart, Nystart - nndraw - 15, ndraw + 15, nndraw + 15)
        narrow = 5
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart + ndraw + 15, Nystart)
        GG.DrawLine(Pbb, Nxstart + ndraw + 15, Nystart, Nxstart + ndraw + 15 - narrow, Nystart - narrow)
        GG.DrawLine(Pbb, Nxstart + ndraw + 15, Nystart, Nxstart + ndraw + 15 - narrow, Nystart + narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart, Nystart - nndraw - 15)
        GG.DrawLine(Pbb, Nxstart, Nystart - nndraw - 15, Nxstart - narrow, Nystart - nndraw - 15 + narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart - nndraw - 15, Nxstart + narrow, Nystart - nndraw - 15 + narrow)
        GG.DrawString("E", drawFont, Brushes.Blue, Nxstart + ndraw - 50, Nystart - 100)
        GG.DrawString("e", drawFontscript, Brushes.Blue, Nxstart + ndraw - 39, Nystart - 94)
        GG.DrawString("=" & Format(Intensity, "####0.##"), drawFont, Brushes.Blue, Nxstart + ndraw - 34, Nystart - 100)
        GG.DrawString(" W/m", drawFont, Brushes.Blue, Nxstart + ndraw + 7, Nystart - 100)
        GG.DrawString("2", drawFontscript, Brushes.Blue, Nxstart + ndraw + 47, Nystart - 100)
        GG.DrawString("0", drawFont, Brushes.Black, Nxstart - 5, Nystart)
        GG.DrawString("f(THz)", drawFont, Brushes.Black, Nxstart + ndraw + 15, Nystart - 15)
        GG.DrawString("K", drawFont, Brushes.Black, Nxstart - 72, Nystart - nndraw - 32)
        GG.DrawString("max", drawFontscript, Brushes.Black, Nxstart - 60, Nystart - nndraw - 25)
        GG.DrawString("(eV)", drawFont, Brushes.Black, Nxstart - 37, Nystart - nndraw - 32)
        GG.DrawString(Format(x(N), "0.####"), drawFont, Brushes.Black, Nxstart + ndraw - 5, Nystart + 5)
        GG.DrawLine(Pens.Black, Nxstart + ndraw, Nystart, Nxstart + ndraw, Nystart + 5)
        nndraw = Ny * yimax
        GG.DrawString(Format(yimax, "####0.0"), drawFont, Brushes.Black, Nxstart - 40, Nystart - nndraw - 10)
        GG.DrawLine(Pens.Black, Nxstart, Nystart - nndraw, Nxstart - 5, Nystart - nndraw)
        If UmaxMethod(1) = 1 Then
            For i = 1 To N
                ndraw = Nx * x(i)
                nndraw = Ny * y1(i)
                nnndraw = Ny * dy1(i)
                GG.DrawLine(Pens.Green, Nxstart + ndraw, Nystart - nndraw + nnndraw, Nxstart + ndraw, Nystart - nndraw - nnndraw)
                GG.DrawLine(Pens.Green, Nxstart + ndraw - 5, Nystart - nndraw + nnndraw, Nxstart + ndraw + 5, Nystart - nndraw + nnndraw)
                GG.DrawLine(Pens.Green, Nxstart + ndraw - 5, Nystart - nndraw - nnndraw, Nxstart + ndraw + 5, Nystart - nndraw - nnndraw)
                GG.FillEllipse(Brushes.Green, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
            Next i
        End If
        If UmaxMethod(2) = 1 Then
            For i = 1 To N
                ndraw = Nx * x(i)
                nndraw = Ny * y2(i)
                nnndraw = Ny * dy2(i)
                GG.DrawLine(Pens.Blue, Nxstart + ndraw, Nystart - nndraw + nnndraw, Nxstart + ndraw, Nystart - nndraw - nnndraw)
                GG.DrawLine(Pens.Blue, Nxstart + ndraw - 5, Nystart - nndraw + nnndraw, Nxstart + ndraw + 5, Nystart - nndraw + nnndraw)
                GG.DrawLine(Pens.Blue, Nxstart + ndraw - 5, Nystart - nndraw - nnndraw, Nxstart + ndraw + 5, Nystart - nndraw - nnndraw)
                GG.FillEllipse(Brushes.Blue, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
            Next i
        End If
        If UmaxMethod(3) = 1 Then
            For i = 1 To N
                ndraw = Nx * x(i)
                nndraw = Ny * y3(i)
                nnndraw = Ny * dy3(i)
                GG.DrawLine(Pens.Red, Nxstart + ndraw, Nystart - nndraw + nnndraw, Nxstart + ndraw, Nystart - nndraw - nnndraw)
                GG.DrawLine(Pens.Red, Nxstart + ndraw - 5, Nystart - nndraw + nnndraw, Nxstart + ndraw + 5, Nystart - nndraw + nnndraw)
                GG.DrawLine(Pens.Red, Nxstart + ndraw - 5, Nystart - nndraw - nnndraw, Nxstart + ndraw + 5, Nystart - nndraw - nnndraw)
                GG.FillEllipse(Brushes.Red, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
            Next i
        End If
    End Sub
    Private Sub DrawLagtimeFrequency(ByVal x() As Double, ByVal y1() As Double, ByVal dy1() As Double, _
                                  ByVal y2() As Double, ByVal dy2() As Double, _
                                  ByVal N As Integer, ByVal GG As Graphics)
        '
        'DRAW DIAGRAM LAGTIME-FREQUENCY
        '
        Dim i, ndraw, nndraw, nnndraw, Nxstart, Nystart, narrow As Integer
        Dim Nx, Ny, ymax, yimax As Double
        Dim Pbb As New Pen(Color.Black, 2)
        Dim Pred As New Pen(Color.Red, 2)
        Dim drawFont As New Font("Arial", 10)
        Dim drawFontscript As New Font("Arial", 7)
        Nxstart = 750
        Nystart = 580
        ymax = y1(1) + dy1(1)
        yimax = y1(1)
        If TlagMethod(1) = 1 Then
            For i = 1 To N
                If y1(i) + dy1(i) > ymax Then
                    ymax = y1(i) + dy1(i)
                    yimax = y1(i)
                End If
            Next i
        End If
        If TlagMethod(2) = 1 Then
            For i = 1 To N
                If y2(i) + dy2(i) > ymax Then
                    ymax = y2(i) + dy2(i)
                    yimax = y2(i)
                End If
            Next i
        End If
        Nx = 200 / x(N)
        Ny = 140 / ymax
        nndraw = Ny * ymax
        ndraw = Nx * x(N)
        GG.FillRectangle(Brushes.Gainsboro, Nxstart, Nystart - nndraw - 15, ndraw + 15, nndraw + 15)
        narrow = 5
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart + ndraw + 15, Nystart)
        GG.DrawLine(Pbb, Nxstart + ndraw + 15, Nystart, Nxstart + ndraw + 15 - narrow, Nystart - narrow)
        GG.DrawLine(Pbb, Nxstart + ndraw + 15, Nystart, Nxstart + ndraw + 15 - narrow, Nystart + narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart, Nystart - nndraw - 15)
        GG.DrawLine(Pbb, Nxstart, Nystart - nndraw - 15, Nxstart - narrow, Nystart - nndraw - 15 + narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart - nndraw - 15, Nxstart + narrow, Nystart - nndraw - 15 + narrow)
        GG.DrawString("E", drawFont, Brushes.Blue, Nxstart + ndraw - 50, Nystart - 30)
        GG.DrawString("e", drawFontscript, Brushes.Blue, Nxstart + ndraw - 39, Nystart - 24)
        GG.DrawString("=" & Format(Intensity, "####0.##"), drawFont, Brushes.Blue, Nxstart + ndraw - 34, Nystart - 30)
        GG.DrawString(" W/m", drawFont, Brushes.Blue, Nxstart + ndraw + 7, Nystart - 30)
        GG.DrawString("2", drawFontscript, Brushes.Blue, Nxstart + ndraw + 47, Nystart - 30)
        GG.DrawString("0", drawFont, Brushes.Black, Nxstart - 5, Nystart)
        GG.DrawString("f(THz)", drawFont, Brushes.Black, Nxstart + ndraw + 15, Nystart - 15)
        GG.DrawString("t", drawFont, Brushes.Black, Nxstart - 46, Nystart - nndraw - 32)
        GG.DrawString("lag", drawFontscript, Brushes.Black, Nxstart - 40, Nystart - nndraw - 25)
        GG.DrawString("(s)", drawFont, Brushes.Black, Nxstart - 25, Nystart - nndraw - 32)
        GG.DrawString(Format(x(N), "0.####"), drawFont, Brushes.Black, Nxstart + ndraw - 5, Nystart + 5)
        GG.DrawLine(Pens.Black, Nxstart + ndraw, Nystart, Nxstart + ndraw, Nystart + 5)
        nndraw = Ny * yimax
        GG.DrawString(Format(yimax, "####0.0"), drawFont, Brushes.Black, Nxstart - 40, Nystart - nndraw - 10)
        GG.DrawLine(Pens.Black, Nxstart, Nystart - nndraw, Nxstart - 5, Nystart - nndraw)
        If TlagMethod(1) = 1 Then
            For i = 1 To N
                ndraw = Nx * x(i)
                nndraw = Ny * y1(i)
                nnndraw = Ny * dy1(i)
                GG.DrawLine(Pens.Green, Nxstart + ndraw, Nystart - nndraw + nnndraw, Nxstart + ndraw, Nystart - nndraw - nnndraw)
                GG.DrawLine(Pens.Green, Nxstart + ndraw - 5, Nystart - nndraw + nnndraw, Nxstart + ndraw + 5, Nystart - nndraw + nnndraw)
                GG.DrawLine(Pens.Green, Nxstart + ndraw - 5, Nystart - nndraw - nnndraw, Nxstart + ndraw + 5, Nystart - nndraw - nnndraw)
                GG.FillEllipse(Brushes.Green, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
            Next i
        End If
        If TlagMethod(2) = 1 Then
            For i = 1 To N
                ndraw = Nx * x(i)
                nndraw = Ny * y2(i)
                nnndraw = Ny * dy2(i)
                GG.DrawLine(Pens.Red, Nxstart + ndraw, Nystart - nndraw + nnndraw, Nxstart + ndraw, Nystart - nndraw - nnndraw)
                GG.DrawLine(Pens.Red, Nxstart + ndraw - 5, Nystart - nndraw + nnndraw, Nxstart + ndraw + 5, Nystart - nndraw + nnndraw)
                GG.DrawLine(Pens.Red, Nxstart + ndraw - 5, Nystart - nndraw - nnndraw, Nxstart + ndraw + 5, Nystart - nndraw - nnndraw)
                GG.FillEllipse(Brushes.Red, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
            Next i
        End If
    End Sub
    Private Sub DrawIcIntensity(ByVal x() As Double, ByVal y() As Double, ByVal dy() As Double, ByVal N As Integer, ByVal GG As Graphics)
        '
        'DRAW DIAGRAM ICURRENT-IRRADIANCE
        '
        Dim i, imax, ndraw, nndraw, nnndraw, Nxstart, Nystart, narrow As Integer
        Dim Nx, Ny, ymax As Double
        Dim Pbb As New Pen(Color.Black, 2)
        Dim Pred As New Pen(Color.Red, 2)
        Dim drawFont As New Font("Arial", 10)
        Dim drawFontscript As New Font("Arial", 7)
        Nxstart = 750
        Nystart = 180
        ymax = y(1) + dy(1)
        imax = 1
        For i = 1 To N
            If y(i) + dy(i) > ymax Then
                ymax = y(i) + dy(i)
                imax = i
            End If
        Next i
        Nx = 200 / x(N)
        Ny = 140 / ymax
        nndraw = Ny * ymax
        ndraw = Nx * x(N)
        GG.FillRectangle(Brushes.Gainsboro, Nxstart, Nystart - nndraw - 15, ndraw + 15, nndraw + 15)
        narrow = 5
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart + ndraw + 15, Nystart)
        GG.DrawLine(Pbb, Nxstart + ndraw + 15, Nystart, Nxstart + ndraw + 15 - narrow, Nystart - narrow)
        GG.DrawLine(Pbb, Nxstart + ndraw + 15, Nystart, Nxstart + ndraw + 15 - narrow, Nystart + narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart, Nystart - nndraw - 15)
        GG.DrawLine(Pbb, Nxstart, Nystart - nndraw - 15, Nxstart - narrow, Nystart - nndraw - 15 + narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart - nndraw - 15, Nxstart + narrow, Nystart - nndraw - 15 + narrow)
        GG.DrawString("f=" & Format(Frequency * 10.0, "###0.###") & " THz", drawFont, Brushes.Blue, Nxstart + ndraw - 50, Nystart - 30)
        GG.DrawString("0", drawFont, Brushes.Black, Nxstart - 5, Nystart)
        GG.DrawString("E", drawFont, Brushes.Black, Nxstart + ndraw + 15, Nystart - 15)
        GG.DrawString("e", drawFontscript, Brushes.Black, Nxstart + ndraw + 26, Nystart - 10)
        GG.DrawString("(W/m", drawFont, Brushes.Black, Nxstart + ndraw + 31, Nystart - 15)
        GG.DrawString("2", drawFontscript, Brushes.Black, Nxstart + ndraw + 72, Nystart - 18)
        GG.DrawString(")", drawFont, Brushes.Black, Nxstart + ndraw + 76, Nystart - 15)
        GG.DrawString("I", drawFont, Brushes.Black, Nxstart - 40, Nystart - nndraw - 32)
        GG.DrawString("c", drawFontscript, Brushes.Black, Nxstart - 35, Nystart - nndraw - 25)
        GG.DrawString("(μΑ)", drawFont, Brushes.Black, Nxstart - 30, Nystart - nndraw - 32)
        GG.DrawString(Format(x(N), "0.####"), drawFont, Brushes.Black, Nxstart + ndraw - 5, Nystart + 5)
        GG.DrawLine(Pens.Black, Nxstart + ndraw, Nystart, Nxstart + ndraw, Nystart + 5)
        nndraw = Ny * y(imax)
        GG.DrawString(Format(y(imax), "###0.0##"), drawFont, Brushes.Black, Nxstart - 48, Nystart - nndraw - 10)
        GG.DrawLine(Pens.Black, Nxstart, Nystart - nndraw, Nxstart - 5, Nystart - nndraw)
        For i = 1 To N
            ndraw = Nx * x(i)
            nndraw = Ny * y(i)
            nnndraw = Ny * dy(i)
            GG.DrawLine(Pbb, Nxstart + ndraw, Nystart - nndraw + nnndraw, Nxstart + ndraw, Nystart - nndraw - nnndraw)
            GG.DrawLine(Pbb, Nxstart + ndraw - 5, Nystart - nndraw + nnndraw, Nxstart + ndraw + 5, Nystart - nndraw + nnndraw)
            GG.DrawLine(Pbb, Nxstart + ndraw - 5, Nystart - nndraw - nnndraw, Nxstart + ndraw + 5, Nystart - nndraw - nnndraw)
            GG.FillEllipse(Brushes.Green, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
        Next i
    End Sub
    Private Sub DrawUmaxIntensity(ByVal x() As Double, ByVal y1() As Double, ByVal dy1() As Double, _
                                  ByVal y2() As Double, ByVal dy2() As Double, ByVal y3() As Double, ByVal dy3() As Double, _
                                  ByVal N As Integer, ByVal GG As Graphics)
        '
        'DRAW DIAGRAM UMAX-IRRADIANCE
        '
        Dim i, ndraw, nndraw, nnndraw, Nxstart, Nystart, narrow As Integer
        Dim Nx, Ny, ymax, yimax As Double
        Dim Pbb As New Pen(Color.Black, 2)
        Dim Pred As New Pen(Color.Red, 2)
        Dim drawFont As New Font("Arial", 10)
        Dim drawFontscript As New Font("Arial", 7)
        Nxstart = 750
        Nystart = 400
        ymax = y1(1) + dy1(1)
        yimax = y1(1)
        If UmaxMethod(1) = 1 Then
            For i = 1 To N
                If y1(i) + dy1(i) > ymax Then
                    ymax = y1(i) + dy1(i)
                    yimax = y1(i)
                End If
            Next i
        End If
        If UmaxMethod(2) = 1 Then
            For i = 1 To N
                If y2(i) + dy2(i) > ymax Then
                    ymax = y2(i) + dy2(i)
                    yimax = y2(i)
                End If
            Next i
        End If
        If UmaxMethod(3) = 1 Then
            For i = 1 To N
                If y3(i) + dy3(i) > ymax Then
                    ymax = y3(i) + dy3(i)
                    yimax = y3(i)
                End If
            Next i
        End If
        Nx = 200 / x(N)
        Ny = 180 / ymax
        nndraw = Ny * ymax
        ndraw = Nx * x(N)
        GG.FillRectangle(Brushes.Gainsboro, Nxstart, Nystart - nndraw - 15, ndraw + 15, nndraw + 15)
        narrow = 5
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart + ndraw + 15, Nystart)
        GG.DrawLine(Pbb, Nxstart + ndraw + 15, Nystart, Nxstart + ndraw + 15 - narrow, Nystart - narrow)
        GG.DrawLine(Pbb, Nxstart + ndraw + 15, Nystart, Nxstart + ndraw + 15 - narrow, Nystart + narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart, Nystart - nndraw - 15)
        GG.DrawLine(Pbb, Nxstart, Nystart - nndraw - 15, Nxstart - narrow, Nystart - nndraw - 15 + narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart - nndraw - 15, Nxstart + narrow, Nystart - nndraw - 15 + narrow)
        GG.DrawString("f=" & Format(Frequency * 10.0, "###0.###") & " THz", drawFont, Brushes.Blue, Nxstart + ndraw - 50, Nystart - 30)
        GG.DrawString("0", drawFont, Brushes.Black, Nxstart - 5, Nystart)
        GG.DrawString("E", drawFont, Brushes.Black, Nxstart + ndraw + 15, Nystart - 15)
        GG.DrawString("e", drawFontscript, Brushes.Black, Nxstart + ndraw + 26, Nystart - 10)
        GG.DrawString("(W/m", drawFont, Brushes.Black, Nxstart + ndraw + 31, Nystart - 15)
        GG.DrawString("2", drawFontscript, Brushes.Black, Nxstart + ndraw + 72, Nystart - 18)
        GG.DrawString(")", drawFont, Brushes.Black, Nxstart + ndraw + 76, Nystart - 15)
        GG.DrawString("K", drawFont, Brushes.Black, Nxstart - 72, Nystart - nndraw - 32)
        GG.DrawString("max", drawFontscript, Brushes.Black, Nxstart - 60, Nystart - nndraw - 25)
        GG.DrawString("(eV)", drawFont, Brushes.Black, Nxstart - 37, Nystart - nndraw - 32)
        GG.DrawString(Format(x(N), "0.####"), drawFont, Brushes.Black, Nxstart + ndraw - 5, Nystart + 5)
        GG.DrawLine(Pens.Black, Nxstart + ndraw, Nystart, Nxstart + ndraw, Nystart + 5)
        nndraw = Ny * yimax
        GG.DrawString(Format(yimax, "####0.0"), drawFont, Brushes.Black, Nxstart - 40, Nystart - nndraw - 10)
        GG.DrawLine(Pens.Black, Nxstart, Nystart - nndraw, Nxstart - 5, Nystart - nndraw)
        If UmaxMethod(1) = 1 Then
            For i = 1 To N
                ndraw = Nx * x(i)
                nndraw = Ny * y1(i)
                nnndraw = Ny * dy1(i)
                GG.DrawLine(Pens.Green, Nxstart + ndraw, Nystart - nndraw + nnndraw, Nxstart + ndraw, Nystart - nndraw - nnndraw)
                GG.DrawLine(Pens.Green, Nxstart + ndraw - 5, Nystart - nndraw + nnndraw, Nxstart + ndraw + 5, Nystart - nndraw + nnndraw)
                GG.DrawLine(Pens.Green, Nxstart + ndraw - 5, Nystart - nndraw - nnndraw, Nxstart + ndraw + 5, Nystart - nndraw - nnndraw)
                GG.FillEllipse(Brushes.Green, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
            Next i
        End If
        If UmaxMethod(2) = 1 Then
            For i = 1 To N
                ndraw = Nx * x(i)
                nndraw = Ny * y2(i)
                nnndraw = Ny * dy2(i)
                GG.DrawLine(Pens.Blue, Nxstart + ndraw, Nystart - nndraw + nnndraw, Nxstart + ndraw, Nystart - nndraw - nnndraw)
                GG.DrawLine(Pens.Blue, Nxstart + ndraw - 5, Nystart - nndraw + nnndraw, Nxstart + ndraw + 5, Nystart - nndraw + nnndraw)
                GG.DrawLine(Pens.Blue, Nxstart + ndraw - 5, Nystart - nndraw - nnndraw, Nxstart + ndraw + 5, Nystart - nndraw - nnndraw)
                GG.FillEllipse(Brushes.Blue, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
            Next i
        End If
        If UmaxMethod(3) = 1 Then
            For i = 1 To N
                ndraw = Nx * x(i)
                nndraw = Ny * y3(i)
                nnndraw = Ny * dy3(i)
                GG.DrawLine(Pens.Red, Nxstart + ndraw, Nystart - nndraw + nnndraw, Nxstart + ndraw, Nystart - nndraw - nnndraw)
                GG.DrawLine(Pens.Red, Nxstart + ndraw - 5, Nystart - nndraw + nnndraw, Nxstart + ndraw + 5, Nystart - nndraw + nnndraw)
                GG.DrawLine(Pens.Red, Nxstart + ndraw - 5, Nystart - nndraw - nnndraw, Nxstart + ndraw + 5, Nystart - nndraw - nnndraw)
                GG.FillEllipse(Brushes.Red, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
            Next i
        End If
    End Sub
    Private Sub DrawLagtimeIntensity(ByVal x() As Double, ByVal y1() As Double, ByVal dy1() As Double, _
                              ByVal y2() As Double, ByVal dy2() As Double, _
                              ByVal N As Integer, ByVal GG As Graphics)
        '
        'DRAW DIAGRAM LAGTIME-IRRADIANCE
        '
        Dim i, ndraw, nndraw, nnndraw, Nxstart, Nystart, narrow As Integer
        Dim Nx, Ny, ymax, yimax As Double
        Dim Pbb As New Pen(Color.Black, 2)
        Dim Pred As New Pen(Color.Red, 2)
        Dim drawFont As New Font("Arial", 10)
        Dim drawFontscript As New Font("Arial", 7)
        Nxstart = 750
        'Nxstart = 1030
        Nystart = 580
        ymax = y1(1) + dy1(1)
        yimax = y1(1)
        If TlagMethod(1) = 1 Then
            For i = 1 To N
                If y1(i) + dy1(i) > ymax Then
                    ymax = y1(i) + dy1(i)
                    yimax = y1(i)
                End If
            Next i
        End If
        If TlagMethod(2) = 1 Then
            For i = 1 To N
                If y2(i) + dy2(i) > ymax Then
                    ymax = y2(i) + dy2(i)
                    yimax = y2(i)
                End If
            Next i
        End If
        Nx = 200 / x(N)
        Ny = 140 / ymax
        nndraw = Ny * ymax
        ndraw = Nx * x(N)
        GG.FillRectangle(Brushes.Gainsboro, Nxstart, Nystart - nndraw - 15, ndraw + 15, nndraw + 15)
        narrow = 5
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart + ndraw + 15, Nystart)
        GG.DrawLine(Pbb, Nxstart + ndraw + 15, Nystart, Nxstart + ndraw + 15 - narrow, Nystart - narrow)
        GG.DrawLine(Pbb, Nxstart + ndraw + 15, Nystart, Nxstart + ndraw + 15 - narrow, Nystart + narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart, Nxstart, Nystart - nndraw - 15)
        GG.DrawLine(Pbb, Nxstart, Nystart - nndraw - 15, Nxstart - narrow, Nystart - nndraw - 15 + narrow)
        GG.DrawLine(Pbb, Nxstart, Nystart - nndraw - 15, Nxstart + narrow, Nystart - nndraw - 15 + narrow)
        GG.DrawString("f=" & Format(Frequency * 10.0, "###0.###") & " THz", drawFont, Brushes.Blue, Nxstart + ndraw - 50, Nystart - 90)
        GG.DrawString("0", drawFont, Brushes.Black, Nxstart - 5, Nystart)
        GG.DrawString("E", drawFont, Brushes.Black, Nxstart + ndraw + 15, Nystart - 15)
        GG.DrawString("e", drawFontscript, Brushes.Black, Nxstart + ndraw + 26, Nystart - 10)
        GG.DrawString("(W/m", drawFont, Brushes.Black, Nxstart + ndraw + 31, Nystart - 15)
        GG.DrawString("2", drawFontscript, Brushes.Black, Nxstart + ndraw + 72, Nystart - 18)
        GG.DrawString(")", drawFont, Brushes.Black, Nxstart + ndraw + 76, Nystart - 15)
        GG.DrawString("t", drawFont, Brushes.Black, Nxstart - 46, Nystart - nndraw - 32)
        GG.DrawString("lag", drawFontscript, Brushes.Black, Nxstart - 40, Nystart - nndraw - 25)
        GG.DrawString("(s)", drawFont, Brushes.Black, Nxstart - 25, Nystart - nndraw - 32)
        GG.DrawString(Format(x(N), "0.####"), drawFont, Brushes.Black, Nxstart + ndraw - 5, Nystart + 5)
        GG.DrawLine(Pens.Black, Nxstart + ndraw, Nystart, Nxstart + ndraw, Nystart + 5)
        nndraw = Ny * yimax
        GG.DrawString(Format(yimax, "####0.0"), drawFont, Brushes.Black, Nxstart - 40, Nystart - nndraw - 10)
        GG.DrawLine(Pens.Black, Nxstart, Nystart - nndraw, Nxstart - 5, Nystart - nndraw)
        If TlagMethod(1) = 1 Then
            For i = 1 To N
                ndraw = Nx * x(i)
                nndraw = Ny * y1(i)
                nnndraw = Ny * dy1(i)
                GG.DrawLine(Pens.Green, Nxstart + ndraw, Nystart - nndraw + nnndraw, Nxstart + ndraw, Nystart - nndraw - nnndraw)
                GG.DrawLine(Pens.Green, Nxstart + ndraw - 5, Nystart - nndraw + nnndraw, Nxstart + ndraw + 5, Nystart - nndraw + nnndraw)
                GG.DrawLine(Pens.Green, Nxstart + ndraw - 5, Nystart - nndraw - nnndraw, Nxstart + ndraw + 5, Nystart - nndraw - nnndraw)
                GG.FillEllipse(Brushes.Green, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
            Next i
        End If
        If TlagMethod(2) = 1 Then
            For i = 1 To N
                ndraw = Nx * x(i)
                nndraw = Ny * y2(i)
                nnndraw = Ny * dy2(i)
                GG.DrawLine(Pens.Red, Nxstart + ndraw, Nystart - nndraw + nnndraw, Nxstart + ndraw, Nystart - nndraw - nnndraw)
                GG.DrawLine(Pens.Red, Nxstart + ndraw - 5, Nystart - nndraw + nnndraw, Nxstart + ndraw + 5, Nystart - nndraw + nnndraw)
                GG.DrawLine(Pens.Red, Nxstart + ndraw - 5, Nystart - nndraw - nnndraw, Nxstart + ndraw + 5, Nystart - nndraw - nnndraw)
                GG.FillEllipse(Brushes.Red, Nxstart + ndraw - 4, Nystart - nndraw - 4, 8, 8)
            Next i
        End If
    End Sub
    Public Function LeastSquares(ByVal x As Double(), ByVal y As Double(), ByVal N As Integer) As Double()
        'Calculates the Least Squares Line y=A+Bx
        'Results(1)=A
        'Results(2)=dA
        'Results(3)=B
        'Results(4)=dB
        'Results(5)=R
        Dim Results(6) As Double
        Dim Sxi, Syi, Sxi2, Sxiyi, D, SyiABxi2, SY As Double
        Dim i As Integer
        i = 1
        Sxi = 0.0
        Syi = 0.0
        Sxi2 = 0.0
        Sxiyi = 0.0
        Do While i <= N
            Sxi = Sxi + x(i)
            Syi = Syi + y(i)
            Sxi2 = Sxi2 + x(i) ^ 2
            Sxiyi = Sxiyi + x(i) * y(i)
            i = i + 1
        Loop
        D = N * Sxi2 - Sxi ^ 2
        Results(1) = (Sxi2 * Syi - Sxi * Sxiyi) / D
        Results(3) = (N * Sxiyi - Sxi * Syi) / D
        i = 1
        SyiABxi2 = 0.0
        Do While i <= N
            SyiABxi2 = SyiABxi2 + (y(i) - Results(1) - Results(3) * x(i)) ^ 2
            i = i + 1
        Loop
        SY = Math.Sqrt(SyiABxi2 / (N - 2))
        Results(2) = SY * Math.Sqrt(Sxi2 / D)
        Results(4) = SY * Math.Sqrt(N / D)
        Results(5) = Math.Sqrt(SyiABxi2)
        Return Results
    End Function
    Public Function LeastSquaresred(ByVal x As Double(), ByVal y As Double(), ByVal N As Integer) As Double
        'Calculates the Least Squares Line y=Bx
        'Results=B
        Dim Results As Double
        Dim Sxi2, Sxiyi As Double
        Dim i As Integer
        i = 1
        Sxi2 = 0.0
        Sxiyi = 0.0
        For i = 1 To N
            Sxi2 = Sxi2 + x(i) ^ 2
            Sxiyi = Sxiyi + x(i) * y(i)
        Next i
        Results = Sxiyi / Sxi2
        Return Results
    End Function
    Public Function MeanValue(ByVal x As Double(), ByVal N As Integer) As Double()
        'Calculates the Mean Value of xi
        'Results(1)=Mean(x)
        'Results(2)=dMean(x)
        Dim Results(2) As Double
        Dim Sxi, Sxixm2 As Double
        Dim i As Integer
        Sxi = 0.0
        For i = 1 To N
            Sxi = Sxi + x(i)
        Next i
        If N >= 1 Then
            Results(1) = Sxi / N
        Else
            Results(1) = 0.0
        End If
        Sxixm2 = 0.0
        For i = 1 To N
            Sxixm2 = Sxixm2 + (x(i) - Results(1)) ^ 2
        Next i
        If N > 1 Then
            Results(2) = Math.Sqrt(Sxixm2 / (N * (N - 1)))
        Else
            Results(2) = 0.0
        End If
        Return Results
    End Function
    Public Sub WriteToFile()
        Dim i As Integer
        Dim umaxii(100), dumaxii(100), tlag1ii(100), dtlag1ii(100) As Double
        Dim FileName As String = "c:\Results.txt"
        My.Computer.FileSystem.WriteAllText(FileName, vbCrLf, True)
        If CheckBox5.Checked = True Then
            My.Computer.FileSystem.WriteAllText(FileName, "CLASSICAL MODEL" & vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "Nrepeat=" & Ncycle & "   Ee=" & Intensity & " W/m2" & vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "f(THz)    " & "escaped e-" & vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "------------------------------" & vbCrLf, True)
            For i = 1 To NmodeFrequency
                If UmaxMethod(1) = 1 Then
                    umaxii(i) = umax1ii(i)
                    dumaxii(i) = dumax1ii(i)
                End If
                If UmaxMethod(2) = 1 Then
                    umaxii(i) = umax2ii(i)
                    dumaxii(i) = dumax2ii(i)
                End If
                If UmaxMethod(3) = 1 Then
                    umaxii(i) = umax3ii(i)
                    dumaxii(i) = dumax3ii(i)
                End If
                If TlagMethod(1) = 1 Then
                    tlag1ii(i) = tlag11ii(i)
                    dtlag1ii(i) = dtlag11ii(i)
                End If
                If TlagMethod(2) = 1 Then
                    tlag1ii(i) = tlag12ii(i)
                    dtlag1ii(i) = dtlag12ii(i)
                End If
            Next i
            For i = 1 To NmodeFrequency
                If yesc(i) = 1 Then
                    My.Computer.FileSystem.WriteAllText(FileName, String.Format("{0,8:###0.000}", Frequencyi(i)) & "  " & "Yes" & vbCrLf, True)
                Else
                    My.Computer.FileSystem.WriteAllText(FileName, String.Format("{0,8:###0.000}", Frequencyi(i)) & "  " & "No" & vbCrLf, True)
                End If
            Next i
            My.Computer.FileSystem.WriteAllText(FileName, vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "f(THz)     " & "Kmax(eV)   " & "dKmax(eV)" & vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "------------------------------" & vbCrLf, True)
            For i = 1 To NmodeFrequency
                My.Computer.FileSystem.WriteAllText(FileName, String.Format("{0,8:###0.000}", Frequencyi(i)) & "  " & String.Format("{0,10:###0.00000}", umaxii(i)) & "  " & String.Format("{0,9:###0.00000}", dumaxii(i)) & vbCrLf, True)
            Next i
            My.Computer.FileSystem.WriteAllText(FileName, vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "f(THz)     " & "tlag(s)   " & "dtlag(s)" & vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "------------------------------" & vbCrLf, True)
            For i = 1 To NmodeFrequency
                My.Computer.FileSystem.WriteAllText(FileName, String.Format("{0,8:###0.000}", Frequencyi(i)) & "  " & String.Format("{0,8:###0.000}", tlag1ii(i)) & "  " & String.Format("{0,8:###0.000}", dtlag1ii(i)) & vbCrLf, True)
            Next i
        Else
            My.Computer.FileSystem.WriteAllText(FileName, "Nrepeat=" & Ncycle & "   f=" & Frequency * 10.0 & " THz" & vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "Ee(W/m2)  " & " Ic(μΑ)     " & "dIc(μΑ)" & vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "------------------------------" & vbCrLf, True)
            For i = 1 To NmodeIntensity
                If UmaxMethod(1) = 1 Then
                    umaxii(i) = umax1ii(i)
                    dumaxii(i) = dumax1ii(i)
                End If
                If UmaxMethod(2) = 1 Then
                    umaxii(i) = umax2ii(i)
                    dumaxii(i) = dumax2ii(i)
                End If
                If UmaxMethod(3) = 1 Then
                    umaxii(i) = umax3ii(i)
                    dumaxii(i) = dumax3ii(i)
                End If
                If TlagMethod(1) = 1 Then
                    tlag1ii(i) = tlag11ii(i)
                    dtlag1ii(i) = dtlag11ii(i)
                End If
                If TlagMethod(2) = 1 Then
                    tlag1ii(i) = tlag12ii(i)
                    dtlag1ii(i) = dtlag12ii(i)
                End If
            Next i
            For i = 1 To NmodeIntensity
                My.Computer.FileSystem.WriteAllText(FileName, String.Format("{0,6:###0.0}", Intensityi(i)) & "  " & String.Format("{0,10:###0.00000}", Icii(i)) & "  " & String.Format("{0,10:###0.00000}", dIcii(i)) & vbCrLf, True)
            Next i
            My.Computer.FileSystem.WriteAllText(FileName, vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "Ee(W/m2) " & " Kmax(eV)   " & "dKmax(eV)" & vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "------------------------------" & vbCrLf, True)
            For i = 1 To NmodeIntensity
                My.Computer.FileSystem.WriteAllText(FileName, String.Format("{0,6:###0.0}", Intensityi(i)) & "  " & String.Format("{0,10:###0.00000}", umaxii(i)) & "  " & String.Format("{0,9:###0.00000}", dumaxii(i)) & vbCrLf, True)
            Next i
            My.Computer.FileSystem.WriteAllText(FileName, vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "Ee(W/m2) " & " tlag(s)   " & "dtlag(s)" & vbCrLf, True)
            My.Computer.FileSystem.WriteAllText(FileName, "------------------------------" & vbCrLf, True)
            For i = 1 To NmodeIntensity
                My.Computer.FileSystem.WriteAllText(FileName, String.Format("{0,6:###0.0}", Intensityi(i)) & "  " & String.Format("{0,8:####0.00}", tlag1ii(i)) & "  " & String.Format("{0,6:##0.00}", dtlag1ii(i)) & vbCrLf, True)
            Next i
        End If
        My.Computer.FileSystem.WriteAllText(FileName, "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" & vbCrLf, True)
    End Sub
    Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
        If istop = 0 Then
            If Button2.Enabled = False Then
                If CheckBox1.Checked = False Then
                    If Timer1.Enabled = True Then
                        Timer1.Enabled = False
                    Else
                        Timer1.Enabled = True
                    End If
                Else
                    If Timer2.Enabled = True Then
                        Timer2.Enabled = False
                    Else
                        Timer2.Enabled = True
                    End If
                End If
            End If
        Else
            Timer1.Enabled = False
            Timer2.Enabled = False
        End If
    End Sub
    Private Sub Button3_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button3.Click
        Timer1.Enabled = False
        Timer2.Enabled = False
        istop = 1
        Button1.Enabled = False
    End Sub
    Private Sub Button4_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button4.Click
        If startlight = 0 Then
            startlight = 1
        Else
            startlight = 0
            lightreach = 0
        End If
        CheckBox1.Enabled = True
    End Sub


End Class