############################################
#  Alessandro MIRONE
#  2001
#  ESRF
########################################################
#
# dabax.py reads x-ray data tables compiled by Manuel Sanchez del Rio
#
# It offers also some extra facilities like KK, mainly used
# in packages like PPM, PyFullerton  ......
#
######################################################
#     CONTRIBUTIONS
#####################################################
#
#  Olivier Dhez,   4 Febraury 2003
#        added
#                  - BrennaCowan Tables  f1f2
#                  - FluorYield class implementig dabax's FluorYield_Krause.dat tables
#                  - class XrayEmission  
#                  - class XrayEmissionWeights  
#
#
#




ISPARALLEL=0
ISMASTER=1
try:
     import mpi
     ISPARALLEL=1
     if(mpi.rank==0):
       ISMASTER=1
     else:
       ISMASTER=0
except:
     pass


import string 
import Numeric 
import arrayfns             
import sys 
import copy
import specfile
import math
import os
import math
import FFT
import arrayfns

USE_KK_M=0


# da mettere a 1 per la diffrazione da cristalli

RECORDER_CONSIDER_THETA=0

try:
  import Minimiser
  import KKpy_c2
  print "  KK features availables "
  SIMPLE_Dabax=0
  print "  WARNING     WARNING      WARNING      WARNING      "
  print "  Dabax with minimising fonctionalities     "
  print "  Output data might be buffered     "
  print "  Disactivate the loading  of Minimiser module for a simpler use    "
  print " packages loaded is "
  print Minimiser.__file__
except:
  SIMPLE_Dabax=1



LAMBDAEV=12398.52
AVOGADRO=6.0221367e23

print " DABAX VERSION 2003-04-2"

if(os.getenv("DABAX_DIR")!=None):
  dabaxpath=os.getenv("DABAX_DIR")
else:
#   dabaxpath="/scisoft/ESRF_sw/linux_i386/PACKAGES/PPM/DATA/"
    dabaxpath="/scisoft/DABAX/data"

def par(a):
    if( hasattr(a,'getvalue' )):
          return 1.0*a.getvalue()
    else:
          return 1.0*a      



######################################################################
#  
#  class to manipulate f0/f12 objects

class Dabax_Scatterer:
# DWsigma, Debye Waller is the sigma of the atom vibration


   is_an_index=0
   is_a_tensor=0

   def __init__(self, f0list=None,f0weights=None, f1f2list=None,f1f2weights=None,
                convention="+iKx" , DWsigma=None ):
      self.DWsigma=DWsigma
      if(f0list!=None):
          self.f0list      = f0list
      else:
          self.f0list      = []

      if(f1f2list is not None):
          self.f1f2list    = f1f2list
      else:
          self.f1f2list      = []

      if(f0weights is not None):
          self.f0weights   = f0weights
      else:
           self.f0weights     = []


      if(f1f2weights is not None):
          self.f1f2weights = f1f2weights
      else:
          self.f1f2weights      = []

      self.classe="Scatterer"
      self.convention= convention

      self.recorder = {}
      self.recorderer_key_list=[]

   def __repr__(self):

     for key in dir(self):
       if(key not in ['self', ]): print "            ",  key, "=", getattr(self,key)
     return ''

   def __add__(self, other):
      if("classe" not in dir(self) or self.classe!= "Scatterer"):
         print " summing wrong class to Scatterer"
         sys.exit(1)

      newf0list= self.f0list +other.f0list
      newf0weights= Numeric.array(list(self.f0weights)+list(other.f0weights))

      newf1f2list= self.f1f2list+other.f1f2list
      newf1f2weights= Numeric.array(list(self.f1f2weights) +list(other.f1f2weights))

      self.recorder = {}
      self.recorderer_key_list=[]

      return Dabax_Scatterer(newf0list,newf0weights,newf1f2list,newf1f2weights)


   def Z(self):
       res=0

       for i  in range(0,len(self.f0list) ):
            res=res+  self.f0list[i].f0Energy( 1000.0 )*self.f1f2weights[i]
       return  res

   def __rmul__(self, fact):

      self.recorder = {}
      self.recorderer_key_list=[]

      return self.__mul__(fact)

   def __mul__(self, fact):
      newf0list  = (self.f0list)
      newf1f2list= (self.f1f2list)

      newf0weights   = fact*self.f0weights
      newf1f2weights = fact*self.f1f2weights

      self.recorder = {}
      self.recorderer_key_list=[]

      return Dabax_Scatterer(newf0list,newf0weights,newf1f2list,newf1f2weights)


   def F_Lambda(self,  Lambda, theta=None ): 
      return self.  F_Energy(  LAMBDAEV/Lambda, theta )


   def F_Energy(self,  En, theta=None ):   

      if(len(self.recorderer_key_list) > 10 ):
       for rem_key in self.recorderer_key_list[10:] :
           del self.recorder[rem_key]
       self.recorderer_key_list =self.recorderer_key_list[:10] 

      if( isinstance(En,Numeric.ArrayType)  ):
        if(theta is None):
           theta=Numeric.ones( len(En)     )*RECORDER_CONSIDER_THETA*math.pi/180.0

        key = En.tostring() + theta.tostring()

      else:
        if(theta is None):
           theta=0*math.pi/180.0
        key = (En, theta)


      if(key not in self.recorder.keys() ):

        res=complex(0,0)

        for i  in range(0,len(self.f0list) ):
            res=res+self.f0weights[i]*self.f0list[i].f0Energy(  En, theta)

        for i  in range(0,len(self.f1f2list) ):
            res=res+self.f1f2weights[i]*self.f1f2list[i].f1f2Energy(  En)

        if(self.DWsigma!=None):
            q=2*Numeric.sin(theta)*LAMBDAEV/En
            res=res*Numeric.exp(- (2*math.pi*math.pi*self.DWsigma*self.DWsigma)* q*q      )

        self.recorder[key]= res
        self.recorderer_key_list.append(key)
      else:
        res= self.recorder[key]


      if(self.convention=="+iKx"): 
         return Numeric.conjugate(res)
      else:
         return res


#
#############################################################################################






class return_value:
  pass
def IndexFromTable(Tablef0, Tablef12, *args):
   """  

     FUNCTION:

     IndexFromTable(Tablef0, Tablef12, *args).

      The arguments are 2*nel in number
     where nel is the number of elements
     entering the composition.

     nome_el=args[2*i]
     density = args[2*i+1]


     RETURN:
        an object containing

   res.DabaxList = DabaxList
   res.DensityList = DensityList
   res.MassList= MassList

   this object is used to define a material
   by its optical indexes as in  object PPM_ElementalLayer
   of library PPM or in object KK.

   """
   res=return_value()
   res.is_an_index=0
   nel=len(args)
   if( nel%2 !=0):
      raise " Wrong Number of arguments in IndexFromTable"
   nel=nel/2

   DabaxList=[]
   DensityList=[]
   MassList=[]

   for i in range(nel):
     nome_el=args[2*i]
     density = args[2*i+1]
     f0  = Tablef0.Element(nome_el)
     f1f2= Tablef12.Element(nome_el)
     mass= AtomicProperties(nome_el, 'AtomicMass')
     scatterer=Dabax_Scatterer((f0,),(1,), (  f1f2,),(1,)       )

     DabaxList.append(scatterer)
     DensityList.append(density)
     MassList.append(mass)
   res.DabaxList = DabaxList
   res.DensityList = DensityList
   res.MassList= MassList
   return res



def  CalculateIndexesFromTable( object  , wavelenght):
      """ 
         given an object in the form returned by IndexFromTable 
       ( an object containing:

          res.DabaxList = DabaxList
          res.DensityList = DensityList
          res.MassList= MassList
        )
        calculates the optical index at wavelenght(s)

      """
      self = object
      indexes=Numeric.zeros(Numeric.shape(wavelenght),Numeric.Complex)
      for i in range(0,len(self.DabaxList)):
         scatterer=self.DabaxList[i]
         density  =par(self.DensityList[i])
         mass     =self.MassList[i]
         if( scatterer.is_an_index ):
           add = 1.0-scatterer.index(wavelenght)
         else:
           if( scatterer.__class__ == return_value ):
             add= 1- CalculateIndexesFromTable( scatterer,wavelenght )
           else:
             add=415.22*(  scatterer.F_Lambda(wavelenght) )*( 
                             wavelenght*wavelenght/12398.52/12398.52/mass)

         indexes=indexes+ add*density


      return 1-indexes  







class KK:
  """
      CLASS:

      it returns index(wavelenghts) as the result of a KK
      transformation.

      *** def __init__(self, nomefile_or_betaObject=None, material=None,
                   E1=None, E2=None, N=None,
                    e1=None, e2=None, 
                    Fact=None, Dx=None):

      --- material is an object like res returned by IndexFromTable :

         res.DabaxList = DabaxList
         res.DensityList = DensityList
         res.MassList= MassList

      --- E1 E2 are the extrema between which the KK has to be done.
          A file containing N energy points and their beta is generated
          from material object
          and given in input to the kk program

      --- e1 e2 specifies an interval contained into E1,E2.
          Doing the KK transforms, beta values inside e1,e2
          are read from file named nomefile.

      --- Fact is a factor multiplying the betas read from nomefile.


  """

  is_a_tensor=0
  is_an_index = 1

  def __init__(self, filename_or_betaObject=None, material=None,
               E1=None, E2=None, N=None,
               e1=None, e2=None, 
               Fact=None, RelativeDensity=1.0, Dx=None, maglia=None, Nmaglia=10 ,idxr1=None, idxr2=None):

    if(filename_or_betaObject is None):
       return

    self.nomefile_or_betaObject = filename_or_betaObject
    self.material = material

    self.E1      = 1.0*E1
    self.E2      = 1.0*E2
    self.N       = N
    self.idxr1 = idxr1
    self.idxr2 = idxr2

    energies = self.E1 + Numeric.arrayrange(self.N)*(self.E2-self.E1)/self.N
    if(maglia is None):
      self.energies_for_inte=energies
      self.Nm = self.N
    else:
      if( hasattr(maglia, "maglia") ):
           maglia=maglia.maglia
      X=[maglia[0].X ]
      for i in range(1,len(maglia) ):
        for k in range(Nmaglia):
          X.append(  ( (Nmaglia-k-1.0)*maglia[i-1].X   +  (k+1) * maglia[i].X   )/ (Nmaglia-0.0)         )

    tmpar=(X + energies.tolist() )
    tmpdic={}
    for fld in tmpar:
         tmpdic[fld]=1
    tmpar=tmpdic.keys()
    tmpar.sort()
    
    self.energies_for_inte = Numeric.array(   tmpar    )
    self.Nm = len(self.energies_for_inte)
      
    self.e1      = 1.0*e1
    self.e2      = 1.0*e2
    self.Fact    = Fact
    self.Dx      = Dx

    if( USE_KK_M):
      dimdata = len(self.energies_for_inte)
      self.data_kk_m = Numeric.zeros([dimdata+1 , dimdata ] ,"d")
      self.kk_m_calculated= 0

    self.oldDx=0
    self.oldFact=0
    self.oldE1=0
    self.oldE2=0
    self.oldN=0
    self.olde1=0
    self.olde2=0

    self.oldWaves = 0
    self.oldRes = 0
    self.colpo=0

    self.DabaxList=[self]
    self.DensityList=[RelativeDensity]
    self.MassList=["does not matter"]

    self.rem_wave={}

  def index(self, wavelenghts):

      fai_kk=1
      sw = wavelenghts.tostring()

      if(type(self.nomefile_or_betaObject)==type("s") ):

         if( (self.oldDx,self.oldFact,self.oldE1,self.oldE2,self.oldN,self.olde1,self.olde2 )==
           (par(self.Dx),par(self.Fact),self.E1,self.E2,self.N,self.e1,self.e2 )):
           fai_kk=0
           if(sw in self.rem_wave.keys()):
              return self.rem_wave[sw]
         else:
           self.rem_wave={}
      else:

         if(hasattr(self.nomefile_or_betaObject,"contribution")==0):
             raise " bad object passed to KK. it should be either a filename or a beta object "

#         print "Minimiser.colpo[0]    ", Minimiser.colpo[0]
#         print "self.colpo    ", self.colpo

         if(Minimiser.colpo[0]== self.colpo ):
             fai_kk=0
             if(sw in self.rem_wave.keys()):
                return self.rem_wave[sw]
         else:
           self.rem_wave={}


      self.colpo = Minimiser.colpo[0]
      indexeses     = Numeric.zeros(len(wavelenghts),"D")
      ZmolTot=0
      DensTot=0

      for i in range(0,len(self.material.DabaxList)):
          scatterer=self.material.DabaxList[i]
          density  =par(self.material.DensityList[i])
          mass     =self.material.MassList[i]

          if( hasattr(scatterer,"Z")==0):
            raise " Z() function missing for object scatterer in function index(), file dabax.py "

          ZmolTot = ZmolTot + scatterer.Z() * density/mass



          add=415.22*(  scatterer.F_Lambda( wavelenghts ) )*( 
                             wavelenghts*wavelenghts /12398.52/12398.52/mass)
          indexeses   =   indexeses  + add*density





      indexeses     = 1 - indexeses


      ## if( fai_kk==0):
      ##   raise " fai_kk==0"
      if(fai_kk):


        indexes_Table = Numeric.zeros(self.Nm,"D")

        energies= self.energies_for_inte

        wavelenghts_Table = 12398.52 / energies

        for i in range(0,len(self.material.DabaxList)):
          scatterer=self.material.DabaxList[i]
          density  =par(self.material.DensityList[i])
          mass     =self.material.MassList[i]

          add=415.22*(  scatterer.F_Lambda( wavelenghts_Table ) )*( 
                             wavelenghts_Table*wavelenghts_Table /12398.52/12398.52/mass)
          indexes_Table   =   indexes_Table  + add*density

        faikk_c=1
        if faikk_c:
          indexes_Table = 1 - indexes_Table
          dati = self.nomefile_or_betaObject.contribution(energies)

          if( USE_KK_M):
            if(self.kk_m_calculated== 0):
              print "chiamo KK_M "
              (emin,emax, imin, imax)=KKpy_c.KKpy_M_c(self.data_kk_m, energies, self.e1, self.e2 , ZmolTot)


              
              print "chiamo KK_M OK"
              self.kk_m_calculated=1
              self. imin=imin
              self.imax =imax
              self. emin=emin
              self.emax =emax
            else:
              emin= self. emin
              emax= self. emax
              imin= self. imin
              imax= self. imax
              
            dim_ = len(self.data_kk_m)-1
            transf=Numeric.zeros( [dim_]   ,"d")
            print " MOLTILICAZIONE "
            transf[imin:imax+1]=Numeric.dot(self.data_kk_m[imin:imax+1 ,0:dim_],
                                            dati  )
            transf[imin:imax+1] = transf[imin:imax+1]+self.data_kk_m[dim_  ][imin:imax+1]

            # transf = Numeric.dot(Numeric.transpose(self.data_kk_m[:-1,:])
            #                                 , dati  )
            # transf = transf+self.data_kk_m[dim_  ]

          else:
            # if(1):
            if(ISPARALLEL==0):
              (transf, emin,emax, imin, imax) = KKpy_c2.KKpy_c(dati, energies, self.e1, self.e2 , ZmolTot)


              tf=open("kkfft2.dat","w")
              for x,y , d in zip(energies,transf, dati ):
                   tf.write(str(x)+"  "+str(y)+" "  +str(d)+"\n")
              tf.close()



              
#              print energies
##               f=open("deltest.del","w")
##               for i in range(len(transf)):
##                    f.write("%e %e  %e\n"%(  energies[i],  transf[i],dati[i]))
##               f.close()


            else:
              inizio =  (self.e1 *(mpi.size-(mpi.rank-0.00001)  )   +  self.e2 *(mpi.rank-0.00001) )/ mpi.size
              fine   =  (self.e1 *(mpi.size-(mpi.rank+1+0.00001)  )   +  self.e2 *(mpi.rank+1+0.00001) )/ mpi.size
              (transf, emin,emax, imin, imax) = KKpy_c.KKpy_c(dati, energies, inizio, fine , ZmolTot)
              tosend = transf[imin:imax+1].tostring()
              ( eminT,emaxT, iminT, imaxT)=( emin,emax, imin, imax)
              for rango in range(mpi.size):
                if(mpi.rank==rango):
                  mpi.bcast(emin,rango)
                  mpi.bcast(emax,rango)
                  mpi.bcast(imin,rango)
                  mpi.bcast(imax,rango)
                  mpi.bcastString(tosend,rango  )
                else:
                  imin_rec=0
                  imax_rec=0
                  emin_rec=0.
                  emax_rec=0.
                  
                  received=""
                  emin_rec=  mpi.bcast(emin_rec,rango)
                  emax_rec=  mpi.bcast(emax_rec,rango)
                  imin_rec=mpi.bcast(imin_rec,rango)
                  imax_rec=mpi.bcast(imax_rec,rango)
                  received=mpi.bcastString(received,rango  )
                  transf[imin_rec:imax_rec+1] = Numeric.fromstring(received,"d")

                  if( emin_rec<eminT): eminT=emin_rec
                  if( emax_rec>emaxT): emaxT=emax_rec
                  if( imin_rec<iminT): iminT=imin_rec
                  if( imax_rec>imaxT): imaxT=imax_rec
              imin=iminT
              imax=imaxT
              emin=eminT
              emax=emaxT
              

            
          x = energies
          # print emin,emax, imin, imax
          # transf[imin:imax+1] =  transf[imin:imax+1]
          if(self.idxr1 is  None ):
               DF1 = transf[imin] - indexes_Table[imin].real
               DF2 = transf[imax] - indexes_Table[imax].real
          else:
               DF1 = transf[imin] - par(self.idxr1)
               DF2 = transf[imax] - par(self.idxr2)
              
          transf[imin:imax+1] = transf[imin:imax+1] - arrayfns.interp(Numeric.array([DF1,DF2]), Numeric.array([emin,emax]), x[imin:imax+1])
          
          # transf=      indexes_Table.real 
          rp = transf
          ip = dati
        else:
          raise " KK with batch file removed "
        self.rem_x=x
        self.rem_rp=rp
        self.rem_ip=ip
        self.rem_emin = emin
        self.rem_emax = emax

      else:
        x  = self.rem_x
        rp = self.rem_rp
        ip = self.rem_ip
        emin = self.rem_emin
        emax = self.rem_emax


      # print " X is ", x
      # print " wavelenghts is ", wavelenghts
      ens = 12398.52/wavelenghts
      masque = 1- Numeric.less(emax , ens) - Numeric.less(  ens, emin)

      nrp = arrayfns.interp(rp,x,ens) * masque
      nip = arrayfns.interp(ip,x,ens) * masque

      nrp = nrp+ indexeses.real  * (1-masque)
      nip = nip+ indexeses.imag  * (1-masque)

      # file=open("tutto","w")
      # for i in range(len(wavelenghts)):
      #   file.write("  %e %e %e  \n" % ( 12398.52/wavelenghts[i], nrp[i] ,nip[i]   )   )
      # file.close()

      p=nrp+complex(0,1)*nip

      (self.oldDx,self.oldFact,self.oldE1,self.oldE2,self.oldN,self.olde1,self.olde2 )= \
          (par(self.Dx),par(self.Fact),self.E1,self.E2,self.N,self.e1,self.e2 )

      self.oldWaves=wavelenghts.tostring()
      self.oldRes=p

      self.rem_wave[sw]=p

      return p



