#!/scisoft/ESRF_sw/linux_i386_03/bin/python
import sys
print sys.argv
batch=0
if not batch:
    import qt
    import qwt
    import iqt
    from qwt.qplt import Plot, Curve, Pen, Red, Green , Yellow, Blue
def write(A, filename):
    f=open(filename,"w")
    for i in range(0,len(A)):
        A[i]=A[i].tolist()
    
    for Al in apply(zip,A):
        for tok in Al:
            if type(tok)==type(1.0):
                f.write("%e  " % tok)
            else:
                f.write( "  %e  %e  " % (tok.real, tok.imag))
            f.write(" ")
        f.write("\n")
    f.close()




def help():
    s="""
 AVAILABLE HXX-COMMANDS ( examples of )
set(system)    # or use directly python syntax to change the attributes of object  system.
               # When using the set(system) function you can also give a filename of a previously
               # saved ( and edited ) file

save(system)   #  You are asked a filename and  you have to quote it like "filename"

scriviFiles(nel3d, prefix="", nhopped=n, spinfixed=0  )  # To generate matrices files. This command
               # creates the Hilbert space and writes the terms composing the Hamiltonian into (sub)directory "datas" which must exist.
               # The above example generates a Hilbert space where nel3d electrons are distributed in the 3d
               # shell in all possible combinations.
               # Giving nhopped=n as argument, besides the configurations with nel3d electrons on 3d also the configurations
               # with 5+i electrons on 3d and 10-i on Ligands orbitals are considered ( 0<=i<=nhopped )
               # If the spinfixed argument is set to 1, the spanned  hilbert space vectors have a total Sz=nel3d/2 for nel3d<=5
               # and Sz = (10-nel3d)/2 for nel3d>5.
            
system.case="datas" # This chooses the system that you have expanded with scriviFiles command.
               #    You can acces this property also with the set(system)    command
  
res=system.GetSpectrum() #   to calculate a spectrum.  res is a list.res[0] is an array of floats : the energies.
               # res[i>=1] are the resonances ( real and imaginary part)
               # for the quadrupole 1s->3d transition for   Mz=-2 up to Mz=2 ( i from 1 up to 5)
               # and the Mixeddipole 1s->LiganD transition for  Mz=-2 up to Mz=2   ( i from 6 up to 10 ) 
               # The mixed dipole contribution is calculated hybridizing a virtual 4p orbital with the oxygen ones.

res=system.GetSpectrum([fm2,  fm1,f0,f1, fp2 ]) # Absorption for a  definite polarisation, where fm2,  fm1,f0,f1, fp2
               # are five coefficients of the polarisatin in  the angular momentum representation.
               # For xy
               # f's = [ 1,0,0,0,-1  ]
               # to calculate instead x2y2
               # f = [ 1,0,0,0,1 , 0   ]
    
write(res,filename)   #  To save res on a file
               #  and the file will have several columns ( 1 for energies and two (real, imag) for each polarisation)

Plot(Curve(res[0], res[n>=1].imag, Pen(Red), "legend")) # To plot. Change n to select the polarisation that you want
               #  The plotting feature uses qt, iqt and qwt. If you run a long job in the background you must
               # desactivate this graphics feature by commenting out, in the file initGDIS.py, the lines where such packages
               # are imported
   
Es, S2s, L2s, SL2, occPs , Szs= system.GetESLcounters()  #  to get statistics: a list of E, the list of 3d  L2 ,
               # the list of 3d S2, the list of 3d SL, a list of ligand expected occupation, a list of 3d Sz expectation
  

re= system.GetRIXS( polarisationIn=[0.0, 1.0,0.0, 1.0,.0] ,polarisationOut=[1.0,0.0], ein=ene,
                        eout1=eout1, eout2=eout2, dout=dout, gammain=0.4, gammaout=  [ 0.5 ,  Echange , 1.0  ]) #
               # calculate the rixs spectra, for an incoming photon at energy ein and quadrupolar polarisation=polarisationIn
               # The outgoing dipolar polarisation is polarisationOut.
               # The spectra of the outgoing radiation is calculated between eout1 and eout2 with an energy step=dout.
               # Gammain is the lorentzian broadening for the absorption of the incoming photon and
               # gammaout for the outgoing one.
               # Gammaout can be either a number or a list. In case of a list the first entry is used for spectra energy
               # below tEchange, the third entry is used above.
         """
    print s

from System import System


def set(system):
    names = system.basenames+system.excinames+system.finnames+system.calcnames
    stop={}
    lungs = [len(system.basenames),len(system.excinames),len(system.finnames),len(system.calcnames)]
    stop[1]="---- BASE HAMILTONIAN ---"
    stop[1+lungs[0]]="---- EXCITED HAMILTONIAN ---"
    stop[1+lungs[0]+lungs[1]]="---- FINAL HAMILTONIAN  ---"
    stop[1+lungs[0]+lungs[1]+lungs[2]]="---- CALCULATION PARAMETERS ---"
    while(1):
        count = 1;
        for name in names:
            if(count in stop.keys()):
                print stop[count]
            num=" "*(3-len(str(count)))+str(count)
            print num+") "+name +" "*(15-len(name)),":    ",
            exec("val=system."+name)
            print getattr(system,name)
            count=count+1
        print "--------------------------------------------------------------------"
        print " select a value to change , 0 to stop, filename to read values"
        user = input(" ")
        if( type(user)==type(1)):
            if( user==0): break
            if(user>0 and user<=len(names)):
                print " insert new value for ", names[user-1]
                newval=input(" ")
                setattr(system,names[user-1], newval)
        elif( type(user)==type("s")):
            s=open(user,"r").read()
            exec(s)
            
def nonint_read(system, filename):
    s=open(filename,"r").read()
    exec(s)
    

def save(system):
    print "--------------------------------------------------------------------"
    print " insert filename to write values"
    user = input(" ")
    nonint_save(system, user)
    
def nonint_save(system, filename):
    names =   system.basenames+system.excinames+system.finnames  +  system.calcnames
    f=open(filename,"w")
    for n in names:
        val = str(getattr(system,n))
        if(type(getattr(system,n))==type("s") ):
            val="\""+val+"\""
        f.write("system."+n+"="+val+"\n")


from matrixcreator import scriviFiles

selectall=1


import parse
bs=parse.Atom("out36_base")
ex=parse.Atom("out36_exci")
fin=parse.Atom("out36_fin")


basekvs= [
    "base_couche1_F0"     , 5.0  , 
    "base_couche1_F2"     , bs.FDict["(3d,3d)2"]  , 
    "base_couche1_F4"     , bs.FDict["(3d,3d)4"]  ,  

    "base_couche0_1_F0"   , 5.5  ,
    "base_couche0_1_G2"   , bs.GDict["(1s,3d)2"]   ,



    "base_couche1_3_F0"   , 5.5  ,
    "base_couche1_3_G1"   , bs.GDict["(3p,3d)1"]   ,
    "base_couche1_3_G3"   , bs.GDict["(3p,3d)3"]   ,


    "base_SO_1"           , bs.zetaDict["3d"]      ,
    "base_Sop_Zero"       , 0.00001  ,
    "base_Sop_Minus"       , 0.0000  ,
    "base_Sop_Plus"       , 0.0000  ,
    "base_counterDL"       , -4,
    ]
excikvs= [
    "exci_couche1_F0"     ,   5.0, 
    "exci_couche1_F2"     , ex.FDict["(3d,3d)2"]  , 
    "exci_couche1_F4"     , ex.FDict["(3d,3d)4"]  ,  
    "exci_couche0_1_F0"   , 5.5  ,
    "exci_couche0_1_G2"   , ex.GDict["(1s,3d)2"]   ,

    "exci_couche1_3_F0"   , 5.5  ,
    "exci_couche1_3_G1"   , ex.GDict["(3p,3d)1"]   ,
    "exci_couche1_3_G3"   , ex.GDict["(3p,3d)3"]   ,


    "exci_SO_1"           , ex.zetaDict["3d"]      ,


    
    "exci_Sop_Zero"       , 0.00001  ,
    "exci_Sop_Minus"       , 0.0000  ,
    "exci_Sop_Plus"       , 0.0000  ,
    "exci_counterDL"       , -4,
   ]
finalkvs= [
    "final_couche1_F0"     ,   5.0, 
    "final_couche1_F2"     , fin.FDict["(3d,3d)2"]  , 
    "final_couche1_F4"     , fin.FDict["(3d,3d)4"]  ,  
    "final_couche0_1_F0"   , 5.5  ,
    "final_couche0_1_G2"   , fin.GDict["(1s,3d)2"]   ,

    "final_couche1_3_F0"   , 5.5  ,
    "final_couche1_3_F2"   , fin.FDict["(3p,3d)2"]   ,
    "final_couche1_3_G1"   , fin.GDict["(3p,3d)1"]   ,
    "final_couche1_3_G3"   , fin.GDict["(3p,3d)3"]   ,



    "final_SO_3"           , fin.zetaDict["3p"]      ,
    
    "final_SO_1"           , fin.zetaDict["3d"]      ,

    
    "final_Sop_Zero"       , 0.00001  ,
    "final_Sop_Minus"       , 0.0000  ,
    "final_Sop_Plus"       , 0.0000  ,
    "final_counterDL"       , -4,
   ]



quadrunames = [
"quadrupole_couches_0_1_Km_-2",
"quadrupole_couches_0_1_Km_-1",
"quadrupole_couches_0_1_Km_0",
"quadrupole_couches_0_1_Km_1",
"quadrupole_couches_0_1_Km_2",

"mixedDipole_couches_0_2_mz_-2",
"mixedDipole_couches_0_2_mz_-1",
"mixedDipole_couches_0_2_mz_0",
"mixedDipole_couches_0_2_mz_1",
"mixedDipole_couches_0_2_mz_2",

"dipole_couches_0_3_Km_-1",
"dipole_couches_0_3_Km_0",
"dipole_couches_0_3_Km_1",

]

calculkvs = [
    "case" ,       "datas",
    "reduc_1",      0.8 ,
    "reduc_0_1",    0.8 ,
    "all1",         0.1  ,
    "El2l3",        700,
    "all2",         0.1 ,
    "shift",        0 ,
    "npunti",       500 , 
    "dxleft",       -0.1 ,
    "dxright",       0.1 , 
    "temp",          0.009 ,
    "erange",        0.1 ,
    "shift_invert",  0,
    "nsearchedeigen",10,
    "Vs", 2.0,
    "Vp", 1.0,
    "VC0",0.0,
    "VC1",0.0,
    "Dips", 1.0,
    "Dipp", 1.0,
    "DREF", 1.0,
    "ALPHAVC", -3.0,
    "ALPHAVSP", -3.0,
    "ALPHADIPO", -2.0,
    "BONDS", [[-1.,0,0], [1.,0,0],[0,-1.,0],[0,1.,0],[0,0,-1.],[0,0,1.]],
    "hoppfactexci",1.0,
    "hoppfactfinal", 1.0,
    "facts_hop",None
    ]


angularnames = [
    "base_Sop_Zero"   ,
    "base_Sop_Minus"  ,
    "base_Sop_Plus"   ,
    "base_Lop_Zero"   ,
    "base_Lop_Minus"  ,
    "base_Lop_Plus" 
    ]

counternames = [
    "base_counterDL",
    "base_Sop_Zero"
    ]

system = System( basekvs[0::2], basekvs[1::2] , excikvs[0::2], excikvs[1::2], finalkvs[0::2], finalkvs[1::2],
                 calculkvs[0::2], calculkvs[1::2] ,   quadrunames, angularnames, counternames)

print """type
help()
to get help
"""
def calcola():
    exec(open("paraVC").read())
    system.case="dataf_m"
    resm=system.GetSpectrum()
    resma = system.GetSpectrum(  [ 1.0 , 0.0 , 0 , 0 ,  -1.0 ] )
    resmb = system.GetSpectrum(  [ 0.0 , 1.0 , 0 , 1.0 , 0 ] )
    system.case="dataf_p"
    resp=system.GetSpectrum()
    respa = system.GetSpectrum(  [ 1.0 , 0.0 , 0 , 0 ,  -1.0 ] )
    respb = system.GetSpectrum( [ 0.0 , 1.0 , 0 , 1.0 , 0 ] )

    system.case="datas"
    ress=system.GetSpectrum([ 0.0 , 1.0 , 0 , 1.0 , 0 ])
    
    dm=(resm[1]+resm[2]+resm[3]).imag
    dp=(resp[1]+resp[2]+resp[3]).imag
    return resm[0], dm, resp[0],dp,resma[1].imag,resmb[1].imag,respa[1].imag,respb[1].imag, ress[1].imag





def calcola2():
    exec(open("paraVC").read())
    system.case="dataf_m2"
    resm=system.GetSpectrum()
    system.case="dataf_p2"
    resp=system.GetSpectrum()
    dm=(resm[1]+resm[2]+resm[3]).imag
    dp=(resp[1]+resp[2]+resp[3]).imag
    return resm[0], dm, resp[0],dp

def calcolaS():
    exec(open("paraVC").read())
    system.case="datas"
    res=system.GetSpectrum()
    d=(res[1]+res[2]+res[3]).imag
    return res[0], d


def S2():
    system.case="datas"
    exec(open("paraVC").read())
    Es, S2s, L2s, SL2, occPs , Szs= system.GetESLcounters()
    print S2s
    print occPs
    print Szs
    


