风螺旋公切线算法详解

好久不见，近来一切可好？2017年最后这段时间里，狂补了一把C#，希望未来能够从软件代码层面实现风螺旋算法的验证与推广。今天跟大家分享的这个话题的底图就是最近一段时间的学习成果：一个基于WPF架构的非常简单的绘图框架，以及对风螺旋的自动化绘制进行的实现。闲话少叙，开始今天的主题。

在掌握了风螺旋切线计算的基础上，就可以开始公切线算法的研究了。公切线的计算是飞行程序模板中非常关键的一项内容，因此，在开始模板算法分享之前，详细回顾一下公切线的算法是非常有必要的。

风螺旋公切线计算首先是基于以下4个基础设定：

/// 1. 同一程序模块内，转弯基础参数相同。无论是程序转弯、还是等待模板，在同一个ICAO标准模板中，转弯参数是相同的，即转弯半径是相同的。

/// 2. 风螺旋的切线与Esita方向相垂直。之前从原理及数学软件上对此项内容进行了证明，这是公切线的一个基础条件。

/// 3. Esita方向与对应的基准圆半径之间的夹角恒等于偏流角DA。ICAO DOC8168中的配图明确标注了此项内容DA= arcsin(w/v)。

/// 4. Esita的大小与sita角度可以互相换算。距离可以换算成旋转角度，反之旋转的角度可以换算为对应的距离。这是本公众号中多次提及的一个基本概念。风螺旋是假想的航空器的飞行轨迹，对这个轨迹整体进行旋转时，它的外扩规律符合按照Esita方向等距离外扩的特性，外扩的大小与旋转过的角度有直接的换算关系。

除以上4个基础设定以外，公切线计算中按照个人偏好使用了以下的标识来表示风螺旋的参数（如图1所示）：

图1 风螺旋参数示意图

 rotation： 从系统的零度旋转至风螺旋起点的角度。航空地图中采用的是磁北为零度，计算机系统中多数以水平向右为零度的起点，顺时针角度增大为默认条件。

offset：表示外扩距离的一个量。对于一个风螺旋来说，外扩和旋转角度可以互相换算，因此，最佳方案是对每一个风螺旋都设定一个初始的系统角度（rotation）和一个初始的外扩距离（offset）。

sita：用来表示从风螺旋起点开始实际的转弯角度。通过该角度可以换算出从起点开始一共外扩了多少距离，或者用来表示风螺旋中的某个具体位置点。

Esita：与sita角相对应的外扩距离。这个Esita的方向与对应半径的延长线之间的夹角为偏流角DA，与该位置处风螺旋切线的夹角为90度。

为了更清晰的表示计算关系，对图1中的关键点进行编号得到图2的内容：

图2  风螺旋公切线的角度关系

    图2中绘制了两条实心的风螺旋线ws1和ws2，以及与它们对应的无外扩的虚线风螺旋。虚线风螺旋与圆周的交点视为风螺旋的实际起始角度。旋转与外扩可以相互换算，因此可知，实线风螺旋与虚线风螺旋之间是“等距平行”关系。

由于Esita与风螺旋的切线相垂直（基础设定2），因此，若两个风螺旋之间存在有公切线，则必定分别有两个Esita与公切线相垂直，如图1中的Esita1和Esita2，它们二者为平行关系。

由于两个Esita与对应的半径Radius之间夹角固定为偏流角DA（基础设定3），因此可知对应半径之间为平行关系。因为，两条风螺旋的基础参数是相同的，因此，两个半径是相等的。于是图2中的abcd是一个平行四边形。从d点向fc边做垂线，得到垂足e点，可知defg是一个矩形。

ws2风螺旋在公切点位置处的Esita可以分成两部分，EsitaA和EsitaB，前者与ws1中的外扩距离是相等的，后者产生的原因是两条风螺旋初始条件的差异所造成。若每弧度可以外扩的距离为Edist，则EsitaB可表示为offset2 - offset1 + (rotation1-rotation2) * Edist，即初始外扩距离之差，加上初始角度之差所换算成的外扩距离，就构成了图2中的EsitaB。图2中的EsitaA可视为两条风螺旋共同旋转过的角度所形成的距离，因为外扩的效率是相同的，相同的外扩角度可以得到相同的外扩距离。rotation1-rotation2这个公式是浪费了两天多的时间才最终纠错纠出来的，之前在Flash系统中不存在此问题，C#中的绘图系统真是复杂许多。rotation1-rotation2代表了以rotation2为基准角度，初始的角度差值。

具体的计算步骤如下：

1、放置风螺旋，计算dc线的角度。

当我们向屏幕中放置一条风螺旋时，它的中心点位置就可视为一个已知条件。设两个中心点坐标分别为c1(x,y)，c2(x,y)，那么，线段ab在软件系统中的角度为: arctan((c2.y-c1.y)/(c2.x-c1.x))，这里需要用到反正切函数，获得一个角度值，称之为基准线角度，用baseLineAngle来表示，圆心点间的距离用distance来表示。dc线与ab线平行且相等，因此，baseLineAngle就是dc线的方向。

2、计算dg线的角度。

图2中⊿dce为一个直角三角形，∠cde的值用angle来表示，它等于arcsin(EsitaB/distance)。图2中逆时针方向角度减小，因此可将线段dg的角度表示为：baseLineAngle – angle - Pi/2 （所有角度值均用弧度来表示），公切线的方向就等于baseLineAngle – angle。

若ws1的外扩距离大于ws2时，angle角的值为负值，位置关系如图3所示：

图3 angle角为负值的情况   

此时，∠cdg的值可表示为Pi/2 + angle（因为angle小于零），线段dg的角度可以表示为baseLineAngle -（Pi/2 + angle），化简后，得到的公式为baseLineAngle – angle- Pi/2。与前面讨论的angle大于零时的公式是完全一致的，也就是说不论ws1或ws2外扩距离谁大，dg线的角度始终可以用baseLineAngle – angle - Pi/2来表示，公切线的角度可以用baseLineAngle– angle来表示。

3、计算dg线的长度。

从dg线的角度反推对应的半径的方向，可采用公式：baseLineAngle - angle - Pi/2 + DA（因为Esita与半径的夹角总等于DA）。得到圆周半径方向之后，分别与各自的rotation参数相减，就可以得到风螺旋真实的旋转角度，用这个旋转角度sita换算成外扩距离，再加上初始的外扩距离，就可以得到dg线的长度，角度、长度都知道以后，可以得到准确的切点位置。连接两个切点的坐标，得到的线段即是两条风螺旋线的公切线。

公切线的计算过程就是这样，您get到了没？

另外再附两张公切线的示意图，因为手工条件下真是很难想像公切线会出现在什么地方（以后这个问题将不是问题　：－））

    这两张图与前面的图的区别在于：初始的 rotation 大小关系对调了一下，结论完全相同。这就是自动化计算的魅力所在，哈哈哈。

        可以预见的未来中，风螺旋线将成为一门新兴的学科，它将与阿基米德螺旋、渐开线螺旋共同纳入 等距螺旋 的框架，成为几何课程的一部分，风螺旋相关计算将成为飞行程序设计入门课程的一部分。本公众号将以此为目标持续努力，全面开源共享。期待在不远的将来风螺旋线精确计算方法将成为DOC8168的一部分，并且成为飞行程序领域送给 几何研究领域的一份礼物！

        祝大家：新年快乐！