Rust

use std::{
    cmp::Reverse,
    collections::{HashMap, HashSet},
    f64::consts::PI,
};

use itertools::Itertools;
use libm::atan2;
use priority_queue::PriorityQueue;

fn l2(i: usize, j: usize) -> usize {
    i * i + j * j
}

fn erupt(data: &[Vec<char>], vi: usize, vj: usize, r: usize) -> usize {
    let r2 = r * r;
    data.iter()
        .enumerate()
        .flat_map(|(i, l)| {
            l.iter()
                .enumerate()
                .filter(move |&(j, v)| *v != '@' && l2(vi.abs_diff(i), vj.abs_diff(j)) <= r2)
        })
        .map(|(_, v)| (*v as u8 - b'0') as usize)
        .sum::<usize>()
}

fn find(data: &[Vec<char>], c: char) -> (usize, usize) {
    data.iter()
        .enumerate()
        .flat_map(|(i, l)| {
            l.iter()
                .enumerate()
                .filter(|&(_, v)| *v == c)
                .map(move |(j, _)| (i, j))
        })
        .exactly_one()
        .unwrap()
}

pub fn solve_part_1(input: &str) -> String {
    let data = input
        .lines()
        .map(|l| l.chars().collect::<Vec<_>>())
        .collect::<Vec<_>>();
    let (vi, vj) = find(&data, '@');
    erupt(&data, vi, vj, 10).to_string()
}

pub fn solve_part_2(input: &str) -> String {
    let data = input
        .lines()
        .map(|l| l.chars().collect::<Vec<_>>())
        .collect::<Vec<_>>();
    let (vi, vj) = find(&data, '@');
    let h = data.len();
    let w = data[0].len();
    let max_r = vi.max(vj).max(h - vi - 1).max(w - vj - 1);
    let mut last_eruption = 0;
    let mut max_eruption = 0;
    let mut max_eruption_r = 0;
    for r in 1..=max_r {
        let eruption = erupt(&data, vi, vj, r);
        let de = eruption - last_eruption;
        if de > max_eruption {
            max_eruption = de;
            max_eruption_r = r;
        }
        last_eruption = eruption;
    }
    (max_eruption_r * max_eruption).to_string()
}

pub fn solve_part_3(input: &str) -> String {
    let data = input
        .lines()
        .map(|l| l.chars().collect::<Vec<_>>())
        .collect::<Vec<_>>();
    let width = data[0].len();
    let height = data.len();
    let (vi, vj) = find(&data, '@');
    let (si, sj) = find(&data, 'S');

    let azimuth = |i: usize, j: usize| atan2(i as f64 - vi as f64, j as f64 - vj as f64);

    let small_rot = |az1: f64, az2: f64| {
        let d = az1 - az2;
        if d > PI {
            d - 2. * PI
        } else if d < -PI {
            d + 2. * PI
        } else {
            d
        }
    };

    let solve = |radius: usize| {
        let r2 = radius * radius;
        let time_limit = ((radius + 1) * 30) as i64;
        let mut queue = PriorityQueue::new();
        let mut rotations = HashMap::new();
        rotations.insert((si, sj, false), 0f64);
        let mut visited = HashSet::new();
        queue.push((si, sj, false), Reverse(0));
        while let Some(((i, j, rotated), Reverse(time))) = queue.pop() {
            if time >= time_limit {
                break;
            }
            visited.insert((i, j, rotated));
            let az = azimuth(i, j);
            let rotation = rotations[&(i, j, rotated)];
            for (di, dj) in [(-1, 0), (1, 0), (0, -1), (0, 1)] {
                let (ni, nj) = (i.wrapping_add_signed(di), j.wrapping_add_signed(dj));
                if ni >= height || nj >= width {
                    continue;
                }
                if l2(ni.abs_diff(vi), nj.abs_diff(vj)) <= r2 {
                    continue;
                }
                let is_rotated = if let Some(previous_rotation) = rotations.get(&(ni, nj, false)) {
                    let rotation = rotation + small_rot(azimuth(ni, nj), az);
                    (rotation - previous_rotation).abs() > 6.
                } else {
                    false
                };
                if (ni, nj, is_rotated) == (si, sj, true) {
                    return Some(time);
                }
                if visited.contains(&(ni, nj, is_rotated)) {
                    continue;
                }
                let new_time: i64 = time + (data[ni][nj] as i8 - '0' as i8) as i64;
                let should_update =
                    match queue.push_increase((ni, nj, is_rotated), Reverse(new_time)) {
                        None => true,
                        Some(Reverse(t)) => t > new_time,
                    };
                if should_update {
                    rotations.insert(
                        (ni, nj, is_rotated),
                        rotation + small_rot(azimuth(ni, nj), az),
                    );
                };
            }
        }
        None
    };
    let (radius, time) = (1..(width.min(height) / 2))
        .map(|radius| (radius, solve(radius)))
        .filter(|(_, s)| s.is_some())
        .min()
        .unwrap();
    (radius as i64 * time.unwrap()).to_string()
}

Python

Just putting the solution to part 3 here, because these solutions are pretty long. Took me a while to figure out that an optimal return path from the bottom of the volcano burn can bleed into the left quadrant.

from collections import deque
import heapq as hq
import re

DIRECTIONS_CARDINAL = [(-1, 0), (0, 1), (1, 0), (0, -1)]
DIRECTIONS_ALL = [
    *DIRECTIONS_CARDINAL,               # cardinal
    (-1, -1), (-1, 1), (1, -1), (1, 1)  # diagonal
]

# yields all 8 possible neighbors
# lava spreads in all 8 directions
def yield_all_neighbors(i, j, m, n):
    for di, dj in DIRECTIONS_ALL:
        ni, nj = i + di, j + dj
        if 0 <= ni < m and 0 <= nj < n:
            yield ni, nj

# yields only the 4 possible cardinal neighbors
# the dragonduck can only move in cardinal directions
def yield_cardinal_neighbors(i, j, m, n):
    for di, dj in DIRECTIONS_CARDINAL:
        ni, nj = i + di, j + dj
        if 0 <= ni < m and 0 <= nj < n:
            yield ni, nj

# finds a special character in the grid
def find_char(data: str, ch: str):
    i = 0
    for row in data.splitlines():
        j = row.find(ch)
        if j != -1:
            return (i, j)
        i += 1
    assert False, f"Character {ch} not found in grid"

# returns squared distance between two points
def get_dist_sq(a, b):
    x1, y1 = a
    x2, y2 = b
    return (x2 - x1) ** 2 + (y2 - y1) ** 2

# Generator to burn cells permanently one step at a time
def bfs_burn(grid, volcano):
    m, n = len(grid), len(grid[0])
    queue = deque([volcano])
    visited = set()

    step = 0
    while queue:
        nq = len(queue)
        for _ in range(nq):
            i, j = queue.popleft()
            if (i, j) in visited:
                continue
            
            # skip cells outside the current step radius
            #   but do not mark them as visited yet
            if get_dist_sq((i, j), volcano) > (step * step):
                continue
            visited.add((i, j))

            # burn cell permanently
            grid[i][j] = -1

            for ni, nj in yield_all_neighbors(i, j, m, n):
                if (ni, nj) in visited:
                    continue
                queue.append((ni, nj))        
        
        step += 1
        # yield here to allow enclosing the volcano after each step
        yield step

def part3(data: str):
    # initialize variables
    start = find_char(data, 'S')
    volcano = find_char(data, '@')
    data = re.sub(r'[@S]', '0', data)
    grid = [[int(c) for c in row] for row in data.splitlines()]
    m, n = len(grid), len(grid[0])

    # initialize lava generator
    lava_spreader = bfs_burn(grid, volcano)

    # Finds shortest path from start to dest within time limit
    # Dijkstra's algorithm with constraints
    def find_shortest_path(start: tuple[int, int], dest: tuple[int, int], volcano_col: int, time_limit, map_half):
        candidates = [(0, start[0], start[1])]
        visited = set()

        while candidates:
            cost, i, j = hq.heappop(candidates)            
            if (i, j) in visited:
                continue
            visited.add((i, j))

            if (i, j) == dest:
                # if we are processing the destination, we have found the shortest path to it
                return cost
            
            if cost > time_limit:
                # if we have exceeded the time limit, prune this path
                continue
            
            # remember: the dragonduck can only move in cardinal directions
            for ni, nj in yield_cardinal_neighbors(i, j, m, n):
                # skip visited or burned cells
                if (ni, nj) in visited or grid[ni][nj] == -1:
                    continue
                # if processing left half, skip paths going to the right of the volcano
                if map_half == 'left' and nj > volcano_col:
                    continue
                # if processing right half, skip paths going to the left of the volcano (ONLY for bottom quadrant)
                # allow moving into the left half in the top quadrant (!!!)
                elif map_half == 'right' and (ni > m // 2 and nj < volcano_col):
                    continue

                # add new candidate path
                hq.heappush(candidates, (cost + grid[ni][nj], ni, nj))

        # we couldn't reach the destination within time limit
        # fail by returning the maximum possible cost
        return 9 * m * n

    # try increasing lava spreads
    while True:
        # spread lava one more step
        step = next(lava_spreader)
        # calculate time limit
        time = step * 30
        # our paths from the start must connect at the bottom of the lava burn
        burn_bottom_boundary = (volcano[0] + step, volcano[1])

        time_taken = find_shortest_path(start, burn_bottom_boundary, volcano[1], time, 'left')
        if time_taken >= time:
            # if we already exceeded time, quit early
            continue
        time_taken += find_shortest_path(burn_bottom_boundary, start, volcano[1], time - time_taken, 'right')
        if time_taken >= time:
            # if we exceeded time, try next lava spread
            continue
        
        # we successfully enclosed the volcano within time
        burn_radius = step - 1
        return time_taken * burn_radius