CN109491639A - 一种光路追踪计算方法及*** - Google Patents

一种光路追踪计算方法及*** Download PDF

Info

Publication number
CN109491639A
CN109491639A CN201811365311.2A CN201811365311A CN109491639A CN 109491639 A CN109491639 A CN 109491639A CN 201811365311 A CN201811365311 A CN 201811365311A CN 109491639 A CN109491639 A CN 109491639A
Authority
CN
China
Prior art keywords
equation
simple cubic
cubic equation
root
laser
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201811365311.2A
Other languages
English (en)
Other versions
CN109491639B (zh
Inventor
任小广
周浩
肖震
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National Defense Technology Innovation Institute PLA Academy of Military Science
Original Assignee
National Defense Technology Innovation Institute PLA Academy of Military Science
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National Defense Technology Innovation Institute PLA Academy of Military Science filed Critical National Defense Technology Innovation Institute PLA Academy of Military Science
Priority to CN201811365311.2A priority Critical patent/CN109491639B/zh
Publication of CN109491639A publication Critical patent/CN109491639A/zh
Application granted granted Critical
Publication of CN109491639B publication Critical patent/CN109491639B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F7/00—Methods or arrangements for processing data by operating upon the order or content of the data handled
    • G06F7/38—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation
    • G06F7/48—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation using non-contact-making devices, e.g. tube, solid state device; using unspecified devices
    • G06F7/483—Computations with numbers represented by a non-linear combination of denominational numbers, e.g. rational numbers, logarithmic number system or floating-point numbers
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F7/00—Methods or arrangements for processing data by operating upon the order or content of the data handled
    • G06F7/38—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation
    • G06F7/48—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation using non-contact-making devices, e.g. tube, solid state device; using unspecified devices
    • G06F7/544—Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation using non-contact-making devices, e.g. tube, solid state device; using unspecified devices for evaluating functions by calculation
    • G06F7/552—Powers or roots, e.g. Pythagorean sums
    • G06F7/5525—Roots or inverse roots of single operands

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Nonlinear Science (AREA)
  • Optical Modulation, Optical Deflection, Nonlinear Optics, Optical Demodulation, Optical Logic Elements (AREA)

Abstract

本发明提供一种光路追踪计算方法及***，包括：采集被测激光的参数值；基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；基于所述激光传播时间计算所述激光的轨迹点。本发明提供的技术方案根据浮点稳定性计算理论，采用高精度计算方法求解光线参数方程，提高了光路追踪的准确性。

Description

一种光路追踪计算方法及***
技术领域
本发明涉及科学研究计算中的高性能计算领域，具体涉及一种光路追踪计算方法。
背景技术
激光在等离子体中的能量沉积问题是激光间接驱动惯性约束聚变(ICF，InertialConfinement Fusion)中的重要物理过程，其模拟精度会影响ICF全过程数值模拟结果的正确性。在ICF整体数值模拟中，目前激光主要采用几何光路建模，在模拟中需要通过求解光线与网格的相交关系来确定光线的运动轨迹。
光路追踪算法有广泛的应用，例如无限通信***、晶体光学、电磁辐射传播及碰撞检测、激光的传播等。在大规模并行机上实现快速高精度光路追踪计算一直是一件有挑战性的任务。计算光路与一系列(3-D或2-D)几何体的交点是光路追踪算法的核心组成部分。光的传播轨迹是由光线与几何体的交点组成的。相交点处及附近区域的物质特性决定了光线的状态和传播方向。如何提高算法的精度和效率往往因应用问题而异。在激光约束聚变问题中，激光作为聚变的源项出现在物理建模之中。具体地，通常采用几何光线近似来模拟激光在等离子体中吸收与传播。激光束被表示成多条激光光线(直线或曲线)。每条光线被赋予一定的激光能量。随着激光的传播，其携带的能量沉积到沿途的等离子体中。在每个计算网格的电子数密度线性分布的假设下，求交算法的核心就是求出光线参数方程与每个网格边线段的交点。在算法的设计上，一方面需要考虑到电子数密度为常数等可能的退化情形，另一方面需要考虑到计算机数字的有限表示位数，理论上相同数值上不同的算法可能会得到不同的数值解。由于浮点计算的局限性，很多常用的算法难以计算出理论上的精确解。而且在计算交点的算法中当判据或阈值达到计算机浮点精度的量级时，数值结果还可能会因计算机不同而不同，数值结果不同可能会引起光路传播方向的改变，或者传播路径长度的改变。
IEEE二进制浮点数算术标准(IEEE 754)是20世纪80年代以来最广泛使用的浮点数运算标准，为许多CPU与浮点运算器所采用，这个标准定义了表示浮点数的格式与反常值，一些特殊数值，例如：无穷(Inf)与非数值(NaN)，以及这些数值的“浮点数运算符”；它也指明了四种数值舍入规则和五种例外状况，包括例外发生的时机与处理方式。
因为IEEE 754标准指明了浮点数的有效存储位数，导致浮点数在计算机中运算结果的不精确性来自于浮点数的有限位数无法精确表达自然界中的一个自然数，进一步导致了一些算法在理论上正确无误，但在计算机上则可能计算结果与真实数值不同，而计算错误。
发明内容
为解决上述技术的不足，本发明提供一种光路追踪计算方法，根据浮点稳定性计算理论，采用扩展精度浮点数计算方法求解光线参数方程，解决了减法相消带来的计算错误，提高了光线参数方程求解的准确性，能精准地实现光路的追踪。
一种光路追踪计算方法，其改进之处在于，所述方法包括：
采集被测激光的参数值；
基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；
基于所述激光传播时间计算所述激光的轨迹点；
所述光线参数方程包括一元四次方程。
优选地，所述基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间包括：
将一元四次光线参数方程简化为一元三次方程；
基于扩展精度浮点数采用三角函数法求解所述一元三次方程得到所述一元三次方程的根；
基于所述一元三次方程的根，计算一元四次光线参数方程的解，得到激光传播时间；
其中，所述扩展精度浮点数包括精度高于双精度的浮点数。
优选地，所述基于扩展精度浮点数采用三角函数法求解所述一元三次方程包括：
计算所述一元三次方程的判别式；
基于所述判别式结果求解所述一元三次方程，得到求解结果；基于所述求解结果确定所述一元三次方程的根。
优选地，所述基于所述判别式结果求解所述一元三次方程，得到求解结果包括：
当所述一元三次方程的判别式≥0时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解；
否则，判断求根公式中是否出现减法相消，当出现减法相消时，采用扩展精度浮点数对所述一元三次方程的求根公式求解；当不出现减法相消时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解。
优选地，所述基于所述求解结果确定所述一元三次方程的根包括：
采用牛顿迭代法对所述求解结果进行迭代计算；
当所述一元三次方程的判别式＜0时，检测迭代计算后的根是否有重根；若存在重根，对其中一个重根采用牛顿迭代法再次进行迭代计算，直至不出现重根为止；
将所述迭代计算结果作为所述一元三次方程的根。
优选地，所述是否出现减法相消按下式进行判断：
式中，m：求根公式中第一浮点数；n：求根公式中第二浮点数，m和n的大小非常接近；ε：检测减法相消的阈值。
优选地，所述一元三次方程的判别式如下式所示：
式中：Δ：判别式；p、q：消去二次项后的一元三次方程的系数；
其中：p按下式进行计算：
q按下式进行计算：
优选地，所述一元三次方程的判别式＜0时，所述一元三次方程的根的个数为3个，求根公式如下式所示：
式中：x1、x2、x3：一元三次方程的根；θ：采用三角函数法求解一元三次方程时的参数；
其中，θ按下式进行计算：
优选地，所述基于所述激光传播时间计算激光的轨迹点包括：
将所述激光传播时间代入到关于光线传播时间的光线参数方程，得到激光的轨迹点。
一种光路追踪计算***，包括参数采集模块、时间获取模块和轨迹点获取模块；
参数采集模块：用于采集被测激光的参数值；
时间获取模块：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；
轨迹点获取模块：基于所述激光传播时间计算所述激光的轨迹点；
所述光线参数方程包括一元四次方程。
与最接近的已有技术比，本发明提供的技术方案具有以下有益效果：
本发明提供的技术方案，根据浮点稳定性计算理论，采用扩展精度浮点数计算方法求解光线参数方程，解决了减法相消带来的计算错误，提高了光线参数方程求解的准确性，精准地实现光路的追踪。
本发明提供的技术方案，采用扩展精度运算、迭代提炼精度和重根检测的方法优化解的精度，提高了光线参数方程求解算法的稳定性。
附图说明
图1是本发明光路追踪计算方法的示意图；
图2是本发明基于扩展精度浮点数采用三角函数法求解一元三次方程的流程图；
图3是本发明光路追踪计算***的示意图。
具体实施方式
为了更好地理解本发明，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。
实施例一、
一种光路追踪计算方法，如图1所示，包括：
步骤1：采集被测激光的参数值；
步骤2：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；
步骤3：基于所述激光传播时间计算所述激光的轨迹点；
所述光线参数方程包括一元四次方程。
步骤1：采集被测激光的参数值。
步骤2：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间。
具体地，所述采用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间包括：
将一元四次光线参数方程简化为一元三次方程；
基于扩展精度浮点数采用三角函数法求解所述一元三次方程得到所述一元三次方程的根；
基于所述一元三次方程的根，计算一元四次光线参数方程的解，得到激光传播时间；
其中，所述扩展精度浮点数包括精度高于双精度的浮点数。
具体地，所述基于扩展精度浮点数采用三角函数法求解所述一元三次方程包括：
计算所述一元三次方程的判别式；
基于所述判别式结果求解所述一元三次方程，得到求解结果；基于所述求解结果确定所述一元三次方程的根。
具体地，所述基于所述判别式结果求解所述一元三次方程，得到求解结果包括：
当所述一元三次方程的判别式≥0时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解；
否则，判断求根公式中是否出现减法相消，当出现减法相消时，采用扩展精度浮点数对所述一元三次方程的求根公式求解；当不出现减法相消时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解。
具体地，所述基于所述求解结果确定所述一元三次方程的根包括：
采用牛顿迭代法对所述求解结果进行迭代计算；
当所述一元三次方程的判别式＜0时，检测迭代计算后的根是否有重根；若存在重根，对其中一个重根采用牛顿迭代法再次进行迭代计算，直至不出现重根为止；
将所述迭代计算结果作为所述一元三次方程的根。
具体地，所述减法相消按下式进行判断：
式中，m：求根公式中第一浮点数；n：求根公式中第二浮点数，m和n的大小非常接近；ε：检测减法相消的阈值。
具体地，所述一元三次方程的判别式如下式所示：
式中：Δ：判别式；p、q：消去二次项后的一元三次方程的系数；
其中：p按下式进行计算：
q按下式进行计算：
具体地，所述一元三次方程的判别式＜0时，所述一元三次方程的根的个数为3个，求根公式如下式所示：
式中：x1、x2、x3：一元三次方程的根；θ：采用三角函数法求解一元三次方程时的参数；
其中，θ按下式进行计算：
步骤3：所述基于所述激光传播时间计算激光的轨迹点。
将所述激光传播时间代入到关于光线传播时间的光线参数方程，得到激光的轨迹点。
实施例二、
一种光路追踪计算方法，如图1所示，所述方法包括：
步骤1：采集被测激光的参数值；
步骤2：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；
步骤3：基于所述激光传播时间计算所述激光的轨迹点；
所述光线参数方程包括一元四次方程。
步骤1：采集被测激光的参数值。
在几何光学近似下，光线传播方程可以描述为
其中是光线的位置矢量，ω和分别是激光的频率和传播速度
t表示激光传播的时间，c是真空中的光速，ne是电子数密度，临界电子数密度nc＝ω2me/4πe2，me是电子质量，e为电子电量。
描述激光在等离子体中传播的是麦克斯韦方程，它的求解要求网格尺度是激光波长量级的，在时间尺度方面也有同样的要求。而目前的模拟程序还远远不能达到。在这种情况下提出了几何光学近似下的光路追踪法。它的基本思想是：将激光束分成若干条光线，每条光线初始时赋予一定的激光能量，在等离子体中传播，其行为满足光线参数方程，沿着光线轨迹，激光能量通过逆韧致过程和其他过程被等离子体吸收。由于激光传播的速度很快，可以假设在一个流体力学步长内能量沉积是瞬时的，在一条光线传播过程中，等离子体的状态是不变的。
为将三维光线应用于以z为对称轴的二维柱几何数值计算中，即
其中(kr，kθ，kz)表示光线在柱坐标系下的速度，可通过坐标变换将直角坐标系中的光线参数方程变换为柱坐标系下的光线参数方程：
通过无量纲化略去关于网格的附标，并且假定在每一个网格上电子数密度的变化是线性的，则可以得到柱坐标系中每个二维网格上的三维光线参数方程：
其中，(r0，z0，θ0)与分别为光线进入点的位置与波矢。求解上述方程可以得到关于t的光线参数方程，具体地，
a)当h*＝0时，光线参数方程是：
b)当Gr＝0时，光线参数方程是：
c)当h*≠0且Gr≠0时，无法精确求解，记令s＝y′(r)，
则光线参数方程可近似表示为：
基于以上光线参数方程的具体形式，可以计算光线参数方程与二维柱几何网格的交点。
具体地，要分别计算光线与每条网格边的交点。把网格边所在的直线方程记作
r＝Az+B，A＝(r2-r1)/(z2-z1)，B＝r1-Az1，
其中(r1，z1)和(r2，z2)分别为任一网格边的两个端点，由此可以相应地得到如下关于t的一元二次方程和一元四次方程。
a)当h*＝0时，
b)当Gr＝0时，
c)当h*≠0且Gr≠0时，
求解出上述一元二次方程或一元四次方程就可以得到参数t的值，代入光线参数方程即可得到交点的坐标值。求解一元二次或四次方程可能会得到多个根，判断真解的条件是t＞10-6且交点在网格边上。如有多个根满足真解条件，应选取最小值作为光线的轨迹点。根据光线的轨迹点就可以求出光线穿过网格的路程，继而计算出这条光线沉积的激光能量。
费拉里算法是求解一元四次方程的一种常用算法，具体如下：
对于一元四次方程
a′x4+b′x3+c′x2+d′x+e＝0
将方程的最高阶次项化为1，执行
将方程移项，化为
x4+b′x3＝-c′x2-d′x-e
方程两边同时加上使左边配成完全平方式，如下所示：
引入参数y并解出y，使得方程左右两边都配成完全平方式。
首先，左右都加上得到
需要将右边配成完全平方式，即右边的判别式为0，如下式所示：
化简为一个三次方程
-y3+c′y2+(4e-b′d′)y+d′2+eb′2-4ec′＝0
然后求解这个一元三次方程，会得到y的三个值，取其中任何一个即可。可以得到原方程等价于
可以得到两个一元二次方程
对这两个一元二次方程化简
求解这两个一元二次方程，即可得到一元四次方程的四个实根。
在上述费拉里算法中，需要求解一个一元三次方程，该方程的解是否精确将直接影响费拉里算法解出的一元四次方程根的正确性。
步骤2：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间。
具体地，所述采用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间包括：
将一元四次光线参数方程简化为一元三次方程；
基于扩展精度浮点数采用三角函数法求解所述一元三次方程得到所述一元三次方程的根；
基于所述一元三次方程的根，计算一元四次光线参数方程的解，得到激光传播时间；
其中，所述扩展精度浮点数包括精度高于双精度的浮点数。
计算机中常用的浮点数是单精度，双精度浮点数，本发明使用了扩展精度浮点数，这些类型浮点数的指数和有效数位数如下：
表1 浮点数的位数分布
类型 占用位数 指数位数 小数位数
单精度 32 8 23
双精度 64 11 52
如图2所示，所述基于扩展精度浮点数采用三角函数法求解所述一元三次方程包括：
计算所述一元三次方程的判别式；
基于所述判别式结果求解所述一元三次方程，得到求解结果；基于所述求解结果确定所述一元三次方程的根。
具体地，所述基于所述判别式结果求解所述一元三次方程，得到求解结果包括：
当所述一元三次方程的判别式≥0时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解；
否则，判断求根公式中是否出现减法相消，当出现减法相消时，采用扩展精度浮点数对所述一元三次方程的求根公式求解；当不出现减法相消时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解。
具体地，所述基于所述求解结果确定所述一元三次方程的根包括：
采用牛顿迭代法对所述求解结果进行迭代计算；
当所述一元三次方程的判别式＜0时，检测迭代计算后的根是否有重根；若存在重根，对其中一个重根采用牛顿迭代法再次进行迭代计算，直至不出现重根为止；
将所述迭代计算结果作为所述一元三次方程的根。
牛顿迭代法求解方程f(x)＝0。给定初始估计解x0，牛顿迭代法是通过以下公式迭代求得下一步解，如下式所示：
牛顿迭代法能够收敛到解的精度的极限取决于能计算到的精确程度。
在浮点计算中，减法相消发生在两个数值上非常接近的浮点数相减，这经常会导致大量的有效数字丢失。考虑一个函数
其值域为[0，0.5]，在x＝1.2×10-5，有效数字为10位的情况下
cosx＝0.9999999999
1-cosx＝0.0000000001
这明显计算错误，这个问题的原因是1-cosx的有效位数只有一位，1-cosx的浮点计算结果是精确的，但减法操作把浮点误差放大了。在这个函数中，利用恒等式cosx＝1-2sin2(x/2)，将方程重写为以下形式从而避免舍入误差
假设有两个相近的浮点数m和n，有和 和是m和n欲代表自然数，Δm和Δn是数据中的舍入误差，或者之前计算的误差，计算x＝m-n，有
当|m-n|＜＜|m|+|n|时，结果和自然界中对应的结果会产生严重的相对误差。
但是，当Δm＝Δn＝0，即浮点数能精确表示自然数时，不会产生减法相消误差。所述减法相消按下式进行判断：
式中，m：求根公式中第一浮点数；n：求根公式中第二浮点数，m和n的大小非常接近；ε：检测减法相消的阈值。
所述一元三次方程如下式所示：
ax3+bx2+cx+d＝0。
具体地，所述一元三次方程的根的个数为3个时，求根公式如下式所示：
式中，p、q：消去二次项后的一元三次方程的系数；
此时，判别式小于零。
其中，采用三角函数法求解一元三次方程时的参数θ按下式进行计算：
p按下式进行计算：
判别式按下式进行计算：
其中：Δ：判别式；
q按下式进行计算：
若判别式等于零：
若p＝q＝0，一元三次方程有一个根，求根公式如下式所示：
否则一元三次方程有两个根：
若q＜0，求根公式如下式所示：
否则求根公式如下式所示：
若判别式＞0，一元三次方程有一个根：
若求根公式如下式所示：
否则求根公式如下式所示：
步骤3：将所述激光传播时间代入光线参数方程，得到光线的轨迹点包括：将所述激光传播时间代入到关于光线传播时间的光线参数方程，得到激光的轨迹点。
表2列出了在光路追踪高精度计算方法中使用判别式扩展精度、求根公式减法相消检测与扩展精度、牛顿迭代法、重根检测与重启牛顿迭代的各种情况，以及各种情况下经过对1300万个一元三次方程进行计算，统计的求解正确率。由此可见使用本发明提供的技术方案可获得100％正确率。而经典三角函数算法，没有使用判别式扩展精度、求根公式减法相消检测与扩展精度、牛顿迭代法、重根检测与重启牛顿迭代的情况下对这些方程的正确率约为85％。
表2 比较表
实施例三、
一种光路追踪计算***，如图3所示，包括参数采集模块、时间获取模块和轨迹点获取模块；
参数采集模块：用于采集被测激光的参数值；
时间获取模块：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；
轨迹点获取模块：基于所述激光传播时间计算所述激光的轨迹点；
所述光线参数方程包括一元四次方程；
具体地，轨迹点获取模块中基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间包括：
将一元四次光线参数方程简化为一元三次方程；
基于扩展精度浮点数采用三角函数法求解所述一元三次方程得到所述一元三次方程的根；
基于所述一元三次方程的根，计算一元四次光线参数方程的解，得到激光传播时间；
其中，所述扩展精度浮点数包括精度高于双精度的浮点数。
具体地，所述基于扩展精度浮点数采用三角函数法求解所述一元三次方程：
计算所述一元三次方程的判别式；
基于所述判别式结果求解所述一元三次方程，得到求解结果；基于所述求解结果确定所述一元三次方程的根。
具体地，所述基于所述判别式结果求解所述一元三次方程，得到求解结果包括：
当所述一元三次方程的判别式≥0时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解；
否则，判断求根公式中是否出现减法相消，当出现减法相消时，采用扩展精度浮点数对所述一元三次方程的求根公式求解；当不出现减法相消时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解。
具体地，所述基于所述求解结果确定所述一元三次方程的根包括：
采用牛顿迭代法对所述求解结果进行迭代计算；
当所述一元三次方程的判别式＜0时，检测迭代计算后的根是否有重根；若存在重根，对其中一个重根采用牛顿迭代法再次进行迭代计算，直至不出现重根为止；
将所述迭代计算结果作为所述一元三次方程的根。
具体地，所述减法相消按下式进行判断：
式中，m：求根公式中第一浮点数；n：求根公式中第二浮点数，m和n的大小非常接近；ε：检测减法相消的阈值。
具体地，所述一元三次方程的判别式如下式所示：
式中：Δ：判别式；p、q：消去二次项后的一元三次方程的系数；
其中：p按下式进行计算：
q按下式进行计算：
具体地，所述一元三次方程的判别式＜0时，所述一元三次方程的根的个数为3个，求根公式如下式所示：
式中：x1、x2、x3：一元三次方程的根；θ：采用三角函数法求解一元三次方程时的参数；
其中，采用三角函数法求解一元三次方程时的参数θ按下式进行计算：
本领域内的技术人员应明白，本申请的实施例可提供为方法、***、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(***)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员依然可以对本发明的具体实施方式进行修改或者等同替换，这些未脱离本发明精神和范围的任何修改或者等同替换，均在申请待批的本发明的权利要求保护范围之内。

