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::contour::ContourPiece;
111
112 const TOL: f64 = 1e-12;
113
114 fn assert_close(a: f64, b: f64) {
115 assert!((a - b).abs() < TOL, "expected {b}, got {a}");
116 }
117
118 fn assert_complex_close(a: Complex<f64>, b: Complex<f64>) {
119 assert_close(a.re, b.re);
120 assert_close(a.im, b.im);
121 }
122
123 #[test]
124 fn indent_replaces_line_through_pole_with_line_arc_line() {
125 let z0 = Complex::new(-1.0, 0.0);
126 let z1 = Complex::new(1.0, 0.0);
127 let pole = Complex::new(0.0, 0.0);
128
129 let contour =
130 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Left, 1e-10);
131
132 assert_eq!(contour.pieces().len(), 3);
133
134 match &contour.pieces()[0] {
135 ContourSegment::Line(line) => {
136 assert_complex_close(line.start(), z0);
137 assert_complex_close(line.end(), Complex::new(-0.1, 0.0));
138 }
139 _ => panic!("expected first piece to be a line"),
140 }
141
142 match &contour.pieces()[1] {
143 ContourSegment::CircularArc(arc) => {
144 assert_complex_close(arc.center(), pole);
145 assert_close(arc.radius(), 0.1);
146 }
147 _ => panic!("expected second piece to be an arc"),
148 }
149
150 match &contour.pieces()[2] {
151 ContourSegment::Line(line) => {
152 assert_complex_close(line.start(), Complex::new(0.1, 0.0));
153 assert_complex_close(line.end(), z1);
154 }
155 _ => panic!("expected third piece to be a line"),
156 }
157 }
158
159 #[test]
160 fn indent_does_nothing_when_pole_is_not_on_line() {
161 let z0 = Complex::new(-1.0, 0.0);
162 let z1 = Complex::new(1.0, 0.0);
163 let pole = Complex::new(0.0, 0.5);
164
165 let contour =
166 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Left, 1e-10);
167
168 assert_eq!(contour.pieces().len(), 1);
169 }
170
171 #[test]
172 fn left_indent_goes_through_upper_half_plane() {
173 let z0 = Complex::new(-1.0, 0.0);
174 let z1 = Complex::new(1.0, 0.0);
175 let pole = Complex::new(0.0, 0.0);
176
177 let contour =
178 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Left, 1e-10);
179
180 let arc = match &contour.pieces()[1] {
181 ContourSegment::CircularArc(arc) => arc,
182 _ => panic!("expected arc"),
183 };
184
185 let midpoint: Complex<f64> = arc.point(0.5);
186
187 assert!(midpoint.im > 0.0, "left indent should pass above the pole");
188 assert_close((midpoint - pole).norm(), 0.1);
189 }
190
191 #[test]
192 fn right_indent_goes_through_lower_half_plane() {
193 let z0 = Complex::new(-1.0, 0.0);
194 let z1 = Complex::new(1.0, 0.0);
195 let pole = Complex::new(0.0, 0.0);
196
197 let contour =
198 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Right, 1e-10);
199
200 let arc = match &contour.pieces()[1] {
201 ContourSegment::CircularArc(arc) => arc,
202 _ => panic!("expected arc"),
203 };
204
205 let midpoint: Complex<f64> = arc.point(0.5);
206
207 dbg!(&midpoint);
208
209 assert!(midpoint.im < 0.0, "right indent should pass below the pole");
210 assert_close((midpoint - pole).norm(), 0.1);
211 }
212
213 struct ComplexFn<F>(F);
214
215 impl<G> crate::Integrable for ComplexFn<G>
216 where
217 G: Fn(&Complex<f64>) -> Complex<f64>,
218 {
219 type Float = f64;
220 type Input = Complex<f64>;
221 type Output = Complex<f64>;
222
223 fn integrand(&self, z: &Complex<f64>) -> Complex<f64> {
224 self.0(z)
225 }
226 }
227
228 #[test]
229 fn upper_semicircle_around_origin_integrates_inverse_z_to_minus_i_pi() {
230 let gk = crate::core::GaussKronrod::<f64>::default();
231
232 let arc = CircularArc::new(Complex::new(0.0, 0.0), 1.0, std::f64::consts::PI, 0.0);
233
234 let f = ComplexFn(|z: &Complex<f64>| Complex::new(1.0, 0.0) / *z);
235
236 let segment = gk
237 .integrate_piece(&f, &arc, crate::core::PathKey::new(0), false)
238 .unwrap();
239
240 assert_complex_close(segment.result, Complex::new(0.0, -std::f64::consts::PI));
241 }
242
243 #[test]
244 fn left_and_right_indents_have_opposite_inverse_z_arc_contributions() {
245 let gk = crate::core::GaussKronrod::<f64>::default();
246 let f = ComplexFn(|z: &Complex<f64>| Complex::new(1.0, 0.0) / *z);
247
248 let z0 = Complex::new(-1.0, 0.0);
249 let z1 = Complex::new(1.0, 0.0);
250 let pole = Complex::new(0.0, 0.0);
251
252 let left =
253 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Left, 1e-10);
254
255 let right =
256 Contour::piecewise_linear(vec![z0, z1]).indent(pole, 0.1, IndentSide::Right, 1e-10);
257
258 let left_arc = match &left.pieces()[1] {
259 ContourSegment::CircularArc(arc) => *arc,
260 _ => panic!("expected arc"),
261 };
262
263 let right_arc = match &right.pieces()[1] {
264 ContourSegment::CircularArc(arc) => *arc,
265 _ => panic!("expected arc"),
266 };
267
268 let left_segment = gk
269 .integrate_piece(&f, &left_arc, crate::core::PathKey::new(0), false)
270 .unwrap();
271 let right_segment = gk
272 .integrate_piece(&f, &right_arc, crate::core::PathKey::new(0), false)
273 .unwrap();
274
275 assert_complex_close(
277 left_segment.result,
278 Complex::new(0.0, -std::f64::consts::PI),
279 );
280
281 assert_complex_close(
283 right_segment.result,
284 Complex::new(0.0, std::f64::consts::PI),
285 );
286 }
287}