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

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



import Numeric
import arrayfns
from Numeric import *
import PPMcore
from Dabax import *
from  arrayfns  import *            
import string
import Minimiser
from Dabax import IndexFromTable
import os

def   shift_add(convoluted,rs,shift,fact):
  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):
  if('getvalue' in dir(x)):
    return x.getvalue()
  else:
    return x

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

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

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


class PPM_Stack:
  interpola=0
  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 calculatescanParallel(self,wavelenghts, angles, nscan):
      if( (  mpi.size % nscan  ) != 0 ):
          print " lenscans = ",nscan , "  mpi.size = ", mpi.size
          raise " PROBLEM (mpi.size  % nscan) != 0  in calculatescanParallel "

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

      nself   = mpi.rank / nscan

      lscan = len(wavelenghts)

      rs=Numeric.zeros(lscan,"d")
      rp=Numeric.zeros(lscan,"d")

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

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

      partialS, partialP = self.calculatescan(wavelenghts[inizio:fine] , angles[inizio:fine]   , 1)
      for i in range(nprocs):
          iproc = i*nscan+nmaster
          if( mpi.rank== nmaster):
              if(i==0):
                  receivedS = partialS
                  receivedP = partialP
              else:
                  (receivedS, receivedP), status = mpi.recv(iproc , mpi.ANY_TAG  )
              inizio = i*lentoken
              fine   = (i+1)*lentoken
              if( fine > lscan ):
                  fine=lscan
              rs[inizio:fine]= receivedS
              rp[inizio:fine]= receivedP
          else:
              if( i == nself ):
                  mpi.send((partialS,partialP), nmaster) 
      return (rs,rp)
  
  def calculatescan(self,wavelenghts, angles, Gives_Intensity=1):
    thicknesslist=[]
    roughnesslist=[]
    indexlist=[]
     

    depth_count = 0
    to_th=0.0
    for layer in self.stack:
        Minimiser.DependentVariable.depth = depth_count
        Minimiser.DependentVariable.to_th = to_th
        
        layerdescr= layer.getThickRoughIndexes( wavelenghts)
        thicknesslist.append( layerdescr[0])
        roughnesslist.append( layerdescr[1])
        indexlist    .append( layerdescr[2])
        depth_count = depth_count + 1
        to_th = to_th + Numeric.sum(to_th)
        
    Minimiser.DependentVariable.depth = 0
    Minimiser.DependentVariable.to_th = 0

    Thickness= concatenate(thicknesslist)
    Roughness= concatenate(roughnesslist)
    Index    = concatenate(indexlist    )


    if(self.interpola):
         print " EFFETTUO INTERPOLAZIONE "
         Thickness, Roughness, Index    = Interpola(  Thickness, Roughness, Index , self.interpola )
         
    
    SubRough=Roughness[0]*1.0
    SubIndex=Index    [0]
    
    # print "XXXXXXXXXXXXXXXXXXXXXXX"
    # print Thickness
    # print Roughness

    Thickness=Thickness[1:]*1.0
    Roughness=Roughness[1:]*1.0
    Index    =Index    [1:]*1.0


    


#    print "SubRough ", SubRough
#    print "Thickness ", Thickness

    (a,b)=PPMcore.PPM_calculatescan(Thickness,Roughness, wavelenghts,
                                angles,Index,SubIndex,SubRough)


    if( Gives_Intensity==1) :
	    rs= (a*conjugate(a)).real
	    rp= (b*conjugate(b)).real
    else:
	    rs= a
	    rp= b
   
    return (rs,rp)


  def calculateFields(self,wavelenghts, angles):

    thicknesslist=[]
    roughnesslist=[]
    indexlist=[]
     
    for layer in self.stack:
        layerdescr= layer.getThickRoughIndexes( wavelenghts)
        thicknesslist.append( layerdescr[0])
        roughnesslist.append( layerdescr[1])
        indexlist    .append( layerdescr[2])

    Thickness= concatenate(thicknesslist)
    Roughness= concatenate(roughnesslist)
    Index    = concatenate(indexlist    )
    
    SubRough=Roughness[0]*1.0
    SubIndex=Index    [0]
    
    # print "XXXXXXXXXXXXXXXXXXXXXXX"
    # print Thickness
    # print Roughness

    Thickness=Thickness[1:]*1.0
    Roughness=Roughness[1:]*1.0
    Index    =Index    [1:]*1.0


#    print "SubRough ", SubRough
#    print "Thickness ", Thickness

    (S,P)=PPMcore.PPM_calculatescan(Thickness,Roughness, wavelenghts,
                                angles,Index,SubIndex,SubRough)

   
    return (S,P)



  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:
  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 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
 

class PPM_ElementalLayer(PPM_LayerBaseClasse):
  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):
      indexes=zeros(shape(wavelenght),Complex)
      for i in range(0,len(self.DabaxList)):
         scatterer=self.DabaxList[i]
         density  =par(self.DensityList[i])
         mass     =self.MassList[i]



         if( scatterer.is_an_index ):
           add = 1.0-scatterer.index(wavelenght)
         else:
           if( scatterer.__class__ == dabax.return_value ):
             add= 1- CalculateIndexesFromTable( scatterer,wavelenght )
           else:
             add=415.22*(  scatterer.F_Lambda(wavelenght) )*( 
               wavelenght*wavelenght/12398.52/12398.52/mass)

         indexes=indexes+ add*density
 

      return ( array(( par(self.thickness),)),
               array(( par(self.roughness),)),
               array((1-indexes  ,)) 
             ) 

      
def Write_Stack(stack,outn , energy):
  if( ISMASTER!=1) :
       return
  count=0


  for l in stack.stack:
    Minimiser.DependentVariable.depth = count
    thickness = l.thickness
    roughness = l.roughness
    res=dabax.return_value()
    res.DabaxList  = l.DabaxList
    res.DensityList= l.DensityList
    res.MassList   = l.MassList

    obj = KK_ForPlot(res,energy,energy ,1 )
    obj.error()
    if count==0:
      mode="w"
    else:
      mode="a"
    obj.write2File(outn, (par(thickness), par(roughness)), mode)

    count=count+1
  Minimiser.DependentVariable.depth = 0



class PPM_SimpleLayer(PPM_ElementalLayer):
  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):
  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):
  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:
  classe="PPM_ComparisonTheoryExperiment"
  def __init__(self,stacks,scan2stack, scanlist, weightlist, normlist=None, noise=None, width=0, printpartial=0, shiftlist=None, CutOffRatio=None, angleshift=None):
      self.CutOffRatio = CutOffRatio
      self.width=width
      self.angleshift=angleshift
      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 write2FileNewPlot(self, filename):
     if( ISMASTER!=1) :
            return

     for count in range(len(self.calculatedscan ) ) :
       # print " stampo scan ", count
       calculated=self.calculatedscan[count]
       scan = self.scanlist[count]
       f=open("%s%d" %(filename, count+1),"w")
       f.write("\n")
       nscanspec=1
       

       if(len(scan)==4):
         expe=1
         scanstring="#N 5\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
              
         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 %e\n"%(calculated[k], scan[0][k], scan[1][k],scan[2][k],scan[3][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))
       if(len(scan)==3):
         for k in range(len(calculated)):
           f.write("%e %e %e %e\n"%(calculated[k], scan[0][k], scan[1][k],scan[2][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 write2FileASCII(self, filename):
     if( ISMASTER!=1) :
            return

     for count in range(len(self.calculatedscan ) ) :
       # print " stampo scan ", count
       calculated=self.calculatedscan[count]
       scan = self.scanlist[count]
       f=open("%s%d" %(filename, count+1),"w")


       if(len(scan)==4):
         for k in range(len(calculated)):
           try: 
             f.write("%e %e %e %e %e\n"%(calculated[k], scan[0][k], scan[1][k],scan[2][k],scan[3][k]  ))
           except:
             f.write("%e %e %e\n"%(calculated[k], scan[0][k], scan[1][k]  ))           
       if(len(scan)==3):
         for k in range(len(calculated)):
           f.write("%e %e %e %e\n"%(calculated[k], scan[0][k], scan[1][k],scan[2][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 GetResForGeFit(self):
       res=[]
       for scan in self.calculatedscan:
            res.append(scan)
       res=Numeric.concatenate(res)
       print res.shape

       return Numeric.log(res)
  
  def GetDataForGeFit(self):
       res=[]
       for scan in self.scanlist:
            res.append(scan[2])
       res=Numeric.concatenate(res)
       return Numeric.log(res)
  
  def writeminimum(self, iiter, res ):
     if( ISMASTER and type(self.printpartial)==type([]) ):
         vs=self.printpartial[0]
         ns=self.printpartial[1]
         f=open("PARTIAL/RESULTA_AT_%d" % iiter  , "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.write("errore= %e\n" %res)
         f.close()




  def error(self, nopartial=0):

     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
         
         print 'DEBUG PPMlayerclass', os.name
         if os.name == 'posix':
                cmd="ln -fs variables%d  PARTIAL/lastvariables " % ( ( self.itercounter % 30 )  )  
                os.system(cmd )
         elif os.name == 'nt':
            pass
         else:
            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
         data=scan[2]
         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]
           if(self.angleshift is not None):
             shift=par(self.angleshift[count])
             angles=scan[1]+shift


           if(len(scan)==4):
             errorweights=scan[3]
           else:
             errorweights=ones(len(data ))
           if( ISPARALLEL ):
               (rs,rp)= self.stacks[self.scan2stack[count]].calculatescanParallel(wavelenghts,angles, lenscans )
           else:
               (rs,rp)= self.stacks[self.scan2stack[count]].calculatescan(wavelenghts,angles)
           if(self.noise!=None):
             rs=rs+par(self.noise[count])
             rp=rp+par(self.noise[count])

           if(self.width!=0):
             # print " CONVOLUTION"
             # cerca i punti di stacco fra gli scan

             convoluted=ones(len(rs))*0.0
             sumcon    =ones(len(rs))*0.0
             ones_     =ones(len(rs))
             frontiera =ones(len(rs))
           
             for ist in range(len(rs)-1):
                  if( wavelenghts[ist]!=wavelenghts[ist+1] and angles[ist]!=angles[ist+1] ):
                       for fi in range(ist- 5*self.width  , ist+5*self.width+1):
                            if(fi>=0 and fi<len(rs)):
                                 frontiera[fi]=0

 
             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)
               if(shift==0):
                    shift_add(convoluted,rs,shift,fact)
                    shift_add(sumcon,ones_,shift,fact)
               else:
                    shift_add(convoluted,rs*frontiera,shift,fact)
                    shift_add(sumcon, ones_*frontiera,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)
         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)

     self.calculatedscan=newscans
     if( ISMASTER):
      if(self.printpartial and nopartial==0):
        self.write2File("PARTIAL/partialresults%d_" % ( self.itercounter % 30 )   )
           #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 )
        print 'DEBUG PPMlayerclass', os.name
        if os.name=='posix':
            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 )
                 print 'DEBUG ', cmd
        elif os.name=='nt':
            pass
        else:
             pass


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

if __name__=='__main__':

  NWAVE=100

  Table_f0= Dabax_f0_Table("f0_WaasKirf.dat")

  Table_f1f2=Dabax_f1f2_Table("f1f2_Windt.dat")


  Sif0  = Table_f0.Element("Si")
  Sif1f2= Table_f1f2.Element("Si")
  Si=Dabax_Scatterer((Sif0,),(1,), (  Sif1f2,),(1,)       )


  Cof0  = Table_f0.Element("Co")
  Cof1f2= Table_f1f2.Element("Co")
  Co=Dabax_Scatterer((Cof0,),(1,), (  Cof1f2,),(1,)       )

 

  Cuf0  = Table_f0.Element("Cu")
  Cuf1f2= Table_f1f2.Element("Cu")
  Cu=Dabax_Scatterer((Cuf0,),(1,), (  Cuf1f2,),(1,)       )

  SiLayer= PPM_ElementalLayer( thickness=0, roughness=0, 
               DabaxList=[Si,], DensityList=[2.3,], MassList=[30.0,])

  CoLayer= PPM_ElementalLayer( thickness=10, roughness=0, 
               DabaxList=[Co,], DensityList=[9.0,], MassList=[56.0,])

  CuLayer= PPM_ElementalLayer( thickness=30, roughness=0, 
               DabaxList=[Cu,], DensityList=[9.0,], MassList=[56.0,])

  
  stacks=[PPM_Stack()]

  stacks[0]=stack[0]+SiLayer

  for i in range(0,14):
     stack[0]=stack[0]+CoLayer+CuLayer


  print stack[0]


  wavelenght= ones(NWAVE)*1.0
  angles    = arrayrange(0.01 ,0.01+NWAVE*0.002,  0.002)*arccos(-1.0)/180.0

  for i in range (0,20):
    (rs,rp)=stack[0].calculatescan(wavelenght,angles)


  f=open("refle","w")
  for i in range(0, len(rs)):
    f.write("%e %e %e\n" %(angles[i]*180.0/arccos(-1.0), rs[i],rp[i]) )

class return_value:
  pass

def ScanReader(filename=None, Np="automatic", wavelenghts_col=1,angles_col=2, 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.
   """  


   print filename

   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))):
      try:
        setattr(res,colonna_dati, data[:,0]*0.0+colonna_n)
      except:
        print " FFF ", filename
        raise " CACA "
    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, res.weight]


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()

def Interpola(  Thickness, Roughness, Index , interpola ):
     print "len(Thickness)"
     print len(Thickness)
     print "len(Index)"
     print len(Index)
     if( interpola% 2 ==0 ):
          raise " interpola deve essere dispari in Interpola(  Thickness, Roughness, Index , interpola ) "
     x=0
     POSlist=[]
     for t in Thickness:
          x=x+t/2.0
          POSlist.append(x)
          x=x+t/2.0
     # per il vuoto
     #POSlist.append(x)

     NEWthickness = [0.]
     for t in Thickness[1:]:
          NEWthickness=NEWthickness+[ t/interpola ]*interpola
     
     NEWposList=[]
     x=0.0
     for t in NEWthickness:
        x=x+t/2.0
        NEWposList.append(x)
        x=x+t/2.0
     # per il vuoto
     # NEWposList.append(x)

     indreal = Numeric.array(Index).real
     indimag = Numeric.array(Index).imag
     print "len(indreal)"
     print len(indreal)
     print "len(POSlist)"
     print len(POSlist)

     shapeind = indreal.shape
     shapeNewInd = [  len(NEWposList), shapeind[1]]
     NEWindreal  = Numeric.zeros(shapeNewInd, 'd')
     NEWindimag  = Numeric.zeros(shapeNewInd, 'd')
     for ilambda in range(shapeind[1]):
          NEWindreal[:,ilambda ] = arrayfns.interp(indreal[:,ilambda ]
                                                   ,POSlist , NEWposList  )
          NEWindimag[:,ilambda ] = arrayfns.interp( indimag[:,ilambda ]
                                                    ,POSlist , NEWposList  )

     NEWrough = [0.0]* len(NEWthickness)

     return Numeric.array(NEWthickness) ,  Numeric.array(NEWrough), NEWindreal + 1.0j*NEWindimag
