Anisotropie Einkristall mit tikz

Hallo,

die folgende Frage wurde bereits im golatex-Forum gepostet, wo man mich hierher verwies.

Ich möchte die Anisotropie eines Einkristalls mittels des E-Moduls im dreidimensionalen Raum in Abhängigkeit der Kristallrichtungen darstellen.

Dazu soll zunächst die Elastizitätsmatrix von Silizium genutzt werden:

$\boldsymbol{C} = 
\begin{bmatrix}
C_{11} & C_{12} & C_{12} &0&0&0 \\ 
C_{12} & C_{11} & C_{12} &0&0&0 \\ 
C_{12} & C_{12} & C_{11} &0&0&0 \\ 
0&0&0&C_{44}&0&0\\  
0&0&0&0&C_{44}&0\\ 
0&0&0&0&0&C_{44}\\ 
\end{bmatrix}$

Elastischen Konstanten:

 $C_{11} = 165,7\,\mathrm{GPa}, \; C_{12} = 63,9\,\mathrm{GPa}, \; C_{44} = 79,6\,\mathrm{GPa}$

Diese Konstanten gelten für ein parallel zu den <100> -Richtungen ausgerichtetes Koordinatensystem. Geht man von einem einachsigen Spannungszustand in einer bestimmten Richtung aus, so kann man den E-Modul für diese Richtung berechnen. Führt man diese Rechnung für beliebige Richtungen durch und trägt die erhaltenen Werte als Abstand um einen Mittelpunkt in ihrer jeweiligen Richtung ab, so entsteht der angehangene Körper.

Zur Umsetzung bin ich auf das Paket tikz-3dplot gestoßen. Dazu hangele ich mich an folgendem Minimalbeispiel entlang.

\documentclass{minimal}
\usepackage{verbatim}
\usepackage{tikz}
\usepackage{3dplot}

\begin{document}
\tdplotsetmaincoords{70}{135}
\begin{tikzpicture}[scale=1,line join=bevel,tdplot_main_coords, fill opacity=.7]
\pgfsetlinewidth{.1pt}
\tdplotsphericalsurfaceplot[parametricfill]{72}{36}%
{sqrt(15/2)*sin(\tdplottheta)*cos(\tdplottheta)}{black}{\tdplotphi}%
{\draw[color=black,-latex] (0,0,0) -- (2,0,0) node[anchor=north east]{$x$};}%
{\draw[color=red,-latex] (0,0,0) -- (0,2,0) node[anchor=north west]{$y$};}%
{\draw[color=black,-latex] (0,0,0) -- (0,0,2) node[anchor=south]{$z$};}%
\end{tikzpicture}
\end{document}

Nun ist jedoch die Frage, wie ich meine Elastizitätsmatrix hier implementiere und diese dann auch zum Rotieren um den Koordinatenursprung bringe.

Kann hier jemand einen Tip geben bzw. gibt es alternative Lösungsmöglichkeiten?

pgfplots kenne ich ganz gut da ich jeden Tag damit arbeite. Bei der vorliegenden Problemstellung war es mir jedoch noch nicht behilflich.

Zur Formel... Komponenten der Steifigkeitstetrade bezogen auf eine Orthonormalbasis:

$\boldsymbol{C} = C_{ijkl} \boldsymbol{e}_i \otimes \boldsymbol{e}_j \otimes \boldsymbol{e}_k \otimes \boldsymbol{e}_l$

Eine Rotation des Koordinatensystems mit einem Winkel $\omega$ um eine beliebige Achse, deren Lage durch den Einheitsvektor $\boldsymbol{e}$ festgelegt ist, kann durch folgenden (orthogonalen) Tensor vorgenommen werden:

$\boldsymbol{Q} = \boldsymbol{1} \mathrm{cos}\omega +\boldsymbol{e} \otimes \boldsymbol{e} (1-\mathrm{cos}\omega) -\boldsymbol{1} \times \boldsymbol{e} \mathrm{sin}\omega$

Einsteinsche Summenkonvention gilt, $\otimes$ - dyadisches Produkt, $\times$ - Kreuzprodukt

In Mathematica hat es schon einmal jemand umgesetzt. Ein Code-Beispiel gibt es dort auch: Anisotropic Elasticity

Ich habe es mal nachgebastelt: Seite 1 Seite 2

Vielleicht wird die Mathematik hierbei klarer.

Eine exemplarische Plot-Funktion sollte lauten:

63.19/(1.*Cos[theta]^4 + 1.51*Cos[phi]^2*Cos[theta]^2*Sin[theta]^2 + 1.*Cos[phi]^4*Sin[theta]^4 + 0.99*Sin[phi]^4*Sin[theta]^4 + 0.37*Sin[2*phi]^2*Sin[theta]^4 + 0.38*Sin[phi]^2*Sin[2*theta]^2)