Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 4 additions & 23 deletions .github/workflows/test.yml
Original file line number Diff line number Diff line change
Expand Up @@ -22,32 +22,13 @@ jobs:
path: src

- name: Install dependencies
run: pip install 'ncrystal-core>=3.9.86' "scikit-build-core>=0.10" "ncrystal-pypluginmgr>=0.0.3" "ncrystal-python>=3.9.86"
run: pip install 'ncrystal>=4.1.4'

- name: Install plugin
run: pip install ./src

- name: Verify plugin loading
shell: python
run: |
import os
os.environ['NCRYSTAL_PLUGIN_RUNTESTS'] = '1'
os.environ['NCRYSTAL_REQUIRED_PLUGINS'] = 'CrysText'
import NCrystal
NCrystal.browsePlugins(dump=True)

- name: Use plugin file
run: nctool -d 'plugins::CrysText/Al_sg225.ncmat'
run: ncrystal-pluginmanager --test CrysText

- name: Load all plugin files
shell: python
run: |
from NCrystal.datasrc import browseFiles
from NCrystal.core import createScatter
for f in browseFiles(factory='plugins'):
if f.name.startswith('CrysText/'):
print('Loading f.fullKey')
createScatter( f'{f.fullKey}'
";dir1=@crys_hkl:0,1,0@lab:0,1,0"
";dir2=@crys_hkl:1,0,0@lab:1,0,0"
";mos=0.1deg;dirtol=180deg" )
- name: Use explicit plugin file
run: nctool -d 'plugins::CrysText/Al_sg225.ncmat'
4 changes: 2 additions & 2 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -23,14 +23,14 @@
name = "ncrystal-plugin-CrysText"
version = "0.0.1"
license = {file = "LICENSE"}
dependencies = [ "ncrystal-core>=3.9.85", "ncrystal-pypluginmgr>=0.0.3" ]
dependencies = [ "ncrystal-core>=4.1.4", "ncrystal-pypluginmgr>=0.0.3" ]
classifiers = [
"Programming Language :: Python :: 3",
"License :: OSI Approved :: Apache Software License",
]

[build-system]
requires = ["scikit-build-core>=0.10", "ncrystal-core>=3.9.85" ]
requires = ["scikit-build-core>=0.10", "ncrystal-core>=4.1.4" ]
build-backend = "scikit_build_core.build"

[tool.scikit-build]
Expand Down
68 changes: 9 additions & 59 deletions src/NCCrystallineTexture.cc
Original file line number Diff line number Diff line change
Expand Up @@ -63,8 +63,7 @@ bool NCP::CrystallineTexture::isApplicable( const NC::Info& info )
return info.countCustomSections(pluginNameUpperCase()) > 0;
}

NCP::CrystallineTexture NCP::CrystallineTexture::createFromInfo( const NC::SCOrientation& sco,
const NC::Info& info,
NCP::CrystallineTexture NCP::CrystallineTexture::createFromInfo( const NC::Info& info,
NC::PlaneProvider * std_pp)
{
//Parse the content of our custom section. In case of syntax errors, we should
Expand Down Expand Up @@ -114,11 +113,10 @@ NCP::CrystallineTexture NCP::CrystallineTexture::createFromInfo( const NC::SCOri
const NCrystal::StructureInfo& struct_info = info.getStructureInfo();

//Parsing done! Create and return our model:
return CrystallineTexture(sco,preferred_orientation1,R1,f1,preferred_orientation2,R2,f2,struct_info,std_pp);
return CrystallineTexture(preferred_orientation1,R1,f1,preferred_orientation2,R2,f2,struct_info,std_pp);
}

NCP::CrystallineTexture::CrystallineTexture( const NC::SCOrientation& sco,
const NCrystal::Vector& preferred_orientation1, double R1, double f1,
NCP::CrystallineTexture::CrystallineTexture( const NCrystal::Vector& preferred_orientation1, double R1, double f1,
const NCrystal::Vector& preferred_orientation2, double R2, double f2,
const NCrystal::StructureInfo& struct_info,
NCrystal::PlaneProvider * plane_provider )
Expand All @@ -129,13 +127,6 @@ NCP::CrystallineTexture::CrystallineTexture( const NC::SCOrientation& sco,
m_R2(R2),
m_f2(f2)
{
//Important note to developers who are using the infrastructure in the
//testcode/ subdirectory: If you change the number or types of the arguments
//for the constructor here, you should make sure to perform a corresponding
//change in three files in the testcode/ directory: _cbindings.py,
//__init__.py, and NCForPython.cc - that way you can still instantiate your
//model directly from your python test code).

nc_assert( preferred_orientation1.mag() > 0.0 );
nc_assert( m_R1 > 0.0 );
nc_assert( m_f1 > 0.0 );
Expand All @@ -144,12 +135,10 @@ NCP::CrystallineTexture::CrystallineTexture( const NC::SCOrientation& sco,
nc_assert( m_f2 > 0.0 );
nc_assert( plane_provider->canProvide() );

m_reclat = getReciprocalLatticeRot( struct_info );
NCrystal::RotMatrix lattice_rot = NC::getLatticeRot( struct_info.lattice_a, struct_info.lattice_b, struct_info.lattice_c,
struct_info.alpha*NC::kDeg, struct_info.beta*NC::kDeg, struct_info.gamma*NC::kDeg );
m_lab2cry = getCrystal2LabRot( sco, m_reclat ).getInv();

//RotMatrix cry2lab = getCrystal2LabRot( sco, m_reclat );

double V0numAtom = struct_info.n_atoms * struct_info.volume;
const double xsectfact = 0.5 / V0numAtom;

Expand All @@ -171,15 +160,9 @@ NCP::CrystallineTexture::CrystallineTexture( const NC::SCOrientation& sco,
}
}

double NCP::CrystallineTexture::calcCrossSection( NC::NeutronEnergy neutron_ekin, const NC::NeutronDirection& ndirlab ) const
double NCP::CrystallineTexture::calcCrossSection( NC::NeutronEnergy neutron_ekin ) const
{
(void)ndirlab;//FIXME: Actually use this!!
//auto ndir = ( m_lab2cry * ndirlab.as<NC::Vector>() ).unit();
//ndir[0],ndir[1],ndir[2]
//auto neutron_HKL = m_reclat * ndir;

double xs_in_barns = 0.0;

const double wl = neutron_ekin.wavelength().dbl();
const double wlsq = NC::ncsquare(wl);
for ( auto& e: m_hklPlanes ) {
Expand All @@ -190,20 +173,14 @@ double NCP::CrystallineTexture::calcCrossSection( NC::NeutronEnergy neutron_ekin
xs_in_barns += e.strength * (P1 * m_f1 + P2 * m_f2);
}
xs_in_barns *= 2.*wlsq; //consideration of the negative hkl

return xs_in_barns;
}

NC::ScatterOutcome NCP::CrystallineTexture::sampleScatteringEvent( NC::RNG& rng, NC::NeutronEnergy neutron_ekin, const NC::NeutronDirection& ndirlab ) const
NC::ScatterOutcomeIsotropic NCP::CrystallineTexture::sampleScatteringEvent( NC::RNG& rng, NC::NeutronEnergy neutron_ekin ) const
{
//Don't do anything:
//return { neutron_ekin, ndirlab };
//return { neutron_ekin, NC::randIsotropicDirection(rng).as<NeutronDirection>() };

NC::NeutronDirection outndirlab = ndirlab; //outgoing neutron direction
const double wl = neutron_ekin.wavelength().dbl();
const double wlsq = NC::ncsquare(wl);
const double xs = calcCrossSection( neutron_ekin, ndirlab ) / (2.*wlsq); //calculate xs
const double xs = calcCrossSection( neutron_ekin ) / (2.*wlsq);
const double rnd = rng.generate(); //random number on [0;1]

double left_bound = 0.;
Expand All @@ -220,38 +197,11 @@ NC::ScatterOutcome NCP::CrystallineTexture::sampleScatteringEvent( NC::RNG& rng,
const double E_hkl = 0.5 * NC::kPiSq * NC::const_hhm / NC::ncsquare(e.d_hkl);
const double mu = 1. - 2 * E_hkl / neutron_ekin.dbl();
nc_assert( NC::ncabs(mu) <= 1.0 );
outndirlab.as<NC::Vector>() = NC::randDirectionGivenScatterMu( rng, mu, ndirlab.as<NC::Vector>() );
break;
return { neutron_ekin, NC::CosineScatAngle{ mu } };
}
else {
left_bound = right_bound;
}
}

return { neutron_ekin, outndirlab };

//if ( ! (neutron_ekin > m_cutoffekin) ) {
//Special case: We are asked to sample a scattering event for a neutron
//energy where we have zero cross section! Although in a real simulation we
//would usually not expect this to happen, users with custom code might
//still generate such calls. The only consistent thing to do when the cross
//section is zero is to not change the neutron state parameters, which means:
//result.ekin_final = neutron_ekin;
//result.mu = 1.0;
//return result;
//}

//Implement our actual model here. Of course it is trivial for the example
//model. For a more realistic or complicated model, it might be that
//additional helper classes or functions should be created and used, in order
//to keep the code here manageable:

//result.ekin_final = neutron_ekin;//Elastic
//result.mu = randIsotropicScatterMu(rng).dbl();

//Same as coherent elastic scattering
//result.ekin_final = neutron_ekin.dbl();
//result.mu = randIsotropicScatterMu(rng).dbl(); // Take isotropic first for test

// return result;
return { neutron_ekin, NC::CosineScatAngle{1.0} };//state unchanged
}
17 changes: 6 additions & 11 deletions src/NCCrystallineTexture.hh
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,6 @@
#include "NCrystal/NCPluginBoilerplate.hh"//Common stuff (includes NCrystal
//public API headers, sets up
//namespaces and aliases)
#include "NCrystal/interfaces/NCSCOrientation.hh"
#include "NCrystal/internal/utils/NCRotMatrix.hh"
#include "NCrystal/internal/utils/NCVector.hh"
#include "NCrystal/internal/extd_utils/NCPlaneProvider.hh"

Expand All @@ -27,9 +25,10 @@ namespace NCPluginNamespace {
//case of syntax errors in the @CUSTOM_ section data):

static bool isApplicable( const NC::Info& );
static CrystallineTexture createFromInfo( const NC::SCOrientation&, const NC::Info&, NC::PlaneProvider * = nullptr );//will raise BadInput in case of syntax errors
static CrystallineTexture createFromInfo( const NC::Info&,
NC::PlaneProvider * = nullptr );//will raise BadInput in case of syntax errors

//The crystalline texuture or preferred orientation correction is taken into account by introducing
//The crystalline texuture or preferred orientation correction is taken into account by introducing
//the cylindrically symmetric Pole-Density Distribution Function (PDDF) or preferred orientation
//distribution function P_hkl(lambda, theta_hkl) which depends on the orientation angle theta_hkl,
//i.e., angle between the preferred orientation and the plan vectors hkl.
Expand All @@ -38,19 +37,18 @@ namespace NCPluginNamespace {

//Constructor gets constant cross section value, and the neutron wavelength
//cutoff:
CrystallineTexture( const NC::SCOrientation&,
const NCrystal::Vector& preferred_orientation1, double R1, double f1,
CrystallineTexture( const NCrystal::Vector& preferred_orientation1, double R1, double f1,
const NCrystal::Vector& preferred_orientation2, double R2, double f2,
const NCrystal::StructureInfo& struct_info,
NC::PlaneProvider * std_pp = nullptr );

//Provide cross sections for a given neutron:
double calcCrossSection( NC::NeutronEnergy, const NC::NeutronDirection&) const;
double calcCrossSection( NC::NeutronEnergy ) const;

//Sample scattering event (rng is random number stream). Results are given
//as the final ekin of the neutron and scat_mu which is cos(scattering_angle).

NC::ScatterOutcome sampleScatteringEvent( NC::RNG&, NC::NeutronEnergy, const NC::NeutronDirection& ) const;
NC::ScatterOutcomeIsotropic sampleScatteringEvent( NC::RNG&, NC::NeutronEnergy ) const;

private:
//Data members:
Expand All @@ -66,9 +64,6 @@ namespace NCPluginNamespace {
double strength; //dspacing*fsq*xsectfact
};
std::vector<HKLPlane> m_hklPlanes;
NCrystal::RotMatrix m_lab2cry;
NCrystal::RotMatrix m_reclat;

};

}
Expand Down
20 changes: 8 additions & 12 deletions src/NCPluginFactory.cc
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ namespace NCPluginNamespace {

using PhysicsModel = CrystallineTexture;

class PluginScatter final : public NC::ProcImpl::ScatterAnisotropicMat {
class PluginScatter final : public NC::ProcImpl::ScatterIsotropicMat {
public:

//The factory wraps our custom PhysicsModel helper class in an NCrystal API
Expand All @@ -40,20 +40,18 @@ namespace NCPluginNamespace {
PluginScatter( PhysicsModel && pm ) : m_pm(std::move(pm)) {}


NC::CrossSect crossSection(NC::CachePtr&,
NC::NeutronEnergy neutron_ekin,
const NC::NeutronDirection& ndirlab ) const override
NC::CrossSect crossSectionIsotropic(NC::CachePtr&,
NC::NeutronEnergy neutron_ekin) const override
{
return NC::CrossSect{ m_pm.calcCrossSection(neutron_ekin, ndirlab) };
return NC::CrossSect{ m_pm.calcCrossSection(neutron_ekin) };
}


NC::ScatterOutcome sampleScatter(NC::CachePtr&,
NC::ScatterOutcomeIsotropic sampleScatterIsotropic(NC::CachePtr&,
NC::RNG& rng,
NC::NeutronEnergy neutron_ekin,
const NC::NeutronDirection& ndirlab) const override
NC::NeutronEnergy neutron_ekin) const override
{
return m_pm.sampleScatteringEvent( rng, neutron_ekin, ndirlab );
return m_pm.sampleScatteringEvent( rng, neutron_ekin );
}

private:
Expand Down Expand Up @@ -100,10 +98,8 @@ NCP::PluginFactory::produce( const NC::FactImpl::ScatterRequest& cfg ) const
//Ok, we are selected as the provider! First create our own scatter model:

auto sc_pp = createStdPlaneProvider( cfg.infoPtr() );
auto sco = cfg.createSCOrientation();
auto sc_ourmodel
= NC::makeSO<PluginScatter>( CrystallineTexture::createFromInfo( sco,
cfg.info(),
= NC::makeSO<PluginScatter>( CrystallineTexture::createFromInfo( cfg.info(),
sc_pp.get() ) );

//Now we just need to combine this with all the other physics. So ask the
Expand Down