储油罐的变位识别与罐容表标定-数学建模毕业论文.doc

《储油罐的变位识别与罐容表标定-数学建模毕业论文.doc》由会员分享，可在线阅读，更多相关《储油罐的变位识别与罐容表标定-数学建模毕业论文.doc（33页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、储油罐的变位识别与罐容表标定摘要 本文用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题，建立了油量容积与油高的函数关系。考虑到算法的实现性，本文把对体积的多重积分通过降维，转化为单重积分，从而避免多重积分边界难确定问题。 问题一：通过MATLAB编程得到油量容积与油高的数值关系，讨论实验数据与理论数据之间的误差，最后得到结果变位后罐容表标定的油量容积大于实际油量容积。通过建立的罐容模型得出了变位后的罐容表。 问题二，对于实际的油罐，采用分割法降低模型建立的难度和编程的难度，把油罐分成三段，利用问题一中的积分方法，分别建立油的体积与h的函数，相加既得。从附件2的数据表中可以得到出油量和显

2、示油高的数据，由出油量的数值逐步累加可以得到累加出油量的值， （表示油罐内原有的油量）将油位高度为时的理论值与实际值进行逐差比较，方差取得最小值时，即可确定和的值。问题转化为以方差最小为目标函数，和为自变量的最优规划问题。得到。模型的合理性分析：把附录中显示油高与显示油量的数据看做是纵向倾斜角度a和横向偏转角度b 值都为0时的实验值。然后用建好的模型对此数据进行求解，如果得到的结果非常接近0则可以充分证明模型的准确性和可靠性。否则模型不可靠不准确缺乏说服力。最后求解的结果为从而模型的准确性和可靠性得证。最后建立了基于数据特征求变位参数模型。随着有高的上升，液面面积有一最大值，利用附表的数据得到

3、最大面积和对应的高度，进而求出同时求解过程对误差进行了分析与检验，说明结果可靠性很高。结果与上个模型的差异可能主要来源于系统误差，如油罐参数偏差，油罐形变，罐内油管的体积及残留油垢等。关键词：变为识别 罐容表标定 优化问题 数值积分 特征分析一.问题的重述通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以

4、下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。图1是一种典型的储油罐尺寸及形状示意图，其主体为圆柱体，两端为球冠体。图2是其罐体纵向倾斜变位的示意图，图3是罐体横向偏转变位的截面示意图。请你们用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。 （1）为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为a=4.10的纵向变位两种情况做了实验，实验数据如附件1所示。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。（2）对于图1所示的实际储油罐，试建立

5、罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据（附件2），根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。附件1：小椭圆储油罐的实验数据附件2：实际储油罐的检测数据二.问题分析 本问题的关键是找到油罐储油量与油位高度的对应关系。 对于问题一，首先无变位时容易根据立体几何的知识得到储油体积与油位高度的函数关系，并且可以得到无变位时的理论罐容表，与给出的无变位时的数据

6、比较，作误差分析；倾斜角为的纵向变位后，由于随着油位高度的取值不同，储油体积的计算不同，于是分段计算得到储油体积与油位高度的函数关系，同样也可以得到变位后的理论罐容表，与所给数据比较，计算误差。 问题二，由于实际油罐形状复杂，建立罐体变位后标定罐容表的数学模型需要一定的处理技巧。为了减小算法实现的难度，计算体积时把对体积的多重积分通过降维，转化为单重积分，从而避免多重积分边界难确定问题。根据附表的数据，利用最小二乘发求解。使当纵向倾斜角度a和横向偏转角度b 值时理论值与所给数据规律方差最小。利用求解出来的模型求出。纵向倾斜角度a和横向偏转角度b 值都是未知数即使求解出来，由于没有与真实值对比没

