from sys import *
from copy import *
# import NumericA
from Numeric import *
from  arrayfns  import *            
# from pgplot import *
from umathextra import *
from Minimiser import *
from Dabax import *


TR=transpose

import profile

def connect(a,b):
   a.next=b
   b.before=a


def   shift_add(convoluted,rs,shift,fact):
  if(abs(shift)>len(convoluted)):
   return
  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 LayerBefore(self):
        if(self.before.classe=="Interface"):
           adjacent=self.before.before
        else:
           adjacent=self.before
        return adjacent

def LayerNext(self):
        if(self.next.classe=="Interface"):
           adjacent=self.next.next
        else:
           adjacent=self.next
        return adjacent




#######################################################################3
class Free_Strain_Diffusion:
    def __init__(self, TopSpacingsD=None,BottomSpacingsD=None,TopInterdiffusion=None,BottomInterdiffusion=None  ):
       self.TopSpacingsD         = TopSpacingsD
       self. BottomSpacingsD     = BottomSpacingsD
       self. TopInterdiffusion   = TopInterdiffusion
       self.BottomInterdiffusion = BottomInterdiffusion

class Fuller_Interface:
  def __init__(self, 
               Separation=None,
               SigmaSeparation=0,
               Interdiffusion=0,
               InterdiffusionExponent=0,
               StearnsWidth=None,
               REdiff0=0,
               REdiff1=0,
               I_Planes=None,
               Perov=None
              ):

    self.Separation=Separation
    self.Interdiffusion=Interdiffusion
    self.InterdiffusionExponent=InterdiffusionExponent
    self.classe="Interface"
    self.SigmaSeparation=SigmaSeparation
    self.StearnsWidth=StearnsWidth
    self.REdiff0=REdiff0
    self.REdiff1=REdiff1
    self.I_Planes=I_Planes

    if(Perov is  None):
      if(self.I_Planes is not None):
         self.Perov_Planes=Perov_Planes(I_Planes)
      else:
         self.Perov_Planes=None
    else:
      self.Perov_Planes=Perov 

  def __repr__(self):
    print "/////////////////////// ____ Interface  _____"
    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  ==   'PeriodM_I'):
       connect(self,other.stack[0])
       connect(other.stack[-1],self)
       newstack=[self] +other.stack
    elif(other.classe in [  'LayerM_I',  'LayerM_IV' ,  'LayerM_XXI' ] ):
       connect(self,other)
       connect(other,self)
       newstack=[self]+[other]
    else:
       print "Error summing a wrong class "
       return 0
    res=Fuller_PeriodM_I()
    res.stack=newstack
    return res
    
###############################################################################3
# a class to specify arbitrary deformations at the surface
#
class InterfacePlanes:
  NpL=0
  NpH=0
  ShiftsL=[]
  ShiftsH=[]
  MixingL=[]
  MixingH=[]

class Perov_Planes:
#            for  (diff, shift)  in self.next.Perov_Planes.L  :
   def __init__(self, I_Planes):
	self.L=[]
  	self.H=[]

        bs=I_Planes.BottomSpacingsD 
        if(bs is None):
            bs=[]
        bi=I_Planes.BottomInterdiffusion
        if(bi is None):
            bi=[]
        bs=bs+[0.0]*(  max(len(bi), len(bs)) -len(bs)     )
        bi=bi+[0.0]*(  max(len(bi), len(bs)) -len(bi)     )
        for i in range (  len (bs)    ):
           self.L.append(  (  bi[i] ,bs[i]  ) ) 
         


        bs=I_Planes.TopSpacingsD 
        if(bs is None):
            bs=[]
        bi=I_Planes.TopInterdiffusion
        if(bi is None):
            bi=[]
        bs=bs+[0.0]*(  max(len(bi), len(bs)) -len(bs)     )
        bi=bi+[0.0]*(  max(len(bi), len(bs)) -len(bi)     )
        for i in range (  len (bs)    ):
           self.H.append(  (  bi[i] ,bs[i]  ) ) 
         


#######################################################################3
 
class Fuller_ElementalPlane:
  def __init__(self, Element=None, InPlaneDensity=1):
    self.classe="ElementalPlane"
    self.Element=Element
    self.InPlaneDensity=InPlaneDensity

  def ScatteringPower(self, Lambda, theta):
    return  conjugate(self.Element.F_Lambda( Lambda, theta))*par(self.InPlaneDensity)

  def VolumeScatteringPower(self, Lambda, dens):
    # dens is the volumic density (Natomes  A**-3)
    return self.ScatteringPower(Lambda, 0.0   )*dens


  def DeltaBeta(self, Lambda, dens):  
    # dens is the volumic density (Natoms  A**-3)
    E=LAMBDAEV/Lambda
    return self.ScatteringPower(Lambda,  Lambda*0  )*415.22*dens/E/E/AVOGADRO*1.0e+24


  def __repr__(self):
    print "    ElementalPlane made of ",  self.Element, "; InPlaneDensity=", \
          self.InPlaneDensity
    return ''
#######################################################################3
 


def Fuller_MixHTc( a,b, f,spacing=None):
   if(a.classe!=b.classe):
      return
   if(a.classe=="HTcPlaneC"):
      res= Fuller_HTcPlaneC()
      res.Oxygen=(1-f)*a.Oxygen+f*b.Oxygen
      res.Copper=(1-f)*a.Copper+f*b.Copper
      res.Barium=(1-f)*a.Barium+f*b.Barium
      res.RE=(1-f)*a.RE+f*b.RE

      res.RE_Cu2  =(1-f)*par(a.RE_Cu2 )+f*par(b.RE_Cu2 ) 
      res.Cu2_Ba  =(1-f)*par(a.Cu2_Ba )+f*par(b.Cu2_Ba  ) 
      res.Ba_Cu1  =(1-f)*par(a.Ba_Cu1 )+f*par(b.Ba_Cu1 ) 
      res.Ba_O4   =(1-f)*par(a.Ba_O4  )+f*par(b.Ba_O4   ) 
      res.Cu2_O2  =(1-f)*par(a.Cu2_O2 )+f*par(b.Cu2_O2  ) 

      return res

   if(a.classe=="PerovskitePlane"):

      res=Fuller_PerovskitePlane ()
      res.Oxygen=a.Oxygen*(1-f)+b.Oxygen*f
      res.Internal=(1-f)*a.Internal+f*b.Internal
      res.External=(1-f)*a.External+f*b.External

      if(spacing is not None):
         res.Spacing=a.Spacing+spacing
      else:
         res.Spacing=(1-f)*a.Spacing+f*b.Spacing

      return res

class Fuller_HTcPlaneC:
  def __init__(self,   Oxygen=None, Copper=None, Barium=None, RE=None,
               RE_Cu2=None, Cu2_Ba=None, Ba_Cu1=None,
               Ba_O4=None, Cu2_O2=None, InPlaneDensity=1):

      self.Oxygen=Oxygen
      self.Copper=Copper
      self.Barium=Barium
      self.RE=RE

      self.RE_Cu2=RE_Cu2
      self.Cu2_Ba=Cu2_Ba
      self.Ba_Cu1=Ba_Cu1
      self.Ba_O4=Ba_O4
      self.Cu2_O2=Cu2_O2
      self.classe="HTcPlaneC"
      self.InPlaneDensity=InPlaneDensity

      
  def setArrays(self):
    Ba_Cu1=par(self.Ba_Cu1)
    Ba_O4 =par(self.Ba_O4)
    Cu2_Ba=par(self.Cu2_Ba)
    Cu2_O2=par(self.Cu2_O2)
    RE_Cu2=par(self.RE_Cu2)

    self.Total= 2*Ba_Cu1 +2*Cu2_Ba+ 2*RE_Cu2
    self.upperThickness=Ba_Cu1
    self.lowerThickness=Ba_Cu1

    self.O_array=array( (0,
                    Ba_Cu1+Ba_O4, 
                    Ba_Cu1 +Cu2_Ba + Cu2_O2 , 
                    Ba_Cu1 +Cu2_Ba+ 2*RE_Cu2-Cu2_O2  ,
                    Ba_Cu1 +2*Cu2_Ba+ 2*RE_Cu2 -Ba_O4  
                   )
                 )

    self.Cu_array=array( (0,
                    Ba_Cu1 +Cu2_Ba , 
                    Ba_Cu1 +Cu2_Ba+ 2*RE_Cu2  ,
                   )
                 )

    self.Ba_array=array( (
                    Ba_Cu1, 
                    Ba_Cu1 +2*Cu2_Ba+ 2*RE_Cu2  
                   )
                 )

    self.RE_array=array( (
                    Ba_Cu1 +Cu2_Ba+ RE_Cu2  
                   )
                 )

    self.O_weights  =array( (1,1,2,2,1)  )
    self.Cu_weights =array( (1,1,1)  )
    self.Ba_weights =array( (1,1) )
    self.RE_weights =array( (1) )

  def totalAbsorption(self, Lambda, theta ):
   pesi=(sum(O_weights), sum(Cu_weights), sum(Ba_weights),sum(RE_weights))
   atoms=(self.Oxygen, self.Copper, self. Barium, self.RE)
   res=0.0
   for i in range(0,4):
     res=res+(pesi[i]*atoms[i].DeltaBeta( Lambda, par(self.InPlaneDensity))).imag
   return exp(-2*res*2*math.pi/Lambda/sin(theta))


  def ScatteringPower(self, Lambda, theta):
   iKperp=complex(0,1)*sin(theta)*4*math.pi/Lambda


   #########################################################
   # ready to sum up
   #########################################################
   F=0 
            
   F=F+    ( self.Oxygen.ScatteringPower( Lambda, theta) *  
                 sum( self.O_weights * exp(iKperp*self.O_array) ))
   F=F+    ( self.Copper.ScatteringPower( Lambda, theta) *  
                 sum( self.Cu_weights * exp(iKperp*self.Cu_array) ))
   F=F+    ( self.Barium.ScatteringPower( Lambda, theta) *  
                 sum( self.Ba_weights * exp(iKperp*self.Ba_array) ))
   F=F+    ( self.RE.ScatteringPower( Lambda, theta) *  
                 sum( self.RE_weights * exp(iKperp*self.RE_array) ))
   return F


  def TotalThickness(self):
    return self.Total 

  def getUpperThickness(self):
    return self.upperThickness

  def getLowerThickness(self):
    return self.lowerThickness

  def __repr__(self):
    print "    ElementalPlane made of ",  self.Element, "; InPlaneDensity=", \
          self.InPlaneDensity
    return ''
 
#######################################################################3
      
class Fuller_LayerM_I:
  def __init__(self, Nplanes=1, SigmaNplanes=0,
               Plane=None,
               PlaneD =None,
               PlaneDD2=0,
               PlaneDD1=0,
               PlaneDDAlpha=0,
               Amorphous=0,
               Thickness=0,
               SigmaThickness=0,
               Density=0
              ):
    self.Plane=Plane
    self.PlaneD=PlaneD
    self.PlaneDD2=PlaneDD2
    self.PlaneDD1=PlaneDD1
    self.PlaneDDAlpha=PlaneDDAlpha
    self.Nplanes =Nplanes
    self.SigmaNplanes=SigmaNplanes
    self.Amorphous=Amorphous
    self.Thickness=Thickness
    self.SigmaThickness=SigmaThickness
    self.Density=Density


    self.classe="LayerM_I"
    self.model="unknown"



  def __repr__(self):
    print "////////////////////// ____ LayerM_I_____"
    for key in dir(self):
      if(key not in ['self', 'before','next','HeightsArrays']): 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  ==   'PeriodM_I'):
       connect(self,other.stack[0])
       connect(other.stack[-1], self)
       newstack=[self] +other.stack
    elif(other.classe in [  'LayerM_I',  'LayerM_IV' ,  'LayerM_XXI' ,   'Interface'] ):
       connect(self ,other)
       connect(other,self )
       newstack=[self]+[other]
    else:
       print "Error summing a wrong class "
       return 0

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

  def getUpperThickness(self):
    return par(self.PlaneD)+par(self.PlaneDD2)


  def getLowerThickness(self):
    return par(self.PlaneD)+par(self.PlaneDD1)



  def SetupSharedArrays(self):


      Nplanes=par(self.Nplanes)
      PlaneD=par(self.PlaneD)
      PlaneDD1=par(self.PlaneDD1)
      PlaneDD2=par(self.PlaneDD2)
      SigmaNplanes=par(self.SigmaNplanes)
      PlaneDDAlpha=par(self.PlaneDDAlpha)


      ####################################################
      #  Calculation of the low interface thickness
      ##################################################
      if(self.before.classe=="Interface" and par(self.before.Separation)!=None   ):
        SeparationLower=par(self.before.Separation)
        SigmaSeparationLower=par( self.before.SigmaSeparation)
      else:
        adjacent=LayerBefore(self)

        SeparationLower=( PlaneD+PlaneDD1+
                             adjacent.getUpperThickness() )*0.5

        if(self.before.classe=="Interface"  and par(self.before.SigmaSeparation)!=None):
          SigmaSeparationLower= par(self.before.SigmaSeparation)

        else:
          SigmaSeparationLower=0






      if(self.Amorphous):
          self.SeparationLower=SeparationLower
          self.SigmaSeparationLower=SigmaSeparationLower
          return

      ##     setup the arrays that depend only on structure and not
      ##      on Wavelenght/angle
      Npmin =int( round ( Nplanes ) -(round(4.243* abs( SigmaNplanes)   )+1))
      Npmax =int( round ( Nplanes )+ (round(4.243* abs( SigmaNplanes )  )+1)+1)


      if (Npmin < 0): Npmin=0
      ##
      # Probabilities calculation
      # as Fullerton, take two gaussians centered on the two
      # integers adjacent to Nplanes
      ##
      NlayersArray=arrayrange(Npmin,Npmax+1,1)
      CenterLow=int(Nplanes)
      CenterHigh=int(Nplanes)+1


      LowArraya=  NlayersArray+0.5   - CenterLow 
      LowArrayb=  clip( (NlayersArray-0.5),0.0, 1.e10 )   - CenterLow
         

      HighArraya=  NlayersArray+0.5   -CenterHigh  
      HighArrayb=   clip( NlayersArray-0.5,0.0, 1.e10 )   -CenterHigh  

      deno=(SigmaNplanes*sqrt(2.0))

      LowArraya=erf(LowArraya/deno)
      LowArrayb=erf(LowArrayb/deno)
      HighArraya=erf(HighArraya/deno)
      HighArrayb=erf(HighArrayb/deno)

      LowArray=LowArraya-LowArrayb
      HighArray=HighArraya-HighArrayb
      LowArray=LowArray/sum(LowArray)
      HighArray=HighArray/sum(HighArray)
      ##
      # we do exactly the same average as Fullerton
      ##
      # LowArray[0]=0 # ??

      Probabilities=(HighArray*(Nplanes-int(Nplanes))+
                     LowArray*(int(Nplanes+1)-Nplanes))
      ###
      # Normalization
      ###
      Probabilities=Probabilities/sum(Probabilities  )

      
      ####################################################
      #  Calculation of Npmin,Npmax set of thickness
      ##################################################
      HeightsArrays={}

      for nl in range(Npmin,Npmax+1):
         nep=nl-1
         dum=ones(max(1,nep))
         dum=dum*PlaneD
         D1=PlaneDD1
         D2=PlaneDD2
         alpha=exp(-PlaneDDAlpha)
         for i in range (0,3):
            l=i
            u=nep-1-i
            if(l>u): break
            if(l!=u):
               dum[l]=dum[l]+D1
               dum[u]=dum[u]+D2
            else:
               dum[l]=dum[l]+0.5*(D1+D2)
            D1=D1*alpha
            D2=D2*alpha
         # print dum 
         
         ####################################################################
         #  EVERYTHING-FREE INTERDIFFUSION
         #
         I_P_L=None
         I_P_H=None
         if(self.before.classe=="Interface" and par(self.before.I_Planes)!=None ):
              I_P_L=self.before.I_Planes

         if(self.next.classe=="Interface" and par(self.next.I_Planes)!=None ):
              I_P_H=self.next.I_Planes

         ###########################################################################
         #  TRAPEZOIDAL INTERDIFFUSION MODEL
         #

         Swa=0
         Swb=0
         Sda=0
         Sdb=0
         if(self.before.classe=="Interface" and par(self.before.StearnsWidth)!=None ):
           Swa=par(self.before.StearnsWidth)
           Sda=par(LayerBefore(self).PlaneD)
         if(self.next.classe=="Interface" and par(self.next.StearnsWidth)!=None ):
           Swb=par(self.next.StearnsWidth)
           Sdb=par(LayerNext(self).PlaneD)

         ####################################################################
         # Check that only one model at a time be used

         if( (Swa>0 or Swb>0) and (I_P_H!=None or I_P_L!=None) ) :
             raise " Use only one interdiffusion model at a time  "

         #######################################################33
         # TRAPEZOIDAL

         if(Swa>0 or Swb>0):
           for i in range (0, nep):

              if(Swa>0): pA1= (0.5*Swa - (1.0 +i))/Swa
              else:       pA1=0
              if(Swb>0): pA2= (0.5*Swb - (1.0 +(nep-1)-i))/Swb
              else      : pA2=0

              pA1=max(pA1,0)
              pA2=max(pA2,0)
 	      pA = pA1+pA2
              if(pA>0.5):
                pA1=pA1*0.5/pA
                pA2=pA2*0.5/pA
              dum[i]=dum[i]+pA1*(Sda-par(self.PlaneD)) 
              dum[i]=dum[i]+pA2*(Sdb-par(self.PlaneD)) 
         #
         # END OF THE TRAPEZOIDAL MODEL
         #  
         #######################################################33
         # EVERYTHING-FREE  MODEL

         if(I_P_L!=None   ):
           for i in range (0, min(I_P_L.NpH, nep) ):
              shft = par (  I_P_L.ShiftsH[i]  )
              dum[i]=dum[i]+ shft

         if(I_P_H!=None   ):
           for i in range (0, min(I_P_H.NpL, nep) ):
              shft = par (  I_P_H.ShiftsL[i]  )
              dum[nep-1-i]=dum[nep-1-i]+ shft






         HeightsArrays[nl-Npmin]= concatenate(((0,),cumsum(dum)))[0:nl]

      Heights=arrayrange(Npmax+1-Npmin)*0.0
      # in case of zeo layers we have to compensate 
      # for the non existing interface
      Heights[0]=-SeparationLower
      for i in range (max(Npmin,1),Npmax+1):
        Heights[i-Npmin]=  (HeightsArrays[i-Npmin]) [-1]

      ###############################################################
      # Calculation of average thickness and variance
      ###############################################################

      means= Probabilities *  Heights
      AveragedThickness=sum(means)
      means2=means*Heights

      ###
      # Things to remember
      ##
      self.Probabilities=Probabilities
      self.NlayersArray =NlayersArray
      self.HeightsArrays=HeightsArrays
      self.Heights=Heights
      self.AveragedThickness=AveragedThickness
      self.SeparationLower=SeparationLower
      self.SigmaSeparationLower=SigmaSeparationLower

      # print " Averaged Thickness= ", AveragedThickness       


  def ScatteringPower(self, Lambda, theta):
   iKperp=complex(0,1)*sin(theta)*4*math.pi/Lambda
   if(self.Amorphous):
     #############################################
     #  *** Amorphous Material A ***
     vtha = par(self.SigmaThickness)
     tha =  par(self.Thickness)
     dens=  par(self.Density)
     faa= self.Plane.VolumeScatteringPower( Lambda, dens)  
     q=sin(theta)*4*math.pi/Lambda

     x1 = exp(-q*q*vtha*vtha/2.)
     AvPh=ada = x1*exp(iKperp*tha)
     AvF=afa = complex(0,1)*faa*( 1 - ada)/q

     AvPhF = adafa = complex(0,1)*conjugate(faa)*(1 - ada)/q
     x2 = 2.*real(faa*conjugate(faa) )/(q*q)
     x3 = x2*(1.-x1*cos(q*tha))
     AvFF=afa2 = x3 
     beta=real(self.Plane.DeltaBeta( Lambda, dens))    
     self.totalAbsorption=exp(-2*tha*beta*2*math.pi/Lambda/sin(theta))

   else:
      beta                =  (self.Plane.DeltaBeta( Lambda, par(self.Plane.InPlaneDensity)) ).imag    
      self.totalAbsorption=  exp(-2* par(self.Nplanes)    *beta*2*math.pi/Lambda/sin(theta)  )  

      #############################################
      # Scattering power of a plane
      #############################################
      f=self.Plane.ScatteringPower( Lambda, theta)  

      ##############################################################
      # Taking into account the interface diffusion if any 
      ##############################################################
      lowerdiffOK=0
      upperdiffOK=0

      if( self.before.classe=='Interface'):
         fl=self.before.before. Plane.ScatteringPower( Lambda, theta)
      else:
         fl=self.before.Plane.ScatteringPower( Lambda, theta)

      if( self.next.classe=='Interface'):
         fu=self.next.next. Plane.ScatteringPower( Lambda, theta)
      else:
         fu=self.next.Plane.ScatteringPower( Lambda, theta)
         
      
      if( self.before.classe=='Interface' and par(self.before.Interdiffusion)>0):
        lowerdiffOK=1
        lowerdiff =par(self.before.Interdiffusion)
        A=exp(-abs(before.InterdiffusionExponent))
        fld=(fl*lowerdiff+(1-lowerdiff)*f,
             fl*(lowerdiff*A)+(1-lowerdiff*A)*f ,
             fl*(lowerdiff*A*A)+(1-lowerdiff*A*A)*f
            )

      if( self.next.classe=='Interface' and par(self.next.Interdiffusion)>0):
        upperdiffOK=1
        upperdiff =par(self.next.Interdiffusion)
        A=exp(-par(self.next.InterdiffusionExponent))
        fud=(fu*upperdiff+(1-upperdiff)*f,
             fu*(upperdiff*A)+(1-upperdiff*A)*f ,
             fu*(upperdiff*A*A)+(1-upperdiff*A*A)*f
            )
      ####################################################################
      #  EVERYTHING-FREE INTERDIFFUSION
      #
      I_P_L=None
      I_P_H=None
      N_P_L=0
      N_P_H=0
      if(self.before.classe=="Interface" and par(self.before.I_Planes)!=None ):
          I_P_L=self.before.I_Planes
          fl=self.before.before.Plane.ScatteringPower( Lambda, theta)        
          N_P_L= self.before.I_Planes.NpL
      if(self.next.classe=="Interface" and par(self.next.I_Planes)!=None ):
          I_P_H=self.next.I_Planes
          fu=self.next.next.Plane.ScatteringPower( Lambda, theta)        
          N_P_H= self.next.I_Planes.NpH

      ###############################################################################3
      #  TRAPEZOIDAL INTERDIFFUSION MODEL
      #
      Swa=0
      Swb=0
      if(self.before.classe=="Interface" and par(self.before.StearnsWidth)!=None ):
        Swa=par(self.before.StearnsWidth)
        fl=self.before.before.Plane.ScatteringPower( Lambda, theta)        

      if(self.next.classe=="Interface" and par(self.next.StearnsWidth)!=None ):
        Swb=par(self.next.StearnsWidth)
        fu=self.next.next.Plane.ScatteringPower( Lambda, theta)        

      ####################################################################
      # Check that only one model at a time be used

      if( (Swa>0 or Swb>0) and (I_P_H!=None or I_P_L!=None) ) :
             raise " Use only one interdiffusion model at a time  "

      ############################################
      # for every thickness of the ensemble we
      # calculate the scattering powers of each layer 
      ############################################   
      scatterings={}
      npmin=self.NlayersArray[0]
      _ones=ones



      for nl in  self.NlayersArray.tolist():
         scatts=outerproduct(  _ones(nl, Complex), f)
         if(upperdiffOK and lowerdiffOK):
             for l in range (0,3):
                u=nl-l
                if (l>u): break
                if(l==u):
                  scatts[l]=0.5*(fld[l]+fud[l])
                else:
                  scatts[l]=fld[l]
                  scatts[u]=fud[l]
         elif( upperdiffOK):
             for l in range (0,3):
                u=nl-l
                if (u<0): break
                scatts[u]=fud[l]
         elif( lowerdiffOK):
             for l in range (0,3):
                if (l>nl): break
                scatts[l]=fld[l]
 
         ##############################################
         #  TRAPEZOIDAL INTERDIFFUSION MODEL
         #
         if(Swa>0 or Swb>0):



           nep=nl-1
           for i in range (0, nep+1):
              if(Swa>0): pA1= (0.5*Swa - (0.5+i))/Swa
              else:       pA1=0
              if(Swb>0): pA2= (0.5*Swb - (0.5+(nep)-i))/Swb
              else      : pA2=0
              pA1=max(pA1,0)
              pA2=max(pA2,0)
 	      pA = pA1+pA2
              if(pA>0.5):
                pA1=pA1*0.5/pA
                pA2=pA2*0.5/pA
              scatts[i]=scatts[i]+pA1*(fl-f) 
              scatts[i]=scatts[i]+pA2*(fu-f) 
         #
         # END OF THE TRAPEZOIDAL MODEL
         #  
         #######################################################
         # EVERYTHING-FREE  MODEL

         if(I_P_L!=None or  I_P_H!=None  ):
           nep=nl-1
           for i in range (0, nep+1) :
              if(N_P_L>0   and   N_P_L> i  ):
                pA1=par(I_P_L.MixingH[i])
              else:
                pA1=0

              if(N_P_H>0 and N_P_H>(nep-i)     ):
                pA2=par(I_P_H.MixingL[nep-i])
              else:
                pA2=0
              if(pA>0.5):
                pA1=pA1*0.5/pA
                pA2=pA2*0.5/pA
              scatts[i]=scatts[i]+pA1*(fl-f) 
              scatts[i]=scatts[i]+pA2*(fu-f) 
 
         #
         ########################################################


         scatterings[nl-npmin]=scatts
         self.scatterings=scatterings

      #####################################################
      # ready to sum up
      #########################################################
             
      AvF=0        # Averaged F
      AvFF=0       # Averaged conj(F)*F
      AvPh=0       # Averaged Phase through thickness (< exp(iKperp*Thickness)>)
      AvPhF=0      # Averaged product of conj(F) X exp(iKperp*Thickness)

      s=scatterings
      ik=iKperp
      h=self.HeightsArrays


      scats= array([sum(  ( s[i]*   exp(outerproduct(h[i], ik )))       )
                        for i in range(0,len(self.NlayersArray))])

      probs=self.Probabilities
      phases= exp( outerproduct(iKperp, self.Heights)     )
     
      scats=swapaxes(scats,0,1)

      AvPh= sum( swapaxes( probs*phases,0,1 )     )
      AvF  =sum( swapaxes( probs*  scats       ,0,1)    )
      AvFF =sum( swapaxes( probs* (scats*conjugate(scats)),0,1)  )
      AvPhF=sum( swapaxes( probs*conjugate(scats)*phases,0,1)   )

