from Numeric import *
import sys
import qt
import qwt
import iqt
from qwt.qplt import Plot, Curve, Pen, Red, Green , Yellow



def write(A, filename):
    A=list(A)
    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()


import string
import Numeric 
def leggi(nome):
    s=open(nome,"r").read()
    sl=string.split(s,"\n")
    ll=[]
    for l in sl:
        ll.append(map(string.atof,string.split(l)))
    del ll[-1]
    return Numeric.array(ll)

## ----------------------------------------------------------------------------------------------------
##     PARAMETERS FROM COWAN
##     When the program starts it search for out36_base and out36_exci files which are the rcn output for
##     the ground and excited configurations.
##     This initial choice is overwritten by following calls to the set(system) function if you
##     change those attributes.
## ---------------------------------------------------------------------------------
  
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(5,prefix="datas", nhopped=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 5 electrons are distributed in the 3d
               # shell in all possible combinations.
               # Giving nhopped=n as argument, besides the configurations with 5 electrons on 3d also the configurations
               # with 5+i electrons on 3d and 10-i on Ligands orbitals are considered ( 0<=i<=nhopped )
            
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). Dipolar transitions are calculated as
               # res[i], where i=1,2,3 and Mz= -1,0,1
    
res=system.GetSpectrum([ fm1,f0,f1 ]) # Absorption for a  definite polarisation. Fm1,f0,f1 are three coefficients.
               # They can be complex. So you can define any polarisation :
               # For Z polarisation [ fm1,f0,f1 ] = [ 0, 1, 0]
               # For X polarisation [ fm1,f0,f1 ] = [-1/sqrt(2.0) , 0 ,1/sqrt(2.0)]
               # For Y polarisation [ fm1,f0,f1 ] = [ 1.0j /sqrt(2.0) , 0 ,1.0j/sqrt(2.0)]
               # ..... and so on
    
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 init.py, the lines where such packages
               # are imported
   
Es, S2s, L2s, SL2, occPs , Szs, Lzs= 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
 
         """
    print s

from  System import System
# from SimplePlot import plot


def set(system):
    names = system.basenames+system.excinames+system.calcnames
    stop={}
    lungs = [len(system.basenames),len(system.excinames),len(system.calcnames)]
    stop[1]="---- BASE HAMILTONIAN ---"
    stop[1+lungs[0]]="---- EXCITED HAMILTONIAN ---"
    stop[1+lungs[0]+lungs[1]]="---- 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 : you must input a \"quoted\" string"
    user = input(" ")
    nonint_save(system, user)
    
def nonint_save(system, filename):
    names =   system.basenames+system.excinames+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")

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_F2"   ,  bs.FDict["(2p,3d)2"] ,
    "base_couche0_1_G1"   , bs.GDict["(2p,3d)1"]   ,
    "base_couche0_1_G3"   , bs.GDict["(2p,3d)3"]   ,

    "base_SO_0"           , bs.zetaDict["2p"]      ,
    "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_F2"     , ex.FDict["(2p,3d)2"]  , 
    "exci_couche0_1_G1"   , ex.GDict["(2p,3d)1"]   ,
    "exci_couche0_1_G3"   , ex.GDict["(2p,3d)3"]   ,
    "exci_SO_0"           , ex.zetaDict["2p"]      ,
    "exci_SO_1"           , ex.zetaDict["3d"]      ,
    "exci_Sop_Zero"       , 0.00001  ,
    "exci_Sop_Minus"       , 0.0000  ,
    "exci_Sop_Plus"       , 0.0000  ,
    "exci_counterDL"       , -4,
   ]



quadrunames = [
"dipole_couches_0_1_Km_-1",
"dipole_couches_0_1_Km_0",
"dipole_couches_0_1_Km_1",   
]

calculkvs = [
    "case" ,       "./",
    "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 ,
    "tolefact", 1.0e-6, 
    "shift_invert",  0,
    "nsearchedeigen",10,
    "NstepsTridiag", 250, 
    "Vs", 2.0,
    "Vp", 1.0,
    "VC0",0.2,
    "VC1",0.0,
    "DREF", 1.0,
    "ALPHAVC", -3.0,
    "ALPHAVSP", -3.0,
    "BONDS", [[-1.,0,0], [1.,0,0],[0,-1.,0],[0,1.,0],[0,0,-1.],[0,0,1.]],
    "factorhopexci",1.0,
    "facts_hop", None

    ]


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

counternames = [
    "base_counterDL",
    ]

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

print """type
help()
to get help
"""