7、法说明其正确合理性和可靠性。为此可以把显示油高与显示油量的数据看做是来源一另一种变位情况，只是纵向倾斜角度a和横向偏转角度b 值都为0。然后用建好的模型进行求解如果得到的结果非常接近0则可以充分证明模型的准确性和可靠性。前面的模型都是建立油量体积与液位高度的函数关系，建立模型。同理我们亦可以通过建立油罐内液面面积与液位高度的函数关系，通过这个关系的特征对部分问题题求解，如求解纵向倾斜角度a和横向偏转角度b 值。三.模型假设1.实验数据真实可靠。2.油罐的参数准确，偏差可忽略。3.不考虑外界因素的干扰。四.主要符号说明油面高度油位高度油体截面积油的体积椭圆油罐的长度油截面与轴成的角随着值的变化的

8、函数球冠体的半径油罐中间圆柱体截面半径注：用到的局部变量在文中直接注释。五.模型的建立与求解5.1问题一 为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为的纵向变位两种情况做了实验。由于小 椭圆型储油罐各参数已经确定，于是可以通过各参数计算得到油罐储油量与油位高度的函数关系，以此得到一个理论罐容表。 为此这里先分别对罐体无变位和倾斜角为的纵向变位两种情况建立理论模型，再利用理论模型计算得到的理论值与实际给的数据比较，作误差分析。 5.1.1 无变位的情况时罐容表理论模型如图建立直角坐标系， 椭圆柱体截面方程为 ，其中 =0.89m，

9、=0.6m。令，, 则油面高度： 无变位时油位高度等于油面高度，所以计算油体截面积 所以，油的体积为 ，其中为椭圆油罐的长度。 代入数值计算得到油的容积与油位高度的关系： 5.1.2倾斜角为的纵向变位情况时罐容表理论模型（1）如图所示建立直角坐标系计算油截面面积 截面的参数方程油位高度表示油截面与轴成的角随着值的变化的函数因此，截面积可以表示为： （2）这里我们采用油罐逐步移动法来进行分析求解油的体积。假设有一高度与题目所给的小椭圆油罐相等，以的倾角斜放在水平地面上，对 角线代表油面且与水平地面平齐的长油罐。如图所示建立坐标系，将原点定义在长油罐的右下端，轴沿油罐方向向左。 保持长油罐固定不变

10、，用椭圆柱油罐从右到左逐步移动移动来模拟往油罐中注入油的情况，从而求出在不同情况下小椭圆油罐中油的体积。我们将这种方法称为油罐逐步移动法。 如图所示，表示小椭圆油罐最右端油面到罐底的高度。表示小椭圆油罐最左端到原点的距离。为长油罐的长度，由几何关系可知： 由于当的取值不同时，油体积的计算不同，所以得分别计算，如图所示，共有五种情况： （1） 当时，。但由于确定油位高度的油浮子在距油罐左端0.4米处，此时油浮子不能测到油位高度，故不能对罐容进行标定。 （2）当时， 。 （3）当时， 。 （4）当时，。 （5）当时， 。综上所述:其中 5.1.3 模型求解由于上面所建模型的被积函数较复杂，不能直接

11、进行符号积分，在此我们确保精度在1e-10内采用数值积分的方法，调用MATLAB函数库中的quadl数值积分函数编程求出数值解。 程序见附录9.3.2。5.1.4误差分析.无变位时理论计算结果与实验结果对比、理论值与实验值的差值如图： 可以看出随着油位高度的增加，油的容积差值基本呈线性增长，可以考虑是由于油罐内进油管道和出油管道的体积占用了油的体积，使得计算的理论值比实际值偏高且线性增长。. 有变位且时理论计算结果与实验结果对比、理论值与实验值的差值如图：由上面数据可得理论值始终大于实验值，该实验存在系统误差。5.1.5变位后罐容表标定油位高度（mm）容积(L)油位高度(mm)容积(L)油位高

12、度(mm)容积(L)油位高度(mm)容积(L)103.531044310630.14626101841.7979103112206.26351320665.58086201885.1319203151.234309.974764330701.52566301928.5139303190.114014.75629340737.95846401971.9319403228.6125020.69084350774.85776502015.3729503266.7226027.85416360812.2036602058.8249603304.4217036.3163370849.9747670210

13、2.2759703341.6918046.14242380888.15376802145.7139803378.5119057.39353390926.72176902189.1259903414.86110070.12695400965.66087002232.510003450.7211084.396764101004.9547102275.82410103486.064120100.25414201044.5847202319.08610203520.87130117.74754301084.5357302362.27310303555.114140136.9234401124.7917

14、402405.37210403588.769150157.81844501165.3367502448.37210503621.808160180.25914601206.1557602491.25910603654.2170203.99944701247.2347702534.0210703685.915180228.90664801288.5577802576.64310803716.917190254.88494901330.1117902619.11510903747.171200281.85775001371.8818002661.42311003776.636210309.7608

15、5101413.8548102703.55211103805.266220338.53875201456.0158202745.49111203833.013230368.14265301498.3528302787.22511303859.819240398.52855401540.8518402828.7411403885.618250429.65675501583.4998502870.02211503910.332260461.49065601626.2838602911.05711603933.858270493.99675701669.198702951.8311703956.05

16、6280527.14385801712.2088802992.32611803976.655290560.90245901755.3238903032.53111903995.537300595.24526001798.5249003072.42712004012.7455.2 问题二 对于图1所示的实际储油罐，先建立罐体变位后标定罐容表的理论数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。通过误差修正后得到修正模型，根据附表的数据求解出纵向倾斜角度a和横向偏转角度b 值，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步利用附件2中的实

17、际检测数据来分析检验模型的正确性与方法的可靠性。这里根据附表的数据特征，在没有利用油体积与液位高度的理论模型情况下得出纵向倾斜角度a和横向偏转角度b值。5.2.1实际储油罐变位后罐容表标定理论模型考虑到纵向倾斜角和横向偏转角，将油罐分为三段分别进行容积的计算，这三段分别为中间的圆柱体和两端的两个球冠体。为了方便计算分别建立油的体积与h的关系函数1.对于中间圆柱体容积的计算按照前面的方法建立直角坐标系则 对于圆柱体沿方向上的截面，如图所示建立另一个直角坐标系。可得其中沿轴方向上截面的面积可表示为：则容积：其中。2.对于左端球冠体容积的计算确定球冠体的半径： 如图建立直角坐标系，则此时左端球冠体的

18、表面方程为：。确定油面方程：这里确定球冠体内油面距平面上轮廓线的距离（如图所示）： 所以可推导出的表达式：沿轴方向上球冠体的截面半径：因此可以求得球冠体沿轴方向上油的截面面积： 所以 3. 对于右端球冠体容积的计算确定球冠体的半径： 建立直角坐标系，则此时右端球冠体的表面方程为：，确定油面方程：沿轴方向上球冠体的截面半径：确定球冠体内油面距平面上轮廓线的距离（如图所示）： 其中因此可以求得球冠体沿轴方向上油的截面面积： 所以综上所述:油的容积： 油位高度与的关系为：因此可以确定油的容积与油位高度之间的函数关系。5.2.2模型求解：算法实现过程：由于被积函数的表达式过于复杂，不好直接求出原函数。

19、我们采用数值积分的方法，在的误差范围内，对于任意一个油位高度值，分别求出上述三段油的容积，叠加后即为总容积。对于中间圆柱体容积的计算，采用函数库中的数值积分函数进行积分求解。对于左右两端球冠体部分的容积计算，采用定步长法进行积分求解。（程序见附录 ）确定实际储油罐变位后的和值：从附件2的数据表中可以得到出油量和显示油高的数据，由出油量的数值逐步累加可以得到出油量的累加值，而累加值与油位高度之间有如下关系： （表示油罐内原有的油量）即为油位高度为时，罐内油的容积。将油位高度为时的理论值与实际值进行逐差比较，当它们的波动最小，即方差取得最小值时，问题转化为优化问题，即可确定和的值。没加油之前的方差

20、加油之后数据的方差为目标函数为通过软件编程进行搜索方差取得最小值时的和的值。结果如下： 0.00105L5.2.3模型的检验：附录中的试验给出了显示油高和显示油量容积数据，经过分析此时的数据应该是建立在的基础上得到的，也就是说这是无变位时候的罐容表。利用这两列数据作为评判模型准确性的标准，代入模型中求得此时的和值，比较求得的和值与之间的误差，用于衡量模型的准确性。在不同的油位高度下，将用模型求得的油量容积与显示油量容积逐次做差，当差值的方差最小时得到的和值即为所求的和。根据上述模型的算法，用软件编程求得的和的值如下表：00.0030L用模型求得油的容积与实际值的图像如下：从上面两条曲线中可以看