#      for i in  range(0,len(self.Probabilities)):
#         thickness=self.Heights[i]
#         phase=  exp( iKperp*thickness )
#         prob=self.Probabilities[i]
#         AvPh =AvPh  + prob*phase
##
#         nl=self.NlayersArray[i]
#         if(nl==0): continue
#
##         heights=self.HeightsArrays[i]
#
##         scat=     sum( scatterings[nl]* exp(iKperp*heights) )
#
#         AvF  =AvF   + prob*scat
##         AvFF =AvFF  + prob*scat*conjugate(scat)
# #        AvPhF=AvPhF + prob*conjugate(scat)*phase
#

# end of the if: else: for amorphous

   InterPh = exp( iKperp*self.SeparationLower +
                     iKperp*iKperp*self.SigmaSeparationLower*self.SigmaSeparationLower/2.0  )



   (self.AvF, self.AvFF,
       self.AvPh, self.AvPhF,
       self.InterPh )=( AvF,AvFF,AvPh, AvPhF, InterPh )


#######################################################################3
      
class Fuller_LayerM_IV:
  def __init__(self, Nplanes=1, SigmaNplanes=0,
               Plane=None,
               PlaneD =None,
               PlaneDD2=0,
               PlaneDD1=0,
               PlaneDDAlpha=0,
               ContinuosRoughness=0
              ):
    self.Plane=Plane
    self.PlaneD=PlaneD
    self.PlaneDD2=PlaneDD2
    self.PlaneDD1=PlaneDD1
    self.PlaneDDAlpha=PlaneDDAlpha
    self.Nplanes=Nplanes
    self.SigmaNplanes=SigmaNplanes
    self.ContinuosRoughness=ContinuosRoughness


    self.classe="LayerM_IV"
    self.model="unknown"


  def getUpperThickness(self):
    return par(self.PlaneD)+par(self.PlaneDD2)

  def getLowerThickness(self):
    return par(self.PlaneD)+par(self.PlaneDD1)

  def __repr__(self):
    print "////////////////////// ____ LayerM_IV_____"
    for key in dir(self):
      if(key not in ['self', 'before','next','HeightsArrays']): 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  ==   'PeriodM_I'):
       connect(self,other.stack[0])
       connect(other.stack[-1], self)
       newstack=[self] +other.stack
    elif(other.classe in [  'LayerM_I',  'LayerM_IV' ,  'LayerM_XXI' ,   'Interface']):
       connect(self ,other)
       connect(other,self )
       newstack=[self]+[other]
    else:
       print "Error summing a wrong class "
       return 0

    res=Fuller_PeriodM_I()
    res.stack=newstack
    return res
    
  def SetupSharedArrays(self):

      Nplanes=par(self.Nplanes)
      PlaneD=par(self.PlaneD)
      PlaneDD1=par(self.PlaneDD1)
      PlaneDD2=par(self.PlaneDD2)
      SigmaNplanes=par(self.SigmaNplanes)
      PlaneDDAlpha=par(self.PlaneDDAlpha)
