Plugin no longer pretends to describe an anisotropic material

Author: tkittel

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]

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<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 ) {
@@ -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.;
@@ -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
 }

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<HKLPlane> m_hklPlanes;
-    NCrystal::RotMatrix m_lab2cry;
-    NCrystal::RotMatrix m_reclat;
-
   };
 
 }

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<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