Integrating Superposition of several sinusoidal waves squared times its derivative

I'm currently finishing some research from a internship I did last summer and part of the project involves calculating the centroids of coronal holes from synoptic maps. If you aren't aware of what these are, they are basically holes or regions of the corona, an upper layer of the suns atmosphere, where it is less dense. The purpose of this is so I can, track how the motion of these coronal holes change over time. The way this is done is by taking a certain latitude band, say 10 degrees south to 10 degrees north, finding the centroids for coronal holes after each carrington rotation, solar month, and finding the slope with all of these points. In general, calculating these centroids is not an easy task because in general the shapes of these regions, on solar maps, are non uniform and can take on any flat 2D shape. However last summer I did find a solution that should work, at least theoretically. Essentially my idea was to generate parametric plots that look identical to the images of each coronal hole and then use green's theorem, with the parametric equations, to calculate the centroids through line integration. Calculus 3 was still fresh in my mind at that point and I thought, if we can use Green's Theorem to calculate area, why not centroids as well? After doing some research I found that these line integrals do in fact exist and I went through the derivations to prove to myself that the line integrals were equivalent to the double integrals used to calculate centroids. If you are wondering what these line integrals are, here they are.

And here are the double integrals that the line integrals are equivalent to.

My mentors had another method that they had planned for me to use on IDL but encouraged me to try applying my method instead. I ended up working on Mathematica instead, because of useful code I found on other forums. From there I added some of my own code and did a bunch of tinkering over the course of several weeks. At the time, I did not know if this was the most efficient method computationally and I still don't but I would still like to get this to be fully functional. I got the method to work for the most part, but at times I'm still having issues getting an accurate answer. Sometimes it gives an answer that looks reasonable but at other times it gives an answer that is obviously completely off just by eyeballing it. I believe that the issue solely lies in convergence issues with the integrals but I will show all of my code, just in case some of the issues lie elsewhere.

param[x_, m_, t_] :=
Module[{f, n = Length[x], nf},
f = Chop[Fourier[x]][[;; Ceiling[Length[x]/2]]];
nf = Length[f];
Total[Rationalize[
2 Abs[f]/Sqrt[n] Sin[
Pi/2 - Arg[f] + 2. Pi Range[0, nf - 1] t], .01][[;;
Min[m, nf]]]]]

tocurve[Line[data_], m_, t_] := param[#, m, t] & /@ Transpose[data]
img = \!$$\* GraphicsBox[ TagBox[RasterBox[CompressedData[" 1:eJzt229sE3Ucx/HP/f/dXa/XrmUmiFl8hCgzMWbBTDfnnH+GmyYVE+IziJMY khG3ZUYhZiISM0WmENQwUNRoGNncEP8skoFREQMD94dtOjIyKZtUxv5gO2XZ 6ukjn5ioafzu4PtK2vTaJ+82l+v3d71ev3pd7HEZQLXw7mKrnrqrqmrV0w+H vI1HKqufWFNZ8VhpZU3Fmoqq21Yr3pMLvds6CfjjcZoxxhhjjDHGGGOMZU7s 7r22uhvYodhNktMCvC9rb1BHMZ95fdu3QIy6ghHLW9oHHIuGxywxCHwXDl2S EAe+iYpG6jTmG0CdJn9RU9tFHZJh5eWtwHvAx0Bjff2V9u4y697ir71DRyR8 2XWms7LStjkDpBw97Wjp667pp65jvrFs2X6B4fz8buqQDAPeCVs/Rs1kUBsU yn7qnHntjqI24CjQD3QBA8CIrkxYxoQuex/j996ShzqQ+QbQrWD4xps/pw7J mA11nd5Mku2OuWIibMRltFIXzXfPbj+7qaF/8yvx+q2Jl172FjhdhjYOTKr2 jOWe9Y4nu976gbqRzXfLlw/YYkRYvwAnnnvhDHVOxtx6i7fA+Sxk9oQCHdQt vuQEhlR9XDaSemhWUfoltfW1LbxaZH+r9J7jQKfrJAN2WlHGluZ+SV2UYSWF zcBK6gpfAj6QtQSQgpGG9psuj0jYRx3F5ilghYrjAZECzgtz3HCTwMmcRV33 l3QW33mw6PZ26kBGTNL6nKw5SYlbkcvAsdwlV85CmGUc0GaK85HorOmmJHVY 038ylOSfZ+EOALtzb/h08WK+2ODqJUmfCOecKk4BTQ/EDlHnsPlOlt82xQVZ vhhwx5yslKqe1kV79qIPc3JaVGwPiUbXraFuZAQCzpshpztg93vfLPmF/IMO +0eAoVBkzs2aAno2bkr89aWCvJ1UVYxQMNhqqwlb8naJ9kdjR6hzmD/UbRyV 5IShz2rWpCQdamgYpC5ixCAaVNEbDU5HrZGVZT3UOcw3gL1BZ9wUSVkZta2p J9f+TF3EyNiRFm8UkcWZcE4KaJWkHdRFzDc2VJ821N6g8GaSSU2ZjthzQNuW zX3UXez/lpfXBOwzjSHLuRRwR4Fda2tPUEcxPwGaLT1uKNOqkpakX01lqKSw kzqKESgu+8owTmrKcLY1xtfPs//gmfUXJPmiok6ZYiYUnPb2ond38lhyNSoo /QjoEPKAqw9u28p/7mP/WvmDRxRlUNVGDRFfsOAcsIe6iNEoKtujKs3Xhvts 6TB1C/OfJTc9DzTZdqdu9GrGKdM8KsuvUkcxGgX3vQjURmxvPmmgbmE+E134 EBATdrOwDgMtmt7hugeB9dRdjDHGGGOMMcYYY4xl2O/DK24g "], {{0, 26}, {367, 0}}, {0, 255}, ColorFunction->RGBColor], BoxFormImageTag["Byte", ColorSpace -> "RGB", Interleaving -> True], Selectable->False], DefaultBaseStyle->"ImageGraphics", ImageSizeRaw->{367, 26}, PlotRange->{{0, 367}, {0, 26}}]$$
img = Binarize[img~ColorConvert~"Grayscale"];
lines = Cases[
Normal@ListContourPlot[Reverse@ImageData[img],
Contours -> {0.5}], _Line, -1];

ParametricPlot[
Evaluate[tocurve[#, 100000000, t] & /@ lines], {t, 0, 1},
Frame -> True, Axes -> False, ImageSize -> Large]
ParametricEquations = Evaluate[tocurve[#, 100000000, t] & /@ lines]
n = Length[ParametricEquations]
img2 = Binarize[img~ColorConvert~"Grayscale"];

bsFs = Cases[
Normal@ListContourPlot[Reverse@ImageData[img2], Contours -> {0.5}],
Line[x_] :> BSplineFunction[x], All];

ParametricPlot[Evaluate[Through@bsFs@t], {t, 0, 1}, Frame -> True,
Axes -> False, ImageSize -> Large, Axes -> false, Frame -> True]
stepsize = 0.0001
list = Table[
N[Evaluate[Through@bsFs@t] /. t -> i], {t, 0, 1, stepsize}]
xmatrix =
Table[Table[list[[i, j, 1]], {i, 1, (1/stepsize) + 1}], {j, 1, n}]
ymatrix =
Table[Table[list[[i, j, 2]], {i, 1, (1/stepsize) + 1}], {j, 1, n}]
Minx = Table[xmatrix[[i, Ordering[xmatrix[[i]], 1]]], {i, 1, n}]
Miny = Table[ymatrix[[i, Ordering[ymatrix[[i]], 1]]], {i, 1, n}]
x = Table[ParametricEquations[[i, 1]], {i, 1, n}]
y = Table[ParametricEquations[[i, 2]], {i, 1, n}]
Mx = Table[NMinValue[{x[[i]], 0 <= t <= 1}, t], {i, 1, n}]
My = Table[NMinValue[{y[[i]], 0 <= t <= 1}, t], {i, 1, n}]
NewParametricEquations =
Table[{x[[i]] - (Mx[[i]] - Minx[[i]]),
y[[i]] - (My[[i]] - Miny[[i]])}, {i, 1, n}]
xCentroids =
Table[NIntegrate[
NewParametricEquations[[i, 1]]^2*
D[NewParametricEquations[[i, 2]], t], {t, 0, 1}]/
NIntegrate[
2*NewParametricEquations[[i, 1]]*
D[NewParametricEquations[[i, 2]], t], {t, 0, 1}], {i, 1, n}]
yCentroids =
Table[NIntegrate[-NewParametricEquations[[i, 2]]^2*
D[NewParametricEquations[[i, 1]], t], {t, 0, 1}]/
NIntegrate[-2*NewParametricEquations[[i, 2]]*
D[NewParametricEquations[[i, 1]], t], {t, 0, 1}], {i, 1, n}]
Centroids = Table[{xCentroids[[i]], yCentroids[[i]]}, {i, 1, n}]
plot = ParametricPlot[NewParametricEquations, {t, 0, 1},
ImageSize -> Large, Axes -> False, Frame -> True, PlotRange -> All];
centroidpoints =
ListPlot[Table[{First[xCentroids[[i]]], First[yCentroids[[i]]]}, {i,
1, n}] -> Range[n], LabelingFunction -> Bottom,
PlotStyle -> {Red, PointSize[Medium]}];
coronalholepoints =
ListPlot[Table[{First[N[NewParametricEquations[[i, 1]] /. t -> 0]],
First[N[NewParametricEquations[[i, 2]] /. t -> 0]]}, {i, 1,
n}] -> Range[n], LabelingFunction -> Left,
PlotStyle -> {Blue, PointSize[Medium]}];
Show[plot, centroidpoints, coronalholepoints]


Here is an example of an image with some coronal holes I attempted to calculate the centroids of. All other colors have been filtered out in order to generate accurate parametric equations for each coronal hole (the blue regions in this case, which indicates that the coronal holes have a positive polarity).

And here is an image of the parametric plot of all the coronal holes, along with red dots which represent the calculated centroids. As you can see, the first two centroids look reasonable but the last one is way off and isn't even inside the coronal hole. As a reference, here is a parametric equations that Mathematica generated for the coronal holes.

{{{24.7056 - 1/45 Sin[1/11 - (3349 t)/13] -
1/77 Sin[3/8 - (754 t)/3] - 1/23 Sin[8/9 - (3217 t)/16] -
1/76 Sin[1/18 - (377 t)/2] - 1/28 Sin[1/4 - (1640 t)/9] -
1/30 Sin[7/13 - (1797 t)/11] - 1/53 Sin[2/9 - (1885 t)/12] -
1/23 Sin[22/15 - (4191 t)/29] - 1/28 Sin[5/9 - (2111 t)/16] -
1/22 Sin[7/5 - (2199 t)/25] - 1/16 Sin[5/4 - (1552 t)/19] -
1/9 Sin[1/4 - (553 t)/8] - 1/11 Sin[8/9 - (377 t)/6] -
1/8 Sin[13/10 - (553 t)/11] - 59/15 Sin[13/12 - (88 t)/7] +
1/5 Sin[(1891 t)/43] + 19/5 Sin[17/5 + (44 t)/7] +
5/4 Sin[103/26 + (132 t)/7] + 2/3 Sin[25/14 + (201 t)/8] +
6/13 Sin[47/14 + (377 t)/12] + 1/7 Sin[13/4 + (377 t)/10] +
1/9 Sin[25/6 + (509 t)/9] + 1/7 Sin[5/14 + (377 t)/5] +
1/24 Sin[23/14 + (377 t)/4] + 1/24 Sin[45/13 + (1307 t)/13] +
1/16 Sin[8/5 + (1175 t)/11] + 1/11 Sin[94/21 + (1131 t)/10] +
1/71 Sin[31/7 + (955 t)/8] + 1/25 Sin[13/5 + (377 t)/3] +
1/27 Sin[3/14 + (1244 t)/9] + 1/19 Sin[5/4 + (754 t)/5] +
1/33 Sin[65/22 + (2287 t)/13] + 1/23 Sin[58/19 + (1753 t)/9] +
1/27 Sin[3/2 + (3525 t)/17] + 1/55 Sin[20/7 + (1709 t)/8] +
1/65 Sin[14/3 + (2419 t)/11] + 1/28 Sin[10/19 + (1131 t)/5] +
1/85 Sin[49/11 + (3952 t)/17] + 1/43 Sin[10/3 + (3104 t)/13] +
1/33 Sin[5/11 + (5391 t)/22] + 1/56 Sin[7/2 + (2375 t)/9] +
1/43 Sin[14/15 + (3041 t)/11] + 1/92 Sin[41/11 + (1131 t)/4] +
1/80 Sin[32/9 + (2953 t)/10]}, {14.8678 -
1/86 Sin[5/4 - (1885 t)/6] - 1/83 Sin[16/15 - (3349 t)/13] -
1/21 Sin[6/5 - (3217 t)/16] - 1/21 Sin[3/5 - (1640 t)/9] -
1/46 Sin[7/11 - (1797 t)/11] - 1/36 Sin[7/8 - (1244 t)/9] -
1/11 Sin[4/7 - (2111 t)/16] - 1/44 Sin[17/18 - (955 t)/8] -
1/34 Sin[7/6 - (1307 t)/13] - 1/15 Sin[1 - (2199 t)/25] -
1/7 Sin[23/15 - (553 t)/11] - 3/13 Sin[23/15 - (377 t)/10] -
2/5 Sin[7/10 - (377 t)/12] - 8/11 Sin[4/3 - (132 t)/7] -
10/7 Sin[5/6 - (88 t)/7] - 145/18 Sin[13/9 - (44 t)/7] +
4/11 Sin[12/11 + (201 t)/8] + 2/9 Sin[1/3 + (1627 t)/37] +
1/8 Sin[14/5 + (509 t)/9] + 1/24 Sin[137/46 + (377 t)/6] +
1/10 Sin[35/12 + (553 t)/8] + 1/66 Sin[5/8 + (377 t)/5] +
1/9 Sin[57/14 + (1307 t)/16] + 1/30 Sin[5/6 + (1175 t)/11] +
1/22 Sin[50/11 + (1131 t)/10] + 1/21 Sin[8/3 + (377 t)/3] +
1/30 Sin[22/5 + (3324 t)/23] + 1/22 Sin[26/11 + (754 t)/5] +
1/27 Sin[34/15 + (1885 t)/12] + 1/24 Sin[19/7 + (1866 t)/11] +
1/69 Sin[23/8 + (2287 t)/13] + 1/23 Sin[1/5 + (377 t)/2] +
1/30 Sin[11/3 + (1753 t)/9] + 1/21 Sin[37/11 + (1709 t)/8] +
1/64 Sin[29/7 + (2419 t)/11] + 1/38 Sin[10/19 + (1131 t)/5] +
1/31 Sin[4/5 + (5391 t)/22] + 1/36 Sin[29/8 + (2375 t)/9] +
1/86 Sin[5/7 + (3041 t)/11] +
1/64 Sin[48/11 + (1131 t)/4]}}, {{217.089 -
1/41 Sin[36/35 - (2375 t)/9] - 1/65 Sin[2/3 - (3952 t)/17] -
1/17 Sin[1/30 - (3217 t)/16] - 1/18 Sin[31/21 - (2287 t)/13] -
1/19 Sin[1 - (1866 t)/11] - 1/22 Sin[3/13 - (1885 t)/12] -
1/12 Sin[3/7 - (4191 t)/29] - 1/25 Sin[11/12 - (377 t)/3] -
1/7 Sin[11/9 - (1552 t)/19] - 1/4 Sin[5/9 - (377 t)/5] -
4/15 Sin[1/9 - (509 t)/9] - 2/11 Sin[3/10 - (377 t)/12] -
9/7 Sin[7/6 - (201 t)/8] - 24/25 Sin[1/8 - (88 t)/7] +
82/15 Sin[51/13 + (44 t)/7] + 7/9 Sin[4 + (132 t)/7] +
3/11 Sin[34/9 + (377 t)/10] + 1/9 Sin[75/19 + (1627 t)/37] +
2/7 Sin[2/9 + (553 t)/11] + 1/8 Sin[29/7 + (377 t)/6] +
3/11 Sin[25/7 + (553 t)/8] + 1/14 Sin[13/3 + (2023 t)/23] +
1/11 Sin[27/11 + (377 t)/4] + 1/33 Sin[9/7 + (1307 t)/13] +
1/10 Sin[55/12 + (1175 t)/11] + 1/42 Sin[47/11 + (1131 t)/10] +
1/28 Sin[71/18 + (955 t)/8] + 1/13 Sin[11/5 + (2111 t)/16] +
1/22 Sin[40/11 + (1244 t)/9] + 1/17 Sin[35/18 + (754 t)/5] +
1/18 Sin[49/16 + (1797 t)/11] + 1/13 Sin[21/8 + (1640 t)/9] +
1/17 Sin[7/4 + (377 t)/2] + 1/37 Sin[5/8 + (1753 t)/9] +
1/50 Sin[42/11 + (1709 t)/8] + 1/28 Sin[23/5 + (2419 t)/11] +
1/60 Sin[5/4 + (3104 t)/13] + 1/70 Sin[49/15 + (5391 t)/22] +
1/36 Sin[25/6 + (754 t)/3] + 1/24 Sin[8/5 + (3349 t)/13] +
1/99 Sin[7/3 + (2972 t)/11] + 1/63 Sin[159/79 + (1885 t)/6] +
1/64 Sin[16/15 + (2375 t)/7] +
1/75 Sin[25/11 + (2463 t)/7]}, {13.1114 -
1/97 Sin[10/9 - (3041 t)/11] - 1/46 Sin[9/8 - (2375 t)/9] -
1/28 Sin[1/20 - (3952 t)/17] - 1/27 Sin[2/5 - (3217 t)/16] -
1/51 Sin[8/7 - (1753 t)/9] - 1/21 Sin[11/9 - (2287 t)/13] -
1/19 Sin[3/7 - (1866 t)/11] - 1/21 Sin[3/5 - (1885 t)/12] -
1/15 Sin[7/10 - (4191 t)/29] - 1/16 Sin[4/5 - (377 t)/3] -
1/12 Sin[16/13 - (1175 t)/11] - 1/34 Sin[4/13 - (1307 t)/13] -
1/11 Sin[4/13 - (377 t)/6] - 1/6 Sin[1/11 - (509 t)/9] -
11/23 Sin[5/16 - (201 t)/8] - 57/8 Sin[22/23 - (44 t)/7] +
27/13 Sin[47/10 + (88 t)/7] + 3/4 Sin[2/5 + (132 t)/7] +
6/11 Sin[25/7 + (377 t)/12] + 1/76 Sin[14/5 + (377 t)/10] +
1/8 Sin[70/23 + (1627 t)/37] + 2/9 Sin[13/3 + (553 t)/11] +
1/6 Sin[25/7 + (553 t)/8] + 1/8 Sin[1/4 + (377 t)/5] +
1/25 Sin[49/12 + (1307 t)/16] + 1/20 Sin[68/15 + (2023 t)/23] +
1/11 Sin[8/3 + (377 t)/4] + 1/46 Sin[10/7 + (1131 t)/10] +
1/48 Sin[41/9 + (955 t)/8] + 1/15 Sin[23/11 + (2111 t)/16] +
1/20 Sin[11/3 + (1244 t)/9] + 1/15 Sin[37/16 + (754 t)/5] +
1/26 Sin[19/6 + (1797 t)/11] + 1/16 Sin[34/13 + (1640 t)/9] +
1/29 Sin[9/5 + (377 t)/2] + 1/43 Sin[17/4 + (1709 t)/8] +
1/35 Sin[19/5 + (2419 t)/11] + 1/36 Sin[9/5 + (1131 t)/5] +
1/65 Sin[7/6 + (3104 t)/13] + 1/43 Sin[27/7 + (5391 t)/22] +
1/27 Sin[37/8 + (754 t)/3] + 1/32 Sin[18/11 + (3349 t)/13] +
1/68 Sin[31/12 + (2972 t)/11] + 1/79 Sin[1/10 + (1131 t)/4] +
1/72 Sin[3/8 + (3594 t)/11]}}, {{287.833 -
1/39 Sin[8/15 - (1552 t)/19] - 1/36 Sin[9/11 - (377 t)/5] -
1/22 Sin[11/8 - (553 t)/8] - 1/32 Sin[5/9 - (377 t)/6] -
1/5 Sin[7/11 - (1891 t)/43] - 1/16 Sin[5/9 - (377 t)/12] -
1/10 Sin[11/10 - (201 t)/8] + 19/7 Sin[19/6 + (44 t)/7] +
2/5 Sin[1/21 + (88 t)/7] + 1/8 Sin[2/5 + (132 t)/7] +
1/9 Sin[1/72 + (377 t)/10] + 1/37 Sin[27/7 + (553 t)/11] +
1/9 Sin[24/11 + (509 t)/9] + 1/31 Sin[16/7 + (2023 t)/23] +
1/60 Sin[11/5 + (377 t)/4] +
1/43 Sin[9/5 + (1175 t)/11]}, {19.1795 - 1/39 Cos[(1307 t)/16] -
1/21 Sin[9/13 - (553 t)/8] - 1/37 Sin[1/11 - (377 t)/6] -
1/11 Sin[14/9 - (1891 t)/43] - 1/21 Sin[8/15 - (377 t)/10] +
23/9 Sin[51/11 + (44 t)/7] + 1/3 Sin[1/5 + (88 t)/7] +
2/7 Sin[2/5 + (132 t)/7] + 1/12 Sin[79/17 + (201 t)/8] +
1/7 Sin[3/4 + (377 t)/12] + 1/15 Sin[9/2 + (553 t)/11] +
1/14 Sin[31/9 + (509 t)/9] + 1/27 Sin[55/12 + (377 t)/5] +
1/39 Sin[13/5 + (2023 t)/23] + 1/47 Sin[5/2 + (1175 t)/11] +
1/72 Sin[51/11 + (1131 t)/10]}}}


It goes without saying but these equations are very long. And on top of that for the integrand you have to square these equations and multiply by the derivative. That's a pretty complicated integral, which on top of the the large number of terms, I imagine, is also highly oscillating. Since the integrals are so complicated, I imagine that's what is causing the centroids to be inaccurate at times. Whenever I run my code, I've also been getting an error message, like this, which seems to support this.

NIntegrate::levmaxord: (287.833 -1/39 Sin[8/15-(1552 t)/19]-1/36 Sin[9/11-(377 t)/5]-1/22 Sin[11/8-(553 t)/8]-1/32 Sin[5/9-(377 t)/6]-1/5 Sin[7/11-(1891 t)/43]-1/16 Sin[5/9-(377 t)/12]-1/10 Sin[11/10-(201 t)/8]+19/7 Sin[19/6+(44 t)/7]+2/5 Sin[1/21+(88 t)/7]+1/8 Sin[2/5+(132 t)/7]+1/9 Sin[1/72+(377 t)/10]+1/37 Sin[27/7+(553 t)/11]+1/9 Sin[24/11+(509 t)/9]+1/31 Sin[16/7+(2023 t)/23]+1/60 Sin[11/5+(377 t)/4]+1/43 Sin[9/5+(1175 t)/11])^2 is a Levin function of differential order 1088 which exceeds value of option "MaxOrder" -> 50. Treating (287.833 -1/39 Sin[8/15-(1552 t)/19]-1/36 Sin[9/11-(377 t)/5]-1/22 Sin[11/8-(553 t)/8]-1/32 Sin[5/9-(377 t)/6]-1/5 Sin[7/11-(1891 t)/43]-1/16 Sin[5/9-(377 t)/12]-1/10 Sin[11/10-(201 t)/8]+19/7 Sin[19/6+(44 t)/7]+2/5 Sin[1/21+(88 t)/7]+1/8 Sin[2/5+(132 t)/7]+1/9 Sin[1/72+(377 t)/10]+1/37 Sin[27/7+(553 t)/11]+1/9 Sin[24/11+(509 t)/9]+1/31 Sin[16/7+(2023 t)/23]+1/60 Sin[11/5+(377 t)/4]+1/43 Sin[9/5+(1175 t)/11])^2 as a non-Levin function.


I've tried a bunch of different integration methods, but nothing seems to work better than the method Mathematica automatically uses. I wanted to have an alternative way to calculate the centroid, so I could compare my results and see how far off some of the more reasonable answers are. So I came up with an alternative method that turns the parametric plot into a region using DiscretizeGraphics and then calculates the centroid using RegionCentroid and I believe that this gives more accurate results. All three x coordinates for each method are different by at least a couple of pixels. This may not seem significant, but based on the size of the image 1 pixel in the x-direction is approximately 1 degree in longitude, which could make a significant difference when I'm calculating the slope. Do you know any way for these integrals to converge more accurately or perhaps some other unforeseen problems that are causing the numbers to be off? Any help would be greatly appreciated. Sorry in advance if my question is obvious. I would still consider myself a novice when it comes to Mathematica, so there's a lot I don't know still. Also, sorry that my post was so long, I may have wrote more than I needed but I wanted to give enough background to make clear the full context of the problem.

• My first guess (I don’t have a computer at hand) is that for the third figure, $x^2 \approx 84000$ and rationalizing the coefficient 2. Pi n to with 0.01 introduces small errors in the integral that when multiplied by $x^2$ are significant. Why not use the exact 2 Pi n for the coefficient? – Michael E2 Jun 14 '20 at 14:28
• The method here could probably be adapted to do your integrals exactly and efficiently. You’d need to add rules for the integrals of cubic products of sine and cosine. – Michael E2 Jun 14 '20 at 14:34
• Actually I don’t see the need to Rationalize[] at all. It’s cute for the “person curves.” But it adds an error, which again is fine for person curves. – Michael E2 Jun 14 '20 at 14:40

img = Import["https://i.stack.imgur.com/VrChT.jpg"];
`