fork(1) download
copy 
  1. import math
  2.  
  3. def total_horizontal_displacement(alpha0, speeds, widths, boundary_normal_angle):
  4.     """
  5.     Given an initial angle alpha0 (radians, measured relative to the layer normal),
  6.     propagate the refraction (using Snell’s law) through the layers and compute
  7.     the total horizontal displacement.
  8.  
  9.     The layers have a width (measured along the normal) given by widths[i].  
  10.     The absolute travel direction in layer i is (boundary_normal_angle - alpha[i]).
  11.     """
  12.     n = len(speeds)
  13.     alpha = [0.0] * n
  14.     alpha[0] = alpha0
  15.     # Propagate using Snell’s law: v_i * sin(alpha_i) = v_(i+1) * sin(alpha_(i+1))
  16.     for i in range(n - 1):
  17.         ratio = (speeds[i] / speeds[i+1]) * math.sin(alpha[i])
  18.         if abs(ratio) > 1:
  19.             return 1e9  # No real solution in this branch.
  20.         alpha[i+1] = math.asin(ratio)
  21.  
  22.     horiz = 0.0
  23.     for i in range(n):
  24.         # Distance traveled in layer i
  25.         L = widths[i] / math.cos(alpha[i])
  26.         # Absolute travel angle relative to horizontal:
  27.         theta = boundary_normal_angle - alpha[i]
  28.         horiz += L * math.cos(theta)
  29.     return horiz
  30.  
  31. def total_time(alpha0, speeds, widths, boundary_normal_angle):
  32.     """
  33.     Compute total travel time given the initial angle alpha0.
  34.     """
  35.     n = len(speeds)
  36.     alpha = [0.0] * n
  37.     alpha[0] = alpha0
  38.     for i in range(n - 1):
  39.         ratio = (speeds[i] / speeds[i+1]) * math.sin(alpha[i])
  40.         if abs(ratio) > 1:
  41.             return 1e9
  42.         alpha[i+1] = math.asin(ratio)
  43.  
  44.     time = 0.0
  45.     for i in range(n):
  46.         L = widths[i] / math.cos(alpha[i])
  47.         time += L / speeds[i]
  48.     return time
  49.  
  50. def solve_snell(speeds, widths, boundary_normal_angle, target):
  51.     """
  52.     Uses binary search on alpha0 so that the total horizontal displacement equals the target.
  53.     Returns the optimal alpha0 (in radians) and the corresponding total travel time.
  54.     """
  55.     def f(alpha0):
  56.         return total_horizontal_displacement(alpha0, speeds, widths, boundary_normal_angle) - target
  57.  
  58.     low = 0.0
  59.     high = math.radians(10)  # Start with a small guess.
  60.     # Increase high until f(high) > 0 (i.e. horizontal displacement exceeds target)
  61.     while f(high) < 0:
  62.         high *= 2
  63.         if high > math.radians(89):  # safeguard
  64.             break
  65.  
  66.     # Binary search for alpha0
  67.     for _ in range(100):
  68.         mid = (low + high) / 2.0
  69.         if f(mid) > 0:
  70.             high = mid
  71.         else:
  72.             low = mid
  73.     alpha_opt = (low + high) / 2.0
  74.     t_opt = total_time(alpha_opt, speeds, widths, boundary_normal_angle)
  75.     return alpha_opt, t_opt
  76.  
  77. def main():
  78.     # Speeds in each layer:
  79.     # Layers 0 and 6 represent terrain outside the marsh (10 leagues/day),
  80.     # layers 1-5 are inside the marsh.
  81.     speeds = [10, 9, 8, 7, 6, 5, 10]
  82.  
  83.     # Each layer has width 10 leagues measured along the normal to the boundaries.
  84.     widths = [10] * 7
  85.  
  86.     # The marsh boundaries are oriented at 45° relative to horizontal.
  87.     # Thus the boundary normal is at 45°.
  88.     boundary_normal_angle = math.radians(45)
  89.  
  90.     # We require that the total horizontal displacement equals 100 leagues.
  91.     target = 100.0
  92.  
  93.     alpha_opt, t_opt = solve_snell(speeds, widths, boundary_normal_angle, target)
  94.  
  95.     print("Optimal initial angle (relative to normal) = {:.6f} degrees".format(math.degrees(alpha_opt)))
  96.     print("Total travel time = {:.13f} days".format(t_opt))
  97.  
  98. if __name__ == "__main__":
  99.     main()
Success #stdin #stdout 0.1s 14292KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Optimal initial angle (relative to normal) = 27.850759 degrees
Total travel time = 15.2880830803470 days
https://ideone.com/iLHwV9
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!