Вища геодезія

% Пiставляємо значення змiнних u та v та спрощуємо отриманi вирази.

% В результатi змiнним r та n присвоєнi вирази, якi залежать лише вiд t

r= simplify(subs(r));

n= simplify(subs(n));

N = cross(cross(dr, d2r), dr);

% Кривизна кривої як функцiя t

k = norm(cross(dr, d2r)) / norm(dr)^3;

% Вектор кривизни як функцiя t

k_vector = k * N / norm(N);

% Нормальна i геодезична кривизна як функцiї t

kn = k * abs(dot(N, n)) / norm(N) / norm(n); kg = sqrt(k^2 - kn^2);

% Значення t, яке буде пiдставлено в аналiтичнi вирази t = 0.5;

% Видаємо результати у вiкно Web Browser

ShowSym(...

3.9. Геодезична лiнiя

Геодезичною лiнiєю називається така крива на поверхнi, у якої геодезична кривизна дорiвнює нулю в усiх її точках. Серед геометричних властивостей геодезичної лiнiї найважливiшою є те, що вона є найкоротшим шляхом мiж двома точками на криволiнiйнiй поверхнi.

Для рiвностi нулю геодезичної кривизни лiнiї на поверхнi необхiдно i

#„

достатньо, щоб вектор кривизни k кривої був в усiх точках направлений по

#„

нормалi до поверхнi, а тому достатньо, щоб вектор N був перпендикулярний до поверхнi всюди вздовж кривої.

Пiд натуральним параметром s для даної точки кривої розумiють довжину дуги кривої мiж даною точкою i деякою точкою кривої, прийнятою

при перемiщеннi точки вздовж кривої змiнюються не тiльки її криволiнiйнi

координати u, v, а може змiнюватись також i кут θ, тобто для кривої на поверхнi величини u, v, θ є функцiями натурального параметру s.

На поверхнi, координатна сiтка якої є ортогональною, геодезична лiнiя описується наступною системою диференцiальних рiвнянь:

Першi двi початковi умови визначають криволiнiйнi координати (u0, v0) точки поверхнi, через яку проходить геодезична в початковiй точцi (при s = 0), а остання умова задає напрямок геодезичної в початковiй точцi. Єдинiсть розв’язку означає, що через точку поверхнi з координатами (u0, v0) в заданому напрямку проходить лише одна геодезична лiнiя.

Обчислення похiдних та спрощення отриманого виразу можна автоматизувати, якщо скористатись можливостями символьної математики системи MATLAB, як це показано в наступному прикладi.

clc, clear

% Оголошуємо символьнi змiннi

syms a e B L

%Задаємо властивостi символьних змiнних − область їх допустимих значень

%(це необхiдно для коректного спрощення символьних виразiв)

assume(a>0 & e>0 & e<1 & B>-pi/2 & B<pi/2)

% Вирази для залежностей N(B), X(B),Y (B), Z(B)

N = a/sqrt(1-e^2*sin(B)^2); X = N * cos(B) * cos(L);

Y = N * cos(B) * sin(L);

Z = N *(1-e^2) * sin(B);

% Знаходимо вирази для перших та других похiдних та заносимо їх в змiннi dX = diff(X, B);

dY = diff(Y, B); dZ = diff(Z, B);

d2X = diff(X, B, 2); d2Y = diff(Y, B, 2); d2Z = diff(Z, B, 2);

%Промiжнi обчислення: окремо обчислюємо та спрощуємо чисельник та

%знаменник в пiдкореневому виразi формули (4.7).

%Застосування функцiй factor (розкладання на множники) та simplify

p1 = simplify(factor(dX^2 + dY^2 + dZ^2));

p2 = simplify((dX*d2Y - d2X*dY)^2 + (dX*d2Z - d2X*dZ)^2 + ...

(dZ*d2Y - d2Z*dY)^2);

В загальному випадку нормальна площина може перетинатися з меридiаном пiд довiльним кутом A. В цьому випадку радiус кривизни нормального перерiзу буде виражатись формулою

Очевидно, пiдставляючи в (4.10) A= 0 та A= π/2, одержимо вiдповiдно радiуси кривизни меридiану та першого вертикалу.

При будь-якому розташуваннi нормального перерiзу, проведеного через точку Q, вiн буде залишатись плоскою кривою, а тому вектор головної нормалi до нього, побудований в точцi Q буде лежати в нормальнiй площинi, а значить, буде перпендикулярний до поверхнi елiпсоїда. Звiдси випливає, що в точцi Q геодезична кривизна нормального перерiзу дорiвнює 0.

Середнiм радiусом кривизни Rсер поверхнi в точцi називають границю середнього арифметичного значення радiусiв кривизни RA всiх можливих нормальних перерiзiв, проведених через дану точку поверхнi в рiзних напрямках:

Повною або гаусовою кривизною поверхнi називається добуток кривизн головних нормальних перерiзiв поверхнi в данiй точцi

X

Рис. 4.3. До поняття першого вертикалу

MeridianArcLength Функцiя призначена для обчислення довжини дуги меридiана на елiпсоїдi. Формат її виклику наступний:

s = MeridianArcLength(B1, B2, a, e)

Тут B1, B2 — геодезичнi широти початкової та кiнцевої точок дуги (в радiанах); a, e — велика пiввiсь та ексцентриситет елiпсоїда. Довжина дуги повертається в тих же одиницях, в яких задана велика пiввiсь. Обчислення реалiзовано за формулою типу (4.15), в якiй взято дещо бiльшу кiлькiсть доданкiв (врахованi степенi e8, e10, e12).

похибка якої складає не бiльше 0,2 мм при довжинi дуги S < 100 км. В (4.16) Mm означає радiус кривизни меридiана, обчислений в точцi з широ-

похибка якої не перевищує 1,5 мм.

Використовуючи формулу Сiмпсона чисельного iнтегрування, отримає-

B1, B2, Bm вiдповiдно. Похибка формули (4.18) не перевищує 0,2 мм при довжинi дуги до 400 км.

Якщо довжина дуги меридiана перевищує 400–500 км та при використаннi формул (4.16), (4.18) необхiдно зберегти вказану точнiсть обчислення, то дугу слiд роздiлити на частини i до кожної частини окремо застосувати формулу (4.16) або (4.18).

Можна наближено розв’язати i обернену задачу: за вiдомою довжиною дуги меридiана знайти широту (передбачається, що дуга починається на

екваторi). По сутi, ця задача зводиться до розв’язання рiвняння

Z B

де C0, C2, C2, C3, . . . — невiдомi коефiцiєнти. Пiдставимо B1 = 0 та B2 =(4.19)

яка наближено розв’язує поставлену задачу:

де позначено ξ = saм , s = sin ξ, c = cos ξ.

B_via_ArcLength Функцiя за вiдомою довжиною s дуги меридiана, що вiдкладається вiд екватора, на основi формули типу (4.20) повертає вiдповiдне значення геодезичної широти B.

Формат її виклику наступний: