import Hilbertxx
import Sparsa
import Numeric
import math
import time
import string
import math
import copy
import random
 
usecomplex = 1
usepython  = 1


import os
place=os.path.dirname(__file__)
import sys

sys.path.append(place+"/../../python")

from matrixop import *


def estraigamma(gamma, enes):
  if type(gamma)==type(1.0):
    return gamma
  else:
    g1,e12,g2 = gamma
    res=g1*Numeric.less(enes, e12)+g2*Numeric.less(e12,enes)
    return res


def changename(name, c):
  if c==1:
    return name
  elif c==2:
    return string.replace(name,"base","exci")
  else:
    return string.replace(name,"base","final")

  
def ReadS( name )  :
  t_base    = Hilbertxx.Hxx_TransitionMatrix()
  t_base.Read_add(name,1.0,0)
  ars_base=t_base   .create_3Array_and_Clean()
  S_base=Sparsa.Sparsa3A( ars_base[2], ars_base[0], ars_base[1] )
  return S_base

def SfromT( t_base )  :
  if t_base is None:
    return None
  ars_base=t_base   .create_3Array_and_Clean()
  S_base=Sparsa.Sparsa3A( ars_base[2], ars_base[0], ars_base[1] )
  return S_base

 

class System:
  def __init__(self , basenames, basepars , excinames, excipars,  finnames, finpars,  calcnames, calcpars,   quadrunames, angularnames, counternames):
    self.basenames=basenames
    self.excinames=excinames
    self.finnames=finnames
    self.calcnames=calcnames
    self.quadrunames = quadrunames
    self.angularnames  = angularnames
    self.counternames=counternames

    
    for nomelist, vallist in zip( [basenames, excinames, finnames, calcnames ], [ basepars, excipars, finpars,calcpars ] ):
      for name,par in zip(  nomelist, vallist ):
        setattr(self,name,par)


  def GetESLcounters(self, composeonly=1, fvect=None):


    if fvect is   None:

      H=self.GetH(composeonly)
      if( not usepython):
        vects,vals = spectrumFunction.Diago(H[0],0.1,self.nsearchedeigen, self.shift_invert)
      else:
        if usecomplex:
          H = TrialMatrixC(*H )
        else:
          H=H[0]
        vals , vects = solveEigenSystem( H , self.nsearchedeigen)
    else:
        vals, vects = [0.0], [fvect]

    if( usepython):
      # class4sparse = H.__class__
      class4sparse = TrialMatrixC
      class4vect   = class4sparse.getClass4Vect()
      if fvect is None:
        dim = H.dim
      else:
        dim = vects[0].len()
      L2s=[]
      S2s=[]
      SLs=[]
      opslist=[]
      for  opsnames in [self.angularnames[:3], self.angularnames[3:]]      :
        ops = [ ReadS( self.case+"/" + changename(name, composeonly) )      for name in opsnames       ]
        opslist.append(ops)
      for l,(As,Bs) in zip(    [L2s, S2s, SLs ] ,           [ [opslist[0]]*2, [opslist[1]]*2,   [opslist[0], opslist[1]] ]):
        for v in vects:
          res=0.0
          dum1 = class4vect(dim)
          dum2 = class4vect(dim)
          for (A,B), fact in zip( zip(As,Bs    )  , [1.0,0.5,0.5]):

            if(usecomplex):
              A=TrialMatrixC(A,None)
              B=TrialMatrixC(B,None)
            dum1.set_to_zero()
            dum2.set_to_zero()
            
            A.Moltiplica(dum1, v)
            B.Moltiplica(dum2, v)
            add=class4vect.scalare(dum1, dum2)
            res=res+add*fact
          l.append(res)
      cntlist=[]
      cntlist = [ ReadS( self.case+"/" + changename(name, composeonly) )      for name in self.counternames      ]
      cresListS=[]
      for cnt in cntlist:
        if(usecomplex):
          cnt=TrialMatrixC(cnt,None)
        cresList=[]
        for v in vects:
          dum1.set_to_zero()
          cnt.Moltiplica(dum1, v)
          res=class4vect.scalare(v, dum1)
          cresList.append(res)
        cresListS.append(cresList)

    else:

  
      L2s=[]
      S2s=[]
      SLs=[]
      opslist=[]
      for  opsnames in [self.angularnames[:3], self.angularnames[3:]]      :
        ops = [ ReadS( self.case+"/" + changename(name, composeonly) )      for name in opsnames       ]
        opslist.append(ops)
      for l,(As,Bs) in zip(    [L2s, S2s, SLs ] ,           [ [opslist[0]]*2, [opslist[1]]*2,   [opslist[0], opslist[1]] ]):
        for v in vects:
          res=0.0
          for (A,B), fact in zip( zip(As,Bs    )  , [1.0,0.5,0.5]):
            v0=Numeric.zeros(len(vects[0]),"d")
            v1=Numeric.zeros(len(vects[0]),"d")
            A.Moltiplica(v0, v)
            B.Moltiplica(v1, v)
            add=Numeric.sum(v1*v0)
            res=res+add*fact
          l.append(res)

      cntlist=[]
      cntlist = [ ReadS( self.case+"/" + changename(name, composeonly) )      for name in self.counternames      ]
      cresListS=[]
      for cnt in cntlist:
        cresList=[]
        for v in vects:
          v0=Numeric.zeros(len(vects[0]),"d")
          cnt.Moltiplica(v0, v)
          res=Numeric.sum(v*v0)
          cresList.append(res)
        cresListS.append(cresList)
    return [ vals, L2s, S2s, SLs ] + cresListS
                          
                          

  def GetSpectrum(self, polarisation=None):
    return self.Get( polarisation=polarisation, composeonly=0, )
  
  def GetH(self,composeonly=1):
    return self.Get( composeonly=composeonly)

  def GetRIXS(self, composeonly=0, polarisationIn=None,polarisationOut=None,
                         ein=None, eout1=None, eout2=None, dout=None, gammain=None, gammaout=None):

    calcpars=[  getattr(self,name) for name in  self.calcnames ]
    basenames = self.basenames
    basecoeffs = [  getattr(self,name) for name in  basenames ]
    Riduci( basecoeffs, basenames , self.calcnames, calcpars )

    excinames = self.excinames
    excicoeffs = [  getattr(self,name) for name in  excinames ]
    Riduci(excicoeffs ,  excinames, self.calcnames, calcpars )



    finnames = self.finnames
    fincoeffs = [  getattr(self,name) for name in  finnames ]
    Riduci(fincoeffs ,  finnames, self.calcnames, calcpars )



    result=self.spectrumRIXS("%s/"%self.case+"/", basecoeffs, basenames,
                         excicoeffs, excinames, fincoeffs, finnames, polarisationIn=polarisationIn,polarisationOut=polarisationOut,
                         ein=ein, eout1=eout1, eout2=eout2, dout=dout, gammain=gammain, gammaout=gammaout,
                         HFe=self.hoppfactexci, HFf=self.hoppfactfinal)
    return result

 
  def Get(self, composeonly=0, polarisation=None):
    calcpars=[  getattr(self,name) for name in  self.calcnames ]    
    basenames = self.basenames
    basecoeffs = [  getattr(self,name) for name in  basenames ]
    Riduci( basecoeffs, basenames , self.calcnames, calcpars )
    
    excinames = self.excinames
    excicoeffs = [  getattr(self,name) for name in  excinames ]
    Riduci(excicoeffs ,  excinames, self.calcnames, calcpars )



    finnames = self.finnames
    fincoeffs = [  getattr(self,name) for name in  finnames ]
    Riduci(fincoeffs ,  finnames, self.calcnames, calcpars )



    if(composeonly):
      
      t_base    = Hilbertxx.Hxx_TransitionMatrix()
      t_baseI = Hilbertxx.Hxx_TransitionMatrix()
      if( composeonly==1):
        self.compose( t_base, t_baseI,        "%s/" % self.case+"/", basecoeffs, basenames, "base"   )
      else:
        self.compose( t_base, t_baseI,        "%s/" % self.case+"/", excicoeffs, excinames, "base"   )
      ars_base=t_base   .create_3Array_and_Clean()

      if( usecomplex or  not usepython ) :
        S_base=Sparsa.Sparsa3A( ars_base[2], ars_base[0], ars_base[1] )
      else:
        S_base=TrialMatrix( ars_base[2], ars_base[0], ars_base[1] )

      
      ars_base=t_baseI   .create_3Array_and_Clean()
      if( usecomplex or  not usepython ) :
        S_baseI=Sparsa.Sparsa3A( ars_base[2], ars_base[0], ars_base[1] )
      else:
        S_baseI=TrialMatrix( ars_base[2], ars_base[0], ars_base[1] )
      return S_base, S_baseI
    
    result=self.spectrum("%s/"%self.case+"/", basecoeffs, basenames,
                         excicoeffs, excinames, fincoeffs, finnames, polarisation=polarisation,
                         HFe=self.hoppfactexci, HFf=self.hoppfactfinal)
    return result
  


  
  def compose(self, t_base  , t_baseI,  repe, coeff, names, prefix , hoppfact=1.0)  :

    self.get_nomi_coeffs()
    for i in range(len(coeff) ):
      if( coeff[i] !=0.0) :
        if type(coeff[i])!=type(1.0j):
          t_base.Read_add( repe+names[i], coeff[i], 0)
        else:
          if coeff[i].real:
            t_base.Read_add( repe+names[i], coeff[i].real, 0)
          if coeff[i].imag:
            t_baseI.Read_add( repe+names[i], coeff[i].imag, 0)



