codigoTest.py

import scipy.special 
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
from mpmath import *

mpl.rcParams["text.usetex"] = True

def riccatiBessel_jn(n, z):
    return z * scipy.special.spherical_jn(n, z, derivative=False)

def riccatiBessel_jn_derivative(n, z):
    return (1+n) * scipy.special.spherical_jn(n, z, derivative=False) - z * scipy.special.spherical_jn((1+n), z, derivative=False)

def spherical_hankel_2k(n, z):
    jn = scipy.special.spherical_jn(n, z, derivative=False)
    yn = scipy.special.spherical_yn(n, z, derivative=False)
    
    return jn - (yn * 1j)


def riccatiBessel_hn(n, z):
    return z * spherical_hankel_2k(n, z)
    

def riccatiBessel_hn_derivative(n, z):
    return (1+n) * spherical_hankel_2k(n, z) - z * spherical_hankel_2k((1+n), z)

def varC_n (u_r, M, ka, n):
    chi_ka = riccatiBessel_hn(n, ka)
    chiD_ka = riccatiBessel_hn_derivative(n, ka)
    psi_ka = riccatiBessel_jn(n, ka)
    psi_Mka = riccatiBessel_jn(n, (M*ka))
    psiD_ka = riccatiBessel_jn_derivative(n, ka)
    psiD_Mka = riccatiBessel_jn_derivative(n, (M*ka))
    
    num = M * u_r * ((chi_ka * psiD_ka) - (chiD_ka * psi_ka))
    den = (u_r * chi_ka * psiD_Mka) - (M * chiD_ka * psi_Mka)
    
    return num/den

def varD_n (u_r, M, ka, n):
    chi_ka = riccatiBessel_hn(n, ka)
    chiD_ka = riccatiBessel_hn_derivative(n, ka)
    psi_ka = riccatiBessel_jn(n, ka)
    psi_Mka = riccatiBessel_jn(n, (M*ka))
    psiD_ka = riccatiBessel_jn_derivative(n, ka)
    psiD_Mka = riccatiBessel_jn_derivative(n, (M*ka))
    
    num = M * M * ((chi_ka * psiD_ka) - (chiD_ka * psi_ka))
    den = (M * chi_ka * psiD_Mka) - (u_r * chiD_ka * psi_Mka)
    
    return num/den

ceilingX = lambda x : math.ceil(x + (4.05 * (x ** (1/3))) + 2)

def Rn (M, x, n):
    Mx = M * x
    psi_n1Mx = riccatiBessel_jn((n+1), Mx)
    psi_Mx = riccatiBessel_jn(n, Mx)
    
    num = M * psi_n1Mx * np.conj(psi_Mx)
    den = M ** 2
    

    return num.imag/den.imag

def Sn (M, x, n):
    M2 = M ** 2
    Mx = M * x
    conjM = np.conj(M)
    rn = Rn(M, x, n)
    rn1 = Rn(M, x, (n+1))
    mod2Psi_Mx = abs(riccatiBessel_jn(n, Mx)) ** 2
    mod2Psin1_Mx = abs(riccatiBessel_jn((n+1), Mx)) ** 2
    
    fristTerm = -(1j / (2 * M2.imag))
    secondTerm = x * ((M * mod2Psi_Mx) + (conjM * mod2Psin1_Mx))
    thirdTerm = (M + ((2 * (n+1)) * (M2.real / M))) * rn
    fourthTerm = ((2*n) + 1) * conjM * rn1
    
    return fristTerm * (secondTerm - thirdTerm + fourthTerm)

def J1(x, M, epislon2l, ur):
    firstTerm = (3 * epislon2l) / ((abs(M) ** 2) * (x ** 3))
    summationResult = 0
    nMax = ceilingX(x)
    etaR = 1 / M
    mod2EtaR = abs(etaR) ** 2
    
    for i in range(1, (nMax+1)):
        cn = varC_n (ur, M, x, i) 
        conjCn = np.conj(cn)
        cn1 = varC_n (ur, M, x, (i+1))
        conjCn1 = np.conj(cn1)
        rn = Rn (M, x, i)
        rn1 = Rn (M, x, (i+1))
        dn = varD_n (ur, M, x, i)
        conjDn = np.conj(dn)
        dn1 = varD_n (ur, M, x, (i+1))
        conjDn1 = np.conj(dn1)
        sn = Sn (M, x, i)
        
        fristTermFor = ((i*(i+2)) / M) * ((cn1*conjCn*rn1) + (mod2EtaR*dn1*conjDn*rn))
        secondTermFor = ((i*(i+2)) / (i+1)) * ((cn*conjCn1) + (mod2EtaR*dn1*conjDn))
        thirdTermFor =  np.conj(etaR) * (((2*i)+1) / (i*(i+1))) * (cn*conjDn)
        
        result = fristTermFor - ((secondTermFor + thirdTermFor) * sn)
        
        summationResult += result
        
    return firstTerm * summationResult.imag

def imagem1():
    m = 1.57 - 0.038j
    x = np.linspace(0.01, 20, 200)
    ur = 1
    eps_r2l = - (m**2).imag
    results = []

    for i in x:
        results.append(J1(i, m, eps_r2l, ur))

    plt.figure(figsize=[7,5])
    plt.plot(x, results, 'g', label= "m = 1.57 - 0.038j")
    plt.xlabel('Size Parameter x')
    plt.ylabel('Asymmetry Factor J1')
    plt.ylim(-0.10, 0.15) 
    plt.xlim(0, 20)
    plt.grid()
    plt.legend(loc='best')
    plt.show()
    
def imagem2():
    m = 1.57 - 0.38j
    x = np.linspace(0.01, 20, 200)
    ur = 1
    eps_r2l = - (m**2).imag
    results = []

    for i in x:
        results.append(J1(i, m, eps_r2l, ur))

    plt.figure(figsize=[7,5])
    plt.plot(x, results, 'g', label= "m = 1.57 - 0.38j")   
    plt.xlabel('Size Parameter x')
    plt.ylabel('Asymmetry Factor J1')
    plt.ylim(-0.50, 0.01)
    plt.xlim(0, 20)
    plt.grid()
    plt.legend(loc='best')
    plt.show()

    
imagem1()
imagem2()