diff --git a/examples/weather_v2.cpp b/examples/weather_v2.cpp
new file mode 100644
index 00000000..dc6095c9
--- /dev/null
+++ b/examples/weather_v2.cpp
@@ -0,0 +1,880 @@
+#include <psemek/app/application_base.hpp>
+#include <psemek/app/default_application_factory.hpp>
+#include <psemek/gfx/painter.hpp>
+#include <psemek/gfx/gl.hpp>
+#include <psemek/math/camera.hpp>
+#include <psemek/math/gauss.hpp>
+#include <psemek/random/generator.hpp>
+#include <psemek/random/device.hpp>
+#include <psemek/random/uniform_ball.hpp>
+#include <psemek/pcg/perlin.hpp>
+#include <psemek/pcg/fractal.hpp>
+#include <psemek/util/ndarray.hpp>
+#include <psemek/util/clock.hpp>
+#include <psemek/log/log.hpp>
+#include <psemek/io/file_stream.hpp>
+
+using namespace psemek;
+
+auto make_perlin(random::generator & rng, int min_octave, int max_octave, float power)
+{
+	std::vector<pcg::perlin<float, 2>> octaves;
+	std::vector<float> weights;
+
+	random::uniform_sphere_vector_distribution<float, 2> random_vector{};
+
+	for (int octave = min_octave; octave < max_octave; ++octave)
+	{
+		int size = 1 << octave;
+		util::ndarray<math::vector<float, 2>, 2> gradients({size + 1, size + 1});
+		for (auto & g : gradients)
+			g = random_vector(rng);
+		octaves.emplace_back(std::move(gradients));
+		weights.push_back(std::pow(power, - octave));
+	}
+
+	float weight_sum = 0.f;
+	for (auto w : weights)
+		weight_sum += w;
+	for (auto & w : weights)
+		w /= weight_sum;
+
+	return pcg::fractal<pcg::perlin<float, 2>>(std::move(octaves), std::move(weights));
+}
+
+void make_random_vector_field(random::generator & rng, util::ndarray<math::vector<float, 2>, 2> & result, int min_octave, int max_octave, float scale)
+{
+	auto noise_1 = make_perlin(rng, min_octave, max_octave, 2.f);
+	auto noise_2 = make_perlin(rng, min_octave, max_octave, 2.f);
+
+	for (int y = 0; y < result.height(); ++y)
+	{
+		for (int x = 0; x < result.width(); ++x)
+		{
+			math::point p{(x + 0.5f) / result.height(), (y + 0.5f) / result.width()};
+			result(x, y) = math::normalized(math::vector{2.f * noise_1(p) - 1.f, 2.f * noise_2(p) - 1.f}) * scale;
+		}
+	}
+}
+
+template <typename T>
+void add(util::ndarray<T, 2> & dst, util::ndarray<T, 2> const & src)
+{
+	auto src_begin = src.begin();
+	auto src_end = src.end();
+	auto dst_begin = dst.begin();
+	for (; src_begin != src_end;)
+		*dst_begin++ += *src_begin++;
+}
+
+template <typename T>
+void scale(util::ndarray<T, 2> & array, float factor)
+{
+	for (auto & v : array)
+		v *= factor;
+}
+
+float do_lerp(float x, float y, float t)
+{
+	return math::lerp(x, y, t);
+}
+
+math::vector<float, 2> do_lerp(math::vector<float, 2> const & x, math::vector<float, 2> const & y, float t)
+{
+	return math::lerp(length(x), length(y), t) * math::slerp(math::normalized(x), math::normalized(y), t);
+}
+
+template <typename T>
+void lerp(util::ndarray<T, 2> & dst, util::ndarray<T, 2> const & src0, util::ndarray<T, 2> const & src1, float t)
+{
+	auto src0_begin = src0.begin();
+	auto src0_end = src0.end();
+	auto src1_begin = src1.begin();
+	auto dst_begin = dst.begin();
+
+	for (; src0_begin != src0_end;)
+		*dst_begin++ = do_lerp(*src0_begin++, *src1_begin++, t);
+}
+
+template <typename T>
+T sample_bilinear(util::ndarray<T, 2> const & array, math::point<float, 2> p)
+{
+	p[0] = math::clamp(p[0] - 0.5f, {0.f, array.width() - 1.f});
+	p[1] = math::clamp(p[1] - 0.5f, {0.f, array.height() - 1.f});
+
+	int ix = std::min<int>(array.width() - 2, std::floor(p[0]));
+	int iy = std::min<int>(array.height() - 2, std::floor(p[1]));
+
+	float tx = p[0] - ix;
+	float ty = p[1] - iy;
+
+	return math::lerp(
+		math::lerp(array(ix + 0, iy + 0), array(ix + 1, iy + 0), tx),
+		math::lerp(array(ix + 0, iy + 1), array(ix + 1, iy + 1), tx),
+		ty
+	);
+}
+
+std::string season_str(int season)
+{
+	switch (season)
+	{
+	case 0: return "Winter";
+	case 1: return "Spring";
+	case 2: return "Summer";
+	case 3: return "Autumn";
+	}
+
+	return "???";
+}
+
+std::string day_night_str(int day_night)
+{
+	switch (day_night)
+	{
+	case 0: return "Night";
+	case 1: return "Day";
+	}
+
+	return "???";
+}
+
+struct climate_snapshot
+{
+	util::ndarray<math::vector<float, 2>, 2> wind_velocity;
+	util::ndarray<float, 2> temperature;
+	util::ndarray<float, 2> humidity;
+};
+
+static constexpr int N = 128;
+
+struct solver
+{
+	util::ndarray<float, 2> const & terrain;
+	random::generator & rng;
+	int season;
+	int day_night;
+
+	climate_snapshot & snapshot;
+
+	int solve_velocity_iterations = 0;
+	int max_solve_velocity_iterations = N / 2;
+
+	int solve_temperature_iterations = 0;
+	int max_solve_temperature_iterations = N / 4;
+
+	int solve_humidity_iterations = 0;
+	int max_solve_humidity_iterations = 3 * N;
+
+	util::ndarray<math::vector<float, 2>, 2> initial_random_wind_velocity;
+	util::ndarray<float, 2> pressure;
+	util::ndarray<float, 2> new_temperature;
+	util::ndarray<float, 2> new_humidity;
+
+	int stage = 0;
+	bool finished = false;
+
+	solver(util::ndarray<float, 2> const & terrain,	random::generator & rng, int season, int day_night, climate_snapshot & snapshot)
+		: terrain(terrain)
+		, rng(rng)
+		, season(season)
+		, day_night(day_night)
+		, snapshot(snapshot)
+	{
+		initial_random_wind_velocity.resize({N, N});
+		snapshot.wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
+		pressure.resize({N, N}, 0.f);
+		snapshot.temperature.resize({N, N}, 0.f);
+		new_temperature.resize({N, N}, 0.f);
+		snapshot.humidity.resize({N, N}, 0.f);
+		new_humidity.resize({N, N}, 0.f);
+
+		make_random_vector_field(rng, initial_random_wind_velocity, 3, 6, 0.005f);
+
+		for (int y = 0; y < N; ++y)
+			for (int x = 0; x < N; ++x)
+				if (terrain(x, y) < 0.f)
+					initial_random_wind_velocity(x, y) *= 5.f;
+				else
+					initial_random_wind_velocity(x, y) *= std::exp(- 2.f * terrain(x, y));
+
+		snapshot.wind_velocity = initial_random_wind_velocity.copy();
+
+		for (int y = 0; y < N; ++y)
+			for (int x = 0; x < N; ++x)
+				snapshot.temperature(x, y) = expected_temperature(y);
+	}
+
+	float expected_temperature(int y)
+	{
+		return math::lerp(40.f, -20.f, (y + 0.5f) / N) // Base latitude temperature in spring/autumn
+			- 15.f * std::cos(season * math::pi * 0.5f) // Seasonal offset of +/-15C
+			- 5.f * std::cos(day_night * math::pi) // Day-night fluctuation of +/-5C
+			;
+	}
+
+	float expected_humidity(int x, int y)
+	{
+		if (terrain(x, y) < 0.f)
+			return 1.f;
+		else
+			return 0.f;
+	}
+
+	void step()
+	{
+		if (stage == 0)
+		{
+			for (int y = 1; y + 1 < N; ++y)
+			{
+				for (int x = 1; x + 1 < N; ++x)
+				{
+					float divergence = 0.f
+						+ (initial_random_wind_velocity(x + 1, y)[0] - initial_random_wind_velocity(x - 1, y)[0]) * 0.5f
+						+ (initial_random_wind_velocity(x, y + 1)[1] - initial_random_wind_velocity(x, y + 1)[1]) * 0.5f
+						;
+
+					// (sum - 4.f * p(x,y)) = divergence
+					// sum - 4.f * p = divergence
+					// p = (sum - divergence) / 4.f
+					pressure(x, y) = (pressure(x - 1, y) + pressure(x + 1, y) + pressure(x, y - 1) + pressure(x, y + 1) - divergence) / 4.f;
+				}
+			}
+
+			for (int y = 1; y + 1 < N; ++y)
+			{
+				for (int x = 1; x + 1 < N; ++x)
+				{
+					math::vector pressure_gradient
+					{
+						(pressure(x + 1, y) - pressure(x - 1, y)) * 0.5f,
+						(pressure(x, y + 1) - pressure(x, y - 1)) * 0.5f,
+					};
+					snapshot.wind_velocity(x, y) = initial_random_wind_velocity(x, y) - pressure_gradient;
+
+					if (terrain(x, y) > 0.f)
+						snapshot.wind_velocity(x, y) *= std::exp(- 2.f * terrain(x, y));
+				}
+			}
+
+			++solve_velocity_iterations;
+			if (solve_velocity_iterations == max_solve_velocity_iterations)
+				++stage;
+			return;
+		}
+
+		if (stage == 1)
+		{
+			for (int y = 0; y < N; ++y)
+			{
+				for (int x = 0; x < N; ++x)
+				{
+					float step_time = 100.f;
+
+					math::point p{x + 0.5f, y + 0.5f};
+					float advection = sample_bilinear(snapshot.temperature, p - snapshot.wind_velocity(x, y) * step_time);
+
+					float heating_speed = 0.001f;
+					float heating_delta = (expected_temperature(y) - snapshot.temperature(x, y)) * (- std::expm1(- heating_speed * step_time));
+
+					new_temperature(x, y) = heating_delta + advection;
+				}
+			}
+			std::swap(snapshot.temperature, new_temperature);
+
+			for (int y = 1; y + 1 < N; ++y)
+			{
+				for (int x = 1; x + 1 < N; ++x)
+				{
+					float average = (0.f
+						+ snapshot.temperature(x - 1, y)
+						+ snapshot.temperature(x + 1, y)
+						+ snapshot.temperature(x, y - 1)
+						+ snapshot.temperature(x, y + 1)
+					) * 0.25f;
+
+					new_temperature(x, y) = math::lerp(snapshot.temperature(x, y), average, 0.5f);
+				}
+			}
+			std::swap(snapshot.temperature, new_temperature);
+
+			++solve_temperature_iterations;
+			if (solve_temperature_iterations >= max_solve_temperature_iterations)
+				++stage;
+			return;
+		}
+
+		if (stage == 2)
+		{
+			for (int y = 0; y < N; ++y)
+			{
+				for (int x = 0; x < N; ++x)
+				{
+					float step_time = 100.f;
+
+					math::point p{x + 0.5f, y + 0.5f};
+					float advection = sample_bilinear(snapshot.humidity, p - snapshot.wind_velocity(x, y) * step_time);
+
+					float evaporation_speed = 0.0001f;
+					float evaporation_delta = (expected_humidity(x, y) - snapshot.humidity(x, y)) * (- std::expm1(- evaporation_speed * step_time));
+
+					new_humidity(x, y) = evaporation_delta + advection;
+				}
+			}
+			std::swap(snapshot.humidity, new_humidity);
+
+			for (int y = 1; y + 1 < N; ++y)
+			{
+				for (int x = 1; x + 1 < N; ++x)
+				{
+					float average = (0.f
+						+ snapshot.humidity(x - 1, y)
+						+ snapshot.humidity(x + 1, y)
+						+ snapshot.humidity(x, y - 1)
+						+ snapshot.humidity(x, y + 1)
+					) * 0.25f;
+
+					new_humidity(x, y) = math::lerp(snapshot.humidity(x, y), average, 0.5f);
+				}
+			}
+			std::swap(snapshot.humidity, new_humidity);
+
+			++solve_humidity_iterations;
+			if (solve_humidity_iterations >= max_solve_humidity_iterations)
+			{
+				++stage;
+				finished = true;
+			}
+			return;
+		}
+	}
+};
+
+struct weather_app
+	: app::application_base
+{
+	random::generator rng{random::device{}};
+	// random::generator rng{0, 0};
+
+	gfx::painter painter;
+
+	math::box<float, 2> simulation_box;
+
+	bool paused = true;
+	bool show_velocity = false;
+	bool show_temperature = false;
+	bool show_humidity = false;
+	bool show_biomes = true;
+
+	gfx::pixmap_rgba biomes_map;
+	util::ndarray<float, 2> terrain;
+
+	// 0 - Winter
+	// 1 - Spring
+	// 2 - Summer
+	// 3 - Autumn
+	int season = 0;
+
+	// 0 - Night
+	// 1 - Day
+	int day_night = 0;
+
+	climate_snapshot snapshots[4][2];
+	climate_snapshot average;
+	climate_snapshot display_snapshot;
+
+	bool need_update_display_snapshot = false;
+
+	std::optional<struct solver> solver;
+
+	float display_season = 0.f;
+
+	util::clock<> frame_clock;
+
+	weather_app(options const &, context const & context)
+	{
+		simulation_box = {{{0.f, N}, {0.f, N}}};
+
+		terrain.resize({N, N}, 0.f);
+		auto heightmap = gfx::read_image<std::uint8_t>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "heightmap_seed_3.png"});
+		for (int y = 0; y < N; ++y)
+			for (int x = 0; x < N; ++x)
+				terrain(x, y) = ((heightmap(x, y) / 255.f) * 2048.f - 512.f) / 1024.f;
+
+		biomes_map = gfx::read_image<gfx::color_rgba>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "biomes.png"});
+
+		context.vsync(false);
+	}
+
+	void on_event(app::key_event const & event) override
+	{
+		application_base::on_event(event);
+
+		if (event.down && event.key == app::keycode::SPACE)
+			paused ^= true;
+
+		if (event.down && event.key == app::keycode::V)
+			show_velocity ^= true;
+
+		if (event.down && event.key == app::keycode::T)
+			show_temperature ^= true;
+
+		if (event.down && event.key == app::keycode::W)
+			show_humidity ^= true;
+
+		if (event.down && event.key == app::keycode::B)
+			show_biomes ^= true;
+	}
+
+	void update() override
+	{
+		float const dt = frame_clock.restart().count();
+
+		if (solver && solver->finished)
+		{
+			solver.reset();
+			++day_night;
+			if (day_night == 2)
+			{
+				day_night = 0;
+				++season;
+			}
+
+			if (season == 4)
+			{
+				compute_average();
+				need_update_display_snapshot = true;
+			}
+		}
+
+		if (!solver && season < 4)
+			solver.emplace(terrain, rng, season, day_night, snapshots[season][day_night]);
+
+		if (!paused && solver)
+			solver->step();
+
+		if (state().key_down.contains(app::keycode::LEFT))
+		{
+			display_season -= 0.5f * dt;
+			need_update_display_snapshot = true;
+		}
+
+		if (state().key_down.contains(app::keycode::RIGHT))
+		{
+			display_season += 0.5f * dt;
+			need_update_display_snapshot = true;
+		}
+
+		while (display_season < 0.f)
+			display_season += 1.f;
+		while (display_season >= 1.f)
+			display_season -= 1.f;
+	}
+
+	void compute_average()
+	{
+		average.wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
+		average.temperature.resize({N, N}, 0.f);
+		average.humidity.resize({N, N}, 0.f);
+
+		for (auto const & season : snapshots)
+			for (auto const & s : season)
+			{
+				add(average.wind_velocity, s.wind_velocity);
+				add(average.temperature, s.temperature);
+				add(average.humidity, s.humidity);
+			}
+
+		scale(average.wind_velocity, 1.f / 8.f);
+		scale(average.temperature, 1.f / 8.f);
+		scale(average.humidity, 1.f / 8.f);
+	}
+
+	void update_display_snapshot()
+	{
+		display_snapshot.wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
+		display_snapshot.temperature.resize({N, N}, 0.f);
+		display_snapshot.humidity.resize({N, N}, 0.f);
+
+		int i = std::min<int>(3, std::floor(4.f * display_season));
+		int j = (i + 1) % 4;
+		float t = (4.f * display_season) - i;
+
+		auto & s0 = snapshots[i][0];
+		auto & s1 = snapshots[j][0];
+
+		lerp(display_snapshot.wind_velocity, s0.wind_velocity, s1.wind_velocity, t);
+		lerp(display_snapshot.temperature, s0.temperature, s1.temperature, t);
+		lerp(display_snapshot.humidity, s0.humidity, s1.humidity, t);
+	}
+
+	void present() override
+	{
+		if (need_update_display_snapshot)
+		{
+			update_display_snapshot();
+			need_update_display_snapshot = false;
+		}
+
+		gl::ClearColor(0.f, 0.f, 0.f, 0.f);
+		gl::Clear(gl::COLOR_BUFFER_BIT);
+
+		float const aspect_ratio = state().size[0] * 1.f / state().size[1];
+		math::box<float, 2> view_box = math::expand(simulation_box, 1.f);
+		if (view_box[0].length() / view_box[1].length() > aspect_ratio)
+			view_box[1] = math::expand(view_box[1], (view_box[0].length() / aspect_ratio - view_box[1].length()) / 2.f);
+		else
+			view_box[0] = math::expand(view_box[0], (view_box[1].length() * aspect_ratio - view_box[0].length()) / 2.f);
+
+		std::optional<math::vector<int, 2>> mouseover_cell;
+		{
+			auto mouse = math::lerp(view_box, math::vector{state().mouse[0] * 1.f / state().size[0], 1.f - state().mouse[1] * 1.f / state().size[1]});
+			int x = std::floor(mouse[0]);
+			int y = std::floor(mouse[1]);
+
+			if (x >= 0 && x < N && y >= 0 && y < N)
+				mouseover_cell = {x, y};
+		}
+
+		[[maybe_unused]] float const pixel_size = view_box[0].length() / state().size[0];
+
+		auto map_color = [](float value, gfx::color_4f const & negative, gfx::color_4f const & positive){
+			return math::lerp(negative, positive, 1.f / (1.f + std::exp(- value)));
+		};
+
+		auto map_color_positive = [](float value, gfx::color_4f const & zero, gfx::color_4f const & positive){
+			return math::lerp(zero, positive, 2.f / (1.f + std::exp(- value)) - 1.f);
+		};
+
+		auto map_temperature = [&](float temperature)
+		{
+			return map_color(0.5f * std::floor(temperature / 10.f), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f});
+		};
+
+		auto map_humidity = [&](float humidity)
+		{
+			return map_color_positive(std::round(10.f * humidity), {0.f, 1.f, 1.f, 0.f}, {0.f, 1.f, 1.f, 0.75f});
+		};
+
+		auto map_biome = [this](float temperature, float precipitation)
+		{
+			auto x = math::clamp<int>(math::unlerp({0.f, 0.5f}, precipitation) * biomes_map.width() , {0, biomes_map.width()  - 1});
+			auto y = math::clamp<int>(math::unlerp({-10.f, 30.f}, temperature  ) * biomes_map.height(), {0, biomes_map.height() - 1});
+
+			return gfx::to_colorf(biomes_map(x, y));
+		};
+
+		auto & avg_snapshot = (season < 4) ? snapshots[season][day_night] : average;
+		auto & snapshot = (season < 4) ? snapshots[season][day_night] : display_snapshot;
+
+		for (int y = 0; y < N; ++y)
+		{
+			for (int x = 0; x < N; ++x)
+			{
+				gfx::color_4f color = gfx::color_4f::zero();
+
+				bool const is_water = terrain(x, y) < 0.f;
+				bool const is_land = !is_water;
+
+				if (is_water)
+					color = map_color(8.f * terrain(x, y), {0.f, 0.f, 0.125f, 1.f}, {0.f, 1.f, 1.5f, 1.f});
+				else if (show_biomes)
+				{
+					float temperature = avg_snapshot.temperature(x, y) - std::max(0.f, terrain(x, y)) * 30.f;
+					color = map_biome(temperature, avg_snapshot.humidity(x, y));
+				}
+				else
+					color = gfx::color_4f{0.3f, 0.5f, 0.1f, 1.f};
+
+				if (is_land && x > 0 && x + 1 < N && y > 0 && y + 1 < N)
+				{
+					math::vector terrain_gradient
+					{
+						(terrain(x + 1, y) - terrain(x - 1, y)) / 2.f,
+						(terrain(x, y + 1) - terrain(x, y - 1)) / 2.f,
+					};
+
+					auto terrain_normal = math::normalized(math::vector{-terrain_gradient[0], -terrain_gradient[1], 0.125f});
+
+					float lightness = 0.5f + 0.5f * math::dot(terrain_normal, math::normalized(math::vector{1.f, 2.f, 3.f}));
+
+					color = gfx::dark(color, 1.f - lightness);
+				}
+
+				if (show_temperature)
+				{
+					color = gfx::blend(color, map_temperature(snapshot.temperature(x, y)));
+				}
+
+				if (show_humidity)
+				{
+					color = gfx::blend(color, map_humidity(snapshot.humidity(x, y)));
+				}
+
+				painter.rect({{{x, x + 1.f}, {y, y + 1.f}}}, gfx::to_coloru8(color));
+			}
+		}
+
+		if (show_velocity)
+		{
+			for (int y = 0; y < N; ++y)
+			{
+				for (int x = 0; x < N; ++x)
+				{
+					math::point center{x + 0.5f, y + 0.5f};
+					auto v = snapshot.wind_velocity(x, y);
+					if (auto l = math::length(v); l > 0.f)
+					{
+						float const magnification = 200.f;
+						float const max_length = 1.5f;
+						v *= 0.5f * max_length * (1.f - std::exp(- magnification * l)) / l;
+					}
+					auto n = math::ort(v) * 0.3f;
+
+					auto color = gfx::color_4f{1.f, 1.f, 1.f, 1.f};
+					if (!show_temperature)
+					{
+						float t = 0.2f * snapshot.temperature(x, y);
+						if (t >= 0.f)
+							color = map_color_positive( t, color, {1.f, 0.5f, 0.125f, 1.f});
+						else
+							color = map_color_positive(-t, color, {0.125f, 0.5f, 1.f, 1.f});
+					}
+
+					painter.triangle(center - v - n, center - v + n, center + v, gfx::to_coloru8(color));
+				}
+			}
+		}
+
+		int row = 0;
+		auto push_text = [&](std::string const & text) mutable
+		{
+			painter.text(view_box.corner(0.f, 1.f) - math::vector{0.f, row * pixel_size * 2.f * 12.f}, text, {.scale = {2.f * pixel_size, - 2.f * pixel_size}, .x = gfx::painter::x_align::left, .y = gfx::painter::y_align::top, .c = {255, 255, 255, 255}});
+			++row;
+		};
+
+		if (solver)
+		{
+			push_text(std::format("{} {}", season_str(season), day_night_str(day_night)));
+
+			if (solver->stage >= 0)
+				push_text(std::format("Velocity {}/{}", solver->solve_velocity_iterations, solver->max_solve_velocity_iterations));
+			else
+				push_text("");
+
+			if (solver->stage >= 1)
+				push_text(std::format("Temperature {}/{}", solver->solve_temperature_iterations, solver->max_solve_temperature_iterations));
+			else
+				push_text("");
+
+			if (solver->stage >= 2)
+				push_text(std::format("Humidity {}/{}", solver->solve_humidity_iterations, solver->max_solve_humidity_iterations));
+			else
+				push_text("");
+
+			push_text("");
+		}
+
+		if (mouseover_cell)
+		{
+			int x = (*mouseover_cell)[0];
+			int y = (*mouseover_cell)[1];
+			painter.rect({{{x, x + 1.f}, {y, y + 1.f}}}, {255, 255, 255, 127});
+
+			push_text(std::format("{} {}", x, y));
+			push_text(std::format("V = {:.3f}", length(snapshot.wind_velocity(x, y)) * 1000.f));
+			push_text(std::format("H = {:.3f}", terrain(x, y) * 1024.f));
+			push_text(std::format("T = {:.3f}", snapshot.temperature(x, y)));
+			push_text(std::format("W = {:.3f}", snapshot.humidity(x, y)));
+			push_text("");
+		}
+
+		if (season == 4)
+		{
+			++row;
+
+			float axis_width = 300.f;
+			float axis_height = 24.f;
+
+			float axis_left = view_box[0].min + 30.f * pixel_size;
+			float axis_bottom = view_box[1].max - row * pixel_size * 2.f * 12.f - axis_height * pixel_size;
+
+			auto line = [&](float x0, float y0, float x1, float y1, float w, gfx::color_rgba const & c)
+			{
+				painter.line(
+					{axis_left + x0 * axis_width * pixel_size, axis_bottom + y0 * axis_height * pixel_size},
+					{axis_left + x1 * axis_width * pixel_size, axis_bottom + y1 * axis_height * pixel_size},
+					pixel_size * w, c, false
+				);
+			};
+
+			line(0.f, 0.5f, 1.f, 0.5f, 2.f, {255, 255, 255, 127});
+			line(display_season, 0.f, display_season, 1.f, 4.f, {255, 255, 255, 255});
+
+			++row;
+			++row;
+		}
+
+		if (mouseover_cell && season == 4)
+		{
+			int x = (*mouseover_cell)[0];
+			int y = (*mouseover_cell)[1];
+
+			float plot_width = 300.f;
+			float plot_height = 150.f;
+
+			float plot_left = view_box[0].min + 30.f * pixel_size;
+			float plot_bottom = view_box[1].max - row * pixel_size * 2.f * 12.f - plot_height * 1.5f * pixel_size;
+
+			auto line = [&](float x0, float y0, float x1, float y1, gfx::color_rgba const & c)
+			{
+				painter.line(
+					{plot_left + x0 * plot_width * pixel_size, plot_bottom + y0 * plot_height * pixel_size},
+					{plot_left + x1 * plot_width * pixel_size, plot_bottom + y1 * plot_height * pixel_size},
+					pixel_size * 2.f, c, false
+				);
+			};
+
+			auto begin_plot = [&]{
+				painter.rect({{{plot_left, plot_left + plot_width * pixel_size}, {plot_bottom, plot_bottom + plot_height * pixel_size}}}, {255, 255, 255, 31});
+
+				for (int x = 0; x <= 4; ++x)
+					line(x / 4.f, 0.f, x / 4.f, 1.f, {255, 255, 255, (x == 0 || x == 4) ? 127 : 31});
+			};
+
+			auto draw_plot = [&](std::vector<float> const & ys, gfx::color_rgba const & c)
+			{
+				std::vector<float> xs(5);
+				for (int i = 0; i <= 4; ++i)
+					xs[i] = i / 4.f;
+
+				auto system = math::matrix<float, 16, 16>::zero();
+				auto coeffs = math::vector<float, 16>::zero();
+
+				for (int i = 0; i < 4; ++i)
+				{
+					float x0 = xs[i];
+					system[2 * i + 0][4 * i + 0] = 1.f;
+					system[2 * i + 0][4 * i + 1] = x0;
+					system[2 * i + 0][4 * i + 2] = x0 * x0;
+					system[2 * i + 0][4 * i + 3] = x0 * x0 * x0;
+					coeffs[2 * i + 0] = ys[i];
+
+					float x1 = xs[i + 1];
+					system[2 * i + 1][4 * i + 0] = 1.f;
+					system[2 * i + 1][4 * i + 1] = x1;
+					system[2 * i + 1][4 * i + 2] = x1 * x1;
+					system[2 * i + 1][4 * i + 3] = x1 * x1 * x1;
+					coeffs[2 * i + 1] = ys[i + 1];
+				}
+
+				for (int i = 0; i < 4; ++i)
+				{
+					int j = (i + 1) % 4;
+					float x0 = xs[i + 1];
+					float x1 = xs[j];
+
+					system[8 + 2 * i + 0][4 * i + 1] = 1.f;
+					system[8 + 2 * i + 0][4 * i + 2] = 2.f * x0;
+					system[8 + 2 * i + 0][4 * i + 3] = 3.f * x0 * x0;
+					system[8 + 2 * i + 0][4 * j + 1] = - 1.f;
+					system[8 + 2 * i + 0][4 * j + 2] = - 2.f * x1;
+					system[8 + 2 * i + 0][4 * j + 3] = - 3.f * x1 * x1;
+
+					system[8 + 2 * i + 1][4 * i + 2] = 2.f;
+					system[8 + 2 * i + 1][4 * i + 3] = 6.f * x0;
+					system[8 + 2 * i + 1][4 * j + 2] = - 2.f;
+					system[8 + 2 * i + 1][4 * j + 3] = - 6.f * x1;
+				}
+
+				math::gauss(system, coeffs);
+
+				auto eval = [&](int i, float x)
+				{
+					return 0.f
+						+ coeffs[4 * i + 0]
+						+ coeffs[4 * i + 1] * x
+						+ coeffs[4 * i + 2] * x * x
+						+ coeffs[4 * i + 3] * x * x * x
+						;
+				};
+
+				for (int i = 0; i < 4; ++i)
+				{
+					int D = 8;
+					for (int j = 0; j < D; ++j)
+					{
+						float x0 = (i + (j + 0.f) / D) / 4.f;
+						float x1 = (i + (j + 1.f) / D) / 4.f;
+						line(x0, eval(i, x0), x1, eval(i, x1), c);
+					}
+
+					if (i == std::min<int>(3, std::floor(4.f * display_season)))
+					{
+						float x = display_season;
+						float y = eval(i, x);
+						painter.circle({plot_left + x * plot_width * pixel_size, plot_bottom + y * plot_height * pixel_size}, 4.f * pixel_size, c);
+					}
+				}
+			};
+
+			// Temperature
+			{
+				begin_plot();
+
+				for (int y = 0; y <= 10; ++y)
+				{
+					line(0.f, y / 10.f, 1.f, y / 10.f, {255, 255, 255, (y == 5) ? 127 : 31});
+				}
+
+				std::vector<float> ys[2];
+
+				for (int d = 0; d < 2; ++d)
+				{
+					for (int s = 0; s <= 4; ++s)
+					{
+						float value = 0.5f + snapshots[s % 4][d].temperature(x, y) / 100.f;
+						ys[d].push_back(value);
+					}
+				}
+
+				draw_plot(ys[0], gfx::color_rgba{0, 127, 255, 255});
+				draw_plot(ys[1], gfx::color_rgba{255, 127, 0, 255});
+			}
+
+			plot_bottom -= (plot_height + 50.f) * pixel_size;
+
+			// Humidity
+			{
+				begin_plot();
+
+				for (int y = 0; y <= 5; ++y)
+				{
+					line(0.f, y / 5.f, 1.f, y / 5.f, {255, 255, 255, (y == 0) ? 127 : 31});
+				}
+
+				std::vector<float> ys;
+
+				for (int s = 0; s <= 4; ++s)
+				{
+					float value = (snapshots[s % 4][0].humidity(x, y) + snapshots[s % 4][1].humidity(x, y)) / 2.f;
+					ys.push_back(value);
+				}
+
+				draw_plot(ys, {0, 255, 255, 255});
+			}
+		}
+
+		painter.render(math::orthographic_camera{view_box}.transform());
+	}
+};
+
+namespace psemek::app
+{
+
+	std::unique_ptr<application::factory> make_application_factory()
+	{
+		return default_application_factory<weather_app>({.name = "Weather simulation test"});
+	}
+
+}