##     self.get_nomi_coeffs()
##     for i in range(len(coeff) ):
##       if( coeff[i] !=0.0) :
##         print " leggo " , repe+names[i]
##         t_base.Read_add( repe+names[i], coeff[i], 0)




    countername = names[0]
    pos = string.find(countername,"_")
    countername=countername[:pos+1]+"counterDL"
    t_counter = Hilbertxx.Hxx_TransitionMatrix()
    t_counter.Read_add(   repe+ countername , 1.0, 0      )
    fc, tc, cc = t_counter.create_3Array_and_Clean()
    cc=(10.0-cc).astype("i")

    if  self.facts_hop is not None and  int( max(cc)) !=len(self.facts_hop)-1:
      print max(cc)

      raise " self.facts_hop non corrisponed in lunghezza al numero di stati di valenza "


    for n,c in zip(self.grapnomi+self.crystalnomi, self.grapamps+self.crystalamps):
       if( c.real !=0.0):
        print " AGGOINGO A PARTE REALE ",  repe+prefix+n, c.real
        if  self.facts_hop is not None :
          t_base.Read_add_forhops( repe+prefix+n, c.real*hoppfact, 0,cc, self.facts_hop  )
        else:
          t_base.Read_add( repe+prefix+n, c.real*hoppfact, 0 )
       if(c.imag!=0.0):
        print " AGGOINGO A PARTE IMMAGINARIA ",  repe+"anti"+prefix+n, c.imag
        if  self.facts_hop is not None :
          t_baseI.Read_add_forhops( repe+"anti"+prefix+n, c.imag*hoppfact, 0,cc, self.facts_hop )
        else:
          t_baseI.Read_add( repe+"anti"+prefix+n, c.imag*hoppfact, 0)



    
  def get_nomi_coeffs(self):
    from miscele import getCoefficients
    (self.BONDS, self.Vs, self.Vp, self.ALPHAVSP, 1.0, 0.5,  self.ALPHADIPO, self.DREF)
    grapamps,self.grapdipolamps = getCoefficients(self.BONDS, self.Vs, self.Vp, self.ALPHAVSP, self.Dips,self.Dipp,  self.ALPHADIPO, self.DREF)
    
    print " #################################### "
    print self.BONDS, self.Vs, self.Vp, self.ALPHAVSP, 1.0, self.ALPHADIPO, self.DREF
    print grapamps
    print self.grapdipolamps
    
    self.grapnomi=[]
    self.grapamps=[]
    for i in range(5):
      for j in range(5):
        self.grapnomi.append(  "_mixer_"+str(-2+i)+"_"+str(-2+j)  )
        self.grapamps.append(grapamps[i][j])
    print self.grapnomi
    print self.grapamps


    from miscele import getCoefficientsCrystal
    crystalamps  = getCoefficientsCrystal(self.BONDS, self.VC0, self.VC1, self.ALPHAVC,  self.DREF)

    self.crystalnomi=[]
    self.crystalamps=[]
    for i in range(5):
      for j in range(i,5):
        self.crystalnomi.append(  "_CRYSTAL_1_%d_%d_"%(j,i) )
        self.crystalamps.append(crystalamps[j][i])
     

  def spectrum(self, repe, coeff,  names,coeff_exc, names_exc,coeff_fin, names_fin,  polarisation=None, HFe = 1.0, HFf=1.0):



    t_base = Hilbertxx.Hxx_TransitionMatrix()
    t_baseI = Hilbertxx.Hxx_TransitionMatrix()
    
    self.compose( t_base, t_baseI,  repe,  coeff, names , "base" )
        
    t_excited = Hilbertxx.Hxx_TransitionMatrix()
    t_excitedI = Hilbertxx.Hxx_TransitionMatrix()
    
    self.compose( t_excited, t_excitedI,  repe, coeff_exc, names_exc ,"exci", hoppfact=HFe) 

##    t_final = Hilbertxx.Hxx_TransitionMatrix()
##    t_finalI = Hilbertxx.Hxx_TransitionMatrix()
    
##    self.compose( t_final, t_finalI,  repe, coeff_fin, names_fin ,"final", hoppfact=HFf) 


        
    ars_base=t_base.create_3Array_and_Clean()
    ars_ecci=t_excited.create_3Array_and_Clean()

    if( usecomplex or not usepython ):
      S_base=Sparsa.Sparsa3A( ars_base[2], ars_base[0], ars_base[1] )
    else:
      S_base=TrialMatrix( ars_base[2], ars_base[0], ars_base[1] )

    if( usecomplex or not usepython ):
      S_ecci=Sparsa.Sparsa3A( ars_ecci[2], ars_ecci[0], ars_ecci[1] )
    else:
      S_ecci=TrialMatrix( ars_ecci[2], ars_ecci[0], ars_ecci[1] )


    ars_base=t_baseI.create_3Array_and_Clean()
    ars_ecci=t_excitedI.create_3Array_and_Clean()

    S_baseI=Sparsa.Sparsa3A( ars_base[2], ars_base[0], ars_base[1] )

    S_ecciI=Sparsa.Sparsa3A( ars_ecci[2], ars_ecci[0], ars_ecci[1] )


    operatorS=[]
    if(usepython  and not polarisation  ):
      nomebase="mixedDipole_couches_0_2_mz_"
      for pol in range(3):
        t_dip = Hilbertxx.Hxx_TransitionMatrix() 
        for finale in range(5):
          print " PER I DIPOLI AGGIUNGO "  , repe+nomebase+str(-2+finale), self.grapdipolamps[finale][pol]
          t_dip.Read_add( repe+nomebase+str(-2+finale), self.grapdipolamps[finale][pol], 0)
          
        ars = t_dip.create_3Array_and_Clean()
        if( usecomplex):
          operatorS.append(TrialMatrixC(Sparsa.Sparsa3A( ars[2], ars[0], ars[1] ), None)  )
        else:
          operatorS.append(Sparsa.Sparsa3A( ars[2], ars[0], ars[1] ) )
        
          
 ##    elif polarisation is not None:
      
##       t_dip = Hilbertxx.Hxx_TransitionMatrix() 
##       for opname, fact  in zip(self.quadrunames, polarisation):
##         t_dip.Read_add( repe+opname, fact, 0)
##       ars = t_dip.create_3Array_and_Clean()
##       if( usecomplex):
##         operatorS.append(TrialMatrixC(Sparsa.Sparsa3A( ars[2], ars[0], ars[1] ), None)  )
##       else:
##         operatorS.append(Sparsa.Sparsa3A( ars[2], ars[0], ars[1] ) )



    elif polarisation is not None:
      
      t_dipR=Hilbertxx.Hxx_TransitionMatrix()
      t_dipI=None

      for opname, fact  in zip(self.quadrunames, polarisation):

        cR=0
        cI=0
        if type(fact)!=type(1.0j):
          cR=fact
        else:
          cR=fact.real
          cI=fact.imag
        if( cR):
          t_dipR.Read_add(repe+opname , cR,0)
        if( cI):
          if(t_dipI is None):
            t_dipI=Hilbertxx.Hxx_TransitionMatrix()
          t_dipI.Read_add(  repe+opname   , cI,0)

      sR=SfromT(t_dipR)
      sI=SfromT(t_dipI)
 
      operatorS.append(TrialMatrixC(sR,sI)  )
      






      
    else:
      
      for opname  in self.quadrunames:
        t_dip = Hilbertxx.Hxx_TransitionMatrix() 
        t_dip.Read_add( repe+opname, 1.0, 0)
        ars = t_dip.create_3Array_and_Clean()
        if( usecomplex):
          operatorS.append(TrialMatrixC(Sparsa.Sparsa3A( ars[2], ars[0], ars[1] ), None)  )
        else:
          operatorS.append(Sparsa.Sparsa3A( ars[2], ars[0], ars[1] ) )
        

    allX=Numeric.array([-100000.0, self.El2l3-0.2-self.shift*0 , self.El2l3+0.2-self.shift*0, 100000.0])
    allY=Numeric.array([(self.all1),(self.all1),(self.all2),(self.all2) ])



    #####
    print " OPERATORS ", operatorS
    if(not usepython ):
     
      result = spectrumFunction.spectrumFunctionC( operatorS, S_base, S_ecci, self.erange,self.temp,
                                                   self.dxleft,self.dxright,self.npunti,  
                                                   allX, allY ,self.shift_invert ,self.nsearchedeigen  )
    else:
      if usecomplex:
        S_base = TrialMatrixC( S_base, S_baseI)
        S_ecci = TrialMatrixC( S_ecci, S_ecciI)
        
      
      evals, evects = solveEigenSystem( S_base , self.nsearchedeigen)


      # evals1, evects1 = solveEigenSystem( S_ecci , self.nsearchedeigen)



      data4spect=faispettri( evals, evects   , S_ecci ,self.erange , operatorS, min(225, S_ecci.dim))
      print data4spect
      print dir(S_ecci)

      factors = getFactors( data4spect["enes"] , self.temp , self.erange ) 

      result=[0]*(len(operatorS)+1)
      energies=Numeric.arrayrange( self.dxleft ,  self.dxright, (-self.dxleft + self.dxright+0.0)/ self.npunti  )
      result[-1]=energies
      energies = energies + self.all1*1.0j * Numeric.less(self.El2l3 , energies ) + self.all2*1.0j * Numeric.less(energies ,  self.El2l3 )
      
      for pol in range(len(operatorS) ):
        result[pol]=0.0
        for i in range(len(factors)):
          result[pol]=result[pol] +  factors[i]*spectra(data4spect["alphas"][i][pol],
                                                        energies+evals[i])*data4spect["norms"][i][pol]

    print " spec calcolato "
    result[len(operatorS)][:]=result[len(operatorS)]+self.shift
    result = [result[len(operatorS)],]+result[:len(operatorS) ]

    evals1, evects1 = solveEigenSystem( S_ecci , self.nsearchedeigen)
    print evals
    print evals1 

    return result






  def spectrumRIXS(self, repe, coeff,  names,coeff_exc, names_exc,coeff_fin, names_fin,  polarisationIn=None,polarisationOut=None,
               ein=None, eout1=None, eout2=None, dout=None, gammain=None, gammaout=None, HFe = 0.0, HFf=1.0):


    raise " OK " 
    t_base = Hilbertxx.Hxx_TransitionMatrix()
    t_baseI = Hilbertxx.Hxx_TransitionMatrix()
    
    self.compose( t_base, t_baseI,  repe,  coeff, names , "base" )
        
    t_excited = Hilbertxx.Hxx_TransitionMatrix()
    t_excitedI = Hilbertxx.Hxx_TransitionMatrix()
    
    self.compose( t_excited, t_excitedI,  repe, coeff_exc, names_exc ,"exci", hoppfact=HFe) 

    t_final = Hilbertxx.Hxx_TransitionMatrix()
    t_finalI = Hilbertxx.Hxx_TransitionMatrix()
    
    self.compose( t_final, t_finalI,  repe, coeff_fin, names_fin ,"final", hoppfact=HFf) 


    ars_base=t_base.create_3Array_and_Clean()
    ars_ecci=t_excited.create_3Array_and_Clean()
    ars_fin =t_final  .create_3Array_and_Clean()





    S_base=Sparsa.Sparsa3A( ars_base[2], ars_base[0], ars_base[1] )
    S_ecci=Sparsa.Sparsa3A( ars_ecci[2], ars_ecci[0], ars_ecci[1] )
    S_fin =Sparsa.Sparsa3A( ars_fin [2], ars_fin [0], ars_fin [1] )

    S_ecciT =Sparsa.Sparsa3A( ars_ecci [2], ars_ecci [1], ars_ecci [0] )

    print " S_fin dim " , S_fin.dim

    ars_base =t_baseI   .create_3Array_and_Clean()
    ars_ecci =t_excitedI.create_3Array_and_Clean()
    ars_fin  =t_finalI  .create_3Array_and_Clean()


    S_baseI = Sparsa.Sparsa3A( ars_base[2], ars_base[0], ars_base[1] )
    S_ecciI = Sparsa.Sparsa3A( ars_ecci[2], ars_ecci[0], ars_ecci[1] )
    S_finI  = Sparsa.Sparsa3A( ars_fin[2], ars_fin[0], ars_fin[1] )


    S_ecciIT  = Sparsa.Sparsa3A( ars_ecci[2], ars_ecci[1], ars_ecci[0] )

    print " S_finI dim " , S_finI.dim

    
    S_base = TrialMatrixC( S_base, S_baseI)
    S_ecci = TrialMatrixC( S_ecci, S_ecciI)
    S_ecciT  = TrialMatrixC( S_ecciT, S_ecciIT)
    S_fin  = TrialMatrixC( S_fin, S_finI)
        
    print S_base.dim
    print S_ecci.dim
    print S_fin.dim
    evals, evects = solveEigenSystem( S_base , self.nsearchedeigen)

    ein=ein+evals[0]


    qsomma=evects[0]
    fv=open("vecto","w")
    for i in range(qsomma.len()):
      print i, "  ", qsomma.vr.get_value(i)
      fv.write("%d %e %e \n"%(i,qsomma.vr.get_value(i),qsomma .vi.get_value(i)))




    # raise " stampato veto  OK %d " % S_base.dim 
    




