Curvilinear coordinate system around body of revolution In Boundary-Layer Theory by Schlichting he gives the boundary-layer equations for a body of revolution according to the paper by Boltze$^1$. Unfortunately, this paper is in German. He apparently uses the following local curvilinear coordinate system to derive the boundary-layer equations:

I'm not specifically interested in deriving the boundary-layer equations for a body of revolution, however, I am working on a problem with a very similar geometry and am having trouble wrapping my head around the coordinate system.
Question:
For the locally orthogonal basis vectors $\mathbf{e}_x$, $\mathbf{e}_y$, and $\mathbf{e}_\theta$, what are the gradient, divergence, and Laplacian operators?
This is what I have tried:
$$x = \hat{x}$$
$$y = \hat{y}$$
$$r(x)\cos\theta = \hat{z}$$
where $\hat{x}$, $\hat{y}$, and $\hat{z}$ are coordinates in a Cartesian system with the same origin. Lamé's coefficients are then:
$$h_x = \left\vert\left(\frac{\partial \hat{x}}{\partial x}, \frac{\partial \hat{y}}{\partial x}, \frac{\partial \hat{z}}{\partial x}\right)\right\vert = 1$$
$$h_y = \left\vert\left(\frac{\partial \hat{x}}{\partial y}, \frac{\partial \hat{y}}{\partial y}, \frac{\partial \hat{z}}{\partial y}\right)\right\vert = 1$$
$$h_\theta = \left\vert\left(\frac{\partial \hat{x}}{\partial \theta}, \frac{\partial \hat{y}}{\partial \theta}, \frac{\partial \hat{z}}{\partial \theta}\right)\right\vert = \left\vert r(x)\sin\theta\right\vert$$
With these coefficients it is straightforward to define the operators since this system is orthogonal (for example here are the formulas).
I do not think that the results above are correct though since it compares to a local Cartesian system at each point. What are the correct transformations?


*

*Boltze, Ernst. Grenzschichten an Rotationskorpern in Flussigkeiten mit kleiner Reibung. Georg-August-Universitat zu Gottingen., 1908.


UPDATE:
For any future Googlers, the correct transformation is: (based on Fig. 4 of @Floris's answer, $\xi$ is aligned with $U_\infty$ in the picture above)
$$\xi = x\cos\phi - y\sin\phi$$
$$\eta = -(r(\xi) + y\cos\phi + x\sin\phi)\sin\theta$$
$$\zeta = (r(\xi) + y\cos\phi + x\sin\phi)\cos\theta$$
where $\phi$ is the local angle of the body with respect to the axis of rotation, i.e. $\phi = arctan(\frac{dr}{dx})$. The paper by Boltze assumes this angle is small. The Jacobian is then:
$$J = r(\xi) + y\cos\phi + x\sin\phi$$
which is equivalent to Boltze's value for small $\phi$ and $y \ll r(\xi)$. And the Lamé coefficients are:
$$h_x = 1$$
$$h_y = 1$$
$$h_\theta = r(\xi) + y\cos\phi + x\sin\phi$$
 A: I found the original text at http://babel.hathitrust.org/cgi/pt?id=uc1.b2619178 (page 5 onwards). The relevant pages and key translation (each time I place the page before the translation - and I will repeat snapshots of equations in the translation where appropriate):

1) Contrary to the work by H Blasius, in this work we treat the problem of a body of revolution with its axis along the direction of the current. It is convenient to use the following coordinate system for the boundary layer (figure 2)

x = arc length along meridian curve
y = normal distance from the surface of the body
r = vertical distance to axis of rotation. r is a function of x and y, 
   but ultimately is only used as a way to describe the radius of the body of rotation.

2) We must transform the basic hydrodynamic equations into this framework. In the Cartesion framework, they are (in vector notation) 

Transformation of the acceleration components $(\upsilon\nabla)\upsilon$ results in components like $u\frac{du}{dx}$, $v\frac{du}{dy}$ etc. as well as centrifugal forces like $\frac{u.v}{r}$ and $\frac{u^2}{r}$. Those components that appear in the expression for x components may be omitted since they are of order $s$ compared to other components like $u\frac{du}{cx}+v\frac{du}{dy} \sim 1$

The resulting shear force, according to figure 3, is

So

Thus, the following is left from the dynamic equation for the x direction:

Estimating the order of magnitude of the Y equation yields $\frac{\partial p}{\partial y} \sim \frac{u^2}{r} \sim 1$
From this we conclude that $\frac{\partial p}{\partial x}$ in the boundary layer differs from $\frac{\partial p}{\partial x}$ in the potential flow field by an order of $s$, which we will neglect since we only preserve components $\sim 1$. Then we can use $\frac{\partial p}{\partial x}$ that was computed in the potential flow field for the boundary layer. With that, the y equation has been fully exhausted for our use - it could only serve for more subtle calculations of $p$ in the boundary layer.


*We want to write the continuity equation:



Where for a moment we consider $\xi, \eta$, and $\zeta$ as the rectangular coordinates. We then take the very general transformation
$$\xi=\xi(x,y,z); \eta=\eta(x,y,z); \zeta=\zeta(x,y,z);$$
By elementary computation, one finds the new continuity equation, in which :

The expression thus found is quite generally valid.
=====
I did the best I could. It doesn't answer your question exactly but I hope you can take it from here. Feel free to ask questions in the comments.
