Программирование перевода координат из одной системы в другую.

stout, chnav, Спасибо за ответы!

Пробовал, написал свои функции вычисления Sin() и Cos() разложением в ряд Тейлора по этим формулам:

В цикле вычислял до десятого члена ряда.. Результат остался тем же. Потому и грешил на тип данных.

--- Сообщения объединены, , Оригинальное время сообщения: ---

Уважаемый stout, а не подскажете, где можно ознакомиться с теорией строгих формул, тех, которые вместо Бурша-Вольфа?
Спасибо.

 

Это здря, ибо смысла в 10 членах нет. Другое дело, если бы ограничились парой членов.

У самого Бурши. Только обратите внимание на знаки углов.

 

stout, Спасибо большое! Буду изучать
По поводу количества членов ряда. Пробовал с различными значениями. Первый вариант функций вообще предусматривал выход из цикла, если величина члена ряда опустится ниже заданной точности. Точность выставил 1.0Е-19, и выход из цикла происходил после первой же итерации, что еще больше укрепило мои подозрения, что результат вычислений выходит за границы мантиссы, т.е. не помещается в область значащих позиций.

ЗЫ: Спасибо за функцию GetPrecisionMode. Раньше как то надобности небыло, а сейчас потянул за "веревочку", а там.. "Там столько интересного! Во!" (© Мультик "Крылья, ноги и хвосты" )

 

В Builder, как одном из вариантов С++, есть макросы LDBL_EPSILON = 1.084202172485504434E-19 и DBL_EPSILON =2.2204460492503131E-16
То есть LDBL_EPSILON = 2^-63 , а DBL_EPSILON = 2^-52
Сам машинный эпсилон определяется как минимальное положительное число, удовлетворяющие неравенству 1 + EPSILON > 1.

Сравнение с 1.0Е-19 чревато бесконечным циклом. Обычно сравнивают с n×EPSILON, где n – от единицы до восьми.
Для Delphi аналогов не припомню.
Для установки точности в gcc использую макросы

Код:

#define SET_PREC_DOUBLE \
  unsigned short __old_cw, __new_cw; \
  asm volatile ("fnstcw %0":"=m" (__old_cw)); \
  __new_cw = (__old_cw & ~0x300) | 0x200; \
  asm volatile ("fldcw %0": :"m" (__new_cw));
 
#define SET_PREC_LONG_DOUBLE \
  unsigned short __old_cw, __new_cw; \
  asm volatile ("fnstcw %0":"=m" (__old_cw)); \
  __new_cw = (__old_cw & ~0x200) | 0x300; \
  asm volatile ("fldcw %0": :"m" (__new_cw));
 
#define RESTORE_PREC  asm volatile ("fldcw %0": :"m" (__old_cw));

 

Спасибо большое! Сижу разбираюсь с числами с плавющей точкой, в первый раз потребовалось так глубоко залезть.
Думаю, Ваша информация также будет очень полезна.

ЗЫ: Просмотрел алгоритм который Вы порекомендовали, именно его я и реализовал в своей программе, только взял его не из первоисточника. Спасибо!

 

Я много о чем болтаю что рекомендую. О каком конкретно алгоритме идёт речь?

 

Вот этот: "Извлечен. cтраницы из... Бурша М. - Основы космической геодезии I- 1971.pdf"

 

Больше 20 назад, ещё на Turbo C и asm реализовал арифметику с плавающей точкой, расширив представление long double до 16 байт, т.е. где-то 32-33 десятичных знака. Алгоритмы брал из Библии Программиста (Кнут).
Число представлялось структурой

Код:

 typedef struct  {
 unsigned char	 tail[6];
 long	 double   head   ;
 } FloatDec33;

Вот тогда пришлось повозиться с внутренним представление чисел, порядком байт в слове и бит в байте.

Оффтоп

