################################################################################
# Alessandro Mirone
# Aprile 2001
#  ESRF

#########################
##Disclaimer
##==========
##
## This  software is provided without warranty of any kind.
## No liability is taken for any loss or damages, direct or
## indirect, that may result through the use of it. No warranty
## is made with respect to this documantation, or the programs
## and functions therein.
## There are no warranties that the programs or their documentation
## are free of error, or that it is consistent with any standard,
## or that it will meet the requirement for a particular application.

## Copyright
## =========
##
## We have adopted the GNU LIBRARY GENERAL PUBLIC LICENSE to apply to
## the this software.
## For more information on the license terms, see
## http://www.gnu.org/copyleft/lgpl.txt

print " PPMlayerclass.py DICROICO VERSIONE 2004-4-19 "

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



import Numeric




from Numeric import *
import PPMcore
from Dabax import *
from  arrayfns  import *            
import string
import Minimiser
import Numeric
import MLab
import LinearAlgebra
import arrayfns

TR=Numeric.transpose

def   shift_add(convoluted,rs,shift,fact):

  """ used to calculate convolutions :
      an array rs is shifted by a given amount (shift)
      and summed to the result ( convoluted)
  """

  if(shift>0):
   insect=convoluted[shift:]
   add(insect,  fact*rs[:-shift],  insect)
  elif(shift<0):
   insect=convoluted[:shift]
   add(insect,  fact*rs[-shift:],  insect)
  else:
   add(convoluted,  fact*rs,  convoluted)


def par(x):
  """ Used to manipulated parameters that could
      be either numerics or of the variable type.
      Variable can be any variable object that can be 
      modified by an external entity ( for example
      a minimisator object ).
      In this case ( if they are of the variable type)
      they have to provide the getvalue function.

      Any object having the getvalue function can be treated
      as a variable
  """
  if('getvalue' in dir(x)):
    return x.getvalue()
  else:
    return x

def connect(a,b):
   """ stackable objects, stacks, layers..
       are piled one on the top of the other
       by pointing next member to the upper 
       contigous object, and before member
       to the previous object.

       next goes to the top, and before 
       points toward the bottom ( substrate)
   """

   a.next=b
   b.before=a

def SumThings(*args):

  """ connects things togheter
      The Plus operator ( + )
      is redefined in the class PPM_stack ( see below)
      and uses the  connect function
  """

  res=PPM_Stack()
  for tok in args:
    res=res+tok
  return res

def MultiplyAThing(N, tok):


  """ sums tok N times to itself to form 
      a repeated periodic  structure
  """

  res=PPM_Stack()
  for i  in range(N):
    # print "Moltiplicazione ", i
    res=res+tok
  return res



def Reflectivity( res, propagazione,   interfaccia):
       """ recursive formula giving the reflectivity
           at an upper interface as function of the reflectivity
           at the lower interface, of the propagation
           through the intermediate material and of the
           upper interface
       """

       propagazione.shape=( propagazione.shape[0], 1,4)

       propa = interfaccia * propagazione

       a = propa[:,0:2,0:2]
       b = propa[:,0:2,2:4]
       c = propa[:,2:4,0:2]
       d = propa[:,2:4,2:4]

       fast = 1

       res= Numeric.array(res)

       if(fast==0):
  
         aR = Numeric.array( map( Numeric.dot, a, res) ) + b
         cR = Numeric.array( map( Numeric.dot, c, res) ) + d       
         # res = Numeric.array(     map( Numeric.dot ,   aR, map(   LinearAlgebra.inverse, cR)  )  )
         res = (     map( Numeric.dot ,   aR, map(   LinearAlgebra.inverse, cR)  )  )
       else:
         aR  = PPMcore.PPM_MAPDOT( a, res)  + b
         cR  = PPMcore.PPM_MAPDOT( c, res)  + d       
         # res = (     map( Numeric.dot ,   aR, map(   LinearAlgebra.inverse, cR)  )  )
         res = PPMcore.PPM_MAPDOT( aR, PPMcore.PPM_MAPINV(cR) )
       return res
    

class PPM_Stack:
  """ Represent a multilayer as a sequence of layers
      connected by next and before pointers ( functions connect manipulate them)
  """ 
  def __init__(self):
    self.stack=[]
    self.classe="PPM_Stack"

  def __add__(self, other):
    newstack=[]

    if( 'classe' not in dir(other)):
       print "Error summing a wrong class "
       return 0

    if(other.classe  ==   'PPM_Stack'):
       if ( len(self.stack)):
         connect( self.stack[-1]    , other.stack[0]       ) 
         connect( other.stack[-1]   , self.stack[0],     ) 
       newstack=self.stack+other.stack

    elif(other.classe  in [  'PPM_LayerBaseClass'] ):
       if ( len(self.stack)):
         connect( self.stack[-1], other      ) 
         connect(other,  self.stack[0]       ) 
       newstack=self.stack+[other]

    else:
       print "Error summing a wrong class "
       return 0

    res=PPM_Stack()
    res.stack=newstack
    return res


  def calculatescanTensor(self,wavelenghts, angles, Polarisation):

     """ calculate two scan, one as usual and the
         other by reflecting the theta with respect to the normal
         to check the dichroism
     """

     res=[0,0]
     res[0]= self.calculatescanTensor_(wavelenghts, angles, Polarisation)
     res[1]= self.calculatescanTensor_(wavelenghts, math.pi-angles, Polarisation)
     return res

  def calculatescanTensorParallel(self,wavelenghts, angles, Polarisation, nscan):

     """ calculate two scan, one as usual and the
         other by reflecting the theta with respect to the normal
         to check the dichroism
     """

     res=[0,0]
     res[0]= self.calculatescanTensor_Parallel(wavelenghts, angles, Polarisation, nscan)
     res[1]= self.calculatescanTensor_Parallel(wavelenghts, math.pi-angles, Polarisation, nscan)
     return res




  def calculatescanTensor_Parallel(self,wavelenghts, angles, Polarisation, nscan):
      if( (  mpi.size % nscan  ) != 0 ):
          print " lenscans = ",nscan , "  mpi.size = ", mpi.size
          raise " PROBLEM (mpi.size  % nscan) != 0  in calculatescanTensor_Parallel "

      nprocs=  mpi.size/ nscan
      
      nmaster = mpi.rank % nscan

      nself   = mpi.rank / nscan

      lscan = len(wavelenghts)




      res=Numeric.zeros(lscan,"d")

      lentoken = int(lscan*1.0/ nprocs  + 0.999999999)

      inizio = nself*lentoken
      fine   = (nself+1)*lentoken
      if( fine > lscan ):
          fine=lscan

      Thickness, Roughness,indexlist,SubIndex,SubRough = self.getDescription(wavelenghts[inizio:fine])

      eigenstuff=PPMcore.PPM_calculatescanTensor(
        Thickness[1:],Roughness[1:],
        wavelenghts[inizio:fine],
        angles[inizio:fine],indexlist[1:],
        SubIndex,SubRough)

      K0=2*pi/wavelenghts[inizio:fine]

      propagazioni = Propagations( eigenstuff, Thickness,K0)


    

      print " calcolo interfacce "
      interfacce= Interfaces(eigenstuff,Roughness ,K0)
      print " calcolo interfacce OK "
      
      NI=len(interfacce)
      NW=len(wavelenghts[inizio:fine])
      
      ###############################################################
      # reflectivity on the substrate
      #
      
      inter = interfacce[0]
      fast=1
      if(fast):
        print " OK  ? "
        reslocal = PPMcore.PPM_MAPDOT( inter[:,0:2,2:4 ]  , PPMcore.PPM_MAPINV(inter[:,2:4,2:4 ])  )    
      else:
        reslocal = Numeric.array(map( Numeric.dot, inter[:,0:2,2:4 ]  ,map(LinearAlgebra.inverse, inter[:,2:4,2:4 ])  )   ) 
      print " OK"
      ###################################################################
      # calculation of reflectivity at every interface
      #
 
      print " calcolo per " , NI , " interfacce "
      for i in range(1,NI):
        reslocal = Reflectivity( reslocal, propagazioni[i],   interfacce[i])

      print " la shape di reslocal  est ", reslocal.shape
      pol = ones([ len(reslocal)  ,1])*Polarisation
      # reslocal=map(Numeric.dot, reslocal, pol)
      print " la shape di  pol  est ", pol.shape
      reslocal = PPMcore.PPM_MAPDOT_V( reslocal, pol)
      reslocal=TR(  sum(TR( reslocal*Numeric.conjugate(reslocal)        )))
      reslocal=reslocal.real

      print " la shape di reslocal est ", reslocal.shape
      # COMMUNICATION
      for i in range(nprocs):
          iproc = i*nscan+nmaster
          if( mpi.rank== nmaster):
              if(i==0):
                  receivedS = reslocal
              else:
                  receivedS, status = mpi.recv(iproc , mpi.ANY_TAG  )                  
              inizio = i*lentoken
              fine   = (i+1)*lentoken
              if( fine > lscan ):
                  fine=lscan
              res[inizio:fine]= receivedS
          else:
              if( i == nself ):
                  mpi.send(reslocal, nmaster) 
      return  res
  




  def calculatescanTensor_(self,wavelenghts, angles, Polarisation):

    """ reflectivity calculated by recursive formula
    """

    Thickness, Roughness,indexlist,SubIndex,SubRough = self.getDescription(wavelenghts)


    print " dopo il substrato ci sono ", len(Thickness[1:]), " layers "
    eigenstuff=PPMcore.PPM_calculatescanTensor(
				Thickness[1:],Roughness[1:],
				wavelenghts,
                                angles,indexlist[1:],
				SubIndex,SubRough)
   
    K0=2*pi/wavelenghts


    propagazioni = Propagations( eigenstuff, Thickness,K0)


    

    print " calcolo interfacce "
    interfacce= Interfaces(eigenstuff,Roughness ,K0)
    print " calcolo interfacce OK "
	
    NI=len(interfacce)
    NW=len(wavelenghts)

    ###############################################################
    # reflectivity on the substrate
    #

    inter = interfacce[0]
    fast=1
    if(fast):
      print " OK  ? "
      res = PPMcore.PPM_MAPDOT( inter[:,0:2,2:4 ]  , PPMcore.PPM_MAPINV(inter[:,2:4,2:4 ])  )    
    else:
      res = Numeric.array(map( Numeric.dot, inter[:,0:2,2:4 ]  ,map(LinearAlgebra.inverse, inter[:,2:4,2:4 ])  )   ) 
    print " OK"
    ###################################################################
    # calculation of reflectivity at every interface
    #
 
    print " calcolo per " , NI , " interfacce "
    for i in range(1,NI):
       res = Reflectivity( res, propagazioni[i],   interfacce[i])

    print " PASSO DI LI "
    # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
    pol = ones([ len(res)  ,1])*Polarisation
    # res=map(Numeric.dot, res, pol)
    res = PPMcore.PPM_MAPDOT_V( res, pol)
    res=TR(  sum(TR( res*Numeric.conjugate(res)        )))

    print " PASSO DI QUA "
    return res.real




  def getDescription(self,wavelenghts):

    """ The stack is described by 
        Thickness, Roughness, indexlist, SubIndex, SubRough
        ( Sub means substrate )
    """

    thicknesslist=[]
    roughnesslist=[]
    indexlist=[]
     

    print " ci sono " , len(self.stack), " layers in stack "
    depth_count = 0.0
    for layer in self.stack:
        Minimiser.DependentVariable.depth = depth_count
        layerdescr= layer.getThickRoughIndexes( wavelenghts)
        thicknesslist.append( layerdescr[0])
        roughnesslist.append( layerdescr[1])
        indexlist    .extend( layerdescr[2])
        depth_count = depth_count + 1.0
        
    Minimiser.DependentVariable.depth = 0

    indexlist.append( array([ [[1.,0,0],[0,1,0],[0,0,1]] ]  *  len(indexlist[0])  )*complex(1.,0) )
    roughnesslist.append(array([0.0]))
    thicknesslist.append(array([0.0]))
 
    Thickness= concatenate(thicknesslist)
    Roughness= concatenate(roughnesslist)
    
    SubRough=Roughness[0]*1.0    
    SubIndex=indexlist    [0]

    return Thickness, Roughness, indexlist, SubIndex, SubRough



  def calculateFields(self,wavelenghts, angles):

    """   returns APLUS, AMOIN, eigenstuff
          APLUS[i], AMOIS[i] with i starting at 0 ( vacuum)
          up to NI ( number of interfaces ) standing for
          substrate, are the coefficients of the upgoing and downgoing waves.

          Electric and magnetic fields can be reconstracted from knowledge of
          Aplus and Amois in the whole stack, by multipling the E and B
          vectors contained in eigenstuff 
    """

    Thickness, Roughness,indexlist,SubIndex,SubRough = self.getDescription(wavelenghts)

    eigenstuff=PPMcore.PPM_calculatescanTensor(
				Thickness[1:],Roughness[1:],
				wavelenghts,
                                angles,indexlist[1:],
				SubIndex,SubRough)

   
    K0=2*pi/wavelenghts

    propagazioni = Propagations( eigenstuff, Thickness,K0)

    interfacce= Interfaces(eigenstuff,Roughness,K0 )
	
    NI=len(interfacce)
    NW=len(wavelenghts)

    ###############################################################
    # reflectivity on the substrate
    #

    inter = interfacce[0]
    # res = map( Numeric.dot, inter[:,0:2,2:4 ]  ,map(LinearAlgebra.inverse, inter[:,2:4,2:4 ])  )    
    res   = PPMcore.PPM_MAPDOT(  inter[:,0:2,2:4 ]  ,PPMcore.PPM_MAPINV(inter[:,2:4,2:4 ] )  )    

    ###################################################################
    # calculation of reflectivity at every interface
    #
    Ref={}
    Ref[0]=res
    for i in range(1,NI):
       Ref[i] = Reflectivity( Ref[i-1], propagazioni[i],   interfacce[i])

    ########################################################################
    # starting from the surface we go back
    #
    APLUS={}
    AMOIN={}

    propagazioni = Propagations( eigenstuff, Thickness,K0, reverse=1)
    interfacce= Interfaces(eigenstuff,Roughness,K0,  reverse=1 )

    APLUS[NI] = TR(ones(len(wavelenghts) )*TR(    array([[[1.,0],[0,1.]]],"D")    ))
    AMOIN[NI] =  Ref[NI-1]
    
    b= zeros([ NW ,4,2]           ,"D")
    b[:,0:2,:] = APLUS
    b[:,2:4,:] = AMOIN

    for i in range(NI-2,-1,-1):
       b = propaga(b, propagazioni[i], interfacce[i+1])
       APLUS[i+1] = b[:,0:2,:]
       AMOIN[i+1] = b[:,2:4,:]
      
    return APLUS, AMOIN, eigenstuff







