معادلة فنسنتي

صيغتا فينسنتي هما طريقتان تكراريتان مترابطتان تُستخدمان في علم المساحة
لحساب المسافة بين نقطتين على سطح إهليلجي، طورهما ثاديوس فينسنتي (1975).
وتستندان إلى افتراض أن شكل الأرض إهليلجي مفلطح، وبالتالي فهي أكثر دقة من
الطرق التي تفترض كروية الأرض، مثل مسافة الدائرة العظمى.

الطريقة الأولى (المباشرة) تحسب موقع نقطة تبعد مسافة وسمتًا محددين عن نقطة أخرى.
أما الطريقة الثانية (العكسية) فتحسب المسافة الجغرافية والسمت بين نقطتين محددتين.
وقد استخدمت هذه الطرق على نطاق واسع في علم المساحة نظرًا لدقتها التي تصل
إلى 0.5 مم (0.02 بوصة) على سطح الأرض الإهليلجي.

## خلفية

كان هدف فينسنتي هو التعبير عن الخوارزميات الحالية للجيوديسيات على شكل
إهليلجي بصيغة تقلل من طول البرنامج (فينسنتي 1975). يذكر تقريره غير
المنشور (فينسنتي 1975) استخدام حاسبة مكتبية من طراز Wang 720،
والتي كانت تحتوي على كيلوبايتات قليلة من الذاكرة.

للحصول على دقة جيدة للخطوط الطويلة، استخدم فينسنتي الحل الكلاسيكي
الذي وضعه ليجاندر (1806) وبيسل (1825) وهيلميرت (1880) والمبني
على الكرة المساعدة.

اعتمد فينسنتي على صياغة هذه الطريقة التي قدمها رينسفورد (1955).

أظهر ليجاندر أنه يمكن رسم خريطة دقيقة لجيوديسية إهليلجية على دائرة
عظمى على الكرة المساعدة عن طريق رسم خريطة لخط العرض الجغرافي
على خط عرض مخفض، وضبط سمت الدائرة العظمى على أنه يساوي
سمت الجيوديسية.

يعطى بعد ذلك خط الطول على الشكل الإهليلجي والمسافة على طول الخط
الجيوديسي بدلالة خط الطول على الكرة وطول القوس على طول الدائرة
العظمى باستخدام تكاملات بسيطة.

وقد قدم بيسل وهيلميرت سلاسل متقاربة سريعة لهذه التكاملات،
مما يسمح بحساب الخط الجيوديسي بدقة متناهية.

من أجل تقليل حجم البرنامج، أخذ فينسينتي هذه السلاسل، وأعاد توسيعها باستخدام الرمز الأول من كل سلسلة كمعامل صغير، وقام بتقصيرها إلى O(f^3) أدى ذلك إلى تعابير مُدمجة لتكاملات الطول والمسافة. وُضعت التعبيرات بصيغة **هورنر** (أو المُتداخلة)، لأن هذا يسمح بتقييم كثيرات الحدود باستخدام سجل مؤقت واحد فقط. وأخيرًا، استُخدمت تقنيات تكرارية بسيطة لحل المعادلات الضمنية بالطريقتين المباشرة والعكسية؛ ورغم بطء هذه الطرق (وفي حالة الطريقة العكسية، لا تتقارب أحيانًا)، إلا أنها تُؤدي إلى أقل زيادة في حجم الكود.

## '''الرموز والمصطلحات'''

## '''حل المسألة العكسية'''
بالنظر إلى إحداثيات النقطتين **( φ₁ ، L₁ )** و**( φ₂ ، L₂ )**، فإن المسألة العكسية تحدد السمتين **α₁** ، **α₂** والمسافة الإهليجية **s**.

أحسب **U₁** و **U₂** و **L** ، ثم اضبط القيمة الأولية لـ **λ** **= L**. ثم قيّم المعادلات التالية بشكل متكرر حتى تتقارب **λ** :

: \sin\sigma = \sqrt{\left(\cos U_2 \sin\lambda\right)^2 + \left(\cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos\lambda\right)^2}
: \cos\sigma = \sin U_1 \sin U_2 + \cos U_1 \cos U_2 \cos\lambda\,
:
: \sigma = \operatorname{arctan2}\left(\sin\sigma, \cos\sigma\right)
: يتم تجنب حساب **σ** (سيجما) مباشرةً مِن دالة جيب **σ** أو جيب تمام **σ** للحفاظ على الدقة العددية بالقرب مِن القطبين وخط الاستواء.
:
: \sin\alpha_0 = \frac{\cos U_1 \cos U_2 \sin\lambda}{\sin\sigma}
: إذا كانت المسافة الزاوية **σ = 0** فإن قيمة جيب **α₀** تصبح غير مُعَرَّفة (غير محددة). وتُمثل هذه الحالة نقطة نهاية مُطابِقة لنقطة البداية، أو مُعاكِسة لها قُطريًا (على الطرف الآخر مِن القطر).
:
: \cos^2\alpha_0 = 1 - \sin^2\alpha_0
:
: \cos\left(2\sigma_\text{m}\right) = \cos\sigma - \frac{2\sin U_1 \sin U_2}{\cos^2\alpha_0} = \cos\sigma - \frac{2\sin U_1 \sin U_2}{1 - \sin^2\alpha_0}
: عندما تقع نقطتا البداية والنهاية على خط الاستواء، فإن **C = 0** ولا يتم استخدام قيمة جيب تمام **(2σₘ)**. في هذه الحالة، تكون القيمة الحدية لجيب تمام **(2σₘ) = -1**.
:
: C = \frac{f}{16} \cos^2\alpha_0 \left[4 + f\left(4 - 3\cos^2 \alpha_0\right)\right]
: \lambda = L + (1 - C)f \sin \alpha_0 \left\{ \sigma + C \sin \sigma \left[\cos\left(2\sigma_\text{m}\right) + C \cos\sigma \left(-1 + 2\cos^2\left(2\sigma_\text{m}\right)\right)\right]\right\}

عندما تتقارب **λ** إلى الدرجة المطلوبة مِن الدقة (10^−12 يتوافق مع ما يقرب من 0.006 مم)، يتم إيقاف التكرار ثم تطبق المعادلات التالية باستخدام **λ** الأخيرة:

: \begin{align}
           u^2 &= \cos^2\alpha_0 \left(\frac{a^2 - b^2}{b^2}\right) \\
             A &= 1 + \frac{u^2}{16384} \left(4096 + u^2 \left[-768 + u^2 \left(320 - 175u^2\right)\right]\right) \\
             B &= \frac{u^2}{1024} \left(256 + u^2 \left[-128 + u^2 \left(74 - 47u^2\right)\right]\right) \\
  \Delta\sigma &= B \sin\sigma \left\{
                    \cos(2\sigma_\text{m}) + \frac{1}{4} B \left(
                      \cos\sigma \left[-1 + 2 \cos^2\left(2\sigma_\text{m}\right)\right] -
                      \frac{1}{6} B \cos\left[2\sigma_\text{m}\right]\left[-3 + 4\sin^2\sigma\right]
                        \left[-3 + 4\cos^2\left(2\sigma_\text{m}\right)\right]
                    \right)
                  \right\}  \\
             s &= b A(\sigma - \Delta \sigma) \, \\
      \alpha_1 &= \operatorname{arctan2}\left(
                    \cos U_2 \sin\lambda, \cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos\lambda
                  \right) \\
      \alpha_2 &= \operatorname{arctan2}\left(
                    \cos U_1 \sin\lambda, -\sin U_1 \cos U_2 + \cos U_1 \sin U_2 \cos\lambda
                  \right)
\end{align}

بين نقطتين متقابلتين تقريبًا، قد تفشل الصيغة التكرارية في التقارب؛ سيحدث هذا عندما يكون التخمين الأول عند **λ** كما تم حسابه بواسطة المعادلة أعلاه أكبر من **π** في القيمة المطلقة؛ وعلاج هذا في الجزء الثاني الذي كتبه فينسنتي بورقته الثانية. 

## '''حل المسألة المباشرة'''
بالنظر إلى نقطة البداية **( φ₁ ، L₁ )** والسمت الأولي **α₁** ، والمسافة **s** ، على طول الخط الجيوديسي، فإن المشكلة هي إيجاد نقطة النهاية **( φ₂ ، L₂ )** والسمت، α₂ .

ابدأ بحساب ما يلي:

\begin{align}
           U_1 &= \arctan\left[(1 - f)\tan\varphi_1\right] \\
      \sigma_1 &= \operatorname{arctan2}\left(\tan U_1, \cos\alpha_1\right) \\
    \sin\alpha_0 &= \cos U_1 \sin\alpha_1 \\
           u^2 &= \cos^2\alpha_0 \left(\frac{a^2 - b^2}{b^2}\right)
                = \left(1 - \sin^2\alpha_0\right)\left(\frac{a^2 - b^2}{b^2}\right) \\
             A &= 1 + \frac{u^2}{16384} \left(4096 + u^2 \left[-768 + u^2(320 - 175u^2)\right]\right) \\
             B &= \frac{u^2}{1024} \left(256 + u^2\left[-128 + u^2\left(74 - 47u^2\right)\right]\right)
\end{align}

ثم باستخدام القيمة الأولية **σ** ، كرر المعادلات التالية:

\begin{align}
  2\sigma_\text{m} &= 2\sigma_1 + \sigma \\
      \Delta\sigma &= B \sin\sigma \left\{
                          \cos\left(2\sigma_\text{m}\right) + \frac{1}{4} B \left(
                            \cos\sigma \left[-1 + 2\cos^2\left(2\sigma_\text{m}\right)\right] -
                            \frac{1}{6} B \cos\left[2\sigma_\text{m}\right]\left[-3 + 4\sin^2\sigma\right]
                              \left[-3 + 4\cos^2\left(2\sigma_\text{m}\right)\right]
                          \right)
                        \right\} \\
            \sigma &= \frac{s}{bA} + \Delta\sigma
\end{align}

إلى الحد الذي لا يكون معه تغيير كبير في **σ** وبمجرد الحصول على **σ** بدقة كافية (فرق 10^-12)، قم بتعيين:

\begin{align}
    \varphi_2 &= \operatorname{arctan2}\left(
                \sin U_1 \cos\sigma + \cos U_1 \sin\sigma \cos\alpha_1,
                (1 - f) \sqrt{\sin^2\alpha_0 + \left(\sin U_1 \sin\sigma - \cos U_1 \cos\sigma \cos\alpha_1\right)^2}
              \right) \\
   \lambda &= \operatorname{arctan2}\left(
                \sin\sigma \sin\alpha_1,
                \cos U_1 \cos\sigma - \sin U_1 \sin\sigma \cos\alpha_1
              \right) \\
         C &= \frac{f}{16} \cos^2\alpha_0 \left[4 + f\left(4 - 3\cos^2\alpha_0\right)\right] \\
         L &= \lambda - (1 - C)f \sin\alpha_0 \left\{\sigma + C\sin\sigma \left(\cos\left[2\sigma_\text{m}\right] + C\cos\sigma\left[-1 + 2\cos^2\left(2\sigma_\text{m}\right)\right]\right)\right\} \\
       L_2 &= L + L_1 \\
  \alpha_2 &= \operatorname{arctan2}\left(
                 \sin\alpha_0,
                -\sin U_1 \sin\sigma + \cos U_1 \cos\sigma \cos\alpha_1
              \right)
\end{align}

إذا كانت نقطة البداية عند القطب الشمالي أو الجنوبي، فإن المعادلة الأولى تكون غير محددة. إذا كان السمت الابتدائي شرقًا أو غربًا، فإن المعادلة الثانية تكون غير محددة. إذا استُخدمت دالة atan2 القياسية ثنائية الوسيطات، فعادةً ما تُعامل هذه القيم بشكل صحيح.

مع ملاحظة أن دالة **atan2** تعطي قيمة عددية متممة (أي مجموعهما 90 درجة) للقيمة التي تعطيها دالة **atan** العادية مع عدم اختلاف الإشارة لو كانت النتيجة سالبة؛

