Attribute VB_Name = "Sine_Stability"
Function SinStabX(ma As Double, q As Double, nmax As Integer) 'call for the function in cell as =SinStabX(a,q,Calc Step)

'ma is Mathieu parameter a
'q is the Mathieu parameter q
'nmax is Calc Step from the spreadsheet

 r = nmax * 1#
 tau = 1 / r 'step increment
  tau2 = 0   'step in the waveform
  trace = 0
If ma = 0 Then
ma = 1E-23          'apparently VBA does not like it to be 0 even if the calc doesn't use it
End If

  Pi = WorksheetFunction.Pi()

     For n = 1 To nmax
         

              If n = 1 Then 'create the identiy matrix and use once
                  
                    a = 1#
                    b = 0#
                    c = 0#
                    d = 1#
                    
                  End If

    
               
               f1 = ma + 2 * q * Sin(tau2 * Pi * 2) 'creation of the approximated sine wave
                                                    'f1 is the waveform potential
              
             
                 If f1 > 0 Then
            
                    'rather than an array or matrix we find the elements
                    
                    e = Cos(Pi * tau * Sqr(f1))
                    f = (1 / Sqr(f1)) * Sin(Pi * tau * Sqr(f1))
                    g = -Sqr(f1) * Sin(Pi * tau * Sqr(f1))
                    h = Cos(Pi * tau * Sqr(f1))
                
                End If

                If f1 < 0 Then
                
                   

                    e = WorksheetFunction.Cosh(Pi * tau * Sqr(-f1))
                    f = (1 / Sqr(-f1)) * WorksheetFunction.Sinh(Pi * tau * Sqr(-f1))
                    g = Sqr(-f1) * WorksheetFunction.Sinh(Pi * tau * Sqr(-f1))
                    h = WorksheetFunction.Cosh(Pi * tau * Sqr(-f1))

                End If
               
            a1 = (a * e + b * g) 'these is an analytical form of any pair of two by two matix multiplication
            b1 = (a * f + b * h)
            c1 = (c * e + d * g)
            d1 = (c * f + d * h)

            trace = a1 + d1     'for stability we need the trace only

            a = a1              'redefine the n = 1 identity as left matrix
            b = b1              'each time through the loop this representation is the left matrix in the next step
            c = c1
            d = d1
            
  tau2 = tau + tau2             'increase the waveform through the period
         
Next

 SinStabX = trace               'send the last trace to the spreadsheet cell
 End Function
 
 Function SinStabY(ma As Double, q As Double, nmax As Integer)


 r = nmax * 1#
 tau = 1 / r
  trace = 0
If ma = 0 Then
ma = 1E-23
End If

  Pi = WorksheetFunction.Pi()

     For n = 1 To nmax
         

              If n = 1 Then
                  
                    a = 1#
                    b = 0#
                    c = 0#
                    d = 1#
                    
                  End If

    
               
               f1 = -ma - 2 * q * Sin(tau2 * Pi * 2) 'Y-axis gets the opposite waveform potential, rest is the same

              
             
                 If f1 > 0 Then
            
            

                    e = Cos(Pi * tau * Sqr(f1))
                    f = (1 / Sqr(f1)) * Sin(Pi * tau * Sqr(f1))
                    g = -Sqr(f1) * Sin(Pi * tau * Sqr(f1))
                    h = Cos(Pi * tau * Sqr(f1))
                
                End If

                If f1 < 0 Then
                
                   

                    e = WorksheetFunction.Cosh(Pi * tau * Sqr(-f1))
                    f = (1 / Sqr(-f1)) * WorksheetFunction.Sinh(Pi * tau * Sqr(-f1))
                    g = Sqr(-f1) * WorksheetFunction.Sinh(Pi * tau * Sqr(-f1))
                    h = WorksheetFunction.Cosh(Pi * tau * Sqr(-f1))

                End If
               
            a1 = (a * e + b * g)
            b1 = (a * f + b * h)
            c1 = (c * e + d * g)
            d1 = (c * f + d * h)

            trace = a1 + d1

            a = a1
            b = b1
            c = c1
            d = d1
            
  tau2 = tau + tau2
         
Next

 SinStabY = trace
 End Function