Click here to Skip to main content
15,886,067 members
Articles / Programming Languages / Visual Basic

Mathemathics Framework

Rate me:
Please Sign up or sign in to vote.
4.76/5 (56 votes)
16 Sep 2008CPOL6 min read 75.3K   6.2K   171  
.NET Mathematical Framework
Imports System.Math
Imports System.Drawing
Imports System.Windows.Forms

#Region "Enumeraciones"

Public Enum PointDescriptor
    pNormal
    pPolo
    pCero
    pIndeterminated
    pUnknow
    pAngleRZeroIZero
    pZeroExpZero
    pOverflow
    pBranch
End Enum

Public Enum TypeOfMathError
    ENone                   'Antes de calcular
    NoError                 'despues de calcular
    Warning
    EUnknownError
    EOverFlow
    ENoImplemented
    EDividedForZero
    ECantBeNegative         '
    EThisIsLowerValue       'cuando devuelve obligadamente el minimo posible
    EThisIsBiggerValue      'cuando se ve obligado a devolver el maximo
    EUndeterminate
End Enum

Public Enum ETypeOfComplex
    Binomical
    Polar
End Enum

#End Region

#Region "Estructuras"

#Region "Complejos"

Public Structure ComplexBinomical
    Public Real As Single
    Public Imaginario As Single
    Public PointDescription As PointDescriptor
    Private Shared objMath As BVMathFunctions

    Shared Sub New()
        objMath = New BVMathFunctions
    End Sub

    Public Sub New(ByVal newReal As Single, ByVal newImag As Single)
        objMath = New BVMathFunctions
        Real = newReal
        Imaginario = newImag
    End Sub

    Public Function ToPolar() As ComplexPolar
        Return objMath.ComplexBin2Pol(Me)
    End Function
End Structure

Public Structure ComplexPolar
    Public Modulo As Single
    Public Angulo As Single
    Public PointDescription As PointDescriptor
    Private Shared objMath As BVMathFunctions

    Shared Sub New()
        objMath = New BVMathFunctions
    End Sub

    Public Sub New(ByVal newModulo As Single, ByVal newAngulo As Single)
        objMath = New BVMathFunctions
        Modulo = newModulo
        Angulo = newAngulo
    End Sub

    Public Function ToBinomical() As ComplexBinomical
        Return objMath.ComplexPol2Bin(Me)
    End Function

End Structure

'clase que permite almacenar complejos en forma mas transparente
Public Structure ComplexUndefinied
    Implements ICloneable

    Private intX As Single
    Private intY As Single
    Private intETypeOfComplex As ETypeOfComplex

    Private Shared objMath As BVMathFunctions

#Region "propiedades"

    Public ReadOnly Property TypeOfComplex() As ETypeOfComplex
        Get
            Return intETypeOfComplex
        End Get
    End Property

    Public Property real() As Single
        Get
            If intETypeOfComplex = ETypeOfComplex.Binomical Then
                Return intX
            Else
                Return objMath.ComplexPol2Bin(intX, intY).Real
            End If
        End Get
        Set(ByVal Value As Single)
            If intETypeOfComplex = ETypeOfComplex.Binomical Then
                intX = Value
            Else
                intX = objMath.ComplexModule(Value, intY)
            End If
        End Set
    End Property

    Public Property imaginario() As Single
        Get
            If intETypeOfComplex = ETypeOfComplex.Binomical Then
                Return intY
            Else
                Return objMath.ComplexPol2Bin(intX, intY).Imaginario
            End If
        End Get
        Set(ByVal Value As Single)
            If intETypeOfComplex = ETypeOfComplex.Binomical Then
                intY = Value
            Else
                intY = objMath.ComplexModule(intX, Value)
            End If
        End Set
    End Property

    Public Property modulo() As Single
        Get
            If intETypeOfComplex = ETypeOfComplex.Polar Then
                Return intX
            Else
                Return objMath.ComplexBin2Pol(intX, intY).Modulo
            End If
        End Get
        Set(ByVal Value As Single)
            If intETypeOfComplex = ETypeOfComplex.Polar Then
                intX = Value
            Else
                intX = objMath.ComplexReal(Value, intY)
            End If
        End Set
    End Property

    Public Property angulo() As Single
        Get
            If intETypeOfComplex = ETypeOfComplex.Polar Then
                Return intY
            Else
                Return objMath.ComplexBin2Pol(intX, intY).Angulo
            End If
        End Get
        Set(ByVal Value As Single)
            If intETypeOfComplex = ETypeOfComplex.Polar Then
                intY = Value
            Else
                intY = objMath.ComplexImaginary(intX, Value)
            End If
        End Set
    End Property

    Public Property binomical() As ComplexBinomical
        Get
            If intETypeOfComplex = ETypeOfComplex.Binomical Then
                Return New ComplexBinomical(intX, intY)
            Else
                Return objMath.ComplexPol2Bin(intX, intY)
            End If
        End Get
        Set(ByVal Value As ComplexBinomical)
            If intETypeOfComplex = ETypeOfComplex.Binomical Then
                intX = Value.Real
                intY = Value.Imaginario
            Else
                Dim z As ComplexPolar
                z = objMath.ComplexBin2Pol(Value)
                intX = z.Modulo
                intY = z.Angulo
            End If
        End Set
    End Property

    Public Property polar() As ComplexPolar
        Get
            If intETypeOfComplex = ETypeOfComplex.Polar Then
                Return New ComplexPolar(intX, intY)
            Else
                Return objMath.ComplexBin2Pol(intX, intY)
            End If
        End Get
        Set(ByVal Value As ComplexPolar)
            If intETypeOfComplex = ETypeOfComplex.Polar Then
                intX = Value.Modulo
                intY = Value.Angulo
            Else
                Dim z As ComplexBinomical
                z = objMath.ComplexPol2Bin(Value)
                intX = z.Real
                intY = z.Imaginario
            End If
        End Set
    End Property

#End Region

#Region "Constructores"

    Shared Sub New()
        Dim x1 As ComplexUndefinied

        x1.objMath = New BVMathFunctions
        x1.intX = 0
        x1.intY = 0
        x1.intETypeOfComplex = ETypeOfComplex.Binomical
    End Sub

    Public Sub New(ByVal initReal As Single, ByVal initImaginary As Single)
        objMath = New BVMathFunctions
        intX = initReal
        intY = initImaginary
        intETypeOfComplex = ETypeOfComplex.Binomical
    End Sub

    Public Sub New(ByVal initBinomic As ComplexBinomical)
        objMath = New BVMathFunctions
        intX = initBinomic.Real
        intY = initBinomic.Imaginario
        intETypeOfComplex = ETypeOfComplex.Binomical
    End Sub

    Public Sub New(ByVal initPolar As ComplexPolar)
        objMath = New BVMathFunctions
        intX = initPolar.Modulo
        intY = initPolar.Angulo
        intETypeOfComplex = ETypeOfComplex.Polar
    End Sub

#End Region

#Region "metodos"

    Public Sub ConvertTo(ByVal TypeOfComplex As ETypeOfComplex, ByVal Convert As Boolean)
        If Convert Then
            If intETypeOfComplex <> TypeOfComplex Then
                If TypeOfComplex = ETypeOfComplex.Binomical Then
                    Dim z As ComplexBinomical
                    z = objMath.ComplexPol2Bin(intX, intY)
                    intX = z.Real
                    intY = z.Imaginario
                Else
                    Dim z As ComplexPolar
                    z = objMath.ComplexBin2Pol(intX, intY)
                    intX = z.Modulo
                    intY = z.Angulo
                End If
            End If
        Else
            intETypeOfComplex = TypeOfComplex
        End If
    End Sub

    Public Function Clone() As Object Implements System.ICloneable.Clone
        Dim z As ComplexUndefinied
        If Me.TypeOfComplex = ETypeOfComplex.Binomical Then
            z = New ComplexUndefinied(Me.binomical)
        Else
            z = New ComplexUndefinied(Me.polar)
        End If
        Return z
    End Function

#End Region

End Structure


'Public Structure ComplexBinomicalD
'    Public Real As Double
'    Public Imaginario As Double
'End Structure

'Public Structure ComplexPolarD
'    Public Modulo As Double
'    Public Angulo As Double
'End Structure


#End Region


#Region "Puntos"

'Structure PointFAdv
'    Public X As Single
'    Public Y As Single

'    Public strX As String
'    Public strY As String

'    Private Shared objstrMath As BVStrMath

'    Shared Sub New()
'        objstrMath = New BVStrMath
'    End Sub

'    Overrides Function toString() As String
'        Return "X: " & objstrMath.Multiplos(X, 3, BVStrMath.TypeOfDecimales.TCifrasSignificativas) & _
'            " Y: " & objstrMath.Multiplos(Y, 3, BVStrMath.TypeOfDecimales.TCifrasSignificativas)
'    End Function

'End Structure


#End Region

Public Structure PointComplexUndefinied2D

    Private Z As ComplexUndefinied
    Private X As Single

    Private Shared objMath As BVMathFunctions

#Region "propiedades"

    Public ReadOnly Property TypeOfComplex() As ETypeOfComplex
        Get
            Return Z.TypeOfComplex
        End Get
    End Property

    Public Property real() As Single
        Get
            Return Z.real
        End Get
        Set(ByVal Value As Single)
            Z.real = Value
        End Set
    End Property

    Public Property imaginario() As Single
        Get
            Return Z.imaginario
        End Get
        Set(ByVal Value As Single)
            Z.imaginario = Value
        End Set
    End Property

    Public Property modulo() As Single
        Get
            Return Z.modulo
        End Get
        Set(ByVal Value As Single)
            Z.modulo = Value
        End Set
    End Property

    Public Property angulo() As Single
        Get
            Return Z.angulo
        End Get
        Set(ByVal Value As Single)
            Z.angulo = Value
        End Set
    End Property

    Public Property binomical() As ComplexBinomical
        Get
            Return Z.binomical
        End Get
        Set(ByVal Value As ComplexBinomical)
            Z.binomical = Value
        End Set
    End Property

    Public Property polar() As ComplexPolar
        Get
            Return Z.polar
        End Get
        Set(ByVal Value As ComplexPolar)
            Z.polar = Value
        End Set
    End Property

#End Region