def kkfft(x, a):

     print " il passo est " ,  (x[-1]- x[0]  )/(len(a)-1)

     period=      (x[-1]- x[0]  )/(len(a)-1)

     conv = Numeric.zeros( len(a), "d")

     conv[1:] = 1.0/ (x[1:]-x[0])  - 1.0/( x[-1]+period-x[1:]) 

     afft=FFT.fft(a)
     convfft=FFT.fft(conv)

     result=1.0/math.pi  * FFT.inverse_fft(afft*convfft )

     return -result*period


class KKFFT:
  """
      CLASS:

      it returns index(wavelenghts) as the result of a KK
      transformation.

      *** def __init__(self, nomefile_or_betaObject=None, material=None,
                   E1=None, E2=None, N=None,
                    e1=None, e2=None, 
                    Fact=None, Dx=None):

      --- material is an object like res returned by IndexFromTable :

         res.DabaxList = DabaxList
         res.DensityList = DensityList
         res.MassList= MassList

      --- E1 E2 are the extrema between which the KK has to be done.
          A file containing N energy points and their beta is generated
          from material object
          and given in input to the kk program

      --- e1 e2 specifies an interval contained into E1,E2.
          Doing the KK transforms, beta values inside e1,e2
          are read from file named nomefile.

      --- Fact is a factor multiplying the betas read from nomefile.


  """

  is_a_tensor=0
  is_an_index = 1

  def __init__(self, filename_or_betaObject=None, material=None,
               E1=None, E2=None, N=None,
               e1=None, e2=None, 
               Fact=None, RelativeDensity=1.0, Dx=None, maglia=None, Nmaglia=10 ,idxr1=None, idxr2=None):

    if(filename_or_betaObject is None):
       return

    self.nomefile_or_betaObject = filename_or_betaObject
    self.material = material

    self.E1      = 1.0*E1
    self.E2      = 1.0*E2
    self.N       = N
    self.idxr1 = idxr1
    self.idxr2 = idxr2
    self.e1      = 1.0*e1
    self.e2      = 1.0*e2
    self.Fact    = Fact
    self.Dx      = Dx

    
    if(maglia is None):
         energies = self.E1 + Numeric.arrayrange(self.N)*(self.E2-self.E1)/self.N
         self.Nm = self.N
    else:
         print " setto la maglia "
         if( hasattr(maglia, "maglia") ):
              maglia=maglia.maglia
         DX = (maglia[-1].X-maglia[0].X   )/Nmaglia/len(maglia)
         Npunti = int((self.E2-self.E1)/DX)
         energies = self.E1 + Numeric.arrayrange(Npunti)*DX
         self.Nm = Npunti
         
    self.energies_for_inte=energies
    
      
    self.oldDx=0
    self.oldFact=0
    self.oldE1=0
    self.oldE2=0
    self.oldN=0
    self.olde1=0
    self.olde2=0

    self.oldWaves = 0
    self.oldRes = 0
    self.colpo=0

    self.DabaxList=[self]
    self.DensityList=[RelativeDensity]
    self.MassList=["does not matter"]

    self.rem_wave={}

  def index(self, wavelenghts):

      print " KKFFT  " 

      fai_kk=1
      sw = wavelenghts.tostring()

      if(type(self.nomefile_or_betaObject)==type("s") ):

         if( (self.oldDx,self.oldFact,self.oldE1,self.oldE2,self.oldN,self.olde1,self.olde2 )==
           (par(self.Dx),par(self.Fact),self.E1,self.E2,self.N,self.e1,self.e2 )):
           fai_kk=0
           if(sw in self.rem_wave.keys()):
              return self.rem_wave[sw]
         else:
           self.rem_wave={}
      else:

         if(hasattr(self.nomefile_or_betaObject,"contribution")==0):
             raise " bad object passed to KK. it should be either a filename or a beta object "

#         print "Minimiser.colpo[0]    ", Minimiser.colpo[0]
#         print "self.colpo    ", self.colpo

         if(Minimiser.colpo[0]== self.colpo ):
             fai_kk=0
             if(sw in self.rem_wave.keys()):
                return self.rem_wave[sw]
         else:
           self.rem_wave={}

      self.colpo = Minimiser.colpo[0]
      indexeses     = Numeric.zeros(len(wavelenghts),"D")
      ZmolTot=0
      DensTot=0

      for i in range(0,len(self.material.DabaxList)):
          scatterer=self.material.DabaxList[i]
          density  =par(self.material.DensityList[i])
          mass     =self.material.MassList[i]

          if( hasattr(scatterer,"Z")==0):
            raise " Z() function missing for object scatterer in function index(), file dabax.py "


          #    lo z andrebbe correto per il relativistico Z* = Z-(Z/82.5)**2.37
          # ZmolTot = ZmolTot + scatterer.Z() * density/mass
          tmp_Z =  scatterer.Z()
          ZmolTot = ZmolTot + ( tmp_Z - ( tmp_Z/82.5 )**2.37   ) * density/mass


          add=415.22*(  scatterer.F_Lambda( wavelenghts ) )*( 
                             wavelenghts*wavelenghts /12398.52/12398.52/mass)
          indexeses   =   indexeses  + add*density





      indexeses     = 1 - indexeses
      if(fai_kk):


        indexes_Table = Numeric.zeros(self.Nm,"D")

        energies= self.energies_for_inte

        wavelenghts_Table = 12398.52 / energies

        for i in range(0,len(self.material.DabaxList)):
          scatterer=self.material.DabaxList[i]
          density  =par(self.material.DensityList[i])
          mass     =self.material.MassList[i]

          add=415.22*(  scatterer.F_Lambda( wavelenghts_Table ) )*( 
                             wavelenghts_Table*wavelenghts_Table /12398.52/12398.52/mass)



          indexes_Table   =   indexes_Table  + add*density

        
        if 1:
          indexes_Table = 1 - indexes_Table
          dati = self.nomefile_or_betaObject.contribution(energies)


          if 1:
            # if(1):
            if(1):


              # quindi trasfromazione di dati/energies a cui si aggiunge -zfatt*ZmolTot/e[i]/e[i]
              # dove 
              # zfatt = 12398.52*12398.52*  0.529/137.036/137.036 *0.6022/2/PI
              # dopo di che emin emax sono due punti della maglia energies
              # che stanno all ' interno di e1 e2
              # imin imax sono i loro indici
              transf = (kkfft(energies, dati)).real


              tf=open("kkfft.dat","w")
              for x,y in zip(energies,transf ):
                   tf.write(str(x)+"  "+str(y)+"\n")
              tf.close()


##               (transf2, emin2,emax2, imin2, imax2) = KKpy_c.KKpy_c(dati, energies, self.e1, self.e2 , ZmolTot)


              
##               raise " OK " 
              zfatt = 12398.52*12398.52*  0.529/137.036/137.036 *0.6022/2/math.pi

              transf=transf - zfatt*ZmolTot / energies/energies

              
              imin = arrayfns.interp(Numeric.arrayrange(len(energies)), energies, [self.e1] )[0]
              imin = int(imin)
              imax = arrayfns.interp(Numeric.arrayrange(len(energies)), energies, [self.e2] )[0]
              imax = int(imax)+1

##               f=open("del.del","w")
##               for i in range(imin, imax):
##                    f.write("%e %e %e %e %e \n"%(energies[i], transf[i], transf2[i], dati[i] , - (zfatt*ZmolTot / energies/energies)[i]))
##               f.close()

              emin=energies[imin]
              emax=energies[imax]

              if( emin> self.e1  or emax < self.e2 ) :
                   raise " PROBLEM WITH INTERVAL e1 e2 "

              
                   
              
              # (transf, emin,emax, imin, imax) = KKpy_c.KKpy_c(dati, energies, self.e1, self.e2 , ZmolTot)

            
          x = energies
          # print emin,emax, imin, imax
          # transf[imin:imax+1] =  transf[imin:imax+1]
          if(self.idxr1 is  None ):
               DF1 = transf[imin] - indexes_Table[imin].real
               DF2 = transf[imax] - indexes_Table[imax].real
          else:
               DF1 = transf[imin] - par(self.idxr1)
               DF2 = transf[imax] - par(self.idxr2)


          transf[imin:imax+1] = transf[imin:imax+1] - arrayfns.interp(Numeric.array([DF1,DF2]), Numeric.array([emin,emax]), x[imin:imax+1])




          
          # transf=      indexes_Table.real 
          rp = transf
          ip = dati
        self.rem_x=x
        self.rem_rp=rp
        self.rem_ip=ip
        self.rem_emin = emin
        self.rem_emax = emax

      else:
        x  = self.rem_x
        rp = self.rem_rp
        ip = self.rem_ip
        emin = self.rem_emin
        emax = self.rem_emax


      
      ens = 12398.52/wavelenghts
      masque = 1- Numeric.less(emax , ens) - Numeric.less(  ens, emin)

      nrp = arrayfns.interp(rp,x,ens) * masque
      nip = arrayfns.interp(ip,x,ens) * masque

      nrp = nrp+ indexeses.real  * (1-masque)
      nip = nip+ indexeses.imag  * (1-masque)

      p=nrp+complex(0,1)*nip

      (self.oldDx,self.oldFact,self.oldE1,self.oldE2,self.oldN,self.olde1,self.olde2 )= \
          (par(self.Dx),par(self.Fact),self.E1,self.E2,self.N,self.e1,self.e2 )

      self.oldWaves=wavelenghts.tostring()
      self.oldRes=p

      self.rem_wave[sw]=p

      return p

try:
     USV = os.getenv("USEKKFFT")
     if USV is not None:
          KK=KKFFT
except:
     pass




def KK_write(KKobj,minE, maxE,n , outn, *args):
   obj = KK_ForPlot(KKobj,minE, maxE,n )
   obj.error()
   obj.write2File(outn, args)


def Beta_write(Bobj,minE, maxE,n , outn, pol=1, dichroism=0):
   a=Numeric.arrayrange(minE, maxE, (maxE-minE+0.0)/n   )
   res = Bobj.contribution(a, pol, dichroism)

   f=open("%s" %(outn),"w")
   for i in range(len(a)):
        f.write("%e %e \n"%( a[i],  res[i] ) )
   f.close()



class KK_ForPlot:

  classe="KK_ForPlot"

  def __init__(self,KKobj,minE, maxE,n ,printpartial=0 ):

      self.KKobj=KKobj
      self.minE=minE
      self.maxE=maxE
      self.n=n
      self.printpartial=printpartial

      energies = minE +  Numeric.arrayrange(n)*(maxE-1.0*minE)/n
      self.wavelenghts=12398.52/energies


  def write2File(self, filename, args , mode="w"):
      f=open("%s" %(filename),mode)
      for a in args:
          f.write("%e "% a)
      for i in range(len(self.wavelenghts)):
        f.write("%e %e %e \n"%(  12398.52/self.wavelenghts[i],  self.calculatedscan[0][i]-1, self.calculatedscan[1][i] ) )
      f.close()


###############################################################
#  the following commented code is taken from the
#  visualisation part of beansgui.py
#
#       calculated=fitobject.calculatedscan[count]
#        scan_1 = fitobject.scanlist[count][0]
#        scan_2 = fitobject.scanlist[count][1]
#        scan = fitobject.scanlist[count][2]
#        if(scan is None):
#          scan=calculated


  def error(self, nopartial=0):

     error=0

     if( hasattr(self, "calculatedscan")==0): self.calculatedscan=[]

     newscans=[]

     wavelenghts= self.wavelenghts

     if( self.KKobj.__class__ == return_value ):
       r=  CalculateIndexesFromTable(self.KKobj ,wavelenghts ) 
     else:
       r = self.KKobj.index(wavelenghts)



     newscans.append(r.real)

     newscans.append(r.imag)

     self.scanlist=[ [   wavelenghts,  wavelenghts,None   ],[wavelenghts,  wavelenghts,None]     ]

     # print " WAVE IS ", wavelenghts

     self.calculatedscan=newscans
     if(self.printpartial and nopartial==0):
         self.write2File("KKresult")

     return 0

###############################################################################################



class Dabax_f0_Table:

  """ Object to retrieve elements f0 tables.
      Creator is Dabax_f0_Table(Name)

      Name can be either  "f0_WaasKirf.dat" ( mind the quotes)
      or "f0_CromerMann.dat" : the two kind of tables implemented


      The function Dabax_f0_Table.Element(self, "name of element") returns
      an object having the methods 
                 f0Lambda(self, Lambda, theta=None)
                 f0Energy(self, Energy, theta=None):              

  """


  def __init__(self, Name=None):
#    self.path="/home/alessandro/data/f0"
    self.path=dabaxpath

    self.dataprocessors   ={   "f0_WaasKirf.dat": self.process_WaasKirf,
                             "f0_CromerMann.dat": self.process_WaasKirf
                           }
    self.elementretrievers={  "f0_WaasKirf.dat":   self.Element_WaasKirf ,
                            "f0_CromerMann.dat":   self.Element_WaasKirf  }
    if(Name!=None):
      self.load(Name)    


  def getDataBaseNames(self):
    return self.elementretrievers.keys()

  def load(self, location):
     self.Location=self.path+"/"+ location
     self.fileType=location
     self.data=specfile.Specfile(self.Location)

     print "Opened file %s. It is of the  %s type"%(
                          self.Location,self.fileType)


     if( (self.fileType in self.dataprocessors.keys()) and
         (self.fileType in self.elementretrievers.keys()) ):

        self.Element= self.elementretrievers[self.fileType]
        self.dataprocessors[self.fileType]()

     else:
        print " I can't process a file of the type %s", self.fileType
        print " Possible types are ", self.fileType in dataprocessors.keys
        exit(1)

  def process_WaasKirf(self):

     self.elementlist=[]
     for scan in self.data:
       self.elementlist.append( string.split(scan.command()) [-1])

  def Element_WaasKirf(self,name):
      index=self.elementlist.index(name)
      if(index==-1):
         print "%s not found in %s" % ( name, self.Location)
         exit(1)
      lista=self.data[index].dataline(1)
      res= f0_WaasKirf( lista[5],lista[0:5], lista[6:11]   )
      res.element=name
      return res


###############################################
#
#  class to get the f0 from WaasKirf (internal use only)

class f0_WaasKirf:
  def __init__(self, c,a,b):
    self.a_array=Numeric.array(a)
    self.b_array=(Numeric.array(b))
    self.c=c
    self.classe="class to get the f0 from WaasKirf (internal use only)"

  def __repr__(self):
     for key in dir(self):
       if(key not in ['self','a_array','b_array','c' ]): print "            ",  key, "=", getattr(self,key)
     return ''

  def f0Lambda(self, Lambda, theta=None):
      if( isinstance(Lambda,Numeric.ArrayType)  ):
        if(theta is None):
           theta=Numeric.ones( len(En)     )*0*math.pi/180.0
        k=  Numeric.sin(theta) / Lambda
        k=-k*k
        return (self.c+Numeric.sum( Numeric.swapaxes(self.a_array* Numeric.exp(  Numeric.reshape(k,[len(k),1])* self.b_array),0,1)   ) )
      else:
        if(theta is None):
           theta=0*math.pi/180.0
        k=  Numeric.sin(theta) / Lambda
        k=-k*k
        return (self.c+ Numeric.sum(self.a_array*Numeric.exp( self.b_array*k) ) )


  def f0Energy(self,  Energy, theta=None):
      return self.f0Lambda(LAMBDAEV /Energy, theta)


#########################################################################################
#
#  f12  tables

class Dabax_f1f2_Table:

  """ Object to retrieve elements f0 tables.
      Creator is Dabax_f1f2_Table(Name)

      Name can be either  "f1f2_Sasaki.dat" ( mind the quotes)
      or "f1f2_Windt.dat" : the two kind of tables implemented


      The function Dabax_f1f2_Table.Element(self, "name of element") returns
      an object having the methods 
                 f1f2Lambda(self, Lambda, theta=None)
                 f1f2Energy(self, Energy, theta=None):              

  """

  def __init__(self, Name=None):
#    self.path="/home/alessandro/data/f12"
    self.path=dabaxpath
    self.dataprocessors   ={"f1f2_Windt.dat": self.process_Windt,"f1f2_Sasaki.dat":self.process_Sasaki, 
                "f1f2_Henke.dat": self.process_Windt,"f1f2_asf_Kissel.dat":self.process_Windt,
                "f1f2_BrennanCowan.dat":self.process_Sasaki}
    self.elementretrievers={"f1f2_Windt.dat": self.Element_Windt,"f1f2_Sasaki.dat":self.Element_Windt ,
                            "f1f2_Henke.dat": self.Element_Windt, "f1f2_asf_Kissel.dat": self.Element_Kissel,
                            "f1f2_BrennanCowan.dat":self.Element_Brennan}
    if(Name!=None):
      self.load(Name)    

  def getDataBaseNames(self):
    return self.elementretrievers.keys()

  def load(self, location):

     self.Location=self.path+"/"+ location
     self.fileType=location
     self.data=specfile.Specfile(self.Location)

     print "Opened file %s. It is of the  %s type"%(
                          self.Location,self.fileType)


     if( (self.fileType in self.dataprocessors.keys()) and
         (self.fileType in self.elementretrievers.keys()) ):

        self.Element= self.elementretrievers[self.fileType]
        self.dataprocessors[self.fileType]()

     else:
        msg="***********************************************************\n"
        msg = msg+ " I can't process a file of the type %s\n" % self.fileType
        msg=msg+" Possible types are %s\n STOP" % self.dataprocessors.keys()
        raise msg


  def  process_Windt(self): 
     self.elementlist=[]
     for scan in self.data: 
       self.elementlist.append( string.split(scan.command()) [-1])

  def  process_Sasaki(self): 
     self.elementlist=[]
     for scan in self.data: 
       self.elementlist.append( string.split(scan.command()) [-1])

  def Element_Windt(self,name):
     self.labelE = 'PhotonEnergy[eV]'
     self.factorE = 1.0
     self.labelf1 = 'f1'
     self.labelf2 = 'f2'
     return self.Element_generic(name)

  def Element_Brennan(self,name):
     self.labelE = ' Energy[eV]'
     self.factorE = 1.0
     self.labelf1 = 'f1'
     self.labelf2 = 'f2'
     return self.Element_generic(name)

  def Element_Kissel(self,name):
     self.labelE = 'PhotonEnergy[KeV]'
     self.factorE = 1000.0
     self.labelf1 = 'f1'
     self.labelf2 = 'f2'
     return self.Element_generic(name)

  def Element_generic(self,name):
      index=self.elementlist.index(name)
      if(index==-1):
         print "%s not found in %s" % (name, self.Location)
         exit(1)

      ene=self.data[index].datacol(self.labelE)*self.factorE
      f1 =self.data[index].datacol(self.labelf1)
      f2 =self.data[index].datacol(self.labelf2)


      hea=self.data[index].header("")
      head=""
      for tok in hea: head=head+tok

      if(string.find(head, "1ADD 0.0")>0):
         # retrieving Z
         Z = AtomicProperties(name, 'Z')
         f1=f1-Z
#                -----
      res= Ef1f2_Interpolated(ene, f1,  f2  )
      res.element=name
      return res

#
#########################################################################################




###############################################
#
#  class to get f1f2 (internal use only)

