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

from Dabax import *

from PyFullerton  import *






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

Table_f0= Dabax_f0_Table("f0_WaasKirf.dat")
# Table_f1f2=Dabax_f1f2_Table("f1f2_Windt.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_Layer  = Fuller_LayerM_I(Plane=Or_Plane, PlaneD = 2.3711,  Nplanes=11.7958 , SigmaNplanes=0.9063,
                              PlaneDD1=0, PlaneDD2= 0, PlaneDDAlpha=0.5)

Ni_Layer  = Fuller_LayerM_I( Plane=Ni_Plane,PlaneD = 2.0665    ,  Nplanes=5.0959,SigmaNplanes=0.0646 ,
                              PlaneDD1= 0.0, PlaneDD2= 0.0, PlaneDDAlpha=0.5)




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


period .SetupCalculations()






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

Qfixe= 2*sin(3.815000e+01*math.pi/180 / 2.0 ) * 8815.0/ 12398.52

lambdalist= 12398.52/ arrayrange(8200,8801,1.0)       

thetalist = arcsin(   Qfixe *lambdalist /2.0)              


res=(period .DoCalculation( lambdalist , thetalist ,21 , L=None, SigmaNs=0, Ds=0, SubPlane=None , sigmatheta=0.0,
                              pol='Nothing', CRough=0 ))




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 %e\n"% (thetalist[i]*180.0/math.pi, 12398.52/lambdalist[i], res[i]  ) )

print "finished  calculations"














