Определите положение ECEF

Поскольку система координат ECEF является геоцентрической, можно использовать сферические уравнения, чтобы определить положение по широте, долготе и геоцентрическому радиусу.

Вычислите геоцентрические радиусы в метрах в массиве широт от +/-90 степеней с помощью geocradius.

lat = -90:90;
r = geocradius( lat );
rlat = convang( lat, 'deg', 'rad');
z = r.*sin(rlat);

Поскольку долгота не влияет на зональную гармоническую гравитационную модель, примите, что положение y равно нулю.

x = r.*cos(rlat);
y = zeros(size(x));

Вычисление зональной гармонической гравитации для Земли

Использование gravityzonal вычислить массив зональной гармонической силы тяжести в координатах ECEF для массива позиций ECEF в метрах в секунду в квадрате.

[gx_zonal, gy_zonal, gz_zonal] = gravityzonal([x' y' z']);

Вычислите WGS84 силу тяжести

Использование gravitywgs84 вычислить WGS84 силу тяжести в нисходящей и северной осях на поверхности Земли, массив геодезических широт в степенях и 0 степенях долготы с помощью точного метода с атмосферой, без центробежных эффектов и без предварительной обработки.

lat_gd = geoc2geod(lat, r);
long_gd = zeros(size(lat));
[gd_wgs84, gn_wgs84] = gravitywgs84( zeros(size(lat)), lat_gd, long_gd, 'Exact', [false true false 0] );

Определите силу тяжести для сферической Земли

Вычислите массив сферической гравитации для массива геоцентрических радиусов в метрах в секунду в квадрате с помощью гравитационного параметра Земли в метрах в кубе в секунду в квадрате.

GM  = 3.986004415e14;
gd_sphere = -GM./(r.*r);

Графики сравнения для различных гравитационных моделей

Чтобы сравнить гравитационные модели, их выходы должны быть в одной системе координат. Можно преобразовать зональную силу тяжести из координат ECEF в координаты NED с помощью матрицы Direction Cosine из dcmecef2ned.

dcm_ef = dcmecef2ned( lat_gd, long_gd );
gxyz_zonal = reshape([gx_zonal gy_zonal gz_zonal]', [3 1 181]);
for m = length(lat):-1:1
    gned_zonal(m,:) = (dcm_ef(:,:,m)*gxyz_zonal(:,:,m))';
end

figure(1)
plot(lat_gd, gned_zonal(:,3), lat, -gd_sphere, '--', lat_gd, gd_wgs84, '-.')
legend('Zonal Harmonic', 'Spherical', 'WGS84','Location','North')
xlabel('Geodetic Latitude (degrees)')
ylabel('Down gravity (m/s^2)')
grid

Фигура 1: Гравитация по оси Вниз в метрах в секунду за секунду

figure(2)
plot( lat_gd, gned_zonal(:,1), lat_gd, gn_wgs84, 'r-.' )
hold on
plot([lat(1) lat(end)], [0 0], 'g--' )
legend('Zonal Harmonic', 'WGS84', 'Spherical','Location','Best')
xlabel('Geodetic Latitude (degrees)')
ylabel('North gravity (m/s^2)')
grid
hold off

Фигура 2: Сила тяжести на северной оси в метрах в секунду за секунду

Вычислите общую силу тяжести для WGS84 и из вектора зональной силы тяжести в метрах в секунду за секунду.

gtotal_wgs84 = gravitywgs84( zeros(size(lat)), lat_gd, zeros(size(lat)), 'Exact', [false true false 0] );
gtotal_zonal = sqrt(sum([gx_zonal.^2 gy_zonal.^2 gz_zonal.^2],2));

figure(3)
plot( lat, gtotal_zonal, lat_gd, gtotal_wgs84, '-.' )
legend('Zonal Harmonic', 'WGS84','Location','North')
xlabel('Geodetic Latitude (degrees)')
ylabel('Total gravity (m/s^2)')
grid

Фигура 3: Общая гравитация в метрах в секунду за секунду

Сравнение гравитационных моделей с центробежными эффектами

Теперь вы видели сравнение гравитации невращающейся Земли. Исследуйте центробежные эффекты от вращения Земли на гравитационных моделях.

Вычисление гравитационных центробежных эффектов для Земли

Использование gravitycentrifugal вычислить массив центробежных эффектов в координатах ECEF для массива позиций ECEF в метрах в секунду в квадрате.

[gx_cent, gy_cent, gz_cent] = gravitycentrifugal([x' y' z']);

Добавьте центробежные эффекты к зональной гармонической гравитации.

gx_cent_zonal = gx_zonal + gx_cent;
gy_cent_zonal = gy_zonal + gy_cent;
gz_cent_zonal = gz_zonal + gz_cent;

Вычислите WGS84 силу тяжести с центробежными эффектами

Использование gravitywgs84 вычислить WGS84 силу тяжести в нисходящей и северной осях на поверхности Земли, массив геодезических широт в степенях и 0 степенях долготы с помощью точного метода с атмосферой, центробежными эффектами и без предварительной обработки.

[gd_cent_wgs84, gn_cent_wgs84] = gravitywgs84( zeros(size(lat)), lat_gd, long_gd, 'Exact', [false false false 0] );

Вычислите общую силу тяжести с центробежными эффектами для WGS84 и из вектора зональной силы тяжести в метрах в секунду за секунду.

gtotal_cent_wgs84 = gravitywgs84( zeros(size(lat)), lat_gd, zeros(size(lat)), 'Exact', [false false false 0] );
gtotal_cent_zonal = sqrt(sum([gx_cent_zonal.^2 gy_cent_zonal.^2 gz_cent_zonal.^2],2));

Графики сравнения для различных гравитационных моделей с центробежными эффектами

Чтобы сравнить гравитационные модели, их выходы должны быть в одной системе координат. Можно преобразовать зональную силу тяжести из координат ECEF в координаты NED с помощью матрицы Direction Cosine из dcmecef2ned. На рисунке 5 можно увидеть некоторое различие между зональной гармонической гравитацией с центробежными эффектами и WGS84 гравитацией с центробежными эффектами. Большинство различий происходит из-за различий между зональной гармонической гравитацией и WGS84 гравитацией.

gxyz_cent_zonal = reshape([gx_cent_zonal gy_cent_zonal gz_cent_zonal]', [3 1 181]);
for m = length(lat):-1:1
    gned_cent_zonal(m,:) = (dcm_ef(:,:,m)*gxyz_cent_zonal(:,:,m))';
end

figure(4)
plot(lat_gd, gned_cent_zonal(:,3), lat_gd, gd_cent_wgs84, '-.')
legend('Zonal Harmonic', 'WGS84','Location','North')
xlabel('Geodetic Latitude (degrees)')
ylabel('Down gravity (m/s^2)')
grid

Фигура 4: Гравитация с центробежными эффектами на оси Вниз в метрах в секунду за секунду

figure(5)
plot( lat_gd, gned_cent_zonal(:,1), lat_gd, gn_cent_wgs84, '--', lat_gd, (gned_zonal(:,1)-gn_wgs84'), '-.' )
axis([-100 100 -0.0002 0.0002])
legend('Zonal Harmonic', 'WGS84', 'Error Between Models w/o Centrifugal Effects', 'Location','Best')
xlabel('Geodetic Latitude (degrees)')
ylabel('North gravity (m/s^2)')
grid

Фигура 5: Сила тяжести на северной оси в метрах в секунду за секунду

figure(6)
plot( lat, gtotal_cent_zonal, lat_gd, gtotal_cent_wgs84, '-.' )
legend('Zonal Harmonic', 'WGS84','Location','North')
xlabel('Geodetic Latitude (degrees)')
ylabel('Total gravity (m/s^2)')
grid

Фигура 6: Общая гравитация с центробежными эффектами в метрах в секунду за секунду