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, # xy
                                         array([0,1j,0,1j,0])/sqrt(2.0),   # yz
                                         array([0,1,0,-1,0])/sqrt(2.0),    # xz
                                         array([1,0,0, 0,1])/sqrt(2.0),    # x2y2
                                         array([0,0,1,0,0 ])               # z2
                                         ]
                                        ),
                                  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 getCoefficientsCrystalF(vectorslmn, V0, V1,V2, alpha, dref  ):

    bondsVsp=[]
    bondsA  =[]

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

    cart2harm = Numeric.swapaxes( array([
                                         array([0,0,0,1,0,0,0])   ,            # 5 z3 - 3 z r2

                                         array([0,0, 1j,0,1j, 0, 0])/sqrt(2.0), # (3z z-r2) y
                                         array([0,0,  1,0,-1 , 0, 0])/sqrt(2.0), #   (3z z-r2) x
                                         
                                         array([0,-1,0,0,0,  1,0])/sqrt(2.0)/1.0j,   # z(xy)
                                         array([0, 1,0,0,0,  1,0])/sqrt(2.0),    # z(x2-y2)
                                         
                                         array([  1j,0,0,0,0,0,  1j])/sqrt(2.0),    # y3-2yx2
                                         array([  1 ,0,0,0,0,0, -1 ])/sqrt(2.0),    # x3-2xy2
                                         
                                         ]
                                        ),
                                  0,1
                                  )

    # (2 z3 - 3 z (x2+y2))**2    ===>  4* 3*2  + 9*( 4  ) =60
    # ( z(x2-y2) )**2            ===>  2*2  = 4
    #  (z(xy))**2                ====>  1
    # ((4 z z-x2-y2) y)             ===> 16*2+2+6   = 40
    #  (y3-3yx2)**2              ====>  6+9*2=24

    
    

    # (5 z3 - 3 z r2)  (5 zp3 - 3 zp r2)  = 25 * 3*2 * <zp,z>**3  -5 *3*3*<z, zp> * 2 *2   +  9 *(6<zp, z> +4*<z, zp>)
    #                                     = 50*3 zzp**3 -90 zzp                                                               

    # ((5 z z-r2 ) y) (5 zp3 - 3 zp r2)    = 25 *3<zpy> *2<zpz>**2 -5*3*<zpy>*2 -15 *<zpy>*2 +3*<zpy> *3*2 + 3*4*yzp
    #                                     = 150 yzp zzp**2 -30  yzp

    #   z(xy)             (5 zp3 - 3 zp r2)   = 5 *3<zpz> *2*<zpx> <zpy>
    #                                          =30 zzp xzp yzp
    
    #  z(x2-y2)           (5 zp3 - 3 zp r2)   = 5*3*<zpz>*2*(<zpx> <zpx>- <zpy> <zpy>)
    #                                         = 30 zzp ( xzp**2 - yzp**2 ) 

    #  x3-3xy2             (5 zp3 - 3 zp r2)   = 5*3*2<zpx>**3   -5*3*3*<zpx> *2<zpy>**2 -3*3*<zpx> *2  +9<zpx> *2 
    #                                          = 30 xzp **3 -90 xzpyzp**2

    #  y3-3yx2             (5 zp3 - 3 zp r2)   =  
    #                                          = 30 yzp **3 -90 yzpxzp**2


    #   (5 z3 - 3 z r2) (5 zp zp-r2) xp   =  25 *3<z xp> *2<zpz>**2 -5*3*<z xp>*2 -15 *<z xp>*2 +3*<z xp> *3*2+3*4*<z xp >
    #                                    =  150  zxp zzp**2         (-30-30+18 +12)zxp= 90 zxpzzp**2   -30 zxp

    #    (5 z z-r2) x    (5 zp zp-r2) xp   =   25 *( 2*<zpz>**2 <x xp>  + 2*<zxp> 2 <x zp> <z zp>)-5*<x xp>*2*2  +<x xp>*3*2+4*<x xp>
    #                                    =    25*(2 zzp**2 xxp+4 zzp zxp xzp   )    + xxp(-20+6+4) 
    
    #    (5z z-r2) y    (5zp zp-r2) xp   =   25 *( 2*<zpz>**2 <y xp>  + 2*<xpz> 2 <y zp> <z zp>)-5*(<y xp>*2)*2  +<y xp>*3*2+4*<y xp>
    #                                    =   25*( 2 zzp**2 yxp +4 yzp zxp zzp) -10 yxp
    
    #    z(xy)          (5 zp zp-r2) xp   =   5*  <z zp> (<x zp>  <y xp> + <y zp> <x xp>)
    #                                        +5* <x zp> (<y zp>  <z xp> + <z zp> <y xp>)
    #                                        +5* <y zp> (<x zp>  <z xp> + <z zp> <x xp>)
    #                                    =   5*2* (<z zp><x zp><y xp>+<x zp><y zp><z xp>+<y zp><z zp><x xp>) 
    
    #    z(x2-y2)       (5zp zp-r2) xp   =   5<z xp> *(2<x zp>**2-2<y zp>**2)  +5*2<z zp>2*(<x zp> < x xp>-<y zp> < y xp>)
    #                                    =  10  zxp *( xzp**2 -yzp**2)  + 20  zzp (  xzp xxp -yzp yxp) 

    
    #    x3-3xy2        (3zp zp-r2) xp   = 5*3<x xp> *2*<x zp>**2 
    #                                      -3<x xp> *2
    #                                      -15<x xp> *2<y zp>**2 -15*2<y xp>*2<x zp> <y zp>
    #                                      +3<x xp>*2           
    #                                    =  30 xxp xzp**2 -30 xxp yzp**2 -60   yxp xzp yzp     

    #    y3-3yx2        (5zp zp-r2) xp   =           
    #                                    =  30 yxp yzp**2 -30 yxp xzp**2 - 60   xxp yzp xzp     





    
    #    (5 z3 - 3 z r2)(5zp zp-r2) yp   =  15 *3<z yp> *2<zpz>**2 -5*3*<z yp>*2 -9 *<z yp>*2 +3*<z yp> *3*2+3*4*<z yp >=0
    #                                    = 0

    
    #    (5z z-r2) x    (5zp zp-r2) yp   =   9 *( 2*<zpz>**2 <x yp>  + 2*<ypz> 2 <x zp> <z zp>)-3*(<x yp>*2)*2  +<x yp>*3*2+4*<x yp>=0
    #                                    = 0
    
    #    (5z z-r2) y    (5zp zp-r2) yp   =   
    #                                    =  25* 2 zzp**2   -10 

    
    #    z(xy)          (5zp zp-r2) yp   =   5*  <z zp> (<x zp>  <y yp> + <y zp> <x yp>)
    #                                        +5* <x zp> (<y zp>  <z yp> + <z zp> <y yp>)
    #                                        +5* <y zp> (<x zp>  <z yp> + <z zp> <x yp>)
    #                                    =   5*2* (<z zp><x zp><y yp>) 

    
    #    z(x2-y2)       (5zp zp-r2) yp   =   5<z yp> *(2<x zp>**2-2<y zp>**2)  +5*2<z zp>2*(<x zp> < x yp>-<y zp> < y yp>)
    #                                    =   20  zzp (   -yzp yyp)

    
    #    x3-3xy2          (5 zp zp-r2) yp =  -60   yyp xzp yzp

    
    #    y3-3yx2        (5zp zp-r2) yp   =           
    #                                    =  30 yyp yzp**2 -30 yyp xzp**2      




    
    #              (5 z3 - 3 z r2)   zp(xpyp)  = 5 *3<z zp> *2*<z xp> <z yp>
    #                                          = 30  zzp zxp zyp  =0
    
    #               (5z z-r2) x    zp(xpyp)  =   5*  <zp z > (<xp z>  <yp x> + <yp z> <xp x>)
    #                                        +5* <xp z> (<yp z>  <zp x> + <zp z> <yp x>)
    #                                        +5* <yp z> (<xp z>  <zp x> + <zp z> <xp x>)
    #                                        =  0
    
    #               (5z z-r2) y    zp(xpyp)  =   5*  <zp z > (<xp z>  <yp y> + <yp z> <xp y>)
    #                                        +5* <xp z> (<yp z>  <zp y> + <zp z> <yp y>)
    #                                        +5* <yp z> (<xp z>  <zp y> + <zp z> <xp y>)
    #                                        = 10  zzp zxp

    
    #                z(x y )       zp(xpy )  = <zpz> <x xp> + <z xp> <x zp>
    #                z(x2-y2)      zp(xpy )  = -<z zp>2<y xp> -<z xp>2<y zp>
    #                 x3-3xy2      zp(xpy )  =-3*2(xzp yxp+xxp yzp )
    #                 y3-3yx2      zp(xpy )  = 3*2*yzp yxp -3*2 xzp xxp





    #              (5 z3 - 3 z r2)   zp(xp2-yp2)  = 5*3*<zpz>*2*(<zxp> <zxp>- <zyp> <zyp>)
    #                                             = 30*zzp*(zxp**2 - zyp**2 )
    
    #              (5z z-r2) x   zp(xp2-yp 2)   =   5<zp x> *(2<xp z>**2-2<yp z>**2)  +5*2<zp z>2*(<xp z> < xp x>-<yp z> < yp x)
    #                                           =10*xzp*(zxp**2-zyp**2)+20*zzp*(zxp*xxp-zyp*xyp) 

    
    #              (5z z-r2) y     zp(xp2-yp2)   =   5<zp y> *(2<xpz>**2-2<yp z>**2)  +5*2<z zp>2*(<xp z> < xp y>-<yp z> < yp y>)
    #                                            =   10*yzp*(zxp**2-zyp**2)+20*zzp*(zxp*yxp-zyp*yyp)          )

    
    #               z(xy )  zp(xp2-y2)           = zzp 2 xxpyxp+2 zxp(xzpyxp+ yzp xzp)
    #                                            
    #               z(x2-y2)         zp(xp2-y2) = zzp2 xxp**2 +4 zxpxzp xxp -(zzp2 yxp**2 +4 zxpyzp yxp)+2 zzp
    #                                           =  zzp*2*xxp**2 +4*zxp*xzp*xxp-(zzp*2*yxp**2 +4*zxp*yzp*yxp)+2*zzp
    
    #               x3-3xy2         zp(xp2-y2)           = 3 xzp 2 xxp**2     -3*2*xzp yxp**2 -3*2*2 xxpyzpyxp  +3 2 x zp
    #                                                    = 6*xzp*xxp**2     -6*xzp*yxp**2 -12*xxp*yzp*yxp  +6*xzp
    #               y3-3yx2         zp(xp2-y2)           = 3 yzp 2 yxp**2 - 3 2 y zp - 3 yzp 2 xxp**2-3 2 2 yxp(xzp xxp  )
    #                                                    =  6 yzp*yxp**2 - 6*yzp- 6*yzp*xxp**2-12*yxp*(xzp*xxp)
    

    

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

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

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


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


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

        zzp=arv[2]
        yzp=arv[1]
        xzp=arv[0]

        zyp=0
        yyp=1
        xyp=0

        zxp=-arv[0]
        yxp=0
        xxp= arv[2]
        
        #       z3                                          z2y                                   z2x                            
        Z3 = [  (50*3*zzp**3-90* zzp)/60.0,  (150.0*yzp*zzp**2-30*yzp)/sqrt(60.0*40.0) , (150.0*xzp*zzp**2-30.0*xzp)/sqrt(60.0*40.0),
        #                zxy                    zx2y2                              y3                     x3         
                (30*zzp*xzp*yzp)/sqrt(60.0),(30*zzp*(xzp**2-yzp**2))/sqrt(60.0*4),
        #                y3                                              x3         
           (30*yzp **3 -90*yzp *xzp**2)/sqrt(60.0*24.0), (30*xzp **3 -90*xzp *yzp**2)/sqrt(60.0*24.0) ]


        norm2=40.0                     
        #       z3                                          z2y                                   z2x
        Z2Y = [  (0.0)/sqrt(norm2*60.0),  (25*(2*zzp**2)-10)/sqrt(norm2*40.0) ,                               (0.0)/sqrt(norm2*40.0),
        #                zxy                                   zx2y2                                  
                (10.0*(zzp*xzp*yyp) )/sqrt(norm2) , (20*zzp*(-yzp* yyp))/sqrt(norm2*4),
        #                y3                                              x3         
           (30*yyp*yzp**2-30*yyp*xzp**2)/sqrt(norm2*24.0), ( -60*yyp*xzp*yzp)/sqrt(norm2*24.0) ]


        norm2=40.0                                      
        #       z3                                          z2y                                                           z2x
        Z2X =[(150*zxp*zzp**2-30*zxp)/sqrt(norm2*60.0),   (25*(2*zzp**2*yxp + 4*yzp*zxp*zzp) -10*yxp)/ sqrt(norm2*40.0),(25*(2*zzp**2*xxp+4*zzp*zxp*xzp)-10*xxp)/sqrt(norm2*40.0),
        #                zxy                                                zx2y2                              y3                     x3         
             (10*(zzp*xzp*yxp+xzp*yzp*zxp + yzp*zzp*xxp ))/sqrt(norm2),(10*zxp* (xzp**2-yzp**2)+20*zzp*(xzp*xxp-yzp*yxp))/sqrt(norm2*4),
        #                y3                                              x3         
           (30*yxp*yzp**2-30*yxp*xzp**2-60*xxp*yzp*xzp )/sqrt(norm2*24.0), (30*xxp*xzp**2-30*xxp*yzp**2-60*yxp*xzp*yzp)/sqrt(norm2*24.0) ]



        
        norm2=1.0
        #       z3                                          z2y                       z2x
        ZXY = [  (0)/sqrt(norm2*60.0),  (10*zzp*zxp)/sqrt(norm2*40.0) , (0)/sqrt(norm2*40.0),
        #                zxy                    zx2y2                              y3                     x3         
                (zzp*xxp + zxp*xzp)/sqrt(norm2),(-zzp*2*yxp -zxp*2*yzp)/sqrt(norm2*4),
        #                y3                                              x3         
           (6*yzp*yxp-6*xzp*xxp)/sqrt(norm2*24.0), (-3*2*(xzp*yxp+xxp*yzp))/sqrt(norm2*24.0) ]


        norm2=4.0
        #       z3                                          z2y                                   z2x
        ZX2Y2 = [  (30*zzp*(zxp**2 - zyp**2 ))/sqrt(norm2*60.0),  (10*yzp*(zxp**2-zyp**2)+20*zzp*(zxp*yxp-zyp*yyp))/sqrt(norm2*40.0) ,
                   (10*xzp*(zxp**2-zyp**2)+20*zzp*(zxp*xxp-zyp*xyp) )/sqrt(norm2*40.0),
        #                zxy                    zx2y2                              y3                     x3         
                ( zzp*2*xxp*yxp+2*zxp*(xzp*yxp+yzp*xzp))/sqrt(norm2),(zzp*2*xxp**2 +4*zxp*xzp*xxp-(zzp*2*yxp**2 +4*zxp*yzp*yxp)+2*zzp)/sqrt(norm2*4),
        #                y3                                              x3         
           ( 6* yzp*yxp**2 - 6*yzp- 6*yzp*xxp**2-12*yxp*(xzp*xxp))/sqrt(norm2*24.0), ( 6*xzp*xxp**2     -6*xzp*yxp**2 -12*xxp*yzp*yxp  +6*xzp)/sqrt(norm2*24.0) ]



        Z3=array(Z3)
        Z2X=array(Z2X)
        Z2Y=array(Z2Y)
        ZXY=array(ZXY)
        ZX2Y2=array(ZX2Y2)