def plotta(d):
    res=Plot(Curve(d[0], d[1], Pen(), "legend"), Curve(d[2], d[3], Pen(Red), "legend" ))
    resb=Plot(Curve(d[0], d[4], Pen(Red), "legend"), Curve(d[2], d[5], Pen(Red), "legend" ),
              Curve(d[0], d[6], Pen(   ), "legend"), Curve(d[0], d[7], Pen(   ), "legend" ))
    return res, resb

def plottaab(d,a,b,shift,fact):
    cs=[Curve(d[0], d[1], Pen(), "legend"), Curve(d[2], d[3], Pen(Red), "legend" )]
##    cs.append(      Curve(a[:,0]-shift, fact*a[:,1], Pen(Red), "")  )                                              
##    cs.append(      Curve(b[:,0]-shift, fact*b[:,1], Pen(), "")                                               )
    res=Plot(*cs)
    resb=Plot(Curve(d[0], d[4], Pen(Red), "legend"), Curve(d[2], d[5], Pen(Red), "legend" ),
              Curve(d[0], d[6], Pen(   ), "legend"), Curve(d[0], d[7], Pen(   ), "legend" ), Curve(d[0], d[8], Pen(   ), "legend" )   )
    return res, resb

def plottaS(d):
    res=Plot(Curve(d[0], d[1], Pen(), "legend"))
    return res

def plottaSl(dl):
    a=[]
    for d in dl:
        a.append(Curve(d[0], d[1], Pen(), "legend"))
    res=Plot(*a)
    return res

from Numeric import *


def leggifile(nome):
    s=open(nome,"r").read()
    sl=string.split(s,"\n")
    res=[]
    for l in sl:
        l=string.split(l)
        if(len(l)!=2):
            break
        res.append(map(string.atof,l))
    return array(res)

    


def fit8():
    import os
    os.system("cp par paraVC")
    d=calcola()
    a=0
    b=0
    return plottaab(d,a,b,6515.5,0.2e-6)

import os
if batch:
    import os
    os.system("cp para8 paraVC")
    exec(open("paraVC").read())
    system.case="datas"
    
    red={}
    import sys

    Eins=arrayrange(24,28+7.5,0.15)

    
    ene = Eins[string.atoi(sys.argv[1] )]

    

    
    shift=ene-Eins[0]
    
    Es, S2s, L2s, SL2, occPs , Szs= system.GetESLcounters()

    E0=Es[0]
    ene=ene+Es[0]

    eout1=630+shift*0
    eout2=790.001+shift*0
    dout=0.15

    
    re= system.GetRIXS( polarisationIn=[1.0, 1.0,0.0, 1.0,.0] ,polarisationOut=[1.0,0.0], ein=ene,
                        eout1=eout1, eout2=eout2, dout=dout, gammain=0.4, gammaout=  [ 0.5 ,  eout1+25 , 1.0  ])
    red[ene]=re
##     re1= system.GetRIXS( polarisationIn=[1.0, 0 ,  0.0 , .0  ,  1.0] ,polarisationOut=[1.0,0.0], ein=ene,
##                         eout1=eout1, eout2=eout2, dout=dout, gammain=0.5, gammaout=  [ 1.0 ,  660.0 + 57.01   ,  4.0  ])

##     re2= system.GetRIXS( polarisationIn=[1.0, 0 ,  0.0 , .0  ,  -1.0] ,polarisationOut=[1.0,0.0], ein=ene,
##                         eout1=eout1, eout2=eout2, dout=dout, gammain=0.5, gammaout=  [ 1.0 ,  660.0 + 57.01   ,  4.0  ])

##     re3= system.GetRIXS( polarisationIn=[.0, 1.0,0.0, -1.0,.0] ,polarisationOut=[1.0,0.0], ein=ene,
##                         eout1=eout1, eout2=eout2, dout=dout, gammain=0.5, gammaout=  [ 1.0 ,  660.0 + 57.01   ,  4.0  ])


##     re[1]=re[1]+re1[1]+re2[1]+re3[1]

    
    file=open("out_"+sys.argv[1],"w")
    for d in re[1]:
        file.write(" %e "% d.imag)
    file.write("\n")
        

    file=open("outX_","w")

    nn=len(re[0])
    for v in Eins:
        file.write((" %e "*nn) % ((v,)*nn) )
        file.write("\n")

    
    file=open("outY_","w")


    for i in range(len(Eins)):
        file.write(" %e  "*len(re[0]) % tuple((re[0].tolist()) ))
        file.write("\n")