def Interfaces(eigenstuff,Roughness,K0 ):
  """ this routine is in charge of creating 
      an array of matrices representing
      the interfaces.
  """
  conditions=[]
  conditions_inv=[]

  print " LA lunghezz di eigenstuff est ", len(eigenstuff)

  for i in range(len(eigenstuff)):
    E= eigenstuff[i][2]
    B= eigenstuff[i][3]
    tok=zeros( E.shape[:2] + (4,) , 'D'  )
    tok[:,:, :2] = E[:,:,:2] #  parallel components
    tok[:,:,2:4] = B[:,:,:2] #  parallel components
    tok=Numeric.swapaxes(tok,1,2)
    conditions.append(tok)
    #conditions_inv.append(  map( LinearAlgebra.inverse,tok )          ) 
    conditions_inv.append(   PPMcore.PPM_MAPINV4(tok)  ) 
  # print  "conditions      " , conditions
  # print  "conditions_inv  " , conditions_inv
  res=[]
  print " PRIMA PARTE OK"
  for i in range(len(eigenstuff)-1  ):
     # TODO : scriver map numerico 4X4
     # toapp =   array(map(Numeric.dot, conditions_inv[i+1], conditions[i]) )
     toapp =   PPMcore.PPM_MAPDOT4( conditions_inv[i+1], conditions[i]) 


    #  ROUGHNESS

     for k1 in range(4):
       for k2 in range(4):
         toapp_inset = toapp[:,k1,k2]
	 dw = K0*Roughness[i]*( eigenstuff[i][1][:,k2] - eigenstuff[i+1][1][:,k1] );
         Numeric.multiply(toapp_inset, Numeric.exp(  -dw*dw/2.), toapp_inset)
     res.append( toapp )
  print " RETURN "
  return res


