( sebist | 2021. 06. 24., cs – 08:32 )

Most nem tudtam megtalálni azt a leírást(valószínűleg a változónevek is azt követik), ami alapján a lenti jó régi C kód született, de talán használható belőle valami.

long double rad ( long double x)
{
    return (x/180.0*M_PI);
}    

long double deg ( long double x)
{
    return (x/M_PI*180.0);
}    

struct WGS84 {
    long double fi,la;
};
    
    
struct WGS84 eov2wgs (long double EOVY, long double EOVX)
{

struct WGS84 ret;
long double a1=6378160,b1=6356774.52,a2=6378137,b2=6356752.31;
long double dX=52.684,dY=-71.194,dZ=-13.975,eX=0.312,eY=0.1063,eZ=0.3729,k=1.0191e-6;
long double FI2,LA2,st,ct,teta,P,ev2,e22,f2,Zv,Yv,Xv;
long double h=0; // magassag
long double Z,Y,X,N,e2,f,fi,la;

long double C1=1.0007197049;
long double E1=rad(19.048571778);

long double G1=rad(47.1);

long double H1=6379296.419;

long double I1=47.0+7.0/60.0+20.0578/3600.0;
long double J1=rad(I1);

long double K1=EOVX-200000;
long double L1=EOVY-650000;

long double M1=2.0*(atan(exp(K1/H1))-M_PI/4.0);
long double N1=L1/H1;

long double O1=47+1.0/6.0;
long double P1=asin(cos(G1)*sin(M1)+sin(G1)*cos(M1)*cos(N1));

long double Q1=asin(sin(N1)*cos(M1)/cos(P1));

//R1
long double T1=rad(O1);

long double U1=6378160;
long double V1=6356774.516;

long double W1=(U1*U1-V1*V1)*cos(T1)*cos(T1)/V1/V1;    
long double X1=180*3600/M_PI;

long double Y1=sqrt(1.0+W1);
long double Z1=1.5*W1*tan(T1)/X1;

long double AA1=0.5*W1*(-1.0+tan(T1)*tan(T1)-W1+5*W1*tan(T1))/Y1/X1/X1;
long double S1=(P1-J1)*X1;

fi=(T1+S1*Y1/X1-S1*S1*Z1/X1+S1*S1*S1*AA1/X1);
la=(E1+Q1/C1);
    
f=(a1-b1)/a1;
e2=2*f-f*f;

N=a1/(sqrt(1-e2*sinl(fi)*sinl(fi)));
X=(N+h)*cos(fi)*cosl(la);

Y=(N+h)*cosl(fi)*sinl(la);
Z=(N*(1-e2)+h)*sinl(fi);
    
Xv=dX+(1+k)*( X+rad(eZ/3600)*Y-rad(eY/3600)*Z);
Yv=dY+(1+k)*(-X*rad(eZ/3600)+Y+Z*rad(eX/3600));
Zv=dZ+(1+k)*( X*rad(eY/3600)-Y*rad(eX/3600)+Z);

f2=(a2-b2)/a2;
e22=2*f2-f2*f2;
ev2=(a2*a2-b2*b2)/b2/b2;    

P=sqrt(Xv*Xv+Yv*Yv);
teta=atan2l(Zv*a2,P*b2);

ct=cosl(teta);
st=sinl(teta);
FI2=atan2l(Zv+ev2*b2*st*st*st,P-e22*a2*ct*ct*ct);
LA2=atan2l(Yv,Xv);

ret.fi=deg(FI2);
ret.la=deg(LA2);
        
return ret;

}