21、出，用模型模拟所得的油的容积与实际值非常接近，且角度的误差很小，能达到以内，由此说明我们所建立的模型准确性很高，且方法较优。5.2.4变位后罐容表标定值 由于和值已经确定，并且经过检验准确度足够高，由此得到罐体变位后的储油体积与油位高度的对应关系，得到如下罐容表：油位高度（mm）容积(L)油位高度（mm）容积(L)100376.168160033053.952001125.172170035863.843002216.817180038654.514003694.801190041410.525005423.54200044116.126007361.309210046755.12700947

22、7.009220049310.6380011745.11230051764.8690014143.48240054098.7100016652.2250056291.21110019252.89260058318.69120021928.27270060153.11130024661.84280061758.89140027437.66290063083.84150030240.14300066338.545.3基于实验数据特征的求解模型我们的目的是找到油体积与油位高度之间的关系，然后得到罐容表。对于问题中的实际储油罐，我们可以建立油体积与油位高度的函数关系。无变位时，油位高度等于油面高度，其

23、对应的油面面积 即为图像的斜率。当储油罐纵向倾斜时，油体积与油位高度的函数关系为，则同理，无变位油面高度与变位后油面高度的关系 因为 对（1）微分得 所以，即 对于争个油罐如图，经分析知道，当油面过油罐几何中心时油面面积最大。此时必有一与之对应的油位高，即 最大。由几何知识易知， 又因为所以，。实际题目只给了部分，实际工作情况时，出油量显示油高的数据。为了方便数据处理，先对出油量进行累加，得到累加出油量。由于在数据采集的中间左右的位置进行了一次性补充进油作为分界，得到两组数据（可以看作为两次实验），分别累加得到累加出油。出油量 ， 假设储油罐初始油量为，则油罐中油量 对上式微分， 所以，可以利

24、用累加出油量的规律来可得最大油面面积和取的最大油面面积时的油位高度。 实际的数据都是离散的，且由于各种偶然误差会导致数据与真实值有偏差。为利用这些离散量得到结果，且是结果准确度和可信度尽可能高，必先对数据进行分析和处理。通过累加出油量得到 分别绘制图像如图，发现与非常接近，定性的说明了误差很小。理论上，但由于由误差，则有 对于可以看作重复试验，由此具有相同的方差，则方差为。由于每个数据对应的不一样，无法进行直接比较分析。为了对比较分析，先对进行三次样条插值后得到两组新的累加出油量函数。则方差，通过 分别得到两组油面面积值，并单独绘制图，如图所示，由于数据波动，不能直接得到的最大值及取得该值时的

25、高度值。图大致呈抛物线。因此可以利用二次曲线进行拟合，但由于图像并非严格的抛物线，所以拟合区域不能过大，但也不能过小，可取1000m 2100下面为三组数据拟合的结果： 拟合参数变位1-0.0000388325 0.1220618572 44.8547865774 变位2-0.0000388608 0.1221931000 44.7185187579 平均-0.0000388466 0.1221274786 44.7866526677 拟合的均方差分别为拟合方差变位10.4856357863 变位20.7869125339 平均0.6362741601 拟合图像 由上述拟合方程可得最大面积和对

26、应的高度为高度体积差面积变位11571.6458162532 140.7737901663 28.1547580333 变位21572.1898744187 140.7738960625 28.1547792125 平均1571.9178453360 140.7738431144 28.1547686229 对于未变位时，最大面积(R为两端球冠的半径R=1.625)计算得所以，同理得到数据6.4798385545 6.4794590690 6.4796488145 由可得到2.0385287339 2.0539970818 2.0462629422 上述求解方法为默认给定油罐的参数是精确的，没