Claims (10)

1.一种光路追踪计算方法，其特征在于，所述方法包括：
采集被测激光的参数值；
基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；
基于所述激光传播时间计算所述激光的轨迹点；
所述光线参数方程包括一元四次方程。
2.如权利要求1所述的计算方法，其特征在于，所述基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间包括：
将一元四次光线参数方程简化为一元三次方程；
基于扩展精度浮点数采用三角函数法求解所述一元三次方程得到所述一元三次方程的根；
基于所述一元三次方程的根，计算一元四次光线参数方程的解，得到激光传播时间；
其中，所述扩展精度浮点数包括精度高于双精度的浮点数。
3.如权利要求2所述的计算方法，其特征在于，所述基于扩展精度浮点数采用三角函数法求解所述一元三次方程包括：
计算所述一元三次方程的判别式；
基于所述判别式结果求解所述一元三次方程，得到求解结果；基于所述求解结果确定所述一元三次方程的根。
4.如权利要求3所述的计算方法，其特征在于，所述基于所述判别式结果求解所述一元三次方程，得到求解结果包括：
当所述一元三次方程的判别式≥0时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解；
否则，判断求根公式中是否出现减法相消，当出现减法相消时，采用扩展精度浮点数对所述一元三次方程的求根公式求解；当不出现减法相消时，采用单精度或双精度浮点数对所述一元三次方程求根公式进行求解。
5.如权利要求3所述的计算方法，其特征在于，所述基于所述求解结果确定所述一元三次方程的根包括：
采用牛顿迭代法对所述求解结果进行迭代计算；
当所述一元三次方程的判别式＜0时，检测迭代计算后的根是否有重根；若存在重根，对其中一个重根采用牛顿迭代法再次进行迭代计算，直至不出现重根为止；
将所述迭代计算结果作为所述一元三次方程的根。
6.如权利要求4所述的计算方法，其特征在于，所述求根公式中是否出现减法相消按下式进行判断：
式中，m：求根公式中第一浮点数；n：求根公式中第二浮点数，m和n的大小非常接近；ε：检测减法相消的阈值。
7.如权利要求2所述的计算方法，其特征在于：所述一元三次方程的判别式如下式所示：
式中：Δ：判别式；p、q：消去二次项后的一元三次方程的系数；
其中：p按下式进行计算：
q按下式进行计算：
8.如权利要求7所述的计算方法，其特征在于，所述一元三次方程的判别式＜0时，所述一元三次方程的根的个数为3个，求根公式如下式所示：
式中：x1、x2、x3：一元三次方程的根；θ：采用三角函数法求解一元三次方程时的参数；
其中，采用三角函数法求解一元三次方程时的参数θ按下式进行计算：
9.如权利要求1所述的计算方法，其特征在于，所述基于所述激光传播时间计算激光的轨迹点包括：
将所述激光传播时间代入到关于光线传播时间的光线参数方程，得到激光的轨迹点。
10.一种光路追踪计算***，其特征在于，所述***包括参数采集模块、时间获取模块和轨迹点获取模块；
参数采集模块：用于采集被测激光的参数值；
时间获取模块：基于所述采集的参数值，利用高精度计算方法求解预先构建的光线参数方程，计算激光传播时间；
轨迹点获取模块：基于所述激光传播时间计算所述激光的轨迹点；
所述光线参数方程包括一元四次方程。
CN201811365311.2A 2018-11-16 2018-11-16 一种光路追踪计算方法及*** Active CN109491639B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811365311.2A CN109491639B (zh) 2018-11-16 2018-11-16 一种光路追踪计算方法及***

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811365311.2A CN109491639B (zh) 2018-11-16 2018-11-16 一种光路追踪计算方法及***

