设阳光方向为 $ \hat{s}=\cos\theta\ \hat{x} - \sin\theta\ \hat{y}$ ($\hat{x}, \hat{y}$ 分别为水平方向与垂直方向的单位向量); 观察者位置为 $(X,Y)$， 设水面高度变化为 $y = h \sin (k x) $。接下来，我们需要计算的是水面的切面方向，列出光线的反射方程，也就是入射阳光方向与水面切面的角度等于反射光线与水面切面的角度。

水面切线方向是 $\vec{t} = \hat{x} + dy/dx \ \hat{y} = \hat{x} + hk \cos (kx) \hat{y}$。（这里不需要计算单位向量，因为 t 长度会在后面的点积中消去）

反射光线向量（从水面反射点到观察点的向量）是 $\vec{r} = (X-x)\ \hat{x} + (Y-y)\ \hat{y}$，单位向量是

 $\hat{r} = \frac{ (X-x)\ \hat{x} + (Y-y)\ \hat{y} }{\sqrt{ (X-x)^2 + (Y-y)^2}}$

反射方程是 入射光线向量与切向的点积等于反射线向量与切向的点积，用 bra-ket 符号表示为，$<\hat{s}|\vec{t}> =  <\hat{r} | \vec{t}>$

 也就是 

$\cos\theta - hk \sin\theta \cos kx = \frac{ (X-x) + hk \ (Y-y)\ \cos kx }{\sqrt{ (X-x)^2 + (Y-y)^2}}$

$ hk (\sin\theta + \frac{ (Y-y)}{\sqrt{ (X-x)^2 + (Y-y)^2}}) \cos kx = \cos \theta - \frac{ (X-x) }{\sqrt{ (X-x)^2 + (Y-y)^2}}$

上面的方程中 $y = h \sin kx $，是随着 x 变化的，这使方程变得相当复杂。假设观察者高度 Y 远远大于水面波动幅度 h，我们可以将 Y-y 近似为 Y, 因此

$hk (\sin\theta + \frac{ Y}{\sqrt{ (X-x)^2 + Y^2}}) \cos kx = \cos \theta - \frac{ (X-x) }{\sqrt{ (X-x)^2 + Y^2}}$

先进行一下 sanity check，如果水面完全平整，也就是  k =0 , 上面的方程 reduce 成 $\cos \theta = \frac{ (X-x) }{\sqrt{ (X-x)^2 + Y^2}}$，这是对的。

最终，我们的反射方程是（如果我没有算错）：

$hk \cos kx = \frac{\sqrt{ (X-x)^2 + Y^2} \cos \theta - (X-x)}{\sqrt{ (X-x)^2 + Y^2}\sin \theta+Y}$

左边是余弦函数，右边是一个与 x 有关的数，存在多个 x 解(也就是多个日影）是可能的。进一步分析（如计算不同距离的日影大小）估计就比较 tedious了，用计算机数值解可能更有效。出于对这个问题的兴趣，我进行了数值计算，设水波动幅度为2厘米，X=3000m, Y = 200m，阳光入射角度 2/30，水波长 30米，结果得到日影位置序列：Found Sun image at: [83.53923367373932, 96.29285813823678, 113.96389798691882, 174.88654492547624, 184.95635957504467, 214.44386039408593, 235.97349038834463, 243.87496449474165, 302.407948894744, 328.6930491317427] 

Number of images:  10

PYTHON CODE如下，鉴于是数字计算，就没有做 Y-y ~ Y 的近似。有兴趣的可以自己调整参数。我试了一下，如果水波幅度 只有1厘米，则日影数量只有两个。