I'm guessing you wish to verify or at least check the results of that paper. This can be done somewhat more generally using methods based on Groebner bases. I'll also use a resultant or two, but later on.

First we can set up the equations, using angles th[1], etc.

Clear[m]
thetas = Array[th, 3];
trigpolys = TrigExpand[Join[{Total[Cos[thetas]] - m}, 
    Table[Total[Cos[j*thetas]], {j, {5, 7}}]]];
trigvars = Variables[trigpolys];

Next we "algebraicize", that is, turn these into polynomials in new variables c[...] instead of Cos[...]. We also add the defining relations involving these c variables and corresponding s variables, in order to account for trig identities.

reprule = {Cos[t_] :> c[t], Sin[t_] :> s[t]};
cosvars = Cases[trigvars, Cos[t_] :> c[t]];
auxpolys = Map[#^2 + s[#[[1]]]^2 - 1 &, cosvars];
newpolys = Join[trigpolys /. reprule, auxpolys];
newvars = Variables[newpolys];

Here is what we have.

In[127]:= newvars

(* Out[127]= {m, c[th[1]], c[th[2]], c[th[3]], s[th[1]], s[th[2]], s[th[3]]}*)

In[128]:= newpolys

(* Out[128]= {-m + c[th[1]] + c[th[2]] + c[th[3]], 
 c[th[1]]^5 + c[th[2]]^5 + c[th[3]]^5 - 10 c[th[1]]^3 s[th[1]]^2 + 
  5 c[th[1]] s[th[1]]^4 - 10 c[th[2]]^3 s[th[2]]^2 + 
  5 c[th[2]] s[th[2]]^4 - 10 c[th[3]]^3 s[th[3]]^2 + 
  5 c[th[3]] s[th[3]]^4, 
 c[th[1]]^7 + c[th[2]]^7 + c[th[3]]^7 - 21 c[th[1]]^5 s[th[1]]^2 + 
  35 c[th[1]]^3 s[th[1]]^4 - 7 c[th[1]] s[th[1]]^6 - 
  21 c[th[2]]^5 s[th[2]]^2 + 35 c[th[2]]^3 s[th[2]]^4 - 
  7 c[th[2]] s[th[2]]^6 - 21 c[th[3]]^5 s[th[3]]^2 + 
  35 c[th[3]]^3 s[th[3]]^4 - 7 c[th[3]] s[th[3]]^6, -1 + c[th[1]]^2 + 
  s[th[1]]^2, -1 + c[th[2]]^2 + s[th[2]]^2, -1 + c[th[3]]^2 + 
  s[th[3]]^2} *)

For any given numeric value of m we can solve for these trig variables, find real solutions, then take arccosines to recover values for the variables. The tool of choice for this example is NSolve, and I set it to use a method I trust to figure out which solutions are real valued.

solveForAngles[mval_] := Module[
  {sol, rsol, rsolUnique, angles, np2 = newpolys /. m -> mval},
  sol = NSolve[np2, Variables[np2], Method -> "CompanionMatrix"];
  rsol = Select[sol, FreeQ[#, Complex] &];
  If[Length[rsol] > 0,
   rsol = Union[cosvars /. rsol, SameTest -> (Norm[#1 - #2] <= .01 &)]];
  rsolUnique = Select[rsol, #[[1]] >= #[[2]] >= #[[3]] &];
  angles = Map[ArcCos, rsolUnique]/Degree
  ]

I show a few examples. The first is in agreement with the paper. The others are in ranges where it is claimed no solutions exist. Ours indeed violate a constraint. I'll return to that later.

AbsoluteTiming[solveForAngles[160/100]]

(* Out[119]= {17.063261, {{39.0176641435, 54.3352645924, 
   76.1130573337}, {19.0061442175, 52.4438553921, 87.4220926268}}} *)

These next three areclaimed not to have viable solutions. They have "solutions" but those all violate the requirement that all angles be no larger than Pi/2 (90 degrees, that is).

AbsoluteTiming[solveForAngles[1]]

(* Out[120]= {1.118488, {{43.2203549792, 71.9462855358, 92.2141079569}}} *)

AbsoluteTiming[solveForAngles[78/100]]

(* Out[121]= {16.871983, {{24.7173755287, 48.9384491915, 
   141.743882101}, {12.9921299006, 82.967292915, 108.471739602}}} *)

AbsoluteTiming[solveForAngles[70/100]]

(* Out[122]= {1.014744, {{29.7307422945, 52.8316106998, 
   140.581186894}, {11.8669526034, 84.7929687949, 111.677542596}}} *)

We can actually find values for m at which the number of real solutions might change. The idea is to determine for which values there might be a multiple root. It is at such points that real root counts can change. The method I show involves finding a polynomial in m that covers what is called the real discriminant variety. This is not really the best method to find what we want, but it works.

--- edit 1 ---

Below is not quite correct. Some of the points it turns up are fine but several needed ones are missing. I'll give a different method in edits below.

--- end edit 1 ---

Step 1: Find a univariate polynomial in a "random" linear combination of the c variables. This uses a GroebnerBasis computation with a term order that eliminates all variables except a new one that defines this combination, and m.

Step 2: Take the discriminant with respect to this new variable. It gives a polynomial in m for which that polynomial will have multiple roots in the new variable.

Step 3: Repeat steps 1 and 2 with a different random combination.

Step 4: Take the GCD of the two discriminants. This removes some (not necessarily all) spurious factors from the discriminants. Factor it to remove multiplicity.

Step 5: Find the roots of this GCD. Order them. In the interval between each pair the number of real roots will be constant.

Step 6: Take a test point in each interval and use the solver above to see how many roots we get.

Here are steps 1 and 2.

keepvars = {new, m}; 
AbsoluteTiming[
 poly1 = First[
    GroebnerBasis[Join[newpolys, {cosvars.{1, 2, -5} - new}], 
     keepvars, Complement[newvars, keepvars], 
     MonomialOrder -> EliminationOrder]];]
AbsoluteTiming[disc1 = Discriminant[poly1, new];]

(* Out[79]= {4.323621, Null}

Out[80]= {27.160578, Null} *)

Step 3:

AbsoluteTiming[
 poly2 = First[
    GroebnerBasis[Join[newpolys, {cosvars.{1, -4, 3} - new}], 
     keepvars, Complement[newvars, keepvars], 
     MonomialOrder -> EliminationOrder]];]
AbsoluteTiming[disc2 = Discriminant[poly2, new];]

(* Out[81]= {4.058456, Null}

Out[82]= {27.027516, Null} *)

Step 4:

disc12 = PolynomialGCD[disc1, disc2];
fdisc12 = FactorList[disc12]

(* Out[86]= {{\
2396842397475203197019749003558903268848355855677200505019017827417895\
6065439121965366330359490883444670117808992944128000000000000000000000\
00000000000000000000000000000000000000000000000, 1}, {m, 
  34}, {5 - 20 m^2 + 16 m^4, 
  15}, {-1715000 + 16078125 m^2 - 59682000 m^4 + 114072000 m^6 - 
   125283200 m^8 + 83417600 m^10 - 34263040 m^12 + 8458240 m^14 - 
   1146880 m^16 + 65536 m^18, 
  6}, {-3828844687500 + 50094051328125 m^2 - 284043552187500 m^4 + 
   923218114000000 m^6 - 1922050921000000 m^8 + 
   2715567495200000 m^10 - 2692924281600000 m^12 + 
   1913068505600000 m^14 - 984907212800000 m^16 + 
   368870268928000 m^18 - 100048214425600 m^20 + 
   19361352908800 m^22 - 2594701312000 m^24 + 227834593280 m^26 - 
   11744051200 m^28 + 268435456 m^30, 3}} *)

It has several factors and that makes one suspect spurious factors remain. A repeat of step 3 did not improve on this situation. We can live with it though since the testing will give root counts in each interval.

--- edit 2 ---

Instead of working with these discriminant polynomials, we can find the set of values for m for which the Jacobian matrix becomes singular.

jac = Det[Grad[newpolys, Complement[newvars, {m}]]];

AbsoluteTiming[
 mpoly = First[
    GroebnerBasis[Join[newpolys, {jac}], {m}, 
     Complement[newvars, {m}], MonomialOrder -> EliminationOrder]];]

(* Out[286]= {69.60784, Null} *)

mvals = N[Sort[m /. NSolve[mpoly, WorkingPrecision -> 100]]];
posmvals = Select[mvals, Im[#] == 0 && Re[#] > 0 &];
posmvals = Join[{0}, posmvals, {3/2*posmvals[[-1]]}];
intervals = Partition[posmvals, 2, 1];
checkpoints = Map[Mean, intervals]
intervals[[-1]] = {intervals[[-1, 1]], Infinity};
intervals

(* Out[347]= {0.255173537989, 0.513999093158, 0.663549719281, \
0.864459502625, 0.925903176502, 0.933391430476, 0.944687127846, \
0.981395292443, 1.10482777725, 1.34069141389, 1.5016490019, \
1.67923612424, 1.84469230783, 2.18921440071, 2.63988589, \
2.76237649831, 2.78108889041, 3.49173418289}

Out[349]= {{0, 0.510347075978}, {0.510347075978, 
  0.517651110338}, {0.517651110338, 0.809448328224}, {0.809448328224, 
  0.919470677026}, {0.919470677026, 0.932335675978}, {0.932335675978, 
  0.934447184975}, {0.934447184975, 0.954927070716}, {0.954927070716, 
  1.00786351417}, {1.00786351417, 1.20179204033}, {1.20179204033, 
  1.47959078745}, {1.47959078745, 1.52370721636}, {1.52370721636, 
  1.83476503212}, {1.83476503212, 1.85461958354}, {1.85461958354, 
  2.52380921788}, {2.52380921788, 2.75596256213}, {2.75596256213, 
  2.7687904345}, {2.7687904345, 
  2.79338734632}, {2.79338734632, \[Infinity]}} *)

--- end edit 2 ---

Steps 5 and 6:

facs = Rest[fdisc12[[All, 1]]]

(* Out[90]= {m, 
 5 - 20 m^2 + 16 m^4, -1715000 + 16078125 m^2 - 59682000 m^4 + 
  114072000 m^6 - 125283200 m^8 + 83417600 m^10 - 34263040 m^12 + 
  8458240 m^14 - 1146880 m^16 + 65536 m^18, -3828844687500 + 
  50094051328125 m^2 - 284043552187500 m^4 + 923218114000000 m^6 - 
  1922050921000000 m^8 + 2715567495200000 m^10 - 
  2692924281600000 m^12 + 1913068505600000 m^14 - 
  984907212800000 m^16 + 368870268928000 m^18 - 
  100048214425600 m^20 + 19361352908800 m^22 - 2594701312000 m^24 + 
  227834593280 m^26 - 11744051200 m^28 + 268435456 m^30} *)

mvals = Sort[Flatten[m /. Map[NSolve, facs]]];
posmvals = Select[mvals, Im[#] == 0 && Re[#] > 0 &];
posmvals = Join[{0}, posmvals, {3/2*posmvals[[-1]]}];
intervals = Partition[posmvals, 2, 1];
checkpoints = Map[Mean, intervals]
intervals[[-1]] = {intervals[[-1, 1]], Infinity};
intervals

(* Out[186]= {0.293892626146, 0.698616790258, 0.864459502647, \
0.935263596683, 0.952991793506, 0.981395292443, 1.26578536525, \
1.6792361242, 1.84469230781, 2.18921440071, 2.63988589001, \
2.77467495424, 3.49173418293}

Out[188]= {{0, 0.587785252292}, {0.587785252292, 
  0.809448328223}, {0.809448328223, 0.919470677071}, {0.919470677071, 
  0.951056516295}, {0.951056516295, 0.954927070716}, {0.954927070716, 
  1.00786351417}, {1.00786351417, 1.52370721632}, {1.52370721632, 
  1.83476503208}, {1.83476503208, 1.85461958354}, {1.85461958354, 
  2.52380921788}, {2.52380921788, 2.75596256213}, {2.75596256213, 
  2.79338734634}, {2.79338734634, \[Infinity]}} *)

AbsoluteTiming[testroots = Table[
   rsolns = solveForAngles[checkpoints[[j]]]; {j, intervals[[j]], 
    Length[rsolns], rsolns}, {j, Length[checkpoints]}]]

(* Out[166]= {61.667555, {{1, {0, 0.587785252292}, 
   3, {{54.7513536803, 81.090782053, 115.982882112}, {34.0546293484, 
     65.8926779525, 160.571705675}, {10.9553501839, 96.2722984079, 
     125.354186923}}}, {2, {0.587785252292, 0.809448328223}, 
   2, {{29.8345781569, 52.8968752584, 140.542849192}, {11.8494348378, 
     84.8279647455, 111.729204153}}}, {3, {0.809448328223, 
    0.919470677071}, 
   3, {{45.549565081, 78.6114802511, 91.9079242571}, {20.1718312387, 
     43.1705109004, 143.467951833}, {14.3841881375, 82.2977867001, 
     103.781313526}}}, {4, {0.919470677071, 0.951056516295}, 
   1, {{44.2920961884, 74.6562073327, 
     92.5869630927}}}, {5, {0.951056516295, 0.954927070716}, 
   1, {{43.9910189292, 73.8622231635, 
     92.5450530153}}}, {6, {0.954927070716, 1.00786351417}, 
   1, {{43.5204099866, 72.6775956991, 
     92.3774266156}}}, {7, {1.00786351417, 1.52370721632}, 
   2, {{39.9684306002, 63.06123127, 87.3434556197}, {20.358204688, 
     63.6962804253, 96.5967057007}}}, {8, {1.52370721632, 
    1.83476503208}, 
   2, {{37.7152160559, 53.888296551, 72.6136388049}, {17.1106438802, 
     49.0989543703, 86.0581993339}}}, {9, {1.83476503208, 
    1.85461958354}, 
   3, {{31.3553777718, 54.9045761749, 65.4305393327}, {8.7057447955, 
     22.7098471888, 93.7990961278}, {7.25406558369, 35.5229322073, 
     87.7756038844}}}, {10, {1.85461958354, 2.52380921788}, 
   1, {{15.1903498592, 39.6103995078, approx
     63.0150638453}}}, {11, {2.52380921788, 2.75596256213}, 
   0, {}}, {12, {2.75596256213, 2.79338734634}, 
   0, {}}, {13, {2.79338734634, \[Infinity]}, 0, {}}}} *)

Upshot: If I did this correctly, there are 3 solutions for m in the range from 0 tp0 approx 0.588, 2 from there to around .809, 3 again to .919, the 1 in the next few intervals, taking us to approx 1.008, 2 in the next two intervals to approx 1.835, 3 from there to 1.855, 1 from there to 2.524, and no solutions thereafter. The "isolated" solutions in the paper seem mostly to correspond to interval endpoints above. This can happen if a pair of complex roots cross the real axis and continue back into complex space.

I cannot say I carefully checked every equation. But several of the interval endpoints seem to agree with points noted in the paper. We simply disagree, in some intervals, on the root counts. It is possible that I misunderstand their criteria for viable roots.

--- edit 3 ---

Here is the corrected interval check.

AbsoluteTiming[
 testroots = 
  Table[Print[j]; 
   rsolns = solveForAngles[checkpoints[[j]]]; {Length[rsolns], 
    intervals[[j]], rsolns}, {j, Length[checkpoints]}]]

(* {98.445928, {
{3, {0, 
    0.510347075978}, {{55.4865182669, 82.3874510418, 
     116.352987938}, {35.7088346767, 67.132149226, 
     160.983891036}, {11.3356453419, 97.409690159, 
     126.609338463}}}, 
{2, {0.510347075978, 
    0.517651110338}, {{49.5663369914, 69.9342373027, 
     118.533114359}, {10.2989741549, 89.887555732, 
     118.154559552}}}, 
{2, {0.517651110338, 
    0.809448328224}, {{32.8512904477, 54.6075533725, 
     139.087037752}, {11.4304778133, 85.7378166282, 
     113.01279735}}}, 
{3, {0.809448328224, 
    0.919470677026}, {{45.5495650814, 78.6114802526, 
     91.9079242565}, {20.1718312399, 43.1705109026, 
     143.467951832}, {14.3841881371, 82.2977866997, 
     103.781313528}}},
 {2, {0.919470677026, 
    0.932335675978}, {{44.4532853355, 75.0970424867, 
     92.5849841136}, {13.0057575463, 29.0882553625, 
     157.267074562}}}, 
{1, {0.932335675978, 
    0.934447184975}, {{44.3242132711, 74.7430698199, 
     92.5880133494}}}, 
{1, {0.934447184975, 
    0.954927070716}, {{44.1313669796, 74.2280629526, 
     92.5715535077}}},
 {1, {0.954927070716, 
    1.00786351417}, {{43.5204099866, 72.6775956991, 
     92.3774266156}}}, 
{3, {1.00786351417, 
    1.20179204033}, {{41.6716287257, 68.2180003377, 
     90.7572336902}, {18.1975797709, 67.5868325642, 
     103.087587897}, {9.46369977719, 43.6717914389, 
     127.219444268}}}, 
{2, {1.20179204033, 
    1.47959078745}, {{39.5686202031, 60.7566121796, 
     85.3362100186}, {20.8724726924, 61.4040188189, 
     94.1469421284}}}, 
{3, {1.47959078745, 
    1.52370721636}, {{39.4242401699, 56.2098791719, 
     80.0358640193}, {20.4372411221, 56.065975505, 
     89.6358329844}, {6.43504243337, 18.9100737994, 
     115.981374633}}}, 
{2, {1.52370721636, 
    1.83476503212}, {{37.715215944, 53.8882965349, 
     72.6136388882}, {17.1106440659, 49.0989543658, 
     86.0581992805}}}, 
{3, {1.83476503212, 
    1.85461958354}, {{31.3553779352, 54.9045764169, 
     65.4305390203}, {8.70575731679, 22.7098478084, 
     93.7990939876}, {7.25408029759, 35.522931441, 
     87.7756024696}}}, 
{1, {1.85461958354, 
    2.52380921788}, {{15.1903493275, 39.6103988393, 
     63.01506448}}}, 
{0, {2.52380921788, 
    2.75596256213}, {}}, 
{1, {2.75596256213, 
    2.7687904345}, {{6.54308173016, 16.1746714205, 
     36.0529819078}}}, 
{0, {2.7687904345, 
    2.79338734632}, {}}, 
{0, {2.79338734632, \[Infinity]}, {}}}} *)

One observation is that where the paper claimed an isolated solution around 2.76 there is actually a (narrow) range containing solutions, specifically, the interval from approx 2.756 to 2.769. Similarly, at around 0.809 we get a third real root in the interval with that as lower bound, and there is a narrow region where that third root has all angles less than Pi/2. Similar no doubt happens at other points where th paper claimed isolated solutions.

I have not figured out where I messed up the prior discrimination variety computation. I may figure it out...eventually...

--- end edit 3 ---

--- edit 4 ---

One constraint in the paper is that angles not exceed 90 degrees. This translates to their cosines not becoming negative. It turns out to be quite simple to find values of m for which this happens. Since the setup is symmetric in all angles we only need find values of m for which one cosine changes signs. To do this we need to find a polynomial in m and one cosine variable. We then set that cosine to zero and solve for m. This gives the st of such values where an angle can cross 90 degrees.

keepvars = {cosvars[[1]], m}; AbsoluteTiming[
 poly1 = First[
    GroebnerBasis[newpolys, keepvars, Complement[newvars, keepvars], 
     MonomialOrder -> EliminationOrder]];]
Union[Select[m /. NSolve[poly1 /. cosvars[[1]] -> 0], 
  Im[#] == 0 && Re[#] >= 0 &]]

(* Out[658]= {0.0385876, Null}

Out[659]= {0., 0.510061, 0.825296, 0.919098, 1.1461, 1.48713, 1.85442} *)

So these points should also be used in constructing the intervals, at which point we can do the counts with the added constraint that all angles are acute (I don't have time to do this right now though). The proximity of some of these new values to the ones above should explain why some of the m-intervals of viable solutions are so narrow.

--- end edit 4 ---

Which brings us to...

--- edit 5 ---

I believe I now understand why the discriminant method failed to give important values in m, as well as why the Jacobian method did find them. The latter finds all points in parameter space where the system is singular. This can happen for reasons other than multiplicity. In particular it can happen when a cosine variable is +-1 and some of the equations in our system are stationary.

As it happens, these (well, =1 at least) are important values for this problem because they show where an angle may go between positive and negative, or where the root to the algebraic system is real but its arccos (giving the angle) is imaginary because the value exceeded 1 in magnitude.

These are not points of multiplicity however, and that's why they do not appear in the discriminant variety. The way to recover them using that method is similar to what we did to determine when angles transition between acute and obtuse. Specifically, we separately solve for m under the assumption that a cosine variable is 1 or -1 and join those points to the discriminant points.

keepvars = {cosvars[[1]], m};
    GroebnerBasis[newpolys, keepvars, Complement[newvars, keepvars], 
     MonomialOrder -> EliminationOrder]];]
Union[Select[m /. NSolve[poly1 /. cosvars[[1]] -> 1], 
  Im[#] == 0 && Re[#] >= 0 &]

(* Out[671]= {0., 0.517651, 0.932336, 1.20179, 1.47959, 2.76879} *)

keepvars = {cosvars[[1]], m};
 poly1 = First[
    GroebnerBasis[newpolys, keepvars, Complement[newvars, keepvars], 
     MonomialOrder -> EliminationOrder]];]
Union[Select[m /. NSolve[poly1 /. cosvars[[1]] -> -1], 
  Im[#] == 0 && Re[#] >= 0 &]

(* Out[684]= {0., 0.510347, 0.934447} *)

These are indeed the values that were missing from the discriminant method but present in the Jacobian computation.

Another thing to note is that, as in the original article, one can recast this purely in terms of cosines. The outcome will be no different but the computation speed might improve considerably from having fewer equations and variables.

--- end edit 5 ---