Publications (2)

Publication Number Publication Date
CN109491639A true CN109491639A (zh) 2019-03-19
CN109491639B CN109491639B (zh) 2019-08-13

Family

ID=65695987

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811365311.2A Active CN109491639B (zh) 2018-11-16 2018-11-16 一种光路追踪计算方法及***

Country Status (1)

Country Link
CN (1) CN109491639B (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN203414169U (zh) * 2013-07-26 2014-01-29 华东师范大学 一种光路装置
CN104238025A (zh) * 2013-06-21 2014-12-24 华为技术有限公司 光路处理方法和装置
CN206584041U (zh) * 2016-11-14 2017-10-24 深圳市镭神智能***有限公司 一种光路***
CN107727084A (zh) * 2017-09-18 2018-02-23 北京卫星环境工程研究所 机器人高精度测量仪对星上立方镜法线方向自动搜索方法
US20180188398A1 (en) * 2016-12-30 2018-07-05 Institute Of Geology And Geophysics, Chinese Academy Of Sciences Gaussian Beam Migration Imaging Method and Device
CN108445469A (zh) * 2018-05-11 2018-08-24 天津大学 一种多线激光雷达的转镜扫描装置及方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104238025A (zh) * 2013-06-21 2014-12-24 华为技术有限公司 光路处理方法和装置
CN203414169U (zh) * 2013-07-26 2014-01-29 华东师范大学 一种光路装置
CN206584041U (zh) * 2016-11-14 2017-10-24 深圳市镭神智能***有限公司 一种光路***
US20180188398A1 (en) * 2016-12-30 2018-07-05 Institute Of Geology And Geophysics, Chinese Academy Of Sciences Gaussian Beam Migration Imaging Method and Device
CN107727084A (zh) * 2017-09-18 2018-02-23 北京卫星环境工程研究所 机器人高精度测量仪对星上立方镜法线方向自动搜索方法
CN108445469A (zh) * 2018-05-11 2018-08-24 天津大学 一种多线激光雷达的转镜扫描装置及方法

Also Published As

Publication number Publication date
CN109491639B (zh) 2019-08-13

Similar Documents

Publication Publication Date Title
Moon et al. Exact charge-conserving scatter–gather algorithm for particle-in-cell simulations on unstructured grids: A geometric perspective
Larsson et al. Stable computation of differentiation matrices and scattered node stencils based on Gaussian radial basis functions
Beilina et al. Numerical linear algebra: Theory and applications
Manocha et al. Algorithms for intersecting parametric and algebraic curves II: Multiple intersections
Brönnimann et al. Computing exact geometric predicates using modular arithmetic with single precision
Verschelde et al. Mixed-volume computation by dynamic lifting applied to polynomial system solving
Teukolsky Formulation of discontinuous Galerkin methods for relativistic astrophysics
Berberich et al. An exact, complete and efficient implementation for computing planar maps of quadric intersection curves
Garrett et al. Optimization and large scale computation of an entropy-based moment closure
Dassios On solutions and algebraic duality of generalised linear discrete time systems
CN109739470A (zh) 一种基于2型双曲cordic任意指数函数的计算***
Lokutsievskiy et al. Liouville integrability of sub-Riemannian problems on Carnot groups of step 4 or greater
Lulli et al. Out-of-equilibrium finite-size method for critical behavior analyses
Cano et al. A new method for numerical integration of singular functions on the plane
Frenod et al. An exponential integrator for a highly oscillatory Vlasov equation
CN109491639B (zh) 一种光路追踪计算方法及***
Duan et al. An efficient ADER discontinuous Galerkin scheme for directly solving Hamilton-Jacobi equation
CN109508172B (zh) 一种光线轨迹点计算方法及***
Todhunter et al. Mathematical approach to the validation of field surface texture parameter software
Han et al. Two uniform tailored finite point schemes for the two dimensional discrete ordinates transport equations with boundary and interface layers
Garrett et al. A fast solver for implicit integration of the Vlasov--Poisson system in the Eulerian framework
Cheng et al. Determining the topology of real algebraic surfaces
Luo A finite volume method based on weno reconstruction for compressible flows on hybrid grids
Li et al. An efficient surface–surface intersection algorithm based on geometry characteristics
Kytmanov An algorithm for calculating power sums of roots for a class of systems of nonlinear equations

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant