Commit f4f01b27 authored by Remi Perenon's avatar Remi Perenon
Browse files

Implementing compute diffraction pattern with lowlevel API

parent 2c230bdb
Pipeline #8060 passed with stages
in 12 minutes and 35 seconds
......@@ -2,24 +2,29 @@ import sys
import os
sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), os.pardir)))
import matplotlib.pyplot as plt
import CFML_api
filename = sys.argv[1]
cif_file = CFML_api.CIFFile(filename)
Cell = cif_file.cell
SpG = cif_file.space_group
A = cif_file.atom_list
cell = cif_file.cell
space_group = cif_file.space_group
atom_list = cif_file.atom_list
Cell.print_description()
SpG.print_description()
A.print_description()
cell.print_description()
space_group.print_description()
atom_list.print_description()
powder_pattern_simulation_conditions = CFML_api.PowderPatternSimulationConditions()
ReflectionList = CFML_api.ReflectionList(
Cell, SpG, True, 0, powder_pattern_simulation_conditions.getSinThetaOverLambdaMax())
ReflectionList.compute_structure_factors(SpG, A)
reflection_list = CFML_api.ReflectionList(
cell, space_group, True, 0, powder_pattern_simulation_conditions.getSinThetaOverLambdaMax())
reflection_list.compute_structure_factors(space_group, atom_list)
diffraction_pattern = CFML_api.DiffractionPattern(powder_pattern_simulation_conditions, reflection_list, cell.get_reciprocal_cell_vol())
Pat = CFML_api.DiffractionPattern(powder_pattern_simulation_conditions, ReflectionList)
plt.plot(diffraction_pattern.x, diffraction_pattern.y)
plt.show()
\ No newline at end of file
......@@ -16,9 +16,18 @@ import CFML_api.crysfml_api
import CFML_api.FortranBindedClass
class DiffractionPattern(CFML_api.FortranBindedClass):
def __init__(self, simulation_conditions, reflection_list):
def __init__(self, simulation_conditions, reflection_list, scale_f):
self._set_fortran_address(CFML_api.crysfml_api.diffraction_patterns_compute_powder_pattern(
simulation_conditions, reflection_list))
simulation_conditions.to_dict(), reflection_list.get_fortran_address(), scale_f)["address"])
def get_x(self):
return CFML_api.crysfml_api.diffraction_patterns_get_x(self.get_fortran_address())["x"]
def get_y(self):
return CFML_api.crysfml_api.diffraction_patterns_get_y(self.get_fortran_address())["y"]
x = property(get_x)
y = property(get_y)
class PseudoVoigt():
......@@ -54,7 +63,7 @@ class PseudoVoigt():
ag= 0.93943727869965133377234032841018/self.fwhm
bg= 2.7725887222397812376689284858327/(self.fwhm**2)
al= 0.63661977236758134307553505349006/self.fwhm
bl= 4.0/self.fwhm**2)
bl= 4.0/self.fwhm**2
gauss = ag* exp(-bg*x2)
lor = al/(1.0+bl*x2)
pv_val = self.eta*lor + (1.0 - self.eta)*gauss
......
......@@ -20,6 +20,7 @@ add_library(${LIBRARY_NAME} SHARED
Extensions/API_IO_Formats.f90
Extensions/API_Reflections_Utilities.f90
Extensions/API_Structure_Factors.f90
Extensions/CFML_API_Calc_Powder_Pattern.f90
Extensions/forpy_mod.f90)
# Forpy must be preprocessed
......
......@@ -21,10 +21,11 @@ module API_Diffraction_Patterns
use CFML_Diffraction_Patterns, only: &
Diffraction_Pattern_Type
use CFML_API_Calc_Powder_Patterm, only: &
Calc_Powder_Pattern
use CFML_API_Calc_Powder_Pattern, only: &
Calc_Powder_Pattern, &
Powder_Pattern_Simulation_Conditions_type
use API_reflection_utilities, only: &
use API_reflections_utilities, only: &
Reflection_List_type_p, &
get_reflection_list_from_arg
......@@ -34,8 +35,76 @@ module API_Diffraction_Patterns
type(Diffraction_Pattern_type), pointer :: p
end type Diffraction_Pattern_type_p
type Powder_Pattern_Simulation_Conditions_type_p
type(Powder_Pattern_Simulation_Conditions_type), pointer :: p
end type Powder_Pattern_Simulation_Conditions_type_p
contains
subroutine get_powpat_from_args(args, pow_pat_sim_c_p, indx)
type(tuple) :: args
type(Powder_Pattern_Simulation_Conditions_type_p), intent(out) :: pow_pat_sim_c_p
integer, optional :: indx
type(object) :: arg_obj
type(dict) :: arg_dict
integer :: ierror
integer :: ii
type(object) :: t
if (present(indx)) then
ierror = args%getitem(arg_obj, indx)
else
ierror = args%getitem(arg_obj, 0)
endif
ierror = cast(arg_dict, arg_obj)
ierror = arg_dict%getitem(pow_pat_sim_c_p%p%job , "job")
ierror = arg_dict%getitem(pow_pat_sim_c_p%p%Lambda, "lambda")
ierror = arg_dict%getitem(pow_pat_sim_c_p%p%U , "u_resolution")
ierror = arg_dict%getitem(pow_pat_sim_c_p%p%V , "v_resolution")
ierror = arg_dict%getitem(pow_pat_sim_c_p%p%W , "w_resolution")
ierror = arg_dict%getitem(pow_pat_sim_c_p%p%X , "x_resolution")
ierror = arg_dict%getitem(pow_pat_sim_c_p%p%Ls , "lorentzian_size")
ierror = arg_dict%getitem(pow_pat_sim_c_p%p%Thmin , "theta_min")
ierror = arg_dict%getitem(pow_pat_sim_c_p%p%Thmax , "theta_max")
ierror = arg_dict%getitem(pow_pat_sim_c_p%p%step , "theta_step")
ierror = arg_dict%getitem(pow_pat_sim_c_p%p%bkg , "background")
end subroutine get_powpat_from_args
subroutine get_diffraction_pattern_type_from_arg(args, diffraction_pattern_type_pointer, indx)
type(tuple) :: args
type(Diffraction_Pattern_type_p), intent(out) :: diffraction_pattern_type_pointer
integer, optional :: indx
type(object) :: arg_obj
type(list) :: arg_list
integer :: diffraction_pattern_type_p12(12)
integer :: ierror
integer :: ii
type(object) :: t
if (present(indx)) then
ierror = args%getitem(arg_obj, indx)
else
ierror = args%getitem(arg_obj, 0)
endif
ierror = cast(arg_list, arg_obj)
do ii=1,12
ierror = arg_list%getitem(t, ii-1)
ierror = cast(diffraction_pattern_type_p12(ii), t)
enddo
diffraction_pattern_type_pointer = transfer(diffraction_pattern_type_p12, diffraction_pattern_type_pointer)
end subroutine get_diffraction_pattern_type_from_arg
function diffraction_patterns_compute_powder_pattern(self_ptr, args_ptr) result(r) bind(c)
type(c_ptr), value :: self_ptr
......@@ -53,6 +122,7 @@ contains
type(Reflection_List_type_p) :: hkl_p
type(Diffraction_Pattern_type_p) :: pattern_p
type(Powder_Pattern_Simulation_Conditions_type_p) :: powpat_p
integer :: pattern_p12(12)
r = C_NULL_PTR ! in case of an exception return C_NULL_PTR
......@@ -61,32 +131,102 @@ contains
! Check if the arguments are OK
ierror = args%len(num_args)
if (num_args /= 2) then
call raise_exception(TypeError, "Calc_Powder_Patter expects exactly 2 arguments")
!Powpat_conditions, Reflectionlist
if (num_args /= 3) then
call raise_exception(TypeError, "Calc_Powder_Pattern expects exactly 3 arguments")
!Powpat_conditions, Reflectionlist, Scalef
call args%destroy
return
endif
!call get_powpat_from_args(args, powpat_p, 0)
allocate(powpat_p%p)
call get_powpat_from_args(args, powpat_p, 0)
call get_reflection_list_from_arg(args, hkl_p, 1)
ierror = args%getitem(arg_obj, 2)
ierror = cast_nonstrict(powpat_p%p%scalef, arg_obj)
allocate(pattern_p%p)
call Calc_Powder_Pattern(powpat_p%p,hkl_p%p,pattern_p%p)
deallocate(powpat_p%p)
pattern_p12 = transfer(pattern_p,pattern_p12)
ierror = list_create(index_obj)
do ii=1,12
ierror = index_obj%append(pattern_p12(ii))
end do
ierror = dict_create(retval)
ierror = retval%setitem("address", index_obj)
r = retval%get_c_ptr()
end function diffraction_patterns_compute_powder_pattern
function diffraction_patterns_get_x(self_ptr, args_ptr) result(r) bind(c)
type(c_ptr), value :: self_ptr
type(c_ptr), value :: args_ptr
type(c_ptr) :: r
type(tuple) :: args
type(dict) :: retval
integer :: num_args
integer :: ierror
type(Diffraction_Pattern_type_p) :: diffraction_pattern_type_pointer
type(ndarray) :: x_array
r = C_NULL_PTR ! in case of an exception return C_NULL_PTR
! use unsafe_cast_from_c_ptr to cast from c_ptr to tuple
call unsafe_cast_from_c_ptr(args, args_ptr)
! Check if the arguments are OK
ierror = args%len(num_args)
! we should also check ierror, but this example does not do complete error checking for simplicity
if (num_args /= 1) then
call raise_exception(TypeError, "diffraction_patterns_get_x expects exactly 1 argument")
call args%destroy
return
endif
! Doing boring stuff
call get_diffraction_pattern_type_from_arg(args, diffraction_pattern_type_pointer)
ierror = ndarray_create(x_array, diffraction_pattern_type_pointer%p%x)
ierror = dict_create(retval)
ierror = retval%setitem("x", x_array)
r = retval%get_c_ptr()
end function diffraction_patterns_get_x
function diffraction_patterns_get_y(self_ptr, args_ptr) result(r) bind(c)
type(c_ptr), value :: self_ptr
type(c_ptr), value :: args_ptr
type(c_ptr) :: r
type(tuple) :: args
type(dict) :: retval
integer :: num_args
integer :: ierror
type(Diffraction_Pattern_type_p) :: diffraction_pattern_type_pointer
type(ndarray) :: y_array
r = C_NULL_PTR ! in case of an exception return C_NULL_PTR
! use unsafe_cast_from_c_ptr to cast from c_ptr to tuple
call unsafe_cast_from_c_ptr(args, args_ptr)
! Check if the arguments are OK
ierror = args%len(num_args)
! we should also check ierror, but this example does not do complete error checking for simplicity
if (num_args /= 1) then
call raise_exception(TypeError, "diffraction_patterns_get_y expects exactly 1 argument")
call args%destroy
return
endif
! Doing boring stuff
call get_diffraction_pattern_type_from_arg(args, diffraction_pattern_type_pointer)
ierror = ndarray_create(y_array, diffraction_pattern_type_pointer%p%y)
ierror = dict_create(retval)
ierror = retval%setitem("y", y_array)
r = retval%get_c_ptr()
end function diffraction_patterns_get_y
end module API_Diffraction_Patterns
......@@ -62,7 +62,9 @@ module API_init
reflections_utilities_hkl_uni_reflist
use API_Diffraction_Patterns, only: &
diffraction_patterns_compute_powder_pattern
diffraction_patterns_compute_powder_pattern, &
diffraction_patterns_get_x, &
diffraction_patterns_get_y
implicit none
......@@ -90,15 +92,25 @@ CONTAINS
integer :: ierror
ierror = forpy_initialize()
call method_table%init(33)
call method_table%init(35)
!--------------------------
! Diffraction Patterns (1)
! Diffraction Patterns (3)
!--------------------------
call method_table%add_method("diffraction_patterns_compute_powder_pattern", & ! method name
"compute the powder diffraction pattern from some experimental conditions and a set of reflections", & !doc-string
METH_VARARGS, & ! this method takes arguments but no keyword arguments
c_funloc(diffraction_patterns_compute_powder_pattern)) ! address of Fortran function to add
call method_table%add_method("diffraction_patterns_get_x", & ! method name
"Get x array", & !doc-string
METH_VARARGS, & ! this method takes arguments but no keyword arguments
c_funloc(diffraction_patterns_get_x)) ! address of Fortran function to add
call method_table%add_method("diffraction_patterns_get_y", & ! method name
"Get y array", & !doc-string
METH_VARARGS, & ! this method takes arguments but no keyword arguments
c_funloc(diffraction_patterns_get_y)) ! address of Fortran function to add
!--------------------------
! Crystallographic Symmetry (7)
!--------------------------
......
module CFML_API_PowderSimulation
module CFML_API_Calc_Powder_Pattern
use CFML_GlobalDeps, only: to_Deg
use CFML_Math_General, only: asind,locate
......@@ -10,13 +10,13 @@ use CFML_PowderProfiles_CW, only: TCH_pVoigt,PseudoVoigt
implicit none
Type :: PowPat_CW_Conditions
Type :: Powder_Pattern_Simulation_Conditions_type
character(len=140) :: title
integer :: job
real :: Lambda, U, V, W, X, Ls, Gs
real :: Thmin, Thmax, step
real :: scalef,bkg
End Type PowPat_CW_Conditions
End Type Powder_Pattern_Simulation_Conditions_type
CONTAINS
......@@ -40,7 +40,7 @@ subroutine compute_powder_pattern(input_filename, mode, PowPat_Conditions, Pat)
! Input/Output
character(len=*) :: input_filename
character(len=3) :: mode
Type(PowPat_CW_Conditions) :: PowPat_Conditions
Type(Powder_Pattern_Simulation_Conditions_type) :: PowPat_Conditions
Type(Diffraction_Pattern_Type) :: Pat
!
......@@ -123,7 +123,7 @@ subroutine Read_Power_Pattern_Simulation_conditions(fich_cfl,PowPat_Conditions)
!
implicit none
Type(file_list_type) :: fich_cfl
Type(PowPat_CW_Conditions) :: PowPat_Conditions
Type(Powder_Pattern_Simulation_Conditions_type) :: PowPat_Conditions
character(len=8) :: radiation
integer :: ier
......@@ -218,7 +218,7 @@ End Subroutine TCH
Subroutine Calc_Powder_Pattern(PowPat_Conditions,Hkl,Pat)
!---- Argument ----!
Type(PowPat_CW_Conditions), intent(in) :: PowPat_Conditions
Type(Powder_Pattern_Simulation_Conditions_type), intent(in) :: PowPat_Conditions
Type(Reflection_List_Type), intent(in) :: hkl
Type(Diffraction_Pattern_Type), intent(out) :: Pat
......@@ -288,4 +288,4 @@ Subroutine Calc_Powder_Pattern(PowPat_Conditions,Hkl,Pat)
return
End Subroutine Calc_Powder_Pattern
end module
end module CFML_API_Calc_Powder_Pattern
......@@ -32,6 +32,21 @@ class PowderPatternSimulationConditions():
self.job = PowderPatternSimulationSource.Neutrons
self.lorentzian_size = 1900.0
self.bkg = 50.0
def to_dict(self):
return {
"lambda":self.lamb,
"u_resolution":self.u_resolution,
"v_resolution":self.v_resolution,
"w_resolution":self.w_resolution,
"x_resolution":self.x_resolution,
"theta_min":self.theta_min,
"theta_step":self.theta_step,
"theta_max":self.theta_max,
"job":self.job,
"lorentzian_size":self.lorentzian_size,
"background":self.bkg
}
def getSinThetaOverLambdaMax(self):
angle_max = min((self.theta_max+10.0)*0.5,90.0) * 180 / np.pi
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment