I believe these equations first appear in Ian Barnes' 1990 PhD thesis Microstructure of bicontinuous phases in surfactant systems.
In section 6.2.1 he derives the Fourier representation of periodic equipotential surfaces:
$$
\varphi(x, y, z) = \frac{8}{\pi}\sum_{j=1}^N q_j \sum_{h,k,l=0}^\infty{\vphantom{\sum}}{\!\!\!'\,\,\,} \frac{\cos(2\pi h(x-x_j))\cos(2\pi k(y-y_j))\cos(2\pi l(z-z_j))}{h^2+k^2+l^2},
$$
where the prime on the last sum indicates that for each index, terms are to be multiplied by $1/2$ if it is zero, and the $ h = k = l = 0$ term is taken to be $0$.
Here there are $N$ charges in $[0, 1]^3$ with values $q_j$ at points $(x_j, y_j, z_j)$, and these charges are then periodically tiled.
Explicit examples are covered in section 6.3 by truncating the above infinite sum, as first done in a paper by Mackay. For the gyroid, in section 6.3.5, there are 16 charges:
8 points with charge $+1$ at $\small \left(\frac{1}{8},\frac{1}{8},\frac{1}{8}\right), \left(\frac{5}{8},\frac{3}{8},\frac{7}{8}\right), \left(\frac{7}{8},\frac{5}{8},\frac{3}{8}\right), \left(\frac{3}{8},\frac{7}{8},\frac{5}{8}\right), \left(\frac{5}{8},\frac{5}{8},\frac{5}{8}\right), \left(\frac{1}{8},\frac{7}{8},\frac{3}{8}\right), \left(\frac{3}{8},\frac{1}{8},\frac{7}{8}\right), \left(\frac{7}{8},\frac{3}{8},\frac{1}{8}\right)$
8 points with charge $-1$ at $\small \left(\frac{7}{8},\frac{7}{8},\frac{7}{8}\right), \left(\frac{3}{8},\frac{5}{8},\frac{1}{8}\right), \left(\frac{1}{8},\frac{3}{8},\frac{5}{8}\right), \left(\frac{5}{8},\frac{1}{8},\frac{3}{8}\right), \left(\frac{3}{8},\frac{3}{8},\frac{3}{8}\right), \left(\frac{7}{8},\frac{1}{8},\frac{5}{8}\right), \left(\frac{5}{8},\frac{7}{8},\frac{1}{8}\right), \left(\frac{1}{8},\frac{5}{8},\frac{7}{8}\right)$
We can verify the well known equation that approximates the gyroid by truncating at $h,k,l\leq1$ in Mathematica:
Do[q[j] = 1, {j, 8}]
Do[q[j + 8] = -1, {j, 8}]
pts = {
{1, 1, 1}, {5, 3, 7}, {7, 5, 3}, {3, 7, 5},
{5, 5, 5}, {1, 7, 3}, {3, 1, 7}, {7, 3, 1},
{7, 7, 7}, {3, 5, 1}, {1, 3, 5}, {5, 1, 3},
{3, 3, 3}, {7, 1, 5}, {5, 7, 1}, {1, 5, 7}
}/8;
Do[{xx[j], yy[j], zz[j]} = pts[[j]], {j, 16}]
gyrapprox[n_] := FullSimplify[
8/π Sum[q[j] (Sum[Piecewise[{
{((1/2)^(Boole[h == 0] + Boole[k == 0] + Boole[l == 0])*
Cos[2π h (x - xx[j])] Cos[2π k (y - yy[j])] Cos[2π l (z - zz[j])])/(h^2 + k^2 + l^2),
h != 0 || k != 0 || l != 0}}], {h, 0, n}, {k, 0, n}, {l, 0, n}]), {j, 16}],
ComplexityFunction -> (LeafCount[#] + 1000 Total[Cases[#, p_Plus :> Length[p], ∞]] &)]
The approximate formula (after dividing both sides by a constant):
g1 = Simplify[π/16 gyrapprox[1] /. {v : x | y | z -> v/(2π)}] == 0
Cos[y]Sin[x] + Cos[z]Sin[y] + Cos[x]Sin[z] == 0
We can take more terms now:
g2 = Simplify[π/16 gyrapprox[2] /. {v : x | y | z -> v/(2π)}] == 0
Cos[y]Sin[x] + Cos[z]Sin[y] + Cos[x]Sin[z] + 2Sin[2x]Sin[2y]Sin[2z]/3 == 0
Omitting the output:
g3 = Simplify[π/16 gyrapprox[3] /. {v : x | y | z -> v/(2π)}] == 0;
Unfortunately only g1, the well known formula, gives an accurate approximation to the gyroid:
opts = {PlotTheme -> "Minimal", Axes -> False, Boxed -> False};
ContourPlot3D[#, {x, 0, 2π}, {y, 0, 2π}, {z, 0, 2π}, opts] & /@ ({g1, g2, g3} /. y -> -y)
