1use super::{CircularArc, Contour, ContourSegment, LineSegment};
2use crate::integrable::ComplexScalar;
3
4use nalgebra::ComplexField;
5use num_traits::FromPrimitive;
6
7#[derive(Clone, Copy, Debug, Eq, PartialEq)]
8pub enum IndentSide {
9 Left,
10 Right,
11}
12
13impl<F> Contour<F>
14where
15 F: ComplexScalar + FromPrimitive,
16{
17 pub fn indent(self, pole: F::Complex, radius: F, side: IndentSide, tolerance: F) -> Self {
18 let mut pieces = Vec::new();
19
20 for piece in self.pieces {
21 match piece {
22 ContourSegment::Line(line) => {
23 if let Some(indented) = line.indent_around_point(pole, radius, side, tolerance)
24 {
25 pieces.extend(indented);
26 } else {
27 pieces.push(ContourSegment::Line(line));
28 }
29 }
30
31 other => pieces.push(other),
32 }
33 }
34
35 Self { pieces }
36 }
37}
38
39impl<C> LineSegment<C> {
40 pub fn indent_around_point<F>(
41 &self,
42 pole: F::Complex,
43 radius: F,
44 side: IndentSide,
45 tolerance: F,
46 ) -> Option<Vec<ContourSegment<F>>>
47 where
48 F: ComplexScalar<Complex = C> + FromPrimitive,
49 C: ComplexField<RealField = F> + Copy,
50 {
51 let start = self.start();
52 let end = self.end();
53
54 let direction = end - start;
55 let length = direction.modulus();
56
57 if length == F::zero() {
58 return None;
59 }
60
61 let unit = direction / F::Complex::from_real(length);
62 let offset = pole - start;
63
64 let t = (offset * direction.conjugate()).real() / direction.modulus_squared();
65
66 if t <= F::zero() || t >= F::one() {
67 return None;
68 }
69
70 let closest = start + direction.scale(t);
71 let distance = (pole - closest).modulus();
72
73 if distance > tolerance * length {
74 return None;
75 }
76
77 if radius <= F::zero() || radius >= length {
78 return None;
79 }
80
81 let tangent_offset = unit * F::Complex::from_real(radius);
82
83 let entry = pole - tangent_offset;
84 let exit = pole + tangent_offset;
85
86 let _normal_angle_shift = match side {
87 IndentSide::Left => F::from_f64(std::f64::consts::FRAC_PI_2).unwrap(),
88 IndentSide::Right => -F::from_f64(std::f64::consts::FRAC_PI_2).unwrap(),
89 };
90
91 let theta0 = (entry - pole).argument();
92 let theta1 = match side {
93 IndentSide::Left => theta0 - F::from_f64(std::f64::consts::PI).unwrap(),
94 IndentSide::Right => theta0 + F::from_f64(std::f64::consts::PI).unwrap(),
95 };
96
97 Some(vec![
98 ContourSegment::Line(LineSegment::new(start, entry)),
99 ContourSegment::CircularArc(CircularArc::new(pole, radius, theta0, theta1)),
100 ContourSegment::Line(LineSegment::new(exit, end)),
101 ])
102 }
103}
104
105#[cfg(test)]
106mod contour_indent_tests {
107 use super::*;
108 use num_complex::Complex;
109
110 use crate::Infallible;
111 use crate::contour::ContourPiece;
112
113 const TOL: f64 = 1e-12;
114
115 fn assert_close(a: f64, b: f64) {
116 assert!((a - b).abs() < TOL, "expected {b}, got {a}");
117 }
118
119 fn assert_complex_close(a: Complex<f64>, b: Complex<f64>) {
120 assert_close(a.re, b.re);
121 assert_close(a.im, b.im);
122 }
123
124 #[test]
125 fn indent_replaces_line_through_pole_with_line_arc_line() {
126 let z0 = Complex::new(-1.0, 0.0);
127 let z1 = Complex::new(1.0, 0.0);
128 let pole = Complex::new(0.0, 0.0);
129
130 let contour =
131 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Left, 1e-10);
132
133 assert_eq!(contour.pieces().len(), 3);
134
135 match &contour.pieces()[0] {
136 ContourSegment::Line(line) => {
137 assert_complex_close(line.start(), z0);
138 assert_complex_close(line.end(), Complex::new(-0.1, 0.0));
139 }
140 _ => panic!("expected first piece to be a line"),
141 }
142
143 match &contour.pieces()[1] {
144 ContourSegment::CircularArc(arc) => {
145 assert_complex_close(arc.center(), pole);
146 assert_close(arc.radius(), 0.1);
147 }
148 _ => panic!("expected second piece to be an arc"),
149 }
150
151 match &contour.pieces()[2] {
152 ContourSegment::Line(line) => {
153 assert_complex_close(line.start(), Complex::new(0.1, 0.0));
154 assert_complex_close(line.end(), z1);
155 }
156 _ => panic!("expected third piece to be a line"),
157 }
158 }
159
160 #[test]
161 fn indent_does_nothing_when_pole_is_not_on_line() {
162 let z0 = Complex::new(-1.0, 0.0);
163 let z1 = Complex::new(1.0, 0.0);
164 let pole = Complex::new(0.0, 0.5);
165
166 let contour =
167 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Left, 1e-10);
168
169 assert_eq!(contour.pieces().len(), 1);
170 }
171
172 #[test]
173 fn left_indent_goes_through_upper_half_plane() {
174 let z0 = Complex::new(-1.0, 0.0);
175 let z1 = Complex::new(1.0, 0.0);
176 let pole = Complex::new(0.0, 0.0);
177
178 let contour =
179 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Left, 1e-10);
180
181 let arc = match &contour.pieces()[1] {
182 ContourSegment::CircularArc(arc) => arc,
183 _ => panic!("expected arc"),
184 };
185
186 let midpoint: Complex<f64> = arc.point(0.5);
187
188 assert!(midpoint.im > 0.0, "left indent should pass above the pole");
189 assert_close((midpoint - pole).norm(), 0.1);
190 }
191
192 #[test]
193 fn right_indent_goes_through_lower_half_plane() {
194 let z0 = Complex::new(-1.0, 0.0);
195 let z1 = Complex::new(1.0, 0.0);
196 let pole = Complex::new(0.0, 0.0);
197
198 let contour =
199 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Right, 1e-10);
200
201 let arc = match &contour.pieces()[1] {
202 ContourSegment::CircularArc(arc) => arc,
203 _ => panic!("expected arc"),
204 };
205
206 let midpoint: Complex<f64> = arc.point(0.5);
207
208 dbg!(&midpoint);
209
210 assert!(midpoint.im < 0.0, "right indent should pass below the pole");
211 assert_close((midpoint - pole).norm(), 0.1);
212 }
213
214 struct ComplexFn<F>(F);
215
216 impl<G> crate::Integrable for ComplexFn<G>
217 where
218 G: Fn(&Complex<f64>) -> Complex<f64>,
219 {
220 type Float = f64;
221 type Input = Complex<f64>;
222 type Output = Complex<f64>;
223
224 fn integrand(&self, z: &Complex<f64>) -> Complex<f64> {
225 self.0(z)
226 }
227 }
228
229 #[test]
230 fn upper_semicircle_around_origin_integrates_inverse_z_to_minus_i_pi() {
231 let gk = crate::core::GaussKronrod::<f64>::default();
232
233 let arc = CircularArc::new(Complex::new(0.0, 0.0), 1.0, std::f64::consts::PI, 0.0);
234
235 let f = Infallible(ComplexFn(|z: &Complex<f64>| Complex::new(1.0, 0.0) / *z));
236
237 let segment = gk
238 .integrate_piece(&f, &arc, crate::core::PathKey::new(0), false)
239 .unwrap();
240
241 assert_complex_close(segment.result, Complex::new(0.0, -std::f64::consts::PI));
242 }
243
244 #[test]
245 fn left_and_right_indents_have_opposite_inverse_z_arc_contributions() {
246 let gk = crate::core::GaussKronrod::<f64>::default();
247 let f = Infallible(ComplexFn(|z: &Complex<f64>| Complex::new(1.0, 0.0) / *z));
248
249 let z0 = Complex::new(-1.0, 0.0);
250 let z1 = Complex::new(1.0, 0.0);
251 let pole = Complex::new(0.0, 0.0);
252
253 let left =
254 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Left, 1e-10);
255
256 let right =
257 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Right, 1e-10);
258
259 let left_arc = match &left.pieces()[1] {
260 ContourSegment::CircularArc(arc) => *arc,
261 _ => panic!("expected arc"),
262 };
263
264 let right_arc = match &right.pieces()[1] {
265 ContourSegment::CircularArc(arc) => *arc,
266 _ => panic!("expected arc"),
267 };
268
269 let left_segment = gk
270 .integrate_piece(&f, &left_arc, crate::core::PathKey::new(0), false)
271 .unwrap();
272 let right_segment = gk
273 .integrate_piece(&f, &right_arc, crate::core::PathKey::new(0), false)
274 .unwrap();
275
276 assert_complex_close(
278 left_segment.result,
279 Complex::new(0.0, -std::f64::consts::PI),
280 );
281
282 assert_complex_close(
284 right_segment.result,
285 Complex::new(0.0, std::f64::consts::PI),
286 );
287 }
288}