#      delta=self.ContinuousRoughness

      ##     setup the arrays that depend only on structure and not
      ##      on Wavelenght/angle
      Npmin =int( round ( Nplanes ) -(round(4.243* abs( SigmaNplanes)   )+1))
      Npmax =int( round ( Nplanes )+ (round(4.243* abs( SigmaNplanes )  )+1)+1)


      if (Npmin < 0): Npmin=0
      ##
      # Probabilities calculation
      # as Fullerton, take two gaussians centered on the two
      # integers adjacent to Nplanes
      ##
      NlayersArray=arrayrange(Npmin,Npmax+1,1)
      CenterLow=int(Nplanes)
      CenterHigh=int(Nplanes)+1


      LowArraya=  NlayersArray+0.5   - CenterLow 
      LowArrayb=  clip( (NlayersArray-0.5),0.0, 1.e10 )   - CenterLow
         

      HighArraya=  NlayersArray+0.5   -CenterHigh  
      HighArrayb=   clip( NlayersArray-0.5,0.0, 1.e10 )   -CenterHigh  

      deno=(SigmaNplanes*sqrt(2.0))

      LowArraya=erf(LowArraya/deno)
      LowArrayb=erf(LowArrayb/deno)
      HighArraya=erf(HighArraya/deno)
      HighArrayb=erf(HighArrayb/deno)




      LowArray=LowArraya-LowArrayb
      HighArray=HighArraya-HighArrayb
      LowArray=LowArray/sum(LowArray)
      HighArray=HighArray/sum(HighArray)
      ##
      # we do exactly the same average as Fullerton
      ##
      # LowArray[0]=0 # ??

      Probabilities=(HighArray*(Nplanes-int(Nplanes))+
                     LowArray*(int(Nplanes+1)-Nplanes))
      ###
      # Normalization
      ###
      Probabilities=Probabilities/sum(Probabilities  )

      
      ####################################################
      #  Calculation of the low interface thickness
      ##################################################
      if(self.before.classe=="Interface" and par(self.before.Separation)!=None   ):
        SeparationLower=par(self.before.Separation)
        SigmaSeparationLower=par( self.before.SigmaSeparation)
      else:
        if(self.before.classe=="Interface"):
           adjacent=self.before.before
        else:
           adjacent=self.before
        SeparationLower=( PlaneD+PlaneDD1+
                             adjacent.getUpperThickness())*0.5
        SigmaSeparationLower=0


      ####################################################
      #  Calculation of Npmin,Npmax set of thickness
      ##################################################
      HeightsArrays={}
      for nl in range(Npmin,Npmax+1):
         nep=nl-1
         dum=ones(max(1,nep))
         dum=dum*PlaneD
         D1=PlaneDD1
         D2=PlaneDD2
         alpha=exp(-PlaneDDAlpha)
         for i in range (0,3):
            l=i
            u=nep-1-i
            if(l>u): break
            if(l!=u):
               dum[l]=dum[l]+D1
               dum[u]=dum[u]+D2
            else:
               dum[l]=dum[l]+0.5*(D1+D2)
            D1=D1*alpha
            D2=D2*alpha
         # print dum 
         HeightsArrays[nl-Npmin]=concatenate( (
                                        (0, ),
                                         cumsum(dum)
                                        ))[0:nl]

      Heights=arrayrange(Npmax+1-Npmin)*0.0
      # in case of zeo layers we have to compensate 
      # for the non existing interface
      Heights[0]=-SeparationLower
      for i in range (max(Npmin,1),Npmax+1):
        Heights[i-Npmin]=  (HeightsArrays[i-Npmin]) [i-1]
     
      ###############################################################
      # Calculation of average thickness and variance
      ###############################################################

      means= Probabilities *  Heights
      AveragedThickness=sum(means)
      means2=means*Heights

        
      ###
      # Things to remember
      ##
      self.Probabilities=Probabilities
      self.NlayersArray =NlayersArray
      self.HeightsArrays=HeightsArrays
      self.Heights=Heights
      self.AveragedThickness=AveragedThickness
      self.SeparationLower=SeparationLower
      self.SigmaSeparationLower=SigmaSeparationLower

      # print " Averaged Thickness= ", AveragedThickness       

  def ScatteringPower(self, Lambda, theta):

      beta                =  (self.Plane.DeltaBeta( Lambda, par(self.Plane.InPlaneDensity)) ).imag    
      self.totalAbsorption=  exp(-2* par(self.Nplanes)    *beta*2*math.pi/Lambda/sin(theta) )


      #############################################
      # Scattering power of a plane
      #############################################
      f=self.Plane.ScatteringPower( Lambda, theta)  


      ##############################################################
      # Taking into account the interface diffusion if any 
      ##############################################################
      lowerdiffOK=0
      upperdiffOK=0
      if( self.before.classe=='Interface' and par(self.before.Interdiffusion)>0):
        lowerdiffOK=1
        lowerdiff =par(self.before.Interdiffusion)
        A=exp(-abs(par(self.before.InterdiffusionExponent)))
        fl=self.before.before. ScatteringPower( Lambda, theta)        
        fld=(fl*lowerdiff+(1-lowerdiff)*f,
             fl*lowerdiff*A+(1-lowerdiff*A)*f ,
             fl*lowerdiff*A*A+(1-lowerdiff*A*A)*f
            )

      if( self.next.classe=='Interface' and par(self.next.Interdiffusion)>0):
        upperdiffOK=1
        upperdiff =par(self.next.Interdiffusion)
        A=exp(-par(self.next.InterdiffusionExponent))
        fu=next.next.ScatteringPower( Lambda, theta)        
        fud=(fu*upperdiff+(1-upperdiff)*f,
             fu*upperdiff*A+(1-upperdiff*A)*f ,
             fu*upperdiff*A*A+(1-upperdiff*A*A)*f
            )
 
      ############################################
      # for every thickness of the ensemble we
      # calculate the scattering powers of each layer 
      ############################################   
      scatterings={}
      for nl in  self.NlayersArray.tolist():
         scatterings[nl]=outerproduct(ones(nl, Complex),f)
         if(upperdiffOK and lowerdiffOK):
             for l in range (0,3):
                u=nl-l
                if (l>u): break
                if(l==u):
                  scatterings[nl][l]=0.5*(fld[l]+fud[l])
                else:
                  scatterings[nl][l]=fld[l]
                  scatterings[nl][u]=fud[l]
         elif( upperdiffOK):
             for l in range (0,3):
                u=nl-l
                if (u<0): break
                scatterings[nl][u]=fud[l]
         elif( lowerdiffOK):
             for l in range (0,3):
                if (l>nl): break
                scatterings[nl][l]=fld[l]
 
      #####################################################
      # ready to sum up
      #########################################################
      iKperp=complex(0,1)*sin(theta)*4*math.pi/Lambda
             
      AvF=0        # Averaged F
      AvFF=0       # Averaged conj(F)*F
      AvPh=0       # Averaged Phase through thickness (< exp(iKperp*Thickness)>)
      AvPhF=0      # Averaged product of conj(F) X exp(iKperp*Thickness)
      for i in  range(0,len(self.Probabilities)):
         thickness=self.Heights[i]
         prob=self.Probabilities[i]
         phase=  exp( iKperp*thickness   )
         nl=self.NlayersArray[i]
         sigma=par(self.ContinuosRoughness)
         # print "sigma est ", sigma
         dum=sigma*sigma*sin(theta)*4*math.pi/Lambda
         dum=-dum*dum/2.0
         AvPh =AvPh  + prob*phase *exp(dum*nl   )


         if(nl==0): continue
         heights=self.HeightsArrays[i]


         if(self.NlayersArray[i]<7):
           sigmaarray=outerproduct( arrayrange(0, nl), dum )
           scat=     sum( scatterings[nl]* exp(outerproduct(heights, iKperp )+sigmaarray) )
           AvF  =AvF   + prob*scat

           sigmaarray=outerproduct(((nl-1)-arrayrange(0, nl)),dum)
           scat=    sum( scatterings[nl]* exp(outerproduct(heights, iKperp )+sigmaarray) )  

           AvPhF=AvPhF + prob*conjugate(scat)*phase

           dumMat=zeros([nl,nl, len(iKperp) ]  , Complex)
           dumN  =zeros([nl,nl, len(iKperp) ]  , Complex)
           for k1 in range(0,nl):
             for k2 in range(0,nl):
                dumMat[k1][k2]= scatterings[nl][k1]*exp(iKperp*heights[k1])*conjugate(scatterings[nl][k2]*exp(iKperp*heights[k2]))
                dumN[k1][k2]  = exp(     dum*abs(k1-k2)   )
           
           dumMat=dumMat* dumN
           AvFF =AvFF  + prob*sum(sum( dumMat     ))

         else:
           scatts=scatterings[nl]
           d= heights[3]-heights[2]
           sumF=(  (1-exp( (nl-6)*(iKperp*d +dum  )  ))/
                  (1-exp( iKperp*d +dum    ))
                 )

           scat=0
           for k in range (0,3):
             scat=scat+( scatts[k]*
               exp(iKperp*heights[k]+k*dum))
             scat=scat+( scatts[nl-1-k]*
               exp(iKperp*heights[nl-1-k]+(nl-1-k)*dum))

           scat=scat+( scatts[3]*
               exp(iKperp*heights[3]+3*dum))*sumF
            
           AvF  =AvF   + prob*scat
          
           scat=0
           for k in range (0,3):
             scat=scat+( scatts[k]*
               exp(iKperp*heights[k]+(nl-1-k)*dum))
             scat=scat+( scatts[nl-1-k]*
               exp(iKperp*heights[nl-1-k]+k*dum))

           sumFbis=( exp( (nl-4)*dum )* (1-exp(+ (nl-6)*(iKperp*d -dum  )  ))/
                  (1-exp( +iKperp*d -dum    ))
                 )
           scat=scat+( scatts[3]*
               exp(iKperp*heights[3]))*sumFbis

           AvPhF  =AvPhF + prob*conjugate(scat)*phase


           for a in [0,1,2,nl-3,nl-2,nl-1]:
             for b in [0,1,2,nl-3,nl-2,nl-1]:
               if(b<a): continue
               add = scatts[b]*conjugate(scatts[a])*exp(  abs(b-a)*dum 
                             +iKperp*(heights[b] -heights[a])  )    
               if(b!=a):
                  add=2*add

               AvFF=AvFF+add.real*prob
           #
           # sum of the 4th termXothers               
           #
           for a in [0,1,2]:
              add=exp((3-a)*dum+iKperp*(heights[3]-heights[a]))*sumF*scatts[3]*conjugate(scatts[a])
              AvFF=AvFF + 2*add.real*prob

           for a in [nl-3,nl-2,nl-1]:
              add= exp(+iKperp*(heights[3]-heights[a]) + (a-(nl-1))*dum )* sumFbis*scatts[3]*conjugate(scatts[a])
              AvFF=AvFF + 2*add.real*prob
           #
           # sum of the 4th termX4               
           #

           # baraggio


           a=exp(iKperp*d+dum)
           a1=1-a

           en=exp((nl-7)*(iKperp*d+dum))
           add=((nl-6)+2*(a*(  (nl-7)/a1  -a*(1-en)/a1/a1
                        )).real)*scatts[3]*conjugate(scatts[3])

           AvFF=AvFF+add.real*prob

  
      InterPh = exp( iKperp*self.SeparationLower +
                     iKperp*iKperp*self.SigmaSeparationLower*self.SigmaSeparationLower/2.0  )

      InterPhi = exp( iKperp*  self.SeparationLower   )
      InterPhr = exp(  iKperp*iKperp*self.SigmaSeparationLower*self.SigmaSeparationLower/2.0  )


      (self.AvF, self.AvFF,
       self.AvPh, self.AvPhF,
       self.InterPh,
       self.InterPhi, self.InterPhr )=( AvF,AvFF,AvPh, AvPhF, InterPh,
                    InterPhi, InterPhr  )

      return  AvF,AvFF,AvPh, AvPhF, InterPh    


#######################################################################3
      
class Fuller_LayerM_XXI:
  def __init__(self, Nplanes=1, SigmaNplanes=0,Plane=None
             ):
    self.Nplanes=Nplanes
    self.SigmaNplanes=SigmaNplanes
    self.Plane=Plane
    self.classe="LayerM_XXI"

  def getUpperThickness(self):
    if(self.Plane.classe=="HTcPlaneC"):
       return self.Plane.getUpperThickness()
    else:
       return 0.0
       

  def getLowerThickness(self):
    if(self.Plane.classe=="HTcPlaneC"):
       return self.Plane.getLowerThickness()
    else:
       return 0.0
       
       

  def __repr__(self):
    print "////////////////////// ____ LayerM_XXI_____"
    for key in dir(self):
      if(key not in ['self', 'before','next','HeightsArrays']): 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  ==   'PeriodM_I'):
       connect(self,other.stack[0])
       connect(other.stack[-1], self)
       newstack=[self] +other.stack
    elif(other.classe in [  'LayerM_I',   'LayerM_IV' ,  'LayerM_XXI' ,   'Interface']):
       connect(self ,other)
       connect(other,self )
       newstack=[self]+[other]
    else:
       print "Error summing a wrong class "
       return 0

    res=Fuller_PeriodM_I()
    res.stack=newstack
    return res
    
  def SetupSharedArrays(self):

      ##     setup the arrays that depend only on structure and not
      ##      on Wavelenght/angle

      Nplanes=par(self.Nplanes)
      SigmaNplanes=par(self.SigmaNplanes)
      Npmin =int( round ( Nplanes ) -(round(4.243* abs( SigmaNplanes)   )+1))
      Npmax =int( round ( Nplanes )+ (round(4.243* abs( SigmaNplanes )  )+1)+1)


      if (Npmin < 0): Npmin=0
      ##
      # Probabilities calculation
      # as Fullerton, take two gaussians centered on the two
      # integers adjacent to Nplanes
      ##
      NlayersArray=arrayrange(Npmin,Npmax+1,1)
      CenterLow=int(Nplanes)
      CenterHigh=int(Nplanes)+1


      LowArraya=  NlayersArray+0.5   - CenterLow 
      LowArrayb=  clip( (NlayersArray-0.5),0.0, 1.e10 )   - CenterLow
         

      HighArraya=  NlayersArray+0.5   -CenterHigh  
      HighArrayb=   clip( NlayersArray-0.5,0.0, 1.e10 )   -CenterHigh  

      deno=(SigmaNplanes*sqrt(2.0))

      LowArraya=erf(LowArraya/deno)
      LowArrayb=erf(LowArrayb/deno)
      HighArraya=erf(HighArraya/deno)
      HighArrayb=erf(HighArrayb/deno)


      LowArray=LowArraya-LowArrayb
      HighArray=HighArraya-HighArrayb
      LowArray=LowArray/sum(LowArray)
      HighArray=HighArray/sum(HighArray)
      ##
      # we do exactly the same average as Fullerton
      ##
      # LowArray[0]=0 # ??

      Probabilities=(HighArray*(Nplanes-int(Nplanes))+
                     LowArray*(int(Nplanes+1)-Nplanes))
      ###
      # Normalization
      ###
      Probabilities=Probabilities/sum(Probabilities  )

      # print " considero  numeri di layers ", NlayersArray
      # print " con probabilit a" , Probabilities

      ####################################################
      #  Calculation of the low interface thickness
      ##################################################

      if(self.before.classe=="Interface" and par(self.before.Separation)!=None   ):
        SeparationLower=par(self.before.Separation)
        SigmaSeparationLower=par( self.before.SigmaSeparation)
      else:
        if(self.before.classe=="Interface"):
           adjacent=self.before.before
        else:
           adjacent=self.before
        SeparationLower=(  self.Plane.getLowerThickness() +  adjacent.getUpperThickness())*0.5
        SigmaSeparationLower=0


      # print  " SeparationLower est ", SeparationLower
      # print  " SigmaSeparationLower  est ", SigmaSeparationLower

      Planes_n=[]
      Planes_Nmn=[]
      self.interf_Planes_L=0
      self.interf_Planes_H=0

      if(self.Plane.classe=="HTcPlaneC"):
        if(self.before.classe=="Interface"    ):
            self.interf_Planes_L=2
            for REdiff in [par(self.before.REdiff0),par(self.before.REdiff1)]:
                Planes_n.append( Fuller_MixHTc( self.Plane, LayerBefore(self).Plane,  REdiff))
        if(self.next.classe=="Interface"    ):
            self.interf_Planes_H=2
            for REdiff in [par(self.next.REdiff0),par(self.next.REdiff1)]:
                Planes_Nmn.append( Fuller_MixHTc( self.Plane, LayerNext(self).Plane,  REdiff))

      self.Plane.setArrays()


      if(self.Plane.classe=="PerovskitePlane"):
        if(self.before.classe=="Interface"    ):
            self.interf_Planes_L=len(self.before.Perov_Planes.H)
            for  (diff, shift)  in self.before.Perov_Planes.H  :
                diff=par(diff)
                shift=par(shift)
                Planes_n.append( Fuller_MixHTc( self.Plane, LayerBefore(self).Plane,  diff, shift))


      if(self.Plane.classe=="WurtzitePlane"):
        if(self.before.classe=="Interface"    ):
            self.interf_Planes_L=len(self.before.Perov_Planes.H)
            for  (diff, shift)  in self.before.Perov_Planes.H  :
                diff=par(diff)
                shift=par(shift)
                Planes_n.append( Fuller_MixHTc( self.Plane, LayerBefore(self).Plane,  diff, shift))

        if(self.next.classe=="Interface"    ):
            self.interf_Planes_L=len(self.before.Perov_Planes.L)
            for  (diff, shift)  in self.next.Perov_Planes.L  :
                diff=par(diff)
                shift=par(shift)
                Planes_Nmn.append( Fuller_MixHTc( self.Plane, LayerNext(self).Plane, diff, shift ))

      self.Plane.setArrays()
      for tok in Planes_n+Planes_Nmn:
         tok.setArrays()

      ####################################################
      #  Calculation of Npmin,Npmax set of thickness
      ##################################################
      HeightsArrays={}



      Heights=arrayrange(Npmax+1-Npmin)*0.0
      # in case of zero layers we have to compensate 
      # for the non existing interface
      Heights[0]=-SeparationLower

      for nl in range(Npmin,Npmax+1):
         nep=nl                     # in composite planes models a layer is also en epaisseur
         dum=ones(max(1,nep))
         dum=dum*self.Plane.TotalThickness()
         upperthickness= self.Plane.getUpperThickness()

         Planes_n  = Planes_n+ [self.Plane] *nep
         Planes_Nmn  = Planes_Nmn+ [self.Plane] *nep
         for i in range (nep):
            l=i
            u=nep-1-i
            if(l>u): break
            if(l!=u):
               dum[l]=Planes_n  [i].TotalThickness()
               dum[u]=Planes_Nmn[i].TotalThickness()
            else:
               dum[l]=0.5*(Planes_n [i].TotalThickness()+
                                  Planes_Nmn[i].TotalThickness() )
               if(i==0):
                  upperthickness=0.5*(Planes_n [i].getUpperThickness()+
                                      Planes_Nmn[i].getUpperThickness() )



         # print "DUM ", dum 
         HeightsArrays[nl-Npmin]=concatenate( (
                                        (0, ),
                                         cumsum(dum)
                                        ))
         # print " HeightsArrays ", HeightsArrays

      Heights=arrayrange(Npmax+1-Npmin)*0.0
      # in case of zeo layers we have to compensate 
      # for the non existing interface
      Heights[0]=-SeparationLower
      for i in range (max(Npmin,1),Npmax+1):
        Heights[i-Npmin]=  (HeightsArrays[i-Npmin]) [-1]
        HeightsArrays[i-Npmin]=HeightsArrays[nl-Npmin][0:i] 
      ###############################################################
      # Calculation of average thickness and variance
      ###############################################################

      means= Probabilities *  Heights
      AveragedThickness=sum(means)
      means2=means*Heights

        
      ###
      # Things to remember
      ##
      self.Probabilities=Probabilities
      self.NlayersArray =NlayersArray
      self.HeightsArrays=HeightsArrays
      self.Heights=Heights
      self.AveragedThickness=AveragedThickness
      self.SeparationLower=SeparationLower
      self.SigmaSeparationLower=SigmaSeparationLower

      self.planesl=Planes_n
      self.planesu=Planes_Nmn

      # print " Averaged Thickness= ", AveragedThickness       

      # print " self.HeightsArrays " , self.HeightsArrays
      # print " self.Heights       ", self.Heights

  def ScatteringPower(self, Lambda, theta):

      iKperp=complex(0,1)*sin(theta)*4*math.pi/Lambda

      self.totalAbsorption =self.Plane.totalAbsorption( Lambda, theta)  

      #############################################
      # Scattering power of a plane
      #############################################
      f=self.Plane.ScatteringPower( Lambda, theta)  

      ##############################################################
      # Taking into account the interface diffusion if any 
      ##############################################################
      
      #  fld= ( self.planesl[0].ScatteringPower( Lambda, theta), self.planesl[1].ScatteringPower( Lambda, theta) )
      #  fud= ( self.planesu[0].ScatteringPower( Lambda, theta), self.planesu[1].ScatteringPower( Lambda, theta) )
 
      ############################################
      # for every thickness of the ensemble we
      # calculate the scattering powers of each layer 
      ############################################   
      scatterings={}
      for nl in  self.NlayersArray.tolist():
         scatterings[nl]=   outerproduct(  ones(nl, Complex), f)
         for l in range (0,max(  self.interf_Planes_L,   self.interf_Planes_H ) ):
            u=nl-l
            if (l>u): break
            if(l==u):
               scatterings[nl][l]=0.0
	       if(     l< self.interf_Planes_L):
                  scatterings[nl][l] += self.planesl[l].ScatteringPower( Lambda, theta)
               else:
                  scatterings[nl][l] += f
	       if(     u< self.interf_Planes_H):
                  scatterings[nl][u] += self.planesh[u].ScatteringPower( Lambda, theta)
               else:
                  scatterings[nl][l] += f
               scatterings[nl][l] *= 0.5
            else:
	       if(     l< self.interf_Planes_L):
                  scatterings[nl][l] = self.planesl[l].ScatteringPower( Lambda, theta)
	       if(     u< self.interf_Planes_H):
                  scatterings[nl][u] = self.planesh[u].ScatteringPower( Lambda, theta)

      #####################################################
      # ready to sum up
      #########################################################

      # print    scatterings         
      AvF=0        # Averaged F
      AvFF=0       # Averaged conj(F)*F
      AvPh=0       # Averaged Phase through thickness (< exp(iKperp*Thickness)>)
      AvPhF=0      # Averaged product of conj(F) X exp(iKperp*Thickness)
       
      for i in  range(0,len(self.Probabilities)):
         thickness=self.Heights[i]
         phase=  exp( iKperp*thickness )
         prob=self.Probabilities[i]
         AvPh =AvPh  + prob*phase

         nl=self.NlayersArray[i]
         if(nl==0): continue

         heights=self.HeightsArrays[i]

         scat=     sum( scatterings[nl]* TR(exp(outerproduct(iKperp,heights)) ))

         AvF  =AvF   + prob*scat
         AvFF =AvFF  + prob*scat*conjugate(scat)
         AvPhF=AvPhF + prob*conjugate(scat)*phase

      InterPh = exp( iKperp*self.SeparationLower +
                     iKperp*iKperp*self.SigmaSeparationLower*self.SigmaSeparationLower/2.0  )

      InterPhi = exp( iKperp*self.SeparationLower  )
      InterPhr = exp(  iKperp*iKperp*self.SigmaSeparationLower*self.SigmaSeparationLower/2.0  )

       

      (self.AvF, self.AvFF,
       self.AvPh, self.AvPhF,
       self.InterPh,
       self.InterPhi, self.InterPhr )=( AvF,AvFF,AvPh, AvPhF, InterPh,
                    InterPhi, InterPhr  )


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