##         print "z3 z3 ",  sum(Z3*Z3)
##         print "Z2X*Z2X ",  sum(Z2X*Z2X)
##         print "Z2Y*Z2Y ",  sum(Z2Y*Z2Y)
##         print " ZXY*ZXY",  sum(ZXY*ZXY)
##         print " ZX2Y2*ZX2Y2",  sum(ZX2Y2*ZX2Y2)
        
##         print " Z3*ZX2Y2",  sum(Z3*ZX2Y2)
##         print " Z3*Z2Y",  sum(Z3*Z2Y)
##         print " Z3*ZXY",  sum(Z3*ZXY)
##         print " Z3*Z3",  sum(Z3*Z3)
##         print " Z3*Z2X ",  sum(Z3*Z2X),"\n\n"

##         print " ZX2Y2*ZX2Y2",  sum(ZX2Y2*ZX2Y2)
##         print " ZX2Y2*Z2Y",  sum(ZX2Y2*Z2Y)
##         print " ZX2Y2*ZXY",  sum(ZX2Y2*ZXY)
##         print " ZX2Y2*Z3",  sum(ZX2Y2*Z3)
##         print " ZX2Y2*Z2X ",  sum(ZX2Y2*Z2X),"\n\n"


##         print " Z2X*ZX2Y2",  sum(Z2X*ZX2Y2)
##         print " Z2X*Z2Y",  sum(Z2X*Z2Y)
##         print " Z2X*ZXY",  sum(Z2X*ZXY)
##         print " Z2X*Z3",  sum(Z2X*Z3)
##         print " Z2X*Z2X ",  sum(Z2X*Z2X),"\n\n"


##         print " Z2Y*ZX2Y2",  sum(Z2Y*ZX2Y2)
##         print " Z2Y*Z2Y",  sum(Z2Y*Z2Y)
##         print " Z2Y*ZXY",  sum(Z2Y*ZXY)
##         print " Z2Y*Z3",  sum(Z2Y*Z3)
##         print " Z2Y*Z2X ",  sum(Z2Y*Z2X),"\n\n"
        


        add =    Z3[:,None]*Z3    *  V0*(dist/dref)**alpha
        add =   add+  Z2X[:,None]*Z2X    *  V1*(dist/dref)**alpha
        add =   add+  Z2Y[:,None]*Z2Y    *  V1*(dist/dref)**alpha
        add =   add+  ZXY[:,None]*ZXY    *  V2*(dist/dref)**alpha
        add =   add+  ZX2Y2[:,None]*ZX2Y2*  V2*(dist/dref)**alpha

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

        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
