import Numeric
from Numeric import array
from  math import *
import math

def get_5_vettori( bonds   ):
    # bonds e una lista dove ogni items contiene l,m,n, Vs, Vp
    vect_list=[]
    nonprojected_amplitudes=[]
    for do in ["xy","yz","xz", "x2y2","z2" ]:
        newvect=[]
        for l, m, n  ,  Vs, Vp  in bonds:
            for po in ["x","y","z"    ]:
                newvect.append(    hoppingPD( po, do, l, m, n  ,  Vs, Vp     )          ) 
        vect_list.append(Numeric.array(newvect))
    return Numeric.array(vect_list)

def get_3_vettori( bonds   ):
    # bonds e una lista dove ogni items contiene l,m,n, A
    vect_list=[]
    nonprojected_amplitudes=[]
    for pA in ["x","y","z"]:
        newvect=[]
        for l, m, n  , As, Ap in bonds:
            for po in ["x","y","z"    ]:
	        print " appendo per ", As, Ap
                newvect.append(    dipolePP( pA, po, l, m, n  ,  As, Ap    )          ) 
        vect_list.append(Numeric.array(newvect))
    return Numeric.array(vect_list)


def GramSchmidt(vectlist_orig):
    vectlist=vectlist_orig.copy()
    nv = len(vectlist)
    for i in range(nv):
        normi = Numeric.sum( vectlist[i]*vectlist[i])
        if normi > 1.0e-10:
            for j in range(i+1, nv):
                normj = Numeric.sum( vectlist[j]*vectlist[j])
                if( normj>1.0e-10):
                    s= Numeric.sum( vectlist[i]*vectlist[j]) /normj
                    vectlist[i] = vectlist[i] -s*vectlist[j]
            norma=Numeric.sum(vectlist[i]*vectlist[i])
            if( norma>1.0e-10):
                print " quadrato di vectlist[",i,"] est ", norma
                vectlist[i] /=math.sqrt(norma)
    return vectlist
        
def get_newamps_and_newdipolamps(vectlist, dipolamp):
    print vectlist
    newvects=GramSchmidt(vectlist)
    print newvects
    newamps = Numeric.dot( newvects, Numeric.swapaxes(vectlist,0,1)    )
    newdipolamps = Numeric.dot( newvects, Numeric.swapaxes(dipolamp,0,1)    )
    return newamps,newdipolamps
    



def hoppingPD( po, do, l, m, n  ,  Vs, Vp     ):
    # giusto per correttezza
    # PD prende vettori che vanno da P a D
    # Per comodita sara meglio usarli da D a P
    l=-l
    m=-m
    n=-n
    
    caso=po+"-"+do
    if caso=="x-xy" :
        res =m*Vp*(1-2*pow(l,2))+(m*Vs*pow(l,2)*sqrt(3))
    elif caso=="x-yz":
        res =-2*l*m*n*Vp+(l*m*n*Vs*sqrt(3))
    elif caso=="x-xz":
        res =n*Vp*(1-2*pow(l,2))+(n*Vs*pow(l,2)*sqrt(3))
        
    elif caso=="x-x2y2":
        res =-(l*Vp*(-2+2*pow(l,2)+pow(n,2)))+(l*Vs*(pow(l,2)-pow(m,2))*sqrt(3))/2.
    elif caso=="x-z2":
        res =(l*Vs*(-1+3*pow(n,2)))/2.-l*Vp*pow(n,2)*sqrt(3)
        
    elif caso=="y-xy" :
        res =l*Vp*(1-2*pow(m,2))+(l*Vs*pow(m,2)*sqrt(3));
    elif caso=="y-yz" :
        res =n*Vp*(1-2*pow(m,2))+(n*Vs*pow(m,2)*sqrt(3));
    elif caso=="y-xz" :
        res =-2*l*m*n*Vp+(l*m*n*Vs*sqrt(3));
        
    elif caso=="y-x2y2" :
        res =m*Vp*(-2+2*pow(m,2)+pow(n,2))+(m*Vs*(pow(l,2)-pow(m,2))*sqrt(3))/2.;
    elif caso=="y-z2" :
        res =(m*Vs*(-1+3*pow(n,2)))/2.-m*Vp*pow(n,2)*sqrt(3);
        
    elif caso=="z-xy" :
        res =-2*l*m*n*Vp+(l*m*n*Vs*sqrt(3));
    elif caso=="z-yz" :
        res =m*Vp*(1-2*pow(n,2))+(m*Vs*pow(n,2)*sqrt(3));
    elif caso=="z-xz" :  
        res =l*Vp*(1-2*pow(n,2))+(l*Vs*pow(n,2)*sqrt(3));
        
    elif caso=="z-x2y2" :  
        res =-((l-m)*(l+m)*n*Vp)+(n*Vs*(pow(l,2)-pow(m,2))*sqrt(3))/2.;
    elif caso=="z-z2" :  
        res =(n*Vs*(-1+3*pow(n,2)))/2.-n*Vp*(-1+pow(n,2))*sqrt(3);
    else:
        print " caso ", caso , " non riconosciuto"
        raise " caso non riconosciuto "
    return res 

def dipolePP( pA, po, l, m, n  ,  Vs, Vp   ):

  case = pA+po

  if   case=="xx" : res =Vp-Vp*pow(l,2)+Vs*pow(l,2);
  elif case=="xy" : res =-(l*m*Vp)+l*m*Vs;
  elif case=="xz" : res =-(l*n*Vp)+l*n*Vs;
  elif case=="yx" : res =-(l*m*Vp)+l*m*Vs;
  elif case=="yy" : res =Vp-Vp*pow(m,2)+Vs*pow(m,2);
  elif case=="yz" : res =-(m*n*Vp)+m*n*Vs;
  elif case=="zx" : res =-(l*n*Vp)+l*n*Vs;
  elif case=="zy" : res =-(m*n*Vp)+m*n*Vs;
  elif case=="zz" : res =Vp*(pow(l,2)+pow(m,2))+Vs*pow(n,2);
  else  :
	  print " non corrisponde a niente "
	  raise " STOP"
  return res

def getCoefficientsCrystal(vectorslmn, V0, V1, alpha, dref  ):

    bondsVsp=[]
    bondsA  =[]

    res= Numeric.zeros([5,5],"d")

    cart2harm = Numeric.swapaxes( array([array([-1,0,0,0,1])/sqrt(2.0)/1j,
                                         array([0,1j,0,1j,0])/sqrt(2.0),
                                         array([0,1,0,-1,0])/sqrt(2.0),
                                         array([1,0,0, 0,1])/sqrt(2.0),
                                         array([0,0,1,0,0 ])
                                         ]
                                        ),
                                  0,1
                                  )

    harm2cart = Numeric.conjugate( Numeric.swapaxes(cart2harm,0,1))

    jmenoi = Numeric.arrayrange(5)
    
    jmenoi = (  -jmenoi  )[ : , None ]  +   jmenoi

    for b in vectorslmn:
        arv=Numeric.array(b,"d")
        dist=Numeric.sqrt(Numeric.sum(arv*arv))
        arv=(arv/dist)

        Z2 = array( [   sqrt(3.0)*arv[0]*arv[1] , sqrt(3.0)*arv[2]*arv[1] , sqrt(3.0)*arv[0]*arv[2] ,
                        sqrt(3.0)/2.0*(arv[0]*arv[0]-arv[1]*arv[1]) ,   0.5*(3*arv[2]*arv[2]-1)            ] )

        add =    Z2[:,None]*Z2    *  V0*(dist/dref)**alpha
        
        add  = Numeric.dot(   add        ,    harm2cart      )
        add  = Numeric.dot(  cart2harm   ,    add            )
        
        res=   res   +   add 


        omega=math.atan2(arv[1],arv[0])


        arv=array([            math.sqrt( arv[0]*arv[0] +  arv[1]*arv[1])  ,    0.0   ,    arv[2]    ]   )

        
        XZ = array( [        0       ,     0           ,              arv[2]*arv[2]  - arv[0]*arv[0]        ,
                          arv[0]*arv[2] ,   -math.sqrt(3.0) *   arv[0]*arv[2]      ] )
        
        YZ = array( [        arv[0]       ,               arv[2]           ,                 0         ,
                               0        ,                0                                                             ] )

        add =         XZ[:,None]*XZ    *  V1*(dist/dref)**alpha
        add =   add+  YZ[:,None]*YZ    *  V1*(dist/dref)**alpha

        
        add  = Numeric.dot(   add        ,    harm2cart      )
        add  = Numeric.dot(  cart2harm   ,    add            )
        add  = add *       Numeric.exp(      jmenoi *1.0j * omega     )

        res=res+add

    return  res

    
def getCoefficients(vectorslmn, Vs, Vp, alphasp, Ps, Pp,  alphaA , dref  ):
    print " Vp est ", Vp
    bondsVsp=[]
    bondsA  =[]
    for b in vectorslmn:
        arv=Numeric.array(b,"d")
        dist=Numeric.sqrt(Numeric.sum(arv*arv))
        arv=(arv/dist).tolist()
        bondsVsp.append(  arv +[  Vs*(dist/dref)**alphasp, Vp*(dist/dref)**alphasp           ]           )
        bondsA  .append(  arv +[  Ps*(dist/dref)**alphaA   ,  Pp*(dist/dref)**alphaA     ]           )

    vettori5 = get_5_vettori(bondsVsp)
    vettori3 = get_3_vettori(bondsA)
    newamps,newdipolamps=get_newamps_and_newdipolamps(vettori5, vettori3 )

    cart2harm = Numeric.swapaxes( array([array([-1,0,0,0,1])/sqrt(2.0)/1j,
                                         array([0,1j,0,1j,0])/sqrt(2.0),
                                         array([0,1,0,-1,0])/sqrt(2.0),
                                         array([1,0,0, 0,1])/sqrt(2.0),
                                         array([0,0,1,0,0 ])
                                         ]
                                        ),
                                  0,1
                                  )

    harm2cart = Numeric.conjugate( Numeric.swapaxes(cart2harm,0,1))

    newamps = Numeric.dot(   newamps,    harm2cart      )

    return  newamps,newdipolamps



if __name__=="__main__":
    vects=[]
    for i in [-1,1]:
        for j in [0,1,2]:
                add=[0,0,0]
                add[j]=i+1.0/10
                vects.append(add)
    Vs=1.0
    Vp=.0
    alphasp=-3.0
    alphaA =-2.0
    A=1
    dref=sqrt(1.0)
    # nomi, newamps,newdipolamps = getCoefficients(vects , Vs, Vp, alphasp, A, alphaA , dref  )
    sposta=[[-1.0, 0, 0], [1.1, 0, 0], [0, -1.0, 0], [0, 1.1, 0], [0, 0, -1.0], [0, 0, 1.1]]
    spostb=[[-1.0, 0, -0.10000000000000001], [1.0, 0, -0.10000000000000001], [0, -1.0, -0.10000000000000001], [0, 1.0, -0.10000000000000001], [0.070000000000000007, 0.070000000000000007, -1.0], [0.070000000000000007, 0.070000000000000007, 1.0]]
    
    newamps,newdipolamps = getCoefficients(spostb,  Vs=1, Vp=sqrt(3.0)/2,    alphasp= -3.0,  Ps=1.0, Pp=0.,alphaA=-2 , dref=1.0 )

    print " ########### " 
    print newamps
    print newdipolamps