27、偏差。但对于一个实际问题，误差是不可避免的，是客观存在的。为了得到一个更准确的解，同时也可以证明模型的准确性和合理性。这里，用上述模型对无变位时油量与油位高度的关系进行研究。同前面处理方法一样，先插值，分别得到两组罐容表，和，以为步长得一纵向量，方差分析得到: 所以理论所得的罐容表非常准确。误差来源，主要是数据的舍入误差。同理，进行逐差，求得 对各面积二项式拟合得到拟合方程参数分别为无变位1-0.0000391858 0.1173328253 53.8727503088 无变位2-0.0000391865 0.1173348929 53.8710814132 平均-0.0000391861 0

28、.1173338591 53.8719158610 拟合方差分别为无变位10.0006971681 无变位20.0007059390 平均0.0007015536 由上表可知方程拟合的非常好 最大面积及取得最大面积时的油位高度：平均相对变化28.3408392855 28.3408451548 28.3408422201 0.0178739409% 1497.1334046651 1497.1359696064 1497.1346871357 0.1910208576 %由于这里只有无变位时的数据，所以即为油罐的半径。由上面计算可得，油罐半径变小 0.191 %，而数据误差因此说明油罐的参数的

29、确有微小的偏差，这一微小的误差充分说明模型的准确性。基于上述分析，进一步修正与的值，则 计算结果如下：平均值相对修正量2.1273356833 1.3637% 6.5691857259 3.823% 微小的相对修正量同时也证明通过这个模型的准确性和可靠性下面为油罐变位与没变位时的体积差的函数图的比较6.结果分析与检验对于问题一通过建立油量容积与油高的数值关系，得到油高对应的理论油量容积，与实验数据比较得到一总是大于0的偏差平均误差约2%，说该实验具有系统误差。通过变位和无变位得到在相同油高下变位的油罐油量容积大。且最大值在油高0.6m左最大。通过建立的罐容模型对罐容表标定结果比较理想。问题二：

30、对于问题二我们分别从两个不同角度建立了两个不同的模型。模型一通过直接推导出油量容积与油高的数学模型，求解出两个参数0.00105L当无变位时油量容积与油高的关系极为显示油量与显示油高的数据的关系，此时可以吧无变位时的情况看做一种特殊的变为情况，把模型待入求其变位参数结果为00.0030L参数非常接近0，所以可以说明模型的可靠性。另一个是基于数据的特征，基于对特征的分析建立模型，结果为2.1273356833 6.5691857259 结果与第一个模型很接近，也进一步证明两个模型都具有合理性。7.模型的优缺点评价及改进7.1模型的优点（1）整个模型建立和求解过程中，都采用了降维的思想，将二重积分

31、转化为一重积分，避免多重积分边界难以确定的问题。（2）在处理问题一的有倾角变位问题时，通过引入一个虚拟油罐，采用了实际油罐在虚拟油罐上逐步移动法来进行分段求解，这一方法有效且直观的解决了分段积分问题。（3）在不能进行符号积分时，在确保精度的条件下采用数值积分的方法，利用数值模型的思想解决问题。7.2模型的缺点（1）在问题一的求解时，没有充分考虑实际值与理论值之间的误差。（2）问题二的编程求解时，和的值采用的是人工搜索，程序设计不足。7.3模型的改进（1）要通过附录中给出的实测数据，充分考虑实际值与理论计算值之间的误差关系，通过曲线拟合的方法进行误差修正，从而得到更加精确的罐容表。（2）为了便于

32、算法的实现，这里还可以尝试其他的积分方法。8.参考文献1同济大学数学系.高等数学(下册-6版).北京:高等教育出版社,2007.6(2008重印).2王绵森,马知恩.工科数学分析基础(上册)(第二版).北京:高等教育出版社,2006.2(2008重印).3王连堂.数学建模.陕西:陕西师范大学出版社,2008.5.4黄永安,李文成,高小科.Matlab7.0/Simulink6.0应用实例仿真与高效算法开发.北京:清华大学出版社,2008.6.9.附录9.1 附件问题A附件问题B附件9.2 示意图9.3 问题一程序yi.mclcclearv=(501003356.913406.913456.91

33、.3506.913556.913606.913656.913706.91+262)./1000;h=159.02176.14192.59.949.05963.80978.91994.431010.431026.991044.251062.371081.591102.331125.321152.361193.49./1000;b=0.6;va=size(h);for i=1:length(h) va(i)=1.3083*(asin(h(i)-b)/b)+(h(i)-b)/b*sqrt(1-(h(i)/b-1)2)+pi/2); d(i)=va(i)-v(i);endva;%m=va-v;%plo

