最近开始研究《mathematica与大学物理计算》，学习光学时遇到一些困难。
光线传播时有一个方程
方程中s是光路的弧长，r是一个矢量，表示光路上的一个位置，n是折射率，它是位置的函数，意思就是光路上任意一点处的折射率与位置矢量对路径弧长的微分的乘积对弧长的微分等于折射率的梯度，写成分量式就是两个或三个微分方程，用NDSolve解出来就能画图了，书上给出了一个例子，下面是代码
sm = 2500; n[y_] := 2 - 0.001 y ;
equ = {D[n[y[s]] y'[s], s] == -0.001,
D[n[y[s]] x'[s], s] == 0,
x[0] == 0, x'[0] == Sqrt[2]/2,
y[0] == 0, y'[0] == Sqrt[2]/2};
sol = NDSolve[equ, {x, y}, {s, 0, sm}];
ParametricPlot[{x[s], y[s]} /. sol[[1]], {s, 0, sm}]
其中n选择了一个随高度增加而减小的函数，模拟的是光在大气中的折射。
以上是背景
现在问题来了，如果折射率不连续，比如光从空气射入水中，这时如何方便地获得光路？简单地在折射率中使用UnitStep好像不行，书中使用的方法比较繁琐，而且对不同形状的界面都要进行推导，每次都要写几十行，计算光线斜率的表达式又很长，一堆方括号，很容易出错，所以希望有一个简单通用的办法。

解决了，非常感谢

用zemax模拟更方便吧？

refWave[{x1_, y1_}, {m1_, m2_}, h_, d_, f_,
x_] := {x,
y1 + (x -
x1) (m1 + (m2 -
m1) With[{t = Clip[Rescale[x - x1, {-h, h}], {0, 1}]},
t^2 (3 - 2 t)])} +
d Sin[2 \[Pi] f x] Normalize[{-Piecewise[{{m1, x - x1 < -h}, {m2,
h < x - x1}}, (m1 + m2)/
2 + (m2 - m1) (x - x1) (3/2 - ((x - x1)/h)^2)/h], 1}]
ParametricPlot[
Table[refWave[
With[{u = 1/2 + t/10}, {-Sqrt[3]/2 (u - 1), (3 u - 1)/2}], {1/10 +
t/20, 1/3 - 3 t/2}, 1/20, 1/50, 8/(5/4 - t), x], {t, 7/10,
1/10, -1/10}] // Evaluate, {x, -1/2, 3/2},
AspectRatio -> Automatic, Axes -> None, Background -> Black,
Epilog -> {White,
Polygon[{{-17/(20 Sqrt[3]), 3/20}, {-6 Sqrt[3]/25,
7/25}, {-3/2, (12 Sqrt[3] - 5)/250}, {-3/2, (17 Sqrt[3] - 45)/
300}}]}, PlotRange -> {{-3/2, 3/2}, {-3/4, 5/4}},
PlotStyle -> (RGBColor /@ {"
#9400D3", "#
4B0082", "#0000FF",
"#00FF00", "
#FFFF00", "#
FF7F00", "#FF0000"}),
Prolog -> {Directive[FaceForm[GrayLevel[1/10]],
EdgeForm[Directive[Gray, Thick]]], RegularPolygon[3]}]