#Region "Constructores"

    Shared Sub New()
        Dim x1 As PointComplexUndefinied2D

        x1.objMath = New BVMathFunctions
        x1.Z.real = 0
        x1.Z.imaginario = 0
        x1.X = 0
    End Sub

    Public Sub New(ByVal initReal As Single, ByVal initImaginary As Single, ByVal initX As Single)
        objMath = New BVMathFunctions
        Z.real = initReal
        Z.imaginario = initImaginary
        X = initX
    End Sub

    Public Sub New(ByVal initBinomic As ComplexBinomical, ByVal initX As Single)
        objMath = New BVMathFunctions
        Z.real = initBinomic.Real
        Z.imaginario = initBinomic.Imaginario
        X = initX
    End Sub

    Public Sub New(ByVal initPolar As ComplexPolar, ByVal initX As Single)
        objMath = New BVMathFunctions
        Z.real = initPolar.Modulo
        Z.imaginario = initPolar.Angulo
        X = initX
    End Sub

    Public Sub New(ByVal initComplexU As ComplexUndefinied, ByVal initX As Single)
        objMath = New BVMathFunctions
        If Z.TypeOfComplex = ETypeOfComplex.Binomical Then
            Z.binomical = initComplexU.binomical
        Else
            Z.polar = initComplexU.polar
        End If
        X = initX
    End Sub


#End Region

#Region "metodos"

    Public Sub ConvertTo(ByVal TypeOfComplex As ETypeOfComplex, ByVal Convert As Boolean)
        Z.ConvertTo(TypeOfComplex, Convert)
    End Sub

#End Region

End Structure


'SIN TERMINAR
Public Structure PointComplexUndefinied3D

    Private Z As ComplexUndefinied
    Private XY As Point

    Private Shared objMath As BVMathFunctions

#Region "propiedades"

    Public ReadOnly Property TypeOfComplex() As ETypeOfComplex
        Get
            Return Z.TypeOfComplex
        End Get
    End Property

    Public Property real() As Single
        Get
            Return Z.real
        End Get
        Set(ByVal Value As Single)
            Z.real = Value
        End Set
    End Property

    Public Property imaginario() As Single
        Get
            Return Z.imaginario
        End Get
        Set(ByVal Value As Single)
            Z.imaginario = Value
        End Set
    End Property

    Public Property modulo() As Single
        Get
            Return Z.modulo
        End Get
        Set(ByVal Value As Single)
            Z.modulo = Value
        End Set
    End Property

    Public Property angulo() As Single
        Get
            Return Z.angulo
        End Get
        Set(ByVal Value As Single)
            Z.angulo = Value
        End Set
    End Property

    Public Property binomical() As ComplexBinomical
        Get
            Return Z.binomical
        End Get
        Set(ByVal Value As ComplexBinomical)
            Z.binomical = Value
        End Set
    End Property

    Public Property polar() As ComplexPolar
        Get
            Return Z.polar
        End Get
        Set(ByVal Value As ComplexPolar)
            Z.polar = Value
        End Set
    End Property

#End Region

#Region "Constructores"

    Shared Sub New()
        Dim x1 As PointComplexUndefinied3D

        x1.objMath = New BVMathFunctions
        x1.Z.real = 0
        x1.Z.imaginario = 0
    End Sub

    Public Sub New(ByVal initReal As Single, ByVal initImaginary As Single)
        objMath = New BVMathFunctions
        Z.real = initReal
        Z.imaginario = initImaginary
    End Sub

    Public Sub New(ByVal initBinomic As ComplexBinomical)
        objMath = New BVMathFunctions
        Z.real = initBinomic.Real
        Z.imaginario = initBinomic.Imaginario
    End Sub

    Public Sub New(ByVal initPolar As ComplexPolar)
        objMath = New BVMathFunctions
        Z.real = initPolar.Modulo
        Z.imaginario = initPolar.Angulo
    End Sub

#End Region

#Region "metodos"

    Public Sub ConvertTo(ByVal TypeOfComplex As ETypeOfComplex, ByVal Convert As Boolean)
        Z.ConvertTo(TypeOfComplex, Convert)
    End Sub

#End Region

End Structure


#End Region

#Region "Type Classes"

'objeto que almacena un numero complejo
'permite obtener rapidamente cualquier valor del objeto guardado
'ya que fueron guardados cuando se actualizo cualquier propiedad
<Obsolete("La clase Complex undefinied ocupa menos memoria")> _
Public Class Complex
    Protected Enum TypeWhatIsUpdated
        TNothing        '
        TBinomical      'se actualizo solo la prte binomica
        TPolar          'solo la polar
        TBoth           'ambas
    End Enum

    Private objBin As ComplexBinomical
    Private objPol As ComplexPolar
    Private Shared objMath As BVMathFunctions
    Protected WhatIsUpdated As TypeWhatIsUpdated

#Region "propiedades"

    Public Property real() As Single
        Get
            If WhatIsUpdated = TypeWhatIsUpdated.TPolar Then
                objBin = objMath.ComplexPol2Bin(objPol)
                WhatIsUpdated = TypeWhatIsUpdated.TBoth
            End If
            Return objBin.Real
        End Get
        Set(ByVal Value As Single)
            objBin.Real = Value
            WhatIsUpdated = TypeWhatIsUpdated.TBinomical
        End Set
    End Property

    Public Property imaginario() As Single
        Get
            If WhatIsUpdated = TypeWhatIsUpdated.TPolar Then
                objBin = objMath.ComplexPol2Bin(objPol)
                WhatIsUpdated = TypeWhatIsUpdated.TBoth
            End If
            Return objBin.Imaginario
        End Get
        Set(ByVal Value As Single)
            objBin.Imaginario = Value
            WhatIsUpdated = TypeWhatIsUpdated.TBinomical
        End Set
    End Property

    Public Property modulo() As Single
        Get
            If WhatIsUpdated = TypeWhatIsUpdated.TBinomical Then
                objPol = objMath.ComplexBin2Pol(objBin)
                WhatIsUpdated = TypeWhatIsUpdated.TBoth
            End If
            Return objPol.Modulo
        End Get
        Set(ByVal Value As Single)
            objPol.Modulo = Value
            WhatIsUpdated = TypeWhatIsUpdated.TPolar
        End Set
    End Property

    Public Property angulo() As Single
        Get
            If WhatIsUpdated = TypeWhatIsUpdated.TBinomical Then
                objPol = objMath.ComplexBin2Pol(objBin)
                WhatIsUpdated = TypeWhatIsUpdated.TBoth
            End If
            Return objPol.Angulo
        End Get
        Set(ByVal Value As Single)
            objPol.Angulo = Value
            WhatIsUpdated = TypeWhatIsUpdated.TPolar
        End Set
    End Property

    Public Property binomical() As ComplexBinomical
        Get
            If WhatIsUpdated = TypeWhatIsUpdated.TPolar Then
                objBin = objMath.ComplexPol2Bin(objPol)
                WhatIsUpdated = TypeWhatIsUpdated.TBoth
            End If
            Return objBin
        End Get
        Set(ByVal Value As ComplexBinomical)
            objBin = Value
            WhatIsUpdated = TypeWhatIsUpdated.TBinomical
        End Set
    End Property

    Public Property polar() As ComplexPolar
        Get
            If WhatIsUpdated = TypeWhatIsUpdated.TBinomical Then
                objPol = objMath.ComplexBin2Pol(objBin)
                WhatIsUpdated = TypeWhatIsUpdated.TBoth
            End If
            Return objPol
        End Get
        Set(ByVal Value As ComplexPolar)
            objPol = Value
            WhatIsUpdated = TypeWhatIsUpdated.TPolar
        End Set
    End Property

#End Region

#Region "Constructores"

    Public Sub New()
        objMath = New BVMathFunctions
        objBin = New ComplexBinomical
        objPol = New ComplexPolar
        WhatIsUpdated = TypeWhatIsUpdated.TNothing
    End Sub

    Public Sub New(ByVal initReal As Single, ByVal initImaginary As Single)
        objMath = New BVMathFunctions
        Me.real = initReal
        Me.imaginario = initImaginary
        WhatIsUpdated = TypeWhatIsUpdated.TBinomical
    End Sub

    Public Sub New(ByVal initBinomic As ComplexBinomical)
        objMath = New BVMathFunctions
        Me.binomical = initBinomic
        WhatIsUpdated = TypeWhatIsUpdated.TBinomical
    End Sub

    Public Sub New(ByVal initPolar As ComplexPolar)
        objMath = New BVMathFunctions
        Me.polar = initPolar
        WhatIsUpdated = TypeWhatIsUpdated.TPolar
    End Sub
#End Region

End Class



#End Region




'Basics Math Functions,Constants...
'faltan implementar funciones como:
'parte real, imaginaria, angulo2modulo, signo
Public NotInheritable Class BVMathFunctions

#Region "Constantes"

    'Public Const Pi = 3.1415926535897931          '180
    Public Const PiSobre2 = PI / 2      '1.5707963267949     '90
    Public Const Pi2 = 2 * PI           '6.28318530717959         '360
    Public Const Pi1 = PI / 180         '0.0174532925199433       '1
    Public Const Pi10 = PI / 18         '0.174532925199433       '0.1
    Public Const Neper = Math.E         ' 2.71828182845905
    Public Const cteLn10 = 2.3025850929940456840179914547D
    Public Const cteLn2 = 0.6931471805599453094172321215D
    Public Const cteInvLn2 = 1.442695040888963407359924681D '1,442695040888963407359924681
    Public Const intOverflow = 1.0E+30
    Public Const intInfinitecimal = 1.0E-30

#End Region

#Region "Variables"
    Public BlnError As TypeOfMathError
    Public PointResult As PointDescriptor
#End Region

    Public Sub New()
        BlnError = TypeOfMathError.ENone
        PointResult = PointDescriptor.pNormal
    End Sub

    'nuevas
#Region "Nuevas"

    'verica si x1 esta dentro del rango x+-Dx
    Public Function CompareRange(ByVal x1 As Single, ByVal x As Single, ByVal Dx As Single) As Boolean
        If x1 <= x + Dx Then
            If x1 >= x - Dx Then
                Return True
            End If
        End If
        Return False
    End Function

    'verica si (x1,y1) esta dentro del circulo con centro en (x,y) y radio Dx
    'por ahora solo verifica dentro del cuadrado!!!!
    Public Function CompareRange(ByVal p1 As PointF, ByVal p As PointF, ByVal r As Single) As Boolean
        If CompareRange(p1.X, p.X, r) Then
            If CompareRange(p1.Y, p.Y, r) Then
                Return True
            End If
        End If
        Return False
    End Function

    Public Function CompareRange(ByVal z1 As ComplexUndefinied, ByVal z As ComplexUndefinied, ByVal r As Single) As Boolean
        If CompareRange(z1.real, z.real, r) Then
            If CompareRange(z1.imaginario, z.imaginario, r) Then
                Return True
            End If
        End If
        Return False
    End Function

#End Region


    'sin terminar
    Public Function complexDerivade(ByVal z1 As ComplexBinomical, ByVal Z1plusDz As ComplexBinomical, ByVal Dz As ComplexBinomical) As ComplexBinomical
        If Dz.Real <> 0 Or Dz.Imaginario <> 0 Then
            If Dz.Real <> 0 Then
                complexDerivade.Real = (z1.Real - Z1plusDz.Real) / Dz.Real
            Else
                complexDerivade.Real = intOverflow
            End If
            If Dz.Imaginario <> 0 Then
                complexDerivade.Imaginario = (z1.Imaginario - Z1plusDz.Imaginario) / Dz.Imaginario
            Else
                complexDerivade.Imaginario = intOverflow
            End If
        Else
            Stop
            Throw New Exception("Dz can't be Zero")
            'complexDerivade.Real = intOverflow
        End If
    End Function

    'sin terminar
    Public Function RealDerivade(ByVal x1 As Single, ByVal x1plusDx As Single, ByVal Dz As ComplexBinomical) As ComplexBinomical
        If Dz.Real <> 0 Or Dz.Imaginario <> 0 Then
            If Dz.Real <> 0 Then
                RealDerivade.Real = x1 / Dz.Real
            Else
                RealDerivade.Real = intOverflow
            End If
            If Dz.Imaginario <> 0 Then
                RealDerivade.Imaginario = (x1 - x1plusDx) / Dz.Imaginario
            Else
                RealDerivade.Imaginario = intOverflow
            End If
        Else
            Stop
            Throw New Exception("Dz can't be Zero")
            'complexDerivade.Real = intOverflow
        End If
    End Function

    'Function Simpson(ByVal a As Double, ByVal b As Double, ByVal N As Integer) As Double
    '    Dim x, x0, x1, x2, h As Double
    '    Dim i As Integer
    '    h = (b - a) / N
    '    x0 = f(a) + f(b)
    '    x1 = 0
    '    x2 = 0
    '    For i = 1 To (N - 1)
    '        x = a + i * h
    '        If Par(i) Then
    '            x2 = x2 + f(x)
    '        Else
    '            x1 = x1 + f(x)
    '        End If
    '    Next i
    '    Simpson = h * (x0 + 2 * x2 + 4 * x1) / 3
    'End Function




#Region "decisiones"


    Public Function IsEqual(ByVal z1 As ComplexBinomical, ByVal z2 As ComplexBinomical) As Boolean
        If z1.Real = z2.Real Then
            If z1.Imaginario = z2.Imaginario Then
                Return True
            Else
                Return False
            End If
        Else
            Return False
        End If
    End Function

    '30/10/04
    Public Function IsEqual(ByVal z1 As ComplexUndefinied, ByVal z2 As ComplexUndefinied) As Boolean
        If z1.real = z2.real Then
            If z1.imaginario = z2.imaginario Then
                Return True
            Else
                Return False
            End If
        Else
            Return False
        End If
    End Function

    Public Function IsModuleMajor(ByVal z1 As ComplexPolar, ByVal z2 As ComplexPolar) As Boolean
        If z1.Modulo > z2.Modulo Then
            Return True
        Else
            Return False
        End If
    End Function

    Public Function IsModuleMinus(ByVal z1 As ComplexPolar, ByVal z2 As ComplexPolar) As Boolean
        If z1.Modulo < z2.Modulo Then
            Return True
        Else
            Return False
        End If

    End Function

    Public Function IsModuleMinusOrEqual(ByVal z1 As ComplexPolar, ByVal z2 As ComplexPolar) As Boolean
        If z1.Modulo <= z2.Modulo Then
            Return True
        Else
            Return False
        End If
    End Function

    Public Function IsModuleMajorOrEqual(ByVal z1 As ComplexPolar, ByVal z2 As ComplexPolar) As Boolean
        If z1.Modulo >= z2.Modulo Then
            Return True
        Else
            Return False
        End If
    End Function



#End Region

#Region "Funciones Especiales"

    Public Function vbCombinatory(ByVal n As Integer, ByVal m As Integer) As Single
        'Stop
        vbCombinatory = 0.5 + Exp(vbFactLn(n) - vbFactLn(m) - vbFactLn(n - m))
        vbCombinatory = Int(vbCombinatory)
    End Function

    Public Function GammaLn(ByVal n As Integer) As Single
        Dim y As Double, x As Double, tmp As Double, ser As Double
        Dim j As Integer
        Static cof(5) As Double

        'cof(0) = Array(76.1800917294715, -86.5053203294168, 24.0140982408309, -1.23173957245015, -1.23173957245015, -5.395239384953E-06)
        cof(0) = 76.1800917294715
        cof(1) = -86.5053203294168
        cof(2) = 24.0140982408309
        cof(3) = -1.23173957245015
        cof(4) = -1.23173957245015
        cof(5) = -0.000005395239384953

        'Stop
        y = n
        x = n
        tmp = x + 5.5
        tmp = tmp - (x + 0.5) * Log(tmp)
        ser = 1.00000000019001
        For j = 0 To 5
            y = y + 1
            ser = ser + cof(j) / y
        Next j
        GammaLn = -tmp + Log(2.506628274631 * ser / x)
    End Function

    Public Function vbFactLn(ByVal n As Integer) As Double
        Static A(101) As Double
        'Stop
        If n < 0 Then
            MsgBox("Negative Factorial in Routine vbFactLn", vbCritical, "Error")
        ElseIf n <= 1 Then
            vbFactLn = 0
        ElseIf n <= 100 Then
            If A(n) <> 0 Then
                vbFactLn = A(n)
            Else
                A(n) = GammaLn(n + 1)
                vbFactLn = A(n)
            End If
        Else
            vbFactLn = GammaLn(n + 1)
        End If
    End Function

    Public Function vbFactorial(ByVal n As Integer) As Single
        Dim j As Integer
        Static ArrayN() As Single = {1, 1, 2, 6, 24}
        Static Ntop As Integer


        If n < 0 Then
            MsgBox("Error", vbCritical)
            Exit Function
        End If
        If n > 32 Then
            vbFactorial = 1
            Exit Function
        End If

        Do While Ntop < n
            j = Ntop
            Ntop = Ntop + 1
            ArrayN(Ntop) = ArrayN(j) * Ntop
        Loop
        vbFactorial = ArrayN(n)
    End Function

#End Region