Правда, уже был опыт работы на больших машинах и на ассемблере в том числе.
Когда делали обращение матрицы методом окаймления при выводе ПЗ-85 на ЕС-1060 (Машина требовала помещение площадью 200 м² и потребляла 80 кВА.),
Саша Санаров нашёл ошибку в машинном микрокоде умножения длинных чисел (PL/I Float Dec (33)) и,
чтобы получить правильный результат, надо было переставить байты в словах,
а после умножения всё вернуть "взад".
И эта (censored) была только на этой машине, но лучшего тогда вообще не было.
Страшно сказать, оперативной памяти был целый реальный мегабайт, то есть просто дофига!

А откуда, если не секрет?

 

Не секрет. Из Ваших советов пользователям этого форума.
Вы упоминали в своих постах алгоритм 7 параметрического преобразования с не упрощенной матрицей, а с полными матрицами разворотов вокруг каждой из осей, с честными синусами и косинусами. Вот я и попробовал его реализовать. Некоторые части для реализации нашел, снова таки по Вашей ссылке в одной из тем форума, здесь: http://vadimchazov.narod.ru, остальное достроил сам. Но прочитав главу из Бурши, который Вы любезно прислали, понял, что у меня получилось реализовать именно этот алгоритм.
Алгоритм работает, но в процессе тестирования сходимости алгоритмов прямого и обратного преобразования, была найдена потеря точности в работе этих цепочек. Она проявляется только на параметрах с разворотами вокруг осей, и составляет расхождение в 9 знаке после запятой в градусах, или 0.2мм при пересчете точки саму в себя в СК-42. Понимаю, мелочь, возможно можно пренебречь, но пренебрегать не хочется.
Как то так.

 

Такого не должно быть. По определению. Могу посмотреть не замыленным взглядом код, или приведите данные для контрольного примера, т.е. координаты самой точки и 7 параметров. Набросаю тестовый пример.

 

Кстати, я тоже сейчас обдумываю аналогичный вариант - заполнение числа Extended побайтно, но не слишком уверен что это поможет решить задачу.

Оффтоп

Во время прохождения службы, был оператором электронно - вычислительной техники стрельбовой РЛС. ЭВМ которую обслуживал имела 2кБ ОЗУ и 4кб ПЗУ. Я своими глазами видел как ПРОШИВАЮТ магнитное ПЗУ: иголка с продетой проволочкой, ее проводят вдоль рядов магнитных сердечников, где нужна 1 - делают виток вокруг сердечника, где 0 - проводят мимо. Сама процедура очень напоминала вышивание крестиком

--- Сообщения объединены, , Оригинальное время сообщения: ---

Был бы очень благодарен за свежий взгляд на код с Вашей стороны. Может лучше попробовать через личные сообщения?

 

В будущем приложении планирую получать координаты gps, однако работать с градусами невозможно. Один градус по широте не равен одному градусу по долготе в переводе на метры да еще и соотношение м/градусы по широте уменьшается с отдалением от экватора.
У меня есть следующая задача: я хочу взять выбранную точку на карте, преобразовать ее в систему координат, которая гарантирует относительное соблюдение пропорций в любой точке на поверхности Земли, добавить +1 км по x и +1 км по y, преобразовать результат обратно углы и отобразить на карте. При этом чтобы Google Maps показывало расстояние между выбранной и найденой точками равное sqrt(2) км не зависимо от положения начальной точки. Желательно иметь точность +-100м.
Вопрос в следующем - существует ли такая система координат, как называется и где можно почитать.
И еще - наткнулся на UTM. Гарантирует ли она соблюдение пропорций? Если да - где можно почитать как над ней выполнять операции, потому что эти зоны просто сбивают с толку.

 

Запутался, помогите найти ошибку.
Пытаюсь по формулам ГОСТа (22), (23) получить значения аналогичные рассчитанным на сайте mapbasic/
Вот код.


