Draft: Source estimation

common/data/element.h

@@ -182,4 +182,39 @@ Data<T_Physic> arg(const Data<T_Physic>& d)
     return out;
 };
 
+/**
+ * @brief Returns a vector ([f][s]) that tells for each shot what coefficient minimizes the L² misfit for that shot.
+ * If d0 is the computed value for a unit shot and d the reference data, we look for a complexe number "alpha" such that 
+ * alpha * d0 is the best approximation to d.
+ * It is given by alpha = <d0, d> / <d0, d0> with < , > the complex inner product
+*/
+template<Physic T_Physic>
+std::vector<std::vector<std::complex<double>>> leastSquareShotIntensity(const Data<T_Physic>& d0, const Data<T_Physic>& d)
+{
+    areCompatible(d0, d);
+    std::vector<std::vector<std::complex<double>>> alpha;
+    alpha.resize(d.nf());
+    for (auto &shots: alpha)
+        shots.resize(d.ns());
+    
+
+    for (unsigned int f = 0; f < d.nf(); f++)
+    {
+        for (unsigned int s = 0; s < d.ns(); s++)
+        {
+            std::complex<double> num = 0.;
+            std::complex<double> den = 0.;
+            
+            for (unsigned int r = 0; r < d.nr(s); r++)
+            {
+                num += conj<T_Physic>(d0.value(f,s,r)) * d.value(f,s,r);
+                den += conj<T_Physic>(d0.value(f,s,r)) * d0.value(f,s,r);
+            }
+
+            alpha[f][s] = num / den;
+        }
+    }
+    return alpha;
+};
+
 #endif // H_COMMON_DATA_ELEMENT