Merge branch 'main' into null_suppression_optim

eisvogel/core/impl/GreensFunction.hh

@@ -71,9 +71,11 @@ public:
 
   // Apply this Green's function to the line current segment `curr_seg` and accumulate the result in `signal`.
   // The signal is calculated starting from time `t_sig_start` with `num_samples` samples using `t_sig_samp` as sampling interval.
-  template <class KernelT, class QuadratureT = Quadrature::TrapezoidalRule>
+  template <class KernelT, typename ResultT = float, class QuadratureT = Quadrature::TrapezoidalRule>
   void apply_accumulate(const LineCurrentSegment& seg, scalar_t t_sig_start, scalar_t t_sig_samp, std::size_t num_samples,
-			std::vector<scalar_t>& signal, Green::OutOfBoundsBehavior oob_mode = Green::OutOfBoundsBehavior::RaiseError);
+			std::vector<ResultT>& signal, 
+			Green::OutOfBoundsBehavior oob_mode = Green::OutOfBoundsBehavior::RaiseError,
+			scalar_t weight = 1.0);
 
   template <class KernelT>
   void fill_array(const RZTCoordVector& start_coords, const RZTCoordVector& end_coords, const RZTVector<std::size_t>& num_samples, chunk_t& array);
@@ -84,9 +86,9 @@ public:
 public: // TODO: make this private, called by unit test at the moment
 
   // Calculate inner product with source current over the time interval [t_start, t_end) and accumulate into `result`
-  template <class KernelT>
+  template <class KernelT, typename ResultT = float>
   void accumulate_inner_product(const RZCoordVectorView coords, scalar_t t_start, scalar_t t_samp, std::size_t num_samples,
-				const RZFieldVectorView source, std::vector<scalar_t>::iterator result, scalar_t weight = 1.0f,
+				const RZFieldVectorView source, std::vector<ResultT>::iterator result, scalar_t weight = 1.0f,
 				Green::OutOfBoundsBehavior oob_mode = Green::OutOfBoundsBehavior::RaiseError);
 
 private:

eisvogel/core/impl/GreensFunction.hxx

@@ -182,9 +182,9 @@ void CylindricalGreensFunction::fill_array(const RZTCoordVector& start_pos, cons
   }  
 }
 
-template <class KernelT, class QuadratureT>
+template <class KernelT, typename ResultT, class QuadratureT>
 void CylindricalGreensFunction::apply_accumulate(const LineCurrentSegment& seg, scalar_t t_sig_start, scalar_t t_sig_samp, std::size_t num_samples,
-						 std::vector<scalar_t>& signal, Green::OutOfBoundsBehavior oob_mode) {
+						 std::vector<ResultT>& signal, Green::OutOfBoundsBehavior oob_mode, scalar_t weight) {
 
   // std::cout << "HH in apply_accumulate HH" << std::endl;
 
@@ -315,8 +315,9 @@ void CylindricalGreensFunction::apply_accumulate(const LineCurrentSegment& seg,
 	assert(block_sample_ind_start + block_num_samples <= signal.size());
 
 	auto block_result = signal.begin() + block_sample_ind_start;	
-	accumulate_inner_product<KernelT>(coords_rz[i_pt], convolution_t_start, t_sig_samp, block_num_samples, source_rz[i_pt], block_result,
-					  quadrature_weights[i_pt] * itgr_step * (-1.0), oob_mode);  // negative sign from how Green's function is defined
+	accumulate_inner_product<KernelT, ResultT>(coords_rz[i_pt], convolution_t_start, t_sig_samp, block_num_samples, source_rz[i_pt], block_result,
+						   quadrature_weights[i_pt] * itgr_step * weight * (-1.0),  // negative sign from how Green's function is defined
+						   oob_mode);
 
 	// std::cout << " . . . . . . . . " << std::endl;
       }
@@ -326,9 +327,9 @@ void CylindricalGreensFunction::apply_accumulate(const LineCurrentSegment& seg,
   }  
 }
 
-template <class KernelT>
+template <class KernelT, typename ResultT>
 void CylindricalGreensFunction::accumulate_inner_product(const RZCoordVectorView rz_coords, scalar_t t_start, scalar_t t_samp, std::size_t num_samples,
-							 const RZFieldVectorView source, std::vector<scalar_t>::iterator result, scalar_t weight,
+							 const RZFieldVectorView source, std::vector<ResultT>::iterator result, scalar_t weight,
 							 Green::OutOfBoundsBehavior oob_mode) {
 
   // Buffer to hold the interpolated values

examples/plot_greens_function/animate_greens_function.py

@@ -19,12 +19,12 @@ def generate_animation(outpath, data, range_x, range_y, range_t, xlabel = "r [m]
     gs = GridSpec(1, 1, figure = fig)
     ax = fig.add_subplot(gs[0])
 
-    linthresh = 1e-5
+    linthresh = 1e-8
     vmin = -2e-3
     vmax = 2e-3
 
-    fieldplot = ax.pcolorfast(range_x, range_y, np.transpose(data[:, :, 0]), cmap = "bwr",
-                              norm = colors.SymLogNorm(linthresh = linthresh, linscale = 1, base = 10, vmin = vmin, vmax = vmax))
+    fieldplot = ax.imshow(np.flip(np.transpose(data[:, :, 0]), axis = 0), cmap = "bwr", extent = [*range_x, *range_y], interpolation = "bicubic",
+                          norm = colors.SymLogNorm(linthresh = linthresh, linscale = 1, base = 10, vmin = vmin, vmax = vmax))
     ax.set_aspect("equal")
     ax.set_xlabel(xlabel, fontsize = fs)
     ax.set_ylabel(ylabel, fontsize = fs)
@@ -46,7 +46,7 @@ def generate_animation(outpath, data, range_x, range_y, range_t, xlabel = "r [m]
 
     def animate(it):
         ind_t = downsample_t * it
-        fieldplot.set_array(np.transpose(data[:, :, ind_t]))
+        fieldplot.set_array(np.flip(np.transpose(data[:, :, ind_t]), axis = 0))
         timestep_text.set_text(r"$t = {:.2f}\,\mathrm{{ns}}$".format(range_t[0] + ind_t * (range_t[1] - range_t[0]) / num_frames))
         print(f"Rendered frame {ind_t} / {num_frames}")
         return fieldplot

extern/corsika/include/C8SignalCalculator.hh

@@ -14,7 +14,7 @@ public:
   void get_antenna_location(std::array<float, 3>& ant_xyz);
 
   void calculate(const std::array<float, 4>& track_start_xyzt, const std::array<float, 4>& track_end_xyzt, float track_charge,
-		 float t_sig_start, float t_sig_samp, std::size_t num_samples, std::vector<float>& signal);
+		 float t_sig_start, float t_sig_samp, std::size_t num_samples, std::vector<double>& signal, float weight = 1.0);
 
 private:
 

extern/corsika/src/C8SignalCalculator.cxx

@@ -16,12 +16,12 @@ void C8SignalCalculator::get_antenna_location(std::array<float, 3>& ant_xyz) {
 }
 
 void C8SignalCalculator::calculate(const std::array<float, 4>& track_start_xyzt, const std::array<float, 4>& track_end_xyzt, float track_charge,
-				   float t_sig_start, float t_sig_samp, std::size_t num_samples, std::vector<float>& signal) {
+				   float t_sig_start, float t_sig_samp, std::size_t num_samples, std::vector<double>& signal, float weight) {
 
   LineCurrentSegment track(XYZCoordVector{track_start_xyzt[0], track_start_xyzt[1], track_start_xyzt[2]},    // track start position
 			   XYZCoordVector{track_end_xyzt[0], track_end_xyzt[1], track_end_xyzt[2]},   // track end position
 			   track_start_xyzt[3], track_end_xyzt[3], track_charge);
 
   Green::OutOfBoundsBehavior oob_mode = Green::OutOfBoundsBehavior::Ignore; 
-  m_gf -> apply_accumulate<Interpolation::Kernel::Linear>(track, t_sig_start, t_sig_samp, num_samples, signal, oob_mode);
+  m_gf -> apply_accumulate<Interpolation::Kernel::Linear, double>(track, t_sig_start, t_sig_samp, num_samples, signal, oob_mode, weight);
 }