SET DECIMALS TO 14
Pi_= 3.14159265358979
ro= 206264.8062
*
aP = 6378245
alP = 1 / 298.3
e2P = 2 * alP - alP*alP
*
aW = 6378137
alW = 1 / 298.257223563
e2W = 2 * alW - alW*alW
*
a = (aP + aW) / 2
e2 = (e2P + e2W) / 2
da = aW - aP
de2 = e2W - e2P
*
dx = 23.57
dy = -140.95
dz = -79.8
wx = 0
wy = -0.35
wz = -0.79
ms = -0.00000022
Bd=56
LD=37
******* wgs -> 42
******* должны получить 56.0,37,0 --> 55.99996418016184, 37.001898209869594
B = Bd * Pi_ / 180
L = Ld * Pi_ / 180
H=0
M = a * (1 - e2) / (1 - e2 * Sin(B) ** 2) ** 1.5
N = a * (1 - e2 * Sin(B) ** 2) ** -0.5

*****ïåðâûå çíà÷åíèÿ
B1=Bd-fdB(B,L,H)/3600
***55.99996419792149
? B1
***37.00189822705185
L1=Ld-fdL(B,L,H)/3600
? L1
H1=H+fdH(B,L,H)
*** 5.364847404449037
? H1
***** âòîðûå çíà÷åíèÿ
Bdi=B1
Ldi=L1
Hi=H1
B = Bdi * Pi_ / 180
L = Ldi * Pi_ / 180
B2=Bd-fdB(B,L,Hi)/3600
L2=Ld-fdL(B,L,Hi)/3600
H2=Hi+fdH(B,L,H)
Bend=(B1+B2)/2
MESSAGEBOX(Bend)
Lend=(L1+L2)/2
Hend=(H1+H2)/2
Function fdB
LPARAMETERS B, L, H
dB = ro / (M + H) * (N / a * e2 * Sin(B) * Cos(B) * da + (N ** 2 / a ** 2 + 1) * N * Sin(B) * Cos(B) * de2 / 2 - (dx * Cos(L) + dy * Sin(L)) * Sin(B) + dz * Cos(B)) - wx * Sin(L) * (1 + e2 * Cos(2 * B)) + wy * Cos(L) * (1 + e2 * Cos(2 * B)) - ro * ms * e2 * Sin(B) * Cos(B)
RETURN dB
ENDFUNC
Function fdL
LPARAMETERS B, L, H
dL = ro / ((N + H) * Cos(B)) * (-dx * Sin(L) + dy * Cos(L)) + Tan(B) * (1 - e2) * (wx * Cos(L) + wy * Sin(L)) - wz
RETURN dL
ENDFUNC
Function fdH
LPARAMETERS B, L, H
dH = -a / N * da + N * Sin(B) ** 2 * de2 / 2 + (dx * Cos(L) + dy * Sin(L)) * Cos(B) + dz * Sin(B) - N * e2 * Sin(B) * Cos(B) * (wx / ro * Sin(L) - wy / ro * Cos(L)) + (a ** 2 / N + H) * ms
RETURN dH
ENDFUNC


Что неправильно? Помогите пожалуйста умные головы.

 

Знаки 7-параметров.

Код:

#include <cstdlib>
#include <fstream>
#include <iostream>
#include <stdint.h>
#include <stdio.h>
#include <cmath>
 
using namespace std;
 
typedef double REAL;
 
const REAL DegPerRad = 45/atan(REAL(1));
const REAL ArcSecPerRad = 3600*DegPerRad;
 
//
// ===========================================================================
//
inline REAL Int(const REAL& Value)
{
    return (Value > 0) ? floor(Value) : ceil(Value);
}
//
// ===========================================================================
//
inline REAL Sqr(const REAL& Value)
{
    return Value*Value;
}
//
// ===========================================================================
//
inline REAL DmsToDeg(const REAL& DMS)
{
    REAL D = Int(DMS/10000);
    return ((DMS/200 - 32*D - Int((DMS - D*10000)/100)/5)/18);
}
//
// ===========================================================================
//
inline REAL DmsToRad(const REAL& DMS)
{
    return DmsToDeg(DMS)/DegPerRad;
}
//
// ===========================================================================
//
inline REAL DegToDms(const REAL& Deg)
{
    REAL D = Int(Deg);
    return (9*Deg + 16*D + Int(60*(Deg - D))/10)*400;
}
//
// ===========================================================================
//
inline REAL RadToDms(const REAL& Rad)
{
    return DegToDms(Rad*DegPerRad);
}
//
// ===========================================================================
//
void BLH2XYZ(const REAL& a // SemimajorAxis
            ,const REAL& f // Flattening
            ,const REAL& B
            ,const REAL& L
            ,const REAL& H
            ,REAL& X
            ,REAL& Y
            ,REAL& Z)
{
 
    REAL ee = f*(2 - f);
    REAL sinB = sin(B);
    REAL cosB = cos(B);
    REAL N = a/sqrt(1 - ee*Sqr(sinB));
    X = (N + H)*cosB*cos(L);
    Y = (N + H)*cosB*sin(L);
    Z = (N*(1 - ee) + H)*sinB;
 
} // end void BLH2XYZ
//
// ===========================================================================
//
REAL Bowring(const REAL& a // SemimajorAxis
            ,const REAL& f // Flattening
            ,const REAL& X // Cartesian coordinates (X , Y, Z)
            ,const REAL& Y
            ,const REAL& Z)
{
    REAL ee = f*(2 - f);
    REAL eea = ee*a;
    REAL RRxy = X*X + Y*Y;
    REAL Rxy = sqrt(RRxy);
    REAL s = sqrt(Z*Z/(1 - ee) + RRxy);
    REAL TanB = Z/Rxy*s/(s - eea);
    REAL TanBeta = (1 - f)*TanB;
    // при таком начальном значении обеспечивается точность, как минимум, 16 знаков мантиссы
    // для tan(B) в диапазоне высот -1 000 км. < H < 20 000 км.
 
    REAL SqrCosBeta = 1/(1 + Sqr(TanBeta));
    REAL cosBeta = sqrt(SqrCosBeta);
    REAL sinBeta = TanBeta*cosBeta;
    REAL SqrSinBeta = 1 - SqrCosBeta;
    TanB = (Z + eea/(1 - f)*SqrSinBeta*sinBeta)/(Rxy - eea*SqrCosBeta*cosBeta);
 
    return TanB;
}
//
// ===========================================================================
//
void Helmert(const REAL XYZA[], const REAL Shift[], const REAL angle[], const REAL& scale, REAL XYZB[])
{
 
    XYZB[0] = Shift[0] + (1 + scale)*(XYZA[0] + angle[2]*XYZA[1] - angle[1]*XYZA[2]);
    XYZB[1] = Shift[1] + (1 + scale)*(XYZA[1] - angle[2]*XYZA[0] + angle[0]*XYZA[2]);
    XYZB[2] = Shift[2] + (1 + scale)*(XYZA[2] + angle[1]*XYZA[0] - angle[0]*XYZA[1]);
}
//
// ===========================================================================
//
void XYZ2BLH(const REAL& a // SemimajorAxis
            ,const REAL& f // Flattening
            ,const REAL& X // Cartesian coordinates (X , Y, Z)
            ,const REAL& Y
            ,const REAL& Z
            ,REAL& B
            ,REAL& L
            ,REAL& H)
{
 
    REAL TanB = Bowring(a, f, X, Y, Z);
    B = atan(TanB);
    REAL SqrCosB = 1/(1 + Sqr(TanB));
    REAL cosB = sqrt(SqrCosB);
    REAL sinB = TanB*cosB;
    REAL Rxy = sqrt(Sqr(X) + Sqr(Y));
    L = 2*atan(Y/(X + Rxy));                                       // формула Ваничека
    H = Rxy*cosB + Z*sinB - a*sqrt(1 - f*(2 - f)*Sqr(sinB)); // не самый лучший подход
}
 
