How to find a basis of the coxeter plane I want to draw a picture that projects the E8 root system to its coxeter plane. 
The coxeter plane is defined as follows: the coxeter element of the weyl group of E8 has a simple eigenvalue $e^{2\pi i/h}$ where $h=30$. The real and imag part of the eigenvector of this eigenvalue span a 2d plane $P$ in $\mathbb{R}^8$. To project the root system to $P$, just find a pair of orthogonal unit vectors $u,v$ in $P$, and for each root $x$, compute the inner product $(u,x)$ and $(v,x)$, this is the 2d point we want.
My problem is: how to find such a basis? Of course one can write down explicitly the matrices of the simple roots, multiply them in any order, this is the matrix of the coxeter element, then compute the eigenvectors, but this way is too tedious.
I learned from http://www.math.ubc.ca/~cass/research/pdf/Element.pdf
that the 2D plane $P$ can also spanned by two eigenvectors  $u',v'$ of the cartan matrix $C$. Here one computes the max and min eigenvaluesof $C$,and $u',v'$ are the eigenvectors of them.
But when I followed this way, the figure produced is symmetric but doesn't display the 8x30 ring pattern. I wonder whether there is something wrong with my work.
----------------------------------
The cartan matrix $C$:(without a factor 2)
$$\begin{pmatrix} 1. & -0.5 & 0. &  0. &  0.  & 0. &  0. &  0. \\
 -0.5 & 1.  &-0.5&  0. &  0. &  0.  & 0.  & 0. \\
 0. & -0.5 & 1.  & -0.5 & 0.&   0.&   0. &  0. \\
  0.  & 0.  &-0.5  &1. & -0.5 & 0. &  0. &  0. \\
  0.  & 0.  & 0.  &-0.5 & 1. & -0.5 & 0. & -0.5\\
 0.  & 0.  & 0.  & 0. & -0.5 & 1. & -0.5 & 0. \\
  0.  & 0.  & 0.  & 0.  & 0. & -0.5  &1.  & 0. \\
 0.  & 0.  & 0. &  0. & -0.5  &0.&   0.  & 1. \end{pmatrix}$$
The resulting figure:

 A: As so often happens, the answer is Steinberg.
R. Steinberg, Finite reflection groups. Trans. Amer. Math. Soc. 91 (1959) 493–504.
He explicitly constructs two lines that span the plane.  A while back, for my paper
Noncrossing partitions, clusters and the Coxeter plane, Sém. Lothar. Combin. 63 (2010), Article B63b.
I worked through this and wrote maple code to construct an orthonormal basis for the plane.  I looked and there are not many clues in my paper, but I still have the code.  (It's ugly...so I won't post it here, but send me an email and I can share it.  I don't know if mathoverflow has a messaging feature, but google my name and email me.)  Or, have more fun and work through it yourself.
A: I found the answer in a Python script that I translated to R. The main point is in the paragraph starting with "Now we split the simple roots into two disjoint sets I and J such that..." but personally I don't understand it.
# --- Step two: compute a basis of the Coxeter plane ---

# A set of simple roots listed by rows of 'delta'
delta <- rbind(
  c(1, -1, 0, 0, 0, 0, 0, 0),
  c(0, 1, -1, 0, 0, 0, 0, 0),
  c(0, 0, 1, -1, 0, 0, 0, 0),
  c(0, 0, 0, 1, -1, 0, 0, 0),
  c(0, 0, 0, 0, 1, -1, 0, 0),
  c(0, 0, 0, 0, 0, 1, 1, 0),
  c(-.5, -.5, -.5, -.5, -.5, -.5, -.5, -.5),
  c(0, 0, 0, 0, 0, 1, -1, 0)
)

# Dynkin diagram of E8:
# 1---2---3---4---5---6---7
#                 |
#                 8
# where vertex i is the i-th simple root.

# The Cartan matrix:
cartan = tcrossprod(delta)

# Now we split the simple roots into two disjoint sets I and J
# such that the simple roots in each set are pairwise orthogonal.
# It's obvious to see how to find such a partition given the
# Dynkin graph above: I = [1, 3, 5, 7] and J = [2, 4, 6, 8],
# since roots are not connected by an edge if and only if they are orthogonal.
# Then a basis of the Coxeter plane is given by
# u1 = sum (c[i] * delta[i]) for i in I,
# u2 = sum (c[j] * delta[j]) for j in J,
# where c is an eigenvector for the minimal
# eigenvalue of the Cartan matrix.
eig <- eigen(cartan)

# The eigenvalues returned by eigen() are in descending order
# and the eigenvectors are listed by columns.
ev <- eig$vectors[, 8L]
u1 <- rowSums(vapply(c(1L, 3L, 5L, 7L), function(i){
  ev[i] * delta[i, ]
}, numeric(8L)))
u2 <- rowSums(vapply(c(2L, 4L, 6L, 8L), function(i){
  ev[i] * delta[i, ]
}, numeric(8L)))

# Gram-Schmidt u1, u2 and normalize them to unit vectors.
u1 <- u1 / sqrt(c(crossprod(u1)))
u2 <- u2 - c(crossprod(u1, u2)) * u1
u2 <- u2 / sqrt(c(crossprod(u2)))

But this works:

