/src/PROJ/src/transformations/molodensky.cpp

Source
/***********************************************************************

                  (Abridged) Molodensky Transform

                    Kristian Evers, 2017-07-07

************************************************************************

    Implements the (abridged) Molodensky transformations for 2D and 3D
    data.

    Primarily useful for implementation of datum shifts in transformation
    pipelines.

    The code in this file is mostly based on

        The Standard and Abridged Molodensky Coordinate Transformation
        Formulae, 2004, R.E. Deakin,
        http://www.mygeodesy.id.au/documents/Molodensky%20V2.pdf



************************************************************************
* Copyright (c) 2017, Kristian Evers / SDFE
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
*
***********************************************************************/

#include <errno.h>
#include <math.h>

#include "proj.h"
#include "proj_internal.h"

PROJ_HEAD(molodensky, "Molodensky transform");

static PJ_XYZ pj_molodensky_forward_3d(PJ_LPZ lpz, PJ *P);
static PJ_LPZ pj_molodensky_reverse_3d(PJ_XYZ xyz, PJ *P);

namespace { // anonymous namespace
struct pj_opaque_molodensky {
    double dx;
    double dy;
    double dz;
    double da;
    double df;
    int abridged;
};
} // anonymous namespace

static double RN(double a, double es, double phi) {
    /**********************************************************
        N(phi) - prime vertical radius of curvature
        -------------------------------------------

        This is basically the same function as in PJ_cart.c
        should probably be refactored into it's own file at some
        point.

    **********************************************************/
    double s = sin(phi);
    if (es == 0)
        return a;

    return a / sqrt(1 - es * s * s);
}

static double RM(double a, double es, double phi) {
    /**********************************************************
        M(phi) - Meridian radius of curvature
        -------------------------------------

        Source:

            E.J Krakiwsky & D.B. Thomson, 1974,
            GEODETIC POSITION COMPUTATIONS,

            Fredericton NB, Canada:
            University of New Brunswick,
            Department of Geodesy and Geomatics Engineering,
            Lecture Notes No. 39,
            99 pp.

            http://www2.unb.ca/gge/Pubs/LN39.pdf

    **********************************************************/
    double s = sin(phi);
    if (es == 0)
        return a;

    /* eq. 13a */
    if (phi == 0)
        return a * (1 - es);

    /* eq. 13b */
    if (fabs(phi) == M_PI_2)
        return a / sqrt(1 - es);

    /* eq. 13 */
    return (a * (1 - es)) / pow(1 - es * s * s, 1.5);
}

static PJ_LPZ calc_standard_params(PJ_LPZ lpz, PJ *P) {
    struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *)P->opaque;
    double dphi, dlam, dh;

    /* sines and cosines */
    double slam = sin(lpz.lam);
    double clam = cos(lpz.lam);
    double sphi = sin(lpz.phi);
    double cphi = cos(lpz.phi);

    /* ellipsoid parameters and differences */
    double f = P->f, a = P->a;
    double dx = Q->dx, dy = Q->dy, dz = Q->dz;
    double da = Q->da, df = Q->df;

    /* ellipsoid radii of curvature */
    double rho = RM(a, P->es, lpz.phi);
    double nu = RN(a, P->es, lpz.phi);

    /* delta phi */
    dphi = (-dx * sphi * clam) - (dy * sphi * slam) + (dz * cphi) +
           ((nu * P->es * sphi * cphi * da) / a) +
           (sphi * cphi * (rho / (1 - f) + nu * (1 - f)) * df);
    const double dphi_denom = rho + lpz.z;
    if (dphi_denom == 0.0) {
        lpz.lam = HUGE_VAL;
        return lpz;
    }
    dphi /= dphi_denom;

    /* delta lambda */
    const double dlam_denom = (nu + lpz.z) * cphi;
    if (dlam_denom == 0.0) {
        lpz.lam = HUGE_VAL;
        return lpz;
    }
    dlam = (-dx * slam + dy * clam) / dlam_denom;

    /* delta h */
    dh = dx * cphi * clam + dy * cphi * slam + dz * sphi - (a / nu) * da +
         nu * (1 - f) * sphi * sphi * df;

    lpz.phi = dphi;
    lpz.lam = dlam;
    lpz.z = dh;

    return lpz;
}

static PJ_LPZ calc_abridged_params(PJ_LPZ lpz, PJ *P) {
    struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *)P->opaque;
    double dphi, dlam, dh;

    /* sines and cosines */
    double slam = sin(lpz.lam);
    double clam = cos(lpz.lam);
    double sphi = sin(lpz.phi);
    double cphi = cos(lpz.phi);

    /* ellipsoid parameters and differences */
    double dx = Q->dx, dy = Q->dy, dz = Q->dz;
    double da = Q->da, df = Q->df;
    double adffda = (P->a * df + P->f * da);

    /* delta phi */
    dphi = -dx * sphi * clam - dy * sphi * slam + dz * cphi +
           adffda * sin(2 * lpz.phi);
    dphi /= RM(P->a, P->es, lpz.phi);

    /* delta lambda */
    dlam = -dx * slam + dy * clam;
    const double dlam_denom = RN(P->a, P->es, lpz.phi) * cphi;
    if (dlam_denom == 0.0) {
        lpz.lam = HUGE_VAL;
        return lpz;
    }
    dlam /= dlam_denom;

    /* delta h */
    dh = dx * cphi * clam + dy * cphi * slam + dz * sphi - da +
         adffda * sphi * sphi;

    /* offset coordinate */
    lpz.phi = dphi;
    lpz.lam = dlam;
    lpz.z = dh;

    return lpz;
}

static PJ_XY pj_molodensky_forward_2d(PJ_LP lp, PJ *P) {
    PJ_COORD point = {{0, 0, 0, 0}};

    point.lp = lp;
    // Assigning in 2 steps avoids cppcheck warning
    // "Overlapping read/write of union is undefined behavior"
    // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
    const auto xyz = pj_molodensky_forward_3d(point.lpz, P);
    point.xyz = xyz;

    return point.xy;
}

static PJ_LP pj_molodensky_reverse_2d(PJ_XY xy, PJ *P) {
    PJ_COORD point = {{0, 0, 0, 0}};

    point.xy = xy;
    point.xyz.z = 0;
    // Assigning in 2 steps avoids cppcheck warning
    // "Overlapping read/write of union is undefined behavior"
    // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
    const auto lpz = pj_molodensky_reverse_3d(point.xyz, P);
    point.lpz = lpz;

    return point.lp;
}

static PJ_XYZ pj_molodensky_forward_3d(PJ_LPZ lpz, PJ *P) {
    struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *)P->opaque;
    PJ_COORD point = {{0, 0, 0, 0}};

    point.lpz = lpz;

    /* calculate parameters depending on the mode we are in */
    if (Q->abridged) {
        lpz = calc_abridged_params(lpz, P);
    } else {
        lpz = calc_standard_params(lpz, P);
    }
    if (lpz.lam == HUGE_VAL) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return proj_coord_error().xyz;
    }

    /* offset coordinate */
    point.lpz.phi += lpz.phi;
    point.lpz.lam += lpz.lam;
    point.lpz.z += lpz.z;

    return point.xyz;
}

static void pj_molodensky_forward_4d(PJ_COORD &obs, PJ *P) {
    // Assigning in 2 steps avoids cppcheck warning
    // "Overlapping read/write of union is undefined behavior"
    // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
    const auto xyz = pj_molodensky_forward_3d(obs.lpz, P);
    obs.xyz = xyz;
}

static PJ_LPZ pj_molodensky_reverse_3d(PJ_XYZ xyz, PJ *P) {
    struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *)P->opaque;
    PJ_COORD point = {{0, 0, 0, 0}};
    PJ_LPZ lpz;

    /* calculate parameters depending on the mode we are in */
    point.xyz = xyz;
    if (Q->abridged)
        lpz = calc_abridged_params(point.lpz, P);
    else
        lpz = calc_standard_params(point.lpz, P);

    if (lpz.lam == HUGE_VAL) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return proj_coord_error().lpz;
    }

    /* offset coordinate */
    point.lpz.phi -= lpz.phi;
    point.lpz.lam -= lpz.lam;
    point.lpz.z -= lpz.z;

    return point.lpz;
}

static void pj_molodensky_reverse_4d(PJ_COORD &obs, PJ *P) {
    // Assigning in 2 steps avoids cppcheck warning
    // "Overlapping read/write of union is undefined behavior"
    // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
    const auto lpz = pj_molodensky_reverse_3d(obs.xyz, P);
    obs.lpz = lpz;
}

PJ *PJ_TRANSFORMATION(molodensky, 1) {
    struct pj_opaque_molodensky *Q = static_cast<struct pj_opaque_molodensky *>(
        calloc(1, sizeof(struct pj_opaque_molodensky)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = (void *)Q;

    P->fwd4d = pj_molodensky_forward_4d;
    P->inv4d = pj_molodensky_reverse_4d;
    P->fwd3d = pj_molodensky_forward_3d;
    P->inv3d = pj_molodensky_reverse_3d;
    P->fwd = pj_molodensky_forward_2d;
    P->inv = pj_molodensky_reverse_2d;

    P->left = PJ_IO_UNITS_RADIANS;
    P->right = PJ_IO_UNITS_RADIANS;

    /* read args */
    if (!pj_param(P->ctx, P->params, "tdx").i) {
        proj_log_error(P, _("missing dx"));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }
    Q->dx = pj_param(P->ctx, P->params, "ddx").f;

    if (!pj_param(P->ctx, P->params, "tdy").i) {
        proj_log_error(P, _("missing dy"));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }
    Q->dy = pj_param(P->ctx, P->params, "ddy").f;

    if (!pj_param(P->ctx, P->params, "tdz").i) {
        proj_log_error(P, _("missing dz"));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }
    Q->dz = pj_param(P->ctx, P->params, "ddz").f;

    if (!pj_param(P->ctx, P->params, "tda").i) {
        proj_log_error(P, _("missing da"));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }
    Q->da = pj_param(P->ctx, P->params, "dda").f;

    if (!pj_param(P->ctx, P->params, "tdf").i) {
        proj_log_error(P, _("missing df"));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }
    Q->df = pj_param(P->ctx, P->params, "ddf").f;

    Q->abridged = pj_param(P->ctx, P->params, "tabridged").i;

    return P;
}

Coverage Report

Created: 2026-09-28 07:23

Line	Count	Source
1		/***********************************************************************
2
3		(Abridged) Molodensky Transform
4
5		Kristian Evers, 2017-07-07
6
7		************************************************************************
8
9		Implements the (abridged) Molodensky transformations for 2D and 3D
10		data.
11
12		Primarily useful for implementation of datum shifts in transformation
13		pipelines.
14
15		The code in this file is mostly based on
16
17		The Standard and Abridged Molodensky Coordinate Transformation
18		Formulae, 2004, R.E. Deakin,
19		http://www.mygeodesy.id.au/documents/Molodensky%20V2.pdf
20
21
22
23		************************************************************************
24		* Copyright (c) 2017, Kristian Evers / SDFE
25		*
26		* Permission is hereby granted, free of charge, to any person obtaining a
27		* copy of this software and associated documentation files (the "Software"),
28		* to deal in the Software without restriction, including without limitation
29		* the rights to use, copy, modify, merge, publish, distribute, sublicense,
30		* and/or sell copies of the Software, and to permit persons to whom the
31		* Software is furnished to do so, subject to the following conditions:
32		*
33		* The above copyright notice and this permission notice shall be included
34		* in all copies or substantial portions of the Software.
35		*
36		* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
37		* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
38		* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
39		* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
40		* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
41		* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
42		* DEALINGS IN THE SOFTWARE.
43		*
44		***********************************************************************/
45
46		#include <errno.h>
47		#include <math.h>
48
49		#include "proj.h"
50		#include "proj_internal.h"
51
52		PROJ_HEAD(molodensky, "Molodensky transform");
53
54		static PJ_XYZ pj_molodensky_forward_3d(PJ_LPZ lpz, PJ *P);
55		static PJ_LPZ pj_molodensky_reverse_3d(PJ_XYZ xyz, PJ *P);
56
57		namespace { // anonymous namespace
58		struct pj_opaque_molodensky {
59		double dx;
60		double dy;
61		double dz;
62		double da;
63		double df;
64		int abridged;
65		};
66		} // anonymous namespace
67
68	0	static double RN(double a, double es, double phi) {
69		/**********************************************************
70		N(phi) - prime vertical radius of curvature
71		-------------------------------------------
72
73		This is basically the same function as in PJ_cart.c
74		should probably be refactored into it's own file at some
75		point.
76
77		**********************************************************/
78	0	double s = sin(phi);
79	0	if (es == 0)
80	0	return a;
81
82	0	return a / sqrt(1 - es * s * s);
83	0	}
84
85	0	static double RM(double a, double es, double phi) {
86		/**********************************************************
87		M(phi) - Meridian radius of curvature
88		-------------------------------------
89
90		Source:
91
92		E.J Krakiwsky & D.B. Thomson, 1974,
93		GEODETIC POSITION COMPUTATIONS,
94
95		Fredericton NB, Canada:
96		University of New Brunswick,
97		Department of Geodesy and Geomatics Engineering,
98		Lecture Notes No. 39,
99		99 pp.
100
101		http://www2.unb.ca/gge/Pubs/LN39.pdf
102
103		**********************************************************/
104	0	double s = sin(phi);
105	0	if (es == 0)
106	0	return a;
107
108		/* eq. 13a */
109	0	if (phi == 0)
110	0	return a * (1 - es);
111
112		/* eq. 13b */
113	0	if (fabs(phi) == M_PI_2)
114	0	return a / sqrt(1 - es);
115
116		/* eq. 13 */
117	0	return (a * (1 - es)) / pow(1 - es * s * s, 1.5);
118	0	}
119
120	0	static PJ_LPZ calc_standard_params(PJ_LPZ lpz, PJ *P) {
121	0	struct pj_opaque_molodensky Q = (struct pj_opaque_molodensky )P->opaque;
122	0	double dphi, dlam, dh;
123
124		/* sines and cosines */
125	0	double slam = sin(lpz.lam);
126	0	double clam = cos(lpz.lam);
127	0	double sphi = sin(lpz.phi);
128	0	double cphi = cos(lpz.phi);
129
130		/* ellipsoid parameters and differences */
131	0	double f = P->f, a = P->a;
132	0	double dx = Q->dx, dy = Q->dy, dz = Q->dz;
133	0	double da = Q->da, df = Q->df;
134
135		/* ellipsoid radii of curvature */
136	0	double rho = RM(a, P->es, lpz.phi);
137	0	double nu = RN(a, P->es, lpz.phi);
138
139		/* delta phi */
140	0	dphi = (-dx * sphi * clam) - (dy * sphi * slam) + (dz * cphi) +
141	0	((nu * P->es * sphi * cphi * da) / a) +
142	0	(sphi * cphi * (rho / (1 - f) + nu * (1 - f)) * df);
143	0	const double dphi_denom = rho + lpz.z;
144	0	if (dphi_denom == 0.0) {
145	0	lpz.lam = HUGE_VAL;
146	0	return lpz;
147	0	}
148	0	dphi /= dphi_denom;
149
150		/* delta lambda */
151	0	const double dlam_denom = (nu + lpz.z) * cphi;
152	0	if (dlam_denom == 0.0) {
153	0	lpz.lam = HUGE_VAL;
154	0	return lpz;
155	0	}
156	0	dlam = (-dx * slam + dy * clam) / dlam_denom;
157
158		/* delta h */
159	0	dh = dx * cphi * clam + dy * cphi * slam + dz * sphi - (a / nu) * da +
160	0	nu * (1 - f) * sphi * sphi * df;
161
162	0	lpz.phi = dphi;
163	0	lpz.lam = dlam;
164	0	lpz.z = dh;
165
166	0	return lpz;
167	0	}
168
169	0	static PJ_LPZ calc_abridged_params(PJ_LPZ lpz, PJ *P) {
170	0	struct pj_opaque_molodensky Q = (struct pj_opaque_molodensky )P->opaque;
171	0	double dphi, dlam, dh;
172
173		/* sines and cosines */
174	0	double slam = sin(lpz.lam);
175	0	double clam = cos(lpz.lam);
176	0	double sphi = sin(lpz.phi);
177	0	double cphi = cos(lpz.phi);
178
179		/* ellipsoid parameters and differences */
180	0	double dx = Q->dx, dy = Q->dy, dz = Q->dz;
181	0	double da = Q->da, df = Q->df;
182	0	double adffda = (P->a * df + P->f * da);
183
184		/* delta phi */
185	0	dphi = -dx * sphi * clam - dy * sphi * slam + dz * cphi +
186	0	adffda * sin(2 * lpz.phi);
187	0	dphi /= RM(P->a, P->es, lpz.phi);
188
189		/* delta lambda */
190	0	dlam = -dx * slam + dy * clam;
191	0	const double dlam_denom = RN(P->a, P->es, lpz.phi) * cphi;
192	0	if (dlam_denom == 0.0) {
193	0	lpz.lam = HUGE_VAL;
194	0	return lpz;
195	0	}
196	0	dlam /= dlam_denom;
197
198		/* delta h */
199	0	dh = dx * cphi * clam + dy * cphi * slam + dz * sphi - da +
200	0	adffda * sphi * sphi;
201
202		/* offset coordinate */
203	0	lpz.phi = dphi;
204	0	lpz.lam = dlam;
205	0	lpz.z = dh;
206
207	0	return lpz;
208	0	}
209
210	0	static PJ_XY pj_molodensky_forward_2d(PJ_LP lp, PJ *P) {
211	0	PJ_COORD point = {{0, 0, 0, 0}};
212
213	0	point.lp = lp;
214		// Assigning in 2 steps avoids cppcheck warning
215		// "Overlapping read/write of union is undefined behavior"
216		// Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
217	0	const auto xyz = pj_molodensky_forward_3d(point.lpz, P);
218	0	point.xyz = xyz;
219
220	0	return point.xy;
221	0	}
222
223	0	static PJ_LP pj_molodensky_reverse_2d(PJ_XY xy, PJ *P) {
224	0	PJ_COORD point = {{0, 0, 0, 0}};
225
226	0	point.xy = xy;
227	0	point.xyz.z = 0;
228		// Assigning in 2 steps avoids cppcheck warning
229		// "Overlapping read/write of union is undefined behavior"
230		// Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
231	0	const auto lpz = pj_molodensky_reverse_3d(point.xyz, P);
232	0	point.lpz = lpz;
233
234	0	return point.lp;
235	0	}
236
237	0	static PJ_XYZ pj_molodensky_forward_3d(PJ_LPZ lpz, PJ *P) {
238	0	struct pj_opaque_molodensky Q = (struct pj_opaque_molodensky )P->opaque;
239	0	PJ_COORD point = {{0, 0, 0, 0}};
240
241	0	point.lpz = lpz;
242
243		/* calculate parameters depending on the mode we are in */
244	0	if (Q->abridged) {
245	0	lpz = calc_abridged_params(lpz, P);
246	0	} else {
247	0	lpz = calc_standard_params(lpz, P);
248	0	}
249	0	if (lpz.lam == HUGE_VAL) {
250	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
251	0	return proj_coord_error().xyz;
252	0	}
253
254		/* offset coordinate */
255	0	point.lpz.phi += lpz.phi;
256	0	point.lpz.lam += lpz.lam;
257	0	point.lpz.z += lpz.z;
258
259	0	return point.xyz;
260	0	}
261
262	0	static void pj_molodensky_forward_4d(PJ_COORD &obs, PJ *P) {
263		// Assigning in 2 steps avoids cppcheck warning
264		// "Overlapping read/write of union is undefined behavior"
265		// Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
266	0	const auto xyz = pj_molodensky_forward_3d(obs.lpz, P);
267	0	obs.xyz = xyz;
268	0	}
269
270	0	static PJ_LPZ pj_molodensky_reverse_3d(PJ_XYZ xyz, PJ *P) {
271	0	struct pj_opaque_molodensky Q = (struct pj_opaque_molodensky )P->opaque;
272	0	PJ_COORD point = {{0, 0, 0, 0}};
273	0	PJ_LPZ lpz;
274
275		/* calculate parameters depending on the mode we are in */
276	0	point.xyz = xyz;
277	0	if (Q->abridged)
278	0	lpz = calc_abridged_params(point.lpz, P);
279	0	else
280	0	lpz = calc_standard_params(point.lpz, P);
281
282	0	if (lpz.lam == HUGE_VAL) {
283	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
284	0	return proj_coord_error().lpz;
285	0	}
286
287		/* offset coordinate */
288	0	point.lpz.phi -= lpz.phi;
289	0	point.lpz.lam -= lpz.lam;
290	0	point.lpz.z -= lpz.z;
291
292	0	return point.lpz;
293	0	}
294
295	0	static void pj_molodensky_reverse_4d(PJ_COORD &obs, PJ *P) {
296		// Assigning in 2 steps avoids cppcheck warning
297		// "Overlapping read/write of union is undefined behavior"
298		// Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
299	0	const auto lpz = pj_molodensky_reverse_3d(obs.xyz, P);
300	0	obs.lpz = lpz;
301	0	}
302
303	1	PJ *PJ_TRANSFORMATION(molodensky, 1) {
304	1	struct pj_opaque_molodensky Q = static_cast<struct pj_opaque_molodensky >(
305	1	calloc(1, sizeof(struct pj_opaque_molodensky)));
306	1	if (nullptr == Q)
307	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
308	1	P->opaque = (void *)Q;
309
310	1	P->fwd4d = pj_molodensky_forward_4d;
311	1	P->inv4d = pj_molodensky_reverse_4d;
312	1	P->fwd3d = pj_molodensky_forward_3d;
313	1	P->inv3d = pj_molodensky_reverse_3d;
314	1	P->fwd = pj_molodensky_forward_2d;
315	1	P->inv = pj_molodensky_reverse_2d;
316
317	1	P->left = PJ_IO_UNITS_RADIANS;
318	1	P->right = PJ_IO_UNITS_RADIANS;
319
320		/* read args */
321	1	if (!pj_param(P->ctx, P->params, "tdx").i) {
322	1	proj_log_error(P, _("missing dx"));
323	1	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
324	1	}
325	0	Q->dx = pj_param(P->ctx, P->params, "ddx").f;
326
327	0	if (!pj_param(P->ctx, P->params, "tdy").i) {
328	0	proj_log_error(P, _("missing dy"));
329	0	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
330	0	}
331	0	Q->dy = pj_param(P->ctx, P->params, "ddy").f;
332
333	0	if (!pj_param(P->ctx, P->params, "tdz").i) {
334	0	proj_log_error(P, _("missing dz"));
335	0	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
336	0	}
337	0	Q->dz = pj_param(P->ctx, P->params, "ddz").f;
338
339	0	if (!pj_param(P->ctx, P->params, "tda").i) {
340	0	proj_log_error(P, _("missing da"));
341	0	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
342	0	}
343	0	Q->da = pj_param(P->ctx, P->params, "dda").f;
344
345	0	if (!pj_param(P->ctx, P->params, "tdf").i) {
346	0	proj_log_error(P, _("missing df"));
347	0	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
348	0	}
349	0	Q->df = pj_param(P->ctx, P->params, "ddf").f;
350
351	0	Q->abridged = pj_param(P->ctx, P->params, "tabridged").i;
352
353	0	return P;
354	0	}