用tikz画球坐标系下的体积微元

最近上传热学，上来第一道作业就是球坐标系下导热微分方程的推导，我想着既然推导嘛，总得从画个好看的图开始😁就在网上找了一下球坐标系下体积微元的图，怎么说呢，大多数都不怎么清晰，而且考虑到可能还要在图中做标注，网上的图果然不是很满意。 

那就自己用咱们万能的绘图工具Tikz绘制一个？想法是可以的，但是3D的我之前完全没有画过啊，不管他，先去TEXample找找。3D的例子是挺多的，但是一个容易懂的代码才是最好的。

反复找了之后看到一个比较容易理解的：3dplot，就决定分享给大家了（可以在CTAN上下载：http://www.ctan.org/tex-archive/graphics/pgf/contrib/tikz-3dplot/）。

这个包的绘图思路包括主要是两种，一种是在主坐标系（tdplot_mian_coords）下的xOy平面绘图，另一种是在旋转后的坐标系（tdplot_rotated_coords）下的xOy平面内绘图。旋转的参数可以用命令\tdplotsetrotatedcoord{alpha}{beta}{gamma}由欧拉角确定（但要注意的是这里的欧拉角是'ZYZ'形式！在计算时要注意）。在坐标系coord下以O为圆心，r为半径画圆弧的命令是\tdplotdrawarc[coords]{(O)}{r}{θ}{θ+Δθ}。

画出来的效果如下图 ，嗯，怎么说，挺舒服的。

 附上画图的latex代码：

\documentclass{article}
\usepackage{verbatim}
\usepackage{tikz}  
\usepackage{3dplot} 
\usepackage[active,tightpage]{preview}
\PreviewEnvironment{tikzpicture}
\setlength\PreviewBorder{2mm}

\begin{document}

\tdplotsetmaincoords{60}{110}

\pgfmathsetmacro{\rvec}{.8}
\pgfmathsetmacro{\thetavec}{40}
\pgfmathsetmacro{\phivec}{45}

\pgfmathsetmacro{\dphivec}{20}
\pgfmathsetmacro{\dthetavec}{20}
\pgfmathsetmacro{\drvec}{.2}

\begin{tikzpicture}[scale=5,tdplot_main_coords]

\coordinate (O) at (0,0,0);

\tdplotsetcoord{A}{\rvec}{\thetavec}{\phivec}
\tdplotsetcoord{B}{\rvec}{\thetavec + \dthetavec}{\phivec}
\tdplotsetcoord{C}{\rvec}{\thetavec + \dthetavec}{\phivec + \dphivec}
\tdplotsetcoord{D}{\rvec}{\thetavec}{\phivec + \dphivec}

\tdplotsetcoord{E}{\rvec + \drvec}{\thetavec}{\phivec}
\tdplotsetcoord{F}{\rvec + \drvec}{\thetavec + \dthetavec}{\phivec}
\tdplotsetcoord{G}{\rvec + \drvec}{\thetavec + \dthetavec}{\phivec + \dphivec}
\tdplotsetcoord{H}{\rvec + \drvec}{\thetavec}{\phivec + \dphivec}

\draw[thick,->] (0,0,0) -- (0.8,0,0) node[anchor=north east]{$x$};
\draw[thick,->] (0,0,0) -- (0,0.8,0) node[anchor=north west]{$y$};
\draw[thick,->] (0,0,0) -- (0,0,0.8) node[anchor=south]{$z$};


\draw[-stealth,color=red] (O) -- (A) node[midway, anchor = south, color=black]{$r$};
\draw[color=red] (O) -- (B);
\draw[dashed, color=red] (O) -- (C);
\draw[dashed, color=red] (O) -- (D);

\draw[dashed, color=red] (O) -- (Axy);
\draw[dashed, color=red] (A) -- (Axy);
\draw[dashed, color=red] (O) -- (Dxy);
\draw[dashed, color=red] (D) -- (Dxy);
\draw[](A) -- (E) node[midway, anchor = south, color=black]{$dr$};
\draw[](B) -- (F);
\draw[dashed](C) -- (G);
\draw[dashed](D) -- (H);

\tdplotdrawarc{(O)}{0.2}{0}{\phivec}{anchor=north}{$\phi _0$}
\tdplotdrawarc{(O)}{0.35}{\phivec}{\phivec + \dphivec}{anchor=north}{$d\phi$}

\tdplotsetthetaplanecoords{\phivec}

\tdplotdrawarc[tdplot_rotated_coords]{(0,0,0)}{0.5}{0}{\thetavec}{anchor=south west}{$\theta _0$}
\tdplotdrawarc[tdplot_rotated_coords]{(0,0,0)}{0.3}{\thetavec}{\thetavec + \dthetavec}{anchor=south west}{$d\theta$}

\tdplotdrawarc[tdplot_rotated_coords]{(O)}{\rvec}{\thetavec}{\dthetavec + \thetavec}{}{}
\tdplotdrawarc[tdplot_rotated_coords]{(O)}{\rvec + \drvec}{\thetavec}{\dthetavec + \thetavec}{}{}

\tdplotsetthetaplanecoords{\phivec + \dphivec}
\tdplotdrawarc[dashed, tdplot_rotated_coords]{(O)}{\rvec}{\thetavec}{\dthetavec + \thetavec}{}{}
\tdplotdrawarc[tdplot_rotated_coords]{(O)}{\rvec + \drvec}{\thetavec}{\dthetavec + \thetavec}{}{}

\tdplotsetrotatedcoords{55}{-50.4313}{-6.4086}
\tdplotdrawarc[dashed, tdplot_rotated_coords]{(O)}{\rvec}{0}{12.8173}{}{}
\tdplotdrawarc[tdplot_rotated_coords]{(O)}{\rvec + \drvec}{0}{12.8173}{}{}

\tdplotsetrotatedcoords{55}{-30.3813}{-8.6492}
\tdplotdrawarc[dashed, tdplot_rotated_coords]{(O)}{\rvec}{0}{17.2983}{}{}
\tdplotdrawarc[tdplot_rotated_coords]{(O)}{\rvec + \drvec}{0}{17.2983}{}{}
\end{tikzpicture}

\end{document}

留心的朋友们应该注意到了，我在进行坐标旋转的时候，欧拉角是魔数！！！

其实这是因为它采用的ZYZ的欧拉角的缘故，这里如果是ZYX的欧拉角我们可以直接口算并用A点方位φ和θ以及D(r, θ, φ+dφ)来表示，但这个ZYZ型的，写解析解过于麻烦（这可能也是为什么我找不到别人贴这个体积微元画法的代码的缘故吧🤣），这里我是用matlab根据φ和θ计算得到的，也一并放上来。

function eul = getEul(theta, phi, dp)
    az=deg2rad(phi);
    el=deg2rad(90-theta);
    daz = deg2rad(dp);
    aa = [az az + daz];
    ee = [el el];
    rr = [1 1];
    [xx,yy,zz]=sph2cart(aa,ee,rr);
    p1 = [xx(1); yy(1); zz(1)];
    p2 = [xx(2); yy(2); zz(2)];
    e23 = cross(p1, p2);
    e23 = e23/norm(e23);
    e22 = cross(e23, p1);
    E0 = [1 0 0;
        0 1 0;
        0 0 1];
    rotM = [p1 e22 e23]/E0;
    eul = rotm2eul(rotM, 'ZYZ');
    eul = rad2deg(eul);
end

本来写的那版matlab代码还附加计算了∠AOD（这个在原来Latex代码里也是魔数）的，但想起来写博客的时候已经找不到了，就只草草写了一个算欧拉角的，毕竟现在知道了A（matlab代码中P1）、D（matlab代码中P2）算夹角也挺简单。

 

参考资料：

欧拉角和四元数：https://www.cnblogs.com/21207-iHome/p/6894128.html

TEXample-3dplot：http://www.texample.net/tikz/examples/the-3dplot-package/