装载船跟踪

加载astKonTikiData。垫文件。对于这个例子，它包含了船的轨迹，速度和航向。这个文件分别在变量“lat”和“long”中存储不同轨迹点的纬度和经度。这些变量包含了从卡亚奥港到拉罗亚环礁每天一个航迹点的足够数据。此外，该文件还存储了每天的船舶速度值，单位为“V”，航向值为“T”。

加载阿斯特康提基达塔

创建一个观察结构

航海还原过程是一系列的步骤，导航下面，以确定他的船的经度和纬度。它是基于在美国实用导航[1]，航海年鉴[2]中所述的理论，和解释性补编天文年历[3]。该过程使用来自六分仪，时钟，指南针，和航海图获得的观测数据。它返回每个所观察到的物体的截距（p）和真方位角（Z）。本例使用的观察结构阵列，观测值，以包含观察数据。用于结构阵列的字段为：

为简单起见，假设所有的测量是在在船上的相同位置，具有相同六分仪，在相同的环境温度和压力：

obs.h = 4;obs.IC = 0;obs.P = 982;obs.T = 15;

探险队于1947年4月28日出发。因此，初始化该日期的观测结构:

obs.year = 1947;obs.month = 4;obs.day = 28;

初始化航迹推算过程

为了开始航迹推算过程，定义船舶位置的初始条件。将纬度和经度的固定解决方案的纬度分别存储在latFix和longFix变量中。在本例中，对于第一个修复位置，使用秘鲁的卡亚奥的纬度和经度:

longFix =零（大小（长））;latFix =零（大小（LAT））;longFix（1）=长（1）;latFix（1）= LAT（1）;

每日航迹推算

对于这个例子，假定一个修复获得日常使用观测数据。因此，本例中使用“for循环”为每个观测。变量“m”个作为表示由于从端口出发每经过一天的计数器：

为了m=1：尺寸（横向，1）-1

增加6月和4月的天数并进行天数调整，这两个月只有30天：

obs.day=obs.day+1；[obs.month，obs.day]=ashelperdaycheck（obs.year，obs.month，obs.day）；

实际纬度和经度

从之前装载的航迹点提取出船舶每天的实际位置。这个例子使用这个值来计算当地时区和所选行星在天空中的位置:

longActual=long（m+1）；latActual=lat（m+1）；

行星选择

在给定的纬度和经度下，选择可以观测的行星。以下代码使用预计算数据:

如果longActual> -90 obs.object = {“土星”；“海王星”};elseiflongActual <=  -  90 && longActual> -95 obs.object = {“土星”；“海王星”；“木星”};elseiflongActual < = -95奥林匹克广播服务公司。对象= {“海王星”；“木星”};结尾

UTC时间计算

根据假设的经度调整本地时间为UTC。在这个例子中，假设所有观察都在每天的同一时间采取在晚上8点当地时间。

obs.hour = 20;

对于航位推算过程，用当前位置的估计来更新观测结构。在这种情况下，位置是使用之前的定位、船的速度V和航向T来估计的。

obs经度=长固定点（m）+（1/60）*V（m）*sind（T（m））/cosd（latFix（m））；obs纬度=latFix（m）+（1/60）*V（m）*cosd（T（m））；

使用helper函数将本地时间调整为UTC，如HelperLongitudeHour。此功能根据船舶的估计经度调整UTC观测时间。

obs.UTC = astHelperLongitudeHour（OBS）;

六分仪高度计算

对于每一颗行星，astHelperNauticalCalculation助手函数计算六分仪的测量值，这是康提基号上的工作人员应该测量的。这个函数模拟了行星的实际行为，同时补偿了当地的条件。该函数使用行星星历和ECI到ECEF变换矩阵。分析不包括行星像差、引力光偏转和光现象像差。

obs.Hs = astHelperNauticalCalculation（OBS，latActual，longActual）;

计算位置

下面的计算代替了航海年鉴的使用。它们包括使用行星星历表和ECI到ECEF转换矩阵。

初始化观测对象的赤纬、格林威治时角(GHA)和高度:

obs.declination =零（大小（obs.Hs））;obs.GHA =零（大小（obs.Hs））;obs.altitude =零（大小（obs.Hs））;

计算测量时间的修改儒略日期:

mjd=mjuliandate（obs.UTC）；

计算UT1和UTC之间的区别：

dUT1=deltaUT1（mjd，“行动”，'没有任何'）;

使用美国海军天文台的TAI-UTC（dAT）值计算ECI到ECEF转换矩阵：

dAT=1.4228180；TM=dcmeci2ecef('IAU-76 / FK5', dAT obs.UTC dUT1);

计算地球时间的儒略日期以近似重心动力时间:

jdTT = juliandate (obs.UTC) + (dAT + 32.184) / 86400;

计算每一个天体的赤纬、格林尼治时角和高度:

为了k = 1:长度(obs.object)

计算每个行星的ECI位置：

posECI=Planethemeris（jdTT，“地球”，obs.object{k}，“405”，“公里”）;

计算每个行星的ECEF位置：

posECEF=TM*posECI'；

使用ECEF位置计算格林尼治小时角（GHA）和偏角：

obs.GHA（K）= -atan2d（posECEF（2），posECEF（1））;obs.declination（K）= atan2d（posECEF（3），SQRT（posECEF（1）^ 2 +．..波塞夫(2)(2);；

利用ECEF to LLA变换函数计算地球表面到行星中心的距离:

posLLA=ecef2lla（1000*posECEF'）；观测高度（k）=波士拉（3）；

结尾

行星物体的视力降低

对观察结构阵列中指定的每一颗行星减少视线:

[p, Z] = astHelperNauticalReduction(观察);

从使用以下公式目前的修复程序的最后一次定方位计算增量的经度和纬度。这些公式是基于航海天文历。

AP =总和（COSD（Z）^ 2）;BP =总和（COSD（Z）*信德（Z））;CP =总和（信德（Z）^ 2）;DP =总和（P * COSD（Z））;EP =总和（P *信德（Z））;GP =鸭*的Cp-BP ^ 2;

根据减少量计算纬度和经度的增量:

deltaLongFix =(美联社* Ep-Bp * Dp) / (Gp * cosd (latFix (m)));deltaLatFix = (Cp * Dp-Bp * Ep) /全科医生;

计算出纬度和经度的增量后，将其加到估计的位置，得到观测时间的固定值:

longFix（m+1）=obs.经度+deltaLongFix；纬度固定点（m+1）=观测纬度+三角洲固定点；

结尾

导航解决方案可视化

下图显示的是实际跟踪和瞄准缩减解决方案：金宝搏官方网站

astHelperVisualization（长，纬度，longFix，latFix，“阴谋”）

您可以使用Mapping Toolbox获得描述美洲大陆和法属波利尼西亚视力缩小解决方案的更详细的图表。金宝搏官方网站

astHelperVisualization（长，纬度，longFix，latFix，“地图”）

在经度和纬度相对误差累积为容器帆从卡亚到拉洛阿。这个错误是由于六分仪高度的测量误差小。对于6月9日，该还原法计算海王星和木星真方位角（Z）。对于海王星和木星的真实方位接近180度分开。这种差异会导致在相对误差小高峰。这个错误，但是，仍然是还原方法的误差范围内。

工具书类

[1]鲍迪奇，N.美国实用导航。国家地理空间情报局，2012。

[2]英国水文局。航海年鉴2012商业版。天堂出版公司2011。

[3]城市，肖恩E.和P.肯尼斯Seidelmann。解释补充天文年历。第三版，大学科技图书，2013。

美国海军天文台。