##     evals, evects = solveEigenSystem( S_fin , 40)

##     qsomma=evects[1]
##     fv=open("vecto","w")
##     for i in range(qsomma.len()):
##       print i, "  ", qsomma.vr.get_value(i)
##       fv.write("%d %e %e \n"%(i,qsomma.vr.get_value(i),qsomma .vi.get_value(i)))




##     raise " OK " 
    


    operatorS=[]

    t_dip  = Hilbertxx.Hxx_TransitionMatrix() 
    t_dipI = Hilbertxx.Hxx_TransitionMatrix() 
    for opname, fact  in zip(self.quadrunames, polarisationIn):
      fact=complex(fact)
      if fact.real !=0.0:
        t_dip.Read_add( repe+opname, fact.real, 0)
      if fact.imag !=0.0:
        t_dipI.Read_add( repe+opname, fact.imag, 0)

      
    ars  = t_dip .create_3Array_and_Clean()
    arsI = t_dipI.create_3Array_and_Clean()

    operatorS.append(TrialMatrixC(Sparsa.Sparsa3A( ars[2], ars[0], ars[1] ),Sparsa.Sparsa3A( arsI[2], arsI[0], arsI[1] ) )  )
    
    class4sparse = S_ecci .__class__
    class4vect   = class4sparse.getClass4Vect()

    dumA = class4vect(S_ecci.dim)
    dumA.set_to_zero()
    operatorS[0].Moltiplica(dumA, evects[0]  )
    norma = (dumA.scalare(dumA)).real

    dumA.normalizzaauto() 

    qsomma = class4vect(S_ecci.dim)
    qsomma.set_to_zero()

    diff  = checkInv( S_ecci, ein, gammain, qsomma, dumA )
    normdiff = class4vect.scalare(diff,diff).real
    print " il modulo della differenza est ", normdiff
    count=0
    print class4vect.scalare(diff,diff).real

    oldnormf=1.0
    normf   =1.0

    oldnormdiff=1.0
    


    bestsolution  = class4vect(S_ecci.dim)
    bestvalue = 1.0e40

    count=0
    
    
    # abbiamo qui finalmente che qsomma moltiplicato per sqrt(norma)
    # e il nostro vettore che ci permette di calcolare lo spettro

    # qsomma=0
    

    S_ecci .shift=  - ein -0*gammain

    
    qsomma=CGsolve(S_ecci, gammain*gammain,  dumA, tol=1.0e-10, nmax = 1000, restart=None)
    bestsolution=CGsolve(S_ecci, gammain*gammain,  dumA, tol=1.0e-20, nmax = 1000, restart=qsomma)
    
    S_ecci .shift=  - ein +  1.0j*gammain
    qsomma.set_to_zero()
    S_ecci.Moltiplica(qsomma, bestsolution)

    print  (qsomma.scalare(qsomma)).real

    raise " OK " 


    fv=open("vectoA","w")
    for i in range(qsomma.len()):
      print i, "  ", qsomma.vr.get_value(i)
      fv.write("%d %e %e \n"%(i,qsomma.vr.get_value(i),qsomma .vi.get_value(i)))


    operatorS=[]
    t_dip  = Hilbertxx.Hxx_TransitionMatrix() 
    t_dipI = Hilbertxx.Hxx_TransitionMatrix() 
    for opname, fact  in zip(self.quadrunames[10:], polarisationOut):
      fact=complex(fact)
      if fact.real !=0.0:
        t_dip.Read_add( repe+opname, fact.real, 0)
      if fact.imag !=0.0:
        t_dipI.Read_add( repe+opname, fact.imag, 0)
    ars  = t_dip .create_3Array_and_Clean()
    arsI = t_dipI.create_3Array_and_Clean()
    ##  da notare la trasposizione !!!!!
    operatorS.append(TrialMatrixC(Sparsa.Sparsa3A( ars[2], ars[0], ars[1] ),Sparsa.Sparsa3A(arsI[2], arsI[0], arsI[1] ) )  )




    if(1):


      dumC = class4vect(S_fin.dim)
      dumC.set_to_zero()


      operatorS[0].Moltiplica(dumC, qsomma)
      qsomma = dumC
      fv=open("vectoD","w")
      for i in range(qsomma.len()):
        print i, "  ", qsomma.vr.get_value(i)
        fv.write("%d %e %e \n"%(i,qsomma.vr.get_value(i),qsomma .vi.get_value(i)))
      raise "  OK " 

      
      data4spect=faispettri( Numeric.array([0.0]), [qsomma]   , S_fin ,0.00001 , operatorS, 200)
      normb = data4spect["norms"][0][0]
      
      
      result = [None, None]
      
      energies=Numeric.arrayrange( eout1 ,  eout2, dout   )
      
      result[0]=energies
      energies = energies + estraigamma(gammaout, energies) *1.0j  
      
      result[1]= spectra(data4spect["alphas"][0][0],energies)*normb*norma
      
 
    return result



  def spectrumRIXS(self, repe, coeff,  names,coeff_exc, names_exc,coeff_fin, names_fin,  polarisationIn=None,polarisationOut=None,
               ein=None, eout1=None, eout2=None, dout=None, gammain=None, gammaout=None, HFe = 0.0, HFf=1.0):


    t_base = Hilbertxx.Hxx_TransitionMatrix()
    t_baseI = Hilbertxx.Hxx_TransitionMatrix()
    
    self.compose( t_base, t_baseI,  repe,  coeff, names , "base" )
        
    t_excited = Hilbertxx.Hxx_TransitionMatrix()
    t_excitedI = Hilbertxx.Hxx_TransitionMatrix()
    
    self.compose( t_excited, t_excitedI,  repe, coeff_exc, names_exc ,"exci", hoppfact=HFe) 

    t_final = Hilbertxx.Hxx_TransitionMatrix()
    t_finalI = Hilbertxx.Hxx_TransitionMatrix()
    
    self.compose( t_final, t_finalI,  repe, coeff_fin, names_fin ,"final", hoppfact=HFf) 


    ars_base=t_base.create_3Array_and_Clean()
    ars_ecci=t_excited.create_3Array_and_Clean()
    ars_fin =t_final  .create_3Array_and_Clean()


    S_base=Sparsa.Sparsa3A( ars_base[2], ars_base[0], ars_base[1] )
    S_ecci=Sparsa.Sparsa3A( ars_ecci[2], ars_ecci[0], ars_ecci[1] )
    S_fin =Sparsa.Sparsa3A( ars_fin [2], ars_fin [0], ars_fin [1] )

    S_ecciT =Sparsa.Sparsa3A( ars_ecci [2], ars_ecci [1], ars_ecci [0] )

    print " S_fin dim " , S_fin.dim

    ars_base =t_baseI   .create_3Array_and_Clean()
    ars_ecci =t_excitedI.create_3Array_and_Clean()
    ars_fin  =t_finalI  .create_3Array_and_Clean()


    S_baseI = Sparsa.Sparsa3A( ars_base[2], ars_base[0], ars_base[1] )
    S_ecciI = Sparsa.Sparsa3A( ars_ecci[2], ars_ecci[0], ars_ecci[1] )
    S_finI  = Sparsa.Sparsa3A( ars_fin[2], ars_fin[0], ars_fin[1] )


    S_ecciIT  = Sparsa.Sparsa3A( ars_ecci[2], ars_ecci[1], ars_ecci[0] )

    print " S_finI dim " , S_finI.dim

    
    S_base = TrialMatrixC( S_base, S_baseI)
    S_ecci = TrialMatrixC( S_ecci, S_ecciI)
    S_ecciT  = TrialMatrixC( S_ecciT, S_ecciIT)
    S_fin  = TrialMatrixC( S_fin, S_finI)
        
    print S_base.dim
    print S_ecci.dim
    print S_fin.dim
    evals, evects = solveEigenSystem( S_base , self.nsearchedeigen)

    ein=ein+evals[0]


    operatorS=[]

    t_dip  = Hilbertxx.Hxx_TransitionMatrix() 
    t_dipI = Hilbertxx.Hxx_TransitionMatrix() 
    for opname, fact  in zip(self.quadrunames, polarisationIn):
      fact=complex(fact)
      if fact.real !=0.0:
        t_dip.Read_add( repe+opname, fact.real, 0)
      if fact.imag !=0.0:
        t_dipI.Read_add( repe+opname, fact.imag, 0)

      
    ars  = t_dip .create_3Array_and_Clean()
    arsI = t_dipI.create_3Array_and_Clean()

    operatorS.append(TrialMatrixC(Sparsa.Sparsa3A( ars[2], ars[0], ars[1] ),Sparsa.Sparsa3A( arsI[2], arsI[0], arsI[1] ) )  )
    
    class4sparse = S_ecci .__class__
    class4vect   = class4sparse.getClass4Vect()

    dumA = class4vect(S_ecci.dim)
    dumA.set_to_zero()
    operatorS[0].Moltiplica(dumA, evects[0]  )
    norma = (dumA.scalare(dumA)).real

    dumA.normalizzaauto() 

    
    S_ecci .shift=  - ein -1.0j*gammain
    S_ecciT.shift=  - ein +1.0j*gammain


    S_ecci .shift=  - ein -0*gammain

    qsomma=CGsolve(S_ecci, gammain*gammain,  dumA, tol=1.0e-10, nmax = 1000, restart=None)
    bestsolution=CGsolve(S_ecci, gammain*gammain,  dumA, tol=1.0e-20, nmax = 1000, restart=qsomma)

    
    S_ecci .shift=  - ein +  1.0j*gammain
    qsomma.set_to_zero()

    S_ecci.Moltiplica(qsomma, bestsolution)


    operatorS=[]
    t_dip  = Hilbertxx.Hxx_TransitionMatrix() 
    t_dipI = Hilbertxx.Hxx_TransitionMatrix() 
    for opname, fact  in zip(self.quadrunames[10:], polarisationOut):
      fact=complex(fact)
      if fact.real !=0.0:
        t_dip.Read_add( repe+opname, fact.real, 0)
      if fact.imag !=0.0:
        t_dipI.Read_add( repe+opname, fact.imag, 0)
    ars  = t_dip .create_3Array_and_Clean()
    arsI = t_dipI.create_3Array_and_Clean()
    operatorS.append(TrialMatrixC(Sparsa.Sparsa3A( ars[2], ars[0], ars[1] ),Sparsa.Sparsa3A(arsI[2], arsI[0], arsI[1] ) )  )

    if(1):