class Ef1f2_Interpolated:
  def __init__(self, ene,f1,f2):
    self.ene=ene
    self.f1=f1
    self.f2=f2
    self.classe="Ef1f2_Interpolated class to get f1f2 (internal use only) "

  def __repr__(self):
     for key in dir(self):
       if(key not in ['self','ene','f1','f2' ]): print "            ",  key, "=", getattr(self,key)
     return ''

  def f1f2Lambda(self, Lambda):

    energy= LAMBDAEV/Lambda
    return self.f1f2Energy(energy)

  def f1f2Energy(self, energy):
    if(isinstance(energy,Numeric.ArrayType)):
      return arrayfns.interp(self.f1,self.ene, energy )+complex(0,1)*arrayfns.interp(self.f2,self.ene, energy ) 
    else:
      return complex(arrayfns.interp(self.f1,self.ene,(energy,))[0],arrayfns.interp(self.f2,self.ene,(energy,))[0])
  
  def f1Energy(self, energy):
        # Olivier Dhez
        return arrayfns.interp(self.f1,self.ene,(energy,))[0]

  def f2Energy(self, energy):
        # Olivier Dhez
        return arrayfns.interp(self.f2,self.ene,(energy,))[0]


class _AtomicProperties:
  databaseA=specfile.Specfile("%s/AtomicConstants.dat" % dabaxpath  )
  databaseB=specfile.Specfile("%s/AtomicDensities.dat" % dabaxpath  )


  elementlistA=[]
  for scan in databaseA: 
     elementlistA.append( string.split(scan.command()) [-1])
  elementlistB=[]
  for scan in databaseB: 
     elementlistB.append( string.split(scan.command()) [-1])
  propA=string.split(databaseA[0].header('L')[0]  )[1:]
  propB=string.split(databaseA[0].header('L')[0]  )[1:]
  __doc__="I can tell you about %s  %s\n call function getProperty with two args: atom name and property name" %( string.join(propA), string.join(propB)) 

def AtomicPropertiesList():
    return _AtomicProperties.propA+_AtomicProperties.propB

def AtomicProperties(atom, property):
     if(property in _AtomicProperties.propA and atom in _AtomicProperties.elementlistA):
          return _AtomicProperties.databaseA[_AtomicProperties.elementlistA.index(atom)].datacol(property)   [0]
     elif(property in _AtomicProperties.propB and atom in _AtomicProperties.elementlistB):
          return _AtomicProperties.databaseB[_AtomicProperties.elementlistB.index(atom)].datacol(property)   [0]
     elif(property == "Z"):
          return  _AtomicProperties.elementlistA.index(atom)+1
     else:
       print _AtomicProperties.elementlistA
       print _AtomicProperties.elementlistB
       print _AtomicProperties.propA
       print _AtomicProperties.propB
       print " property %s not found "% property
       exit(0)

AtomicProperties.__doc__="I can tell you about %s  %s , Z\n call function  getProperty with two args: atom name and property name" %( string.join(_AtomicProperties.propA), string.join(_AtomicProperties.propB)) 



class FluorYield:
    """
        Atomic Radiative Yields for K and L Shells (M.O. Krause)

        This data set contains data calculated by using mucal (A subroutine to 
        calculate x-ray crossections) by Pathikrit Badyopadhyay.

        data is available for all the elements from 1 to 94.

        REFERENCE

        M.O. Krause "Atomic Radiative and Raduiattionless Yields for K and
        L Shells", J. Phys. Chem. Ref. Data, Vol. 8, No 2, 1979, pag 307-327

        The mucal program can be found at: 
            http://ixs.csrri.iit.edu/database/programs/mcmaster.html
        and a Periodic Table WWW interface is at:
            http://www.csrri.iit.edu/periodic-table.html

        Column description:
          1:FluorescenceY(K)
          2:FluorescenceY(L1)
          3:FluorescenceY(L2)
          4:FluorescenceY(L3)

        The scan #S 1000 contains the information for all the elements, where the first 
        column is the Atomic Number Z. This is useful to make plots of the values ofa 
        given variable as a function of Z.

        Author: O. Dhez (olivier.dhez@esrf.fr)

        USE:
            first initialize with a = FluorYield()
            Get value: a.Element(element, k or l1 or l2 or l3 or none)
                none give a tuple with k, l1, l2, l3
    """

    def __init__(self):
        print 'Start Fluor Yield'
        self.path = dabaxpath
        self.load('FluorYield_Krause.dat')    

    def load(self, location):
        self.Location=self.path+"/"+ location
        self.data=specfile.Specfile(self.Location)

        print "Opened file %s." %(self.Location)

        self.process()

    def process(self):
        self.elementlist=[]
        for scan in self.data: 
            self.elementlist.append( string.split(scan.command()) [-1])



    def Element(self, name, FY=None):
        try:
            index = self.elementlist.index(name)
        except ValueError:
            print "\n %s not found in %s\n" % (name, self.Location) 
            exit(1)
        temp = self.data[index].dataline(1)
        if (FY == 'k'):
            return temp[0]
        elif (FY == 'l1'):
            return temp[1]
        elif (FY == 'l2'):
            return temp[2]
        elif (FY == 'l3'):
            return temp[3]
        elif (FY == None):
            return (temp[0], temp[1],temp[2],temp[3])
        else:
            print '%s not found'%FY
            exit(0)



class XrayEmission:
    """
        XREmission.dat 
         X-Ray Emission Lines, K-level and L-level emission lines in KeV
         
          Values are from J. A. Bearden, "X-Ray Wavelengths", Review of Modern
          Physics, (January 1967) pp. 86-99, unless otherwise noted.
          
        Author: O. Dhez (olivier.dhez@esrf.fr)

        USE:
            first initialize with a = XrayEmission()
            Get value: 
                    a.EmiNRJ(element, level or none)
                        none give a tuple with with all the levels
            Name of the level:
                    a.ElementLevelList(element)
    """

    def __init__(self):
        self.path=dabaxpath

        self.load('XREmission.dat')    

    def load(self, location):
        self.Location=self.path+"/"+ location
        self.data=specfile.Specfile(self.Location)

        print "Opened file %s." %(self.Location)

        self.process()

    def process(self):
        self.elementlist=[]
        for scan in self.data: 
            self.elementlist.append( string.split(scan.command()) [-1])



    def EmiNRJ(self, name, level=None):
        try:
            index_element  = self.elementlist.index(name)
        except ValueError:
            print "\n %s not found in %s\n" % (name, self.Location)
            exit(1)
            
        level_label=self.data[index_element].alllabels()
        if (level != None):
            try:
                level_label.index(level)
            except ValueError:
                print "\n The %s emission line is not found for %s element\n" % (level, name)
                print " The exiting emission lines are %s for this element\n"%level_label
                exit(1)
                
        if (level):
            return float(self.data[index_element].datacol(level)[0])
        elif (level == None):
            return tuple(self.data[index_element].dataline(1))

    def ElementLevelList(self, name):
        index_element  = self.elementlist.index(name)
        if(index_element==-1):
            print "%s not found in %s" % (name, self.Location)
            exit(1)
        return self.data[index_element].alllabels()

class XrayEmissionWeights:
    """
        Weights for X-Ray Emission Lines 
        Sheet_Stand 
         This file belongs to the DABAX library. More information on
         DABAX can be found at:
         http://www.esrf.fr/computing/scientific/dabax/
        Weights for X-Ray Emission Lines

        Weights for fluorescence lines.
        Reference: Unknown origin. Data from a file given by P. Mangiagalli (ESRF)
 
        Row Titles: Siegbahn description of the transition (e.g., Ka1,Lb1)
         and IUPAC description (between brackets (e.g., K-L3,L2-M4)
        Row Data: Weights 

        Notes: 
        1) The corresponding IUPAC transition notation is under the #UIUPAC keyword
        2) The corresponding Siegbahn notation is under the #USIEGBAHN keyword
        3) The corresponding line energies [keV] are under the #UENERGY keyword

        Author: O. Dhez (olivier.dhez@esrf.fr)

        USE:
            first initialize with a = XrayEmissionWeights()
            Get value: 
                    a.EmiWeights(element, level or none)
                        none give a tuple with with all the levels
                    a.IUPAC(name)
                        give the IUPAC transition levels notation name for the given element
                    a.SIEGBAHN(name)
                        give the SIEGBAHN transition levels notation name for the given element
    """

    def __init__(self):
        self.path=dabaxpath

        self.load('XREmissionWeights.dat')    

    def load(self, location):
        self.Location=self.path+"/"+ location
        self.data=specfile.Specfile(self.Location)

        print "Opened file %s." %(self.Location)

        self.process()

    def process(self):
        self.elementlist=[]
        for scan in self.data: 
            self.elementlist.append( string.split(scan.command()) [-1])



    def Element(self, name):
        try:
            index_element  = self.elementlist.index(name)
        except ValueError:
            print "\n %s not found in %s\n" % (name, self.Location)
            exit(1)
        return index_element
        
    def EmiWeights(self, name, level=None):
        try:
            index_element  = self.elementlist.index(name)
        except ValueError:
            print "\n %s not found in %s\n" % (name, self.Location)
            exit(1)
            
        level_label=self.SIEGBAHN(name)
        if (level != None):
            try:
                index_label=level_label.index(level)
            except ValueError:
                print "\n The %s emission line is not found for %s element\n" % (level, name)
                print " The exiting emission lines are %s for this element\n"%level_label
                exit(1)
#            level_name = level_label[index_label]+'('+self.IUPAC(name)[index_label]+')'
            return self.data[index_element].datacol(index_label+1)[0]
        elif (level == None):
            return tuple(self.data[index_element].dataline(1))
        
    def IUPAC(self, name):
        try:
            index_element  = self.elementlist.index(name)
        except ValueError:
            print "\n %s not found in %s\n" % (name, self.Location)
            exit(1)
        iupac_list =[]
        iupac_list.append(string.split(self.data[index_element].header('UIUPAC')[0]))
        return iupac_list[0][1:]

    def SIEGBAHN(self, name):
        try:
            index_element  = self.elementlist.index(name)
        except ValueError:
            print "\n %s not found in %s\n" % (name, self.Location)
            exit(1)
        seig_list =[]
        seig_list.append(string.split(self.data[index_element].header('USIEGBAHN')[0]))
        return seig_list[0][1:]




class IndexFromObject:
    """ filters an index object  adding a multiplicative
        factor for the density. Handy to filter object like KK 
        It has the method index(self, wavelenghts)
    """  


    is_a_tensor = 0
    is_an_index = 0 

    def __init__(self,   Object=None,   RelativeDensity=1.0  ):

        """
          Gli argomenti sono: 
              Object: index object

        """
        self.data  = Object

        if( type(Object)==type([]) ):
          self.DabaxList=Object
          self.DensityList=RelativeDensity
          self.MassList=["does not matter"]*len(Object)
        else:        
          self.DabaxList=[Object]
          self.DensityList=[RelativeDensity]
          self.MassList=["does not matter"]



    def index(self, wavelenghts):
      indexes =  CalculateIndexesFromTable(self ,wavelenghts ) 
      return indexes



class IndexFromObjects:
    """ filters an index object  adding a multiplicative
        factor for the density. Handy to filter object like KK 
        It has the method index(self, wavelenghts)
    """  


    is_a_tensor = 0
    is_an_index = 0

    def __init__(self,    *args ):

        """
          Gli argomenti sono: 
             

        """
        nel=len(args)
        if( nel%2 !=0):
          raise " Wrong Number of arguments in IndexFromTable"
        nel=nel/2
        Object=[]

        self.DabaxList=[]
        self.DensityList=[]
        self.MassList=[]
        
        for i in range(nel):

          self.DabaxList.append(args[2*i])
          self.DensityList.append(args[2*i+1])
          self.MassList.append( "does not matter")

        self.data=self.DabaxList


    def index(self, wavelenghts):
      indexes =  CalculateIndexesFromTable(self ,wavelenghts ) 
      return indexes






class IndexFromFile:
    """ Reads index data from one file. 
        It has the method index(self, wavelenghts)
    """  


    is_a_tensor = 0
    is_an_index = 1

    def __init__(self,   NameFile =None,   RelativeDensity=1.0  ):

        """
          Gli argomenti sono: 
              NameFile: nome del file per gli indici . (Energy, Real(n), Imm(n))

        """
        self.data  = Dabax_IndexReader( NameFile  )


        self.DabaxList=[self]
        self.DensityList=[RelativeDensity]
        self.MassList=["does not matter"]



    def index(self, wavelenghts):

      nr = arrayfns.interp(self.data[:,1], self.data[:,0], 12398.52/ wavelenghts )
      ni = arrayfns.interp(self.data[:,2], self.data[:,0], 12398.52/ wavelenghts  )


      indexes=nr+complex(0.0,1.0)*ni

      return indexes



def Dabax_IndexReader(filename, Np="automatic"  ):

   """ given an optical index filename in input it reads it.

       If Np is 'automatic' the number of points will be determined
          automatically.
       If it is an integer, that will be. Finally if it is 'first line'
         it will be read from the first line.

   """  

   f=open(filename,"r")
   datas=f.read()
   datalines=string.split(datas,"\n")

   if(Np=="automatic"):
       Np=len(datalines)

   elif(Np=="first line"):
       Np=string.atoi(datalines[0])
       datalines=datalines[1:Np+1]
   elif( isinstance(Np,type(1) ) ):
       pass
   else:
       raise " PROBLEM with Np in IndexReader\n"

   data=map(string.split,datalines)

   dim=len(data[0])
   for i in range(Np):
     data[i]=map(string.atof,data[i])

   newdata=[]
   for i in range(Np):
     if(len(data[i])==dim):
         newdata.append(data[i])
   data=Numeric.array(newdata)


   return data













#####
#### BETAMANIPULATORS

""" The list of scan may be composed of a list of objects like this one :
   [ energies, values, [ (weighta, pola),(weightb, polb)....   ], weightscan]
   where energies and values are arrays,  weigth are numbers, and 
   and pol is an integer going from 1 to 3, indicating the polarisation.

  Weigth scan is the statistical wight to be given to the scan
  when doing fits.

   return [res.energies, res.data, Polarisation, res.weight]

"""

class return_value:
  pass


#########################################################


def BetaManipulator_ScanReader(filename=None, Np="automatic", Polarisation=None, energies_col=1, data_col=None, weight_col=None, applyPol="applyPol"):
   """ given a filename in input it reads the scan.

       If Np is 'automatic' the number of points will be determined
          automatically.
       If it is an integer, that will be. Finally if it is 'first line'
         it will be read from the first line.

       The other entries, if integer, will tell the column to read for
       that property. Otherwise one can specify a float number, and that will be.

       The return of this function is :

           return [res.energies, res.data, Polarisation, res.weight]   
   """  

   if(filename is not None  and type(filename)==type("s") ):

     synthetic = 0

     f=open(filename,"r")
     datas=f.read()
     datalines=string.split(datas,"\n")

     if(Np=="automatic"):
         Np=len(datalines)
     elif(Np=="first line"):
         Np=string.atoi(datalines[0])
         datalines=datalines[1:Np+1]
     elif( isinstance(Np,type(1) ) ):
         pass
     else:
         raise " PROBLEM with Np in ScanReader\n"

     data=map(string.split,datalines)
     dim=len(data[0])
     for i in range(Np):
       data[i]=map(string.atof,data[i])
     # print data
     newdata=[]
     for i in range(Np):
       if(len(data[i])==dim):
           newdata.append(data[i])
     data=Numeric.array(newdata)

   else:
      newdata=[]
      columns=[]

      synthetic = 1

      for  (item_col, name_col) in [ (energies_col, "energies" ) ]:
        if( type(item_col)==type([])):
          newdata.append(      Numeric.arrayrange( item_col[0], item_col[1], item_col[2]   )       )
          exec( "%s=len(newdata)"% name_col  )

      if( filename is not None):
          betaobject = filename
          betadata = 0
          if(applyPol=="applyPol"):
            for polFP in Polarisation:
              fact = polFP[0]
              pol  = polFP[1]
              betadata = betadata +  betaobject.contribution(newdata[0],pol) * fact
          else:
            betadata = betadata +  betaobject.contribution(newdata[0]) * fact

          newdata.append(    betadata       )
          data_col=len(newdata)
      else:
        if( weight_col is not None):
          raise " weight_col should be None"
        if( data_col is not None):
          raise " refle_col should be None"

      data=Numeric.array(newdata)


      data=Numeric.transpose(data)



   res=return_value()

   for (colonna_n, colonna_dati) in [ (energies_col, "energies" ),
                                     (data_col, "data" ),(weight_col, "weight" ) ]:

    if( isinstance(colonna_n , type(1))):
      if( colonna_n> Numeric.shape(data)[1]):
          message= " PROBLEM with %s_col > dimension \n"% colonna_dati
          raise message
      setattr(res,colonna_dati, data[:,colonna_n-1]*1.0)
    elif( isinstance( colonna_n , type(1.2))):
      setattr(res,colonna_dati, data[:,0]*0.0+colonna_n)
    elif( colonna_n is None):
      setattr(res,colonna_dati, None)
    elif(colonna_dati=="energies" and synthetic ):
      setattr(res,colonna_dati, data[:,0])
    else:
      message= " PROBLEM with %s_col type\n"% colonna_dati
      raise " PROBLEM with wavelenghts_col type\n"

   return [res.energies, res.data, Polarisation, res.weight]





class BetaManipulator_BetaJoin:
    """
        joins an absorption spectra
        to tabulated optical index ( imaginary part, beta)
        to obtained a beta on the wider range given by
        the tabulated values, but having more precised
        data, given by betaobject, in the shorter range given by betaobject
    """
    def __init__(self,  betaobject, material, trim_left=0, trim_right=0):
       self.betaobject = betaobject
       self.material = material
       self.trim_left  = trim_left
       self.trim_right = trim_right


    def BetaJoined_(self,energies):
        # calculate the index from material
        material=self.material
        betaobject = self.betaobject
        if( type(betaobject)!=type([])):
             betaobject=[betaobject]
        N = len(energies)
        indexes_Table=Numeric.zeros( N , "d")

        wavelenghts_Table = 12398.52 / energies


        for i in range(0,len(self.material.DabaxList)):
           scatterer= material.DabaxList[i]
           density  = par(material.DensityList[i])
           mass     = material.MassList[i]


           # print " Scatterer est ", scatterer

           add=-(415.22*(  scatterer.F_Lambda(wavelenghts_Table) )*( 
                             wavelenghts_Table*wavelenghts_Table/12398.52/12398.52/mass)).imag

           indexes_Table=indexes_Table+ add*density


        result =   indexes_Table
        
        for bobj in betaobject:
             index_object = bobj.contribution(energies)  
        

             min,max = bobj.extrema()

             min = min + self.trim_left
             max = max - self.trim_right




             diff1 = arrayfns.interp( indexes_Table- index_object, energies, [min]) [0]
             diff2 = arrayfns.interp( indexes_Table- index_object, energies, [max]) [0] 
             baseline = diff1 + (energies-min)*((diff2-diff1)/(max-min))
             index_object = index_object + baseline

             result = (result * (Numeric.less(energies,min) + Numeric.less(max,energies) ) + 
                       index_object * (Numeric.less_equal(energies,max) * Numeric.less_equal(min,energies) ) )

        return energies, result

    def contribution(self,  energies, pol=1, dichroism=0):
        return     self.BetaJoined_(energies)[1]









class BetaManipulator_BetaFilter:
    """
        Allows to access multicolumn beta object ( polarisation dependent )
        with a pre-selected polarisation
    """
    def __init__(self,  betaobject, Polarisation = None):
       self.betaobject = betaobject
       self.Polarisation = Polarisation


    def contribution(self,  energies, pol=1, dichroism=0):
        result = 0