def Propagations( eigen, Thickness,K0):

  """ the propagations matrices are diagonal matrices obtained from
      knowledge of Kz eigenvalues.
      The return value is a list of NW X 4 arrays.
      Where is dimension having lenght=4 
      correspond just the diagonal.
      (NW X4 X 4 would be useless ) 
  """

  res=[]
  Thickness=Thickness*complex(0,1)
  for i in range(len(Thickness)):
    kz=TR( TR(eigen[i][1])*K0)
 
    t=Thickness[i]
    espos = kz*t
    res.append(exp(espos))
    
  return res

def __repr__(self):
    print " ************* PPM_Stack ******************"
    count=0
    for layer in self.stack:
       print " _____________ LAYER #%d __________________" % count
       print layer
       count=count+1
    return ''




     
class PPM_LayerBaseClasse:
  """ Base class for a layer.

      It as the __add__ function that overloads + operatore
      ( layer! + layer2  creates a stack )
      
      the function  getThickRoughIndexes is reimplemented by all the 
      other classes.

  """

  def __init__(self, 
              ):
    self.classe="PPM_LayerBaseClass"
 

  def __repr__(self):

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


  def __add__(self, other):
    newstack=[]
    if( 'classe' not in dir(other)):

       print "Error summing a wrong class "
       return 0

    if(other.classe  ==   'PPM_Stack'):

       connect(self,other.stack[0])
       connect(other.stack[-1], self)
       newstack=[self] +other.stack

    elif(other.classe in [  'PPM_LayerBaseClass'] ):

       connect(self ,other)
       connect(other,self )
       newstack=[self]+[other]

    else:

       print "Error summing a wrong class "
       return 0

    res=PPM_Stack()
    res.stack=newstack
    return res
    
  def getThickRoughIndexes( self, wavelenght):
      return ( array((1.0,)), array((0.0,)), array( (ones(shape(wavelenght), Complex),)) ) 


class PPM_ElementalLayer(PPM_LayerBaseClasse):

  """ This by now the most widely used classes
      The creation of a new object is done following 
      this paradigma :

         __init__(self, thickness=0, roughness=0, 
               DabaxList=[], DensityList=[], MassList=[])
        
      DabaxList is a list of Dabax scatterers and/or
      index objects and/or tensor objects

       The object in DabaxList must have the members :
         is_a_tensor = 1/0
         is_an_index = 1/0

       in order for the function to know how to deal with them.

       DensityList may contain number or Variables. 
       MassList just numbers or void ( None ) 
       if mass does not apply ( index files for example ) 
       
   """

  def __init__(self, thickness=0, roughness=0, 
               DabaxList=[], DensityList=[], MassList=[]):

      PPM_LayerBaseClasse.__init__(self)

      self.thickness=thickness
      self.roughness=roughness
   
      self.DabaxList  = DabaxList
      self.DensityList= DensityList
      self.MassList   = MassList



  def getThickRoughIndexes( self, wavelenght):

      """  the returned result containes
           an array of thickness ( just one in this case because the layer is simple : not sliced )
           an array of roughness ( just one for the above reason )
           and array of dielectric matrices OR scalars
 
           Either scalr or matrix depends on the nature od the object contained in DabaxList.
           if one or more objects are of the tensorial type, matrixes will be returned
      """

      # print " ####################################################################################### "
      #  print " calcolo uno strato "

      Fs_s = [None,]*len(self.DabaxList)

      is_tensor_present = 0



      
       
      for i in range(0,len(self.DabaxList)):
         scatterer=self.DabaxList[i]
         mass     =self.MassList[i]
         print self.DensityList[i]
         if(hasattr(scatterer,"is_a_tensor") and scatterer.is_a_tensor):
            # print " ho trovato un tensore "
            is_tensor_present = 1
            #  print " Ottengo questo tensore "
            Fs_s[i] = scatterer.Eps_Lambda(wavelenght)
	 elif(hasattr(scatterer,"is_an_index") and scatterer.is_an_index):
	    Fs_s[i] = (1.0-scatterer.index(wavelenght))
         else:
           if( scatterer.__class__ == dabax.return_value ):
             Fs_s[i] = 1- CalculateIndexesFromTable( scatterer,wavelenght )
           else:
             Fs_s[i] =TR(415.22*(  TR( scatterer.F_Lambda(wavelenght)) )*( 
                                  wavelenght*wavelenght/12398.52/12398.52/mass))



      if( is_tensor_present == 1):
        # print " un tensore e presente"
        dielectric=zeros(shape(wavelenght)+(3,3),Complex)
        ONES   =ones(shape(wavelenght)+( 1,1 ),Complex) * MLab.eye(3)
      else:
        dielectric=zeros(shape(wavelenght) ,Complex)
        ONES   =ones(shape(wavelenght),Complex)

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

         if(Fs_s[i].shape[1:]==(3,3)):
           added_shape=()
           shapematrix = 1
         else:
           if(is_tensor_present == 1):
             added_shape=(1,1)
             shapematrix = MLab.eye(3)
           else:
             added_shape=()
             shapematrix = 1

         # print " be allora?" 
         if((not hasattr(scatterer,"is_a_tensor")) or  scatterer.is_a_tensor==0):
            add = (ONES-Numeric.reshape( Fs_s[i] , Fs_s[i].shape + added_shape  ) *shapematrix     )
            add=add*add-ONES
            dielectric=dielectric+ add * density
         else:
          dielectric=dielectric+(  Numeric.reshape( Fs_s[i] , Fs_s[i].shape + added_shape  )*shapematrix -ONES )*density

      
      dielectric = ONES + dielectric

      return ( array(( par(self.thickness),)),
               array(( par(self.roughness),)),
               array( (dielectric,) ) 
             ) 

      
class PPM_SimpleLayer(PPM_ElementalLayer):

  """  This class is just an utility one.
       It is used to simplify the use of PPM_ElementalLayer.

       The annoying lists DabaxList=[], DensityList=[], MassList=[]
       are passed as member of the instance material

       An instance containing those lists can be created by
       ancillary routines like :
         ComposedIndex
         IndexFromTable
  """
  def __init__(self, thickness=0, roughness=0,  material=None):
    apply(PPM_ElementalLayer.__init__, (self, thickness, roughness, material.DabaxList  
  , material.DensityList   , material.MassList  ))



class PPM_SlicedLayer(PPM_LayerBaseClasse):

  """  This class is not yet fully implemented
       but the idea is to provide a layer
       with vertical gradients, represented as a sequence of
       sublayers
  """

  def __init__(self, thickness=0, Nlayers=0,IndexFunction=None ):

      PPM_LayerBaseClasse.__init__(self)

      self.thickness=thickness
      self.Nlayers=Nlayers
      self.IndexFuntion=IndexFunction

      
  def getThickRoughIndexes( self, wavelenght):
      indexesarraylist=[]

      for i in range(0, self.Nlayers):
         indexesarraylist.append(self.IndexFunction(i, wavelenght  ))

      return ( ones(self.Nlayers)*par(self.thickness)/self.Nlayers,
               zeros(self.Nlayers, Double),
               array(indexesarraylist) 
             )




def writeFit2File(fit,filename):
  """ calculate the error functions
      and then dumps calculated and experimental data on a file
  """
  if( hasattr(fit,"classe")==0 or  fit.classe!="PPM_ComparisonTheoryExperiment"):
     raise " the object that you passed to the writeFit2File function is not of the good class"
  fit.error(nopartial=1)
  fit.write2File(filename)


class PPM_ComparisonTheoryExperiment:

  """ 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="PPM_ComparisonTheoryExperiment"
  def __init__(self,stacks,scan2stack, scanlist, weightlist, normlist=None, noise=None, width=0, printpartial=0,
               meritfunction="abslogdiff", shiftlist=None,CutOffRatio=None):
       
      self.CutOffRatio = CutOffRatio
      self.meritfunction=meritfunction
      self.width=width
      
      self.stacks=stacks
      self.scan2stack=scan2stack
      
      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
      self.noise=noise
      self.shiftlist=shiftlist

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

  def getitercounter(self):
      return self.itercounter

  def write2File(self, filename, tipo="ASCII"):
       if tipo=="ASCII":
            self.write2FileASCII( filename)
       else:
            self.write2FileNewPlot( filename)
            
  def write2FileASCII(self, filename):
    for count in range(len(self.calculatedscan ) ) :
       calculated=self.calculatedscan[count]
       scan = self.scanlist[count]
       f=open("%s%d" %(filename, count+1),"w")
       if(len(scan) in [3,4,5] ):

         for k in range(len(calculated)):

          if( scan[2] is not None):
               f.write("%e %e %e %e\n"%(calculated[k], scan[0][k], scan[1][k],scan[2][k]  ))
          else:
               f.write("%e %e %e\n"%(calculated[k], scan[0][k], scan[1][k] ))
       elif(len(scan)==2):
         for k in range(len(calculated)):
           f.write("%e %e %e\n"%(calculated[k], scan[0][k],scan[1][k] ))
       elif(len(scan)==1):
         for k in range(len(calculated)):
           f.write("%e %e\n"%(calculated[k], scan[0][k] ))
       f.close()
 
  def write2FileNewPlot(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")

       f.write("\n")
       nscanspec=1
       
       if(len(scan) in [3,4,5] ):
         expe=1
         scanstring="#N 4\n#L Yc  Wa  An  Ye  Bo\n"
         try:
              tmp=scan[2][0]
         except:
              scanstring="#N 3\n#L Yc  Wa  An\n"
              expe=0


         f.write("#S%d Scan\n%s"%(nscanspec,scanstring) )
         if(  scan[0][0]!=  scan[0][1]  and scan[1][0]!=  scan[1][1]):
              leading=-1
         elif(  scan[0][0]==  scan[0][1]  and scan[1][0]!=  scan[1][1]):
              leading=0
         else:
              leading=1

         print " leading est " , leading
         print scan[0][0],  scan[0][1],   scan[1][0],   scan[1][1]
         

         for k in range(len(calculated)):
           if(leading>=0 and k):
                if( scan[leading][k] != scan[leading][k-1]):
                     nscanspec+=1
                     f.write("#S%d\n%s"%(nscanspec,scanstring) )
           if(expe):
              f.write("%e %e %e %e\n"%(calculated[k], scan[0][k], scan[1][k],scan[2][k] ))
           else:
             f.write("%e %e %e\n"%(calculated[k], scan[0][k], scan[1][k]  ))
             
       f.write("#S%d Dummy\n#N 2\n#L X  Y\n"%(nscanspec+1))
       f.close()
 	 
  def GetResForGeFit(self):
 	        res=[]
 	        for scan in self.calculatedscan:
 	             res.append(scan)
 	        res=Numeric.concatenate(res)
 	        print res.shape
 	 
 	        # return Numeric.log(res)
		return res
 	   
  def GetDataForGeFit(self):
 	        res=[]
 	        for scan in self.scanlist:
 	             res.append(scan[2])
 	        res=Numeric.concatenate(res)
 	        # return Numeric.log(res)
 	        return res
 	 
  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.

     """

     if( ISMASTER and type(self.printpartial)==type([]) ):
         vs=self.printpartial[0]
         ns=self.printpartial[1]
         f=open("PARTIAL/variables%d" %  ( self.itercounter % 30)  , "w")
	 for k in range(len(vs) ):
           f.write("%s = %e \n" %(ns[k], vs[k].value )   )
           print vs[k].value ,
         print "\n"
         f.write("#"*80+"\n")
	 for k in range(len(vs) ):
           f.write("_%s_ = %e \n" %(ns[k], (vs[k].max-vs[k].min)*(vs[k].max-vs[k].min)/(vs[k].value-vs[k].min)/(vs[k].max-vs[k].value)   ) )
         f.close()



         try:
           import os
           cmd="ln -fs variables%d  PARTIAL/lastvariables " % ( ( self.itercounter % 30 )  )  
           os.system(cmd )
         except:
           pass
     
     error=0
     count=-1
      
     csc=0
     numtot=0

     # '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
     #   stores  the calculated values for eventual retrieval in writefit2file
     #   
     # '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
     if( hasattr(self, "calculatedscan")==0): self.calculatedscan=[]
     newscans=[]
     if ISPARALLEL:
       lenscans= len(self.scanlist)
       if( ( mpi.size  % lenscans )!=0 ):
         print " lenscans = ", lenscans, "  mpi.size = ", mpi.size
         raise " nproc not multiple of len(scanlist)"
     rs=Numeric.zeros(2,"d")
     
     rsstore=[]
     scancount=0
     for scan in self.scanlist:
       count=count+1
       if  (ISPARALLEL==0 or (count) == ( mpi.rank %  lenscans)):
        wavelenghts= scan[0]

        if(self.shiftlist is not None):
           shift=par(self.shiftlist[count])
           wavelenghts = 12398.52/wavelenghts + shift
           wavelenghts = 12398.52/wavelenghts 

        angles=scan[1]
        data=scan[2]
        polar=scan[3]
        if(len(scan)==5):
           errorweights=scan[4]
        else:
           errorweights=ones(len(data ))

        rs=zeros(len(wavelenghts),"d")
        for tok in polar:
          rsadd=0
          if( ISPARALLEL ):
            if(tok[1]==0):
              rsadd= self.stacks[self.scan2stack[count+1]].calculatescanTensor_Parallel(  wavelenghts ,  angles    ,tok[2] , lenscans)
            if(tok[1]==1):
              rsadd= self.stacks[self.scan2stack[count+1]].calculatescanTensor_Parallel(  wavelenghts ,  pi-angles    ,tok[2], lenscans )
          else:
            if(tok[1]==0):
              rsadd= self.stacks[self.scan2stack[count+1]].calculatescanTensor_(  wavelenghts ,  angles    ,tok[2] )
            if(tok[1]==1):
              rsadd= self.stacks[self.scan2stack[count+1]].calculatescanTensor_(  wavelenghts ,  pi-angles    ,tok[2] )
              
          if(self.noise!=None):
            rsadd=rsadd+par(self.noise[count])
        

          rs+=tok[0]*rsadd


        if(self.width!=0):
            print " CONVOLUTION"
            convoluted=ones(len(rs))*0.0
            sumcon    =ones(len(rs))*0.0
            ones_     =ones(len(rs))
 
            nmax= int(  4*self.width     )+1
            for shift in range(-nmax,nmax+1):
               x=shift
               x=x*x/2.0/self.width/self.width
               fact=exp(-x)

               shift_add(convoluted,rs,shift,fact)
               shift_add(sumcon,ones_,shift,fact)
            rs=convoluted/sumcon
        

        if(self.CutOffRatio is not None):
          if(type(self.CutOffRatio)==type([]) ):
               cor = self.CutOffRatio[count]
          else:
               cor = self.CutOffRatio
          CutOffRescale = Numeric.clip(angles/par( cor ), 0.0,1.0)
          
          rs=rs*CutOffRescale
        rsstore.append(rs)


     count=-1
     if( ISPARALLEL ):
          own_rs=rs.tostring()
     for scan in self.scanlist:
          count=count+1
          data=scan[2]
          norm=par(self.normlist[count])
          if( ISPARALLEL ):
               rs = mpi.bcastString(own_rs, count)
               rs=Numeric.fromstring(rs,"d")
          else:
               rs=rsstore[count]
          if(data is None):
               data=norm*rs
          if(len(scan)==4 and  scan[3] is not None):
               errorweights=scan[3]
          else:
               errorweights=ones(len(data ),'d')
              
          newscans.append(rs*norm)
          if(self.meritfunction=="abslogdiff"):
               difference=abs(log(norm*rs )-log(data))
               difference=difference*errorweights

          elif( self.meritfunction=="diffroot" ): 
               difference=  norm*rs-data
               difference=difference*difference*errorweights

          elif(self.meritfunction=="sin4" ): 
               difference=  norm*rs-data
               sangles=sin(angles)
               sangles=sangles/sangles[0]
               difference=abs(difference*sangles*sangles*sangles*sangles)
         
          error=error+sum(difference)*self.weightlist[count]
          numtot=numtot+len(difference)




     self.calculatedscan=newscans
     if(self.printpartial and nopartial==0):
          self.write2File("partialresults")
     if( ISMASTER):
          if(self.printpartial and nopartial==0):
               self.write2File("PARTIAL/partialresults%d_" % ( self.itercounter % 30 )   )
               try:
                    import os
                    for i in range(len(newscans)):
                         cmd="ln -fs partialresults%d_%d  PARTIAL/lastresults_%d " % ( ( self.itercounter % 30 )  , i+1, i+1   )  
                         os.system(cmd )
               except:
                    pass
         
     self.itercounter=self.itercounter+1
     return error/numtot


class MagScatterer:
    """ Reads index data from two files. One for polarization +.
        The other for polarization -. It has a method F_Lambda
        that returns an array of tensors ( matrices)
    """  

    default_F=complex(4.0,0.0)

    is_a_tensor = 1
    is_an_index = 0

    def __init__(self, NameFilePlus=None, NameFileMinus=None, versor = Numeric.array( [0,1.0,0.0] ),
	 RelativeDensity=1.0, Saturation = 1 , ismagnetico=1   ):

        """
          Gli argomenti sono: 
              NameFilePlus/Minus: nomi dei file per gli indici relativi alle due polarizzazioni. (Energy, Real(n), Imm(n))
              saturation: fattore di scala per il dicroismo
              versor: direzione della magnetizzazione
        """
        self.RelativeDensity=RelativeDensity
        if( type(NameFilePlus) == type("s") ):
          self.fromfile  = 1
  	  self.dataPlus  = IndexReader( NameFilePlus   )
	  self.dataMinus = IndexReader( NameFileMinus  )
        else:
          self.fromfile  =  0
          self.dataPlus  =  NameFilePlus  
	  self.dataMinus =  NameFileMinus  


   	self.ismagnetico=ismagnetico

        if(self.ismagnetico):
          antisymm = Numeric.zeros([3,3,3])
          antisymm[0,1,2]=antisymm[1,2,0]=antisymm[2,0,1]=1.0
          antisymm[0,2,1]=antisymm[1,0,2]=antisymm[2,1,0]=-1.0
 
       
          self.rotazione = Numeric.dot(antisymm,versor)
        else:
          self.rotazione = Numeric.array(versor)[:,None]*Numeric.array(versor)


        self.identity  = MLab.eye(3)  

        self.recorder={}

        self.saturation=Saturation

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



    def F_Lambda(self, Lambda, theta=0):
         raise " It is a tensor !!!!!"
    def F_Energy(self, Energy, theta=0):
         raise " It is a tensor !!!!!"


        
    def Eps_Lambda(self, Lambda, theta=0):

      """
           descrizionde di F_Lambda
      """
      Lambda=Numeric.array(Lambda)
      return self.Eps_Energy( 12398.52/Lambda, theta  )

    def Eps_Energy(self, Energy, theta=0):

      """
         Interpola la parte reale ed immaginaria degli indici  sui valori di energia passati come argomento.
      """
      if( isinstance(Energy,Numeric.ArrayType) or  isinstance(Energy,type([]) ) ):
           interpolateAt=Numeric.array(Energy)
      else:
           interpolateAt=Numeric.array([Energy])
     

      key = interpolateAt.tostring()


      if(self.fromfile):

        if(key not in self.recorder.keys() ):
          nrP = arrayfns.interp(self.dataPlus[:,1], self.dataPlus[:,0], interpolateAt  )
          niP = arrayfns.interp(self.dataPlus[:,2], self.dataPlus[:,0], interpolateAt  )
          nrM = arrayfns.interp(self.dataMinus[:,1], self.dataMinus[:,0], interpolateAt  ) 
          niM = arrayfns.interp(self.dataMinus[:,2], self.dataMinus[:,0], interpolateAt  ) 
          indexPlus  = nrP+complex(0,1)*niP
          indexMinus = nrM+complex(0,1)*niM

          epsPlus  = indexPlus  * indexPlus        
          epsMinus = indexMinus * indexMinus        

          semisomma  = 0.5*(epsPlus+epsMinus )
          differenza =      epsMinus - epsPlus

          tensorDiag= Numeric.reshape( semisomma, semisomma.shape+(1,1) ) * self.identity   

          if(self.ismagnetico) :
	     	tensorRot = Numeric.reshape(differenza , differenza.shape+(1,1) )* self.rotazione*complex(0.0,1.0)
	  else:
		tensorRot = Numeric.reshape(differenza , differenza.shape+(1,1) )* self.rotazione

          self.recorder[key]=( tensorDiag ,tensorRot )
        else:
          ( tensorDiag ,tensorRot )= self.recorder[key]

      else:

        
