Attribute VB_Name = "Module2"
Option Explicit

'======================================================================
'實現不同方程的牛頓迭代法 和 一些基本復變函數
'======================================================================

Public Const PI As Double = 3.14159265358979
Public Const e  As Double = 2.718281828


'+++++++++++實現不同方程的牛頓迭代法，并返回方程的各種基本性質+++++++++++++
Function MMi(x0 As Double, y0 As Double, nN As Long, M As Long, _
              nM As Long, Hx As Double, Hy As Double, dL1 As Double, _
              dL2 As Double, dL3 As Double, dL4 As Double) As Long
  '調用時MMi(x0, y0, Int(SeData(0, 13)), M, nM, dX, dY, dL1, dL2, dL3, dL4)

    ' 牛頓迭代法解方程原理：
    ' 不失一般性設方程為： f(Z) = 0  （關于復數Z的函數）
    ' 對 f(Z) 求導函數得： f'(Z)
    ' 對任意復數Z0可以有Z1   , Z1 = Z0 - f(Z0)/f'(Z0)
    ' 同樣，對復數Z1可以有Z2 , Z2 = Z1 - f(Z1)/f'(Z1)
    ' …… …… ……
    ' 則，有迭代式：Z(n+1)=Z(n)-f(Z(n))/f'(Z(n))
    ' 對于選定的起始點，迭代大多都會收斂于方程f(z) = 0 的某個根，
    ' 這就是牛頓迭代法解方程的基本方式；
    ' 但也可能存在許多點，使迭代根本就不會收斂，
    ' 甚至可能出現混沌的狀態。
    
    'dX 第一次迭代x軸的變化率
    'dY 第一次迭代y軸的變化率
    'dL1  第一次迭代移動距離
    'dL2  第nM次迭代移動距離
    'dL3  第nM次迭代距離（0，0）點的距離
    'dL4  迭代得到解以后距離解的大概距離
    
    Dim x1 As Double, x2 As Double, y1 As Double, y2 As Double
    Dim m0 As Double, i As Long, k As Long
    Dim SeTa1 As Double, SeTa2 As Double, P1 As Double
    Dim P2 As Double, A As Double, B As Double
           
    Dim temp1 As Double, temp2 As Double, temp3 As Double, temp4 As Double, temp5 As Double
    Dim temp6 As Double, temp7 As Double, temp8 As Double, temp9 As Double, temp0 As Double
    Dim R1 As Double, R2 As Double
    Dim Sum As Long  'Sum值用來判斷迭代是否已經到達根（方程的解），這是一種比較好的方式
    Dim bx As Double, by As Double, xb As Double, yb As Double
           
    Dim N As Long, N7 As Double
    
    On Error GoTo aaa:
    
    x1 = x0: y1 = y0
    N = Int(SeData(0, 14))
    N7 = 0.01
    Sum = 0
    Select Case nN
       Case 1  '方程 Z^N-1=0
           For i = 1 To M
               P1 = Sqr(x1 * x1 + y1 * y1)
               SeTa1 = ZArg(x1, y1, 0)
               x2 = ((N - 1) * P1 * Cos(SeTa1) + P1 ^ (1 - N) * Cos((1 - N) * SeTa1)) / N
               y2 = ((N - 1) * P1 * Sin(SeTa1) + P1 ^ (1 - N) * Sin((1 - N) * SeTa1)) / N
               temp0 = Abs(Abs(x1) - Abs(x2)) ^ 2 + Abs(Abs(y1) - Abs(y2)) ^ 2
               If temp0 < N7 Then
                   Sum = Sum + 1
                   If Sum > 2 Then
                       dL4 = (temp0) / N7
                       If i > nM Then
                           Exit For
                       End If
                   End If
               End If
               If i = 1 Then
                   Hx = Abs(x1 - x2): Hy = Abs(y1 - y2)
                   dL1 = (Hx ^ 2 + Hy ^ 2) ^ 0.5
               End If
               If i = nM Then
                   dL2 = ((x1 - x2) ^ 2 + (y1 - y2) ^ 2) ^ 0.5
                   dL3 = (x2 ^ 2 + y2 ^ 2) ^ 0.5
               End If
               x1 = x2: y1 = y2
           Next i
           
           MMi = (k - 1) * M / N + i
    
       Case 2, 0 '方程 Z^3-1=0 的特解
           x1 = -x1  '其實沒有必要，這里做了一下水平翻轉
           For i = 1 To M
               x2 = x1: y2 = y1: m0 = (x1 * x1 + y1 * y1) ^ 2
               If (x1 * x1 * x1 - x1 * y1 * y1 * 3 - 1) ^ 2 _
                   + (x1 * x1 * y1 * 3 - y1 * y1 * y1) ^ 2 < N7 Then
                       dL4 = Abs(Sqr(x1 * x1 + y1 * y1) - 1) / N7
                       If i > nM Then
                           Exit For
                       End If
               End If
               x1 = x2 - (x2 + y2) / 6 + (x2 * x2 - y2 * y2 - x2 * y2 * 2) / 6 / m0
               y1 = y2 + (x2 - y2) / 6 + (y2 * y2 - x2 * x2 - x2 * y2 * 2) / 6 / m0
               
               
               If i = 1 Then
                   Hx = Abs(x1 - x2): Hy = Abs(y1 - y2)
                   dL1 = (Hx ^ 2 + Hy ^ 2) ^ 0.5
               End If
               If i = nM Then
                   dL2 = ((x1 - x2) ^ 2 + (y1 - y2) ^ 2) ^ 0.5
                   dL3 = (x1 ^ 2 + y1 ^ 2) ^ 0.5
               End If
           Next i
            If x1 > 0.9 Then
                MMi = i + 1 * 0.33 * M
            ElseIf y1 > 0.8 Then
                MMi = i + (0.33 + 1 * 0.33) * M
            ElseIf y1 < -0.8 Then
                MMi = i + (0.66 + 0.33 * 1) * M
            End If
       Case 3  '方程 Z*(1+Z^A)/(1-Z^A)=R
            bx = x0: by = y0
            A = SeData(0, 15) ': B = SeData(0, 16)
            R1 = SeData(0, 17): R2 = SeData(0, 18)
            For i = 1 To M
               temp1 = bx: temp2 = by
               Call Zfang(temp1, temp2, A, temp3, temp4)
               Call Zji(temp3, temp4, temp1, temp2, temp5, temp6)
               Call Zji(temp3, temp4, R1, R2, temp7, temp8)
               temp5 = temp5 + temp7 - R1 + temp1
               temp6 = temp6 + temp8 - R2 + temp2
               Call Zshang(temp3, temp4, temp1, temp2, temp7, temp8)
               temp7 = temp7 + 1 + (A + 1) * temp3
               temp8 = temp8 + (A + 1) * temp4
               Call Zshang(temp5, temp6, temp7, temp8, temp3, temp4)
               bx = temp1 - temp3
               by = temp2 - temp4
               temp0 = Abs(Abs(bx) - Abs(temp1)) ^ 2 + Abs(Abs(by) - Abs(temp2)) ^ 2
               If temp0 < N7 Then
                   Sum = Sum + 1
                   If Sum > 2 Then
                       dL4 = (temp0) / N7
                       If i > nM Then
                           Exit For
                       End If
                   End If
               End If
               If i = 1 Then
                   Hx = Abs(bx - temp1): Hy = Abs(by - temp2)
                   dL1 = (Hx ^ 2 + Hy ^ 2) ^ 0.5
               End If
               If i = nM Then
                   dL2 = ((bx - temp1) ^ 2 + (by - temp2) ^ 2) ^ 0.5
                   dL3 = (bx ^ 2 + by ^ 2) ^ 0.5
               End If
               
            Next i
            MMi = i
       End Select
    Exit Function