34、t(h,v,r,h,va,b);plot(h.*1000,d.*1000)dyibd3.mclcclear%容积与长度的关系format longh=0:0.01:1.2;a=0.89;b=0.6;r=4.1*pi/180;l1=2.45;lm=2*b/tan(r);v=size(h);xm=size(h);for i=1:length(h) if h(i)=lm*tan(r) xm(i)=h(i)/tan(r)+0.4; endendxm;for i=1:length(h) q=i; s=inline(0.89*0.6*(asin(x.*tan(4.1*pi/180)/0.6-1)+0.5*

35、sin(2*asin(x.*tan(4.1*pi/180)/0.6-1)+pi/2);%s=inline(a*b*(asin(x*tan(r)/b-1)+0.5*sin(2*asin(x*tan(r)/b-1)+pi/2);if xm(i)0 v(i)=quadl(s,0,xm(i),1e-12);elseif xm(i)2.45 v(i)=quadl(s,xm(i)-l1,xm(i),1e-12);elseif xm(i)lm v(i)=quad(s,xm(i)-l1,lm,1e-12)+pi*a*b*(xm(i)-lm);endendv=v*1000;vs=747.86797.86847.

36、86897.86947.86997.861047.861097.791147.79.3097.733147.733197.733247.733297.733299.74+215;j=v-vsvar(j)%plot(h,v,-r,h,vs,-b)%figure 2plot(h,j,-o)9.4问题二程序fangcha.mclearclch1=60.00 2632.23 60448.88 149.09 2624.30 60311.43 68.45 2620.67 60248.03 199.27 2610.29 60065.11 。149.83 437.81 5477.45 118.46 430.8

37、6 5347.56 65.81 426.96 5275.11 115.30 420.01 5146.78 57.09 416.53 5082.90 43.13 413.98 5036.26 ;h1(:,1)=cumsum(h1(:,1);h2(:,1)=cumsum(h2(:,1);h3=-h1(:,1)./1000;h4=-h2(:,1)./1000;h5=h3;h4;r=2.12*pi/180;n=1.5;R=1.625;bt=4.18*pi/180;h=(h1(:,2);h2(:,2)./1000-n).*cos(bt)+n+2*tan(r);l1=8;lm=2*n/tan(r);x=h

38、./tan(r);s=inline(1.52*(asin(x*tan(2.12*pi/180)/1.5-1)+0.5*sin(2*asin(x*tan(2.12*pi/180)/1.5-1)+pi/2);%s=inline(n2*(asin(x.*tan(r)/n-1)+0.5*sin(2*asin(x.*tan(r)/n-1)+pi/2);for i=1:length(h) if x(i)0 v(i)=quadl(s,0,x(i);elseif x(i)8 v(i)=quad(s,x(i)-l1,x(i);elseif x(i)lm v(i)=quad(s,x(i)-l1,lm,1e-12)

39、+pi*R2*(x(i)-lm); endendfor j=1:length(h) d(j)=h(j)-(R-1)*tan(r)-sqrt(2*R-1);for i=1:999 t1=(R-1)+i/1000; mt1=asin(t1*tan(r)+d(j)/sqrt(R2-t12);ht1(i)=t1*tan(r)+d(j)+sqrt(R2-t12);Rt1(i)=sqrt(R2-t12);if ht1(i)0&ht1(i)=2*Rt1(i) st1(i)=pi*(R2-t12);elseif ht1(i)=0 st1(i)=0;endv1(i)=st1(i)/1000;endvz(j)=sum(v1);endfor k=1:length(h) if h(k)0&ht2(i)=2*Rt2(i) st2(i)=pi*(R2-t22);elseif ht2(i)=