I would change to an algebraic system by making new variables for the trigs and adding new equations to account for the trig identities. I've shown this in past but cannot seem to find links today so I'll just redo from scratch.
--- edit ---
Found a link:
http://community.wolfram.com/groups/-/m/t/1324461?p_p_auth=tZRhJF5X
--- end edit ---
Here is the setup, with code corrected as needed.
decz = 50*Degree;
raz = 50 Degree;
xyip[decx_, rax_, decy_, ray_] :=
Cos[decx]*Cos[rax]*Cos[decy]*Cos[ray] +
Cos[decx]*Sin[rax]*Cos[decy]*Sin[ray] + Sin[decx]*Sin[decy]
toZip[dec_, ra_] :=
Cos[dec]*Cos[ra]*Cos[decz]*Cos[raz] +
Cos[dec]*Sin[ra]*Cos[decz]*Sin[raz] + Sin[dec]*Sin[decz]
norm[dec_, ra_] := (Cos[dec]*Cos[ra])^2 + (Cos[dec]*Sin[ra])^2 + Sin[dec]^2
Here is the system with Cos
and Sin
suitably replaced.
trigpolys = {xyip[decx, rax, decy, ray], toZip[decx, rax],
toZip[decy, ray], norm[decx, rax] - 1, norm[decy, ray] - 1};
cosvars =
Union[Select[Cases[trigpolys, _Cos, Infinity], ! NumericQ[#] &]];
auxpolys = Map[#^2 + Sin[#[[1]]]^2 - 1 &, cosvars];
reprule = {Cos -> c, Sin -> s};
newpolys = Join[trigpolys, auxpolys] /. reprule
(* {c[decx] c[decy] c[rax] c[ray] + s[decx] s[decy] +
c[decx] c[decy] s[rax] s[ray],
c[40 \[Degree]] s[decx] + c[decx] c[rax] s[40 \[Degree]]^2 +
c[decx] c[40 \[Degree]] s[40 \[Degree]] s[rax],
c[40 \[Degree]] s[decy] + c[decy] c[ray] s[40 \[Degree]]^2 +
c[decy] c[40 \[Degree]] s[40 \[Degree]] s[ray], -1 +
c[decx]^2 c[rax]^2 + s[decx]^2 + c[decx]^2 s[rax]^2, -1 +
c[decy]^2 c[ray]^2 + s[decy]^2 + c[decy]^2 s[ray]^2, -1 +
c[decx]^2 + s[decx]^2, -1 + c[decy]^2 + s[decy]^2, -1 + c[rax]^2 +
s[rax]^2, -1 + c[ray]^2 + s[ray]^2} *)
Feed it to NSolve
. It turns out to be underdetermined (even though there were more equations than unknowns). All the same there will be usable solutions. We extract the ones that give rise to real angles.
sols = NSolve[newpolys];
realsols = Select[sols, FreeQ[Variables[newpolys] /. #, Complex] &];
goodsols =
Select[realsols,
AllTrue[Variables[newpolys] /. #, -1 <= # <= 1 &] &];
(* NSolve::infsolns: Infinite solution set has dimension at least 1. Returning intersection of solutions with (168848 c[decx])/118839-(135062 c[decy])/118839+(113347 c[rax])/118839-(58857 c[ray])/39613-(139306 s[decx])/118839-(133381 s[decy])/118839-(134851 s[rax])/118839+(35780 s[ray])/39613 == 1. *)
Now map the arc cosine function over the cosine variables to recover possible solutions. In some cases one might need to add multiples of Pi, in others there is no hope because the system can give parasite solutions due to multivalued-ness of the inverse mappings in use.
origvals =
Table[Map[ArcCos, cosvars /. reprule /. solj], {solj, goodsols}];
N[origvals]
(* Out[465]= {{0.592129, 0.32452, 1.51345, 2.85598}, {2.96426, 2.47526,
0.482896, 1.91231}, {0.19045, 2.48012, 2.67529, 1.88646}, {0.698131,
3.14081, 0.874258, 0.697196}, {0.189018, 0.662022, 0.468125,
1.88927}, {0.298065, 0.608675, 1.07298, 0.282036}, {2.54908,
2.81766, 2.90861, 2.03178}, {2.49168, 2.92305, 0.435177,
2.17559}, {0.69111, 0.0832024, 2.43752, 2.34391}, {2.59195, 2.75767,
0.120306, 1.94117}, {0.480135, 2.67816, 1.3677,
0.0601801}, {0.36032, 2.57417, 2.90912, 1.58083}} *)
Three of these give actual solutions to the original system.
Chop[
Map[trigpolys /. Thread[{decx, decy, rax, ray} -> #] &, N@origvals]]
(* Out[460]= {{0.355923, 0.855103, 0, 0, 0}, {0.218072, -0.45012,
0.217554, 0, 0}, {-0.429713, 0, 0.202429, 0, 0}, {-0.753563,
0.984807, -0.632316, 0, 0}, {0.231002, 0.724315, 0.737708, 0,
0}, {0.71937, 0.827124, 0.876019, 0, 0}, {0.680725, 0.666973, 0, 0,
0}, {0, 0, 0, 0, 0}, {0.817483, 0.491215, 0.127326, 0, 0}, {0, 0, 0,
0, 0}, {0, 0.855505, -0.0529827, 0, 0}, {0, 0, 0, 0, 0}} *)