We can ask sage for all subfields, take then those that are totally real among them.
This first try gives:

```
sage: var('x');
sage: L.<a> = NumberField(x^8-7)
sage: [F_data for F_data in L.subfields() if F_data[0].is_totally_real()]
[(Number Field in a0 with defining polynomial x,
Ring morphism:
From: Number Field in a0 with defining polynomial x
To: Number Field in a with defining polynomial x^8 - 7
Defn: 0 |--> 0,
None),
(Number Field in a1 with defining polynomial x^2 - 7,
Ring morphism:
From: Number Field in a1 with defining polynomial x^2 - 7
To: Number Field in a with defining polynomial x^8 - 7
Defn: a1 |--> -a^4,
None)]
```

So the totally real subfields in the list come with their embeddings. I.e. the `F_data`

is a tuple with three components, the pythonically 0.th component is the field.

The first field listed above is $\Bbb Q=\Bbb Q(0)$, where $0$ is the root of $x$, the second field is the field $\Bbb Q(\sqrt 7)$, where $\sqrt 7$ (considered algebraically) is the root of $x^2-7$. We want the maximal degree, so an idea is to sort w.r.t. the degrees and take the maximal degree. The code doing this is:

```
sage: var('x');
sage: L.<a> = NumberField(x^8-7)
sage: L_totally_real_subfields = [(F, emb) for (F, emb, _) in L.subfields() if F.is_totally_real()]
sage: L_totally_real_subfields.sort(key=lambda data: data[0].degree())
sage: L_totally_real_subfields[-1]
(Number Field in a1 with defining polynomial x^2 - 7,
Ring morphism:
From: Number Field in a1 with defining polynomial x^2 - 7
To: Number Field in a with defining polynomial x^8 - 7
Defn: a1 |--> -a^4)
```

Let's "do the same" with some nice cyclotomic field...

```
sage: K.<z> = CyclotomicField(17)
sage: K_totally_real_subfields = [F for (F, emb, _) in K.subfields() if F.is_totally_real()]
sage: len(K_totally_real_subfields)
4
sage: K_totally_real_subfields.sort(key=lambda F: F.degree())
sage: K_totally_real_subfields[-1]
Number Field in z3 with defining polynomial x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1 with z3 = 1.864944458808712?
sage: _.defining_polynomial().roots(ring=QQbar, multiplicities=False)
[-1.965946199367804?,
-1.700434271459229?,
-1.205269272758513?,
-0.5473259801441657?,
0.1845367189266040?,
0.8914767115530766?,
1.478017834441319?,
1.864944458808712?]
```