class Fuller_PeriodM_I:
  def __init__(self):
    self.stack=[]
    self.classe="Period"

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

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

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

    elif(other.classe  in [  'LayerM_I',   'LayerM_IV' ,  'LayerM_XXI' ,   'Interface'] ):
       connect( self.stack[-1], other      ) 
       connect(other,  self.stack[0]       ) 
       newstack=self.stack+[other]

    else:
       print "Error summing a wrong class "
       return 0
    res=Fuller_PeriodM_I()
    res.stack=newstack
    return res


  def SetupCalculations(self):
    self.LayerList=[]
    for tok in self.stack:
      if(tok.classe in [ "LayerM_I","LayerM_IV",  'LayerM_XXI' ]):
          self.LayerList=self.LayerList+[tok]
          tok.SetupSharedArrays()
     


  def DoCalculation(self, Lambda, Theta, M=1, L=None, SigmaNs=0, Ds=0,
                     SubPlane=None, pol='S', sigmatheta=0, CRough=0, SubSep=0.0):
 
    self.SetupCalculations()
    Tprod=ones(len(Theta))
    eprod=ones(len(Theta))
    totabs=ones(len(Theta))


    for tok in self.LayerList:
       tok.ScatteringPower(Lambda,Theta)
       Tprod=Tprod*tok.AvPh
       eprod=eprod*tok.InterPh
       totabs=totabs*tok.totalAbsorption

 
    #################################################
    # The  geometric series give rise to
    # the factor g1 and g2
    ##################################################

    a=Tprod*eprod
    g1=a*(1-exp(log(a)*(M-1)))/(1-a)
    g2=a*(1  - M*exp( log(a)*(M-1) ) + (M-1)*exp( log(a)*M )   )/(1-a)/(1-a)

    intra=0.0
    extra=0.0
    if(SubPlane!=None): 
         fatavscat=1
         avscat=0.0  # averaged scattering. To be used with the scattering from substrate

    lista=self.LayerList
    for i in  range(0,len(lista) ):
      toki=lista[i]
      intra=intra+toki.AvFF
      if(SubPlane!=None):
         if(i): fatavscat=fatavscat*toki.InterPh
         avscat=avscat+toki.AvF*fatavscat
         fatavscat=fatavscat*toki.AvPh


      for j in range(i+1, len(lista) ):
         tokj=lista[j]
         fat=1
         for k in range(i+1,j): fat=fat*lista[k].AvPh
         for k in range(i+1  ,j+1): fat=fat*lista[k].InterPh
         #
         ## power from couples of the same layer
         #
         intra=intra+2*( fat*toki.AvPhF*tokj.AvF).real

      for j in range(0,len(lista)):
         fat=1
         tokj=lista[j]
         if(i+1<j):
             for k in range(i+1,j): fat=fat*lista[k].AvPh
         elif(i+1> j):
             for k in range(j,i+1): fat=fat/lista[k].AvPh
         if(i<j):
             for k in range(i+1  ,j+1): fat=fat*lista[k].InterPh
         elif(j<i):
             for k in range(j+1  ,i+1): fat=fat/lista[k].InterPh
         #
         ## power from couples of different layers
         #
         extra=extra+toki.AvPhF*tokj.AvF*fat


    if(SubPlane!=None):
         avscat=avscat*(1-exp(log(a)*M))/(1-a)

    q= sin(Theta)*4*math.pi/Lambda
    iKperp= complex(0.0,1.0)*q
    if(L!=None):
    #
    ## in plane averaging
    #
         # factorization of the sum 1->M
         g3=(1-exp(log(a)*M))/(1-a)
         g3=g3*conjugate(g3)
         extraL=0.0
         fr=1
         for toki in lista:
           fc=1
           for tokj in lista:
               extraL=extraL+ (toki.InterPh*toki.AvF*
                                 conjugate(tokj.InterPh*tokj.AvF)
                                *fr*fc)
               fc=fc*conjugate(tokj.InterPh*tokj.AvPh )  
           fr=fr*toki.InterPh*toki.AvPh

         xx3=xx2=xx1=0
         if(SigmaNs*SigmaNs>0.05): xx3 = exp( -0.5*(9.0/(SigmaNs*SigmaNs )) )
         if(SigmaNs*SigmaNs>0.01): xx2 = exp( -0.5*(4.0/(SigmaNs*SigmaNs )) )
         if(SigmaNs*SigmaNs>0.0025): xx1 = exp( -0.5*(1.0/(SigmaNs*SigmaNs )) )
         xxtot = 1. + 2.*(xx1+xx2+xx3)
         fat =xx3* exp(-3.*iKperp*Ds ) + xx3*exp(3.*iKperp*Ds )
         fat =fat +xx2*exp(-2.*iKperp*Ds ) + xx2*exp(2.*iKperp*Ds)
         fat =fat +xx1*exp(-iKperp*Ds  ) + xx1*exp(iKperp*Ds )
         fat = (fat + 1.)/xxtot

         if(SubPlane!=None):
            add1= ( L*(M*intra+2*( (M*g1-g2)*extra ).real)+ 
                    (L*(L-1)*fat*g3*extraL).real
                  )/L/L 


            subscat = SubPlane.ScatteringPower(Lambda,Theta)
            subscat = exp(-iKperp*SubSep)*subscat / (1- exp(-iKperp*Ds))

            add2=2*( avscat*
                     conjugate(subscat)
                   ).real+( subscat*conjugate(subscat)).real

            avscat=avscat+0.5*subscat
            add3=2*( avscat*conjugate(subscat)   ).real*fat
            result=  add1+(L*add2 +L*(L-1)*add3)/L/L    
         else:
           result=  ( L*(M*intra+2*( (M*g1-g2)*extra ).real)+ 
                    (L*(L-1)*fat*g3*extraL).real
                  )/L/L          
    else:
      if(SubPlane!=None):
         subscat = SubPlane.ScatteringPower(Lambda,Theta)
         subscat = exp(-iKperp*SubSep)*subscat / (1- exp(-iKperp*Ds))
         res= M*intra+2*( (M*g1-g2)*extra ).real
         res=res+2*( avscat*
                     conjugate(subscat)
                   ).real +( subscat*conjugate(subscat)).real
         result= res

      else:
         result= M*intra+2*( (M*g1-g2)*extra ).real         

#    totabs=1-exp(log(totabs)*M)
    totabs=1
    result=result/sin(2*Theta)*totabs 

    if(pol=='P'):
       cos2t=cos(2*Theta)
       cos2t2=cos2t*cos2t
       result=result*cos2t2
    if(pol=='M'):
       cos2t=cos(2*Theta)
       cos2t2=cos2t*cos2t
       result=result*(1+cos2t2)

    if(pol=='Nothing'):
       result=result*sin(2*Theta)


    if(sigmatheta!=0):
      convoluted=ones(len(result))*0.0
      deltatheta=Theta[1]-Theta[0]
      nmax= int(  5*sigmatheta/ deltatheta       )+1
      sum=0
      for shift in range(-nmax,nmax+1):
        x=shift*deltatheta
        x=x*x/2.0/sigmatheta/sigmatheta
        fact=exp(-x)
        sum=sum+fact  
        shift_add(convoluted,result.real,shift,fact)

      convoluted=convoluted/sum
      result=convoluted+complex(0.0,0.0)
    
    if(CRough!=0):
      result=result*exp( -q*q*CRough*CRough )
    return result.real

  def __repr__(self):
      print self.stack
      return ''




class Fuller_ComparisonTheoryExperiment:
  def __init__(self,stack, scanlist, weightlist, normlist=None, noise=None,
                 **args):
      self.args=args
      self.stack=stack
      self.scanlist=scanlist
      self.weightlist=weightlist
      if( normlist is None):
         self.normlist=[1,]*len(scanlist)
      else:
         self.normlist=normlist
      self.passes=0
      self.itercounter=0
      self.noise=noise

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

  def getitercounter(self):
      return self.itercounter


  def error(self):
     error=0
     count=0
      
     if(0 and self.passes==0):
       self.devices={}
       csc=0
       for scan in self.scanlist:
         self.devices[csc]=pgopen('/xw')
         csc+=1
         pgask(0) 
     self.passes+=1 


     if(0 and (self.passes % 10)==1 and self.passes>1    ): pgpage()

     csc=0
     numtot=0

     # '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
     #   stores  the calculated values for eventual retrieval in writefit2file
     #   
     # '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
     if( hasattr(self, "calculatedscan")==0): self.calculatedscan=[]
     newscans=[]

     for scan in self.scanlist:
        wavelenghts= scan[0]
        angles=scan[1]
        data=scan[2]

        rs= apply(self.stack.DoCalculation, ( wavelenghts, angles),
                    self.args
                   )
 
        if(self.noise!=None):
          rs=rs+par(self.noise[count])
        norm=par(self.normlist[count])

        newscans.append(rs*norm)
     
        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')


        difference=abs(log(norm*rs)-log(data))
        # difference=difference*difference*errorweights
        difference=difference*errorweights

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

        if(0 and (self.passes % 10)==1):
          pgslct(self.devices[csc])
          pgpage()
          maxI=max(scan[2])
          minI=min(scan[2])
          pgsci(5)
          pgenv(0.0,   len(scan[0] ) ,  log(minI),  log(maxI),  0,  1)
          pgsci(1)
          pgline(len(scan[0] ), arrayrange(0,len(scan[0] )), log(scan[2]) ) 
          pgsci(6)
      
          pgline(len(scan[0] ), arrayrange(0,len(scan[0] )), log(norm*rs) )       
        count=count+1
        csc+=1

     self.calculatedscan=newscans

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

class return_value:
   pass

