diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 02831fe..83f8d12 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -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' diff --git a/pyproject.toml b/pyproject.toml index afafef4..172a401 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -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] diff --git a/src/NCCrystallineTexture.cc b/src/NCCrystallineTexture.cc index 156bfc3..ada1e45 100644 --- a/src/NCCrystallineTexture.cc +++ b/src/NCCrystallineTexture.cc @@ -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 @@ -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 ) @@ -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 ); @@ -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; @@ -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() ).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 ) { @@ -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() }; - - 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.; @@ -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::randDirectionGivenScatterMu( rng, mu, ndirlab.as() ); - 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 } diff --git a/src/NCCrystallineTexture.hh b/src/NCCrystallineTexture.hh index f8d3046..80a416d 100644 --- a/src/NCCrystallineTexture.hh +++ b/src/NCCrystallineTexture.hh @@ -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" @@ -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. @@ -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: @@ -66,9 +64,6 @@ namespace NCPluginNamespace { double strength; //dspacing*fsq*xsectfact }; std::vector m_hklPlanes; - NCrystal::RotMatrix m_lab2cry; - NCrystal::RotMatrix m_reclat; - }; } diff --git a/src/NCPluginFactory.cc b/src/NCPluginFactory.cc index 4236e85..301c746 100644 --- a/src/NCPluginFactory.cc +++ b/src/NCPluginFactory.cc @@ -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 @@ -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: @@ -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( CrystallineTexture::createFromInfo( sco, - cfg.info(), + = NC::makeSO( CrystallineTexture::createFromInfo( cfg.info(), sc_pp.get() ) ); //Now we just need to combine this with all the other physics. So ask the