aaa:
    MMi = i
End Function



'一些復變函數
'============================

'Z1^Z2  （復數的復數次方）
Public Sub ZZFang(x1 As Double, y1 As Double, x2 As Double, y2 As Double, k As Double, x As Double, y As Double)
    Dim T As Double, TT As Double
    Dim P As Double, Fai As Double
    
    On Error Resume Next
    
    P = Sqr(x1 * x1 + y1 * y1)
    If P = 0 Then
        x = 0: y = 0
        Exit Sub
    End If
    Fai = ZArg(x1, y1, k)
    T = P ^ x2 * e ^ (-y2 * Fai)
    TT = Log(P) * y2 + x2 * Fai
    x = T * Cos(TT)
    y = T * Sin(TT)
End Sub

'Z1*Z2  （復數乘積）
Public Sub Zji(x1 As Double, y1 As Double, x2 As Double, y2 As Double, x As Double, y As Double)
    x = x1 * x2 - y1 * y2
    y = x1 * y2 + y1 * x2
End Sub

'Z1/Z2  （復數商）
Public Sub Zshang(x1 As Double, y1 As Double, x2 As Double, y2 As Double, x As Double, y As Double)
    Dim T As Double
    T = x2 * x2 + y2 * y2
    If T = 0 Then
        If x1 = 0 Then
            x = 1
            y = 1
        Else
            x = Sgn(x1) * 1E+50
            y = Sgn(y1) * 1E+50
        End If
    Else
        x = (x1 * x2 + y1 * y2) / T
        y = (-x1 * y2 + y1 * x2) / T
    End If
End Sub

'Z^N  （復數的實數次方）
Public Sub Zfang(x1 As Double, y1 As Double, N As Double, x As Double, y As Double)
    Dim T As Double, TT As Double, AtnYX As Double
    T = Sqr(x1 * x1 + y1 * y1)
    AtnYX = Atn(y1 / x1)
            If x1 < 0 Then
                TT = PI + AtnYX
            ElseIf y1 > 0 Then
                TT = AtnYX
            Else 'if y1<=0 then
                TT = PI * 2# + AtnYX
            End If
            
    T = T ^ N
    TT = TT * N
    x = T * Cos(TT)
    y = T * Sin(TT)
End Sub

'Arg(Z)  （復數的輻角）
Public Function ZArg(x As Double, y As Double, k As Double) As Double
    If x = 0 Then
        If y = 0 Then
            ZArg = 0
        ElseIf y > 0 Then
            ZArg = PI / 2#
        Else 'if y<0 then
            ZArg = PI * 1.5
        End If
    ElseIf x > 0 Then
        ZArg = Atn(y / x)
        If ZArg < 0 Then ZArg = PI * 2 + ZArg
    Else 'if x<0 then
        ZArg = Atn(y / x) + PI
    End If
    ZArg = ZArg + PI * 2 * k
End Function

'為了實現 Mandelbrot 特效定義的函數 （其實就是Mandelbrot函數迭代）
Public Sub fz2(x0 As Double, y0 As Double, xx As Double, yy As Double, x As Double, y As Double, N As Long)
    Dim t1 As Double, t2 As Double, t3 As Double, t4 As Double
    Dim i As Long
    
    t1 = x0: t2 = y0
    For i = 1 To N
        t3 = t1 * t1 - t2 * t2 + xx
        t4 = 2 * t1 * t2 + yy
        t1 = t3: t2 = t4
    Next i
    x = t1: y = t2
End Sub


'=======================================================