def ScanReader(filename=None, Np="automatic", wavelenghts_col=1,energy_col=None, angles_col=2, refle_col=None, weight_col=None,
               angle_factor=1.0, TwoTheta=None):
   """ 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])
       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"), (energy_col, "energy_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()
 
   
   if( energy_col != None):
      wavelenghts_col=energy_col
   
   if(TwoTheta is None):
      optional=[(angles_col, "angles" )]
   else:
      optional=[]
      

   for (colonna_n, colonna_dati) in [ (wavelenghts_col, "wavelenghts" ),
                                     (refle_col, "refle" ),(weight_col, "weight" ) ] + optional:
    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"
    # print  getattr(res,colonna_dati, )


   if(TwoTheta!=None):
      Qfixe= 2*sin(TwoTheta*math.pi/180. / 2.0 ) * angles_col/ LAMBDAEV
      if( energy_col!=None):
         res.wavelenghts =  12398.52/res.wavelenghts      
      res.angles = arcsin(   Qfixe *res.wavelenghts /2.0)              
   else:
      if( energy_col!=None):
         res.wavelenghts =  12398.52/res.wavelenghts
      res.angles*=angle_factor

   # print " RETURN VALUE "
   # print [res.wavelenghts, res.angles, res.refle, res.weight]
   return [res.wavelenghts, res.angles, res.refle, res.weight]


def DW2sigma(DW):
   return math.sqrt( DW  /8.0/math.pi/math.pi)



def ScattererFromTable(tf0,tf12,Element,DWsigma=0,convention="+iKx"):
   f0   = tf0 .Element(Element)
   f1f2 = tf12.Element(Element)


   res =   Dabax_Scatterer([f0],  array([1.0]),   [f1f2]  ,  array([1.0])  , convention= convention, 
                       DWsigma= DWsigma  )
   return res

def SumThings(*args):
  res=None 
  for tok in args:
    if(res is None):
       res= tok
    else:
       res=res+tok
  res.SetupCalculations()
  return res

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







#######################################################################3
 



class Fuller_PerovskitePlane:
  def __init__(self,   Oxygen=None, External=None, Internal=None,
               Spacing=None, InPlaneDensity=1, PesiOxygen=[0.5,2.0,0.5]):

      self.Oxygen=Oxygen
      self.External=External
      self.Internal=Internal

      self.Spacing=Spacing

      self.classe="PerovskitePlane"

      self.PesiOxygen=PesiOxygen

      self.InPlaneDensity=InPlaneDensity

      self.upperThickness=0
      self.lowerThickness=0

      
  def setArrays(self):
    Spacing=par(self.Spacing)	


    self.Total= Spacing
    self.upperThickness=0
    self.lowerThickness=0

    self.O_array=array( (0.0,
                    Spacing/2.0, 
                    Spacing,
                   )
                 )

    self.Internal_array=array( (Spacing/2.0, 
		   )
                 )

    self.External_array=array( (0.0,
                    Spacing,
                   )
                 )

    self.O_weights        =  self.PesiOxygen
    self.Internal_weights =  array( (1.0)  )
    self.External_weights =  array( (0.5,0.5) )


  def totalAbsorption(self, Lambda, theta ):
   pesi=(sum(self.O_weights), sum(self.Internal_weights), sum(self.External_weights))
   atoms=(self.Oxygen, self.Internal, self. External)
#   res=0.0
#   for i in range(0,3):
#     res=res+(pesi[i]*atoms[i].DeltaBeta( Lambda, par(self.InPlaneDensity))).imag
   dens=self.InPlaneDensity
   E=LAMBDAEV/Lambda
   res=self.ScatteringPower(Lambda,  Lambda*0  )*415.22*dens/E/E/AVOGADRO*1.0e+24
   res=res.imag

     # dens is the volumic density (Natoms  A**-3)
     # return self.ScatteringPower(Lambda,  Lambda*0  )*415.22*dens/E/E/AVOGADRO*1.0e+24


   return exp(-2*res*2*math.pi/Lambda/sin(theta))



  def ScatteringPower(self, Lambda, theta):

   iKperp=complex(0,1)*sin(theta)*4*math.pi/Lambda
   #########################################################
   # ready to sum up
   #########################################################
   F=0 
            
   F=F+    ( self.Oxygen.F_Lambda( Lambda, theta) *  
                 sum( TR(self.O_weights * exp(outerproduct(iKperp,self.O_array)) ) ) )

   F=F+    ( self.Internal.F_Lambda( Lambda, theta) *  
                 sum( TR(self.Internal_weights * exp(    outerproduct(iKperp,self.Internal_array)    ) ) )) 

   F=F+    ( self.External.F_Lambda( Lambda, theta) *  
                 sum( TR(self.External_weights * exp(      outerproduct(iKperp,self.External_array  )  ) ) ) )

   return F*self.InPlaneDensity


  def TotalThickness(self):
    return self.Total 

  def getUpperThickness(self):
    return self.upperThickness

  def getLowerThickness(self):
    return self.lowerThickness

  def __repr__(self):
    print "    perovkite plane made of ",  self.Oxygen, self.Internal, self.External, "; InPlaneDensity=", \
          self.InPlaneDensity
    return ''
 
#######################################################################


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

class Fuller_WurtzitePlane:
  def __init__(self,   Atom1=None, Atom2=None,
               Spacing=None, InPlaneDensity=1):

      self.Atom1=Atom1
      self.Atom2=Atom2

      self.Spacing=Spacing

      self.classe="WurtzitePlane"

      self.InPlaneDensity=InPlaneDensity

      self.upperThickness=0
      self.lowerThickness=0

      
  def setArrays(self):
    Spacing=par(self.Spacing)	


    self.Total= Spacing
    self.upperThickness=0
    self.lowerThickness=0

    self.Atom1_array=array( (0.0,
                    Spacing/2.0, 
                    Spacing, 
                   )
                 )

    self.Atom2_array=array( (0.375*Spacing,
                             0.875*Spacing, 
                   )
                 )


    self.Atom1_weights =  array( (0.5,1.0,0.5)  )
    self.Atom2_weights =  array( (1.0,1.0) )


  def totalAbsorption(self, Lambda, theta ):
   pesi=(sum(self.Atom1_weights), sum(self.Atom2_weights))
   atoms=(self.Atom1, self.Atom2)
#   res=0.0
#   for i in range(0,3):
#     res=res+(pesi[i]*atoms[i].DeltaBeta( Lambda, par(self.InPlaneDensity))).imag
   dens=self.InPlaneDensity
   E=LAMBDAEV/Lambda
   res=self.ScatteringPower(Lambda,  Lambda*0  )*415.22*dens/E/E/AVOGADRO*1.0e+24
   res=res.imag

     # dens is the volumic density (Natoms  A**-3)
     # return self.ScatteringPower(Lambda,  Lambda*0  )*415.22*dens/E/E/AVOGADRO*1.0e+24


   return exp(-2*res*2*math.pi/Lambda/sin(theta))



  def ScatteringPower(self, Lambda, theta):

   iKperp=complex(0,1)*sin(theta)*4*math.pi/Lambda
   #########################################################
   # ready to sum up
   #########################################################
   F=0 
            
   F=F+    ( self.Atom1.F_Lambda( Lambda, theta) *  
                 sum( TR(self.Atom1_weights * exp(outerproduct(iKperp,self.Atom1_array)) ) ) )

   F=F+    ( self.Atom2.F_Lambda( Lambda, theta) *  
                 sum( TR(self.Atom2_weights * exp(    outerproduct(iKperp,self.Atom2_array)    ) ) )) 

   return F*self.InPlaneDensity


  def TotalThickness(self):
    return self.Total 

  def getUpperThickness(self):
    return self.upperThickness

  def getLowerThickness(self):
    return self.lowerThickness

  def __repr__(self):
    print "    wurtzite plane made of ",  self.Atom1, self.Atom2, "; InPlaneDensity=", \
          self.InPlaneDensity
    return ''
 
#######################################################################




#######################################################################3



if(__name__=='__main__'):

  print LAMBDAEV
  print dir("__main__")

#############################################################
#  Read the tables
#
#  Table_f0        = Dabax_f0_Table("f0_WaasKirf.dat")
#  Table_f1f2      = Dabax_f1f2_Table("f1f2_Windt.dat")
#
#

  con1=convenient()
  con1.value0=0.01
  con1.value12=complex(0.0,0.01)
  Mocon=Dabax_Scatterer([con1], [1],[con1],[1])
  Nicon=Dabax_Scatterer([con1], [2],[con1],[1])

######################################################333
#
#  Compose the elements
#
#  Mo_f0     = Table_f0     .Element("Mo")
#  Mo_f1f2   = Table_f1f2   .Element("Mo")

#  Mo        = Dabax_Scatterer((Mo_f0,),(1,), (  Mo_f1f2,),(1,)       )

#  Ni_f0     = Table_f0    .Element("Ni")
#  Ni_f1f2   = Table_f1f2   .Element("Ni")

#  Ni        = Dabax_Scatterer((Ni_f0,),(1,), (  Ni_f1f2,),(1,)       )



######################################################
#
#  define  the layers
#


  Mo_Plane  = Fuller_ElementalPlane(Element=Mocon,  InPlaneDensity=1)
  Ni_Plane  = Fuller_ElementalPlane(Element=Nicon,  InPlaneDensity=1)

  Mo_Layer  = Fuller_LayerM_I(Plane=Mo_Plane, PlaneD = 2.2,  Nplanes=7.5    , SigmaNplanes=0.3,
                              PlaneDD1=-0.0, PlaneDD2= 0.0, PlaneDDAlpha=0.5)
  Ni_Layer  = Fuller_LayerM_I(Plane=Ni_Plane, PlaneD = 2.,  Nplanes=7.00   , SigmaNplanes=0.3 ,
                              PlaneDD1= 0.0, PlaneDD2= 0.0, PlaneDDAlpha=0.5)


  

  period    = Mo_Layer+Ni_Layer



  print "start calculations"
  f=open("curva","w")
  thetalist=array([ theta for theta in arrayrange(5,25,0.0133)*math.pi/180.0 ])
  lambdalist=array( [1.546 for theta in arrayrange(5,25,0.0133)*math.pi/180.0 ])
  sin(thetalist)
  res=(period .DoCalculation( lambdalist , thetalist , 60  , L=None, SigmaNs=0, Ds=0, SubPlane=None , sigmatheta=0.0,
                              pol='M', CRough=0 ))

  
  thetalist=thetalist*2.0
  for i in range (0,len(thetalist)):
    f.write("%e %e\n"% (thetalist[i]*180.0/math.pi,res[i]  ))


  print "finished  calculations"