**والجدير بالذكر** أن البعض يفضل تقديم العنصر الثاني في الدالة على العنصر الأول للحصول على نتيجة صحيحة ومباشرة بهذا الشكل **(x,y)atan2(cos,sin)** ، **atan2** بينما الصورة التقليدية هي **(y,x)atan(sin,cos)** ، **atan** بمعنى أن عند التطبيق في برنامج الإكسيل وما يشابهه يفضل عكس الترتيب السابق في المعادلات أعلاه للحصول على نتيجة مباشرة وبدون إشارات سالبة.== **تعديل فينسنتي** ==
في رسالته إلى مجلة **Survey Review** في عام 1976، اقترح فينسينتي استبدال تعبيراته المتسلسلة لـ **A** و **B** بصيغ أبسط باستخدام معامل التوسع لـ **Helmert** **k₁** :

: A = \frac{1 + \frac{1}{4} {k_1}^2}{1 - k_1}
: B = {k_1}\left(1 - \frac{3}{8}{k_1}^2\right)
: حيث
: k_1 = \frac{\sqrt{1 + u^2} - 1}{\sqrt{1 + u^2} + 1}

## '''نقاط متقابلة تقريبًا'''
كما هو مذكور أعلاه، فإن الحل التكراري بالصيغة القياسية للمشكلة العكسية يفشل في التقارب أو يتقارب ببطء بالنسبة للنقاط المتقابلة تقريبًا. ومن الأمثلة على التقارب البطيء **(φ₁, L₁)** = (0°, 0°) و **(φ₂, L₂)** = (0.5°, 179.5°) للقطع الناقص **WGS84**. ويتطلب هذا حوالي 130 تكرارًا لإعطاء نتيجة دقيقة تصل إلى 1 مم. وبناءً على كيفية تنفيذ الطريقة العكسية، قد تعيد الخوارزمية النتيجة الصحيحة (19936288.579 مترًا) أو نتيجة غير صحيحة أو مؤشر خطأ. ويتم توفير مثال على نتيجة غير صحيحة بواسطة الأداة المساعدة عبر الإنترنت **NGS** ، والتي تعيد مسافة أطول بحوالي 5 كم. واقترح فينسينتي طريقة لتسريع التقارب في مثل هذه الحالات **(راب، 1993)**.

من أمثلة فشل الطريقة العكسية في التقارب: **(φ₁, L₁)** = (0°, 0°) و **(φ₂, L₂)** = (0.5°, 179.7°) للقطع الناقص **WGS84**. في تقرير غير منشور، قدم **فينسينتي (1975b)** مخططًا تكراريًا بديلًا للتعامل مع مثل هذه الحالات. يُقارب هذا المخطط النتيجة الصحيحة 19944127.421 مترًا بعد حوالي 60 تكرارًا؛ ومع ذلك، في حالات أخرى، يلزم آلاف التكرارات.

قام **كارني (2013)** بإعادة صياغة المشكلة العكسية كمشكلة إيجاد الجذر أحادية البعد؛ ويمكن حلها بسرعة باستخدام طريقة نيوتن لجميع أزواج نقاط الإدخال.

## '''تفصيل الحل العكسي الجيوديسي بين النقاط المتقابلة (فينسنتي 1975ب)'''
**مقدمة**:

تعتمد الطرق الكلاسيكية لحل المسألة الجيوديسية العكسية على منهج "**بيسيل**" (يشار إليها فيما بعد بالطرق القياسية) وتتضمن تحديدًا تكراريًا لـ **λ** (فرق خط الطول على الكرة المساعدة).

من المعروف أن هذه الطرق تفشل في الحالات المتقابلة أو شبه المتقابلة، أي عندما تكون القيمة المطلقة لـ **λ** أكبر من **π** (مثال ذلك الفرق بين خطيّ طول 91° غربًا و 91° شرقًا، |L_2 - L_1| = |-91^\circ - 91^\circ| = 182^\circ)؛

قدم **ثين (1967)** حلاً للمسألة المتقابلة يقتصر على النقاط الواقعة على خط الاستواء، كما أشار إليه **راب (1975)** ؛

أما طريقة **سايتو (1970)**، والتي صُممت للاستخدام العام، فهي صالحة فقط عندما تكون نقاط نهاية الخط في خطوط عرض متقابلة (\phi_2 = -\phi_1) ؛

الطريقة المقترحة، والمستمدة من إعادة ترتيب معادلات وردت في بحث سابق **(فينسنتي، 1973)**، هي طريقة عامة تمامًا. سيظهر أن النهج الذي يعتمد عليه الحل المتقابل المقترح يمكن استخدامه أيضًا في الحالات المشابهة للحل حول الجانب الخلفي للقطع الناقص بين نقاط متقاربة.

الرموز المستخدمة:

المحوران الأكبر والأصغر للقطع الناقص: a, b

التفلطح: f = \frac{a - b}{a}

مربع الانحراف الثاني: \epsilon = \frac{a^2 - b^2}{b^2}

ويرمز للرمز \epsilon بالرمز e'^2 أيضًا وهذا ما نجده مثلًا في معادلات **كارني**، حيث أن e^2 = \frac{a^2 - b^2}{a^2}،

خط العرض الجيوديسي (موجب شمال خط الاستواء): \phi_1, \phi_2

فرق خط الطول الجيوديسي (موجب شرقًا): L

المسافة الجيوديسية (طول الخط الجيوديسي): s

السمت الجيوديسي باتجاه عقارب الساعة من الشمال؛ \alpha_1, \alpha_2 في اتجاه امتداد (P_1 P_2):

سمت الخط الجيوديسي عند خط الاستواء: \alpha_0

خط العرض المخفض: U يُعرّف بـ \tan U = (1 - f) \tan \phi

فرق خطيّ الطول على كرة مساعدة: \lambda

المسافة الزاوية بين P_1 و P_2 على الكرة: \sigma

المسافة الزاوية على الكرة من منتصف الخط إلى تقاطعه مع خط الاستواء: \sigma_m

نحن ندرس خطًا جيوديسيًا ينطلق من النقطة P_1 في نصف الكرة الشمالي عند خط عرض \phi_1 وسمت ابتدائي \alpha_1 = \frac{\pi}{2}. أثناء تتبع هذا الخط، يتغير سمت الخط وفقًا لمعادلة كليرو: \cos U_i \sin \alpha_i = \text{ثابت}

وفي هذه الحالة الخاصة: \sin \alpha_i = \frac{\cos U_1}{\cos U_i} = 1