//
// ===========================================================================
//
void Molodensky(const REAL& SemiMajorAxisA
                ,const REAL& SemiMajorAxisB
                ,const REAL& RecipFlatteningA
                ,const REAL& RecipFlatteningB
                ,const REAL Shift[]
                ,const REAL angle[] // radian
                ,const REAL& scale
                ,const REAL& BA // radian
                ,const REAL& LA // radian
                ,const REAL& HA
                ,REAL& BB // radian
                ,REAL& LB // radian
                ,REAL& HB)
{
    REAL deltaB = 0;
    REAL deltaL = 0;
    REAL deltaH = 0;
    REAL dX = Shift[0];
    REAL dY = Shift[1];
    REAL dZ = Shift[2];
    REAL Wx = angle[0];
    REAL Wy = angle[1];
    REAL Wz = angle[2];
 
    cout << "============= proc Molodensky ===============" << '\n';
 
    REAL fA = 1/RecipFlatteningA; // flattening System A
    REAL fB = 1/RecipFlatteningB; // flattening System B
    REAL sqr_excentricityA = fA*(2 - fA);
    REAL sqr_excentricityB = fB*(2 - fB);
    REAL dA = SemiMajorAxisB - SemiMajorAxisA;
    REAL a = 0.5*(SemiMajorAxisB + SemiMajorAxisA);
    REAL sqr_e = 0.5*(sqr_excentricityA + sqr_excentricityB);
    REAL temp = 1/(RecipFlatteningA*RecipFlatteningB);
    REAL d_sqr_e = (RecipFlatteningA - RecipFlatteningB)*temp*(2 - temp*(RecipFlatteningA + RecipFlatteningB));
    //    d_sqr_e = sqr_excentricityB - sqr_excentricityA;
    int iter_count = 3;
    for(int k = iter_count; k > 0; k--) {
        REAL B = BA + 0.5*deltaB;
        REAL L = LA + 0.5*deltaL;
        REAL H = HA + 0.5*deltaH;
        REAL sinB = sin(B);
        REAL cosB = cos(B);
        REAL cos2B = (cosB - sinB)*(cosB + sinB);
        REAL sinL = sin(L);
        REAL cosL = cos(L);
        REAL W = 1/sqrt(1 - sqr_e*Sqr(sinB));
        REAL N = a*W;
        REAL M = a*(1 - sqr_e)*W*W*W;
        deltaB = (N/a*sqr_e*sinB*cosB*dA + 0.5*(Sqr(N/a) + 1)*N*sinB*cosB*d_sqr_e -
                 (dX*cosL + dY*sinL)*sinB + dZ*cosB)/(M + H) +
                 (Wy*cosL - Wx*sinL)*(1 + sqr_e*cos2B) - scale*sqr_e*sinB*cosB;
        deltaL = (dY*cosL - dX*sinL)/(cosB*(N + H)) + sinB/cosB*(1 - sqr_e)*(Wx*cosL + Wy*sinL) - Wz;
        deltaH = 0.5*N*Sqr(sinB)*d_sqr_e - a/N*dA + (dX*cosL + dY*sinL)*cosB + dZ*sinB -
            N*sqr_e*sinB*cosB*(Wx*sinL - Wy*cosL) + scale*(Sqr(a)/N + H);
 
        cout << " ======= iteration: " << iter_count + 1 - k << '\n';
        cout << "   deltaB = " << deltaB << '\n';
        cout << "   deltaL = " << deltaL << '\n';
        cout << "   deltaH = " << deltaH << '\n';
    }
    BB = BA + deltaB;
    LB = LA + deltaL;
    HB = HA + deltaH;
}
 
