#!/usr/bin/python

from scipy.special import j1 , hyp1f1
from scipy.integrate import simps
from scipy.stats import norm , gamma
import numpy as np
import sys
import math
pi=math.pi
e=math.e

sys.path.append("/mntdirect/_users/semeraro/python_tools") 
import FitTools

N_OrAv = 90*2
N_CapCyl = 50

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

""" SPHERICAL BESSEL FUNCTION j0 & j1 	""" 	

def spherical_j0_np(x):
	return np.where(x==0, 1, np.sin(x)/x )	
	
def spherical_j1_np(x):
	return np.where(x==0, 1./3, (np.sin(x)-x*np.cos(x))/x**2 )	

###########################################################################################
###########################################################################################
""" NORMAL DISTRIBUTION 	""" 	
def Normal(x,xo,w):
	return np.exp(- (x-xo)**2/(2*w**2) ) / (w*np.sqrt(2*pi))

###########################################################################################
###########################################################################################
"""  FORM FACTORS 	""" 

################################	SPHERE
def Amp_Sphere(q,R):
	if R==0 or q==0:
		A=Vol_Sphere(R)
	else:
		A=Vol_Sphere(R) * 3 * spherical_j1(q*R) / (q*R)
	return A

################################	SPHERES 
def Sphere( q, PAR ):
	[n,rho,R,C] = PAR
					
	amp = (4*np.pi*R**3/3) * 3 * spherical_j1_np(q*R) / (q*R)
	
	return n * (rho * amp)**2 + C
	
################################	POLYDISPERSE SPHERES (NORMAl PDF)
def Sphere_Normal (q,PAR):
	[n,rho,R,w,C] = PAR

	N = 21
	R_array = np.linspace(R-3*w,R+3*w,N)
	Normal = norm.pdf(R_array, R, w)
				
	I = np.zeros(q.shape[0], dtype=float)
			
	for r in range(R_array.shape[0]):
		if r==0 or r==N-1 : I+= Sphere(q,[1,1,R_array[r],0]) * Normal[r] / 2
		else : 				I+= Sphere(q,[1,1,R_array[r],0]) * Normal[r]
	I*= 6*w/(N-1)
	
	return n*rho**2*I + C