Mercurial > repos > pointsource_algs / file diff

--- a/src/kernels/gaussian.rs	Thu Dec 01 23:07:35 2022 +0200
+++ b/src/kernels/gaussian.rs	Fri Oct 07 13:11:08 2022 +0300
@@ -16,6 +16,8 @@
     GlobalAnalysis,
     Weighted,
     Bounded,
+    Taylor2Model,
+    Taylor2ModelParams,
 };
 use alg_tools::mapping::Apply;
 use alg_tools::maputil::array_init;
@@ -224,6 +226,49 @@
     }
 }
 
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, R, C, S, const N : usize> Taylor2Model<F, N>
+for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
+where R : Constant<Type=F>,
+      C : Constant<Type=F>,
+      S : Constant<Type=F>,
+      Loc<Loc<F, 3>, N> : Product2Taylor<F, N> {
+
+    fn taylor2model(&self, y : Loc<F, N>) -> Taylor2ModelParams<F, N> {
+        let Convolution(ref ind,
+                        SupportProductFirst(ref cut,
+                                            ref gaussian)) = self;
+        let a = cut.r.value();
+        let b = ind.r.value();
+        let σ = gaussian.variance.value().sqrt();
+        let π = F::PI;
+        let t = F::SQRT_2 * σ;
+        let c = σ * (8.0/π).sqrt();
+        let sc = gaussian.scale();
+        
+        // This is just a product of one-dimensional versions
+        y.map(|x| {
+            let c1 = -(a.min(b + x)); //(-a).max(-x-b);
+            let c2 = a.min(b - x);
+            if c1 >= c2 {
+                Loc([0.0, 0.0, 0.0])
+            } else {
+                let s1 = c1 / t;
+                let s2 = c2 / t;
+                let e1 = F::cast_from(erf(s1.as_()));
+                let e2 = F::cast_from(erf(s2.as_()));
+                let d1 = F::FRAC_2_SQRT_PI * (-s1*s1).exp();
+                let d2 = F::FRAC_2_SQRT_PI * (-s2*s2).exp();
+                let dd1 = -2.0 * s1 * d1;
+                let dd2 = -2.0 * s2 * s2;
+                debug_assert!(e2 >= e1);
+                Loc([c * (e2 - e1), c * (d2 - d1), c * (dd2 - dd1)])
+            }
+        }).product2taylor_scaled(sc)
+    }
+
+}
+
 impl<F : Float, R, C, S, const N : usize>
 Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
 where R : Constant<Type=F>,