averaging_kernels/
averaging_kernels.rs

1use beambook::{geometry::{point,rectangle,ORIGIN},Color,Book,Page,Plot};
2use beamplot::{
3    plot::*,
4    font::Font,
5    homography::Homography,
6    maps::{Map,LinearMap,ExponentialMap}
7};
8
9fn main() {
10    let font = Font::new();
11    let mut bk = Book::new();
12    let mut pg = Page::new();
13
14    let origin = ORIGIN;
15
16    let dx = 100.0;
17    // let pressures = Array1::linspace(200.0,700.0,50);
18    // for &pressure0 in pressures.iter() {
19    let pressure0 = 350.0;
20    let mut pl = Plot::new();
21
22    let size = 10.0;
23    let p1 = point(0.0,200.0);
24    let p3 = dx*point(1.0,0.0);
25    let p2 = p1 + p3;
26    let y0 = 0.0;
27    let y1 = 1013.25;
28    let x0 = -0.5;
29    let x1 = 1.5;
30
31    let x_map = LinearMap::new(0.0,p3.x,x0,x1);
32    let y_map = ExponentialMap::new(0.0,p1.y,y0,y1);
33    let tick_spacing = size/4.0;
34    let ticks_y = try_ticks(10,&y_map,
35			    |pos,y| if pos == Position::Last { format!("hPa {}",y) } else { format!("{}",y) },
36			    tick_spacing).unwrap();
37    println!("Ticks Y: {:?}",ticks_y);
38    let ticks_x = try_ticks(5,&x_map,|_,x| format!("{:.1}",x),tick_spacing).unwrap();
39    println!("Ticks X: {:?}",ticks_x);
40
41    ruler(&font,size,origin,p1,true,false,&ticks_y).plot(&mut pl);
42    ruler(&font,size,p1,p2,false,true,&ticks_x).plot(&mut pl);
43    // ruler(&font,x0,x1,p1,p2,size,false,true,
44    // 	  |_,y| format!("{:.4}",y),
45    // 	  |x0,x1,dl| {
46    // 	      let m = ((1.0/dl).floor() as usize).max(2);
47    // 	      let mut ticks = Array2::zeros((m,2));
48    // 	      ticks.slice_mut(s![..,0]).assign(&Array1::linspace(0.0,1.0,m));
49    // 	      ticks.slice_mut(s![..,1]).assign(&Array1::linspace(x0,x1,m));
50    // 	      ticks
51    // 	  }).plot(&mut pl);
52    // let dx = size*point(20.0,0.0);
53    // ruler(&font,0.0,101325.0,ORIGIN+dx,p1+dx,size).plot(&mut pl);
54    // .plot(&mut pl);
55
56    if x0*x1 < 0.0 { pl.line(0x888,origin+x_map.inverse(0.0)*point(1.0,0.0),origin+p1+x_map.inverse(0.0)*point(1.0,0.0)); }
57    pl.rectangle(0xfff,rectangle(origin,p2));
58    let mut f1 = |p:f64| (-sq((p - 0.2*pressure0)/150.0)).exp() + 0.05*(p/25.0).cos();
59    let mut f2 = |p:f64| (-sq((p - 1.0*pressure0)/150.0)).exp() + 0.10*(p/30.0).cos();
60    let mut f3 = |p:f64| (-sq((p - 1.5*pressure0)/150.0)).exp() + 0.15*(p/35.0).cos();
61    let mut g = |f:&mut dyn Fn(f64)->f64,color:Color| {
62	curve(origin,p1,p2,
63	      &y_map,
64	      &x_map,
65	      1.0,
66	      false,
67	      color,
68	      f).plot(&mut pl)
69    };
70    g(&mut f1,0xf00);
71    g(&mut f2,0x0f0);
72    g(&mut f3,0x0ff);
73    bottom_align(&text_lines(&font, size, 0xfff, &["T"]).rc()).transformed(&Homography::translation(origin)).plot(&mut pl);
74
75    let origin = origin + point(size,0.0) + p3;
76    if x0*x1 < 0.0 { pl.line(0x888,origin+x_map.inverse(0.0)*point(1.0,0.0),origin+p1+x_map.inverse(0.0)*point(1.0,0.0)); }
77    pl.rectangle(0xfff,rectangle(origin,origin+p2));
78    let mut f1 = |p:f64| (-sq((p - 0.2*pressure0)/10.0)).exp() + 0.05*(p/25.0).cos();
79    let mut f2 = |p:f64| (-sq((p - 1.0*pressure0)/10.0)).exp() + 0.10*(p/30.0).cos();
80    let mut f3 = |p:f64| (-sq((p - 1.5*pressure0)/10.0)).exp() + 0.15*(p/35.0).cos();
81    let mut g = |f:&mut dyn Fn(f64)->f64,color:Color| {
82	curve(origin,origin+p1,origin+p2,
83	      &y_map,
84	      &x_map,
85	      1.0,
86	      false,
87	      color,
88	      f).plot(&mut pl)
89    };
90    g(&mut f1,0xf00);
91    g(&mut f2,0x0f0);
92    g(&mut f3,0x0ff);
93    ruler(&font,size,origin+p1,origin+p2,false,true,&ticks_x).plot(&mut pl);
94    bottom_align(&text_lines(&font, size, 0xfff, &["H_2O"]).rc()).transformed(&Homography::translation(origin)).plot(&mut pl);
95
96    let origin = origin + point(size,0.0) + p3;
97    if x0*x1 < 0.0 { pl.line(0x888,origin+x_map.inverse(0.0)*point(1.0,0.0),origin+p1+x_map.inverse(0.0)*point(1.0,0.0)); }
98    pl.rectangle(0xfff,rectangle(origin,origin + p2));
99    let mut f1 = |p:f64| (-sq((p - 0.2*pressure0)/50.0)).exp() + 0.15*(p/15.0).cos();
100    let mut f2 = |p:f64| (-sq((p - 1.0*pressure0)/50.0)).exp() + 0.05*(p/20.0).cos();
101    let mut f3 = |p:f64| (-sq((p - 1.5*pressure0)/50.0)).exp() + 0.10*(p/25.0).cos();
102    let mut f4 = |p:f64| (-sq((p - 1.2*pressure0)/30.0)).exp() + 0.12*(p/50.0).cos();
103    let mut g = |f:&mut dyn Fn(f64)->f64,color:Color| {
104	curve(origin,origin + p1,origin + p2,
105	      &y_map,
106	      &x_map,
107	      1.0,
108	      false,
109	      color,
110	      f).plot(&mut pl)
111    };
112    g(&mut f1,0xf00);
113    g(&mut f2,0x0f0);
114    g(&mut f3,0x0ff);
115    g(&mut f4,0xff0);
116    ruler(&font,size,origin+p1,origin+p2,false,true,&ticks_x).plot(&mut pl);
117    bottom_align(&text_lines(&font, size, 0xfff, &["CH_4"]).rc()).transformed(&Homography::translation(origin)).plot(&mut pl);
118
119    let subtitle =
120	left_align(&top_align(&legend(&font,size,
121				      &[(0xf00,"IASI"),
122					(0x0f0,"TROPOMI"),
123					(0x0ff,"IASI+L1(TROPOMI)"),
124					(0xff0,"IASI+L2(TROPOMI)")],
125				      0xfff).rc())).transformed(&Homography::translation(origin + p3));
126    subtitle.plot(&mut pl);
127
128    bottom_align(
129	&text_lines(&font,
130		    1.5*size,
131		    0xfff,
132		    &[
133			&format!("Averaging kernels, p_0={:.3} hPa",
134				 pressure0),
135			"Note that the kernels have been normalized by pressure",
136			"Error amplification: 3x,1.5x, Q=σAT^4"
137		    ]).rc())
138	.transformed(&Homography::translation(point(0.0,subtitle.area.a.y - 2.0*size))).plot(&mut pl);
139    
140    pg.plot(pl);
141    // }
142    bk.page(pg);
143    bk.save_to_file("traj.mpk").unwrap();
144}
145
146fn sq(x:f64)->f64 { x*x }