##       dumC = class4vect(S_fin.dim)
##       dumD = class4vect(S_fin.dim)
##       dumC.set_to_zero()
##       operatorS[0].Moltiplica(dumC, qsomma)
##       fv=open("vectoD","w")
##       for i in range(dumC.len()):
##         print i, "  ", dumC.vr.get_value(i)
##         fv.write("%d %e %e \n"%(i,dumC.vr.get_value(i),dumC .vi.get_value(i)))
##       dumC.set_to_zero()
##       dumC.set_value(457,5.0)
##       for i in [462,468,492,576,876]:
##         dumC.set_value(i,-1)
##       dumD.set_to_zero()
##       S_fin.Moltiplica(dumD,dumC)
##       dumD.normalizzaauto()
##       dumC.normalizzaauto()
##       print self.GetESLcounters(3,dumD)
##       raise " OK " 
##       data4spect=faispettri( Numeric.array([0.0]), [dumC]   , S_fin ,0.00001 , [None], 300)

      data4spect=faispettri( Numeric.array([0.0]), [qsomma]   , S_fin ,0.00001 , operatorS, 60)
      #data4spect=faispettri( Numeric.array([0.0]), [dumC]   , S_fin ,0.00001 , [None], 60)
      normb = data4spect["norms"][0][0]
      
      
      result = [None, None]
      
      energies=Numeric.arrayrange( eout1 ,  eout2, dout   )
      
      result[0]=energies
      energies = energies + estraigamma(gammaout, energies) *1.0j  
      
      result[1]= spectra(data4spect["alphas"][0][0],energies)*normb*norma
      
 
    return result





def checkInv( S_ecci, ein, gammain, qsomma, evect):
      class4sparse = S_ecci .__class__
      class4vect   = class4sparse.getClass4Vect()

      dum1 = class4vect( evect.len() )
      dum1.set_to_zero()

      S_ecci.Moltiplica(dum1, qsomma )

      dum1.add_from_vect_with_fact( qsomma  , -ein-1.0j*gammain   )
      dum1.add_from_vect_with_fact( evect  , -1.0   )

      return dum1


def sommaqs( qs      , fattori  ):
  class4sparse = TrialMatrixC
  class4vect   = class4sparse.getClass4Vect()
  
  qsomma=class4vect(   qs[0].len() )
  qsomma.set_to_zero()
  for i in range( len(fattori)):
    qsomma.add_from_vect_with_fact( qs[i]  ,  fattori[i]  )
  return qsomma 
 
def inverfirstcolumn( alphas, betas):
    N=len(alphas)
    terms = Numeric.zeros(len(alphas),"d")
    terms[N-1]=1.0/alphas[N-1]
    for i in range(N-2,-1,-1):
        terms[i]=1.0/(alphas[i]-betas[i]*betas[i]*terms[i+1])


    column = Numeric. zeros(len(alphas),"d")

    column[0]=terms[0]
    for i in range(1,N):
        column[i]= -column[i-1]*betas[i-1]*terms[i]

    return column



