For such problems, I try to visualize the points in QGIS to see where they are placed.
From your parameters, I created a custom CRS with the definition:
+proj=sterea +lat_0=54.4353877827032 +lon_0=18.4514121640352 +k=0.999790760649094 +x_0=41614.2107651061 +y_0=17150.1692507701 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
and created points in WGS84 CRS and the Stereo projection. Unfortunately, they do not align. Transforming the WGS84 coordinates to your Stereo projection I get:
X: 52504.84
Y: 5709.89
Visualizing your dataset, I get the following picture:
So what went wrong? You simply swapped X and Y. "Normally", we think of Easting as X and Northing as Y. Your parameters are right, but you have to swap X0 and Y0:
+proj=sterea +lat_0=54.4353877827032 +lon_0=18.4514121640352 +k=0.999790760649094 +x_0=17150.1692507701 +y_0=41614.2107651061 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
I added the points once in WGS84 and once in the Stereo projection. Looking closer at the picture, you see an offset between some points of around 100 metres. This is the point where usually the rotations from the +towgs84 parameters come into effect. Or the central meridian was not the right one.
So now we have a proj definition which you can use with GDAL cs2cs to convert the point coordinates using a batch with the command line:
cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=sterea +lat_0=54.4353877827032 +lon_0=18.4514121640352 +k=0.999790760649094 +x_0=17150.1692507701 +y_0=41614.2107651061 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs in.txt >>out.txt
for as many points as you want.
Doing some linear regression on my own, I found this projection string best fitting for your sample data:
+proj=tmerc +lat_0=54 +lon_0=17.02 +k=1 +x_0=-75731 +y_0=-7791 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
The cs2cs command line now is:
cs2cs +init=epsg:4326 +to +proj=tmerc +lat_0=54 +lon_0=17.02 +k=1 +x_0=-75731 +y_0=-7791 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs in.txt >out.txt
Don't bother that it is not a stereographic projection anymore: For the size of your area, the difference between projection methods is much smaller than finding the right central meridian.
For the linear regression, I used a simple LibreOffice Calc sheet:
Imported your sample data, and made an input file for cs2cs from the first two columns with easting first.
As first try, I chose lat_0 and lon_0 rounded to full degrees of the sample data setting x0 and y0 to zero, ran the cs2cs command and filled the output in column E to G. The next columns compute the difference between target and cs2cs values and the absolute offset.
Below, I computed mean and standard deviation for deltaX deltaY and radius.
After that, I played a bit with the cs2cs parameters, and found that lat_0, x0 and y0 have no effect on the std dev value. So I changed lon_0 until the std dev was less than 1 meter. I copied the std dev row for each run as protocol.
The mean values with inverted sign then give the values for x0 and y0 for the final step where delta X and deltaY is less than 2 metres for all points.
You can do the same with different (but similar) projections, values for lat_0, k, ellipsoid or +towgs84. The sample points should cover your area of interest, as extrapolation will give worse results.
Final Edit
Taking the tangent point from polish UKLAD 65 Zone 3 as mentioned in the comments, I got best results with the following command line:
cs2cs +init=epsg:4326 +to +proj=sterea +lat_0=53.58333333333334 +lon_0=17.00833333333333 +k=0.9998 +x_0=-76625.57 +y_0=-54048.732 +ellps=krass +towgs84=33.4,-146.6,176,0,0,0,0 +units=m +no_defs in.txt >out.txt
Offset is less than 1 meter for all sample points, and resulting standard dev 0.32m. You could even try to use all 7 Helmert parameters to get it better.
I guess you stumbled over this note on page 53:
subsequent civilian adoption of the systems usually ignores the
zone prefix to easting. Where this is the case the formulas below do
not apply: use the standard TM formula separately for each zone
The formulas on page 54 should only be used if the zone number is written before the false Easting, which is not the case here. So I took only the formulas from p 50ff into Libre Office and got these results:
M0 0
e1 0.0016649922
mue1 0.4772587465
M1 3038553.85538554
T1 0.2701106629
C1 0.0059292372
esq 0.0066822021
D 0.0183352698
phi1 0.4793002154
nue1 6381782.36093598
rho1 6348382.77011411
phirad 0.4792124095 phi 27.4568485565
lamrad 1.4867384817 lam 85.1838402406
which returns the same lat and lon on the Everest ellipsoid as
cs2cs +proj=tmerc +lon_0=84 +lat_0=0 +k=0.9999 +a=6377276.345 +b=6356075.41314024 +x_0=500000 +units=m +towgs84=282,726,254 +no_defs +to +proj=longlat +a=6377276.345 +b=6356075.41314024 +towgs84=282,726,254 +no_defs -f "%%.6f"<MUTM88.txt >WGS84out.txt
85.183840 27.456849 1230.000000
Transforming to WGS84 returns:
85.181615 27.457128 1188.383496
which is still not exactly what you expected, but more close than your result.
Best Answer
You can't use the Helmert-transformation (with Bursa-Wolf parameter) with geographic coordinates. First transform the geographic coordinates into geocentric coordinates. see
That's 3D coordinates, so you can easily integrate the altitude value.
But consider the precision using this transformation is 1 to 8 m depending on area and parameters.
But why the altitude value is the difference to NN ?