.. _program_listing_file_src_pdfs_PDFSegmentSchechter.hpp:

Program Listing for File PDFSegmentSchechter.hpp
================================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_pdfs_PDFSegmentSchechter.hpp>` (``src/pdfs/PDFSegmentSchechter.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   #ifndef PDFSEGMENTSCHECHTER_HPP
   #define PDFSEGMENTSCHECHTER_HPP
   
   #include "../utils/RandomDist.hpp"
   #include "../utils/RngThread.hpp"
   #include "PDFCommons.hpp"
   #include "PDFSegment.hpp"
   #include <algorithm>
   #include <cmath>
   #include <fstream>
   #include <gsl/gsl_sf_gamma.h>
   #include <toml.hpp>
   
   namespace pdfs {
   
       class PDFSegmentSchechter : public PDFSegment {
       public:
   
           PDFSegmentSchechter(double sMin, double sMax, double sStar, double alpha) :
               PDFSegment(sMin, sMax), 
               sStar_(sStar), 
               alpha_(alpha),
               norm_(1.0 / (
                       std::pow(sStar_, alpha_ + 1) * 
                       (gsl_sf_gamma_inc(alpha_ + 1, sMin_ / sStar_) -
                       gsl_sf_gamma_inc(alpha_ + 1, sMax_ / sStar_)
                       ))
                   )
               { }
           PDFSegmentSchechter(std::ifstream& file,
               FileFormats fmt,
               double &sMin,
               double &sMax,
               double &wgt);
           PDFSegmentSchechter(toml::node_view<const toml::node> node,
               FileFormats fmt,
               double &sMin,
               double &sMax,
               double &wgt);
           PDFSegmentSchechter(const PDFSegmentSchechter&) = default;
           auto operator=(const PDFSegmentSchechter&) -> PDFSegmentSchechter& = default;
           PDFSegmentSchechter(PDFSegmentSchechter&&) = default;
           auto operator=(PDFSegmentSchechter&&) -> PDFSegmentSchechter& = default;
           ~PDFSegmentSchechter() override = default;
   
           // Evaluation functions
           [[nodiscard]] auto operator()(double x) const -> double override {
               if (x < sMin_ || x > sMax_) {
                   return 0.0;
               }
               return norm_ * std::pow(x, alpha_) * std::exp(-x / sStar_);
           }
           [[nodiscard]] auto expectationValue(const double a, const double b) const -> double override {
               if (a >= b || a > sMax_ || b < sMin_)
               {
                   return 0.0; // Invalid range for expectation value calculation
               }
               if (a == sMax_)
               {
                   return sMax_; // Handle edge cases
               }
               if (b == sMin_)
               {
                   return sMin_; // Handle edge cases
               }
               const double aClamped = std::max(a, sMin_);
               const double bClamped = std::min(b, sMax_);
               return sStar_ * (
                   gsl_sf_gamma_inc(alpha_ + 2, aClamped / sStar_) -
                   gsl_sf_gamma_inc(alpha_ + 2, bClamped / sStar_)
               ) / (
                   gsl_sf_gamma_inc(alpha_ + 1, aClamped / sStar_) -
                   gsl_sf_gamma_inc(alpha_ + 1, bClamped / sStar_)
               );
           }
           [[nodiscard]] auto expectationValue() const -> double override {
               return expectationValue(sMin_, sMax_); // Expectation value over the full range of the segment
           }
           [[nodiscard]] auto integral(const double a, const double b) const -> double override {
               if (a >= b || a >= sMax_ || b <= sMin_) {
                   return 0.0; // Invalid range for integral calculation
               }
               const double aClamped = std::max(a, sMin_);
               const double bClamped = std::min(b, sMax_);
               return norm_ * std::pow(sStar_, alpha_ + 1) * (
                   gsl_sf_gamma_inc(alpha_ + 1, aClamped / sStar_) -
                   gsl_sf_gamma_inc(alpha_ + 1, bClamped / sStar_)
               );
           }
   
           // Drawing function
           [[nodiscard]] auto draw(const double a, const double b) const -> double override {
               const double aClamped = std::max(a, sMin_);
               const double bClamped = std::min(b, sMax_);
               if (aClamped >= bClamped) {
                   return 0.0; // Invalid range for drawing
               }
               const double u = utils::canonical(utils::rng()()); // Uniform random number in [0, 1)
               // Find the value of the deviate y by numerically solving
               // u = \int_a^y x^alpha exp(-x / sStar) dx /
               //     \int_a^b x^alpha exp(-x / sStar) dx
               // using a simple bisection search.
               const double gammaA = gsl_sf_gamma_inc(alpha_ + 1, aClamped / sStar_);
               const double gammaB = gsl_sf_gamma_inc(alpha_ + 1, bClamped / sStar_);
               const double denom = gammaA - gammaB;
               double y = aClamped;
               double step = (bClamped - aClamped) / 2.0;
               while (true) {
                   y += step;
                   const double gammaY = gsl_sf_gamma_inc(alpha_ + 1, y / sStar_);
                   const double resid = u - ((gammaA - gammaY) / denom);
                   if (fabs(resid) < 1e-8) {
                       break;  // Converged
                   }
                   if (resid < 0) {
                       y -= step;
                       step /= 2.0;
                       if (step < 1e-8) {
                           break;
                       }
                   }
               }
               return y;
           }
           [[nodiscard]] auto draw() const -> double override {
               return draw(sMin_, sMax_);
           }
   
           // Serialization
           void toToml(toml::table& tbl) const override;
   
       private:
           double sStar_;  
           double alpha_;  
           double norm_;   
       };
   
   } // namespace pdfs
   
   #endif // PDFSEGMENTSCHECHTER_HPP