ProcessSequence.h

#ifndef _include_ProcessSequence_h_
#define _include_ProcessSequence_h_

#include <corsika/process/ProcessReturn.h>

#include <cmath>
#include <iostream>
#include <typeinfo>

namespace corsika::process {  
  
  /**
     \class BaseProcess

     The structural base type of a process object in a
     ProcessSequence. Both, the ProcessSequence and all its elements
     are of type BaseProcess<T>

   */

  template <typename derived>
  struct BaseProcess {
    derived& GetRef() { return static_cast<derived&>(*this); }
    const derived& GetRef() const { return static_cast<const derived&>(*this); }
  };

  /**
     \class ProcessSequence

     A compile time static list of processes. The compiler will
     generate a new type based on template logic containing all the
     elements.

     \comment Using CRTP pattern,
     https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
   */

  template <typename T1, typename T2>
  class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> > {
  public:
    const T1& A;
    const T2& B;

    ProcessSequence(const T1& in_A, const T2& in_B)
        : A(in_A)
        , B(in_B) {}

    template <typename Particle, typename Trajectory, typename Stack>
    inline EProcessReturn DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
      EProcessReturn ret = EProcessReturn::eOk;
      /*ret |=*/ A.DoContinuous(p, t, s);
      /*ret |=*/ B.DoContinuous(p, t, s);
      return ret;
    } // add trajectory

    template <typename D>
    inline double MinStepLength(D& d) const {
      return std::min(A.MinStepLength(d), B.MinStepLength(d));
    }

    template <typename Particle, typename Trajectory>
    inline Trajectory Transport(Particle& p, double& length) const {
      A.Transport(p, length); // todo: maybe check (?) if there is more than one Transport process implemented?? 
      return B.Transport(p, length); // need to do this also to decide which Trajectory to return!!!!
    }

    template <typename Particle, typename Stack>
    void DoDiscrete(Particle& p, Stack& s) const {
      A.DoDiscrete(p, s);
      B.DoDiscrete(p, s);
    }

    /// TODO the const_cast is not nice, think about the constness here
    inline void Init() const {
      const_cast<T1*>(&A)->Init();
      const_cast<T2*>(&B)->Init();
    }
  };

  /// the + operator that assembles more BaseProcess objects into a ProcessSequence
  template <typename T1, typename T2>
  inline const ProcessSequence<T1, T2> operator+(const BaseProcess<T1>& A,
                                                 const BaseProcess<T2>& B) {
    return ProcessSequence<T1, T2>(A.GetRef(), B.GetRef());
  }

  /*
    template <typename T1>
    struct depth_lhs
    {
    static const int num = 0;
    };



    // terminating condition
    template <typename T1, typename T2>
    struct depth_lhs< Sequence<T1,T2> >
    {
    // try to expand the left node (T1) which might be a Sequence type
    static const int num = 1 + depth_lhs<T1>::num;
    };
  */

  /*
    template <typename T1>
    struct mat_ptrs
    {
    static const int num = 0;

    inline static void
    get_ptrs(const Process** ptrs, const T1& X)
    {
    ptrs[0] = reinterpret_cast<const Process*>(&X);
    }
    };


    template <typename T1, typename T2>
    struct mat_ptrs< Sequence<T1,T2> >
    {
    static const int num = 1 + mat_ptrs<T1>::num;

    inline static void
    get_ptrs(const Process** in_ptrs, const Sequence<T1,T2>& X)
    {
    // traverse the left node
    mat_ptrs<T1>::get_ptrs(in_ptrs, X.A);
    // get address of the matrix on the right node
    in_ptrs[num] = reinterpret_cast<const Process*>(&X.B);
    }
    };
  */

  /*
    template<typename T1, typename T2>
    const Process&
    Process::operator=(const Sequence<T1,T2>& X)
    {
    int N = 1 + depth_lhs< Sequence<T1,T2> >::num;
    const Process* ptrs[N];
    mat_ptrs< Sequence<T1,T2> >::get_ptrs(ptrs, X);
    int r = ptrs[0]->rows;
    int c = ptrs[0]->cols;
    // ... check that all matrices have the same size ...
    set_size(r, c);
    for(int j=0; j<r*c; ++j)
    {
    double sum = ptrs[0]->data[j];
    for(int i=1; i<N; ++i)
    {
    sum += ptrs[i]->data[j];
    }
    data[j] = sum;
    }
    return *this;
    }
  */

} // namespace corsika::process

#endif