الخط الجيوديسي سيعبر خط الاستواء بسمت \alpha_0 محدد بالعلاقة: \sin \alpha_0 = \cos U_1

حتى يصل إلى المنطقة المقابلة عند النقطة P_a حيث: \phi_a = -\phi_1 (خط عرض معاكس)

\lambda = \sigma = \pi (فرق الطول والمسافة الزاوية = 180°)

(السمت النهائي) \alpha_a = \frac{\pi}{2}

حتى هذه النقطة، يمثل الخط أقصر مسافة على سطح القطع الناقص، مما يجعله خطًا جيوديسيًا ويمكن إظهار أنه حتى حدود (f^3)

ويمكن التعبير عن خط الطول الجيوديسي عند هذه النقطة:

L_a = \pi - f \cos U \, \pi \left\{1 - \frac{1}{16} f \sin^2 U \left[4 + f (4 - 3 \sin^2 U)\right]\right\}

حيث: |U| = U_1 = U_2

هذا أكبر خط طول جيوديسي بنهاية الخط الجيوديسي الذي يبدأ بـ \alpha_1 = \frac{\pi}{2}

الآن سنختار النقطة P_2 بحيث (\phi_2 = \phi_a أو \phi_2 \approx \phi_a) و \pi > L > L_a،

الخط الجيوديسي الأقصر بين P_1 و P_2 لابد أن يبدأ بسمت مختلف عن (\frac{\pi}{2}) بسبب تفلطح الإهليج وهناك ثلاث احتمالات:

إذا كان \phi_1 + \phi_2 > 0 فإن: \frac{\pi}{2} > \alpha_1 > 0

إذا كان \phi_1 + \phi_2  فإن: \pi > \alpha_1 > \frac{\pi}{2}

إذا كان \phi_1 + \phi_2 = 0 فإن: \lambda = \sigma = \pi

وهناك حلّان للسمت:

إذا تتبعنا خطًّا يبدأ من \phi_1 > 0 وبسمت \frac{\pi}{2} > \alpha_1 > 0 فهذا المسار سيبلغ خط عرضه الأكبر

وسيعبر خط الاستواء باتجاه الجنوب الشرقي ويتجه نحو خط العرض المقابل حيث (-\phi_a = \phi_1).

يمكن حساب إحداثيات النقاط على طول الخط والتحقق من النتائج عن طريق **الحلول العكسية (Inverse Solutions)** باستخدام الطرق القياسية حتى الوصول إلى المنطقة المقابلة عند النقطة P_0