##            indexPlus = self.dataPlus.index(12398.52/interpolateAt)
##            indexMinus = self.dataMinus.index(12398.52/interpolateAt)
           indexPlus = Mag_Aux_Index(  self.dataPlus , 12398.52/interpolateAt)
           indexMinus = Mag_Aux_Index( self.dataMinus, 12398.52/interpolateAt)

           epsPlus  = indexPlus  * indexPlus        
           epsMinus = indexMinus * indexMinus        
           semisomma  = 0.5*(epsPlus+epsMinus )
           differenza =      epsMinus - epsPlus
           tensorDiag= Numeric.reshape( semisomma, semisomma.shape+(1,1) ) * self.identity   
           if(self.ismagnetico) :
              tensorRot = Numeric.reshape(differenza , differenza.shape+(1,1) )* self.rotazione*complex(0.0,1.0)
	   else:
	      tensorRot = Numeric.reshape(differenza , differenza.shape+(1,1) )* self.rotazione

      tensor = tensorDiag + Minimiser.par(  self.saturation )*tensorRot*0.5

      # '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
      # we need dielectric constant, here why we multiply by 2
      #
      # print " tensor  ", tensor 

      return (tensor)


def Mag_Aux_Index(material, wavelenght):

      Fs_s = [None,]*len(material.DabaxList)

      for i in range(0,len(material.DabaxList)):
         scatterer=material.DabaxList[i]
         mass     =material.MassList[i]
         if(scatterer.is_an_index):
	    Fs_s[i] = (1.0-scatterer.index(wavelenght))
         else:
           if( scatterer.__class__ == dabax.return_value ):
             Fs_s[i] = 1- CalculateIndexesFromTable( scatterer,wavelenght )
           else:
             Fs_s[i] =TR(415.22*(  TR( scatterer.F_Lambda(wavelenght)) )*( 
                                  wavelenght*wavelenght/12398.52/12398.52/mass))


      index=ones(shape(wavelenght) ,Complex)
      ONES   =ones(shape(wavelenght),Complex)

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



         index = index -  Fs_s[i]  * density


      return index




 

def 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



















class return_value:
  pass


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

def ScanReader(filename=None, Np="automatic", wavelenghts_col=1,angles_col=2,Polarisation=None,  refle_col=None, weight_col=None,
               angle_factor=1.0):
   """ 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 ti 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.
   """  


 

   if(filename is not None):
     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=array(newdata)

   else:
      newdata=[]
      columns=[]

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

      if( weight_col is not None):
              raise " weight_col should be None"
      if(refle_col is not None):
              raise " refle_col should be None"
      data=array(newdata)
      data=transpose(data)


   res=return_value()

   for (colonna_n, colonna_dati) in [ (wavelenghts_col, "wavelenghts" ),(angles_col, "angles" ),
                                     (refle_col, "refle" ),(weight_col, "weight" ) ]:
    if( isinstance(colonna_n , type(1))):
      if( colonna_n> 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)
    else:
      message= " PROBLEM with %s_col type\n" % colonna_dati
      raise " PROBLEM with wavelenghts_col type\n"
    
   res.angles*=angle_factor
   return [res.wavelenghts, res.angles, res.refle, Polarisation, res.weight]












class Material:
    def __init__(self):
       self.DabaxList=[]
       self.DensityList=[]
       self.MassList=[]

    def append(self, DabaxList, DensityList, MassList):
       self.DabaxList+=DabaxList
       self.DensityList+=DensityList
       self.MassList+=MassList
          

def ComposedIndex(*args):
   res=Material()
   for tok in args:
     res.append( tok.DabaxList, tok.DensityList,  tok.MassList    )
   return res
 
def AppendIndexFromTable(res, Tablef0, Tablef12, *args):
   nel=len(args)
   if( nel%2 !=0):
      raise " Wrong Number of arguments in IndexFromTable"
   nel=nel/2
   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,)       )
     res.append([scatterer],[density], [mass])


def IndexFromTable(Tablef0, Tablef12, *args):
   res=Material()
   apply( AppendIndexFromTable, (res, Tablef0, Tablef12) + args )
   return res

#   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,)       )
#
#     res.append(DabaxList, DensityList,  MassList  )

def Minimise(fit = None, list_of_variables = None, tol=0.001, temperature=".05*exp(-0.2*x)",
             max_refusedcount=100  , max_isthesame=10,centershiftFact=0 ):
   minimizzatore=Minimiser.minimiser(fit, list_of_variables)
   ss="foo=lambda x: "+temperature
   exec(ss)
   minimizzatore.amoebaAnnealed( tol, foo , max_refusedcount,
                                 max_isthesame=max_isthesame,centershiftFact=centershiftFact )



def MinimiseGeFit(fit = None, list_of_variables = None ):

   wpf = Minimiser.wrapperForFit(fit,list_of_variables )
   xs  = Numeric.array(wpf.get_angles())
   data= Numeric.array(wpf.getData())
   data=Numeric.array([data,data])
   data=Numeric.transpose(data)

   print xs.shape
   print data.shape
   
   fittedpar, chisq, sigmapar = Minimiser.Gefit.LeastSquaresFit(wpf.theory,xs,data)
   wpf.apply_angles(fittedpar)
   fit.error()








