import numpy as np

"""
Explanation: FD_1D_stability 1-D FD stability calculation
GNU General Public License v3.0
Author: Florian Wittkamp
Calculate stability limit of FD-simulations
Stability limit is calculated in terms of the CFL-number, which is
defined as: CFL=v_(max)DT/DX
You get the maximum DT by DT=CFLDX/v_(max)
Theory:
Fei, X., & Xiaohong, T. (2006).
Stability and numerical dispersion analysis of a fourth-order accurate FDTD method. Antennas and Propagation,
IEEE Transactions on, 54(9), 2525-2530.
Initialisation
End of explanation
"""


spatial_order=4

"""
Explanation: Input Parameter
End of explanation
"""


def coeff(order):
    ## Check some conditions
    if int(order)%2!=0:
        print("Error: coeff")
        print("Order has to be an integer multiple of 2!")
        return
    if order==2:
        print("Error: coeff")
        print("Order has to be at least 4!")
        return
    ## Calculation
    c=np.transpose(np.hstack((1, np.zeros(int(order/2)-1))))
    M=np.zeros((int(order/2),int(order/2)))
    # Condition 1: \sum^{N/2}_{k=1} b_k(2k-1)=1
    for n in range(1,int(order/2+1)):
        M[0,n-1]=(2*n-1)
    # Condition 2:  \sum^{N/2}_{k=1} b_k(2k-1)^(2j-1)=0; j=2,3...N/2
    for j in range(2,int(order/2+1)):
        for n in range(1,int(order/2+1)):
            M[j-1,n-1]=(2*n-1)**(2*j-1)
    coeff=np.transpose(np.dot(np.linalg.inv(M),c))
    return(coeff)

"""
Explanation: Calculate Taylor coefficient
Calculate the Taylor coefficient in an arbitrary order for the
second-order derivative.
Usage:
Lets say you want to calculate the 4th order accurate second-order
FD-stencil. Then you have to set order=4 and the result would be used
as follow:
p_x = 1/DH * ( coeff(1) * (p(x+1)-p(x)) + coeff(2) * (p(x+2)-p(x-1)) )
where p_x is the derivative.
End of explanation
"""


sum_fd_stencil=np.sum(np.abs(coeff(spatial_order)))
theta=np.arange(0,2*np.pi,0.01)
print("You choose a spatial order of ",spatial_order)

"""
Explanation: Calculate spatial derivation impact
Assumption: Spatial sampling is at the Nyquist condition
End of explanation
"""


f=-1j*(((np.exp(-1j*theta))**(2.5)-(np.exp(-1j*theta))**1.5))
f2=(2*sum_fd_stencil)
f_ges=f/f2
print("2. Temporal Order has the stability limit CFL=",np.max(f_ges.real))


"""
Explanation: Temporal 2 order (Leapfrog)
End of explanation
"""


f=-1j*(((np.exp(-1j*theta))**2.5-(np.exp(-1j*theta))**1.5)*1);
f2=(2*sum_fd_stencil*(25.0/24.0*(np.exp(-1j*theta))**2
    -1./12.*(np.exp(-1j*theta))**1+1./24.))
f_ges=f/f2
print("3. Temporal Order (ABS-method) has the stability limit CFL=",np.max(f_ges.real))

"""
Explanation: Temporal 3 order (Adams-Bashforth method)
End of explanation
"""


f=-1j*(((np.exp(-1j*theta))**3.5-(np.exp(-1j*theta))**2.5))
f2=(2*sum_fd_stencil*(13./12.*(np.exp(-1j*theta))**3
    -5./24.*(np.exp(-1j*theta))**2+1./6.*(np.exp(-1j*theta))**1-1./24.))
f_ges=f/f2
print("4. Temporal Order (ABS-method) has the stability limit CFL=",np.max(f_ges.real))

"""
Explanation: Temporal 4 order (Adams-Bashforth method)
End of explanation
"""