##struct Array {
##    Array( int size ) ;
##    Array( int fd, int sd ) ;

##    ~Array() ;
##    Array * get_array(int index);



##    double get_value( int index ) const;
##    void set_to_zero(  )        ;
##    void set_to_one(  ) ;
##    void set_value( int index, double value );
##    void set_all_random(  double value );

##    Array* copy() const;

##    Array& operator+=( const Array& a ) ;
##    int len();

##     static void copy_to_a_from_b(Array &a, Array&b);
##     static void mat_mult(Array & b, double * ArrayFLOAT2dNOCHECK ,Array & a);


##     void add_from_vect( Array&vect);
##     void add_from_vect_with_fact( Array&vect, double fact);



##  double *dataAddress() ;
##  void  dumptofile(char * nome);
##  static double scalare(Array &a, Array &b) ;
##  static double sqrtscalare(Array &a, Array &b);
##  void normalizza(double norm);
##  void normalizzaauto();
##  void mult_by_fact(double fact);


##};                                                                                                                                  



#       factors = getFactors( data4spect["enes"] , self.temp , self.erange ) 

#      result=[0]*(len(operatorS)+1)
#      energies=Numeric.arrayrange( self.dxleft ,  self.dxright, (-self.dxleft + self.dxright+0.0)/ self.npunti  )
#      result[-1]=energies
#      energies = energies + self.all1*1.0j * Numeric.less(self.El2l3 , energies ) + self.all2*1.0j * Numeric.less(energies ,  self.El2l3 )
      
#      for pol in range(len(operatorS) ):
#        result[pol]=0.0
#        for i in range(len(factors)):
#          result[pol]=result[pol] +  factors[i]*spectra(data4spect["alphas"][i][pol],
#                                                        energies+evals[i])*data4spect["norms"][i][pol]

#    print " spec calcolato "
#    result[len(operatorS)][:]=result[len(operatorS)]+self.shift
#    result = [result[len(operatorS)],]+result[:len(operatorS) ]


#    return result
     


def    Riduci(coeffs ,  names, calcnames, calcpars ):
  for i in range(len(names)):
    name=names[i]
    reduction = searchReduction(name, calcnames, calcpars   )
    if( reduction!=1.0):
      coeffs[i]=coeffs[i]*reduction
    
def searchReduction(name, calcnames, calcpars ):
  import re
  if(name[-2:]=="F0"):
    return 1
  tags = { "couche\d_\d_[FG]":[6,9], "couche\d_[FG]":[6,7],"MixCouche\d_\d_K":[9,12]}
  namepar = None
  for tag in tags.keys():
    searchobject = re.search(tag,name)
    if searchobject is not None:
      range1,range2= tags[tag]
      postfix= searchobject.group(0)[range1:range2]
      namepar = "reduc_"+postfix
      break
  print " per ", name ," la riduzione est data dal parametro ", namepar
  if namepar is not None:
    try:
      pos = calcnames.index(namepar)
    except ValueError:
      pos=-1
    if pos != -1:
      print "       che e nella lista  = ", calcpars[pos]
      return calcpars[pos]
    else:
      print "       che non e nella lista  = "
      return 1
  return 1

  

  
def getFactors(enes, temperature, tole):
    ne = len(enes.keys())
    E = Numeric.zeros(ne,"d")
    for i in range(ne):   E[i] = enes[i]
    factors = Numeric.exp(-(E-E[0])/temperature)
    factors = factors/Numeric.sum(factors)
    Sum=0
    imax=ne
    for i in range(ne):
        Sum=Sum+factors[i]
        if( factors[i] < tole ):
            imax=i
            break
    factors=factors[:imax]/Sum
    return factors
