<jupyter_start><jupyter_text>Development of spectral solution for Teukolksy's radial equation1. Given a, w, n, m and s, solve for the separation constant A(n) using the minimal solution constraint for the radial equaiton 2. Given A(n), evaluate the series solution for the zero spin radial equation in $\xi$ coordinate space3. Validate calculation of zero spin radial functions on $\xi$ by plotting convergence of the series, and various cases in m and n<jupyter_code>%load_ext autoreload
%autoreload 2
%matplotlib inline
# The Important Things
from matplotlib.pyplot import *
from numpy import *
from positive import *
error(' This notebook  is for development purposees only. The related code and results are known to have significant limitations. The question remains regarding whether there exists a polynomial basis capable of accurately representing the radial functions. ')
# Calculate reference QNM frequency and separation constant
a = 0.9
s = -2
l,m = 2,2
M = 1
n = 0
#
cw_ref,sc_ref = leaver(a,l,m,n,s=s)
print(cw_ref,sc_ref)
# Try calling rnspec_helper
Q_r,vals_r,vecs_r,nrange     = rnspec_legendre_helper(  a, cw_ref, m, n, s=-2, span=20 )
from scipy.special import gamma, factorial
myfactorial = lambda z: gamma(z+1)
z = -3
print(factorial(z))
print(myfactorial(z))
# Try calling rnspec_helper
Q_r2,vals_r2,vecs_r2,nrange_r2     = rnspec_jacobi_helper(  a, cw_ref, m, n, s=-2, span=20 )

Q_a,vals_a_raw,vecs_a,lrange = slmcg_helper( a*cw_ref, s, l, m, span=6 )
vals_a = -vals_a_raw

Q_a2,vals_a2,vecs_a2,lrange_a2 = slspec_jacobi_helper( a*cw_ref, s, l, m, span=6 )


print()

print( sc_ref )
print( vals_a[lrange.index(l)] )
print( vals_a2[lrange_a2.index(l)] )
print( vals_r[nrange.index(n)] )
print( vals_r2[nrange_r2.index(n)] )
print()

print( abs(sc_ref-vals_r[nrange.index(n)]) )
print( abs(sc_ref-vals_a[lrange.index(l)]) )
print( abs(sc_ref-vals_a2[lrange_a2.index(l)]) )
print()

vals_r2
plot( vals_r.real, vals_r.imag, ls='none', marker='s', mfc='none' )
plot( vals_r2.real, vals_r2.imag, ls='none', marker='.', mfc='none', c='k' )

figure( figsize=2*figaspect(0.618) )

clr = rgb(len(vals_a), jet=True, reverse=not True)
for nn,v in enumerate(vals_a):
    plot( v.real, v.imag, marker='.', ls='none', label='from angular eigen-problem' if nn==0 else None, ms=10, mew=2, c=clr[nn] )

clr = rgb(len(vals_r), jet=True, reverse=True)
for nn,v in enumerate(vals_r):
    plot( v.real, v.imag, marker='o', ls='none', label='from radial eigen-problem' if nn==0 else None, ms=10, mfc='none', mew=2, c=clr[nn] )



An = vals_r[n]
Aa = vals_a[l-abs(m)]

plot( An.real, An.imag, marker='o', color='k', ms=20, mfc='none' )
plot( Aa.real, Aa.imag, marker='x', color='k', ms=14, mfc='none' )    

    
# #
# K = 10
# xlim(-K,K)
# ylim(-1,K*K)

# xlim(-250,250)
# ylim(-40,300)

#
x0,x1 = xlim()
y0,y1 = ylim()
plot( (x0,x1), (0,0), lw=1, color='k', zorder=-100, alpha=0.5 )
plot( (0,0), (y0,y1), lw=1, color='k', zorder=-100, alpha=0.5 )

# #
# xlim(x0,x1)
# ylim(y0,y1)

#
xlabel('$\Re\, A_{\ell m n}$')
ylabel('$\Im\, A_{\ell m n}$')

# yscale('log')

xlim(-250,250)
ylim(-40,300)

#
title('$(\ell,m,n)=(%i,%i,%i)$, a=%f'%(l,m,n,a))

#
legend(frameon=True)<jupyter_output><empty_output><jupyter_text>Algorithm for solving for QNMs1. define $\ell,m,n$ and $a$2. define some trial w (this is an initial guess, but the method ahead may not be sensitive to it)3. call rnspec_helper and slmcg_helper4. use n to choose the n'th eigenvalue from rnspec_helper5. use l, m and s to choose the correct eigenvalue from slmcg_helper6. we define the work function to be the distance between these two uniquely defined eigenvalues --<jupyter_code>def eigqnm(a,l,m,n,wstart=None,span=6,leaver_refine=False,verbose=False,method=None,full_output=False):

    #
    from scipy.optimize import minimize
    from numpy import real
    
    #
    if method is None:
        method = 'Powell'
    
    #
    if wstart is None:
        wstart = 0.2+0.1j
    #
    wguess = [abs(wstart),angle(wstart)] # [wstart.real,wstart.imag]
    
    #
    alert('Welcome to %s'%magenta('eigqnm'),header=True,verbose=verbose)
    
    #
    def action( X ):
        
        #
        trial_cw = X[0] * exp(1j*X[1]) # X[0]+1j*X[1]
        # _,vals_r,_,    nrange = rnspec_legendre_helper(  a, trial_cw, m, n, s=-2, span=span )
        _,vals_r,_,    nrange = rnspec_jacobi_helper(  a, trial_cw, m, n, s=-2, span=span )
        _,vals_a_raw,_,lrange = slmcg_helper( a*trial_cw, s, l, m, span=6 )
        vals_a = -vals_a_raw
        
        #
        sca = vals_a[lrange.index(l)]
        scr = vals_r[nrange.index(n)]
        res = real(abs(sca-scr)/( abs(sca)+abs(scr) ))
        
        #alert(res,verbose=verbose)
        
        #
        return res*res
    
    #
    alert('method = %s'%blue(method),verbose=verbose)
    
    #
    options = { 'ftol' : 1e-15}
    foo = minimize( action, wguess, method=method, options=options, tol=1e-15 )
    
    #
    if 'nit' in foo.keys():
        alert('nit : %i'%foo.nit, verbose=verbose)
    if 'nfev' in foo.keys():
        alert('nfev: %i'%foo.nfev,verbose=verbose)
    if 'fun' in foo.keys():
        alert('fun : %g'%foo.fun,verbose=verbose)
    if 'message' in foo.keys():
        alert('note: %s'%(foo.message),verbose=verbose)
    
    #
    cw = foo.x[0]*exp(1j*foo.x[1]) # foo.x[0] + 1j*foo.x[1]
    
    #
    _,vals_a_raw,_,lrange = slmcg_helper( a*cw, s, l, m, span=6 )
    sc = -vals_a_raw[lrange.index(l)]
    er = foo.fun
    
    # refine with continued fraction method
    if leaver_refine:
        #
        guess = array([cw.real,cw.imag,sc.real,sc.imag],dtype=float)
        refined_cw,refined_sc,refined_lvrwrk,retry = lvrsolve(a,l,m,guess,tol=1e-8,s=s)
        #
        cw = refined_cw
        sc = refined_sc
        er = refined_lvrwrk
    
    #
    if full_output:
        foo.sc = sc
        foo.cw = cw
        foo.er = er
        return foo
    else:
        return sc,cw,er
    
a = 0.0
l = 2
m = 2
n = 0

if n<4:
    cw_ref,sc_ref = leaver(a,l,m,n,s=s)
else:
    cw_ref,sc_ref = 0,0
    
#
foo = eigqnm(a,l,m,n,span=13,verbose=True,full_output=True) 

#
test_sc = foo.sc
test_cw = foo.cw
err = foo.er
#
foo = eigqnm(a,l,m,n,span=6,wstart=foo.cw,verbose=True,full_output=True) 

#
test_sc = foo.sc
test_cw = foo.cw
err = foo.er
print()
print( test_sc )
print( sc_ref )
print()
print( test_cw )
print( cw_ref )
print()
print(err)
print( abs(test_cw-cw_ref) )
print()
print( leaver_workfunction( a,l,m, [cw_ref.real,cw_ref.imag,sc_ref.real,sc_ref.imag], s=s ) )
print( leaver_workfunction( a,l,m, [test_cw.real,test_cw.imag,test_sc.real,test_sc.imag], s=s ) )
print()
guess = array([test_cw.real,test_cw.imag,test_sc.real,test_sc.imag],dtype=float)
refined_cw,refined_sc,refined_lvrwrk,retry = lvrsolve(a,l,m,guess,tol=1e-8,s=s)
print(refined_cw,refined_sc,refined_lvrwrk)
print(cw_ref,sc_ref)
# nspace = list(range(6))
# cw = []
# sc = []
# er = []
# for nn in nspace:
#     print(nn)
#     test_sc,test_cw,err = eigqnm(a,l,m,nn,wstart = 0.5-0.26j if nn==0 else test_cw, span=6 ) 
#     cw.append(test_cw)
#     sc.append(test_sc)
#     er.append(err)
# cw = array(cw)
# print(array(er))
# print(array(sc))
# print(array(cw))
# sc = array(sc)
# er = array(er)
# figure()
# plot( sc.real, sc.imag, marker='o', ls='none' )
# figure()
# plot( cw.real, -cw.imag, marker='o', ls='none' )
# figure()
# plot( log10(er), marker='o', ls='none' )
from scipy.optimize import minimize,root
?root<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16209342792129566, 'fraction_numerical': 0.02221924775799759, 'mean_word_length': 3.717171717171717, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '214536', 'license': 'MIT', 'n_tokens_mistral': 3255, 'n_tokens_neox': 3162, 'n_words': 670, 'path': 'issues/open/dev_rnspec_helper.ipynb', 'repo_name': 'llondon6/koalas', 'size': '7471'}