garcalc Public

Watch 0 Fork 0 Star 0
Rust · 35557 bytes Raw Blame History
1 //! 3D plotting engine
2 //!
3 //! Provides surface plotting with:
4 //! - Explicit surfaces: z = f(x, y)
5 //! - Parametric surfaces: (x(u,v), y(u,v), z(u,v))
6 //! - Wireframe and filled rendering
7 //! - Height-based colormaps
8 //! - Rotation and zoom
9
10 use cairo::Context;
11 use garcalc_cas::{Evaluator, Expr};
12
13 use crate::Color;
14
15 /// 3D camera state (spherical coordinates around target)
16 #[derive(Debug, Clone, Copy)]
17 pub struct Camera3D {
18 /// Horizontal angle (radians)
19 pub azimuth: f64,
20 /// Vertical angle (radians)
21 pub elevation: f64,
22 /// Distance from target
23 pub distance: f64,
24 /// Look-at point
25 pub target: (f64, f64, f64),
26 }
27
28 impl Default for Camera3D {
29 fn default() -> Self {
30 Self {
31 azimuth: 45.0_f64.to_radians(),
32 elevation: 30.0_f64.to_radians(),
33 distance: 25.0,
34 target: (0.0, 0.0, 0.0),
35 }
36 }
37 }
38
39 /// Camera preset positions
40 #[derive(Debug, Clone, Copy)]
41 pub enum CameraPreset {
42 Front,
43 Side,
44 Top,
45 Isometric,
46 }
47
48 impl Camera3D {
49 /// Apply a camera preset
50 pub fn apply_preset(&mut self, preset: CameraPreset) {
51 match preset {
52 CameraPreset::Front => {
53 self.azimuth = 0.0;
54 self.elevation = 0.0;
55 }
56 CameraPreset::Side => {
57 self.azimuth = std::f64::consts::FRAC_PI_2;
58 self.elevation = 0.0;
59 }
60 CameraPreset::Top => {
61 self.azimuth = 0.0;
62 self.elevation = 89.0_f64.to_radians();
63 }
64 CameraPreset::Isometric => {
65 self.azimuth = 45.0_f64.to_radians();
66 self.elevation = 30.0_f64.to_radians();
67 }
68 }
69 }
70
71 /// Rotate the camera by delta angles
72 pub fn rotate(&mut self, d_azimuth: f64, d_elevation: f64) {
73 self.azimuth += d_azimuth;
74 self.elevation =
75 (self.elevation + d_elevation).clamp(-89.0_f64.to_radians(), 89.0_f64.to_radians());
76 }
77
78 /// Zoom by a factor
79 pub fn zoom(&mut self, factor: f64) {
80 self.distance = (self.distance / factor).clamp(0.5, 2000.0);
81 }
82
83 /// Get camera position in world coordinates
84 pub fn position(&self) -> (f64, f64, f64) {
85 let cos_elev = self.elevation.cos();
86 let sin_elev = self.elevation.sin();
87 let cos_azim = self.azimuth.cos();
88 let sin_azim = self.azimuth.sin();
89
90 (
91 self.target.0 + self.distance * cos_elev * cos_azim,
92 self.target.1 + self.distance * cos_elev * sin_azim,
93 self.target.2 + self.distance * sin_elev,
94 )
95 }
96 }
97
98 /// 3D viewport bounds
99 #[derive(Debug, Clone, Copy)]
100 pub struct Viewport3D {
101 pub x_min: f64,
102 pub x_max: f64,
103 pub y_min: f64,
104 pub y_max: f64,
105 pub z_min: f64,
106 pub z_max: f64,
107 }
108
109 impl Default for Viewport3D {
110 fn default() -> Self {
111 Self {
112 x_min: -5.0,
113 x_max: 5.0,
114 y_min: -5.0,
115 y_max: 5.0,
116 z_min: -5.0,
117 z_max: 5.0,
118 }
119 }
120 }
121
122 /// Colormap for surface coloring
123 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
124 pub enum Colormap {
125 Viridis,
126 Plasma,
127 Coolwarm,
128 Grayscale,
129 }
130
131 impl Colormap {
132 /// Map a value in [0, 1] to a color
133 pub fn color(&self, t: f64) -> Color {
134 let t = t.clamp(0.0, 1.0);
135 match self {
136 Colormap::Viridis => viridis(t),
137 Colormap::Plasma => plasma(t),
138 Colormap::Coolwarm => coolwarm(t),
139 Colormap::Grayscale => {
140 let v = (t * 255.0) as u8;
141 Color {
142 r: v,
143 g: v,
144 b: v,
145 a: 255,
146 }
147 }
148 }
149 }
150
151 /// Cycle to the next colormap
152 pub fn next(self) -> Self {
153 match self {
154 Colormap::Viridis => Colormap::Plasma,
155 Colormap::Plasma => Colormap::Coolwarm,
156 Colormap::Coolwarm => Colormap::Grayscale,
157 Colormap::Grayscale => Colormap::Viridis,
158 }
159 }
160
161 /// Display name
162 pub fn name(&self) -> &'static str {
163 match self {
164 Colormap::Viridis => "Viridis",
165 Colormap::Plasma => "Plasma",
166 Colormap::Coolwarm => "Coolwarm",
167 Colormap::Grayscale => "Grayscale",
168 }
169 }
170 }
171
172 /// A 3D plottable surface
173 #[derive(Debug, Clone)]
174 pub enum Surface3D {
175 /// z = f(x, y)
176 Explicit {
177 expr: Expr,
178 x_var: String,
179 y_var: String,
180 },
181 /// Parametric: (x(u,v), y(u,v), z(u,v))
182 Parametric {
183 x_expr: Expr,
184 y_expr: Expr,
185 z_expr: Expr,
186 u_var: String,
187 v_var: String,
188 u_range: (f64, f64),
189 v_range: (f64, f64),
190 },
191 /// Spherical: r = f(theta, phi)
192 Spherical {
193 expr: Expr,
194 theta_var: String,
195 phi_var: String,
196 theta_range: (f64, f64),
197 phi_range: (f64, f64),
198 },
199 /// Cylindrical: r = f(theta, z)
200 Cylindrical {
201 expr: Expr,
202 theta_var: String,
203 z_var: String,
204 theta_range: (f64, f64),
205 z_range: (f64, f64),
206 },
207 /// Level surface: f(x,y,z) = c, rendered as z-plane contour slices
208 LevelSurface {
209 expr: Expr,
210 x_var: String,
211 y_var: String,
212 z_var: String,
213 level: f64,
214 },
215 }
216
217 /// Render mode for surfaces
218 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
219 pub enum RenderMode {
220 Wireframe,
221 Filled,
222 FilledWithWireframe,
223 }
224
225 impl RenderMode {
226 /// Cycle to the next render mode
227 pub fn next(self) -> Self {
228 match self {
229 RenderMode::FilledWithWireframe => RenderMode::Filled,
230 RenderMode::Filled => RenderMode::Wireframe,
231 RenderMode::Wireframe => RenderMode::FilledWithWireframe,
232 }
233 }
234
235 /// Display name
236 pub fn name(&self) -> &'static str {
237 match self {
238 RenderMode::Wireframe => "Wireframe",
239 RenderMode::Filled => "Filled",
240 RenderMode::FilledWithWireframe => "FilledWire",
241 }
242 }
243 }
244
245 /// Configuration for 3D plot appearance
246 #[derive(Debug, Clone)]
247 pub struct Plot3DConfig {
248 pub background: Color,
249 pub axis_color: Color,
250 pub wireframe_color: Color,
251 pub colormap: Colormap,
252 pub render_mode: RenderMode,
253 pub grid_lines: usize,
254 pub show_axes: bool,
255 pub show_labels: bool,
256 /// Surface transparency (0.0 = fully transparent, 1.0 = fully opaque)
257 pub surface_alpha: f64,
258 /// Whether to show coordinate plane grids
259 pub show_coord_planes: bool,
260 }
261
262 impl Default for Plot3DConfig {
263 fn default() -> Self {
264 Self {
265 background: Color {
266 r: 30,
267 g: 30,
268 b: 46,
269 a: 255,
270 },
271 axis_color: Color {
272 r: 166,
273 g: 173,
274 b: 200,
275 a: 255,
276 },
277 wireframe_color: Color {
278 r: 100,
279 g: 100,
280 b: 120,
281 a: 255,
282 },
283 colormap: Colormap::Viridis,
284 render_mode: RenderMode::FilledWithWireframe,
285 grid_lines: 40,
286 show_axes: true,
287 show_labels: true,
288 surface_alpha: 0.85,
289 show_coord_planes: false,
290 }
291 }
292 }
293
294 /// 3D graph state and renderer
295 pub struct Graph3D {
296 pub camera: Camera3D,
297 pub viewport: Viewport3D,
298 pub config: Plot3DConfig,
299 pub surfaces: Vec<Surface3D>,
300 }
301
302 impl Default for Graph3D {
303 fn default() -> Self {
304 Self {
305 camera: Camera3D::default(),
306 viewport: Viewport3D::default(),
307 config: Plot3DConfig::default(),
308 surfaces: Vec::new(),
309 }
310 }
311 }
312
313 impl Graph3D {
314 pub fn new() -> Self {
315 Self::default()
316 }
317
318 /// Add an explicit surface z = f(x, y)
319 pub fn add_explicit(&mut self, expr: Expr) {
320 self.surfaces.push(Surface3D::Explicit {
321 expr,
322 x_var: "x".to_string(),
323 y_var: "y".to_string(),
324 });
325 }
326
327 /// Add a parametric surface (x(u,v), y(u,v), z(u,v))
328 pub fn add_parametric(
329 &mut self,
330 x_expr: Expr,
331 y_expr: Expr,
332 z_expr: Expr,
333 u_range: (f64, f64),
334 v_range: (f64, f64),
335 ) {
336 self.surfaces.push(Surface3D::Parametric {
337 x_expr,
338 y_expr,
339 z_expr,
340 u_var: "u".to_string(),
341 v_var: "v".to_string(),
342 u_range,
343 v_range,
344 });
345 }
346
347 /// Cycle the render mode
348 pub fn cycle_render_mode(&mut self) {
349 self.config.render_mode = self.config.render_mode.next();
350 }
351
352 /// Add a spherical surface r = f(theta, phi)
353 pub fn add_spherical(&mut self, expr: Expr) {
354 self.surfaces.push(Surface3D::Spherical {
355 expr,
356 theta_var: "theta".to_string(),
357 phi_var: "phi".to_string(),
358 theta_range: (0.0, std::f64::consts::TAU),
359 phi_range: (0.0, std::f64::consts::PI),
360 });
361 }
362
363 /// Add a cylindrical surface r = f(theta, z)
364 pub fn add_cylindrical(&mut self, expr: Expr) {
365 self.surfaces.push(Surface3D::Cylindrical {
366 expr,
367 theta_var: "theta".to_string(),
368 z_var: "z".to_string(),
369 theta_range: (0.0, std::f64::consts::TAU),
370 z_range: (-5.0, 5.0),
371 });
372 }
373
374 /// Add a level surface f(x,y,z) = c
375 pub fn add_level_surface(&mut self, expr: Expr, level: f64) {
376 self.surfaces.push(Surface3D::LevelSurface {
377 expr,
378 x_var: "x".to_string(),
379 y_var: "y".to_string(),
380 z_var: "z".to_string(),
381 level,
382 });
383 }
384
385 /// Clear all surfaces
386 pub fn clear_surfaces(&mut self) {
387 self.surfaces.clear();
388 }
389
390 /// Reset camera to default
391 pub fn reset_camera(&mut self) {
392 self.camera = Camera3D::default();
393 }
394
395 /// Project 3D point to 2D screen coordinates
396 fn project(&self, x: f64, y: f64, z: f64, width: u32, height: u32) -> (f64, f64) {
397 let cam_pos = self.camera.position();
398
399 // View direction (camera looks at target)
400 let dx = self.camera.target.0 - cam_pos.0;
401 let dy = self.camera.target.1 - cam_pos.1;
402 let dz = self.camera.target.2 - cam_pos.2;
403 let dist = (dx * dx + dy * dy + dz * dz).sqrt();
404
405 // Normalize
406 let forward = (dx / dist, dy / dist, dz / dist);
407
408 // Up vector (world Z)
409 let world_up = (0.0, 0.0, 1.0);
410
411 // Right vector = forward x up
412 let right = (
413 forward.1 * world_up.2 - forward.2 * world_up.1,
414 forward.2 * world_up.0 - forward.0 * world_up.2,
415 forward.0 * world_up.1 - forward.1 * world_up.0,
416 );
417 let right_len = (right.0 * right.0 + right.1 * right.1 + right.2 * right.2).sqrt();
418 let right = (
419 right.0 / right_len,
420 right.1 / right_len,
421 right.2 / right_len,
422 );
423
424 // Actual up = right x forward
425 let up = (
426 right.1 * forward.2 - right.2 * forward.1,
427 right.2 * forward.0 - right.0 * forward.2,
428 right.0 * forward.1 - right.1 * forward.0,
429 );
430
431 // Point relative to camera
432 let px = x - cam_pos.0;
433 let py = y - cam_pos.1;
434 let pz = z - cam_pos.2;
435
436 // Project onto camera plane
437 let cam_x = px * right.0 + py * right.1 + pz * right.2;
438 let cam_y = px * up.0 + py * up.1 + pz * up.2;
439 let cam_z = px * forward.0 + py * forward.1 + pz * forward.2;
440
441 // Perspective projection
442 let scale = if cam_z > 0.1 {
443 self.camera.distance / cam_z
444 } else {
445 self.camera.distance / 0.1
446 };
447
448 // Screen coordinates
449 let aspect = width as f64 / height as f64;
450 let fov_scale = 0.8;
451 let sx = width as f64 / 2.0 + cam_x * scale * height as f64 * fov_scale / aspect;
452 let sy = height as f64 / 2.0 - cam_y * scale * height as f64 * fov_scale;
453
454 (sx, sy)
455 }
456
457 /// Render the 3D graph to a Cairo context
458 pub fn render(&self, ctx: &Context, width: u32, height: u32) {
459 // Background
460 set_color(ctx, self.config.background);
461 ctx.rectangle(0.0, 0.0, width as f64, height as f64);
462 let _ = ctx.fill();
463
464 // Draw axes
465 if self.config.show_axes {
466 self.draw_axes(ctx, width, height);
467 }
468
469 // Draw coordinate plane grids
470 if self.config.show_coord_planes {
471 self.draw_coord_planes(ctx, width, height);
472 }
473
474 // Draw surfaces
475 for surface in &self.surfaces {
476 self.draw_surface(ctx, surface, width, height);
477 }
478 }
479
480 fn draw_axes(&self, ctx: &Context, width: u32, height: u32) {
481 set_color(ctx, self.config.axis_color);
482 ctx.set_line_width(1.5);
483
484 let len = 6.0;
485
486 // X axis (red tint)
487 ctx.set_source_rgba(0.9, 0.4, 0.4, 1.0);
488 let (x0, y0) = self.project(0.0, 0.0, 0.0, width, height);
489 let (x1, y1) = self.project(len, 0.0, 0.0, width, height);
490 ctx.move_to(x0, y0);
491 ctx.line_to(x1, y1);
492 let _ = ctx.stroke();
493
494 // Y axis (green tint)
495 ctx.set_source_rgba(0.4, 0.9, 0.4, 1.0);
496 let (x1, y1) = self.project(0.0, len, 0.0, width, height);
497 ctx.move_to(x0, y0);
498 ctx.line_to(x1, y1);
499 let _ = ctx.stroke();
500
501 // Z axis (blue tint)
502 ctx.set_source_rgba(0.4, 0.4, 0.9, 1.0);
503 let (x1, y1) = self.project(0.0, 0.0, len, width, height);
504 ctx.move_to(x0, y0);
505 ctx.line_to(x1, y1);
506 let _ = ctx.stroke();
507
508 // Labels
509 if self.config.show_labels {
510 set_color(ctx, self.config.axis_color);
511 ctx.set_font_size(12.0);
512
513 let (lx, ly) = self.project(len + 0.5, 0.0, 0.0, width, height);
514 ctx.move_to(lx, ly);
515 let _ = ctx.show_text("X");
516
517 let (lx, ly) = self.project(0.0, len + 0.5, 0.0, width, height);
518 ctx.move_to(lx, ly);
519 let _ = ctx.show_text("Y");
520
521 let (lx, ly) = self.project(0.0, 0.0, len + 0.5, width, height);
522 ctx.move_to(lx, ly);
523 let _ = ctx.show_text("Z");
524 }
525 }
526
527 fn draw_surface(&self, ctx: &Context, surface: &Surface3D, width: u32, height: u32) {
528 match surface {
529 Surface3D::Explicit { expr, x_var, y_var } => {
530 self.draw_explicit_surface(ctx, expr, x_var, y_var, width, height);
531 }
532 Surface3D::Parametric {
533 x_expr,
534 y_expr,
535 z_expr,
536 u_var,
537 v_var,
538 u_range,
539 v_range,
540 } => {
541 self.draw_parametric_surface(
542 ctx, x_expr, y_expr, z_expr, u_var, v_var, *u_range, *v_range, width, height,
543 );
544 }
545 Surface3D::Spherical {
546 expr,
547 theta_var,
548 phi_var,
549 theta_range,
550 phi_range,
551 } => {
552 self.draw_spherical_surface(
553 ctx, expr, theta_var, phi_var, *theta_range, *phi_range, width, height,
554 );
555 }
556 Surface3D::Cylindrical {
557 expr,
558 theta_var,
559 z_var,
560 theta_range,
561 z_range,
562 } => {
563 self.draw_cylindrical_surface(
564 ctx, expr, theta_var, z_var, *theta_range, *z_range, width, height,
565 );
566 }
567 Surface3D::LevelSurface {
568 expr,
569 x_var,
570 y_var,
571 z_var,
572 level,
573 } => {
574 self.draw_level_surface(ctx, expr, x_var, y_var, z_var, *level, width, height);
575 }
576 }
577 }
578
579 fn draw_explicit_surface(
580 &self,
581 ctx: &Context,
582 expr: &Expr,
583 x_var: &str,
584 y_var: &str,
585 width: u32,
586 height: u32,
587 ) {
588 let mut evaluator = Evaluator::new();
589 let n = self.config.grid_lines;
590
591 let x_range = self.viewport.x_max - self.viewport.x_min;
592 let y_range = self.viewport.y_max - self.viewport.y_min;
593 let dx = x_range / n as f64;
594 let dy = y_range / n as f64;
595
596 // Sample the surface
597 let mut points: Vec<Vec<Option<(f64, f64, f64)>>> = Vec::with_capacity(n + 1);
598 let mut z_min = f64::INFINITY;
599 let mut z_max = f64::NEG_INFINITY;
600
601 for i in 0..=n {
602 let mut row = Vec::with_capacity(n + 1);
603 let x = self.viewport.x_min + i as f64 * dx;
604
605 for j in 0..=n {
606 let y = self.viewport.y_min + j as f64 * dy;
607
608 evaluator.set_var(x_var, Expr::Float(x));
609 evaluator.set_var(y_var, Expr::Float(y));
610
611 if let Ok(result) = evaluator.eval(expr) {
612 if let Ok(z) = expr_to_f64(&result) {
613 if z.is_finite() && z >= self.viewport.z_min && z <= self.viewport.z_max {
614 z_min = z_min.min(z);
615 z_max = z_max.max(z);
616 row.push(Some((x, y, z)));
617 } else {
618 row.push(None);
619 }
620 } else {
621 row.push(None);
622 }
623 } else {
624 row.push(None);
625 }
626 }
627 points.push(row);
628 }
629
630 // Avoid division by zero
631 if (z_max - z_min).abs() < 1e-10 {
632 z_max = z_min + 1.0;
633 }
634
635 let quads = self.build_quads(&points, n, z_min, z_max);
636 self.draw_quads(ctx, &quads, width, height);
637 }
638
639 fn draw_parametric_surface(
640 &self,
641 ctx: &Context,
642 x_expr: &Expr,
643 y_expr: &Expr,
644 z_expr: &Expr,
645 u_var: &str,
646 v_var: &str,
647 u_range: (f64, f64),
648 v_range: (f64, f64),
649 width: u32,
650 height: u32,
651 ) {
652 let mut evaluator = Evaluator::new();
653 let n = self.config.grid_lines;
654
655 let du = (u_range.1 - u_range.0) / n as f64;
656 let dv = (v_range.1 - v_range.0) / n as f64;
657
658 // Sample the surface
659 let mut points: Vec<Vec<Option<(f64, f64, f64)>>> = Vec::with_capacity(n + 1);
660 let mut z_min = f64::INFINITY;
661 let mut z_max = f64::NEG_INFINITY;
662
663 for i in 0..=n {
664 let mut row = Vec::with_capacity(n + 1);
665 let u = u_range.0 + i as f64 * du;
666
667 for j in 0..=n {
668 let v = v_range.0 + j as f64 * dv;
669
670 evaluator.set_var(u_var, Expr::Float(u));
671 evaluator.set_var(v_var, Expr::Float(v));
672
673 let x_result = evaluator.eval(x_expr);
674 let y_result = evaluator.eval(y_expr);
675 let z_result = evaluator.eval(z_expr);
676
677 if let (Ok(xr), Ok(yr), Ok(zr)) = (x_result, y_result, z_result) {
678 if let (Ok(x), Ok(y), Ok(z)) =
679 (expr_to_f64(&xr), expr_to_f64(&yr), expr_to_f64(&zr))
680 {
681 if x.is_finite() && y.is_finite() && z.is_finite() {
682 z_min = z_min.min(z);
683 z_max = z_max.max(z);
684 row.push(Some((x, y, z)));
685 } else {
686 row.push(None);
687 }
688 } else {
689 row.push(None);
690 }
691 } else {
692 row.push(None);
693 }
694 }
695 points.push(row);
696 }
697
698 if (z_max - z_min).abs() < 1e-10 {
699 z_max = z_min + 1.0;
700 }
701
702 let quads = self.build_quads(&points, n, z_min, z_max);
703 self.draw_quads(ctx, &quads, width, height);
704 }
705
706 fn draw_spherical_surface(
707 &self,
708 ctx: &Context,
709 expr: &Expr,
710 theta_var: &str,
711 phi_var: &str,
712 theta_range: (f64, f64),
713 phi_range: (f64, f64),
714 width: u32,
715 height: u32,
716 ) {
717 let mut evaluator = Evaluator::new();
718 let n = self.config.grid_lines;
719
720 let dtheta = (theta_range.1 - theta_range.0) / n as f64;
721 let dphi = (phi_range.1 - phi_range.0) / n as f64;
722
723 let mut points: Vec<Vec<Option<(f64, f64, f64)>>> = Vec::with_capacity(n + 1);
724 let mut z_min = f64::INFINITY;
725 let mut z_max = f64::NEG_INFINITY;
726
727 for i in 0..=n {
728 let mut row = Vec::with_capacity(n + 1);
729 let theta = theta_range.0 + i as f64 * dtheta;
730
731 for j in 0..=n {
732 let phi = phi_range.0 + j as f64 * dphi;
733
734 evaluator.set_var(theta_var, Expr::Float(theta));
735 evaluator.set_var(phi_var, Expr::Float(phi));
736
737 if let Ok(result) = evaluator.eval(expr) {
738 if let Ok(r) = expr_to_f64(&result) {
739 if r.is_finite() {
740 let x = r * phi.sin() * theta.cos();
741 let y = r * phi.sin() * theta.sin();
742 let z = r * phi.cos();
743 z_min = z_min.min(z);
744 z_max = z_max.max(z);
745 row.push(Some((x, y, z)));
746 } else {
747 row.push(None);
748 }
749 } else {
750 row.push(None);
751 }
752 } else {
753 row.push(None);
754 }
755 }
756 points.push(row);
757 }
758
759 if (z_max - z_min).abs() < 1e-10 {
760 z_max = z_min + 1.0;
761 }
762
763 let quads = self.build_quads(&points, n, z_min, z_max);
764 self.draw_quads(ctx, &quads, width, height);
765 }
766
767 fn draw_cylindrical_surface(
768 &self,
769 ctx: &Context,
770 expr: &Expr,
771 theta_var: &str,
772 z_var: &str,
773 theta_range: (f64, f64),
774 z_range: (f64, f64),
775 width: u32,
776 height: u32,
777 ) {
778 let mut evaluator = Evaluator::new();
779 let n = self.config.grid_lines;
780
781 let dtheta = (theta_range.1 - theta_range.0) / n as f64;
782 let dz = (z_range.1 - z_range.0) / n as f64;
783
784 let mut points: Vec<Vec<Option<(f64, f64, f64)>>> = Vec::with_capacity(n + 1);
785 let mut z_min_val = f64::INFINITY;
786 let mut z_max_val = f64::NEG_INFINITY;
787
788 for i in 0..=n {
789 let mut row = Vec::with_capacity(n + 1);
790 let theta = theta_range.0 + i as f64 * dtheta;
791
792 for j in 0..=n {
793 let z = z_range.0 + j as f64 * dz;
794
795 evaluator.set_var(theta_var, Expr::Float(theta));
796 evaluator.set_var(z_var, Expr::Float(z));
797
798 if let Ok(result) = evaluator.eval(expr) {
799 if let Ok(r) = expr_to_f64(&result) {
800 if r.is_finite() {
801 let x = r * theta.cos();
802 let y = r * theta.sin();
803 z_min_val = z_min_val.min(z);
804 z_max_val = z_max_val.max(z);
805 row.push(Some((x, y, z)));
806 } else {
807 row.push(None);
808 }
809 } else {
810 row.push(None);
811 }
812 } else {
813 row.push(None);
814 }
815 }
816 points.push(row);
817 }
818
819 if (z_max_val - z_min_val).abs() < 1e-10 {
820 z_max_val = z_min_val + 1.0;
821 }
822
823 let quads = self.build_quads(&points, n, z_min_val, z_max_val);
824 self.draw_quads(ctx, &quads, width, height);
825 }
826
827 fn draw_level_surface(
828 &self,
829 ctx: &Context,
830 expr: &Expr,
831 x_var: &str,
832 y_var: &str,
833 z_var: &str,
834 level: f64,
835 width: u32,
836 height: u32,
837 ) {
838 let mut evaluator = Evaluator::new();
839 let num_slices: usize = 30;
840 let grid_size: usize = 60;
841
842 let z_step = (self.viewport.z_max - self.viewport.z_min) / num_slices as f64;
843 let dx = (self.viewport.x_max - self.viewport.x_min) / grid_size as f64;
844 let dy = (self.viewport.y_max - self.viewport.y_min) / grid_size as f64;
845
846 ctx.set_line_width(1.5);
847
848 for slice in 0..num_slices {
849 let z = self.viewport.z_min + (slice as f64 + 0.5) * z_step;
850 let t = (z - self.viewport.z_min) / (self.viewport.z_max - self.viewport.z_min);
851 let color = self.config.colormap.color(t.clamp(0.0, 1.0));
852 ctx.set_source_rgba(
853 color.r as f64 / 255.0,
854 color.g as f64 / 255.0,
855 color.b as f64 / 255.0,
856 self.config.surface_alpha,
857 );
858
859 evaluator.set_var(z_var, Expr::Float(z));
860
861 // Evaluate f(x,y,z) - c on a 2D grid at this z
862 let mut values = vec![vec![0.0f64; grid_size + 1]; grid_size + 1];
863 for i in 0..=grid_size {
864 let x = self.viewport.x_min + i as f64 * dx;
865 evaluator.set_var(x_var, Expr::Float(x));
866 for j in 0..=grid_size {
867 let y = self.viewport.y_min + j as f64 * dy;
868 evaluator.set_var(y_var, Expr::Float(y));
869 if let Ok(result) = evaluator.eval(expr) {
870 if let Ok(v) = expr_to_f64(&result) {
871 values[i][j] = if v.is_finite() { v - level } else { f64::NAN };
872 } else {
873 values[i][j] = f64::NAN;
874 }
875 } else {
876 values[i][j] = f64::NAN;
877 }
878 }
879 }
880
881 // Marching squares on this z-plane
882 for i in 0..grid_size {
883 for j in 0..grid_size {
884 let x0 = self.viewport.x_min + i as f64 * dx;
885 let y0 = self.viewport.y_min + j as f64 * dy;
886 let x1 = x0 + dx;
887 let y1 = y0 + dy;
888
889 let v00 = values[i][j];
890 let v10 = values[i + 1][j];
891 let v01 = values[i][j + 1];
892 let v11 = values[i + 1][j + 1];
893
894 if v00.is_nan() || v10.is_nan() || v01.is_nan() || v11.is_nan() {
895 continue;
896 }
897
898 let s00 = v00 >= 0.0;
899 let s10 = v10 >= 0.0;
900 let s01 = v01 >= 0.0;
901 let s11 = v11 >= 0.0;
902
903 let case = (s00 as u8) | ((s10 as u8) << 1) | ((s01 as u8) << 2) | ((s11 as u8) << 3);
904
905 let interp = |va: f64, vb: f64| -> f64 {
906 if (va - vb).abs() < 1e-15 { 0.5 } else { va / (va - vb) }
907 };
908
909 let e_bottom = || { let t = interp(v00, v10); (x0 + t * dx, y0) };
910 let e_top = || { let t = interp(v01, v11); (x0 + t * dx, y1) };
911 let e_left = || { let t = interp(v00, v01); (x0, y0 + t * dy) };
912 let e_right = || { let t = interp(v10, v11); (x1, y0 + t * dy) };
913
914 let draw_line = |p1: (f64, f64), p2: (f64, f64)| {
915 let (sx1, sy1) = self.project(p1.0, p1.1, z, width, height);
916 let (sx2, sy2) = self.project(p2.0, p2.1, z, width, height);
917 ctx.move_to(sx1, sy1);
918 ctx.line_to(sx2, sy2);
919 };
920
921 match case {
922 0 | 15 => {}
923 1 | 14 => draw_line(e_bottom(), e_left()),
924 2 | 13 => draw_line(e_bottom(), e_right()),
925 3 | 12 => draw_line(e_left(), e_right()),
926 4 | 11 => draw_line(e_left(), e_top()),
927 5 | 10 => {
928 draw_line(e_bottom(), e_left());
929 draw_line(e_top(), e_right());
930 }
931 6 | 9 => draw_line(e_bottom(), e_top()),
932 7 | 8 => draw_line(e_top(), e_right()),
933 _ => {}
934 }
935 }
936 }
937 let _ = ctx.stroke();
938 }
939 }
940
941 fn draw_coord_planes(&self, ctx: &Context, width: u32, height: u32) {
942 ctx.set_line_width(0.5);
943 ctx.set_source_rgba(0.5, 0.5, 0.6, 0.2);
944
945 let step = 1.0;
946 let range = 5.0;
947
948 // XY plane (z=0)
949 let mut x = -range;
950 while x <= range {
951 let (sx0, sy0) = self.project(x, -range, 0.0, width, height);
952 let (sx1, sy1) = self.project(x, range, 0.0, width, height);
953 ctx.move_to(sx0, sy0);
954 ctx.line_to(sx1, sy1);
955 x += step;
956 }
957 let mut y = -range;
958 while y <= range {
959 let (sx0, sy0) = self.project(-range, y, 0.0, width, height);
960 let (sx1, sy1) = self.project(range, y, 0.0, width, height);
961 ctx.move_to(sx0, sy0);
962 ctx.line_to(sx1, sy1);
963 y += step;
964 }
965 let _ = ctx.stroke();
966 }
967
968 /// Build quads from a grid of sampled points (shared by explicit, parametric, spherical, cylindrical)
969 fn build_quads(
970 &self,
971 points: &[Vec<Option<(f64, f64, f64)>>],
972 n: usize,
973 z_min: f64,
974 z_max: f64,
975 ) -> Vec<Quad> {
976 let mut quads = Vec::new();
977 for i in 0..n {
978 for j in 0..n {
979 if let (Some(p00), Some(p10), Some(p11), Some(p01)) = (
980 points[i][j],
981 points[i + 1][j],
982 points[i + 1][j + 1],
983 points[i][j + 1],
984 ) {
985 let avg_z = (p00.2 + p10.2 + p11.2 + p01.2) / 4.0;
986 let t = (avg_z - z_min) / (z_max - z_min);
987 let color = self.config.colormap.color(t);
988
989 let cam_pos = self.camera.position();
990 let cx = (p00.0 + p10.0 + p11.0 + p01.0) / 4.0;
991 let cy = (p00.1 + p10.1 + p11.1 + p01.1) / 4.0;
992 let cz = (p00.2 + p10.2 + p11.2 + p01.2) / 4.0;
993 let depth = (cx - cam_pos.0).powi(2)
994 + (cy - cam_pos.1).powi(2)
995 + (cz - cam_pos.2).powi(2);
996
997 quads.push(Quad {
998 corners: [p00, p10, p11, p01],
999 color,
1000 depth,
1001 });
1002 }
1003 }
1004 }
1005 quads.sort_by(|a, b| {
1006 b.depth.partial_cmp(&a.depth).unwrap_or(std::cmp::Ordering::Equal)
1007 });
1008 quads
1009 }
1010
1011 /// Draw sorted quads with the current render mode and alpha
1012 fn draw_quads(&self, ctx: &Context, quads: &[Quad], width: u32, height: u32) {
1013 let alpha = self.config.surface_alpha;
1014 for quad in quads {
1015 let corners: Vec<(f64, f64)> = quad
1016 .corners
1017 .iter()
1018 .map(|p| self.project(p.0, p.1, p.2, width, height))
1019 .collect();
1020
1021 if self.config.render_mode != RenderMode::Wireframe {
1022 ctx.set_source_rgba(
1023 quad.color.r as f64 / 255.0,
1024 quad.color.g as f64 / 255.0,
1025 quad.color.b as f64 / 255.0,
1026 alpha,
1027 );
1028 ctx.move_to(corners[0].0, corners[0].1);
1029 for c in &corners[1..] {
1030 ctx.line_to(c.0, c.1);
1031 }
1032 ctx.close_path();
1033 let _ = ctx.fill();
1034 }
1035
1036 if self.config.render_mode != RenderMode::Filled {
1037 set_color(ctx, self.config.wireframe_color);
1038 ctx.set_line_width(0.5);
1039 ctx.move_to(corners[0].0, corners[0].1);
1040 for c in &corners[1..] {
1041 ctx.line_to(c.0, c.1);
1042 }
1043 ctx.close_path();
1044 let _ = ctx.stroke();
1045 }
1046 }
1047 }
1048 }
1049
1050 /// A quad face for depth sorting
1051 struct Quad {
1052 corners: [(f64, f64, f64); 4],
1053 color: Color,
1054 depth: f64,
1055 }
1056
1057 /// Set Cairo source color
1058 fn set_color(ctx: &Context, color: Color) {
1059 ctx.set_source_rgba(
1060 color.r as f64 / 255.0,
1061 color.g as f64 / 255.0,
1062 color.b as f64 / 255.0,
1063 color.a as f64 / 255.0,
1064 );
1065 }
1066
1067 /// Convert Expr to f64
1068 fn expr_to_f64(expr: &Expr) -> Result<f64, ()> {
1069 match expr {
1070 Expr::Integer(n) => Ok(*n as f64),
1071 Expr::Float(x) => Ok(*x),
1072 Expr::Rational(r) => Ok(r.to_f64()),
1073 _ => Err(()),
1074 }
1075 }
1076
1077 // Colormap implementations
1078
1079 fn viridis(t: f64) -> Color {
1080 // Approximation of viridis colormap
1081 let r = (0.267 + t * (0.329 + t * (1.421 + t * (-1.685 + t * 0.668)))).clamp(0.0, 1.0);
1082 let g = (0.004 + t * (1.260 + t * (-0.569 + t * 0.305))).clamp(0.0, 1.0);
1083 let b = (0.329 + t * (1.440 + t * (-2.814 + t * (2.768 - t * 0.723)))).clamp(0.0, 1.0);
1084 Color {
1085 r: (r * 255.0) as u8,
1086 g: (g * 255.0) as u8,
1087 b: (b * 255.0) as u8,
1088 a: 255,
1089 }
1090 }
1091
1092 fn plasma(t: f64) -> Color {
1093 // Approximation of plasma colormap
1094 let r = (0.050 + t * (2.735 + t * (-2.811 + t * 0.826))).clamp(0.0, 1.0);
1095 let g = (0.029 + t * (-0.278 + t * (2.149 + t * (-1.100)))).clamp(0.0, 1.0);
1096 let b = (0.533 + t * (0.667 + t * (-2.519 + t * 1.319))).clamp(0.0, 1.0);
1097 Color {
1098 r: (r * 255.0) as u8,
1099 g: (g * 255.0) as u8,
1100 b: (b * 255.0) as u8,
1101 a: 255,
1102 }
1103 }
1104
1105 fn coolwarm(t: f64) -> Color {
1106 // Blue (cool) to red (warm)
1107 let r = if t < 0.5 { 0.2 + t * 1.6 } else { 1.0 };
1108 let g = if t < 0.5 {
1109 0.2 + t * 1.2
1110 } else {
1111 0.8 - (t - 0.5) * 1.6
1112 };
1113 let b = if t < 0.5 { 1.0 } else { 1.0 - (t - 0.5) * 1.6 };
1114 Color {
1115 r: (r.clamp(0.0, 1.0) * 255.0) as u8,
1116 g: (g.clamp(0.0, 1.0) * 255.0) as u8,
1117 b: (b.clamp(0.0, 1.0) * 255.0) as u8,
1118 a: 255,
1119 }
1120 }
1121
1122 #[cfg(test)]
1123 mod tests {
1124 use super::*;
1125
1126 #[test]
1127 fn test_camera_default() {
1128 let cam = Camera3D::default();
1129 assert!((cam.azimuth - 45.0_f64.to_radians()).abs() < 0.01);
1130 assert!((cam.elevation - 30.0_f64.to_radians()).abs() < 0.01);
1131 }
1132
1133 #[test]
1134 fn test_camera_rotate() {
1135 let mut cam = Camera3D::default();
1136 cam.rotate(0.1, 0.1);
1137 assert!(cam.azimuth > 45.0_f64.to_radians());
1138 assert!(cam.elevation > 30.0_f64.to_radians());
1139 }
1140
1141 #[test]
1142 fn test_colormap_bounds() {
1143 let cm = Colormap::Viridis;
1144 let c0 = cm.color(0.0);
1145 let c1 = cm.color(1.0);
1146 assert!(c0.r < 100); // Dark at 0
1147 assert!(c1.g > 200); // Yellowish at 1
1148 }
1149 }
1150