e2plus.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Jan 14 23:02:31 2022

@author: noumanbutt
"""

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jan 12 11:12:50 2022

@author: noumanbutt
"""

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Oct 11 14:12:26 2021
@author: noumanbutt
"""

import numpy as np
import matplotlib.pyplot as plt
import pickle
from scipy.special import spherical_jn
from scipy import integrate
from itertools import product,permutations

from scipy import arange, pi, sqrt, zeros
from scipy.special import jv, jvp
from scipy.optimize import brentq
from sys import argv
import time
#import matplotlib.pyplot as plt
from scipy.sparse import  isspmatrix

from angular_integrals import angular_integral_values_2p



import sympy
from sympy import S
from sympy.physics.wigner import wigner_3j
from sympy.physics.wigner import clebsch_gordan
from sympy.physics.wigner import gaunt

# print out sample integrals

### recursive method: computes zeros ranges of Jn(r,n) from zeros of Jn(r,n-1)
### (also for zeros of (rJn(r,n))')
### pros : you are certain to find the right zeros values;
### cons : all zeros of the n-1 previous Jn have to be computed;
### note : Jn(r,0) = sin(r)/r




def Jn(r,n):
  return (sqrt(pi/(2.*r))*jv(n+0.5,r))
def Jn_zeros(n,nt):
  zerosj = zeros((n+1, nt), dtype=float)
  zerosj[0] = arange(1,nt+1)*pi
  points = arange(1,nt+n+1)*pi
  racines = zeros(nt+n, dtype=float)
  for i in range(1,n+1):
    for j in range(nt+n-i):
      foo = brentq(Jn, points[j], points[j+1], (i,))
      racines[j] = foo
    points = racines
    zerosj[i][:nt] = racines[:nt]
  return (zerosj)

def rJnp(r,n):
  return (0.5*sqrt(pi/(2*r))*jv(n+0.5,r) + sqrt(pi*r/2)*jvp(n+0.5,r))
def rJnp_zeros(n,nt):
  zerosj = zeros((n+1, nt), dtype=float)
  zerosj[0] = (2.*arange(1,nt+1)-1)*pi/2
  points = (2.*arange(1,nt+n+1)-1)*pi/2
  racines = zeros(nt+n, dtype=float)
  for i in range(1,n+1):
    for j in range(nt+n-i):
      foo = brentq(rJnp, points[j], points[j+1], (i,))
      racines[j] = foo
    points = racines
    zerosj[i][:nt] = racines[:nt]
  return (zerosj)


def integrand(knp,kn,lp,l,i,j):
    
     
     f = lambda x:(x**(2 + 2*i + 2*j))*Jn(kn*x,l)*Jn(knp*x,lp)
     return f
 
def normed_integrand(knp,kn,lp,l,i,j,norm_lp,norm_l):
    
     
     f = lambda x:(x**(2 + 2*i + 2*j))*(Jn(kn*x,l)/norm_l)*(Jn(knp*x,lp)/norm_lp)
     return f    
    
def radial_integral(f):
    value, err = integrate.quadrature(f,0.0,np.pi)
    return value,err









'''
n = 1 #int(argv[1])  # n'th spherical bessel function
nt = 10 #int(argv[2]) # number of zeros to be computed

dr = 0.01
eps = dr/1000

jnz = Jn_zeros(n,nt)[n]
print(jnz/pi)
r1 = arange(eps,jnz[len(jnz)-1],dr)
jnzp = rJnp_zeros(n,nt)[n]
print(jnzp/pi)
r2 = arange(eps,jnzp[len(jnzp)-1],dr)


plt.plot(r1,Jn(r1,n),'b', r2,rJnp(r2,n),'r')
plt.title((str(nt)+' first zeros'))
plt.legend((r'$j_{'+str(n)+'}(r)$', r'$(rj_{'+str(n)+'}(r))\'$'))
plt.plot(jnz,zeros(len(jnz)),'bo', jnzp,zeros(len(jnzp)),'rd')
plt.axhline(y=0.0, color='black', linestyle='-')
#plt.gca().xaxis.set_minor_locator(MultipleLocator(1))
#plt.xaxis.set_minor_formatter(FormatStrFormatter('%d'))
plt.show()

'''


#k0_val=np.zeros([lmax+1,nmax])
#k1_val=np.zeros([lmax+1,nmax])


#for l in range(0,lmax+1):
    
#    k0_val[l] = Jn_zeros(l,nmax)[l]
#    k1_val[l] = rJnp_zeros(l,nmax)[l]
    
#print(k1_val/np.pi)    




    


    


def radial_integral_values_1p(lmax,nmax,e,norms):
    time0 = time.time()
    vals_i1_j0 = np.zeros((nmax,nmax,lmax+1))  # r_1^2 r_2^2
    vals_i0_j1 = np.zeros((nmax,nmax,lmax-1))  
    vals_i2_j0 = np.zeros((nmax,nmax,lmax+1))  # r_1^4
    vals_i3_j0 =  np.zeros((nmax,nmax,lmax+1)) # r_1^6
    
    
    for m,mp in product(range(0,nmax),repeat=2):
        #if(abs(m-mp) <= 1):
          for l in range(0,lmax+1):
            
            vals_i1_j0[m,mp,l], err = radial_integral(normed_integrand(e[l,mp],e[l,m],l,l,1,0,norms[l,mp],norms[l,m]))
        #else:
        #  continue  
    for m,mp in product(range(0,nmax),repeat=2):
        #if(abs(m-mp) <= 1):
          for l in range(0,lmax+1):
            
            vals_i2_j0[m,mp,l], err = radial_integral(normed_integrand(e[l,mp],e[l,m],l,l,2,0,norms[l,mp],norms[l,m]))
        #else:
        #  continue
    
    
    
    for m,mp in product(range(0,nmax),repeat=2):
        #if(abs(m-mp) <= 1):
          for l in range(0,lmax+1):
            
            vals_i3_j0[m,mp,l], err = radial_integral(normed_integrand(e[l,mp],e[l,m],l,l,3,0,norms[l,mp],norms[l,m]))
        #else:
        #  continue
    
    
    
    
    
    
    
    
    for m,mp in product(range(0,nmax),repeat=2):
        
      #if(abs(m-mp) <= 1):  
        for l in range(0,lmax-1):
        
            vals_i0_j1[m,mp,l], err = radial_integral(normed_integrand(e[l+2,mp],e[l,m],l+2,l,0,1,norms[l+2,mp],norms[l,m]))
            
      #else:
       # continue
            
    time1 = time.time()
    print("Time for computing radial integrals: ",time1-time0)                
    return vals_i1_j0,vals_i0_j1,vals_i2_j0,vals_i3_j0   
        
 

#kinetic term v1

def kinetic_term(l1,l2,l3,n1,n2,n3,e0,e1,g):
    
    if( g >= 1.2):
     return 0.5 *(e1[l1,n1]**2 + e1[l2,n2]**2 + e0[l3,n3]**2)
    else: 
     return (e1[l1,n1] + e1[l2,n2] + e0[l3,n3])


#radial term v1


def radial_terms1p(l1,l2,l3,n1,n2,n3,n1p,n2p,n3p,vals_i1_j0,vals_i2_j0,vals_i3_j0,g):
    
    alpha_1 = 2.1810429e-2 ## type B
    val = 0.0
    if(n2 == n2p  and n3 == n3p):
      
      if( g < 1.2):  
        
        val += -(1./(1./(g**2) + alpha_1)) * 0.5 * (1.5**2) * vals_i1_j0[n1p,n1,l1]  
      
      val += (- 0.30104661 )  * vals_i1_j0[n1p,n1,l1]
      
      val += (-1.4488847e-3 ) * vals_i2_j0[n1p,n1,l1]
      
      val += (4.9676959e-5 ) * vals_i3_j0[n1p,n1,l1]
      
    return val 

def radial_terms2p(l1,l2,l3,n1,n2,n3,n1p,n2p,n3p,vals_e0_i1_j0,vals_e1_i1_j0,vals_e0_i2_j0,vals_e1_i2_j0):
    
    val = 0.0
    if(n3 == n3p):
        
      val +=  (1.2790086e-2 ) * vals_e0_i1_j0[n1p,n1,l1]*vals_e1_i1_j0[n2p,n2,l2]
    
      val +=  (-5.5172502e-5) * (vals_e0_i2_j0[n1p,n1,l1]*vals_e1_i1_j0[n2p,n2,l2] + vals_e1_i2_j0[n2p,n2,l2]*vals_e0_i1_j0[n1p,n1,l1])
    
    return val
    
def radial_term(l1,l2,l3,n1,n2,n3,n1p,n2p,n3p,vals_e0_i1_j0,vals_e1_i1_j0):
    
    if(n3 == n3p):
        
      return vals_e0_i1_j0[n1p,n1,l1]*vals_e1_i1_j0[n2p,n2,l2]
    
    else:
      return 0.0

#angular term v1

 
            
def angular_term_zsq(l1,n1,l2,n2,l3,n3,l1p,n1p,l2p,n2p,l3p,n3p,vals_e0_i1_j0,vals_e1_i1_j0,vals_e0_i0_j1,vals_e1_i0_j1,vals_zsq):
    
    #z^2 term
      val = 0.0
    
      
      if((l3p == l3) and (n3p == n3)): #sanity check
          
       if((abs(l1p-l1) == 2 or l1p == l1) and (abs(l2p-l2) == 2 or l2p == l2)):
           
           val = vals_zsq[(l1p,l1,l2p,l2,l3)]
           
           if( l1p == l1):
               val *= vals_e0_i1_j0[(n1p,n1,l1)]
           elif(l1p - l1 == 2):
               val *= vals_e0_i0_j1[(n1,n1p,l1)]
           elif(l1p -l1 == -2):  
               val *= vals_e0_i0_j1[(n1p,n1,l1p)]
                   
               
               
               
           if( l2p == l2):
               val *= vals_e1_i1_j0[(n2p,n2,l2)]
           elif(l2p - l2 == 2):
               val *= vals_e1_i0_j1[(n2,n2p,l2)]   
           elif(l2p -l2 == -2):  
               val *= vals_e1_i0_j1[(n2p,n2,l2p)]
                   
       else:
           val *= 0
       '''  
       if((l1p-l1 == 2) and (l2p-l2 == 2)):
          
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_zsq[(l1p,l1,l2p,l2,l3)]   
      
       elif((l1p-l1 == -2) and (l2p-l2 == -2)):
          
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_zsq[(l1p,l1,l2p,l2,l3)]   
                            
       elif((l1p-l1 == 2) and (l2p-l2 == -2)):
          
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_zsq[(l1p,l1,l2p,l2,l3)]   
      
       elif((l1p-l1 == -2) and (l2p-l2 == 2)):
          
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_zsq[(l1p,l1,l2p,l2,l3)]   
           
        
       elif((l1p-l1 == 2) and (l2p == l2)):
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i1_j0[l2p,l2,n2p,n2]*vals_zsq[(l1p,l1,l2p,l2,l3)]                           
                             
       elif((l1p-l1 == -2) and (l2p == l2)):
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i1_j0[l2p,l2,n2p,n2]*vals_zsq[(l1p,l1,l2p,l2,l3)]                           
      
       elif((l2p-l2 == 2) and (l1p == l1)):
        val += vals_i1_j0[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_zsq[(l1p,l1,l2p,l2,l3)]                           
                             
       elif((l2p-l2 == -2) and (l1p == l1)):
        val += vals_i1_j0[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_zsq[(l1p,l1,l2p,l2,l3)]                           
      
       elif((l2p == l2) and (l1p == l1)):
        val += vals_i1_j0[l1p,l1,n1p,n1]*vals_i1_j0[l2p,l2,n2p,n2]*vals_zsq[(l1p,l1,l2p,l2,l3)] 

       else:
        val = 0.0                          
      
                              
       '''
    
      else: 
        val = 0.0
        print("l3' != l3 and n3' != n3")                     
                              
                                 
     
                                  
                                   
                               
                              
                                
                                
      return val

def angular_term_xysq(l1,n1,l2,n2,l3,n3,l1p,n1p,l2p,n2p,l3p,n3p,vals_e0_i1_j0,vals_e1_i1_j0,vals_e0_i0_j1,vals_e1_i0_j1,vals_xysq):
    
    #z^2 term
      val = 0.0
    
      
      if((l3p == l3) and (n3p == n3)): #sanity check
          
          
          
       if((abs(l1p-l1) == 2 or l1p == l1) and (abs(l2p-l2) == 2 or l2p == l2)):
           
           val = vals_xysq[(l1p,l1,l2p,l2,l3)]
           
           if( l1p == l1):
               val *= vals_e0_i1_j0[(n1p,n1,l1)]
           elif(l1p - l1 == 2):
               val *= vals_e0_i0_j1[(n1,n1p,l1)]
           elif(l1p -l1 == -2):  
               val *= vals_e0_i0_j1[(n1p,n1,l1p)]
               
           if( l2p == l2):
               val *= vals_e1_i1_j0[(n2p,n2,l2)]
           elif(l2p - l2 == 2):
               val *= vals_e1_i0_j1[(n2,n2p,l2)] 
           elif(l2p - l2 == -2):
               val *= vals_e1_i0_j1[(n2p,n2,l2p)] 
               
       else:
           val *= 0
          
          
          
      
       ''' 
       if((l1p-l1 == 2) and (l2p-l2 == 2)):
          
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_xysq[(l1p,l1,l2p,l2,l3)]   
      
       elif((l1p-l1 == -2) and (l2p-l2 == -2)):
          
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_xysq[(l1p,l1,l2p,l2,l3)]   
                            
       elif((l1p-l1 == 2) and (l2p-l2 == -2)):
          
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_xysq[(l1p,l1,l2p,l2,l3)]   
      
       elif((l1p-l1 == -2) and (l2p-l2 == 2)):
          
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_xysq[(l1p,l1,l2p,l2,l3)]   
           
        
       elif((l1p-l1 == 2) and (l2p == l2)):
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i1_j0[l2p,l2,n2p,n2]*vals_xysq[(l1p,l1,l2p,l2,l3)]                           
                             
       elif((l1p-l1 == -2) and (l2p == l2)):
        val += vals_i0_j1[l1p,l1,n1p,n1]*vals_i1_j0[l2p,l2,n2p,n2]*vals_xysq[(l1p,l1,l2p,l2,l3)]                           
      
       elif((l2p-l2 == 2) and (l1p == l1)):
        val += vals_i1_j0[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_xysq[(l1p,l1,l2p,l2,l3)]                           
                             
       elif((l2p-l2 == -2) and (l1p == l1)):
        val += vals_i1_j0[l1p,l1,n1p,n1]*vals_i0_j1[l2p,l2,n2p,n2]*vals_xysq[(l1p,l1,l2p,l2,l3)]                           
      
       elif((l2p == l2) and (l1p == l1)):
        val += vals_i1_j0[l1p,l1,n1p,n1]*vals_i1_j0[l2p,l2,n2p,n2]*vals_xysq[(l1p,l1,l2p,l2,l3)] 

       else:
        val = 0.0                          
       '''
                              
     
    
      else: 
        val = 0.0
        print("l3' != l3 and n3' != n3")                     
                              
                                 
     
                                  
                                   
                               
                              
                                
                                
      return val




                         



def get_gauge_invariant_ham(H,M,g,indices,e0,e1,vals_e0_i1_j0,vals_e0_i0_j1,vals_e0_i2_j0,vals_e0_i3_j0,vals_e1_i1_j0,vals_e1_i0_j1,vals_e1_i2_j0,vals_e1_i3_j0,vals_zsq,vals_xysq):

   time0 =time.time() 
   alpha_1 =  2.1810429e-2
   alpha_2 =  7.5714590e-3

   
   for i in range(0,M):
       
       
      
       multi_indexi= indices[i]
      
       lps = [multi_indexi[0],multi_indexi[1],multi_indexi[2]]
       nps = [multi_indexi[3],multi_indexi[4],multi_indexi[5]]
       
       l1p = multi_indexi[0]
       l2p = multi_indexi[1]
       l3p = multi_indexi[2]
       n1p = multi_indexi[3]
       n2p = multi_indexi[4]
       n3p = multi_indexi[5] 
       
      
       if(tuple([l1p,n1p]) == tuple([l2p,n2p]) == tuple([l3p,n3p])):

         norm_i = 1./2.
        
       elif((tuple([l1p,n1p]) == tuple([l2p,n2p]) )):
           
         norm_i = 1./2.
       
       elif((tuple([l1p,n1p]) == tuple([l3p,n3p])) or (tuple([l2p,n2p]) == tuple([l3p,n3p]))):  
         norm_i = 1./(np.sqrt(2.))
      
       else: 
         norm_i = 1./np.sqrt(2.)
      
      
       for j in range(0,M):
      
          val= 0.0 
          multi_indexj = indices[j]
         
          ls = [multi_indexj[0],multi_indexj[1],multi_indexj[2]]
          ns = [multi_indexj[3],multi_indexj[4],multi_indexj[5]]
         
          l1 = multi_indexj[0]
          l2 = multi_indexj[1]
          l3 = multi_indexj[2]
          n1 = multi_indexj[3]
          n2 = multi_indexj[4]
          n3 = multi_indexj[5]
       
          if(tuple([l1,n1]) == tuple([l2,n2]) == tuple([l3,n3])):

              norm_j = 1./2.
        
          elif((tuple([l1,n1]) == tuple([l2,n2]) )):
          
              norm_j = 1./2.
          
          elif ((tuple([l1,n1]) == tuple([l3,n3])) or (tuple([l2,n2]) == tuple([l3,n3]))):  
              norm_j = 1./(np.sqrt(2.))
      
          else: 
              norm_j = 1./np.sqrt(2.)   
         
         
          for [a,b] in list(permutations([0,1])):
              
                  l1p = lps[a]
                  l2p = lps[b]
                  l3p = lps[2]
                  n1p = nps[a]
                  n2p = nps[b]
                  n3p = nps[2]
                  
                  for [k,l] in list(permutations([0,1])):
              
                     l1 = ls[k]
                     l2 = ls[l]
                     l3 = ls[2]
                     n1 = ns[k]
                     n2 = ns[l]
                     n3 = ns[2]
                  
                  
              
         
         
        
         
                     #if( multi_indexi == multi_indexj):
                     if( (tuple([l1,n1]) == tuple([l1p,n1p])) and (tuple([l2,n2]) == tuple([l2p,n2p])) and (tuple([l3,n3]) == tuple([l3p,n3p]))):
                       val +=   (1./(1./(g**2) + alpha_1)) * kinetic_term(l1,l2,l3,n1,n2,n3,e0,e1,g)
                   
                 
                     if((l1 == l1p) and (l2 == l2p) and (l3 ==l3p)):  
                         
                       val += 0.5 * ((1./(g**2)) + alpha_2) * radial_term(l1,l2,l3,n1,n2,n3,n1p,n2p,n3p,vals_e1_i1_j0,vals_e1_i1_j0)
                   
                       val += 0.5 * ((1./(g**2)) + alpha_2) * radial_term(l1,l3,l2,n1,n3,n2,n1p,n3p,n2p,vals_e1_i1_j0,vals_e0_i1_j0)
                       val += 0.5 * ((1./(g**2)) + alpha_2) * radial_term(l2,l3,l1,n2,n3,n1,n2p,n3p,n1p,vals_e1_i1_j0,vals_e0_i1_j0)
                       
                       # two particle terms
                       val += radial_terms2p(l1,l2,l3,n1,n2,n3,n1p,n2p,n3p,vals_e1_i1_j0,vals_e1_i1_j0,vals_e1_i2_j0,vals_e1_i2_j0)
                       val += radial_terms2p(l1,l3,l2,n1,n3,n2,n1p,n3p,n2p,vals_e1_i1_j0,vals_e0_i1_j0,vals_e1_i2_j0,vals_e0_i2_j0)
                       val += radial_terms2p(l2,l3,l1,n2,n3,n1,n2p,n3p,n1p,vals_e1_i1_j0,vals_e0_i1_j0,vals_e1_i2_j0,vals_e0_i2_j0)
                       
                       
                      
                       # single particle terms
                       val += radial_terms1p(l1,l2,l3,n1,n2,n3,n1p,n2p,n3p,vals_e1_i1_j0,vals_e1_i2_j0,vals_e1_i3_j0,g)
                       val += radial_terms1p(l2,l1,l3,n2,n1,n3,n2p,n1p,n3p,vals_e1_i1_j0,vals_e1_i2_j0,vals_e1_i3_j0,g)
                       val += radial_terms1p(l3,l1,l2,n3,n1,n2,n3p,n1p,n2p,vals_e0_i1_j0,vals_e0_i2_j0,vals_e0_i3_j0,g)
                       
                       
                       val += (-1.2423581e-3 ) * vals_e1_i1_j0[n1p,n1,l1]*vals_e1_i1_j0[n2p,n2,l2]*vals_e0_i1_j0[n3p,n3,l3]
                       

                 
                     if((l3 == l3p) and (n3p == n3)):
                       val -= 0.5 * (1./(g**2) + alpha_2) * 3.*angular_term_zsq(l1,n1,l2,n2,l3,n3,l1p,n1p,l2p,n2p,l3p,n3p,vals_e1_i1_j0,vals_e1_i1_j0,vals_e1_i0_j1,vals_e1_i0_j1,vals_zsq)  
                       val -= 0.5 * (1./(g**2) + alpha_2) * 6.*angular_term_xysq(l1,n1,l2,n2,l3,n3,l1p,n1p,l2p,n2p,l3p,n3p,vals_e1_i1_j0,vals_e1_i1_j0,vals_e1_i0_j1,vals_e1_i0_j1,vals_xysq)
         
                     if((l2 == l2p) and (n2p == n2)):
                       val -= 0.5 * (1./(g**2) + alpha_2) * 3.*angular_term_zsq(l1,n1,l3,n3,l2,n2,l1p,n1p,l3p,n3p,l2p,n2p,vals_e1_i1_j0,vals_e0_i1_j0,vals_e1_i0_j1,vals_e0_i0_j1,vals_zsq)  
                       val -= 0.5 * (1./(g**2) + alpha_2) * 6.*angular_term_xysq(l1,n1,l3,n3,l2,n2,l1p,n1p,l3p,n3p,l2p,n2p,vals_e1_i1_j0,vals_e0_i1_j0,vals_e1_i0_j1,vals_e0_i0_j1,vals_xysq)
        
                     if((l1 == l1p) and (n1p == n1)):
                       val -= 0.5 * (1./(g**2) + alpha_2) * 3.*angular_term_zsq(l2,n2,l3,n3,l1,n1,l2p,n2p,l3p,n3p,l1p,n1p,vals_e1_i1_j0,vals_e0_i1_j0,vals_e1_i0_j1,vals_e0_i0_j1,vals_zsq)  
                       val -= 0.5 * (1./(g**2) + alpha_2) * 6.*angular_term_xysq(l2,n2,l3,n3,l1,n1,l2p,n2p,l3p,n3p,l1p,n1p,vals_e1_i1_j0,vals_e0_i1_j0,vals_e1_i0_j1,vals_e0_i0_j1,vals_xysq)
                  
                 
                 
                     else:
                       continue  
                  
                
          H[i,j] += val * (norm_i * norm_j)
               
   time1 = time.time()
   print("Time for constructing Hamiltonian: ",time1-time0)      
   return H   




def e1p_indices(lmax,nmax):
    Ls=list()
    Ns=list()
    l_vals=[]
    n_vals=[]
    
    time0 = time.time()
    for l1,l2,l3 in product(range(0,lmax+1),repeat= 3):
    
            
            if (abs(l1-l2) <= l3 <= (l1+l2) and (l1%2 == 0) and (l2%2 == 0) and (l3%2== 0)):
               #print("L1 L2 L3 =    ",l1,l2,l3 )
               l_vals=[l1,l2,l3]
               
               Ls.append(l_vals)

    
    
    
    
    for n1,n2,n3 in product(range(0,nmax),repeat= 3):
        n_vals = [n1,n2,n3]
        
        Ns.append(n_vals)
     
     
    
    indices_list=list()
    index=[]
    for ls in Ls:
      
     for ns in Ns:
         
        
         index = [ls[0],ls[1],ls[2],ns[0],ns[1],ns[2]]   
         
        
         
         
         
         
         
         for [a,b] in list(permutations([0,1])):
           
           
            perm_index = [ls[a],ls[b],ls[2],ns[a],ns[b],ns[2]]
            
            if(perm_index in indices_list):
               break
       
         else:
            indices_list.append(index)
        
        
            
       
    
    time1 = time.time()
    print("Time for computing unpermuted e1+ indices ",time1-time0)
    return indices_list











# build a dictionary of zeros
def zeros_dict(lmax,nmax,flag):
    
    z = np.zeros((lmax+1,nmax))
    if(flag == 1):
      for l in range(0,lmax+1):
        z[l] = Jn_zeros(l,nmax)[l]
    
    elif(flag == 0):
        for l in range(0,lmax+1):
          z[l] = rJnp_zeros(l,nmax)[l]
    return z/np.pi

def radial_func_norms(lmax,nmax,e):
    
    norms = np.zeros((lmax+1,nmax))
    
    for l in range(0,lmax+1):
        for n in range(0,nmax):
            
            val ,err = radial_integral(integrand(e[l,n],e[l,n],l,l,0,0))
            norms[l,n] = np.sqrt(val)
    return norms

def energies(indices_list,e0, e1,g):
    energies =list()
     
    for i in range(0,len(indices_list)):
        multi_indexi = indices_list[i]
       
        l1=multi_indexi[0]
        l2=multi_indexi[1]
        l3=multi_indexi[2]
        n1=multi_indexi[3]
        n2=multi_indexi[4]
        n3=multi_indexi[5]
        
        if( g < 1.2):
         E_i = (e1[l1,n1] + e1[l2,n2] + e0[l3,n3])
        else: 
         E_i = 0.5 * (e1[l1,n1]**2 + e1[l2,n2]**2 + e0[l3,n3]**2)
        #print("Index Energy    ", multi_indexi,E_i)
        energies.append(E_i)
         
        
        
    return energies  

def sorted_indices(Es,indices_list,M):
    sorted_indices=list()
    sort_index = np.argsort(Es)
    #sort_energies = np.sort(Es)
    #print(sort_energies[:M])
   
    sorted_indices = list()
    for index in sort_index:
        
       
       sorted_indices.append(indices_list[index])
    
     
    
    
    return sorted_indices[:M]


def simplify(H,M):
    c = 0.2
    for i in range(M):
        for j in range(M):
            if( i != j):
             if( abs(H[i][j])/abs(H[i][i]- H[j][j]) <  c):
                H[i][j] = 0.0
             else:
                 continue
            else:
                continue
    return H            
                

          
 
def get_eigenvalues(g,M,lmax,nmax):
    #indices=e1p_indices(lmax, nmax)
    
    
   
    with open('e1+_indices', 'rb') as fp:
     indices = pickle.load(fp)
    
    
    if ( g >= 1.2):
     with open('epsilon_e0','rb') as fp:
      e0 = pickle.load(fp)
     
     with open('epsilon_e1','rb') as fp:
      e1 = pickle.load(fp) 
   
    else:
     with open('epsilon_e0_typeB', 'rb') as fp:
         
         e0  = pickle.load(fp).reshape(21,10)
         
         
     e1  = np.loadtxt("epsilon_typeB_e_1_omega_1.5.csv", delimiter=',').reshape((21,10))
    
    
    Es = energies(indices, e0, e1,g)
    
    indices_sorted = sorted_indices(Es, indices, M)
    
    
    
    
    #print(indices_sorted)
    
    print("------------------------------------")
    #Ls=remove_perms(indices_sorted)
    
    #print(Ls)
   
    
   
    # loading 1-particle radial integrals
    if( g >= 1.2):
     with open('F_i1_j0_e0', 'rb') as fp:
        vals_e0_i1_j0 = pickle.load(fp)
     with open('F_i0_j1_e0', 'rb') as fp:
        vals_e0_i0_j1 = pickle.load(fp)
     with open('F_i2_j0_e0', 'rb') as fp:
        vals_e0_i2_j0 = pickle.load(fp)    
     with open('F_i3_j0_e0', 'rb') as fp:
        vals_e0_i3_j0 = pickle.load(fp)        
   
     with open('F_i1_j0_e1', 'rb') as fp:
        vals_e1_i1_j0 = pickle.load(fp)
     with open('F_i0_j1_e1', 'rb') as fp:
        vals_e1_i0_j1 = pickle.load(fp)
     with open('F_i2_j0_e1', 'rb') as fp:
        vals_e1_i2_j0 = pickle.load(fp)    
     with open('F_i3_j0_e1', 'rb') as fp:
        vals_e1_i3_j0 = pickle.load(fp)
    
    else:
      with open('F_i1_j0_e0_typeB', 'rb') as fp:
         
         vals_e0_i1_j0 = pickle.load(fp)
        
      with open('F_i0_j1_e0_typeB', 'rb') as fp:
        
         
        vals_e0_i0_j1 = pickle.load(fp)
        
      with open('F_i2_j0_e0_typeB', 'rb') as fp:
         
         vals_e0_i2_j0 = pickle.load(fp)  
        
   
      vals_e0_i3_j0  = np.loadtxt("flatFi3j0List.csv", delimiter=',').reshape((10, 10, 21))
    
    
      vals_e1_i0_j1  = np.loadtxt("flat_F_i0_j1_e1_typeB_omega1.5.csv", delimiter=',').reshape((10, 10, 19))
      vals_e1_i1_j0  = np.loadtxt("flat_F_i1_j0_e1_typeB_omega1.5.csv", delimiter=',').reshape((10, 10, 21))
      vals_e1_i2_j0  = np.loadtxt("flat_F_i2_j0_e1_typeB_omega1.5.csv", delimiter=',').reshape((10, 10, 21))
      vals_e1_i3_j0  = np.loadtxt("flat_F_i3_j0_e1_typeB_omega1.5.csv", delimiter=',').reshape((10, 10, 21))
    
    
    
    f1=open("Vals_Zsq2p","rb")
    
    f2=open("Vals_XYsq2p","rb")
    
    vals_zsq = pickle.load(f1)
    
    vals_xysq = pickle.load(f2)
    
    
    
    
    H = np.zeros((M,M))
    
    H = get_gauge_invariant_ham(H,M,g,indices_sorted,e0,e1,vals_e0_i1_j0,vals_e0_i0_j1,vals_e0_i2_j0,vals_e0_i3_j0,vals_e1_i1_j0,vals_e1_i0_j1,vals_e1_i2_j0,vals_e1_i3_j0,vals_zsq,vals_xysq)
    
    H = simplify(H,M)
    
    evals, evecs = np.linalg.eigh(H)
    
    print(isspmatrix(H))
    
    return H, evals
       

#lmax=12
#nmax=6
#M=1024
#g=0.9


#H,evals = get_eigenvalues(g,M,lmax,nmax)
#print(H.shape)
#print("# of non-zero elements: ",np.count_nonzero(H))
#print("Ground state energy :",evals[0])
#print("1st excited state energy :",evals[1])