Calculating Second Moment of Area of NACA Aerofoils I've been trying to code a Python program to calculate $I_x$, $I_y$ and $I_{xy}$ for a four digit NACA Aerofoil whose outline can be calculated as follows:
The camberline is given by:
$y_c(x)=
\begin{cases}
\frac{M}{P^2}(2Px-x^2),  & 0\le x< P \\
\frac {M}{(1-P)^2}(1-2P+2Px-x^2), & P\le x\le 1
\end{cases}$
The gradient is given by:
$\frac{dy_c}{dx}(x)=
\begin{cases}
\frac{2M}{P^2}(P-x),  & 0\le x< P \\
\frac {2M}{(1-P)^2}(P-x), & P\le x\le 1
\end{cases}$
A thickness distribution is generated by:
$y_t(x)=\frac{T}{0.2}(a_0\sqrt x+a_1x+a_2x^2+a_3x^3+a_4x^4)$
where $a_n$ is a constant.
The equation of the upper surface of the aerofoil is given by:
$x_u=x_c-y_t\sin(\theta), y_u=y_c+y_t\cos(\theta)$
and for the lower surface:
$x_l=x_c+y_t\sin(\theta), y_l=y_c-y_t\cos(\theta)$
Where $\theta=\frac{dy_c}{dx}(x)$. See here for more detail on how this works.
For calculating $I_y$ I found an MIT PDF here but I don't understand how this would translate for calculating $I_y$ and $I_{xy}$. For reference the method they have used is:
$Area (A) = \int_{0}^{c}[y_u-y_l]dx$
$\overline {y} = \frac{1}{A}\int_{0}^{c}\frac{1}{2}[{y_u}^2-{y_l}^2]dx$
$I_y = \int_{0}^{c}\frac{1}{3}[({y_u-\overline {y}})^3-({y_l-\overline {y}})^3]dx$
My best guess for calculating $I_x$ is:
$\overline {x} = \frac{1}{A}\int_{y_{min}}^{y_{max}}\frac{1}{2}[{x_u}^2-{x_l}^2]dx$
$I_x = \int_{0}^{c}\frac{1}{3}[({x_u-\overline {x}})^3-({x_l-\overline {x}})^3]dx$
For $I_{xy}$ I have no idea. As a final question using this method are the second moments of area relative to the centroid of the aerofoil?
 A: Assuming a constant density, the definitions of the moments of inertia about the origin of the coordinate system are $$I_x = \int y^2\ dA\\I_y = \int x^2\ dA\\I_{xy} = \int xy\ dA$$
Where the last one follows the more common sign convention (the inertia tensor has $-I_{xy}$ as its off-diagonal components). 
In the first version of this post, I went through the derivation of the formulas in the MIT paper. But that paper is uses a different method of expressing its curves than here. The curves here are actually parametric, and that works better with a different set of formulas. First, let me recast the curves:
Note that the curves are given by $\alpha(x) = (x_l(x), y_l(x))$ on bottom, and $\beta(x) = (x_u(x), y_u(x))$ on top. Now the "$x$" upon which they vary is the x-coordinate along the camberline, not these, which makes this notation confusing. For that reason, I am changing the variable to $t$, since it is not the $x$ value for the curves we are interested in. For our purposes it is just a parameter. And similarly, to highlight the "$y$" of interest in this calculation, let me just rename $y_c$ as $u$, and $y_t$ as $v$. So:
$$ 0 \le t \le 1\\ \ \\u(t) =\begin{cases}
\frac M{P^2}(2Pt-t^2),  & 0\le t< P \\
\frac M{(1-P)^2}(1-2P+2Pt-t^2), & P\le t\le 1
\end{cases}\\
\dot u(t) = \frac{du}{dt}(t)=
\begin{cases}
\frac{2M}{P^2}(P-t),  & 0\le t< P \\
\frac {2M}{(1-P)^2}(P-t), & P\le t\le 1
\end{cases}\\
v(t)=\frac{T}{0.2}(a_0\sqrt t+a_1t+a_2t^2+a_3t^3+a_4t^4)$$
It will be helpful also to calculate $$w(t) = \frac 1{1+\dot u^2}$$
because then $\cos \theta = w(t)$ and $\sin \theta = \dot u(t)w(t)$, and so the parametric curves forming the boundaries of the airfoil are
$$\alpha(t) = (t + v(t)\dot u(t)w(t), u(t) - v(t)w(t))$$
below and
$$\beta(t) = (t - v(t)\dot u(t)w(t), u(t) + v(t)w(t))$$
above.
The full curve in the counterclockwise direction can be defined by $$\gamma : [0,2] \to \Bbb R^2$$ with
$$\gamma(t) = \begin{cases} \alpha(t) & 0 \le t \le 1\\\beta(2 - t) & 1\le t \le 2\end{cases}$$ 
With this parametric boundary we can use Green's theorem to simplify the various integrals (see formulas 10 through 16), Letting $(x(t), y(t)) = \gamma(t)$,
$$A = \int_0^2 x\dot y\ dt\\
\bar x = -\frac 1{2A}\int_0^2 y^2\dot x\ dt\\
\bar y = \frac 1{2A}\int_0^2 x^2\dot y\ dt\\
I_x = -\frac 13\int_0^2 (y-\bar y)^3\dot x\ dt\\
I_y = \frac 13\int_0^2 (x-\bar x)^3\dot y\ dt\\
I_{xy}= -\frac 12 \int_0^2 (y-\bar y)^2(x-\bar x)\dot x\ dt$$
Where $\dot x = \frac {dx}{dt}, \dot y = \frac {dy}{dt}$, which can be calculated from the formulas above. For the integrals, you'll want to do numerical integration rather than trying to simplify them directly. I've modified the moment integrals to be around the centroid, rather than the origin as in the Wolfram link.
