Solution (source code)

= Solution

Put the slit centres at $x_j=\pm D/2,\pm3D/2$. Each slit has the same width $a$ and illumination. Its contribution to the <Fraunhofer diffraction> <complex amplitude> is
$$
\int_{x_j-a/2}^{x_j+a/2}e^{-iqx}\,dx
=a\operatorname{sinc}(qa/2)e^{-iqx_j}.
$$
The four <phase factors> pair into cosines:
$$
\sum_j e^{-iqx_j}=e^{iqD/2}+e^{-iqD/2}+e^{3iqD/2}+e^{-3iqD/2}
=2\cos(qD/2)+2\cos(3qD/2).
$$
Therefore the <four-slit Fraunhofer diffraction> field is
$$
\boxed{\psi(q)\propto2a[\cos(qD/2)+\cos(3qD/2)]\operatorname{sinc}(qa/2),\qquad
q=\frac{2\pi\sin\theta}{\lambda}.}
$$
The common factor $2a$ can be absorbed in the proportionality constant, giving the requested form. In particular $\psi(0)\propto4a$, as expected from four equal in-phase slits. Normalizing by the central intensity gives
$$
\boxed{\frac{I(q)}{I(0)}=\operatorname{sinc}^2(qa/2)\cos^2(qD/2)\cos^2(qD).}
$$
The last form follows from $\cos x+\cos3x=2\cos2x\cos x$ and avoids apparent singularities in the usual sine-ratio array formula.