int main(int argc, char** argv)
{
 
    REAL SemiMajorAxisA = 6378137;
    REAL SemiMajorAxisB = 6378245;
    REAL RecipFlatteningA = 298.257223563;
    REAL RecipFlatteningB = 298.3;
    REAL fA, fB;
    REAL Shift[3] = { -23.57, +140.95, +79.8 };
    REAL angle[3] = { 0, +0.35, +0.79 };
    REAL scale = +0.00000022;
    REAL BA = 560000.0000;
    REAL LA = 370000.0000;
    REAL HA;
    REAL BB;
    REAL LB;
    REAL HB;
    REAL X, Y, Z;
    REAL XYZA[3];
    REAL XYZB[3];
    char ans;
 
    cout.precision(12);
    cout << fixed << endl;
 
    cout << "\n manual input parametr (Y/N)?";
    cin >> ans;
 
    if(toupper(ans) == 'Y') {
        cout << "Please enter for System A" << endl << "  Semimajor axis in meter = ";
        cin >> SemiMajorAxisA;
        cout << endl << "reciprocal flattening (1/f) = ";
        cin >> RecipFlatteningA;
 
        cout << "Please enter for System B" << endl << "  Semimajor axis in meter = ";
        cin >> SemiMajorAxisB;
        cout << endl << "reciprocal flattening (1/f) = ";
        cin >> RecipFlatteningB;
 
        cout << endl << "Please enter Helmert param:";
        cout << endl << "Shift X in meter = ";
        cin >> Shift[0];
        cout << endl << "Shift Y in meter = ";
        cin >> Shift[1];
        cout << endl << "Shift Z in meter = ";
        cin >> Shift[2];
        cout << endl << "angle of rotation Wx in arc. sec. = ";
        cin >> angle[0];
        cout << endl << "angle of rotation Wy in arc. sec. = ";
        cin >> angle[1];
        cout << endl << "angle of rotation Wz in arc. sec. = ";
        cin >> angle[2];
        cout << endl << "(1 - scale) = ";
        cin >> scale;
 
        cout << endl << "Please enter Geodetic:";
        cout << endl << "Latitude  in format ddmmss.ss... = ";
        cin >> BA;
        cout << endl << "Longitude in format ddmmss.ss... = ";
        cin >> LA;
        cout << endl << "H in meter = ";
        cin >> HA;
    }
 
    for(int k = 0; k < 3; k++)
        angle[k] /= ArcSecPerRad;
    BA = DmsToRad(BA);
    LA = DmsToRad(LA);
    fA = 1/RecipFlatteningA;
    fB = 1/RecipFlatteningB;
 
    cout << '\n' << "System A:" << '\n';
    cout << "  B(dms) = " << RadToDms(BA) << '\n';
    cout << "  L(dms) = " << RadToDms(LA) << '\n';
    cout << "  B(degree) = " << DegPerRad*BA << '\n';
    cout << "  L(degree) = " << DegPerRad*LA << '\n';
    cout << "  H      = " << HB << '\n';
 
 
    Molodensky(SemiMajorAxisA, SemiMajorAxisB, RecipFlatteningA, RecipFlatteningB, Shift, angle, scale, BA, LA, HA, BB,
        LB, HB);
 
    cout << '\n' << "System B:" << '\n';
    cout << "  B(dms) = " << RadToDms(BB) << '\n';
    cout << "  L(dms) = " << RadToDms(LB) << '\n';
    cout << "  B(degree) = " << DegPerRad*BB << '\n';
    cout << "  L(degree) = " << DegPerRad*LB << '\n';
    cout << "  H      = " << HB << '\n';
 
    BLH2XYZ(SemiMajorAxisA, fA, BA, LA, HA, X, Y, Z);
 
    XYZA[0] = X;
    XYZA[1] = Y;
    XYZA[2] = Z;
    cout << '\n' << "System A:" << '\n';
    cout << " X = " << X << '\n';
    cout << " Y = " << Y << '\n';
    cout << " Z = " << Z << '\n';
 
    Helmert(XYZA, Shift, angle, scale, XYZB);
 
    X = XYZB[0];
    Y = XYZB[1];
    Z = XYZB[2];
    cout << '\n' << "System B:" << '\n';
    cout << " X = " << X << '\n';
    cout << " Y = " << Y << '\n';
    cout << " Z = " << Z << '\n';
 
 
    XYZ2BLH(SemiMajorAxisB, fB, X, Y, Z, BB, LB, HB);
 
    cout << '\n' << "Helmert transform" << '\n';
    cout << "  B(dms) = " << RadToDms(BB) << '\n';
    cout << "  L(dms) = " << RadToDms(LB) << '\n';
    cout << "  B(degree) = " << DegPerRad*BB << '\n';
    cout << "  L(degree) = " << DegPerRad*LB << '\n';
    cout << "  H      = " << HB << '\n';
 
    system("PAUSE");
    return EXIT_SUCCESS;
}

 

Я бы заменил все эти вводы с клавы на:

Код:

...
if (argc > 4) SemiMajorAxisA = atoi(argv[4]);
...

Или вообще бы читал их из файла.