KFF (Все сообщения пользователя)

прямая задача
[CODE]function TrueGeoTask(Bo, Lo, Azim, Dist: double): T3dPoint;
var sinAzim, cosAzim  : double;
   tanU1,cosU1,sinU1 : double;
   sinA,cosSqA,A,B    : double;
   sg1,sig,dsg,psig,uSq  : double;
   cos2sM ,sinSg, cosSg: double;

begin
 sinAzim := sin(Azim);
 cosAzim := cos(Azim);
 tanU1   := (1-eWGS84.rf) * tan(Bo);
 cosU1   := 1 / sqrt(1 + tanU1*tanU1);
 sinU1   := tanU1 * cosU1;
 sg1     := arctan2(tanU1, cosAzim);
 sinA    := cosU1 * sinAzim;
 cosSqA  := 1 - sinA*sinA;
 uSq     := cosSqA*(eWGS84.a*eWGS84.a-eWGS84.b*eWGS84.b)/(eWGS84.b*eWGS84.b);
 A       := 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
 B       := uSq/1024*(256+uSq*(-128+uSq*(74-47*uSq)));
 sig     := Dist/(eWGS84.b * A); psig := 0;
 repeat
   cos2sM := cos(2*sg1 + sig);
   sinSg  := sin(sig);
   cosSg  := cos(sig);
   dsg    := B*sinSg*(cos2sM+0.25*B*(cosSg*(-1+2*cos2sM*cos2sM)-
             B/6*cos2sM*(-3+4*sinSg*sinSg)*(-3+4*cos2sM*cos2sM)));
   psig   := sig;
   sig    := Dist / (eWGS84.b*A) + dsg;
 until (abs(sig-psig) <= 1E-12);
 A   := eWGS84.rf/16*cosSqA*(4+eWGS84.rf*(4-3*cosSqA));
 B   := sinU1*sinSg - cosU1*cosSg*cosAzim;
 result.B := arctan2(sinU1*cosSg + cosU1*sinSg*cosAzim, (1-eWGS84.rf)*sqrt(sinA*sinA + B*B));
 result.L  := Lo+ arctan2(sinSg*sinAzim, cosU1*cosSg - sinU1*sinSg*cosAzim) -
              (1-A) * eWGS84.rf * sinA * (sig + A*sinSg*(cos2sM+A*cosSg*(-1+2*cos2sM*cos2sM)));
 result.Z := arctan2(sinA, -B);
end;
[/CODE]

Иван, предлагаю второй вариант решения обратной геодезической задачи.
Код полностью переписан
[CODE]function ReversGeoTask(const P1,P2 : T3dPoint): T3dPoint;
var L,lam, dlam : double;
   uSq, A, B, C, sig, dsig : double;
   sinsig, cossig, sinSqsig : double;
   cos2sigM  : double;
   iterationLimit : integer;
   cosU1, cosU2 : double;
   sinU1, sinU2 : double;
begin
 L :=  P2.Y - P1.Y;
 A := (1-eWGS84.rf) * tan(P1.B);
 B := (1-eWGS84.rf) * tan(P2.B);
 cosU1 := 1 / sqrt((1 + A*A));
 cosU2 := 1 / sqrt((1 + B*B));
 sinU1 := A * cosU1;
 sinU2 := B * cosU2;
 lam := L;  dlam := 0;
 iterationLimit  := 250;
 repeat
   sinSqsig := (cosU2*sin(lam)) * (cosU2*sin(lam)) + (cosU1*sinU2-sinU1*cosU2*cos(lam)) * (cosU1*sinU2-sinU1*cosU2*cos(lam));
   if (sinSqsig = 0) then exit;
   sinsig   := sqrt(sinSqsig);
   cossig   := sinU1*sinU2 + cosU1*cosU2*cos(lam);
   sig      := arctan2(sinsig, cossig);
   B     := cosU1 * cosU2 * sin(lam) / sinsig;
   A        := 1 - B*B;
   if (A <> 0) then
     cos2sigM := cossig - 2*sinU1*sinU2/A
   else
     cos2sigM := 0;
   C    := eWGS84.rf/16*A*(4+eWGS84.rf*(4-3*A));
   dlam := lam;
   lam  := L + (1-C) * eWGS84.rf * B * (sig + C*sinsig*(cos2sigM+C*cossig*(-1+2*cos2sigM*cos2sigM)));
   dec(iterationLimit);
 until (abs(lam-dlam)<=1E-12) or (iterationLimit<1);

 if (iterationLimit=0) then  exit;

 uSq  := (1 - B*B) * (eWGS84.a*eWGS84.a - eWGS84.b * eWGS84.b) / (eWGS84.b*eWGS84.b);
 A    := 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
 B    := uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
 dsig := B*sinsig*(cos2sigM+B/4*(cossig*(-1+2*cos2sigM*cos2sigM)- B/6*cos2sigM*(-3+4*sinsig*sinsig)*(-3+4*cos2sigM*cos2sigM)));

 result.X := eWGS84.b*A*(sig-dsig);
 result.Y := arctan2(cosU2*sin(lam),  cosU1*sinU2-sinU1*cosU2*cos(lam));
 result.Z := arctan2(cosU1*sin(lam), -sinU1*cosU2+cosU1*sinU2*cos(lam));
 if result.Y<0 then result.Y := result.Y+2*pi;
 if result.Z<0 then result.Z := result.Z+2*pi;
end;

[/CODE]

[QUOTE]P.S. Извините, что беспокою - делаю программу для своего хобби, радиоприема.[/QUOTE]
Чесно говоря работу функции на расстояниях больших чем 1000 км я не проверял, поэтому могут быть какие то неточности с дальностью

Хотя я могу пересмотреть работу алгоритма для разных полушарий и больших расстояний и выложить сюда новый код..
Наверное на выходных так и сделаю

Здравствуйте !

Согласен, в коде была ошибка которую я  исправил очень давно но тут на форуме не обновлял

В сообщении выше я внёс правку. Заклбчаетс яона в следующем
строчку в функции [B]ReversGeoTask[/B][CODE]aLen:=D.X/(P0.b*cos(dA12));[/CODE]нужно исправить вот так [CODE]aLen:=Abs(D.X/(P0.b*cos(dA12))); // <<< + Abs[/CODE]

[QUOTE]1. В структуре T3DPoint используются еще и некие поля L и B. (Пример - последняя строчка функции TrueGeoTask) - я их просто объявил дополнительно к полям X,Y,Z. Они имеют какой-то особый смысл?[/QUOTE]
Структуру лучше объявлять вот так[CODE]typedef struct T3dPoint
{
   union { double X; double x; double lat; double B;};
   union { double Y; double y; double lng; double L;};
   union { double Z; double z; double H;   double elev;};
   T3dPoint() {memset(this, 0, sizeof(T3dPoint));};
   int  operator  == (T3dPoint& P) {return (x == P.x && y == P.y && z == P.z);}
   int  operator  != (T3dPoint& P) {return (x != P.x || y != P.y || z != P.z);}
}
 T3dPoint;[/CODE][QUOTE]2. В Вашем коде так же используется вызов функции adjlon[/QUOTE]
Функция нормализации угла в радианах, ну чтобы у вас не получалось 745 градусов. Вариант на Паскале выглядит так:
[CODE] //* reduce argument to range +/- PI *//
function adjlon(lon:double):;
begin
result:=lon;
while Abs(result)>pi do
case ord(result>0) of
  1: result:=-(2*pi-result);
  0: result:=result+2*pi;
end;
end;

function adjlon(Value:double; Latitude : boolean):double;
var limit : double;
begin
result:=Value;
case Latitude of
  true : limit:=pi;
  false: limit:=2*pi;
end;
while Abs(result)>limit/2 do
case ord(result>0) of
 1: if Latitude then
      result:=limit-result
    else
      result:=-(limit-result);
 0: result:=result+limit;
end;
end; [/CODE]

[QUOTE]Задано СК 42 (Прямоугольная в мерах) перевести в Геодезическую СК (международный эллипсоид WGS 1984).
Напишите метод как это делается  конкретный метод, не нужны ссылки...или пример кода, как в обычных  форумах пишут[/QUOTE]
функция номер 1, преобразование Ваши координат в радианы (без создания карты)
[CODE] // Преобразование координат в метрах на местности из заданной зоны
// в геодезические координаты в системе 42г.
// zone - номер исходной зоны системы 1942г
// x;y  - преобразуемые координаты
// на входе метры в одной зоне 42г.;на выходе - радианы.
// При ошибке возвращает 0

function  mapPlaneToGeo42ByZone(zone : integer; var x,y : double) : integer;
{$IFNDEF LINUXAPI} stdcall {$ELSE} cdecl {$ENDIF} external sGisAcces;[/CODE]результаты функции номер один отправляете в функцию номер 2

[CODE] // Преобразование геодезических координаты в радианах из системы 1942г
// (эллипсоид Красовского) в геодезические координаты в радианах
// (общеземной эллипсоид WGS84) (поддерживается не для всех карт !)
// hmap - идентификатор открытых данных
// Bx;Ly  - преобразуемые координаты
// на входе радианы в 42г.; на выходе - радианы в WGS84
// H     - высота в точке (метры)
// Параметр hmap может быть равен нулю

procedure mapGeo42ToGeoWGS84Ex(var Bx, Ly,H : double);
{$IFNDEF LINUXAPI} stdcall {$ELSE} cdecl {$ENDIF}
external sGisAcces;
function  mapGeo42ToGeoWGS84(Map : HMap;var Bx, Ly : double) : integer;
{$IFNDEF LINUXAPI} stdcall {$ELSE} cdecl {$ENDIF}
external sGisAcces;
function  mapGeo42ToGeoWGS843D(Map : HMap; var Bx, Ly, H : double) : integer;
{$IFNDEF LINUXAPI} stdcall {$ELSE} cdecl {$ENDIF}
external sGisAcces;[/CODE]
например у вас есть
X = +06 283 314 М, Y = 09 307 605 М, N_зоны = 09-08, H = +0122 М. [CODE]
var x, y, h : double;
begin
 x := 6283314;
 y := 307605;  
 h := 122;
 mapPlaneToGeo42ByZone (9, x,y);
 mapGeo42ToGeoWGS843D(0, x,y,h); // в комментарии написано, что "Параметр hmap может быть равен нулю"
 
 ИТОГО в x-широта в y-долгота а WGS84 а в h - пересчитанная высота

end;[/CODE]как то так =)

форум Линукс отвечаю на Си

T3dPoint ReversGeoTask(TEllpsType ellipsoid, T3dPoint p1, T3dPoint p2);

где
struct T3dPoint
{
  double X,
  double Y,
  double Z
} T3dPoint ;

const double M_PI = 3.1415926533;

Вызов:
T3dPoint start;
start.X = 54.125454 * M_PI / 180.0;
start.Y = 46.125454 * M_PI / 180.0;

T3dPoint finish;
finish.X = 55.87452 * M_PI / 180.0;
finish.Y = 47.542112 * 1M_PI / 180.0;

T3dPoint ret = ReversGeoTask(ellWGS84,start,  finish)

ret.X - содержит дистанцию в метрах
ret.Y - содержит прямой азимут в радианах, чтобы перевести в градусы нужно   ret.Y *= 180.0 / M_PI;
ret.Z - содержит обратный азимут в радианах, чтобы перевести в градусы нужно   ret.Я *= 180.0 / M_PI;

как то так

[QUOTE]Олег Касьянов написал:
было предусмотрено поле для хранения метрики[/QUOTE]
К слову. СУБД Postgres уникальна и своеобразна тем, что у неё есть следующие  встроенные типы данных

[IMG WIDTH=463 HEIGHT=291]http://s019.radikal.ru/i628/1705/ad/230e0948741e.png[/IMG]

Было бы неплохо, если бы Дарья Лунченко расписала как работать с метрикой
на базе типов показанных на рисунке =)

[CODE]function CreateGisTask(aMapDoc : HMAPDOC; var aParm : TTASKPARM; const aPathName : PAnsiChar) : HMAPTASK; stdcall;
begin
 result:= 0;
 .....
 КОД
 ......
end;[/CODE]где TTASKPARM[CODE] {$IFDEF CPUX64}
 type
   TTASKPARM = packed record   // ПАРАМЕТРЫ ПРИКЛАДНОЙ ЗАДАЧИ
     Language : longint;       // Код языка диалогов (1 -ENGLISH, 2 - RUSSIAN, ...)
     Reserv   : longint;       // Выравнивание           // 20/07/2016
     Resource : Int64;         // Модуль ресурсов приложения
     HelpName : PAnsiChar;         // Полное имя файла ".hlp"
     IniName  : PAnsiChar;         // Полное имя файла ".ini" приложения
     PathShell: PAnsiChar;         // Каталог приложения (exe,dll,...)
     AppName  : PAnsiChar;         // Имя приложения
     Handle   : HWND;                // Идентификатор главного окна приложения
   end;

 {$ELSE}
   TTASKPARM = packed record     // ПАРАМЕТРЫ ПРИКЛАДНОЙ ЗАДАЧИ
     Language  : longint;      // Код языка диалогов (1 -ENGLISH, 2 - RUSSIAN, ...)
     Resource  : longint;      // Модуль ресурсов приложения
     HelpName  : PAnsiChar;        // Полное имя файла ".hlp"
     IniName   : PAnsiChar;        // Полное имя файла ".ini" приложения
     PathShell : PAnsiChar;        // Каталог приложения (exe,dll,...)
     AppName   : PAnsiChar;        // Имя приложения
     Handle    : HWND;         // Идентификатор окна приложения
   end;
{$ENDIF}[/CODE]

[QUOTE]MapFind1.MapSelect.Excode[-1,OL_SQUARE] := false;
MapFind1.MapSelect.Excode[10000017,OL_SQUARE] := true;[/QUOTE]
Вы выбрали и перебираете ТОЛЬКО площадные объекты с кодом 10000017
Уточняющий вопрос, граница (которая не ищется) тоже площадной объект с кодом 10000017 ?