Skip to content

Instantly share code, notes, and snippets.

@chakpongchung

chakpongchung/master_project.ipynb

Created November 5, 2015 17:07
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save chakpongchung/614b1a0606364d160e26 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save chakpongchung/614b1a0606364d160e26 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
master_project.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Derivation for the inner product of 3D rotational gradient operator acted on arbitrary spherical harmonics \n",
"\n",
"Chak-Pong CHUNG, [email protected]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To the best of our knowledge, it is the first time this dot product is derived. As you can compare, now the computation is much more simpler since the number of integral within the unit sphere(the alpha terms) are only 1/3 of the original"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with application in object/scenes response to a lighting environment, often represented using spherical harmonics,including complex global illumination effects.Other application includes simulation of crystal growth in material science."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first paper that\n",
"has been used extensively in games deals with using Spherical Harmonics to represent\n",
"irradiance environment maps efficiently, allowing for interactive rendering of diffuse objects\n",
"under distant illumination"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Reference about spherical harmonics\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"http://mathworld.wolfram.com/SphericalHarmonic.html"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#Spherical harmonics in Sympy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To get back the well known expressions in spherical coordinates we use full expansion:"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"sqrt((2*n + 1)*factorial(-m + n)/factorial(m + n))*exp(I*m*phi)*assoc_legendre(n, m, cos(theta))/(2*sqrt(pi))"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from sympy import Ynm, Symbol, expand_func\n",
"from sympy.abc import n,m\n",
"theta = Symbol(\"theta\")\n",
"phi = Symbol(\"phi\")\n",
"expand_func(Ynm(n, m, theta, phi))"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Relevant page\n",
"http://docs.sympy.org/latest/modules/functions/special.html?highlight=spherical%20harmonics#sympy.functions.special.spherical_harmonics.Ynm"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"3d plot of spherical harmonics\n",
"\n",
"http://stsdas.stsci.edu/download/mdroe/plotting/entry1/index.html"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"# Derivation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We try to derive the dot product of rotational gradient of spherical harmonics with spherical harmonics expansion."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"$$\\vec{R} = \\begin{bmatrix}\n",
" R_{x}\\\\\n",
" R_{y} \\\\\n",
" R_{z}\n",
" \\end{bmatrix}\n",
" =\\begin{bmatrix}\n",
" - \\frac{cos\\theta}{sin\\theta}cos\\phi \\frac{\\partial }{\\partial \\phi} - sin\\phi \\frac{\\partial }{\\partial \\theta}\n",
" \\\\\n",
" - \\frac{cos\\theta}{sin\\theta}cos\\phi \\frac{\\partial }{\\partial \\phi} + cos\\phi \\frac{\\partial }{\\partial \\theta} \n",
" \\\\\n",
" \\frac{\\partial }{\\partial \\phi} \n",
" \\end{bmatrix}\n",
" $$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"R_{\\alpha} = i L_{\\alpha}, where \\, \\alpha = x,y,z\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" $$ \\vec{R} = \\begin{bmatrix}\n",
" R_{x} \\\\\n",
" R_{y} \\\\\n",
" R_{z}\n",
" \\end{bmatrix}\n",
" =\n",
" \\begin{bmatrix}\n",
" i L_{x} \\\\\n",
" i L_{y} \\\\\n",
" i L_{z}\n",
" \\end{bmatrix}\n",
" =\n",
" \\begin{bmatrix}\n",
" \\frac{i}{2}(L_{+}+L_{-} ) \\\\\n",
" \\frac{1}{2}(L_{+}-L_{-} ) \\\\\n",
" i L_{z}\n",
" \\end{bmatrix}\n",
" $$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"L_{+} Y{_l^{m}}= c_{l,m,1} Y{_l^{m+1}}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"$$\n",
"L_{-} Y{_l^{m}}= c_{l,m,2} Y{_l^{m-1}}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"$$\n",
"L_{z} Y{_l^{m}}= m Y{_l^{m}}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"after some manipulations of the above identities, it can be shown that (details will be added later)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ \\vec{R}Y{_j^{p}} \\cdot \\vec{R}Y{_l^{m}} = - \\frac{1}{2} [c_{j,p,1}c_{l,m,2}Y{_j^{p+1}} Y{_l^{m-1}}+c_{j,p,2}c_{l,m,1}Y{_j^{p-1}}Y{_l^{m+1}}] - mp Y{_j^{p}} Y{_l^{m}} $$ "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ =(-1)*(-1)^{m+p} \\sum\\limits_{k=|l-j|}^{l+j}Y{_k^{m+p}}[\\frac{1}{2} (c_{j,p,1}c_{l,m,2} \\alpha_{l,m-1,j,p+1,k} +c_{j,p,2}C_{l,m,1}\\alpha_{l,m+1,j,p-1,k}) + mp \\alpha_{l,m,j,p,k}] $$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $$ \\alpha_{l,m,j,p,k} = \\int_{0}^{2\\pi} \\int_{0}^{\\pi} Y{_j^{p}} Y{_l^{m}} Y{_k^{-m-p}} \\sin \\theta d\\theta d\\phi$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ =[\\frac{(2l+1)(2j+1)(2k+1)}{4\\pi}]^{1/2} * \\left(\\begin{array}{clcr} j & l & k\\\\ 0& 0 & 0 \\end{array}\\right) * \\left(\\begin{array}{clcr} j & l & k\\\\ p& m & -(m+p) \\end{array}\\right) $$ "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and \n",
"$$\n",
"\\left(\\begin{array}{clcr} j & l & k\\\\ 0& 0 & 0 \\end{array}\\right)\n",
"$$\n",
"\n",
"is Wigner 3j symbol"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This formula is computationally expensive since it calculates different integral in the unit sphere(the alpha term) many times. The heaviness of the computation comes from evaluating the Wigner 3j symbol, which involves calculating the factorial, a notoriusly expensive calculations. It also introduces rounding errors."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Wigner 3j symbol"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"reference for Wigner 3j symbol(strange fraction here)\n",
"\n",
"http://docs.sympy.org/dev/modules/physics/wigner.html"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"sqrt(715)/143"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from sympy.physics.wigner import wigner_3j\n",
"wigner_3j(2, 6, 4, 0, 0, 0)"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.18698939800169145"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import math\n",
"math.sqrt(S(5)/143)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"wigner_3j(2, 6, 4, 0, 0, 1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Main ideas of the derivation"
]
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"To simplify the whole derivation, these notations are introduced"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
" Y{_j^{p}} = Y{_j^{p}}(\\theta, \\phi)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"c_{l,m,1}=\\sqrt{(l-m)(l+m+1)}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"c_{l,m,2}=\\sqrt{(l+m)(l-m+1)}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It can be shown that (details will be added later) ,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" $$ \\frac{1}{2} (c_{j,p,1}c_{l,m,2} \\alpha_{l,m-1,j,p+1,k} +c_{j,p,2}C_{l,m,1}\\alpha_{l,m+1,j,p-1,k}) + mp \\alpha_{l,m,j,p,k} $$ "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ =\\frac{1}{2}[\\frac{(2l+1)(2j+1)(2k+1)}{4\\pi}]^{1/2} * \\left(\\begin{array}{clcr} j & l & k\\\\ 0& 0 & 0 \\end{array}\\right) * [ c_{j,p,1}c_{l,m,2} \\left(\\begin{array}{clcr} j & l & k\\\\ p+1& m-1 & -(m+p) \\end{array}\\right) + c_{j,p,2}c_{l,m,1} \\left(\\begin{array}{clcr} j & l & k\\\\ p-1& m+1 & -(m+p) \\end{array}\\right) + 2mp \\left(\\begin{array}{clcr} j & l & k\\\\ p-1& m+1 & -(m+p) \\end{array}\\right) ] $$ "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ = \\frac{1}{2}[\\frac{(2l+1)(2j+1)(2k+1)}{4\\pi}]^{1/2} * \\left(\\begin{array}{clcr} j & l & k\\\\ 0& 0 & 0 \\end{array}\\right) * \\left(\\begin{array}{clcr} k & j & l\\\\ -(m+p)& p & m \\end{array}\\right) *(k^2-j^2-l^2+k-j-(-2pm+2pm)) $$ "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" $$ = \\frac{1}{2}[\\frac{(2l+1)(2j+1)(2k+1)}{4\\pi}]^{1/2} * \\left(\\begin{array}{clcr} k & j & l\\\\ 0& 0 & 0 \\end{array}\\right) * \\left(\\begin{array}{clcr} k & j & l\\\\ -(m+p)& p & m \\end{array}\\right) *(k^2-j^2-l^2+k-j-l) $$ "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ = \\alpha_{l,m,j,p,k} * \\frac{1}{2} (k^2-j^2-l^2+k-j-l) $$ "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"so, we have \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ \\vec{R}Y{_j^{p}} \\cdot \\vec{R}Y{_l^{m}} = (-1)^{m+p} \\sum\\limits_{k=|l-j|}^{l+j}Y{_k^{m+p}} * \\alpha_{l,m,j,p,k} * \\frac{1}{2} (k^2-j^2-l^2+k-j-l) $$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To the best of our knowledge, it is the first time this dot product is derived. As you can compare, now the computation is much more simpler since the number of integral within the unit sphere(the alpha terms) are only 1/3 of the original"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Testing with SymPy"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0\n"
]
}
],
"source": [
"from sympy import *\n",
"from sympy.physics.wigner import wigner_3j\n",
"init_printing(use_latex='mathjax')\n",
"\n",
"def dot_rota_grad_SH(j, p, l, m, theta, phi):\n",
" temp=0\n",
" for k in range(abs(l-j), l+j +1):\n",
" temp+=Ynm(k, m+p, theta,phi)*alpha(l,m,j,p,k)/2*(k**2-j**2-l**2+k-j-l)\n",
"\n",
" return (-S(1))**(m+p)*temp\n",
"\n",
"def alpha(l,m,j,p,k):\n",
" return sqrt((2*l+1)*(2*j+1)*(2*k+1)/(4*pi))*wigner_3j(j, l, k, 0, 0, 0)*wigner_3j(j, l, k, p, m, -m-p)\n",
"\n",
"\n",
"print dot_rota_grad_SH(2,6,2,2,0,0)"
]
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"This dot_rota_grad_SH() and the following are now incorporated into SymPy"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from sympy import Ynm, init_printing, Expr, var, sympify, Dummy, S, Sum, sqrt, pi\n",
"from sympy.physics.wigner import wigner_3j\n",
"# init_printing(use_latex='mathjax')"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"class Wigner3j(Expr):\n",
" def doit(self, **hints):\n",
" num = True\n",
" for i in range(6):\n",
" if not self.args[i].is_number:\n",
" num = False\n",
" if num:\n",
" return wigner_3j(*self.args)\n",
" else:\n",
" return self\n",
"\n",
"def alpha(l,m,j,p,k):\n",
" return sqrt((2*l+1)*(2*j+1)*(2*k+1)/(4*pi)) * Wigner3j(j, l, k, S(0), S(0), S(0))*Wigner3j(j, l, k, p, m, -m-p)\n",
"\n",
"def dot_rota_grad_SH(j, p, l, m, theta, phi):\n",
" j = sympify(j)\n",
" p = sympify(p)\n",
" l = sympify(l)\n",
" m = sympify(m)\n",
" theta = sympify(theta)\n",
" phi = sympify(phi)\n",
" k = Dummy(\"k\")\n",
" return (-S(1))**(m+p) * Sum(Ynm(k, m+p, theta, phi)*alpha(l,m,j,p,k)/2 *(k**2-j**2-l**2+k-j-l), (k, abs(l-j), l+j))"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"-Sum(sqrt(50*_k + 25)*(_k**2 + _k - 12)*Ynm(_k, 1, theta, phi)*Wigner3j(2, 2, _k, 0, 0, 0)*Wigner3j(2, 2, _k, 0, 1, -1)/(4*sqrt(pi)), (_k, 0, 4))"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"var(\"j p l m theta phi\")\n",
"#dot_rota_grad_SH(1, 5, 1, 1, 1, 2)\n",
"dot_rota_grad_SH(2,0,2,1,theta,phi)"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-3*sqrt(5)*Ynm(2, 1, theta, phi)/(14*sqrt(pi)) + 2*sqrt(30)*Ynm(4, 1, theta, phi)/(7*sqrt(pi))\n"
]
}
],
"source": [
"print _.doit()"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"-Sum(_k**2*sqrt(50*_k + 25)*sqrt(2*_k*factorial(_k - 1)/factorial(_k + 1) + factorial(_k - 1)/factorial(_k + 1))*exp(I*phi)*assoc_legendre(_k, 1, cos(theta))*Wigner3j(2, 2, _k, 0, 0, 0)*Wigner3j(2, 2, _k, 0, 1, -1)/(8*pi) + _k*sqrt(50*_k + 25)*sqrt(2*_k*factorial(_k - 1)/factorial(_k + 1) + factorial(_k - 1)/factorial(_k + 1))*exp(I*phi)*assoc_legendre(_k, 1, cos(theta))*Wigner3j(2, 2, _k, 0, 0, 0)*Wigner3j(2, 2, _k, 0, 1, -1)/(8*pi) - 3*sqrt(50*_k + 25)*sqrt(2*_k*factorial(_k - 1)/factorial(_k + 1) + factorial(_k - 1)/factorial(_k + 1))*exp(I*phi)*assoc_legendre(_k, 1, cos(theta))*Wigner3j(2, 2, _k, 0, 0, 0)*Wigner3j(2, 2, _k, 0, 1, -1)/(2*pi), (_k, 0, 4))"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"_.expand(func=True)"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0 0 0 0 0\n",
"0 0 1 -1 0\n",
"0 0 1 0 0\n",
"0 0 1 1 0\n",
"1 -1 0 0 0\n",
"1 -1 1 -1 sqrt(30)*exp(-4*I*phi)*Ynm(2, 2, theta, phi)/(10*sqrt(pi))\n",
"1 -1 1 0 -sqrt(15)*exp(-2*I*phi)*Ynm(2, 1, theta, phi)/(10*sqrt(pi))\n",
"1 -1 1 1 Ynm(0, 0, theta, phi)/sqrt(pi) + sqrt(5)*Ynm(2, 0, theta, phi)/(10*sqrt(pi))\n",
"1 0 0 0 0\n",
"1 0 1 -1 -sqrt(15)*exp(-2*I*phi)*Ynm(2, 1, theta, phi)/(10*sqrt(pi))\n",
"1 0 1 0 -Ynm(0, 0, theta, phi)/sqrt(pi) + sqrt(5)*Ynm(2, 0, theta, phi)/(5*sqrt(pi))\n",
"1 0 1 1 sqrt(15)*Ynm(2, 1, theta, phi)/(10*sqrt(pi))\n",
"1 1 0 0 0\n",
"1 1 1 -1 Ynm(0, 0, theta, phi)/sqrt(pi) + sqrt(5)*Ynm(2, 0, theta, phi)/(10*sqrt(pi))\n",
"1 1 1 0 sqrt(15)*Ynm(2, 1, theta, phi)/(10*sqrt(pi))\n",
"1 1 1 1 sqrt(30)*Ynm(2, 2, theta, phi)/(10*sqrt(pi))\n"
]
}
],
"source": [
"for j in range(2):\n",
" for p in range(-j, j+1):\n",
" for l in range(2):\n",
" for m in range(-l, l+1):\n",
" print j, p, l, m, dot_rota_grad_SH(j, p, l, m, theta, phi).doit()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment