#-*- coding: utf-8 -*-
# Comparação do espectros de espalhamento baseados na função modificada de Bessel do primeiro tipo de ordem 1 
# com o modelo Gaussiano. 
# Incluída convolução com forma de linha originária da colisão frontal
#
# Referência: R. P. Pezzi, et al. Applied Physics Letters, 92 164102 (2008)

from pylab import *
from mpmath import *

PI=math.pi

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


#### Parâmetros iniciais ####
# Feixe:
E0=100000 #eV
mi=1.
Zi=1.

# Alvo:
dedx=192. # eV/nm
dw2dx=20740. # eV^2/nm
Theta_in=0.
Theta_out=70.
Theta_s = PI - (Theta_in + Theta_out)*PI/180
espessura=2

# Cálculo:
deltaE=20000
passoE=20
passox=0.02
e = arange(E0-deltaE, E0, passoE)
besselcamada= arange(E0-deltaE, E0, passoE)*0

## TODO: - Otimizar os cálculos de cada componente: evitar cálculos (de zeros) em energias desnecessários
##       - Definir o intervalo de plotagem do espectro independente do deltaE
#############################


## Define a composição da amostra ##
elementos=[[72.,178.49,217.],[14.,28.,121.],[8.,16.,102]] # Hf, Si e O
#elementos=[[72.,178.49,217.]]#,[14.,28.,121.],[8.,16.,102]] # Hf, Si e O

# Esta lista deve incluir a distribuição
# em profundidade do elemento 
####################################

mt=178
Zt=72

k = ((sqrt (mt**2+mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2	
dedxef=dedx*(k/cos(Theta_in*PI/180.) + 1/cos(Theta_out*PI/180.))
dw2dxef=dw2dx*(k**2./cos(Theta_in*PI/180.)+1./cos(Theta_out*PI/180.))



alpha=dedxef*(2./dw2dxef)
m=alpha*dedxef
lbd = limit(lambda x: m*x*alpha, 0)
#bespecial = besseli(1, 2.*sqrt(lbd*(k*E0-e)) )
#i1(2.*sqrt(lbd*(k*E0-e)))
print lbd
for i in e:
    besselcamada[int(i)] = lbd*exp(-lbd/alpha-alpha*(k*E0-int(i)))*besseli(1, 2.*sqrt(lbd*(k*E0-int(i))) ) /(sqrt(lbd*(k*E0-int(i)))) 

       ### Removendo Not A Number do vetor ###
ondeestaoNaN = isnan(besselcamada);
besselcamada[ondeestaoNaN] = 0

plot(e, besselcamada)

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

#print totalbessel
show() 
##################