От условных прямоугольных координат переходим прямоугольным проекции Г

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
1. От условных прямоугольных координат переходим к прямоугольным в проекции Гаусса-Крюгера.
2. От прямоугольных координат в проекции Гаусса-Крюгера переходим к геодезическим координатам (широта и удаление от осевого меридиана) Здесь можно применять любые формулы, хоть гостовские, хоть из Морозова, или из библиотеки proj.4 (proj_etmerc.c или proj_etmerc_optimized.c) или из  GeographicLib
у меня это реализовано так:
см. модуль uMercator
3. Зная долготу осевого меридиана, вычисляем долготу от начального меридиана.
4. По B,L,H вычисляем декартовы геоцентрические координаты X,Y,Z.
5. С помощью 7 параметров переходим от [X,Y,Z] 42  к [X,Y,Z] WGS-84
6. 6. От [X,Y,Z] WGS-84  переходим к B,L,H в WGS-84
см. модуль uGGConvert
>5. С помощью 7 параметров переходим от [X,Y,Z]42 к [X,Y,Z]WGS-84 На Паскале этот кусочек выглядит так
type TajsDatum = record
DX, DY, DZ: Extended;
scale: Extended;
Rx, Ry, Rz: Extended;
end;
ToWGS84: TajsDatum;
with ToWGS84 do
begin
X3D := scale * (+X3D_42 + Rz * Y3D_42 - Ry * Z3D_42) + Dx;
Y3D := scale * (-Rz * X3D_42 + Y3D_42 + Rx * Z3D_42) + Dy;
Z3D := scale * (+Ry * X3D_42 - Rx * Y3D_42 + Z3D_42) + Dz;
end;
Ниже дано строгое аналитическое решение для преобразования трёхмерных прямоугольных координат в геодезические.
Этот алгоритм опубликовал в 80-х годах польский астроном Казимир Барковский.
Потом этот алгоритм был описан в стандартах IERS 1996 года и стал его частью.
Его можно использовать для контроля составления программы.
procedure XYZ2BLH(X,Y,Z,SemiMajorAxis,flattening: Extended; out B,L,H: Extended);
var
R, BB, EE, F, P, Q, D, V, S, G, T: extended;
begin
R := sqrt( X*X + Y*Y );
BB := SemiMajorAxis * ( 1.0 - flattening );
if Z<0 then BB := -BB;
EE := ( (Z+BB)*BB/SemiMajorAxis-SemiMajorAxis ) / R;
F := ( (Z-BB)*BB/SemiMajorAxis+SemiMajorAxis ) / R;
P := ( EE*F+1.0 ) * 4.0/3.0;
Q := ( Sqr(EE) - Sqr(F) ) * 2.0;
D := P*Sqr(P) + Q*Q;
if D < 0
then V := 2.0 * Sqrt(-P) *
Cos( ArcCos( Q/P/Sqrt(-P) ) / 3.0 )
else begin
S := Q + Sqrt(D);
V := Exp( Ln( Abs(S) ) / 3.0 );
if V<0 then S := -V else S := V;
V := P/S - S;
V := - ( V*V*V + Q+Q ) / (P*3.0);
end;
G := ( EE + Sqrt(Sqr(EE)+V) ) / 2.0;
T := Sqrt( Sqr(G) + (F-V*G)/(G+G-EE) ) - G;
B := ArcTan( (1.0-Sqr(T)) * SemiMajorAxis/(BB*T*2.0) );
H := (R-SemiMajorAxis*T) * Cos(B)+(Z-BB)*Sin(B);
L := ArcTan2(Y,X);
end;
п.с. 
>только вот заметил,что Вы писали про систему СК42, а не СК63 (мск50), или Вы это опечатались?
Нет, это не описка. Один и тот же пункт имеет одни и те же геодезические координаты (широту, долготу и геодезическую высоту относительно отсчётного эллипсоида) во всех этих системах отсчёта.
Проекция для этих систем тоже одна и та же — Гаусса-Крюгера, отличаются только параметры проекции. Так как СК-63 и МСК-50 имеют один и тот же осевой меридиан, то координаты отличаются только на постоянный сдвиг по осям.