from sys import *
from copy import *
import NumericA
from Numeric import *
from  arrayfns  import *            
import math
import time
# from umathextra import *

import  Minimiser 
from Dabax import *
# from dabax import *

from PyFullerton  import *


LAMBDAEV=12398.52
DeuxTheta=38.15


#######################################################################
# specify here the scatterers
#

Table_f0= Dabax_f0_Table("f0_WaasKirf.dat")
Table_f1f2=Dabax_f1f2_Table("f1f2_Sasaki.dat")

Au_f0  = Table_f0.Element("Au")
Au_f1f2 = Table_f1f2.Element("Au")

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

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


Or =   Dabax_Scatterer([Au_f0],  [1.0],   [Au_f1f2]  ,  [1.0]  , convention="-iKx" , 
                       DWsigma=DW2sigma(0.5)    )
Ni=    Dabax_Scatterer([Ni_f0 ], [1.0],   [Ni_f1f2 ]  ,  [1.0] , convention="-iKx" , 
                        DWsigma= DW2sigma(0.23)  )


print Or.F_Energy(8815.0,0.0)
print Ni.F_Energy(8815.0,0.0)


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

Or_Plane  = Fuller_ElementalPlane(Element=Or,  InPlaneDensity=1*0.139)
Ni_Plane  = Fuller_ElementalPlane(Element=Ni,  InPlaneDensity=0.91*0.186)


Or_PlaneD=Variable(2.3711, 1.0,3.0)
Or_Nplanes=Variable(11.7958, 1.0,13.0  )
Or_SigmaNplanes=Variable(0.9063,0.2 , 1.8)

#########################################################################
# la on cree comme d'habitude une couche avec ses parametres
# Avant on mettait des chiffres partout. Maintenant ce que on veut varier
# ce sont des objects avec la fonctions getvalue. Ils sont appeller par getvalue
# si ils ont cette function si non, si ils sont des chiffres, on utilise
# la valeur numerique 
#
#

Or_Layer  = Fuller_LayerM_I(Plane=Or_Plane, PlaneD = Or_PlaneD ,  Nplanes= Or_Nplanes , SigmaNplanes=Or_SigmaNplanes,
                              PlaneDD1=0, PlaneDD2= 0, PlaneDDAlpha=0.5)

Ni_PlaneD=Variable(2.0665 , 1.0,3.0)
Ni_Nplanes=Variable(5.095, 1.0,13.0  )
Ni_SigmaNplanes=Variable(0.0646,0.00001 , 1.8)


Ni_Layer  = Fuller_LayerM_I( Plane=Ni_Plane,PlaneD = Ni_PlaneD ,  Nplanes=Ni_Nplanes,SigmaNplanes= Ni_SigmaNplanes,
                              PlaneDD1= 0.0, PlaneDD2= 0.0, PlaneDDAlpha=0.5)




SigSep=Variable(0.14,0.001,2)

period    = (Fuller_Interface(Separation=None, SigmaSeparation=SigSep, StearnsWidth=0.00)+
             Or_Layer+
             Fuller_Interface(Separation=None, SigmaSeparation=SigSep, StearnsWidth=4.6)+
             Ni_Layer
            )


period .SetupCalculations()






# thetalist=array([ theta for theta in arrayrange(5,25,0.0133)*math.pi/180.0 ])

# thetalist = arrayrange(22,43.25,0.05)*0.5*math.pi/180.0
# lambdalist= 1.4065*ones(len(thetalist) ) 
# res=(period .DoCalculation( lambdalist , thetalist ,M=21 , L=None, SigmaNs=0, Ds=0, SubPlane=None , sigmatheta=0.0,
#                              pol='Nothing', CRough=0 ))

fd=open("dafs_p1.txt","r")
nwave=string.atoi(fd.readline())

wavelenghts1=zeros(nwave,Float)
angles1     =zeros(nwave,Float)
refle1      =zeros(nwave,Float)

for i in range(0, nwave ):
  line=fd.readline()
  line=string.split(line)
  wavelenghts1[i]= string.atof(line[0])
  wavelenghts1[i]=LAMBDAEV/wavelenghts1[i]
  refle1 [i]=string.atof(line[1])
  


Qfixe= 2*sin(DeuxTheta*math.pi/180 / 2.0 ) * 8815.0/ LAMBDAEV
thetalist = arcsin(   Qfixe *wavelenghts1 /2.0)              




print " wavelenghts1 ", wavelenghts1
print " thetalist ", thetalist
print " refle1 ", refle1


scanlist=[[wavelenghts1, thetalist,refle1 ]]
weightlist=[1.0]


norm=Variable(0.0001,0.00001,0.001   )

fit=Fuller_ComparisonTheoryExperiment(period, scanlist, weightlist,   normlist=[norm], noise=None,
              M=21 , L=None, SigmaNs=0, Ds=0, SubPlane=None , sigmatheta=0.0,
                              pol='Nothing', CRough=0   )



ListeVariables= [Ni_PlaneD, Ni_Nplanes , Ni_SigmaNplanes, Or_PlaneD, Or_Nplanes,
                 Or_SigmaNplanes, SigSep, norm]


mini=Minimiser.minimiser(fit, ListeVariables )
mini.amoebaAnnealed( 0.001, lambda x: .2/(1.0+x*0.3), 100,
                               max_isthesame=10, arret=1.0e-6)



print dir(Ni_Layer)
print Ni_Layer.scatterings[3][ : ,0   ]
print Ni_Layer.HeightsArrays[3]
print " ################################# "
print Or_Layer.scatterings[3][ : ,0   ]
print Or_Layer.HeightsArrays[3]
print " ################################# "

hgt=Ni_Layer.HeightsArrays[3]
print hgt[1:]-hgt[0:-1]

hgt=Or_Layer.HeightsArrays[3]
print hgt[1:]-hgt[0:-1]




# print Ni_Layer.scatterings

res=res*0.12 +44
thetalist=thetalist*2.0

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

print "finished  calculations"

time.sleep (10)














