diff --git a/examples/biomes.png b/examples/biomes.png index 00437a7e..36b3c9bb 100644 Binary files a/examples/biomes.png and b/examples/biomes.png differ diff --git a/examples/weather.cpp b/examples/weather.cpp index cb16ac71..8f19d7c7 100644 --- a/examples/weather.cpp +++ b/examples/weather.cpp @@ -43,8 +43,8 @@ auto make_perlin(random::generator & rng, int min_octave, int max_octave, float void make_force_field(random::generator & rng, util::ndarray, 2> & result, float scale) { - auto noise_1 = make_perlin(rng, 2, 6, 2.f); - auto noise_2 = make_perlin(rng, 2, 6, 2.f); + auto noise_1 = make_perlin(rng, 3, 6, 2.f); + auto noise_2 = make_perlin(rng, 3, 6, 2.f); for (int y = 0; y < result.height(); ++y) { @@ -61,7 +61,7 @@ struct weather_app { static constexpr int N = 128; - const float dt = 4.f; + const float dt = 20.f; const float viscosity = 0.f; const float advection_magnification = 1.f; const float temperature_diffusion = 0.001f; @@ -71,13 +71,13 @@ struct weather_app const float coriolis = 0.f; const float coriolis_bands = 2.f; const float friction = 0.f; - const float slope_force = 0.04f; + const float slope_force = 4.f; const float vorticity_confinement = 0.f; - const float evaporation = 0.1f; - const float max_humidity_factor = 0.001f; + const float evaporation = 0.01f; + const float max_humidity_factor = 100.f; const float precipitation_factor = 0.0003f; const float force_field_amplitude = 0.00005f; - const float random_forces = 0.f; + const float random_forces = 1.f; const float force_field_switch_duration = 720.f * 7.5f; // 7.5 days const int force_field_switch_frames = std::round(force_field_switch_duration / dt); // const float friction_factor = 1.f - std::exp(- friction * dt); @@ -89,7 +89,7 @@ struct weather_app gfx::pixmap_rgba biomes_map; - float expected_temperature_at(int y) + float expected_temperature_at(int y) const { // float latitude = (y - N * 0.5f) * 2.f / N; // return std::cos(latitude * float(math::pi)); @@ -97,7 +97,7 @@ struct weather_app return temperature_income_at(y) * dt / (std::exp(cooling * dt) - 1.f); } - float temperature_income_at(int y) + float temperature_income_at(int y) const { // float latitude = (y - N * 0.5f) * 2.f / N; float latitude = y * 1.f / N; @@ -105,6 +105,11 @@ struct weather_app return heating * math::lerp(0.8f, 1.f, 1.f - std::abs(latitude)); } + int wrap(int i) const + { + return (i + N) % N; + } + weather_app(options const &, context const &) { simulation_box_ = {{{0.f, N}, {0.f, N}}}; @@ -141,6 +146,9 @@ struct weather_app float latitude = (N * 0.5f - y) * 2.f / N; velocity_(x, y) = random_velocity(rng) * 0.f + 0.f * math::vector{-std::cos(0.5f * float(math::pi) * latitude * coriolis_bands), 0.f}; temperature_(x, y) = expected_temperature_at(y); + + float max_humidity = std::max(0.f, temperature_(x, y) - 223.f) * max_humidity_factor; + humidity_(x, y) = max_humidity + dt * evaporation * std::max(0.f, temperature_(x, y) - 273.f) * (1.f - precipitation_factor * dt) / precipitation_factor / dt; } } @@ -214,6 +222,9 @@ struct weather_app } float const force_field_t = ((frame_ % force_field_switch_frames) + 0.5f) / force_field_switch_frames; + int xmin = periodic_x ? 0 : 1; + int xmax = periodic_x ? N : N - 1; // exclusive + // Temperature source for (int y = 0; y < N; ++y) { @@ -234,22 +245,14 @@ struct weather_app humidity_(x, y) += dt * evaporation * std::max(0.f, temperature_(x, y) - 273.f); // float discharge = std::min(humidity_(x, y), precipitation_factor * dt); - float discharge = humidity_(x, y) * precipitation_factor * dt; - // float max_humidity = std::max(0.f, temperature_(x, y) - 223.f) * max_humidity_factor; - // float discharge = std::max(0.f, humidity_(x, y) - max_humidity) * precipitation_factor * dt; + // float discharge = humidity_(x, y) * precipitation_factor * dt; + float max_humidity = std::max(0.f, temperature_(x, y) - 223.f) * max_humidity_factor; + float discharge = std::max(0.f, humidity_(x, y) - max_humidity) * precipitation_factor * dt; humidity_(x, y) -= discharge; precipitation_(x, y) = discharge; } } - auto wrap = [](int i) - { - return (i + N) % N; - }; - - int xmin = periodic_x ? 0 : 1; - int xmax = periodic_x ? N : N - 1; // exclusive - // Velocity & temperature advection for (int i = 0; i < N; ++i) { @@ -366,10 +369,14 @@ struct weather_app (std::max(0.f, terrain_(x, y + 1)) - std::max(0.f, terrain_(x, y - 1))) / 2.f, }; - [[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * math::dot(math::normalized(velocity_(x, y)), terrain_gradient)); - velocity_(x, y) *= std::min(1.f, slope_factor); + // [[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * math::dot(math::normalized(velocity_(x, y)), terrain_gradient)); + // velocity_(x, y) *= std::min(1.f, slope_factor); // velocity_(x, y) -= terrain_gradient * slope_force * dt; + [[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * std::pow(math::length(terrain_gradient), 4.f)); + // [[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * std::pow(std::max(0.f, terrain_(x, y)), 1.f)); + velocity_(x, y) *= slope_factor; + // Directional external force // velocity_(x, y)[1] += 0.001f * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands); // velocity_(x, y) += math::direction(frame_ * dt * 2.f * float(math::pi) / 10080.f) * 0.001f * dt; @@ -528,6 +535,68 @@ struct weather_app mouseover_cell = {x, y}; } + if (mouseover_cell && (state().mouse_button_down.contains(app::mouse_button::left) || state().mouse_button_down.contains(app::mouse_button::right))) + { + float delta = state().mouse_button_down.contains(app::mouse_button::left) ? 1.f : -1.f; + + int R = 4; + + for (int iy = -R; iy <= R; ++iy) + { + for (int ix = -R; ix <= R; ++ix) + { + int x = (*mouseover_cell)[0] + ix; + int y = (*mouseover_cell)[1] + iy; + + if (x >= 0 && x < N && y >= 0 && y < N) + { + auto d = math::vector{ix, iy} / 2.f; + terrain_(x, y) += delta * 0.05f * std::exp(- math::dot(d, d)); + } + } + } + } + + if (mouseover_cell && state().mouse_button_down.contains(app::mouse_button::middle)) + { + int R = 4; + + float average = 0.f; + int count = 0; + + for (int iy = -R; iy <= R; ++iy) + { + for (int ix = -R; ix <= R; ++ix) + { + int x = (*mouseover_cell)[0] + ix; + int y = (*mouseover_cell)[1] + iy; + + if (x >= 0 && x < N && y >= 0 && y < N) + { + average += terrain_(x, y); + count += 1; + } + } + } + + average /= count; + + for (int iy = -R; iy <= R; ++iy) + { + for (int ix = -R; ix <= R; ++ix) + { + int x = (*mouseover_cell)[0] + ix; + int y = (*mouseover_cell)[1] + iy; + + if (x >= 0 && x < N && y >= 0 && y < N) + { + auto d = math::vector{ix, iy} / 2.f; + terrain_(x, y) += (average - terrain_(x, y)) * 0.05f * std::exp(- math::dot(d, d)); + } + } + } + } + [[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){ @@ -604,10 +673,34 @@ struct weather_app color = gfx::blend(color, map_color(10000.f * pressure_(x, y), {0.f, 0.f, 1.f, 0.75f}, {1.f, 0.f, 0.f, 0.75f})); if (show_water_vapor_) - color = gfx::blend(color, map_color(humidity_(x, y), {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f})); + color = gfx::blend(color, map_color(humidity_(x, y) * 0.001f, {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f})); if (show_precipitation_) - color = gfx::blend(color, map_color(precipitation_(x, y), {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f})); + { + // color = gfx::blend(color, map_color(precipitation_(x, y), {1.f, 1.f, 1.f, -1.f}, {1.f, 1.f, 1.f, 1.f})); + + float alpha = 2.f / (1.f + std::exp(- precipitation_(x, y))) - 1.f; + + gfx::color_4f cloud_color{1.f, 1.f, 1.f, alpha * 0.875f}; + + if (y > 0 && y + 1 < N) + { + if (periodic_x || (x > 0 && x + 1 < N)) + { + math::vector gradient + { + (precipitation_(wrap(x + 1), y) - precipitation_(wrap(x - 1), y)) / 2.f, + (precipitation_(x, y + 1) - precipitation_(x, y - 1)) / 2.f, + }; + + auto normal = math::normalized(math::vector{-gradient[0], -gradient[1], 2.f}); + auto lightness = 0.5f + 0.5f * math::dot(normal, math::normalized(math::vector{1.f, 2.f, 3.f})); + cloud_color = gfx::dark(cloud_color, 1.f - lightness); + } + } + + color = gfx::blend(color, cloud_color); + } if (show_average_precipitation_) color = gfx::blend(color, map_color(average_precipitation_(x, y), {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f}));