fork download
copy 
  1. import math
  2. import numpy as np
  3.  
  4. def total_horizontal_displacement(alpha0, speeds, widths, boundary_normal_angle):
  5.     """
  6.     Given an initial angle alpha0 (measured relative to the boundary normal),
  7.     propagate angles through each layer via Snell's Law and compute
  8.     the total horizontal displacement from start to finish.
  9.  
  10.     Parameters:
  11.     -----------
  12.     alpha0 : float
  13.         Initial angle (radians) in the first layer relative to the boundary normal.
  14.     speeds : list of float
  15.         speeds[i] = travel speed in layer i
  16.     widths : list of float
  17.         widths[i] = thickness of layer i in the boundary-normal direction
  18.     boundary_normal_angle : float
  19.         Angle (radians) of the boundary normal relative to the horizontal (east direction).
  20.         For example, if the boundary lines themselves are at 45° to the horizontal,
  21.         the boundary normal might also be ±45°, depending on orientation.
  22.  
  23.     Returns:
  24.     --------
  25.     float
  26.         The total horizontal (eastward) distance traveled.
  27.         If Snell's Law fails (sin(...) out of range), returns a large number.
  28.     """
  29.     n = len(speeds) - 1  # Number of boundaries
  30.     alpha = [0.0]*(n+1)
  31.     alpha[0] = alpha0
  32.  
  33.     # Propagate angles using Snell's Law: v_i * sin(theta_i) = v_(i+1) * sin(theta_(i+1))
  34.     for i in range(n):
  35.         ratio = (speeds[i] / speeds[i+1]) * math.sin(alpha[i])
  36.         if abs(ratio) > 1:
  37.             # Invalid angle => no real solution
  38.             return 1e9
  39.         alpha[i+1] = math.asin(ratio)
  40.  
  41.     # Compute total horizontal displacement
  42.     horiz = 0.0
  43.     for i in range(n+1):
  44.         # Path length in layer i = widths[i] / cos(alpha[i])
  45.         # because alpha[i] is the angle to the normal, so the
  46.         # "normal" component of the path is widths[i].
  47.         seg_length = widths[i] / math.cos(alpha[i])
  48.  
  49.         # The angle of travel relative to the horizontal is
  50.         # (boundary_normal_angle + alpha[i]).
  51.         seg_angle = boundary_normal_angle + alpha[i]
  52.  
  53.         # Horizontal component
  54.         horiz += seg_length * math.cos(seg_angle)
  55.  
  56.     return horiz
  57.  
  58. def total_time(alpha0, speeds, widths, boundary_normal_angle):
  59.     """
  60.     Given the same setup as total_horizontal_displacement, compute
  61.     the total travel time if the initial angle is alpha0.
  62.     """
  63.     n = len(speeds) - 1
  64.     alpha = [0.0]*(n+1)
  65.     alpha[0] = alpha0
  66.  
  67.     # Again apply Snell's Law across boundaries
  68.     for i in range(n):
  69.         ratio = (speeds[i] / speeds[i+1]) * math.sin(alpha[0 + i])
  70.         if abs(ratio) > 1:
  71.             return 1e9
  72.         alpha[i+1] = math.asin(ratio)
  73.  
  74.     # Sum travel times in each layer
  75.     total_t = 0.0
  76.     for i in range(n+1):
  77.         seg_length = widths[i] / math.cos(alpha[i])
  78.         total_t += seg_length / speeds[i]
  79.  
  80.     return total_t
  81.  
  82. def solve_snell_for_100_leagues(speeds, widths, boundary_normal_angle, target=100.0):
  83.     """
  84.     Uses a simple binary search on alpha0 to achieve a desired
  85.     total horizontal displacement (default 100 leagues).
  86.     Returns (alpha0_opt, min_time).
  87.     """
  88.     def f(alpha0):
  89.         # difference from desired horizontal displacement
  90.         return total_horizontal_displacement(alpha0, speeds, widths, boundary_normal_angle) - target
  91.  
  92.     # Search range for alpha0: from 0 to something less than pi/2
  93.     # so that cos(alpha0) is positive.
  94.     low, high = 0.0, math.radians(80.0)  # somewhat arbitrary upper bound
  95.     for _ in range(100):
  96.         mid = 0.5*(low + high)
  97.         val = f(mid)
  98.         if val > 0:
  99.             # If we overshoot the horizontal distance, reduce angle
  100.             high = mid
  101.         else:
  102.             # If we undershoot, increase angle
  103.             low = mid
  104.  
  105.     alpha_opt = 0.5*(low + high)
  106.     # Compute the time corresponding to that alpha
  107.     t_opt = total_time(alpha_opt, speeds, widths, boundary_normal_angle)
  108.     return alpha_opt, t_opt
  109.  
  110. def main():
  111.     # Example speeds in each region:
  112.     #   - Region 0: normal terrain (10 leagues/day)
  113.     #   - Regions 1..5: inside the marsh, each with distinct speeds
  114.     #   - Region 6: normal terrain again (10 leagues/day)
  115.     speeds = [10, 9, 8, 7, 6, 5, 10]
  116.  
  117.     # Example widths (in leagues) in the boundary-normal direction:
  118.     # Suppose each region is 10 leagues thick. Adjust as needed.
  119.     widths = [10, 10, 10, 10, 10, 10, 10]
  120.  
  121.     # Suppose the marsh boundaries run at 45° to the horizontal,
  122.     # so the boundary normal is also at ±45°. Let's pick +45°:
  123.     boundary_normal_angle = math.radians(45.0)
  124.  
  125.     # Solve for the path that yields 100 leagues total horizontal displacement
  126.     alpha_solution, time_solution = solve_snell_for_100_leagues(
  127.         speeds, widths, boundary_normal_angle, target=100.0
  128.     )
  129.  
  130.     print(f"Optimal initial angle (relative to boundary normal) = {math.degrees(alpha_solution):.3f} degrees")
  131.     print(f"Total travel time = {time_solution:.5f} days")
  132.  
  133. if __name__ == "__main__":
  134.     main()
Success #stdin #stdout 0.98s 41684KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Optimal initial angle (relative to boundary normal) = 30.000 degrees
Total travel time = 134217737.75181 days
https://ideone.com/raX4iP
language:
Python 3 (python  3.12)
created:
10 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!