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

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

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

 خلفية 
كان هدف فينسينتي هو التعبير عن الخوارزميات الحالية للجيوديسيات على شكل إهليجي بصيغة تُقلل من طول البرنامج (فينسينتي 1975). يذكر تقريره غير المنشور (فينسنتي 1975) استخدام حاسبة مكتبية من طراز وانغ 720، والتي كانت تحتوي على كيلوبايتات قليلة من الذاكرة. للحصول على دقة جيدة للخطوط الطويلة، استخدم فينسنتي الحل الكلاسيكي الذي وضعه ليجاندر (1806) وبيسل (1825) وهيلمرت (1880) والمبني على الكرة المساعدة. اعتمد فينسينتي على صياغة هذه الطريقة التي قدمها رينسفورد (1955)؛ أظهر ليجاندر أنه يمكن رسم خريطة دقيقة لجيوديسية إهليجية على دائرة عظمى على الكرة المساعدة عن طريق رسم خريطة لخط العرض الجغرافي على خط عرض مُخفض، وضبط سمت الدائرة العظمى على أنه يساوي سمت الجيوديسية. يُعطى بعد ذلك خط الطول على الشكل الإهليجي والمسافة على طول الخط الجيوديسي بدلالة خط الطول على الكرة وطول القوس على طول الدائرة العظمى باستخدام تكاملات بسيطة. وقد قدّم بيسل وهيلمرت سلاسل متقاربة سريعة لهذه التكاملات، مما يسمح بحساب الخط الجيوديسي بدقة متناهية.

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

 الرموز والمصطلحات 
الرموز المستخدمةتعريفقيمةaطول المحور الأكبر للقطع الناقص (نصف القطر عند خط الاستواء).(6378137.0 مترًا في WGS-84 )ƒتسطيح القطع الناقص.(1/298.257223563 في WGS-84 )طول المحور الأصغر للقطع الناقص (نصف القطر عند القطبين).(6356752.314245 مترًا في WGS-84)φ₁ ، φ₂خطا عرض الموقعين.خطا العرض المنخفضين (خط العرض على الكرة المساعدة).L₁, L₂خطا الطول للموقعين.الفرق بين خطيّ الطول.λالفرق بين خطيّ الطول بين على الكرة المساعدة.α₁ ، α₂السمتان الأمامي والخلفي بداية ونهاية المسار بين الموقعين.α₀السمت الأمامي للخط الجيوديسي عند خط الاستواء، إذا تم تمديده إلى هذا الحد؛sالمسافة الإهليلجية بين النقطتين؛σالفصل الزاوي بين الموقعين وهي المسافة بالدرجات القوسيةσ₁الفصل الزاوي بين الموقع الأول وخط الاستواءσₘالمسافة الزاوية بين منتصف الخط وخط الاستواء

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

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





 يتم تجنب حساب σ (سيجما) مباشرةً مِن دالة جيب σ أو جيب تمام σ للحفاظ على الدقة العددية بالقرب مِن القطبين وخط الاستواء.


 إذا كانت المسافة الزاوية σ = 0 فإن قيمة جيب α₀ تصبح غير مُعَرَّفة (غير محددة). وتُمثل هذه الحالة نقطة نهاية مُطابِقة لنقطة البداية، أو مُعاكِسة لها قُطريًا (على الطرف الآخر مِن القطر).




 عندما تقع نقطتا البداية والنهاية على خط الاستواء، فإن C = 0 ولا يتم استخدام قيمة جيب تمام (2σₘ). في هذه الحالة، تكون القيمة الحدية لجيب تمام (2σₘ) = -1.




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



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

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

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

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

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

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

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

والجدير بالذكر أن البعض يفضل تقديم العنصر الثاني في الدالة على العنصر الأول للحصول على نتيجة صحيحة ومباشرة بهذا الشكل (x,y)atan2(cos,sin) ، atan2 بينما الصورة التقليدية هي (y,x)atan(sin,cos) ، atan بمعنى أن عند التطبيق في برنامج الإكسيل وما يشابهه يفضل عكس الترتيب السابق في المعادلات أعلاه للحصول على نتيجة مباشرة وبدون إشارات سالبة.

 تعديل فينسنتي 
في رسالته إلى مجلة Survey Review في عام 1976، اقترح فينسينتي استبدال تعبيراته المتسلسلة لـ A و B بصيغ أبسط باستخدام معامل التوسع لـ Helmert k₁  :



 حيث


 نقاط متقابلة تقريبًا 
كما هو مذكور أعلاه، فإن الحل التكراري بالصيغة القياسية للمشكلة العكسية يفشل في التقارب أو يتقارب ببطء بالنسبة للنقاط المتقابلة تقريبًا. ومن الأمثلة على التقارب البطيء ( φ₁ ، 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₂ - L₁│ = │-91 - 91│ = 182│)؛

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

أما طريقة سايتو (1970)، والتي صُممت للاستخدام العام، فهي صالحة فقط عندما تكون نقاط نهاية الخط في خطوط عرض متقابلة (φ₂ = −φ₁) ؛

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

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

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

التفلطح: f = (a − b) / a

مربع الانحراف الثاني: ϵ = (a² − b²) / b²

ويرمز للرمز ϵ بالرمز e'² أيضًا وهذا ما نجده مثلا في معادلات كارني، حيث أن  e² = (a² − b²) / a² ،

خط العرض الجيوديسي (موجب شمال خط الاستواء): φ₁

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

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

السمت الجيوديسي باتجاه عقارب الساعة من الشمال؛ α₂ في اتجاه امتداد (P₁P₂): α₁, α₂

سمت الخط الجيوديسي عند خط الاستواء: α₀

خط العرض المخفض: U يُعرّف بـ tan U = (1 − f) * tan φ

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

المسافة الزاوية بين P₁ و P₂ على الكرة: σ

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

 نحن ندرس خطاً جيوديسياً ينطلق من النقطة P₁ في نصف الكرة الشمالي عند خط عرض φ₁ وسمت ابتدائي α₁ = π/2 أثناء تتبع هذا الخط، يتغير سمت الخط وفقاً لمعادلة كليرو: Cos Uᵢ * sin αᵢ = ثابت

وفي هذه الحالة الخاصة:

 Sin αᵢ = cos U₁ / cos Uᵢ = 1 

الخط الجيوديسي سيعبر خط الاستواء بسمت α₀ محدد بالعلاقة: Sin α₀ = cos U₁

حتى يصل إلى المنطقة المقابلة عند النقطة Pₐ حيث: φₐ = -φ₁ (خط عرض معاكس)

λ = σ = π (فرق الطول والمسافة الزاوية = 180°)

(السمت النهائي) αₐ = π/2

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

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

Lₐ = π - f * cos U * π * {1 - (1 / 16) * f * sin² U*[ 4 + f * (4 - 3 * sin² U ) ] }

حيث: U₂-│ = U₁ = U│

هذا أكبر خط طول جيوديسي بنهاية الخط الجيوديسي الذي يبدأ بـ α₁ = π/2

الآن سنختار النقطة P₂ بحيث ( φ₂ = φₐ أو φ₂ ≈ φₐ) و π > L > Lₐ ،

الخط الجيوديسي الأقصر بين P₁ و P₂ لابد أن يبدأ بسمت مختلف عن (π/2) بسبب تفلطح الإهليج وهناك ثلاث احتمالات:

إذا كان φ₁ + φ₂ > 0 فإن: π/2 > α₁ > 0

إذا كان φ₁ + φ₂ < 0 فإن: π > α₁ > π/2

  إذا كان φ₁ + φ₂ = 0  فإن :  λ = σ = π

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

 إذا تتبعنا خطًّا يبدأ من φ₁ > 0 وبسمت π/2 > α₁ > 0 فهذا المسار سيبلغ خط عرضه الأكبر

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

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

إذا قمنا بتمديد هذا المسار إلى خط عرض φₐ واخترنا أي نقطة P₂ عليه بحيث φ₂ ≥ φₐ، نجد من الحل المباشر أن (λ < π) فرق الطول على الكرة المساعدة أقل من 180° و (σ < π) المسافة الزاوية أقل من 180°، في الحل العكسي λ يتذبذب بين قيم أكبر وأصغر من π بدون تقارب إلى أي قيمة.

 يمكن استنتاج وجود منطقة تعتبر أي نقطة فيها "مضادة" بالنسبة لـ P₁ حتى لو لم يكن هذا صحيحًا هندسياً، في نصف الكرة الشرقي المقابل لـ P₁ يجب أن تمتد طولياً من Lₐ إلى  π كما هو موضح بالمعادلة:

Lₐ = π - f * cos U * π * {1 - (1 / 16) * f * sin² U*[ 4 + f * (4 - 3 * sin² U ) ] }

، ويجب أن يحددها منحنيان متماثلان حول φₐ متوافقين مع

 π/2 > α₁ > 0 و π > α₁ > π/2 

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

 حل عكسي 
أثناء الحسابات الأولية إذا كانت λ أكبر من π كقيمة مطلقة، يجب استخدام طريقة مختلفة؛

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

نستخدم L' = π – L ، عندما L > 0

أو L' = -π - L ، عندما L < 0

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

 λ' = 0 و cos² α₀ = 0.5 و cos 2σₘ = 0 و │σ = π -│U₁ + U₂
ثم تطبق المعادلات التالية أكثر من مرة بشكل تكراري

C = (1/16) * f * cos² α₀ * [4 + f * (4 - 3 * cos² α₀)]

Cos 2σₘ = cos σ - (2 * sin U₁ * sin U₂) / cos²α₀

{D = (1 - C) * f * {σ + C * sin σ * [cos 2σₘ + C * cos σ * (-1 + 2 * cos² 2σₘ)]

Sin α₀ = (L' - λ') / D

sin λ' = (sin α₀ * sin σ) / (cos U₁ * cos U₂)

sin² σ = (cos U₂ * sin λ')² + (cos U₁ * sin U₂ + sin U₁ * cos U₂ * cos λ')²

يتوقف التكرار عندما يصبح التغير في sin α₀ أصغر مِن القيمة المحددة سابقا.

نحصل الآن على:

 Sin α₁ = sin α₀ / cos U₁  

Cos α₁ = ±√(1 - sin² α₁)

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

Cos U₁ * sin U₂ + sin U₁ * cos U₂ * cos λ' < 0

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

1-   Sin α₁ = sin α₀ / cos U₁ 

2-  Cos α₁ = ±√(1 - sin² α₁) 

3-  Tan α₁ = sin α₁ / cos α₁

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

Tan α₂ = sin α₀ / (-sin U₁ * sin σ + cos U₁ * cos σ * cos α₁)

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

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

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

E = √(1 + ϵ * cos² α₀)

F = (E - 1) / (E + 1)

A = (1 + ¼F²) / (1 - F)

B = F * (1 - ⅜F²)

Δσ = B * sin σ * [cos 2σₘ+ ¼B * [cos σ * (-1 + 2 * cos² 2σₘ)] -⅙B * cos 2σₘ* (-3 + 4 * sin² σ) * (-3 + 4 * cos² 2σₘ)]

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

s = b * A * (σ - Δσ) = (1 - f) * a * A* (σ - Δσ)

 الحالة الخاصة للنقاط المتقابلة تماماً 
 (φ₂ = -φ₁) ، لا تتطلب هذه الحالة استخدام طرق تكرارية للحل

في هذه الحالة تكون المسافة الزاوية σ = π وفرق الطول المساعد λ' = 0

لذا لدينا من المعادلة:

 Sin α₀ = (L' - λ') / D  

L′ / (f * π) = sin α₀ * (1 - C) = Q

بالتعويض عن cos² α₀ = 1 - sin² α₀ ، يمكن كتابة المعادلة الأخيرة كمتسلسلة:

Q = a₁ * sin α₀ + a₃ * sin³ α₀ + a₅ * sin⁵ α₀

حيث:

a₁ = 1 - f / 4 - f² / 16

a₃ = f / 4 - f²/ 8

a₅ = 3 * f² / 16

بعكس المتسلسلة نحصل على المعادلة:

Sin α₀ = b₁ * Q + b₃ * Q³ + b₅ * Q⁵

حيث:

b₁ = 1 + f / 4 + f² / 8

b₃ = 1 - b₁

b₅ = 0

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

s = b * A * π 

حيث أنَّ

b = (1 - f) * a

A = (1 + ¼F²) / (1 - F)

F = (E - 1) / (E + 1)

E = √(1 + ϵ * cos² α₀)

cos² α₀ = 1 - sin² α₀

السمتان α₁ و α₂ يحددان باستخدام المعادلات التالية:

Sin α₁ = sin α₀ / cos U₁

Cos α₁ = ±√(1 - sin² α₁)

Tan α₁ = sin α₁ / cos α₁

Tan α₂ = sin α₀ / (-sin U₁ * sin σ + cos U₁ * cos σ * cos α₁)

هذا يوضح أنه في الحالات الخاصة مِن التقابل عندما تكون φ₂ = -φ₁ فإن قيم السمت α₀ والمسافة الجيوديسية s لا تعتمد على خط العرض بل تعتمد فقط على فرق الطول بين النقطتين.

ملاحظة: في حالات التقابل التامّ يكون أقصر مسار بين النقطتين عبر القطبين بمعنى أن كل من α₀ و Q = صفرًا ودائما يكون السمتان الأمامي والخلفي α₁ و α₂ عند النقطتين P₁ و P₂ إلى الشمال وإلى الجنوب فلو كانت النقطة الأولى P₁ شمال النقطة الثانية P₂ أي أنَّ (φ₁ > φ₂) يكون السمت الأول α₁ إلى الجنوب تمامًا أي 180 درجة، وبالتالي يكون السمت عند النقطة الثانية α₂ إلى الشمال تمامًا أي صفر درجة والعكس صحيح.

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

 Cos 2σₘ = cos σ - (2 * sin U₁ * sin U₂) / cos²α₀
 Sin α₀ = (L' - λ') / D  

Cos α₁ = ±√(1 - sin² α₁)

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

 المعادلة:

 Cos 2σₘ = cos σ - (2 * sin U₁ * sin U₂) / cos²α₀ 

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

إذا كان cos α₀ ≈ 0 قريبًا جدًا من الصفر فأي خطأ في cos² α₀ يصبح غير مهم لأنه سيختفي عند ضربه بالقيمة الصغيرة جدًا لـ f² * cos² α₀/4 وستعطي المعادلة:

 D = (1 - C) * f * [σ + C * sin σ * (cos 2σₘ + C * cos σ * (-1 + 2 * cos² 2σₘ)]

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

( في الحالات غير المتقابلة، يكفي استبعاد القسمة على الصفر مثلاً، عن طريق تعيين هذا الحدّ مسبقًا إلى الصفر واستبداله بقيمة محسوبة عندما تكون cos² α₀ غير صفرية)

ملاحظة: 

تكون قيمة cos² α₀ = صفرًا عندما تكون قيمة α₀ مساوية لـ 90 درجة وحينها تكون النقطتان P₁ و P₂ على الاستواء ويكون السمتان الأمامي والخلفي α₁ و α₂ إلى الشرق وإلى الغرب على حسب موقع كل نقطة من الأخرى فلو كانت النقطة الأولى P₁ شرق النقطة الثانية P₂ أي أنَّ (λ₁ > λ₂) سيكون سمت النقطة الأولى α₁ للغرب 270 درجة ويكون سمت النقطة الثانية α₂ للشرق أي 90 درجة والعكس صحيح.

هذا ومِن المزايا التي يمكن استغلالها في دالة atan2(cos,sin) التغلب على نتيجة القسمة على الصفر في الدالة التقليدية atan(sin,cos) وما يشابهها وذلك بوضع قيمة المقام كحد أول في الدالة ووضع قيمة البسط كحد ثاني فلو رمزنا للمقام والبسط بالرمزين x و y نكتب المعادلة هكذا atan2(x,y) وبعد ذلك نحسب ظل الناتج النهائي فتصير المعادلة هكذا (atan2(x,y))tan ومثال رقمي على هذه الطريقة لو كانت قيمة البسط = 3 وقيمة المقام = 0 تكتب هكذا (atan2(0,3))tan والنتيجة 16324552277619100وبالتالي يمكن وضع أي معادلة من بسط ومقام بهذه الصورة لتجنب النتيجة الغير معرَّفة من القسمة على صفر.

المعادلة:

 Sin α₀ = (L' - λ') / D    

يمكن كتابتها بشكل مكافئ رياضيًا كالتالي:

Sin α₀ = (λ - L) / D

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

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

ثم λ' ، كما تم حسابه بالمعادلة:

 sin λ' = (sin α₀ * sin σ) / (cos U₁ * cos U₂)  

 سيكون معطى بمزيد من الأرقام المعنوية والمنازل العشرية مقارنة بـ L' ، وسيتم كسب رقمين معنويين أو أكثر في sin α₀ .

أمَّا المعادلة:

 Cos α₁ = ±√(1 - sin² α₁) 

فتكون غير مستقرة عندما يقترب α₁ من ±(π/2)

هذا يمكن أن يحدث فقط إذا كانت P₂ تتطابق تقريبًا مع Pₐ ؛

 φₐ = -φ₁  ،  Lₐ محدد بمعادلة 

Lₐ = π - f * cos U * π * {1 - (1 / 16) * f * sin² U * [ 4 + f * (4 - 3 * sin² U ) ] }

في هذه الحالة، أي تغيير طفيف في موقع P₂ يسبب تغييرًا كبيرًا في α₁

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

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

لذلك، عدم استقرار المعادلة :

 Cos α₁ = ±√(1 - sin² α₁)  

غير ضار عمليًا.

 الحل قرب الجانب الخلفي للإهليج 
يبدأ الخط الجيوديسي من النقطة P₁ عند خط عرض φ₁ وسمت α₁ ويصل إلى خط العرض

φₐ = -φ₁ عند Pₐ في السمت حيث يكون Sin αₐ = sin α₁ و Cos αₐ = -cos α₁

عند هذه النقطة σₐ = λₐ = π وخط الطول الجيوديسي يُعطى بالمعادلة:

Lₐ = π * [1 – f * (1 - C) * sin α₀]

= π * [1 – f * (1 - C) * cos U₁ * sin α₁]

عند تمديد الخط بسمت αₐ مسافة جيوديسية P₁Pₐ فإنه سيصل إلى النقطة Pₜ حيث يكون

φₜ = φ₁  و Lₜ = 2Lₐ و αₜ = α₁ و σₜ = λₜ = 2π

وإذا تم اختيار P₂ على هذا الخط بحيث تكون المسافة الكرويةP₁PₐP₂ أكبر من π فإن الخط الجيوديسي P₁PₐP₂ ليس خطًا جيوديسيًا ولكنه خط حول الجانب الخلفي من القطع الناقص.

وأقصر مسافة P₁P₂ هي خط جيوديسي يبدأ بالسمت وعلى العموم هو مختلف عن α₁ ± π

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

الصيغة العكسية القياسية Vincenty 1975 بعد تعديلها بجعل sin σ تعطي إشارةً سالبة

(2π ≥ σ > π)

وتتغير زاوايتا السَّمت بمقدار ±π

تفشل هذه الطريقة عندما تكون النقاط قريبةً من بعضها أي عندما تكون داخل منطقة مجاورة معينة، بسبب تذبذب فرق الطول λ بين قيم موجبة وسالبةفي التكرارات، وحينئذٍ يتم الحل باستخدام المعادلات المقدمة في هذه الورقة مع بعض التعديلات كما هو الحال في حالة النقاط المتقابلة، فإن الهدف هو الحصول على sin α₀ عن طريق تطبيق المعادلات:

C = (1/16) * f * cos² α₀ * [4 + f * (4 - 3 * cos² α₀)]

Cos 2σₘ = cos σ - (2 * sin U₁ * sin U₂) / cos²α₀

D = (1 - C) * f * [σ + C * sin σ * (cos 2σₘ + C * cos σ * (-1 + 2 * cos² 2σₘ)]

Sin α₀ = (L' - λ') / D

sin λ' = (sin α₀ * sin σ) / (cos U₁ * cos U₂)

sin² σ = (cos U₂ * sin λ')² + (cos U₁ * sin U₂ + sin U₁ * cos U₂ * cos λ')²

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

مع استخدام المعادلة

 Sin α₀ = (λ - L) / D 

بدلًا مِن

 Sin α₀ = (L' - λ') / D  

لأن L و λ سيكونان أقرب إلى π والحصول على قيمة سيجما σ من المعادلة:

sin² σ = (cos U₁*sin λ)² + (cos U₂*sin U₂ - sin U₁*cos U₂*cos λ)²

وتصبح إشارة cos α₁ سالبةً عندما:

Cos U₁* sin U₂ - cos U₁* cos U₂* cos λ > 0

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

λ = 0  ؛  cos² α₀ = 0.5  ؛  cos (2σₘ) = 0  ؛  σ = 2π - |U₁ - U₂|

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

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

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

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

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

أمثلة حسابية

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

 تم تقريب السّمُوت إلى 3 منازل عشرية فقط من الثانية، حيث أن:

هذه الدقة كافية لإعادة إنتاج المواضع بدقة تصل إلى 0.5 مم عند استخدام الحل المباشر،

هذا ينطبق بشكل خاص عندما تكون σ تساوي π أو 2π أو قريبة منها

الجدول 1: حلول النقاط المتقابلة (Antipodal Solutions)
المتغيرالمثال 1المثال 2المثال 3المثال 4φ₁41∘41′45.880∘00′00.00030∘00′00.0060∘00′00.00φ₂−41∘41′46.20∘00′00.000−30∘00′00.0−59∘59′00.0L179∘59′59. 44179∘41′49.78063179∘40′00.0′179∘50′00.0α₁179∘58′49.16330∘00′00.0039∘24′51.80629∘11′51.07α₂0∘01′10.838150∘00′00.00140∘35′08.194′′150∘49′06.868′′s  متر20,004,566.722819,996,147.416919,994,364.606920,000,433.9629
 تم حسابه باستخدام: لغة برمجة  FORTRAN ، حاسوب IBM 7040/94 

 الجدول 2: حلول الجانب الخلفي (Back Side Solutions)
المتغيرالمثال 1المثال 2المثال 3المثال 4φ₁41∘41′45.880′′0∘00′00.000′′30∘00′00.000′′60∘00′00.000′′φ₂41∘41′46.200′′0∘00′00.000′′30∘00′00.000′′59∘59′00.000′′L0∘00′00.560′′0∘18′10.219′′0∘20′00.000′′0∘10′00.000′′α₁180∘00′35.423′′194∘28′47.448′′198∘30′47.488′′344∘56′31.727′′α₂180∘00′35.423′′194∘28′47.448′′198∘30′47.488′′344∘56′59.622′′s  متر40,009,143.320840,004,938.272240,004,046.711440,006,087.0024
تم حسابه بواسطة برنامج الآلة الحاسبة المكتبية Wang 720

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

 مصادر 

 مراجع 

  English translation of Astron. Nachr. 4, 241–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 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

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