إذا قمنا بتمديد هذا المسار إلى خط عرض \phi_a واخترنا أي نقطة P_2 عليه بحيث \phi_2 \ge \phi_a، نجد من الحل المباشر أن (\lambda  فرق الطول على الكرة المساعدة أقل من 180° و (\sigma  المسافة الزاوية أقل من 180°، في الحل العكسي \lambda يتذبذب بين قيم أكبر وأصغر من \pi بدون تقارب إلى أي قيمة.

يمكن استنتاج وجود منطقة تعتبر أي نقطة فيها "مضادة" بالنسبة لـ P_1 حتى لو لم يكن هذا صحيحًا هندسيًا، في نصف الكرة الشرقي المقابل لـ P_1 يجب أن تمتد طوليًا من L_a إلى \pi كما هو موضح بالمعادلة:L_a = \pi - f \cos U \, \pi \left\{1 - \frac{1}{16} f \sin^2 U \left[4 + f (4 - 3 \sin^2 U)\right]\right\}

ويجب أن يحددها منحنيان متماثلان حول \phi_a متوافقين مع

\frac{\pi}{2} > \alpha_1 > 0 و \pi > \alpha_1 > \frac{\pi}{2}

إنها صورة مرآوية في نصف الكرة الغربي تتبع منطقيًا قيم \alpha_1 في الربعين الثالث والرابع؛ لم يتم تحديد الحدود الدقيقة للمنطقة المضادة بعد.

### '''حل عكسي'''
أثناء الحسابات الأولية إذا كانت \lambda أكبر من \pi كقيمة مطلقة، يجب استخدام طريقة مختلفة؛

فيما يلي وصف لطريقة فعالة بشكل عام:

نستخدم L' = \pi - L، عندما L > 0

أو L' = -\pi - L، عندما L 

استخدام L' و \lambda' بدلاً من L و \lambda يهدف إلى تجنب فقدان الأرقام المعنوية في المعادلات ونعتمد ما يلي كتقديرات أولية:

\lambda' = 0 و \cos^2 \alpha_0 = 0.5 و \cos 2\sigma_m = 0 و |\sigma| = \pi - |U_1 + U_2|

ثم تطبق المعادلات التالية أكثر من مرة بشكل تكراري:

\sin \lambda' = \frac{\sin \alpha_0 \sin \sigma}{\cos U_1 \cos U_2}

\cos 2\sigma_m = \cos \sigma - \frac{2 \sin U_1 \sin U_2}{\cos^2 \alpha_0}

\sin^2 \sigma = (\cos U_2 \sin \lambda')^2 + (\cos U_1 \sin U_2 + \sin U_1 \cos U_2 \cos \lambda')^2

C = \frac{1}{16} f \cos^2 \alpha_0 \left[4 + f (4 - 3 \cos^2 \alpha_0)\right]

D = (1 - C) f \left[\sigma + C \sin \sigma \left(\cos 2\sigma_m + C \cos \sigma (-1 + 2 \cos^2 2\sigma_m)\right)\right]

\sin \alpha_0 = \frac{L' - \lambda'}{D}

يتوقف التكرار عندما يصبح التغير في \sin \alpha_0 أصغر من القيمة المحددة سابقًا.

#### شرح تفصيلي لخطوات التكرار
 **تجنبًا للغموض في تطبيق المعادلات السابقة فإن الفكرة الأساسية للتكرار بهذا التسلسل:**

 
 **<big>1 المعطيات</big>**
لدينا نقطتان على الإهليلج:\phi_1 , \lambda_1 ، \phi_2 , \lambda_2

وثابت التفلطح: f
 **<big>2 حساب فرق خط الطول</big>**
نحسب: L = \lambda_2 - \lambda_1 ، ولأننا في حالة شبه متقابلة نحسب: L' = \pi - L
 **<big>3 حساب خطوط العرض المختزلة</big>**
U_1 = \arctan((1-f)\tan\phi_1)

U_2 = \arctan((1-f)\tan\phi_2)
 **<big>4 القيم الابتدائية</big>**
نضع القيم الابتدائية للتكرار: \lambda'_{init} = 0

التقدير الأولي للقوس الجيوديسي: \sigma_{init} = \pi - |U_1 + U_2|

كما نضع تقديرًا أوليًا: \cos^2\alpha_{0\ init} = 0.5

 **<big>5 التكرار الأول</big>**

 
 
 **أ. <math>\sigma_{init}</math>**

 
 
  \sigma_{init} = \pi - |U_1 + U_2|

 
 **ب. حساب <math>D_{new1}</math>**
C_{init}
= \frac{f}{16}
\cos^2\alpha_{0\ init}
(4 + f(4 - 3\cos^2\alpha_{0\ init}))

ثم: D_{init}
= (1 - C_{init})f
\left[
\sigma_{init}
+
C_{init}\sin\sigma_{init}
\right]
 **<big>جـ. حساب <math>\sin\alpha_{0\ new1}</math></big>**
\sin\alpha_{0\ new1}
= \frac{L' - \lambda'_{init}}{D_{init}}

**<big>د. حساب <math>\cos^2\alpha_{0\ new1}</math></big>**

\cos^2\alpha_{0\ new1}
= 1 - \sin^2\alpha_{0\ new1}
 **<big>هـ.شرط التقارب بعد التكرار الأول</big>**
**<big>نحسب: <math>|\sin\alpha_{0\ new1} - \sin\alpha_{0\ init}|</math></big>**

إذا كان أصغر من حد معين (مثلاً 10^{-12}) نتوقف.
وإلا ننتقل للتكرار الثاني.

 **<big>6 التكرار الثاني</big>**

 
 
 **<big>أ. تحديث <math>\lambda'</math></big>**

\sin\lambda'_{new1}
= \frac{
\sin\alpha_{0\ new1}
\sin\sigma_{init}
}{\cos U_1 \cos U_2}

ثم: \lambda'_{new1}
= \arcsin(\sin\lambda'_{new1})
 **<big>ب. حساب <math>\sigma_{new1}</math></big>**
باستخدام: \lambda'_{new1}

\sin^2\sigma_{new1}
= (\cos U_2 \sin\lambda'_{new1})^2
+ (\cos U_1 \sin U_2
+ \sin U_1 \cos U_2
\cos\lambda'_{new1})^2

ثم: \sin\sigma_{new1}
= \sqrt{\sin^2\sigma_{new1}}

ثم: \cos\sigma_{new1}
= \sin U_1 \sin U_2
+ \cos U_1 \cos U_2
\cos\lambda'_{new1}

ثم: \sigma_{new1}
= \operatorname{atan2}
(\sin\sigma_{new1},\cos\sigma_{new1}) ،

**ملاحظة جانبية**: في بعض التطبيقات البرمجية مثل الإكسيل تقدم دالة جيب التمام على دالة الجيب للحصول على نتيجة مباشرة صحيحة لـدالة **<math>
\operatorname{atan2}
</math>**

 **<big>جـ. حساب <math>D_{new1}</math></big>**
C_{new1}
= \frac{f}{16}
\cos^2\alpha_{0\ new1}
(4 + f(4 - 3\cos^2\alpha_{0\ new1}))

D_{new1}
= (1 - C_{new1})f
\left[
\sigma_{new1}
+
C_{new1}\sin\sigma_{new1}
\right]
 **<big>هـ. حساب <math>\sin\alpha_{0\ new2}</math></big>**
\sin\alpha_{0\ new2}
= \frac{L' - \lambda'_{new1}}{D_{new2}}

 **<big>و. حساب <math>\cos^2\alpha_{0\ new2}</math></big>**
\cos^2\alpha_{0\ new2}
= 1 - \sin^2\alpha_{0\ new2}
 **<big>هـ. شرط التقارب بعد التكرار الثاني</big>**
|\sin\alpha_{0\ new2} - \sin\alpha_{0\ new1}|

إذا كان صغيرًا نتوقف.
وإلا ننتقل للتكرار الثالث.

 **<big>7 التكرار الثالث</big>**

 **<big>أ. تحديث <math>\lambda'</math></big>**

\sin\lambda'_{new2}
= \frac{
\sin\alpha_{0\ new2}
\sin\sigma_{new1}
}{\cos U_1 \cos U_2}

ثم: \lambda'_{new2}
= \arcsin(\sin\lambda'_{new2})
 **<big>ب. حساب <math>\sigma_{new2}</math></big>**
باستخدام: \lambda'_{new2}

بنفس المعادلات السابقة.
 **<big>جـ. حساب <math>D_{new2}</math></big>**
باستخدام: \cos^2\alpha_{0\ new2} ، **<big><math>\sigma_{new2}</math> ،</big>** 
 **<big>د. حساب <math>\sin\alpha_{0\ new3}</math></big>**
\sin\alpha_{0\ new3}
= \frac{L' - \lambda'_{new2}}{D_{new2}}
 **<big>هـ. شرط التقارب بعد التكرار الثالث</big>**
|\sin\alpha_{0\ new3} - \sin\alpha_{0\ new2}|

إذا كان الفرق صغيرًا جدًا يتوقف التكرار.
 **<big>و. تحديث <math>\lambda'</math></big>**
\lambda'_{new3}

\sin\lambda'_{new3}
= \frac{
\sin\alpha_{0\ new3}
\sin\sigma_{new2}
}{\cos U_1 \cos U_2}

ثم: \lambda'_{new3}
= \arcsin(\sin\lambda'_{new3})

بعد تحقق شرط التقارب نستخدم آخر \sin\alpha_{0} ونحسب منها السمت الأمامي \alpha_{1} :

\sin \alpha_1 = \frac{\sin \alpha_0}{\cos U_1}

\cos \alpha_1 = \pm\sqrt{1 - \sin^2 \alpha_1}

وتكون إشارة \cos \alpha_1 سالبةً إذا تحقق التالي:

\cos U_1 \sin U_2 + \sin U_1 \cos U_2 \cos \lambda' 

يُحدد السمت عند P_1 باستخدام المعادلات:

1- \sin \alpha_1 = \frac{\sin \alpha_0}{\cos U_1}

2- \cos \alpha_1 = \pm\sqrt{1 - \sin^2 \alpha_1}

3- \tan \alpha_1 = \frac{\sin \alpha_1}{\cos \alpha_1}

وحساب السمت عند النقطة P_2:

\tan \alpha_2 = \frac{\sin \alpha_0}{-\sin U_1 \sin \sigma + \cos U_1 \cos \sigma \cos \alpha_1}

بالمعادلة الأخيرة يُحدد الربع بناءً على إشارات البسط والمقام.

يمكن لمعظم أنظمة الحاسوب تحديد الظل في الربع المناسب بناءً على هذه الأسس وبالتالي ليس ضروريًا القيام بذلك من خلال تعليمات البرمجة.

يتم حساب المسافة الجيوديسية باستخدام **صيغة هلمرت (Helmert, 1880)**، وتمت إعادة كتابتها هنا بصياغة أكثر ملاءمة للاستخدام مع أجهزة الكمبيوتر:

E = \sqrt{1 + \epsilon \cos^2 \alpha_0}

F = \frac{E - 1}{E + 1}

A = \frac{1 + \frac{1}{4}F^2}{1 - F}

B = F \left(1 - \frac{3}{8}F^2\right)

\Delta\sigma = B \sin \sigma \left[\cos 2\sigma_m + \frac{1}{4}B \left[\cos \sigma(-1 + 2\cos^2 2\sigma_m)\right] - \frac{1}{6}B \cos 2\sigma_m (-3 + 4\sin^2 \sigma)(-3 + 4\cos^2 2\sigma_m)\right]

حساب المسافة الجيوديسية s:

s = b A (\sigma - \Delta\sigma) = (1 - f) a A (\sigma - \Delta\sigma)

### '''الحالة الخاصة للنقاط المتقابلة تمامًا'''
(\phi_2 = -\phi_1)، لا تتطلب هذه الحالة استخدام طرق تكرارية للحل

في هذه الحالة تكون المسافة الزاوية \sigma = \pi وفرق الطول المساعد \lambda' = 0

لذا لدينا من المعادلة: \sin \alpha_0 = \frac{L' - \lambda'}{D}

\frac{L'}{f \pi} = \sin \alpha_0 (1 - C) = Q

بالتعويض عن \cos^2 \alpha_0 = 1 - \sin^2 \alpha_0، يمكن كتابة المعادلة الأخيرة كمتسلسلة:

Q = a_1 \sin \alpha_0 + a_3 \sin^3 \alpha_0 + a_5 \sin^5 \alpha_0

حيث:

a_1 = 1 - \frac{f}{4} - \frac{f^2}{16}

a_3 = \frac{f}{4} - \frac{f^2}{8}

a_5 = \frac{3f^2}{16}

بعكس المتسلسلة نحصل على المعادلة: \sin \alpha_0 = b_1 Q + b_3 Q^3 + b_5 Q^5

حيث:

b_1 = 1 + \frac{f}{4} + \frac{f^2}{8}

b_3 = 1 - b_1

b_5 = 0

حساب المسافة الجيوديسية بواسطة المعادلة: s = b A \pi

حيث أن: b = (1 - f) a

A = \frac{1 + \frac{1}{4}F^2}{1 - F}

F = \frac{E - 1}{E + 1}

E = \sqrt{1 + \epsilon \cos^2 \alpha_0}

\cos^2 \alpha_0 = 1 - \sin^2 \alpha_0

السمتان \alpha_1 و \alpha_2 يحددان باستخدام المعادلات التالية:

\sin \alpha_1 = \frac{\sin \alpha_0}{\cos U_1}

\cos \alpha_1 = \pm\sqrt{1 - \sin^2 \alpha_1}

\tan \alpha_1 = \frac{\sin \alpha_1}{\cos \alpha_1}

\tan \alpha_2 = \frac{\sin \alpha_0}{-\sin U_1 \sin \sigma + \cos U_1 \cos \sigma \cos \alpha_1}

هذا يوضح أنه في الحالات الخاصة من التقابل عندما تكون \phi_2 = -\phi_1 فإن قيم السمت \alpha_0 والمسافة الجيوديسية s لا تعتمد على خط العرض بل تعتمد فقط على فرق الطول بين النقطتين.

**ملاحظة**: في حالات التقابل التام يكون أقصر مسار بين النقطتين عبر القطبين بمعنى أن كل من \alpha_0 و Q = صفرًا ودائمًا يكون السمتان الأمامي والخلفي \alpha_1 و \alpha_2 عند النقطتين P_1 و P_2 إلى الشمال وإلى الجنوب فلو كانت النقطة الأولى P_1 شمال النقطة الثانية P_2 أي أن (\phi_1 > \phi_2) يكون السمت الأول \alpha_1 إلى الجنوب تمامًا أي 180 درجة، وبالتالي يكون السمت عند النقطة الثانية \alpha_2 إلى الشمال تمامًا أي صفر درجة والعكس صحيح.

#### '''تحليل الاستقرارات'''
قد تنشأ أسئلة حول استقرار المعادلات التالية:

\cos 2\sigma_m = \cos \sigma - \frac{2 \sin U_1 \sin U_2}{\cos^2 \alpha_0}
\sin \alpha_0 = \frac{L' - \lambda'}{D}
\cos \alpha_1 = \pm\sqrt{1 - \sin^2 \alpha_1}

وفيما يلي التحقق من ذلك:

المعادلة: \cos 2\sigma_m = \cos \sigma - \frac{2 \sin U_1 \sin U_2}{\cos^2 \alpha_0}

الحد الثاني فيها لا يمكن أن يصبح غير محدد، لعدم وجود خطوط استوائية في حلول النقاط المتقابلة.

إذا كان \cos \alpha_0 \approx 0 قريبًا جدًا من الصفر فأي خطأ في \cos^2 \alpha_0 يصبح غير مهم لأنه سيختفي عند ضربه بالقيمة الصغيرة جدًا لـ f^2 \cos^2 \alpha_0/4 وستعطي المعادلة:

D = (1 - C) f \left[\sigma + C \sin \sigma \left(\cos 2\sigma_m + C \cos \sigma (-1 + 2 \cos^2 2\sigma_m)\right)\right]

نتيجةً صحيحةً.

(في الحالات غير المتقابلة، يكفي استبعاد القسمة على الصفر مثلًا، عن طريق تعيين هذا الحد مسبقًا إلى الصفر واستبداله بقيمة محسوبة عندما تكون \cos^2 \alpha_0 غير صفرية)

### '''ملاحظة:'''
تكون قيمة \cos^2 \alpha_0 = 0 عندما تكون قيمة \alpha_0 مساوية لـ 90 درجة وحينها تكون النقطتان P_1 و P_2 على خط الاستواء ويكون السمتان الأمامي والخلفي \alpha_1 و \alpha_2 إلى الشرق وإلى الغرب على حسب موقع كل نقطة من الأخرى فلو كانت النقطة الأولى P_1 شرق النقطة الثانية P_2 أي أن (\lambda_1 > \lambda_2) سيكون سمت النقطة الأولى \alpha_1 للغرب 270 درجة ويكون سمت النقطة الثانية \alpha_2 للشرق أي 90 درجة والعكس صحيح.

**ملاحظة عملية لتجنب القسمة على صفر، من خارج الورقة ب لفينسنتي:**

"من المزايا التي يمكن استغلالها في دالة \operatorname{atan2}(\sin, \cos) التغلب على نتيجة القسمة على الصفر في الدالة التقليدية \tan^{}\left(\frac{\sin}{\cos}\right) وما يشبِهُها وذلك بوضع قيمة المقام كحد أول في الدالة ووضع قيمة البسط كحد ثاني فلو رمزنا للمقام بالرمز y والبسط بالرمز x نكتب المعادلة هكذا \operatorname{atan2}(y, x) وبعد ذلك نحسب ظل القيمة الناتجة فتصير المعادلة هكذا \tan(\operatorname{atan2}(y, x)).

ومثال رقمي على هذه الطريقة لو كانت قيمة البسط = 3 وقيمة المقام = 0 فإننا نكتب المعادلة هكذا \tan(\operatorname{atan2}(0, 3)) والنتيجة مُعرَّفة 1.63246E+16
 من برنامج إكسيل

ولكن الدالة \operatorname{atan2}(0, 3) تعطي نتيجةً مُعرَّفةً في حالة \frac{\pi}{2} وبالتالي يمكن وضع أي معادلة مِن بسطٍ ومقامٍ جائز أن يكون صفرًا بهذه الصورة \tan(\operatorname{atan2}(y, x)) لتجنب النتيجة غير المُعرّفة من القسمة على صفر ."

المعادلة: \sin \alpha_0 = \frac{L' - \lambda'}{D}

يمكن كتابتها بشكل مكافئ رياضيًا كالتالي: \sin \alpha_0 = \frac{\lambda - L}{D}

لكن هذا قد يؤدي إلى فقدان أرقام معنوية في النتيجة، لأن أول رقمين على الأقل في قيمتي \lambda و L متشابهان، وبالتالي سيتم فقدان نفس العدد من الأرقام في عملية الطرح.

بالتالي لو حوّلنا L إلى L'، ستكون القيمة الجديدة صحيحة بنفس عدد المنازل العشرية المستخدمة في L.

ثم \lambda'، كما تم حسابه مِن المعادلة: \sin \lambda' = \frac{\sin \alpha_0 \sin \sigma}{\cos U_1 \cos U_2}

سيكون مُعطَى بمزيدٍ مِن الأرقام المعنوية والمنازل العشرية مقارنة بـ L'، وسيتم كسب رقمين معنويين أو أكثر في \sin \alpha_0.

أما المعادلة: \cos \alpha_1 = \pm\sqrt{1 - \sin^2 \alpha_1}

فتكون غير مستقرة عندما يقترب \alpha_1 من \pm\frac{\pi}{2}

هذا يمكن أن يحدث فقط إذا كانت P_2 تتطابق تقريبًا مع P_a؛\phi_a = -\phi_1،

L_a محدد بمعادلة L_a = \pi - f \cos U \, \pi \left\{1 - \frac{1}{16} f \sin^2 U \left[4 + f (4 - 3 \sin^2 U)\right]\right\}

في هذه الحالة، أي تغيير طفيف في موقع P_2 يسبب تغييرًا كبيرًا في \alpha_1

بمعنى أن أي خطأ في \alpha_1 بسبب فقدان الدقة في حساب جيب التمام باستخدام الصيغة الفيثاغورسية لن يسبب خطأً في الموضع أكبر مما يمكن توقعه في أي حالة أخرى، أي أن الخطأ ناتج فقط عن العدد المحدود من الأرقام التي يمكن للحاسوب معالجتها.

لذلك، عدم استقرار المعادلة: \cos \alpha_1 = \pm\sqrt{1 - \sin^2 \alpha_1} غير ضار عمليًا.

## '''الحل قرب الجانب الخلفي للإهليج'''
يبدأ الخط الجيوديسي من النقطة P_1 عند خط عرض \phi_1 وسمت \alpha_1 ويصل إلى خط العرض

\phi_a = -\phi_1 عند P_a في السمت حيث يكون \sin \alpha_a = \sin \alpha_1 و \cos \alpha_a = -\cos \alpha_1

عند هذه النقطة \sigma_a = \lambda_a = \pi وخط الطول الجيوديسي يُعطى بالمعادلة:

= \pi [1 - f (1 - C) \cos U_1 \sin \alpha_1] L_a = \pi [1 - f (1 - C) \sin \alpha_0]

عند تمديد الخط بسمت \alpha_a مسافة جيوديسية P_1 P_a فإنه سيصل إلى النقطة P_t حيث يكون

\phi_t = \phi_1 و L_t = 2L_a و \alpha_t = \alpha_1 و \sigma_t = \lambda_t = 2\pi

وإذا تم اختيار P_2 على هذا الخط بحيث تكون المسافة الكروية P_1 P_a P_2 أكبر من \pi فإن الخط الجيوديسي P_1 P_a P_2 ليس خطًا جيوديسيًا ولكنه خط حول الجانب الخلفي من القطع الناقص.

وأقصر مسافة P_1 P_2 هي خط جيوديسي يبدأ بالسمت وعلى العموم هو مختلف عن \alpha_1 \pm \pi

يمكن حساب طول الخط الجيوديسي حول الجانب الخلفي للقطع الناقص والسموت عند نهايته باستخدام:

الصيغة العكسية القياسية (**Vincenty 1975**) بعد تعديلها بجعل \sin \sigma تعطي إشارة سالبة (2\pi \ge \sigma > \pi) وتتغير زاويتا السمت بمقدار \pm \pi

تفشل هذه الطريقة عندما تكون النقاط قريبة من بعضها أي عندما تكون داخل منطقة مجاورة معينة، بسبب تذبذب فرق الطول \lambda بين قيم موجبة وسالبة في التكرارات، وحينئذٍ يتم الحل باستخدام المعادلات المُقدَّمة في هذه الورقة مع بعض التعديلات كما هو الحال في حالة النقاط المتقابلة، فإن الهدف هو الحصول على \sin \alpha_0 عن طريق تطبيق المعادلات:

C = \frac{1}{16} f \cos^2 \alpha_0 \left[4 + f (4 - 3 \cos^2 \alpha_0)\right]

\cos 2\sigma_m = \cos \sigma - \frac{2 \sin U_1 \sin U_2}{\cos^2 \alpha_0}

D = (1 - C) f \left[\sigma + C \sin \sigma \left(\cos 2\sigma_m + C \cos \sigma (-1 + 2 \cos^2 2\sigma_m)\right)\right]

\sin \alpha_0 = \frac{L' - \lambda'}{D}

\sin \lambda' = \frac{\sin \alpha_0 \sin \sigma}{\cos U_1 \cos U_2}

\sin^2 \sigma = (\cos U_2 \sin \lambda')^2 + (\cos U_1 \sin U_2 + \sin U_1 \cos U_2 \cos \lambda')^2

ويتوقف التكرار عندما يصبح التغير في \sin \alpha_0 أصغر من القيمة المحددة سابقًا.

مع استخدام المعادلة \sin \alpha_0 = \frac{\lambda - L}{D}

بدلًا من \sin \alpha_0 = \frac{L' - \lambda'}{D}

لأن L و \lambda سيكونان أقرب إلى \pi والحصول على قيمة سيجما \sigma من المعادلة:

\sin^2 \sigma = (\cos U_1 \sin \lambda)^2 + (\cos U_2 \sin U_2 - \sin U_1 \cos U_2 \cos \lambda)^2

وتصبح إشارة \cos \alpha_1 سالبة عندما: \cos U_1 \sin U_2 - \cos U_1 \cos U_2 \cos \lambda > 0

يمكن استخدام تقريبات ابتدائية لهذا الحلِّ على النحو التالي:

\lambda = 0 ؛ \cos^2 \alpha_0 = 0.5 ؛ \cos 2\sigma_m = 0 ؛ \sigma = 2\pi - |U_1 - U_2|

بشكل عام فَمِن الملاحظ حسب **هيلمرت (1880)** يوجد عدد لا نهائي من الخطوط الجيوديسية التي تصل بين أي نقطتين، ربما يمكن تحديد عدد اللفات حول القطع الناقص بشكل خاص عند: 2n\pi \ge \sigma > (2n-1)\pi

حيث n عدد اللفات الكاملة حول القطع الناقص، في هذا السياق نستخدم n=1 لفة واحدة حول الجانب الخلفي؛

ويمكن التحقق من صحة الحل العكسي حول الجانب الخلفي للإهليج باستخدام أي صيغة مباشرة قياسية.

من المثير للاهتمام ملاحظة أنه في الحل المباشر لا يلزم تقسيم الخط إلى قسمين أو أكثر بحيث يلبي كل منهما الشرط \sigma 

لا يوجد حدٌّ نظريٌّ لطول الخطِّ الذي يمكن استخدامه في الحلِّ المباشر، القيود الوحيدة هي رقمية (تتعلق بدقة الحسابات العددية).

أمثلة حسابية

تم تقديم أمثلة على نتائج **الحلول المضادة (Antipodal)** وحلول الجانب الخلفي في الجدولين 1، 2

تم تقريب السموت إلى 3 منازل عشرية فقط مِن الثانية، حيث أنَّ هذه الدِّقَّة كافية لإعادة إنتاج المواضع بدقة تصل إلى 0.5 مم عند استخدام الحل المباشر،

هذا ينطبق بشكل خاص عندما تكون \sigma تساوي \pi أو 2\pi أو قريبة من تلك القيمة.

**الجدول 1**: **حلول النقاط المتقابلة (Antipodal Solutions)**

 **تم حسابه باستخدام: لغة برمجة**  **FORTRAN** **، حاسوب** **IBM 7040/94** 

 **الجدول 2**: **حلول الجانب الخلفي (Back Side Solutions)**

تم حسابه بواسطة برنامج الآلة الحاسبة المكتبية **Wang 720**

معلمات القطع الناقص المستخدمة: المحور الأكبر **a = 6,378,388.000** مترًا ، ومقلوب التفلطح **(1/f) = 297**

## مصادر

## مراجع

*  English translation of Astron. Nachr. **4**, 241&ndash;254 (1825).
*  English translation of *Die Mathematischen und Physikalischen Theorieen der Höheren Geodäsie*, Vol. 1 (Teubner, Leipzig, 1880).
*  Addenda.
* 
* 
* 
* 
* 
* 
* 

## وصلات خارجية
* Online calculators from Geoscience Australia:
** Vincenty Direct (destination point)
** Vincenty Inverse (distance between points)
* Calculators from the U.S. National Geodetic Survey:
** Online and downloadable PC-executable calculation utilities, including forward (direct) and inverse problems, in both two and three dimensions (accessed 2011-08-01).
* Online calculators with JavaScript source code by Chris Veness (Creative Commons Attribution license):
** Vincenty Direct (destination point)
** Vincenty Inverse (distance between points)
* GeographicLib provides a utility GeodSolve (with MIT/X11 licensed source code) for solving direct and inverse geodesic problems. Compared to Vincenty, this is about 1000 times more accurate (error = 15&nbsp;nm) and the inverse solution is complete. Here is an online version of GeodSolve.
* Complete Vincenty's direct and inverse formulae implementation with source code, Excel VBA implementation by Tomasz Jastrzębski

تصنيف:مسح الأرض
تصنيف:مسافة