#        for polFP in self.Polarisation:
#           result = result + self.betaobject.contribution(self,  energies, polFP[0] ) * polFP[1]

        result = result + self.betaobject.contribution(self,  energies, self.Polarisation ) 

        return   result

class bidon:
  def __init__(self,X):
    self.X=X

class BetaManipulator_ContributionsFromFile:
  """
      reads a file containing up to three polarisations.
      Initialised with shift, factor, saturation for 
      dichroism ( they can be variables ).

      It provides the function contributions(energies, pol)


      Rescale X lambda is set to one when one has raw data
      for absorption. In that case beta is proportional to absorption
      time lambda/lambda0 where lambda0 is the
      middle of the scan.

      The first column of filename is the energy 

  """
  def __init__( self, filename, shift, factor,  rescaleXlambda=1, saturation=1, Three_Files=0, File2=None, File3=None ):

     self.charge_correction=0
     self.data = Dabax_IndexReader(filename, Np="automatic"  )
     self.data=Numeric.transpose(self.data)
     self.Three_Files = Three_Files
     if(self.Three_Files):
       self.data2 = Dabax_IndexReader(File2, Np="automatic"  )
       self.data3 = Dabax_IndexReader(File3, Np="automatic"  )

       self.data2=Numeric.transpose(self.data2)
       self.data3=Numeric.transpose(self.data3)


       self.datas_=[ self.data, self.data2, self.data3 ]

     self.shift = shift
     self.factor = factor
     self.saturation = saturation

     if(rescaleXlambda==1) :
       ene0 = 0.5*(self.data[0][0]+self.data[0][-1])
       self.data[1:,:] = self.data[1:,:] * ene0 / self.data[0,:]
       if(self.Three_Files):
         self.data2[1:,:] = self.data2[1:,:] * ene0 / self.data2[0,:]
         self.data3[1:,:] = self.data3[1:,:] * ene0 / self.data3[0,:]

     self.maglia =[ bidon(X) for X in  self.data[0] ]
  def extrema(self):
   """
       gives the energy minimum and maximum between which
       data are available
   """
   shift = par(self.shift)
   if(self.Three_Files):
      return (max(self.data[0][0]+shift, self.data2[0][0]+shift, self.data3[0][0]+shift),
              min(self.data[0][-1]+shift, self.data2[0][-1]+shift, self.data3[0][-1]+shift ) )   
   else:
      return (self.data[0][0]+shift, self.data[0][-1]+shift)    


  def charge_correct(data):
    pass


  def contribution(self,  energies, pol=1, dichroism=0):
     """
         Fait un'interpolation pour donner
         le beta en fonction des energies.
         Energies must be in the same units as the
         first column of the file given in the
         initialisation.

         pol is the number of the column to retrieve.
         However, if dichroism==1, the saturation factor
         is used :
              the filename is supposed to contain 4 columns
              ( energie D+ Dz D-) and the difference between
              D+ and D- will be rescaled by saturation
     """
     shift = par(self.shift)
     factor= par(self.factor)

     for i in range(1, len(self.data[0]) ):
          if self.data[0][i]==self.data[0][i-1]:
               print self.data[0][i]
               raise " uguali"

     if(self.Three_Files):
       result =  arrayfns.interp( self.datas_[pol-1][1]   , self.datas_[pol-1][0]+shift ,energies)

     else:
      if(len(self.data)>pol):
        result =  arrayfns.interp( self.data[pol]   , self.data[0]+shift ,energies) 
      else:
        result =  arrayfns.interp( self.data[1]   , self.data[0]+shift ,energies) 

     if(self.charge_correction):
        result=self.charge_correct(result)

     saturation = par(self.saturation)

     # if( dichroism==1 and pol in [1,3] ) :
     if( ( self.Three_Files  or  len(self.data)==4 ) and pol in [1,3] ) :
        # saturation = par(saturation)
        if(self.Three_Files):
          result2 = arrayfns.interp( self.datas_[4-pol-1][1]   , self.datas_[4-pol-1][0]+shift ,energies)            
        else:
          result2 = arrayfns.interp( self.data[4-pol]   , self.data[0]+shift ,energies)  

        if(self.charge_correction):
           result2=self.charge_correct(result2)

        result = ((1+saturation)*result + ( 1-saturation)*result2)/2
     
     return factor*result

## p[2]*atan(p[0]*(x-p[1]))+p[3]+(p[4]*(p[5]/2)^2 /((x-p[6])^2+(p[7]/2)^2))^2
## +(p[8]*(p[9]/2)^2 /((x-p[10])^2+(p[11]/2)^2))^2



class BetaManipulator_ContributionsFromContinuum:
     """
           A function specified by :
              E0
              step
              pente
           it is 0 for E<E0

           step+pente*(E-E0) 

     """
     def __init__(self, E0=None, step=None, pente=0,arctanfact=None,  min = -1.0e20 , max =  1.0e20):
         self.E0=E0
         self.step=step
         self.pente=pente
         self.min = min
         self.max = max
         self.arctanfact =  arctanfact

     def contribution(self,  energies, *args):
         """
         """
         step  = par(self.step)
         E0    = par(self.E0)
         pente = par(self.pente)

         result = step + (energies - E0 ) * pente
         if( self.arctanfact is None ) :
           result =result*Numeric.less(E0, energies)
         else:
           result =result*(0.5+Numeric.arctan( ( energies-E0)* par(self.arctanfact) ) / Numeric.pi)

         return result

     def extrema(self):
        """
          gives the energy minimum and maximum between which
          data are available
        """
        return ( self.min , self.max )    



class BetaManipulator_ContributionsFromInterpolation:
     """
           A function specified by :
             - a Beta object for initialisation
             - one or more array of variables like those created by 
                 CreateVariableArray.
           Each variable in the arrays has an energy (X in CreateVariableArray)
           an initial value ( generally 1 ) a min ( generally 1-tol)
           and a max.

           The initialisation process converts 1 to beta(X) etc. etc.
     """
     def __init__(self, BetaObject=None, shift=0, arrs=[], nameForView="x_y_beta_interp.dat"):
        self.arr=[]
        self.shift=shift
        self.X=[]
        self.nameForView=nameForView
        if(type(arrs[0])!=type([])):
             arrs=[arrs]
        for arr in arrs:
          for v in arr:

           X=v.X
           self.X.append(X)

           b = BetaObject.contribution([X])[0]
           v.value = v.value* b
           v.min   = v.min  * b 
           v.max   = v.max  * b 

           self.arr.append(v)
           # print v
           # print self.X
           # print v.value

        # raise " OK " 
        self.X=Numeric.array(self.X)

     def contribution(self,  energies, *args):

         shift = par(self.shift)
         # print " mappo par su ", self.arr
         # raise " OK"
         Y=map(par, self.arr)

         if( self.nameForView is not None ) :
           f=open(self.nameForView,"w")
           for i in range(len(self.X)):
             f.write("%e %e \n"%(self.X[i],Y[i]) )
           f.close()
           
         result =  arrayfns.interp( Y   , self.X ,energies) 

         return result

     def extrema(self):
        """
          gives the energy minimum and maximum between which
          data are available
        """
        return ( self.X[0] , self.X[-1] )    



class BetaManipulator_ContributionsFromLorentz:
     """
           A function specified by :
              E0
              gamma
              height
     """
     def __init__(self, E0=None,height=None , gammaL=None, gammaR=None, min = -1.0e20 , max =  1.0e20 ):

         self.E0=E0
         self.height=height
         self.gammaL=gammaL
         self.gammaR=gammaR
         self.min = min
         self.max = max

     def contribution(self,  energies, *args):
         """
         """
         height  = par(self.height)
         E0    = par(self.E0)
         gammaL = par(self.gammaL)
         gammaR = par(self.gammaR)
         mask = Numeric.less(energies, E0)
         gamma2 = (gammaL*gammaL)*mask + (gammaR*gammaR)*(1-mask)
         l1= gamma2 / ( (energies - E0 ) *(energies - E0 ) + gamma2  )
         result = height *l1 * l1

         return result

     def extrema(self):
        """
          gives the energy minimum and maximum between which
          data are available
        """
        return (  self.min , self.max   )    


class BetaManipulator_SumOfContributions:
  """
       just memorizes a series of contribution object
       Provides a contribution function that returns
       the sum of all contributions
  """
  def __init__(self, *args):
    self.contribss = args

  def contribution(self,  energies, pol=1, dichroism=0):
     result=0
     if( len(self.contribss)%2 == 0  and not
               hasattr(self.contribss[1],"contribution")):
      for i in range(len(self.contribss)/2):
       tok=self.contribss[i*2]
       peso=par(self.contribss[i*2+1])
       result=result + peso*tok.contribution(energies,pol,dichroism)
     else:
      for tok in self.contribss:
       # print type(tok)
       # print tok.__class__
       result=result + tok.contribution(energies,pol,dichroism)
     return result

  def extrema(self):
        """
          gives the energy minimum and maximum between which
          data are available
        """
        min_=1.0e-20
        max_=1.0e+20
        for c in self.contribss:
         if(hasattr(c,"contribution")):
          cmin,cmax = c.extrema()
          min_=max(min_,cmin)
          max_=min(max_,cmax)

        return (min_,max_)        

class BetaManipulator_Comparison:

  """ this object is aware of a stack model and of a set of experimental data,
      (eventually void data, where just the scan are specified)

  """
  classe="BetaManipulator_Comparison"
  def __init__(self,stack, scanlist, weightlist, normlist=None, printpartial=0, meritfunction="diffroot"):

      self.meritfunction=meritfunction
      self.stack=stack
      self.scanlist=scanlist
      self.weightlist=weightlist
      self.printpartial=printpartial

      if( normlist==None):
         self.normlist=[1,]*len(scanlist)
      else:
         self.normlist=normlist
      self.passes=0
      self.itercounter=0

  def setitercounter(self,n):
      self.itercounter=n

  def getitercounter(self):
      return self.itercounter

  def write2File(self, filename):
     """ The function error, calculates the error,
         and stores the model scan data in the variable
         self.calculatedscan

     """ 
     for count in range(len(self.calculatedscan ) ) :
       calculated=self.calculatedscan[count]
       scan = self.scanlist[count]
       f=open("%s%d"%(filename, count+1),"w")
       for k in range(len(calculated)):
         f.write("%e %e %e\n"%(calculated[k], scan[0][k], scan[1][k]  ))
       f.close()




  def error(self, nopartial=0):
     """ For each scan calculated the experimental observable 
         as given from the model.

         The merit function is calculated :
            -- if abslogdiff is choosed as the sum of the asolutes 
            of the differences of the logarithmes multiplied by the wheigths

            -- If diffroot is choosed in the usual way. the weigths
               multiplies the squares

         if nopartial=0 the partial results are saved on files using function write2File
         and name= partialresults.

     """
     error=0
     count=0

     csc=0
     numtot=0

     # '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
     #   stores  the calculated values for eventual retrieval in writefit2file
     #   
     # '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
     if( hasattr(self, "calculatedscan")==0): self.calculatedscan=[]
     if( hasattr(self, "secondmodel"   )==0): self.secondmodel=[]
     newscans=[]
     if( type(self.stack)==type([]) ):
        stack=self.stack[0]
        stack1 = self.stack[1]
     else:
        stack  = self.stack
        stack1 = None


     new_secondmodel=[]

     for scan in self.scanlist:
        energies = scan[0]
        data=scan[1]

        polar=scan[2]

        if(len(scan)==4):
           errorweights=scan[3]
        else:
           errorweights=Numeric.ones(len(data ))

        rs=Numeric.zeros(len(energies),"d")

        for tok in polar:

          rsadd=0

          rsadd= stack.contribution(  energies  ,tok[1] )

          rs+=tok[0]*rsadd


        norm=par(self.normlist[count])
        newscans.append(rs*norm)


        if(data is None):
            if(stack1 is not None):
                 data=0
                 for tok in polar:
                     data=data+ tok[0]*stack1.contribution(  energies  ,tok[1] )
                 new_secondmodel.append(norm*data)

            else:
               data=norm*rs



        if(len(scan)==4 and  scan[3] is not None):
            errorweights=scan[3]
        else:
            errorweights=Numeric.ones(len(data ),'d')

        if(self.meritfunction=="abslogdiff"):
          difference=abs(log(norm*rs)-log(data))
          # difference=difference*difference*errorweights
          difference=difference*errorweights

        elif( self.meritfunction=="diffroot" ): 

          difference=  norm*rs-data
          difference=difference*difference*errorweights


        error=error+Numeric.sum(difference)*self.weightlist[count]
        numtot=numtot+len(difference)

     self.calculatedscan=newscans
     self.secondmodel   =new_secondmodel

     if(self.printpartial and nopartial==0):
         self.write2File("partialresults")

     self.itercounter=self.itercounter+1
     return error/numtot




################################################################################
class Element_ComptonMcMaster:
  """
	 Ffit[ene] = exp{ c0 + c1*[log(ene)]^1 + c2*[log(ene)]^2 + c3*[log(ene)]^3}
  """
  def __init__( self, lista  ):
     self.c=lista
  def CrossSection(self, Ene):

     Ene=1000.0*Ene

     lgs = Numeric.log(Ene)

     res = lista[0] + lgs * (  lista[1] + lgs * ( lista[2] + lgs*lista[3]  )  )

     return res


class CrossSecComptonMcMaster:
  """	
	#F  CrossSec-Compton_McMaster.dat
	#D Fri Sep 27 11:51:46 2002
	#UT Parametrization of incoherently scattered X-ray intensities vs energy.
	#C This file has been created using CrossSec-Compton_McMaster.pro on Fri Sep 27 11:51:46 2002
	#UD  
	#UD  Parametrization of incoherently scattered X-ray intensities vs energy.
	#UD  
	#UD  This file contains the tabulated coefficients for calculation
	#UD  of Compton Cross Section as a function of the photon energy.
	#UD  
	#UD  REFERENCE: 
	#UD    W.H.McMaster, N. Kerr Del Grande, J. H. Mallett and J. H. Hubbell,
	#UD    UCRL-50174 SEC II REV I, available from Nartional Tecnical Information
	#UD    Service, U.S. Dept. of Commerce.
	#UD  
	#UD  The analytical function is: 
	#UD    Ffit[ene] = exp{ c0 + c1*[log(ene)]^1 + c2*[log(ene)]^2 + c3*[log(ene)]^3}
	#UD 
	#UD  where ene is the Photon Energy [in keV] and c0-c3 are the tabulated coefficients.
	#UD 
	#UD  The data in this file has been extracted from the file raycomin.f 
	#UD  of the library by Brennan and Cowan (S. Brennan and P.L. Cowan (1992)
	#UD  Rev. Sci. Instrum. 63,1, 850) 
	#UD  These data contain parameters for for elements with 1<=Z<=92,
	#UD  thus, the data in this DABAX file goes from Z=1 to Z=92.
	#UD  The file raycomin.f is available by ftp from the authors from
	#UD  http://www-ssrl.slac.stanford.edu/absorb.html or
	#UD  ftp://ftpa.aps.anl.gov/pub/cross-section_codes/sb-1.0.tar.Z
	#UD 
	#UD  Column description:  c0 c1 c2 c3
	#UD Note: The coefficients in this file are slightly different from the ones in
	#UD       the file mucal.f which has been used for creating the DABAX file
	#UD       CrossSec_McMaster.dat. We believe that the present set from Brennan's
	#UD       package are more refined than the ones in mucal.f because its reference
	#UD       is more recent.
  """
	
  def __init__(self, Name=None):
    self.path=dabaxpath

    self.dataprocessors   ={   "CrossSec-Compton_McMaster": self.process_Compton, }

    self.elementretrievers={  "CrossSec-Compton_McMaster":   self.Element_Compton , }
    if(Name!=None):
      self.load(Name)    


  def getDataBaseNames(self):
    return self.elementretrievers.keys()

  def load(self, location):
     self.Location=self.path+"/"+ location
     self.fileType=location
     self.data=specfile.Specfile(self.Location)

     print "Opened file %s. It is of the  %s type"%(
                          self.Location,self.fileType)


     if( (self.fileType in self.dataprocessors.keys()) and
         (self.fileType in self.elementretrievers.keys()) ):

        self.Element= self.elementretrievers[self.fileType]
        self.dataprocessors[self.fileType]()

     else:
        print " I can't process a file of the type %s", self.fileType
        print " Possible types are ", self.fileType in dataprocessors.keys
        exit(1)


  def process_Compton(self):
     self.elementlist=[]
     for scan in self.data:
       self.elementlist.append( string.split(scan.command()) [-1])


  def Element_Compton(self,name):
      index=self.elementlist.index(name)
      if(index==-1):
         print "%s not found in %s" % ( name, self.Location)
         exit(1)
      lista=self.data[index].dataline(1)
      res= Element_ComptonMcMaster( lista  )
      res.element=name
      return res







#######################################################
# ESEMPI DI UTILIZZAZIONE
#
if (__name__=="__main__"):

  ## LOW LEVEL USE
  ## --------------------------------------------------------

  Table_f0= Dabax_f0_Table("f0_WaasKirf.dat")

  Table_f1f2=Dabax_f1f2_Table("f1f2_Windt.dat")


  Bef0  = Table_f0.Element("O2-")
  Bef1f2= Table_f1f2.Element("Be")

  print Bef0.f0Energy(1000.0)
  print Bef1f2.f1f2Energy(1000.0)


  Be=Dabax_Scatterer((Bef0,),(1,), (  Bef1f2,),(1,)       )

  print Be.F_Energy(1000)

  print "#######################"



  Wf0  = Table_f0.Element("W")
  Wf1f2= Table_f1f2.Element("W")

  print Wf0  .f0Energy(1000.0,0)
  print Wf1f2.f1f2Energy(1000.0)


  W=Dabax_Scatterer((Wf0,),(1,), (  Wf1f2,),(1,)       )

  print W.F_Energy(1000)

  print "#######################"

  print (1*W+Be*1).F_Energy( 2000  )         
  print (1*W+Be*1).F_Energy( Numeric.array((1000,2000,80000))         )


  ## HIGH LEVEL USE

  melange =  IndexFromTable(Table_f0, Table_f1f2, "Si", 2.0,"O", 1.0)
  print CalculateIndexesFromTable( melange , 10.0)


  ## LOW LEVEL USE

  table_f1f2 = Dabax_f1f2_Table()
  tabname_list = table_f1f2.elementretrievers.keys()
  tab_list=[]

  for name in tabname_list:
       tab_list.append(Dabax_f1f2_Table(name) )

  elname_list = tab_list[0].elementlist

  enelist = [15816., 17794.  , 21747. ]
  print " USING TABLES ", tabname_list
  print "  AT ENERGIES ",  enelist

  for elname in elname_list:
     print "%3s     " %  elname,
     rho = AtomicProperties(elname, "Density[g/ccm]")
     for ene in enelist :
         for tab in tab_list:
           if(elname in tab.elementlist ):
             substance = IndexFromTable ( Table_f0, tab, elname, rho)
             beta = CalculateIndexesFromTable(substance , 12398.52/ene).imag
             mu   = 2*(beta *2*math.pi *ene/12398.52 ) 
             print "%6.4e " % mu,
           else:
             print "------     ",
         print "     ",
     print ""          





#
#
#######################################################