#Region "Funciones b�sicas de los generadores"

    'Escalon unitario, forma binomica
    Public Function ComplexEsc(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        If Z1.Real >= 0 Then
            ComplexEsc.Real = 1
        Else
            ComplexEsc.Real = 0
        End If
    End Function

    'Escalon unitario, forma polar
    Public Function ComplexEsc(ByVal Z1 As ComplexPolar) As ComplexPolar
        Dim z2 As ComplexPolar

        If Z1.Angulo <= PiSobre2 And Z1.Angulo >= -PiSobre2 Then
            ComplexEsc.Modulo = 1
        ElseIf Z1.Angulo <= PI Or Z1.Angulo >= -PI Then
            ComplexEsc.Modulo = 0
        Else
            z2.Angulo = ShrinkAngle(Z1.Angulo, False)
            z2.Modulo = Z1.Modulo
            ComplexEsc = ComplexEsc(z2)
        End If
    End Function

    'Rampa unitaria
    Public Function ComplexRamp(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        If Z1.Real >= 0 Then
            ComplexRamp.Real = Z1.Real
        Else
            ComplexRamp.Real = 0
        End If
    End Function

    'Rampa unitaria
    Public Function ComplexRamp(ByVal Z1 As ComplexPolar) As ComplexBinomical
        ComplexRamp = ComplexRamp(ComplexPol2Bin(Z1))
    End Function

    'pulso unitario
    'falta arreglar, se debe utilizar el paso de calculo como limites de comparacion 
    Public Function ComplexPulse(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        If Z1.Real > -0.5 And Z1.Real < 0.5 Then
            ComplexPulse.Real = 1
        Else
            ComplexPulse.Real = 0
        End If
    End Function

    'pulso unitario
    'falta arreglar, se debe utilizar el paso de calculo como limites de comparacion 
    Public Function ComplexPulse(ByVal Z1 As ComplexPolar) As ComplexBinomical
        ComplexPulse = ComplexPulse(ComplexPol2Bin(Z1))
    End Function

#End Region

#Region "Utilidades"



    '2/2/04
    Public Function ComplexInverse(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim M2 As Single

        M2 = Z1.Real ^ 2 + Z1.Imaginario ^ 2
        If M2 <> 0 Then
            ComplexInverse.Real = Z1.Real / M2
            ComplexInverse.Imaginario = -Z1.Imaginario / M2
        Else
            ComplexInverse.Real = intOverflow
            ComplexInverse.Imaginario = intOverflow
            ComplexInverse.PointDescription = PointDescriptor.pOverflow
        End If
    End Function

    Public Function ComplexReal(ByVal Z1 As ComplexPolar) As ComplexBinomical
        ComplexReal.Real = Z1.Modulo * Math.Cos(Z1.Angulo)
        ComplexReal.Imaginario = 0
    End Function

    Public Function ComplexReal(ByVal Modulo As Single, ByVal Angulo As Single) As Single
        ComplexReal = Modulo * Math.Cos(Angulo)
    End Function

    Public Function ComplexImaginary(ByVal Modulo As Single, ByVal Angulo As Single) As Single
        ComplexImaginary = Modulo * Math.Sin(Angulo)
    End Function

    'determina el modulo de un complejo
    '16/8/03
    Public Function ComplexModule(ByVal Z1 As ComplexBinomical) As Single
        If Z1.Real = 0 Then
            ComplexModule = Abs(Z1.Imaginario)
        ElseIf Z1.Imaginario = 0 Then
            ComplexModule = Abs(Z1.Real)
        ElseIf Z1.Real > Z1.Imaginario Then
            ComplexModule = Abs(Z1.Real) * Sqrt(1 + (Z1.Imaginario / Z1.Real) ^ 2)
        Else
            ComplexModule = Abs(Z1.Imaginario) * Sqrt(1 + (Z1.Real / Z1.Imaginario) ^ 2)
        End If
    End Function

    'determina el modulo de un complejo
    '16/8/03
    Public Function ComplexModule(ByVal real As Single, ByVal imaginario As Single) As Single
        If real = 0 Then
            ComplexModule = Abs(imaginario)
        ElseIf imaginario = 0 Then
            ComplexModule = Abs(real)
        ElseIf real > imaginario Then
            ComplexModule = Abs(real) * Sqrt(1 + (imaginario / real) ^ 2)
        Else
            ComplexModule = Abs(imaginario) * Sqrt(1 + (real / imaginario) ^ 2)
        End If
    End Function

    'determina el modulo entre dos complejos
    '9/5/05
    Public Function ComplexModule(ByVal Z1 As ComplexBinomical, ByVal Z2 As ComplexBinomical) As Single
        ComplexModule = ComplexModule(ComplexSub(Z1, Z2))
    End Function


    'permite pasar a un numero complejo de la forma binomica a la exponencial
    '16/8/03
    Public Function ComplexBin2Pol(ByVal Z1 As ComplexBinomical) As ComplexPolar
        ComplexBin2Pol.Modulo = ComplexModule(Z1)
        ComplexBin2Pol.Angulo = CalcAngle(Z1.Real, Z1.Imaginario, True)
    End Function

    Public Function ComplexBin2Pol(ByVal real As Single, ByVal imaginario As Single) As ComplexPolar
        ComplexBin2Pol.Modulo = ComplexModule(real, imaginario)
        ComplexBin2Pol.Angulo = CalcAngle(real, imaginario, True)
    End Function


    'permite pasar a un numero complejo de la forma binomica a la exponencial
    '16/8/03
    Public Function ComplexPol2Bin(ByVal Z1 As ComplexPolar) As ComplexBinomical
        ComplexPol2Bin.Real = Z1.Modulo * Math.Cos(Z1.Angulo)
        ComplexPol2Bin.Imaginario = Z1.Modulo * Math.Sin(Z1.Angulo)
    End Function

    Public Function ComplexPol2Bin(ByVal Modulo As Single, ByVal Angulo As Single) As ComplexBinomical
        ComplexPol2Bin.Real = Modulo * Math.Cos(Angulo)
        ComplexPol2Bin.Imaginario = Modulo * Math.Sin(Angulo)
    End Function
#End Region

#Region "Complex Functions"

#Region "Corregidas"

#Region "Trigonometricas"

    '16/8/03
    Public Function ComplexSin(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        If Z1.Imaginario = 0 Then
            ComplexSin.Real = Math.Sin(Z1.Real)
            ComplexSin.Imaginario = 0
        Else
            ComplexSin.Real = Math.Sin(Z1.Real) * Math.Cosh(Z1.Imaginario)
            ComplexSin.Imaginario = Math.Cos(Z1.Real) * Math.Sinh(Z1.Imaginario)
        End If
    End Function

    '21/01/04
    Public Function ComplexSin(ByVal Z1 As ComplexPolar) As ComplexBinomical
        Dim z2 As ComplexBinomical
        z2 = ComplexPol2Bin(Z1)
        ComplexSin = ComplexSin(z2)
    End Function

    '19/1/04
    Public Function ComplexCos(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        If Z1.Imaginario = 0 Then
            ComplexCos.Real = Math.Cos(Z1.Real)
            ComplexCos.Imaginario = 0
        Else
            ComplexCos.Real = Math.Cos(Z1.Real) * Math.Cosh(Z1.Imaginario)
            ComplexCos.Imaginario = -Math.Sin(Z1.Real) * Math.Sinh(Z1.Imaginario)
        End If
    End Function

    '21/01/04
    Public Function ComplexCos(ByVal Z1 As ComplexPolar) As ComplexBinomical
        Dim z2 As ComplexBinomical
        z2 = ComplexPol2Bin(Z1)
        ComplexCos = ComplexCos(z2)
    End Function

    '21/01/04
    'Parece OK
    Public Function ComplexTan(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim cosx As ComplexBinomical, sinx As ComplexBinomical

        If Z1.Imaginario = 0 Then
            cosx.Real = Math.Cos(Z1.Real)
            If cosx.Real = 0 Then
                ComplexTan.Real = intOverflow
                ComplexTan.Imaginario = 0
                ComplexTan.PointDescription = PointDescriptor.pOverflow
            Else
                sinx.Real = Math.Sin(Z1.Real)
                ComplexTan.Real = sinx.Real / cosx.Real
            End If
            ComplexTan.Imaginario = 0
        Else
            cosx = ComplexCos(Z1)
            If cosx.Real = 0 And cosx.Imaginario = 0 Then
                ComplexTan.Real = intOverflow
                ComplexTan.Imaginario = intOverflow
                ComplexTan.PointDescription = PointDescriptor.pOverflow
            Else
                sinx = ComplexSin(Z1)
                ComplexTan = ComplexDivision(sinx, cosx)
            End If
        End If
    End Function

    '21/01/04
    Public Function ComplexTan(ByVal Z1 As ComplexPolar) As ComplexPolar
        Dim z2 As ComplexBinomical
        z2 = ComplexPol2Bin(Z1)
        z2 = ComplexTan(z2)
        ComplexTan = ComplexBin2Pol(z2)
    End Function

    '??? viejo metodo
    Public Function ComplexTanOld(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim dx As Single, cosx As Single, sinx As Single
        Dim coshy As Single, sinhy As Single

        If Z1.Imaginario = 0 Then
            ComplexTanOld.Real = Math.Tan(Z1.Real)
            ComplexTanOld.Imaginario = 0
        Else
            cosx = Math.Cos(Z1.Real)
            sinx = Math.Sin(Z1.Real)
            coshy = Math.Cosh(Z1.Imaginario)
            sinhy = Math.Sinh(Z1.Imaginario)
            dx = (cosx * coshy) ^ 2 + (sinx * sinhy) ^ 2
            ComplexTanOld.Real = cosx * sinx / dx
            ComplexTanOld.Imaginario = coshy * sinhy / dx
        End If
    End Function

#Region "Sin verificar"

    '2/2/04
    Public Function ComplexCosec(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim sinx As ComplexBinomical

        If Z1.Imaginario = 0 Then
            sinx.Real = Math.Sin(Z1.Real)
            If sinx.Real = 0 Then
                ComplexCosec.Real = intOverflow
                ComplexCosec.Imaginario = 0
                ComplexCosec.PointDescription = PointDescriptor.pOverflow
            Else
                ComplexCosec.Real = 1 / sinx.Real
            End If
            ComplexCosec.Imaginario = 0
        Else
            sinx = ComplexSin(Z1)
            If sinx.Real = 0 And sinx.Imaginario = 0 Then
                ComplexCosec.Real = intOverflow
                ComplexCosec.Imaginario = intOverflow
                ComplexCosec.PointDescription = PointDescriptor.pOverflow
            Else
                ComplexCosec = ComplexInverse(sinx)
            End If
        End If
    End Function


    '2/2/04
    Public Function ComplexSec(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim cosx As ComplexBinomical

        If Z1.Imaginario = 0 Then
            cosx.Real = Math.Cos(Z1.Real)
            If cosx.Real = 0 Then
                ComplexSec.Real = intOverflow
                ComplexSec.Imaginario = 0
                ComplexSec.PointDescription = PointDescriptor.pOverflow
            Else
                ComplexSec.Real = 1 / cosx.Real
            End If
            ComplexSec.Imaginario = 0
        Else
            cosx = ComplexCos(Z1)
            If cosx.Real = 0 And cosx.Imaginario = 0 Then
                ComplexSec.Real = intOverflow
                ComplexSec.Imaginario = intOverflow
                ComplexSec.PointDescription = PointDescriptor.pOverflow
            Else
                ComplexSec = ComplexInverse(cosx)
            End If
        End If
    End Function


    '2/2/04
    Public Function ComplexCoTan(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim cosx As ComplexBinomical, sinx As ComplexBinomical

        If Z1.Imaginario = 0 Then
            sinx.Real = Math.Sin(Z1.Real)
            If sinx.Real = 0 Then
                ComplexCoTan.Real = intOverflow
                ComplexCoTan.Imaginario = 0
                ComplexCoTan.PointDescription = PointDescriptor.pOverflow
            Else
                cosx.Real = Math.Cos(Z1.Real)
                ComplexCoTan.Real = cosx.Real / sinx.Real
            End If
            ComplexCoTan.Imaginario = 0
        Else
            sinx = ComplexSin(Z1)
            If sinx.Real = 0 And sinx.Imaginario = 0 Then
                ComplexCoTan.Real = intOverflow
                ComplexCoTan.Imaginario = intOverflow
                ComplexCoTan.PointDescription = PointDescriptor.pOverflow
            Else
                cosx = ComplexCos(Z1)
                ComplexCoTan = ComplexDivision(cosx, sinx)
            End If
        End If
    End Function


#End Region


#End Region

#Region "Hiperbolicas"

    '21/01/04
    'tanh(z)=-i*tan(i*z) 
    'aparentemente OK
    Public Function ComplexTanh(ByVal Z1 As ComplexBinomical) As ComplexBinomical

        If Z1.Imaginario = 0 Then
            ComplexTanh.Real = Math.Sinh(Z1.Real) / Math.Cosh(Z1.Real)
        Else
            ComplexTanh = ComplexDivision(ComplexSinh(Z1), ComplexCosh(Z1))
        End If
    End Function



    '19/11/03
    'lo que hay que hacer
    'Ztmp.Real = -Z1.Imaginario
    'Ztmp.Imaginario = Z1.Real
    'Ztmp2.Real = Sin(Ztmp.Real) * Math.cosh(Ztmp.Imaginario)
    'Ztmp2.Imaginario = Math.cos(Ztmp.Real) * Math.sinh(Ztmp.Imaginario)
    'OK??
    Public Function ComplexSinh(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        If Z1.Imaginario = 0 Then
            ComplexSinh.Real = Math.Sinh(Z1.Real)
            ComplexSinh.Imaginario = 0
        Else
            ComplexSinh.Real = Math.Cos(Z1.Imaginario) * Math.Sinh(Z1.Real)
            ComplexSinh.Imaginario = -Math.Sin(-Z1.Imaginario) * Math.Cosh(Z1.Real)
        End If
    End Function


    'sin(z)/z
    '21/01/04
    'OK, corregido
    Public Function ComplexSinc(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim dx As Single, A As Single, B As Single

        If Z1.Real <> 0 Then
            If Z1.Imaginario <> 0 Then
                A = Math.Sinh(Z1.Imaginario) * Math.Cos(Z1.Real)
                B = Math.Sin(Z1.Real) * Math.Cosh(Z1.Imaginario)
                dx = Z1.Real ^ 2 + Z1.Imaginario ^ 2
                ComplexSinc.Real = (Z1.Real * B + Z1.Imaginario * A) / dx
                ComplexSinc.Imaginario = (Z1.Real * A - Z1.Imaginario * B) / dx
            Else
                'A = 0
                'B = Sin(Z1.Real)
                'dx = Z1.Real ^ 2
                ComplexSinc.Real = Math.Sin(Z1.Real) / Z1.Real
                ComplexSinc.Imaginario = 0
            End If

        ElseIf Z1.Imaginario <> 0 Then
            'A = Math.sinh(Z1.Imaginario)
            'B = 0
            'dx = Z1.Imaginario ^ 2
            ComplexSinc.Real = Math.Sinh(Z1.Imaginario) / Z1.Imaginario
            ComplexSinc.Imaginario = 0
        Else
            ComplexSinc.Real = 1
            ComplexSinc.Imaginario = 0
        End If
    End Function


    '21/01/04
    'OK??? TODO indica que OK
    Public Function ComplexCosh(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        If Z1.Imaginario = 0 Then
            ComplexCosh.Real = Math.Cosh(Z1.Real)
            ComplexCosh.Imaginario = 0
        Else
            ComplexCosh.Real = Math.Cos(Z1.Imaginario) * Math.Cosh(Z1.Real)
            ComplexCosh.Imaginario = Math.Sin(Z1.Imaginario) * Math.Sinh(Z1.Real)
        End If
    End Function

#End Region

#Region "Exponencial"

    'Aparentemente anda bien!!!!!!!
    '21/1/04
    Public Function ComplexExp(ByVal Z1 As ComplexBinomical, ByVal Z2 As ComplexBinomical) As ComplexBinomical
        Dim polTmp As ComplexPolar
        Dim U As Single, V As Single
        Dim M As Single, N As Single

        If Z1.Imaginario = 0 Then
            If Z1.Real = 0 Then
                If Z2.Imaginario = 0 Then
                    If Z2.Real = 0 Then
                        '0^0    lo hago = 1, pero es indeterminado
                        ComplexExp.Real = 1
                        ComplexExp.Imaginario = 1
                        ComplexExp.PointDescription = PointDescriptor.pZeroExpZero
                    Else
                        ComplexExp.Real = 0
                        ComplexExp.Imaginario = 0
                    End If
                Else
                    ComplexExp.Real = 0
                    ComplexExp.Imaginario = 0
                End If

            Else 'Z1=x
                If Z2.Imaginario = 0 Then
                    If Z2.Real = 0 Then
                        ComplexExp.Real = 1
                        ComplexExp.Imaginario = 1
                    Else
                        ComplexExp.Real = Z1.Real ^ Z2.Real
                        ComplexExp.Imaginario = 0
                    End If
                Else
                    If Z2.Imaginario = 0 And Z2.Real = 0 Then
                        'Z1^0
                        ComplexExp.Real = 1
                        ComplexExp.Imaginario = 1
                    Else
                        'completo
                        V = 0.5 * Log(Z1.Real ^ 2 + Z1.Imaginario ^ 2)
                        N = CalcAngle(Z1, True)
                        U = Z2.Real * V - Z2.Imaginario * N
                        If U > 80 Then
                            ComplexExp.Real = intOverflow
                            ComplexExp.Imaginario = intOverflow
                            ComplexExp.PointDescription = PointDescriptor.pOverflow
                            Exit Function
                        Else
                            polTmp.Modulo = Math.Exp(U)
                        End If
                        polTmp.Angulo = Z2.Imaginario * V + Z2.Real * N
                        ComplexExp = ComplexPol2Bin(polTmp)
                    End If
                End If
            End If
        Else
            'Z1=x+jy
            If Z2.Imaginario = 0 And Z2.Real = 0 Then
                'Z1^0
                ComplexExp.Real = 1
                ComplexExp.Imaginario = 1
            Else
                'completo
                V = 0.5 * Math.Log(Z1.Real ^ 2 + Z1.Imaginario ^ 2)
                'N = CalcAngle(Z1, False)
                N = CalcAngle(Z1, True)
                U = Z2.Real * V - Z2.Imaginario * N
                If U > 80 Then
                    ComplexExp.Real = intOverflow
                    ComplexExp.Imaginario = intOverflow
                    ComplexExp.PointDescription = PointDescriptor.pOverflow
                    Exit Function
                Else
                    polTmp.Modulo = Math.Exp(U)
                End If
                polTmp.Angulo = Z2.Imaginario * V + Z2.Real * N
                ComplexExp = ComplexPol2Bin(polTmp)
            End If
        End If
    End Function

#End Region

#Region "Logaritmos"

    'logaritmo natural de un complejo
    '18/01/04 OK!!!
    Public Function ComplexLn(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim Phi, Rho As Single

        If Abs(Z1.Real) > intInfinitecimal Or Abs(Z1.Imaginario) > intInfinitecimal Then
            ComplexLn.Real = 0.5 * Math.Log(Z1.Real ^ 2 + Z1.Imaginario ^ 2)
            ComplexLn.Imaginario = CalcAngle(Z1, True)
        Else
            ComplexLn.Real = -intOverflow
            ComplexLn.PointDescription = PointDescriptor.pOverflow
        End If

    End Function


    'logaritmo natural de un complejo
    '18/01/04  OK!!!!
    Public Function ComplexLn(ByVal Z1 As ComplexPolar) As ComplexBinomical
        If Z1.Modulo > 0 Then
            ComplexLn.Real = Math.Log(Z1.Modulo)
        Else
            ComplexLn.Real = -intOverflow
            ComplexLn.PointDescription = PointDescriptor.pOverflow
        End If
        ComplexLn.Imaginario = Z1.Angulo
    End Function


    '2/2/04
    'sin verificar
    Public Function ComplexLog10(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim Phi, Rho As Single

        If Z1.Real <> 0 Or Z1.Imaginario <> 0 Then
            ComplexLog10.Real = 0.5 * Math.Log(Z1.Real ^ 2 + Z1.Imaginario ^ 2) / cteLn10
        Else
            ComplexLog10.Real = -intOverflow
            ComplexLog10.PointDescription = PointDescriptor.pOverflow
        End If
        ComplexLog10.Imaginario = CalcAngle(Z1, True)

    End Function

    Public Function ComplexLog10(ByVal Z1 As ComplexPolar) As ComplexBinomical
        If Z1.Modulo > 0 Then
            ComplexLog10.Real = Math.Log(Z1.Modulo) / cteLn10
        Else
            ComplexLog10.Real = -intOverflow
            ComplexLog10.PointDescription = PointDescriptor.pOverflow
        End If
        ComplexLog10.Imaginario = Z1.Angulo
    End Function



#End Region

#End Region

#Region "Sin corregir"

    'Z1^Z2, Z1 y Z2 complejos
    'deberia ser igual a: Z^c=e^(c*ln(z)), c complejo y z complejo???
    '7/1/04
    'aparentemente anda mas o menos
    Public Function ComplexExp(ByVal Z1 As ComplexPolar, ByVal Z2 As ComplexBinomical) As ComplexPolar
        Dim U As Single, V As Single

        If Z1.Modulo > 0 Then
            If Z2.Imaginario = 0 Then
                If Z2.Real = 0 Then
                    ComplexExp.Modulo = 1
                Else
                    U = Z2.Real * Log(Z1.Modulo)
                    If U > 80 Then
                        ComplexExp.Modulo = intOverflow
                        ComplexExp.PointDescription = PointDescriptor.pOverflow
                    Else
                        ComplexExp.Modulo = Math.Exp(U)
                    End If
                    ComplexExp.Angulo = Z2.Real * Z1.Angulo
                End If
            Else
                If Z2.Real = 0 Then
                    U = -Z2.Imaginario * Z1.Angulo
                    If U > 80 Then
                        ComplexExp.Modulo = intOverflow
                        ComplexExp.PointDescription = PointDescriptor.pOverflow
                    Else
                        ComplexExp.Modulo = Math.Exp(U)
                    End If
                    ComplexExp.Angulo = Z2.Imaginario * Log(Z1.Modulo)
                Else
                    'forma completa, donde z1 y z2 son complejos con partes real e imaginaria <>0
                    V = Log(Z1.Modulo)
                    U = Z2.Real * V - Z2.Imaginario * Z1.Angulo
                    If U > 80 Then
                        ComplexExp.Modulo = intOverflow
                        ComplexExp.PointDescription = PointDescriptor.pOverflow
                    Else
                        ComplexExp.Modulo = Math.Exp(U)
                    End If
                    ComplexExp.Angulo = Z2.Imaginario * V + Z2.Real * Z1.Angulo
                End If
            End If
        Else 'Z1.Modulo=0
            If Z2.Real <> 0 Or Z2.Imaginario <> 0 Then
                ComplexExp.Modulo = 0
            Else
                ComplexExp.Modulo = 1       'le asigno 1, pero es indeterminado
                ComplexExp.PointDescription = PointDescriptor.pZeroExpZero
            End If
        End If
    End Function



    'Z1^Z2, Z1 y Z2 complejos
    '7/1/04
    '????????????vERIFICAR
    Public Function ComplexExp(ByVal Z1 As ComplexPolar, ByVal Z2 As ComplexPolar) As ComplexPolar
        Stop
        If Z1.Modulo > 0 Then
            If Z2.Modulo > 80 Then
                ComplexExp.Modulo = intOverflow
            Else
                ComplexExp.Modulo = Z1.Modulo ^ Z2.Modulo
            End If
            ComplexExp.Angulo = Z1.Angulo * Z2.Modulo
        Else
            ComplexExp.Modulo = 0
        End If
    End Function



    '19/11/03
    '???? no anda!
    Public Function ComplexArcTan(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim U As ComplexBinomical, V As ComplexBinomical, Ztmp As ComplexBinomical

        U.Real = 1 - Z1.Imaginario
        U.Imaginario = Z1.Real
        V.Real = 1 + Z1.Imaginario
        V.Imaginario = Z1.Real
        Ztmp = ComplexDivision(U, V)
        V = ComplexLn(Ztmp)
        U.Real = 0
        U.Imaginario = -0.5
        ComplexArcTan = ComplexProduct(V, U)
    End Function


    '19/11/03
    'sin verificar ???
    Public Function ComplexArcTanh(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim U As ComplexBinomical, V As ComplexBinomical, Ztmp As ComplexBinomical

        U.Real = 1 + Z1.Real
        U.Imaginario = Z1.Imaginario
        V.Real = 1 - Z1.Real
        V.Imaginario = -Z1.Imaginario
        Ztmp = ComplexDivision(U, V)
        V = ComplexLn(Ztmp)
        U.Real = 0.5
        U.Imaginario = 0
        ComplexArcTanh = ComplexProduct(V, U)
    End Function


#End Region

#End Region

#Region "Operadores Aritmeticos complejos"

    '16/8/03
    Public Function ComplexAdd(ByVal Z1 As ComplexBinomical, ByVal Z2 As ComplexBinomical) As ComplexBinomical
        ComplexAdd.Real = Z1.Real + Z2.Real
        ComplexAdd.Imaginario = Z1.Imaginario + Z2.Imaginario
    End Function


    '21/1/04
    Public Function ComplexAdd(ByVal Z1 As ComplexPolar, ByVal Z2 As ComplexPolar) As ComplexPolar
        Dim Xtmp, X2Tmp As ComplexBinomical

        Xtmp = ComplexPol2Bin(Z1)
        X2Tmp = ComplexPol2Bin(Z2)

        X2Tmp.Real = Xtmp.Real + X2Tmp.Real
        X2Tmp.Imaginario = Xtmp.Imaginario + X2Tmp.Imaginario

        ComplexAdd = ComplexBin2Pol(X2Tmp)
    End Function

    '9/5/04
    Public Function ComplexAdd(ByVal Z1 As ComplexUndefinied, ByVal Z2 As ComplexUndefinied) As ComplexUndefinied
        ComplexAdd.real = Z1.real + Z2.real
        ComplexAdd.imaginario = Z1.imaginario + Z2.imaginario
    End Function

    '9/5/04
    Public Function ComplexSub(ByVal Z1 As ComplexBinomical, ByVal Z2 As ComplexBinomical) As ComplexBinomical
        ComplexSub.Real = Z1.Real - Z2.Real
        ComplexSub.Imaginario = Z1.Imaginario - Z2.Imaginario
    End Function

    '9/5/04
    Public Function ComplexSub(ByVal Z1 As ComplexUndefinied, ByVal Z2 As ComplexUndefinied) As ComplexUndefinied
        ComplexSub.real = Z1.real - Z2.real
        ComplexSub.imaginario = Z1.imaginario - Z2.imaginario
    End Function

    'forma Binomica
    '16/8/03
    Public Function ComplexProduct(ByVal Z1 As ComplexBinomical, ByVal Z2 As ComplexBinomical) As ComplexBinomical
        ComplexProduct.Real = Z1.Real * Z2.Real - Z1.Imaginario * Z2.Imaginario
        ComplexProduct.Imaginario = Z1.Imaginario * Z2.Real + Z1.Real * Z2.Imaginario
    End Function

    'forma polar
    '16/8/03
    Public Function ComplexProduct(ByVal Z1 As ComplexPolar, ByVal Z2 As ComplexPolar) As ComplexPolar
        ComplexProduct.Modulo = Z1.Modulo * Z2.Modulo
        ComplexProduct.Angulo = Z1.Angulo + Z2.Angulo
    End Function

    'forma polar
    '16/8/03
    Public Function ComplexDivision(ByVal Z1 As ComplexPolar, ByVal Z2 As ComplexPolar) As ComplexPolar
        If Z2.Modulo <> 0 Then
            ComplexDivision.Modulo = Z1.Modulo / Z2.Modulo
            ComplexDivision.Angulo = Z1.Angulo - Z2.Angulo
        Else
            ComplexDivision.Modulo = intOverflow
            ComplexDivision.Angulo = intOverflow
        End If
    End Function

    'forma Binomica, z1/z2
    '16/8/03
    Public Function ComplexDivision(ByVal Z1 As ComplexBinomical, ByVal Z2 As ComplexBinomical) As ComplexBinomical
        Dim dc As Single, divisor As Single
        If Z2.Real > Z2.Imaginario Then
            If Z2.Real <> 0 Then
                dc = Z2.Imaginario / Z2.Real
            Else
                dc = intOverflow
            End If
            divisor = Z2.Real + Z2.Imaginario * dc
            If divisor = 0 Then
                ComplexDivision.Real = intOverflow
                ComplexDivision.Imaginario = intOverflow
            Else
                ComplexDivision.Real = (Z1.Real + Z1.Imaginario * dc) / divisor
                ComplexDivision.Imaginario = (Z1.Imaginario - Z1.Real * dc) / divisor
            End If
        Else
            If Z2.Imaginario <> 0 Then
                dc = Z2.Real / Z2.Imaginario
            Else
                dc = intOverflow
            End If
            divisor = Z2.Real * dc + Z2.Imaginario
            If divisor = 0 Then
                ComplexDivision.Real = intOverflow
                ComplexDivision.Imaginario = intOverflow
            Else
                ComplexDivision.Real = (Z1.Real * dc + Z1.Imaginario) / divisor
                ComplexDivision.Imaginario = (Z1.Imaginario * dc - Z1.Real) / divisor
            End If
        End If
    End Function


    '21/1/04
    Public Function ComplexSqr(ByVal Z1 As ComplexPolar) As ComplexPolar
        ComplexSqr.Modulo = Math.Sqrt(Z1.Modulo)
        ComplexSqr.Angulo = Z1.Angulo / 2
    End Function


    '16/8/03
    Public Function ComplexSqr(ByVal Z1 As ComplexBinomical) As ComplexBinomical
        Dim w As Single, cd As Single

        If Z1.Real = 0 And Z1.Imaginario = 0 Then
            w = 0
        ElseIf Abs(Z1.Real) >= Abs(Z1.Imaginario) Then
            cd = Z1.Imaginario / Z1.Real
            w = Sqrt(Abs(Z1.Real)) * Sqrt((1 + Sqrt(1 + cd ^ 2)) / 2)
        Else
            cd = Z1.Real / Z1.Imaginario
            w = Sqrt(Abs(Z1.Imaginario)) * Sqrt((Abs(cd) + Sqrt(1 + cd ^ 2)) / 2)
        End If

        If w = 0 Then
            ComplexSqr.Real = 0
            ComplexSqr.Imaginario = 0
        ElseIf Z1.Real >= 0 Then
            ComplexSqr.Real = w
            ComplexSqr.Imaginario = Z1.Imaginario / (2 * w)
        ElseIf Z1.Real < 0 And Z1.Imaginario >= 0 Then
            ComplexSqr.Real = Abs(Z1.Imaginario) / (2 * w)
            ComplexSqr.Imaginario = w
        Else 'If Z1.Imaginario >= 0 Then
            ComplexSqr.Real = Abs(Z1.Imaginario) / (2 * w)
            ComplexSqr.Imaginario = -w
            'Else
            '   Stop
        End If
    End Function

#End Region

#Region "Funciones para vectores (incluye complejos)"

    Public Function Modulo(ByVal SValue As ComplexBinomical) As Single
        Modulo = Sqrt(SValue.Real ^ 2 + SValue.Imaginario ^ 2)
    End Function

    Public Function Modulo(ByVal Real As Single, ByVal Imag As Single) As Single
        Modulo = Sqrt(Real ^ 2 + Imag ^ 2)
    End Function

    Public Function CalcAngle(ByVal Real As Single, ByVal Imag As Single, ByVal Bln360 As Boolean) As Single
        Dim Angle As Single
        'calcula el angulo de un numero complejo, lo devuelve en radianes

        If Real <> 0 Then
            Angle = Math.Atan(Abs(Imag / Real))
        ElseIf Imag <> 0 Then
            Angle = PI / 2
        Else
            'indeterminado
            CalcAngle = 0
            BlnError = TypeOfMathError.EUndeterminate
            PointResult = PointDescriptor.pAngleRZeroIZero
            Exit Function
        End If

        If Bln360 Then
            'devuelve un angulo entre 0 y 360
            If Real >= 0 And Imag >= 0 Then
                CalcAngle = Angle
            ElseIf Real < 0 And Imag >= 0 Then
                CalcAngle = PI - Angle
            ElseIf Real >= 0 And Imag < 0 Then
                CalcAngle = Pi2 - Angle
            ElseIf Real < 0 And Imag < 0 Then
                CalcAngle = Angle + PI
            End If

        Else
            'devuelve un angulo entre +-180
            If Real >= 0 And Imag >= 0 Then
                CalcAngle = Angle
            ElseIf Real < 0 And Imag >= 0 Then
                CalcAngle = PI - Angle
            ElseIf Real >= 0 And Imag < 0 Then
                CalcAngle = -Angle
            ElseIf Real < 0 And Imag < 0 Then
                CalcAngle = Angle - PI
            End If
        End If
    End Function

    Public Function CalcAngle(ByVal SValue As ComplexBinomical, ByVal Bln360 As Boolean) As Single
        Dim Angle As Single
        'calcula el angulo de un numero complejo, lo devuelve en radianes

        If SValue.Real <> 0 Then
            Angle = Math.Atan(Abs(SValue.Imaginario / SValue.Real))
        Else
            Angle = PI / 2
        End If

        If Bln360 Then
            'devuelve un angulo entre 0 y 360
            If SValue.Real >= 0 And SValue.Imaginario >= 0 Then
                CalcAngle = Angle
            ElseIf SValue.Real < 0 And SValue.Imaginario >= 0 Then
                CalcAngle = PI - Angle
            ElseIf SValue.Real >= 0 And SValue.Imaginario < 0 Then
                CalcAngle = Pi2 - Angle
            ElseIf SValue.Real < 0 And SValue.Imaginario < 0 Then
                CalcAngle = Angle + PI
            End If

        Else
            'devuelve un angulo entre +-180
            If SValue.Real >= 0 And SValue.Imaginario >= 0 Then
                CalcAngle = Angle
            ElseIf SValue.Real < 0 And SValue.Imaginario >= 0 Then
                CalcAngle = PI - Angle
            ElseIf SValue.Real >= 0 And SValue.Imaginario < 0 Then
                CalcAngle = -Angle
            ElseIf SValue.Real < 0 And SValue.Imaginario < 0 Then
                CalcAngle = Angle - PI
            End If
        End If
    End Function

    Public Function Degree2Rad(ByVal DegreeAngle As Single, ByVal BlnBetween360 As Boolean) As Single
        Dim Angle2 As Single, n As Integer

        Angle2 = ShrinkAngleDegree(DegreeAngle, True)
        Degree2Rad = Angle2 * PI / 360
    End Function

    'convirte de radiane a grados
    Public Function Rad2Degree(ByVal RadAngle As Single, ByVal BlnBetween360 As Boolean) As Single
        Dim Angle2 As Single, n As Integer

        Angle2 = ShrinkAngle(RadAngle, True)
        Rad2Degree = Angle2 * 360 / PI
    End Function

    'achica el angulo de un numero entre +-pi
    Public Function ShrinkAngleDegree(ByVal DegreeAngle As Single, ByVal BlnBetween360 As Boolean) As Single
        Dim Angle2 As Single, n As Integer

        If Abs(DegreeAngle) > 360 Then
            n = Fix(DegreeAngle / 360)
            ShrinkAngleDegree = DegreeAngle - (n * 360)
        Else
            ShrinkAngleDegree = DegreeAngle
        End If
        If BlnBetween360 = False Then
            If DegreeAngle > 180 Then
                ShrinkAngleDegree = DegreeAngle - 180
            ElseIf DegreeAngle < -180 Then
                ShrinkAngleDegree = DegreeAngle + 180
            End If
        End If
    End Function

    'achica el angulo de un numero entre +-pi
    Public Function ShrinkAngle(ByVal Angle As Single, ByVal BlnBetween360 As Boolean) As Single
        Dim Angle2 As Single, n As Integer

        If Abs(Angle) > Pi2 Then
            n = Fix(Angle / Pi2)
            ShrinkAngle = Angle - (n * Pi2)
        Else
            ShrinkAngle = Angle
        End If
        If BlnBetween360 = False Then
            If Angle > PI Then
                ShrinkAngle = Angle - PI
            ElseIf Angle < -PI Then
                ShrinkAngle = Angle + PI
                'Else
                '    ShrinkAngle = Angle
            End If
        End If
    End Function


#End Region

#Region "Funciones Logaritmicas"

    Public Function Log10(ByVal Valor As Single) As Single
        If Valor > 0 Then
            Log10 = Log(Valor) / cteLn10
        Else
            Log10 = -intOverflow
            BlnError = TypeOfMathError.ECantBeNegative
        End If
    End Function

    'logaritmo en base 2
    Public Function Log2(ByVal Value As Single) As Single
        If Value > 0 Then
            Log2 = Log(Value) * cteInvLn2
        Else
            Log2 = -intOverflow
            BlnError = TypeOfMathError.ECantBeNegative
        End If
    End Function

    Public Function DeciBell(ByVal Value As Single) As Single
        If Value > 0 Then
            'Dim SBin As ComplexBinomical
            'DeciBell = Int(2000 * Log(Value) / Log(10)) / 100
            DeciBell = 20 * Math.Log10(Value)
        Else
            DeciBell = -intOverflow
            BlnError = TypeOfMathError.ECantBeNegative
        End If
    End Function

#End Region

#Region "Miscelaneus Methods"

    'redondea un numero con la cantidad de digitos indicados
    'OK para valores menores que 1
    Public Function RoundNumber(ByVal Valor As Single, ByVal NroDec As Integer) As Single
        Dim intTmp As Single, intTmp2 As Single, ValorTmp As Single
        Dim StrTmpa As String

        If Valor = 0 Then
            RoundNumber = 0
        Else
            ValorTmp = Math.Abs(Valor)

            intTmp = 10 ^ NroDec
            intTmp2 = Log10(ValorTmp)
            If intTmp2 < 0 Then
                'para numeros < 1
                intTmp2 = -Fix(intTmp2)
                intTmp2 = 10 ^ intTmp2
                RoundNumber = (Int(ValorTmp * intTmp2 * intTmp) / intTmp2) / intTmp
                'ElseIf valorTmp >= 1 And valorTmp < 2 Then
                '    RoundNumber = valorTmp '\ intTmp2) * intTmp2
            Else
                'numeros >1
                intTmp2 = Int(intTmp2)
                intTmp2 = 10 ^ intTmp2

                'el int() redondea mal algunos numeros!!??
                RoundNumber = Int((ValorTmp * intTmp) / intTmp2) * intTmp2 / intTmp
                'StrTmpa = ValorTmp * intTmp / intTmp2
                'RoundNumber = InStr(1, StrTmpa, ".", vbTextCompare)
                'If RoundNumber = 0 Then
                'RoundNumber = ValorTmp
                'Else
                '   StrTmpa = Left((ValorTmp * intTmp) / intTmp2, RoundNumber) * intTmp2 / intTmp
                '   RoundNumber = StrTmpa
                'End If
            End If

            If Valor < 0 Then
                RoundNumber = -RoundNumber
            End If
        End If
    End Function

    'retorna el porcentaje(modulo 1) del valor ingresado dentro de la decada 
    'en la que se halla, por ej: v=2,5*E3
    'LimInferior=1E3,   LimSuperior=1E4,   BetweenDecadePercent=log10(2,5)=0,4 (40%)
    '19/7/04
    Public Function BetweenDecadePercent(ByVal v As Single) As Single
        If v > 0 Then
            Dim m, n As Single

            m = Log10(v)
            n = Math.Floor(m)
            Return m - n
        Else
            Stop
        End If
    End Function

    'recorta un numero con la cantidad de digitos decimales indicados
    Public Function Decimales(ByVal Valor As Single, ByVal NroDec As Integer) As Single
        Dim intTmp As Single

        intTmp = 10 ^ NroDec
        'Decimales = (Valor * intTmp) \ intTmp
        Decimales = Int((Valor * intTmp) / intTmp)
    End Function

    'retorna el orden de magnitud tipo (ej, de 2.5*10E3 retorna 3) del valor
    '1/8/03
    Public Function Order(ByVal Valor As Single) As Single
        Dim ValorTmp As Single

        If Valor <> 0 Then
            ValorTmp = Abs(Valor)
            ValorTmp = Math.Log10(ValorTmp)
            Order = Int(ValorTmp)
        Else
            Order = 1
        End If
    End Function

    Public Function IncrementEsc(ByVal Valor As Single, ByVal blnInc As Boolean) As Single
        Dim intTmp As Single
        'incrementa o decrementa el valor utilizando escala logaritmica
        '1/8/03

        intTmp = RoundLog((Valor), (blnInc))
        If blnInc Then
            intTmp = 2 * intTmp
        Else
            intTmp = 0.5 * intTmp
        End If
        IncrementEsc = RoundLog((intTmp), True)
    End Function



    'redondea un numero pero en la escala mas
    'cercana, tipo 1,2 5 10, 20, etc
    '1/8/03
    Public Function RoundMoreClose(ByVal Value As Single, ByVal NumberDigits As Integer) As Single
        Dim tmp As Single, Tmp2 As Single

        tmp = Abs(Value)
        Tmp2 = Order(tmp)
        tmp = Int(tmp / (10 ^ (Tmp2 * NumberDigits)))

        RoundMoreClose = tmp * (10 ^ (Tmp2 * NumberDigits))
        If Value < 0 Then
            RoundMoreClose = -RoundMoreClose
        End If
    End Function


    'redondea un numero pero en la escala exponencial mas
    'cercana, tipo 1,2 5 10, 20, etc
    '1/8/03
    Public Function RoundLog(ByVal Value As Single, ByVal blnRoundUp As Boolean) As Single
        Dim tmp As Single, Tmp2 As Single

        tmp = Value
        Tmp2 = Order(tmp)
        tmp = tmp / (10 ^ Tmp2)

        If blnRoundUp Then
            If tmp <= 0.5 Then
                RoundLog = 0.5 * (10 ^ Tmp2)
            ElseIf tmp <= 1 Then
                RoundLog = 1 * (10 ^ Tmp2)
            ElseIf tmp <= 2 Then
                RoundLog = 2 * (10 ^ Tmp2)
            ElseIf tmp <= 5 Then
                RoundLog = 5 * (10 ^ Tmp2)
            ElseIf tmp <= 10 Then
                RoundLog = 10 * (10 ^ Tmp2)
            Else
                Stop
                RoundLog = 1000
            End If
        Else
            If tmp < 0.5 Then
                RoundLog = 0.2 * (10 ^ Tmp2)
            ElseIf tmp < 1 Then
                RoundLog = 0.5 * (10 ^ Tmp2)
            ElseIf tmp < 2 Then
                RoundLog = 1 * (10 ^ Tmp2)
            ElseIf tmp < 5 Then
                RoundLog = 2 * (10 ^ Tmp2)
            ElseIf tmp < 10 Then
                RoundLog = 5 * (10 ^ Tmp2)
            Else
                Stop
                RoundLog = 1000
            End If
        End If
    End Function
#End Region

#Region "Numerical Recipes"

    'Linear equation solution by Gauss-Jordan elimination, equation (2.1.1) above. a[1..n][1..n]
    'is the input matrix. b[1..n][1..m] is input containing the m right-hand side vectors. On
    'output, a is replaced by its matrix inverse, and b is replaced by the corresponding set of solution
    'vectors.
    Public Sub gaussj(ByVal a(,) As Single, ByVal n As Integer, ByVal b()() As Single, ByVal m As Integer)
        Dim indxc(n) As Integer
        Dim indxr(n) As Integer
        Dim ipiv(n) As Integer
        Dim i, icol, irow, j, k, l, ll As Integer
        Dim big, dum, pivinv As Single

        icol = 0
        irow = 0
        For j = 0 To n - 1
            ipiv(j) = 0
        Next j
        For i = 0 To n - 1
            big = 0
            For j = 0 To n - 1
                If ipiv(j) <> 1 Then
                    For k = 0 To n - 1
                        If ipiv(k) = 0 Then
                            If Math.Abs(a(j, k)) >= big Then
                                big = Math.Abs(a(j, k))
                                irow = j
                                icol = k
                            End If
                        End If
                    Next k
                End If
            Next j
            ipiv(icol) += 1
            If irow <> icol Then
                For l = 0 To n - 1
                    'for (l=1;l<=n;l++)
                    SWAP(a(irow, l), a(icol, l))
                Next l 'for (l=1;l<=m;l++)
                For l = 0 To m - 1
                    SWAP(b(irow)(l), b(icol)(l))
                Next l
            End If
            indxr(i) = irow
            indxc(i) = icol
            If a(icol, icol) = 0 Then
                MessageBox.Show("gaussj: Singular Matrix") '
            End If
            pivinv = 1 / a(icol, icol)
            a(icol, icol) = 1
            For l = 0 To n - 1 '
                a(icol, l) *= pivinv
            Next l
            For l = 0 To m - 1
                b(icol)(l) *= pivinv
            Next l
            For ll = 0 To n - 1
                If ll <> icol Then
                    dum = a(ll, icol)
                    a(ll, icol) = 0
                    For l = 0 To n - 1
                        a(ll, l) -= a(icol, l) * dum
                    Next
                    For l = 0 To m - 1
                        b(ll)(l) -= b(icol)(l) * dum
                    Next l
                End If
            Next ll
        Next i
        For l = n - 1 To 0 Step -1
            If indxr(l) <> indxc(l) Then
                For k = 0 To n - 1
                    SWAP(a(k, indxr(l)), a(k, indxc(l)))
                Next k
            End If
        Next l
    End Sub 'gaussj

    Sub SWAP(ByVal a As Object, ByVal b As Object)
        Dim temp As Object = a
        a = b
        b = temp
    End Sub 'SWAP

#End Region

End Class


Public Class FFTcls

    Private Shared objMath As BVMathFunctions
    Public BlnError As TypeOfMathError


    Public Const cte2Pow63 = 9223372036854775808D '1,442695040888963407359924681

    Sub New()
        objMath = New BVMathFunctions

    End Sub

#Region "FFT Section"

#Region "FFT Utilities"

    Function NumberOfBitsNeeded(ByVal PowerOfTwo As Long) As Byte
        Dim I As Byte
        For I = 0 To 16
            If (PowerOfTwo And (2 ^ I)) <> 0 Then
                NumberOfBitsNeeded = I
                Exit Function
            End If
        Next
    End Function

    Function IsPowerOfTwo(ByVal X As Long) As Boolean
        If (X < 2) Then IsPowerOfTwo = False : Exit Function
        If (X And (X - 1)) = False Then IsPowerOfTwo = True
    End Function

    'devuelve el siguiente valor potencia de 2
    '-1 si hubo un error
    Public Function GetNextPowerOfTwo(ByVal Value As Single) As Long
        Try
            If Value > 0 Then
                Dim n As Long
                n = Math.Floor(objMath.Log2(Value))
                If n < 63 Then
                    GetNextPowerOfTwo = 2 ^ (n + 1)
                Else
                    GetNextPowerOfTwo = cte2Pow63
                    BlnError = TypeOfMathError.EOverFlow
                End If
            Else
                BlnError = TypeOfMathError.ECantBeNegative
            End If
        Catch ex As Exception
            GetNextPowerOfTwo = -1
        End Try
    End Function

    Public Function GetPreviousPowerOfTwo(ByVal Value As Single) As Long
        Try
            If Value > 0 Then
                If Value > 1 Then
                    Dim n As Long
                    n = Math.Ceiling(objMath.Log2(Value))

                    If n > 1 Then
                        If n < 64 Then
                            GetPreviousPowerOfTwo = 2 ^ (n - 1)
                        Else
                            GetPreviousPowerOfTwo = cte2Pow63
                            BlnError = TypeOfMathError.EOverFlow
                        End If
                    Else
                        GetPreviousPowerOfTwo = 2
                    End If
                Else
                    GetPreviousPowerOfTwo = 2
                    BlnError = TypeOfMathError.ECantBeNegative
                End If
            Else
                BlnError = TypeOfMathError.EThisIsLowerValue
            End If
        Catch ex As Exception
            GetPreviousPowerOfTwo = -1
        End Try
    End Function

    Function ReverseBits(ByVal Index As Long, ByVal NumBits As Byte) As Long
        Dim I As Byte, Rev As Long

        For I = 0 To NumBits - 1
            Rev = (Rev * 2) Or (Index And 1)
            Index = Index \ 2
        Next

        ReverseBits = Rev
    End Function

    Function FrequencyOfIndex(ByVal NumberOfSamples As Long, ByVal Index As Long) As Double
        'Based on IndexToFrequency().  This name makes more sense to me.

        If Index >= NumberOfSamples Then
            FrequencyOfIndex = 0.0#
            Exit Function
        ElseIf Index <= NumberOfSamples / 2 Then
            FrequencyOfIndex = CDbl(Index) / CDbl(NumberOfSamples)
            Exit Function
        Else
            FrequencyOfIndex = -CDbl(NumberOfSamples - Index) / CDbl(NumberOfSamples)
            Exit Function
        End If
    End Function


    Sub CalcFrequency(ByVal NumberOfSamples As Long, ByVal FrequencyIndex As Long, ByVal RealIn() As Double, ByVal ImagIn() As Double, ByVal RealOut As Double, ByVal ImagOut As Double)
        Dim K As Long
        Dim Cos1 As Double, Cos2 As Double, Cos3 As Double, Theta As Double, Beta As Double
        Dim Sin1 As Double, Sin2 As Double, Sin3 As Double

        Theta = 2 * PI * FrequencyIndex / CDbl(NumberOfSamples)
        Sin1 = Sin(-2 * Theta)
        Sin2 = Sin(-Theta)
        Cos1 = Cos(-2 * Theta)
        Cos2 = Cos(-Theta)
        Beta = 2 * Cos2

        For K = 0 To NumberOfSamples - 2
            'Update trig values
            Sin3 = Beta * Sin2 - Sin1
            Sin1 = Sin2
            Sin2 = Sin3

            Cos3 = Beta * Cos2 - Cos1
            Cos1 = Cos2
            Cos2 = Cos3

            RealOut = RealOut + RealIn(K) * Cos3 - ImagIn(K) * Sin3
            ImagOut = ImagOut + ImagIn(K) * Cos3 + RealIn(K) * Sin3
        Next
    End Sub

#End Region

#Region "FFT Methods"

    Public Sub FourierTransform(ByVal NumSamples As Long, ByVal RealIn() As Double, ByVal ImageIn() As Double, ByVal RealOut() As Double, ByVal ImagOut() As Double, Optional ByVal InverseTransform As Boolean = False)
        Dim AngleNumerator As Double
        Dim NumBits As Byte, I As Long, j As Long, K As Long, n As Long, BlockSize As Long, BlockEnd As Long
        Dim DeltaAngle As Double, DeltaAr As Double
        Dim Alpha As Double, Beta As Double
        Dim TR As Double, TI As Double, AR As Double, AI As Double

        If InverseTransform Then
            AngleNumerator = -2.0# * PI
        Else
            AngleNumerator = 2.0# * PI
        End If

        If (IsPowerOfTwo(NumSamples) = False) Or (NumSamples < 2) Then
            Call MsgBox("Error in procedure Fourier:" + vbCrLf + " NumSamples is " + CStr(NumSamples) + ", which is not a positive integer power of two.", , "Error!")
            Exit Sub
        End If

        NumBits = NumberOfBitsNeeded(NumSamples)
        For I = 0 To (NumSamples - 1)
            j = ReverseBits(I, NumBits)
            RealOut(j) = RealIn(I)
            ImagOut(j) = ImageIn(I)
        Next

        BlockEnd = 1
        BlockSize = 2

        Do While BlockSize <= NumSamples
            DeltaAngle = AngleNumerator / BlockSize
            Alpha = Sin(0.5 * DeltaAngle)
            Alpha = 2.0# * Alpha * Alpha
            Beta = Sin(DeltaAngle)

            I = 0
            Do While I < NumSamples
                AR = 1.0#
                AI = 0.0#

                j = I
                For n = 0 To BlockEnd - 1
                    K = j + BlockEnd
                    TR = AR * RealOut(K) - AI * ImagOut(K)
                    TI = AI * RealOut(K) + AR * ImagOut(K)
                    RealOut(K) = RealOut(j) - TR
                    ImagOut(K) = ImagOut(j) - TI
                    RealOut(j) = RealOut(j) + TR
                    ImagOut(j) = ImagOut(j) + TI
                    DeltaAr = Alpha * AR + Beta * AI
                    AI = AI - (Alpha * AI - Beta * AR)
                    AR = AR - DeltaAr
                    j = j + 1
                Next

                I = I + BlockSize
            Loop

            BlockEnd = BlockSize
            BlockSize = BlockSize * 2
        Loop

        If InverseTransform Then
            'Normalize the resulting time samples...
            For I = 0 To NumSamples - 1
                RealOut(I) = RealOut(I) / NumSamples
                ImagOut(I) = ImagOut(I) / NumSamples
            Next
        End If
    End Sub

    'transformad de Fourier, recibe dos arrays complexBinomical, uno sale modificado
    Public Sub FourierTransform(ByVal NumSamples As Long, ByVal fIn() As ComplexBinomical, ByVal fOut() As ComplexBinomical, Optional ByVal InverseTransform As Boolean = False)
        Dim AngleNumerator As Double
        Dim NumBits As Byte
        Dim I, j, K, n, BlockSize, BlockEnd As Long
        Dim DeltaAngle, DeltaAr As Double
        Dim Alpha, Beta As Double
        Dim TR, TI, AR, AI As Double

        If InverseTransform Then
            AngleNumerator = -2.0# * PI
        Else
            AngleNumerator = 2.0# * PI
        End If

        If (IsPowerOfTwo(NumSamples) = False) Or (NumSamples < 2) Then
            Call MsgBox("Error in procedure Fourier:" + vbCrLf + " NumSamples is " + CStr(NumSamples) + ", which is not a positive integer power of two.", , "Error!")
            Exit Sub
        End If

        NumBits = NumberOfBitsNeeded(NumSamples)
        For I = 0 To (NumSamples - 1)
            j = ReverseBits(I, NumBits)
            fOut(j).Real = fIn(I).Real
            fOut(j).Imaginario = fIn(I).Imaginario
        Next

        BlockEnd = 1
        BlockSize = 2

        Do While BlockSize <= NumSamples
            DeltaAngle = AngleNumerator / BlockSize
            Alpha = Sin(0.5 * DeltaAngle)
            Alpha = 2.0# * Alpha * Alpha
            Beta = Sin(DeltaAngle)

            I = 0
            Do While I < NumSamples
                AR = 1.0#
                AI = 0.0#

                j = I
                For n = 0 To BlockEnd - 1
                    K = j + BlockEnd
                    TR = AR * fOut(K).Real - AI * fOut(K).Imaginario
                    TI = AI * fOut(K).Real + AR * fOut(K).Imaginario
                    fOut(K).Real = fOut(j).Real - TR
                    fOut(K).Imaginario = fOut(j).Imaginario - TI
                    fOut(j).Real = fOut(j).Real + TR
                    fOut(j).Imaginario = fOut(j).Imaginario + TI
                    DeltaAr = Alpha * AR + Beta * AI
                    AI = AI - (Alpha * AI - Beta * AR)
                    AR = AR - DeltaAr
                    j = j + 1
                Next

                I = I + BlockSize
            Loop

            BlockEnd = BlockSize
            BlockSize = BlockSize * 2
        Loop

        If InverseTransform Then
            'Normalize the resulting time samples...
            For I = 0 To NumSamples - 1
                fOut(I).Real = fOut(I).Real / NumSamples
                fOut(I).Imaginario = fOut(I).Imaginario / NumSamples
            Next
        End If
    End Sub

#End Region

#End Region

End Class

By viewing downloads associated with this article you agree to the Terms of Service and the article's licence.

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


Written By
Engineer Universidad Tecnológica Nacional
Argentina Argentina
This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.

Comments and Discussions