---
title: "Spacekime Analytics (Time Complexity and Inferential Uncertainty)"
author: 'SOCR Team '
date: "`r format(Sys.time(),'%m/%d/%Y')`"
output:
  html_document:
    theme: spacelab
    highlight: tango
    includes:
      before_body: TCIU_header.html
    toc: yes
    number_sections: yes
    toc_depth: 2
    toc_float:
      collapsed: no
      smooth_scroll: yes
    code_folding: hide
  word_document:
    toc: yes
    toc_depth: '2'
  pdf_document:
    toc: yes
    toc_depth: '2'
always_allow_html: yes
subtitle: "Chapter 6: Circular Kime-Phase and Electron Orbit Densities (Circular Distributions)"
---

```{r error=F, message=F, warning=F, }
# install.packages("circular")
# library(crosstalk)
# crosstalkLibs()
library(htmlwidgets)
library(circular)
library(animation)
library(ggplot2)
library(tidyr)
library(plotly)
```

In this [TCIU/Spacekime section](https://tciu.predictive.space/), we show some of the parallels between circular kime-phase densities and [electron orbital densities](https://chem.libretexts.org/Bookshelves/General_Chemistry/Map%3A_Chemistry_-_The_Central_Science_(Brown_et_al.)/06%3A_Electronic_Structure_of_Atoms/6.06%3A_3D_Representation_of_Orbitals) depicted as [2D circular or 3D spherical distributions](https://winter.group.shef.ac.uk/orbitron).

# Random sampling from a given probability distribution

By definition, experimental science and confirmatory evidence decision-making are based on observations, which are typically assumed to be identically-distributed and independent (IID) *random* samples from an unknown distribution. Observable phenomena (or processes) may be modeled by simple *primary distributions*, or by composite superpositions (or mixtures) of several primary (base) distributions.

The animation below illustrates the [process of drawing a *random sample* from any univariate distribution](https://wiki.socr.umich.edu/index.php/SOCR_EduMaterials_Activities_GeneralCentralLimitTheorem). Analogous IID sampling applies for drawing random observations from multiple (multivariate) distributions, see [SOCR bivariate and trivariate distribution webapp](https://socr.umich.edu/HTML5/BivariateNormal/BVN2/).

<video autosize:true controls>
  <source src="SOCR_CLT_App_Video.mp4" type="video/mp4">
</video>

Fundamentally, all *data, information*, and *evidence* is always finite, quantized
as individual or (order-agnostic) combinations of countable observations that can be indexed, tracked, and referenced one-by-one or as subsets. 

In contrast to observed samples, *distributions* correspond to *observable processes* that can be sampled to obtain *instances of data*. However, distributions are holistic analytic models that cannot be themselves observed as data. Only finite samples drawn from probability distributions are observable. Hence, all evidence-based inference is based on *proxy data* aiming to elucidate the properties of the underlying probability distributions. The first two laws of probability theory explain why functions of *random (IID) samples* (statistics!) yield potent *proxy measures* of the intrinsic properties of the source distribution, e.g., average $\leftrightarrow$ mean, sample variance $\leftrightarrow$ population variance, etc.

 - *First fundamental law of probability theory*: The [Law of Large Numbers
 (LLN)](https://doi.org/10.1080/10691898.2009.11889499) reflects inference based on repeated IID sampling of a controlled experiment. For instance, if we are interested in the *relative frequency* of occurrence of one event whose *probability* to be observed at each experiment is $p$. Then, the ratio of the observed *sample frequency* of that event relative to the total number of repetitions converges towards $p$ as the number of (identical and independent) experiments increases.
 - *Second fundamental law of probability theory*: The [Central Limit Theorem (CLT)](https://www.tandfonline.com/doi/full/10.1080/10691898.2008.11889560) yields that the (arithmetic) average (of random variables sampled from a process with well-defined and finite first two moments) will follow approximately a normal distribution, as the sample-size increases.


<!-- **Note**: Perhaps connect with the *sufficiency & exchangeability* paper notions of random sampling, IID and De Finetti theorem ... -->


# Circular Electron Orbit Densities 

In many physics systems and chemistry applications, it's important to model the [spatial distribution of elementary particles, such as electrons](https://chem.libretexts.org/Bookshelves/General_Chemistry/Map%3A_Chemistry_-_The_Central_Science_(Brown_et_al.)/06%3A_Electronic_Structure_of_Atoms/6.06%3A_3D_Representation_of_Orbitals). For instance, in photoelectron spectroscopy experiments, we may measure the intensities of electrons to track their angular dependencies, which may be described as electric dipole approximations. These estimates may subsequently be used to study the electron emissions resulting from their interaction with photons. Suppose we measure the electron emission intensity, $I(\theta)$, the [Cooper-Zare's *Angular Distribution of Photoelectrons* model](https://doi.org/10.1063/1.1668742) suggests that $I(\theta)$ can be expressed as a differential of the partial cross section:

$$I(\theta)\equiv \frac{d\sigma}{d \Omega}=
\frac{\sigma}{4\pi}\left ( 1+ \frac{\beta}{2}(3\cos^2(\theta)-1) \right ).$$

This ordinary differential equation describes the angular distribution of the emitted electrons from
a gas-phase sample in a random orientation, excited by $100\%$ linearly polarized light, see [this reference](https://www.researchgate.net/publication/356920774_Dynamique_electronique_dans_les_molecules_organiques_par_la_spectroscopie_Core-Hole_Clock). The $\theta$ dependence correspond to spherical harmonics $Y_2^0$.

In the above *Cooper-Zare model* equation, 

 - $-\pi\leq \theta<\pi$ is the angle between the electric field vector and the vector that describes the *direction* (phase) of the electron ejected after the collision of the particle (electron) with a photon. In other words, $\theta$ measures the *angle between the direction of the ejected electron and the polarization of the incident light*.
 - $\sigma$ is the angle integrated *total cross-section* of a given particle state, and
 - $-1\le \beta\le 2$ is the [photoelectron asymmetry parameter](https://journals.aps.org/rmp/pdf/10.1103/RevModPhys.54.389), which is independent of $\theta$, but can depend on the excitation energy or the particle's (electrons’) kinetic energy. The $\beta$ term derivations are provided in [Introduction to Photoelectron Angular Distributions](https://link.springer.com/book/10.1007/978-3-031-08027-2). For a level $l$ (e.g., $s,p,d$), the range of $\beta$ is from $-1$, for $\frac{d\sigma}{d\Omega}\approx \sin^2\theta$, to $+2$, for $\frac{d\sigma}{d\Omega}\approx \cos^2\theta$:

$$\beta=\frac{l(l-1)\sigma_{l-1}^2 +(l+1)(l+2)\sigma_{l+1}^2 -
6l(l+1)\sigma_{l+1}\sigma_{l-1}\cos(\delta_{l+1}-\delta_{l-1})}
{(2l+1)(l \sigma _{l-1}^2+(l+1)\sigma _{l+1}^2 )}\ ,$$

 - the *phase-shift* of the $l^{th}$ partial wave is $\delta_l$, 
 - the *radial momentum* $R_{n\ l}(r)$ is normalized by $\int_{0}^{\infty}{[R_{n\ l}(r)]^2r^2\ dr} =1$, 
 - the direction of the *photoelectron momentum* is $G_{k\ l}(r)=\frac{e^{i\delta_l} \sin(kr-l\pi/2) +e^{-i(kr+i\pi/2)} e^{-2i\delta_l}\sin(\delta_l)}{kr}$,
 - the *radial dipole integral* is $R_{n}=\int_{0}^{\infty}{R_{n\ l}(r)\cdot G_{k\ l}(r)\cdot r^3\ dr}$, and
 - the usual *dipole radial matrix element* is 
 
 $$\sigma_{l\pm 1}=\int_{0}^{\infty}{R_{n\ l}(r) R_{k\ l\pm 1}(r)\ r\ dr} .$$

The graph below shows the probabilistic photoelectron differential cross-section showing
the photoionization of an atomic $1s$ orbital. An electron in this orbital has angular momentum $l=0$. Upon the electron-photon collision, the angular momentum of the electron increases to $l = 1$, according to the dipole selection rules ($\Delta l = \pm 1$). The ejected electron leaves the system as a *p-wave*, corresponding to the case of $\beta = 2$. In this case, the intensity in the polarization vector plane is maximized and close to zero in the plane perpendicular to the polarization vector. Similar arguments hold for the orbital distributions of electrons with different characters and their relations to the corresponding $\beta$ parameters.

Note the connections between all four angular distributions (corresponding to $\beta\in \{-1, 0, 1, 2 \}$) in four symmetrical points, where the angular distributions are independent of the
$\beta$-values. For instance, the first point corresponding to $\theta=54.7^o$ (see black line in the first quadrant, or $0.954695$ radians) is called the *magic angle* of the polarization vector of electromagnetic radiation. This *magic* is due to the intensity distribution being independent of the wave-character of the emitted electron. Assuming a perfect, $100\%$, linear polarization, for an electron energy detector mounted at the *magic angle*, the angular distribution vanishes and the differential partial cross-section in
the above equation is proportional to the integral partial cross-section. When the particle analyzing detector is focused at a different angle, the correction factors (dependencies) need to be accounted for due to angular distribution effects.

The significance of the $\beta$ parameter can be explicated in terms of an isolated carbon ($C$) atom. At all energies, an idealized and isolated $C\ 1s$ photoelectron exhibits angular distribution corresponding to $\beta = 2$ (green density). The polar angle (phase) represents the angle between the polarization vector (horizontal axis) and the emitted electron. The (green) plot shows that the electrons are preferentially emitted along the polarization vector. At one extreme, when $\beta = 0$, the angular distribution is completely isotropic (circular, as it's independent of $\theta$). In the other extreme, when $\beta = -1$, the maximum of the distribution is mostly vertical, i.e., perpendicular to the horizontal polarization vector. This interpretation is valid only under the dipole approximation assumption when the wavelength of light is much larger than the dimension of the atoms (the wave is not varying inside the atom). For instance, experiments frequently use small radiation wavelengths, around $2\ nm$, whereas the entire atom may be $10$ times smaller, at about $0.2\ nm$.

The *dipole approximation* refers to the *first order* approximation of the complex exponential when the wavelength of electromagnetic radiation, which transitions between different atomic energy levels, is much larger than the typical size of a light atom. Specifically,

$$e^{i(\omega/c)\ n\cdot x} = 1 + i\frac{\omega}{c} \ n\cdot x + \cdots 
\approx  1  ,$$

since $\frac{\omega}{c} =2\frac{\pi}{\lambda}$, i.e., 
the *wavelength* $\lambda\uparrow$ and the *frequency* $\omega\downarrow$ are inversely proportional, $\lambda=2\pi\frac{c}{\omega}$.

This electron distribution map suggests that in general, experiments are expected to yield more observed electrons when performing a measurement at $0^o$ angle and much fewer electrons when measuring at $90^o$. However, in practice, the difference in the observed electron densities compared to model predictions between the $0^o$ and $90^o$ angles is smaller than expected. This difference is a function of the photon energy. See [this publication demonstrating how the observed intensities can be used to estimate the $\beta$ parameter](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.106.193009#fulltext) for the *chlorine* substituted carbon and for the *methyl* carbon in a series of chloroethanes. The more chlorine atoms are present in a molecule, the [lower $\beta$ parameter is, and the deviation $2-\beta$ is larger at the lowest ionization energies](https://journals.aps.org/rmp/pdf/10.1103/RevModPhys.54.389).

Note that the angular distributions of these particles (electrons) are very structured.
They tend to be symmetric and do not come out independently all over the place. Rather, the distributions emerge in just a few directions, as narrow sprays or jets. About $90\%$ of the time these distributions represent just two polar opposite jets, about $10\%$ of the time there are three jets, $1\%$ four jets, and so on progressively less likely.

```{r electron_densities, warning=F, error=F, message=F}
# Display the electron densities for
theta <- seq(-pi, pi, length.out=200)

# beta = is the asymmetry parameter, independent of \theta, can depend on the excitation
# energy or the electrons’ kinetic energy
# theta is the angle integrated cross-section of a given state
f <- function (theta, beta = -1, sigma=1) {
    return ( (sigma/(4*pi))*(1 + (beta/2) * (3*cos(theta) * cos(theta) - 1)) )
}

df <- data.frame(cbind(theta=theta,
                       beta_1 = f(theta=theta, beta=-1, sigma=1),
                       beta0 = f(theta=theta, beta=0, sigma=1),
                       beta1 = f(theta=theta, beta=1, sigma=1),
                       beta2 = f(theta=theta, beta=2, sigma=1)))

df_long <- gather(df, beta, value, c(beta_1,beta0, beta1,beta2), factor_key=TRUE)

plot1 <- ggplot(df_long, aes(x = theta, y = value, colour = beta)) +
    geom_line() +
    geom_vline(xintercept=0.954695) +   # 54.7^o = 0.954695 radians
    labs(title = "", x = "", y = "")+
    # theme_bw() +
    geom_label(label="54.7 degrees", x=1.0, y=0.2,
               label.padding = unit(0.55, "lines"), # Rectangle size around label
               label.size = 0.35,  color = "black", fill="lightblue") +
    ggtitle("Electron angular density") +
    coord_polar(start = pi/2, direction = -1) +
    scale_color_manual(labels = c("beta=-1", "beta=0", "beta=1", "beta=2"),
                       values = c("red", "gray", "blue", "green")) +
    theme(axis.text.y = element_blank(), plot.title = element_text(hjust=0.5),
          axis.ticks.y = element_blank(), legend.position="bottom") +
    guides(color=guide_legend("Beta\nValues"))
# ggplotly(plot1)
plot1
```

# Kime-Phases Circular Distributions 

```{r kime_phases, eval=F, echo=T, warning=F, error=F, message=F}
library(animation)
saveHTML(
  {
  oopt = ani.options(interval = 0.2)
  set.seed(1234)
  x <- rvonmises(n=1000, mu=circular(pi/5), kappa=3)
  y <- rvonmises(n=1000, mu=circular(-pi/3), kappa=5)
  z <- rvonmises(n=1000, mu=circular(0), kappa=10)
  r <- runif(n=1000, min=0, max=2)
  vectorX = as.vector(x)
  vectorY = as.vector(y)
  vectorZ = as.vector(z)
  plotX = c(vectorX[1])
  plotY = c(vectorY[1])
  plotZ = c(vectorZ[1])
  for (i in 2:1000){
    plotX = c(plotX,vectorX[i])
    plotY = c(plotY,vectorY[i])
    plotZ = c(plotZ,vectorZ[i])
    
    # Try to "stack=T the points ....
    #r1 = sqrt((resx$x)^2+(resx$y)^2)/2; 
    #resx$x = r1*cos(resx$data)
    #resx$x = r1*cos(resx$data)
    
    tempX = circular(plotX)
    tempY = circular(plotY)
    tempZ = circular(plotZ)
    resx <- density(tempX, bw=25, xaxt='n', yaxt='n')
    resy <- density(tempY, bw=25, xaxt='n', yaxt='n')
    resz <- density(tempZ, bw=25, xaxt='n', yaxt='n')
    res <- plot(resx, points.plot=TRUE, xlim=c(-1.1,1.5), ylim=c(-1.5, 1.5), 
                main="Trivariate random sampling of\n kime-magnitudes (times) and kime-directions (phases)", 
                offset=0.9, shrink=1.0, ticks=T, lwd=3, axes=F, stack=TRUE, bins=150)
    lines(resy, points.plot=TRUE, col=2, points.col=2, plot.info=res, offset=1.0, shrink=1.45, lwd=3, stack=T)
    lines(resz, points.plot=TRUE, col=3, points.col=3, plot.info=res, offset=1.1, shrink=1.2, lwd=3, stack=T)
    segments(0, 0, r[i]*cos(vectorX[i]), r[i]*sin(vectorX[i]), lwd=2, col= 'black')
    segments(0, 0, r[i]*cos(vectorY[i]), r[i]*sin(vectorY[i]), lwd=2, col= 'red')
    segments(0, 0, r[i]*cos(vectorZ[i]), r[i]*sin(vectorZ[i]), lwd=2, col= 'green')
    ani.pause()
  }
  
  
  # To save animation as MP4 video
  saveVideo({ 
    oopt = ani.options(interval = 0.2)
    set.seed(1234)
    x <- rvonmises(n=1000, mu=circular(pi/5), kappa=3)
    y <- rvonmises(n=1000, mu=circular(-pi/3), kappa=5)
    z <- rvonmises(n=1000, mu=circular(0), kappa=10)
    r <- runif(n=1000, min=0, max=2)
    vectorX = as.vector(x)
    vectorY = as.vector(y)
    vectorZ = as.vector(z)
    plotX = c(vectorX[1])
    plotY = c(vectorY[1])
    plotZ = c(vectorZ[1])
    for (i in 2:1000){
      plotX = c(plotX,vectorX[i])
      plotY = c(plotY,vectorY[i])
      plotZ = c(plotZ,vectorZ[i])
      
      # Try to "stack=T the points ....
      #r1 = sqrt((resx$x)^2+(resx$y)^2)/2; 
      #resx$x = r1*cos(resx$data)
      #resx$x = r1*cos(resx$data)
      
      tempX = circular(plotX)
      tempY = circular(plotY)
      tempZ = circular(plotZ)
      resx <- density(tempX, bw=25, xaxt='n', yaxt='n')
      resy <- density(tempY, bw=25, xaxt='n', yaxt='n')
      resz <- density(tempZ, bw=25, xaxt='n', yaxt='n')
      res <- plot(resx, points.plot=TRUE, xlim=c(-1.1,1.5), ylim=c(-1.5, 1.5), 
                  main="Trivariate random sampling of\n kime-magnitudes (times) and kime-directions (phases)", 
                  offset=0.9, shrink=1.0, ticks=T, lwd=3, axes=F, stack=TRUE, bins=150)
      lines(resy, points.plot=TRUE, col=2, points.col=2, plot.info=res, offset=1.0, shrink=1.45, lwd=3, stack=T)
      lines(resz, points.plot=TRUE, col=3, points.col=3, plot.info=res, offset=1.1, shrink=1.2, lwd=3, stack=T)
      segments(0, 0, r[i]*cos(vectorX[i]), r[i]*sin(vectorX[i]), lwd=2, col= 'black')
      segments(0, 0, r[i]*cos(vectorY[i]), r[i]*sin(vectorY[i]), lwd=2, col= 'red')
      segments(0, 0, r[i]*cos(vectorZ[i]), r[i]*sin(vectorZ[i]), lwd=2, col= 'green')
      ani.pause()
    }
    }, 
    video.name = "C:/Users/Dinov/Desktop/KimeRandomSamplingAnimation.mp4", 
    other.opts = "-pix_fmt yuv420p -b:v 500K",
    ffmpeg = "E:/Ivo.dir/Ivo_Tools/ffmpeg-20200218-ebee808-win64-static/bin/ffmpeg.exe", 
    ani.width = 600, ani.height = 600
  )
  },
  htmlfile  = "Chapter6_Kime_Phases.html"
  
)
```

<video autosize:true controls>
  <source src="KimeRandomSamplingAnimation.mp4" type="video/mp4">
</video>

# Kime-Phases Helical (Spherical) Distributions 

We can use a spherical representation to visualize a more general case of infinitely-supported kime-phase distributions as helices. This approach does not require transforming the infinite distribution support to $[-\pi,\pi)$ and displays the native distribution as a helix embedded in $\mathbb{R}^3$.


```{r kimePhaseHelixDistribution, echo=T, warning=F, error=F, message=F}
# see these examples: 
#     https://plotly.com/r/3d-line-plots/
#     https://rpubs.com/daryl-hnl/plotly-3d
# parabola (as an example of a density model)
a <- 0; b <- 10 # Z-space vertical range of helix
per <- 2  # helical period
resol <- 100 # bin points density (resolution)
s <- seq(from=a, to=b, length.out=resol)
R <- 50 - 2*(s - (a+b)/2)^2  # parabola open down/inward # plot(s,R)
baseline <- 9.5  # radial offset

# A length, s, parameterized Helix equation ( Kime-phase density model)
# x = R*cos(s), y = R*sin(t), z = s
x <- R * cos(per*s) + baseline/2
y <- R * sin(per*s) + baseline/2
z <- s

# A length, s, parameterized Helix equation (baseline helix)
# x = cos(s), y = sin(t), z = s
xb <- 4*cos(per*s)
yb <- 4*sin(per*s)

# pair start & end indices i <--> i+50
pair_id <- rep(c(1:length(x)), times=2)
  
df <- data.frame(   # concat the x-start and x-end X coordinates
  X = c(x,xb), Y = c(y,yb), Z = c(z,z), R=R, phi=per*s,
  color_col = "lightgray",pair_id=pair_id)

# plot_ly() %>%
#   add_trace(x = df$X, y = df$Y, z = df$Z, color = df$color_col, 
#             type = "scatter3d", mode = "markers") %>%
#   add_trace(x = df$X, y = df$Y, z = df$Z, split = df$pair_id, 
#             line = list(color = "red"), type = "scatter3d", 
#             mode = "lines", showlegend = FALSE, inherit = FALSE)

plot_ly(x=x, y=y, z=z, type="scatter3d", line=list(width=10, color="navy blue"),
        mode='markers+lines', name="Phase Distribution") %>%
  # trace the main Z-axis
  add_trace(x=0, y=0, z=s, mode="lines",
              line = list(width = 3, color="navy blue"), name="core") %>%
              # , hoverinfo="none", showlegend=F) %>%
  # trace the baseline-helix
  add_trace(x=xb, y=yb, z=s, mode="markers+lines", 
            marker=list(size = 5, color="gray"),
            line = list(width = 5, color="gray"), name="Phase Domain") %>%
  # trace the (radial) density heights
  add_trace(x=df$X, y=df$Y, z=df$Z, split=df$pair_id, name="Phase Density", 
            line = list(color = df$color_col), type = "scatter3d",
            hovertext = paste("PhaseDensity(phi=", round(df$phi,3), 
                              ")=<br>  ", round(df$R, 3)),
            mode = "lines", showlegend = FALSE, inherit = FALSE) %>%
  layout(title="Helical Kime-Phase Density",
         scene = list(xaxis = list(title = 'k1'),
                      yaxis = list(title = 'k2'),
                      zaxis = list(title = 'Azimuthal Period')),
         legend = list(title="Functions", orientation = 'h'))

# test
# library(plotly)
# m_x = seq(-2, 2, .01)
# m_y = seq(-2, 2, .01)
# df = expand.grid(m_x, m_y)
# df['matrix'] = exp(-(df$Var1^2+df$Var2^2))
# m_z = matrix(df$matrix, nrow = length(m_x), ncol = length(m_y))
# m_df = list(m_x, m_y, m_z)
# 
# x1 = seq(-2, 0, by=0.0202)
# y1 = runif(100, min=-0.03,max=0.03)
# z1 = exp(-(x1^2+y1^2))
# df = data.frame(x1, y1, z1)
# names(df) <- c("df_x", "df_y", "df_z")
# 
# colors = c( "Blue",  "Cyan", "Green", "Yellow", "Orange", "Red")
# p1 <- plot_ly(x = m_df[[1]], y = m_df[[2]], z = m_df[[3]],
#               colors = colors, color = m_df[[3]]) %>% add_surface()  %>%
#   add_trace(x = ~df_x, y = ~df_y, z = ~df_z, data = df, 
#             type="scatter3d", mode="lines",
#             line = list(color = "blue", width = 4)) %>%
#   layout(title = "Hike_Example",
#        scene = list(aspectratio = list(x = 4, y = 4, z = 1)))
# p1
```

# Diffraction Pattern Interpretation of Kime-Phases

This section offers a *loose analogy* intended to provide *intuition* for the problem of the enigmatic kime-phase recovery.

The figure below illustrates one simple example of why the kime-phases may 
appear to be stochastic. This uncertainty is due to inabilities to visually disambiguate the depth perception
($3^{rd}$ dimension) from observing projections as proxy measurements tracking the underlying kime-phase density.

Look at the 3D scene from the top down, without using the mouse to change the position of the 3D scene. *Can you determine the actual kime-phase just from the angular appearance of the radial spokes (density heights)?* 

This is not possible, unless you navigate the 3D scene and utilize *zoom, pan*, and *rotation* positioning to associate each radial spoke to a specific spherical coordinate pair $(r, \theta, \phi)$. The reason is that the fixed *radius* plays the role of *time*, whereas the different *gray radial projections* correspond to *repeated samples* from the same controlled experiment under identical conditions (e.g., fixed spatiotemporal localization and fixed experimental conditions). In this example, we intentionally use random irregular angular-space sampling, $\phi\sim Uniform(a,b)$, which makes the polygonal density curve somewhat irregular, but also illustrates a more pragmatic model for temporal dynamics, $t_o\leq t\leq t_1$. Looking from the top, without a depth perception, we can't untangle the axial (transverse plane) containing the gray spoke (density value). Use the mouse to hover over the *gray density radial lines* and read off the corresponding phase values. These angular values should appear very irregular, perhaps not random from this small sample, but certainly oddly unusual. Lower and higher azimuthal planes correspond to phases, $\phi$, in the lower and higher spectrum of the domain, $\phi\in [0.1474376, 38.53021]$. However, in the top down view, they are all superimposed in a 2D flat projection of the 3D scene. This inability to untangle the *azimuthal axial plane* introduces an apparent single degree of freedom. In this case, the uncertainty can be easily resolved simply by exploring the actual 3D space. However, in the complex-time representation, that uncertainty manifests as random sample draws from kime-phase distributions, and is related to the Heisenberg uncertainty principle.

In practice, we can't directly observe the entire 5D spacekime universe. However, we can measure a number of observable process characteristics in the 4D Minkowski projection space.

```{r kimePhaseDistributionInterpret, echo=T, warning=F, error=F, message=F}
# this is a simple modification of the earlier example
# parabola (as an example of a density model)
a <- 0; b <- 10 # Z-space vertical range of helix
per <- 4  # helical period
resol <- 100 # bin points density (resolution)
s <- sort(runif(resol, min=a, max=b))
sb <- sort(runif(resol, min=a, max=b))
R <- 50 # - 2*(s - (a+b)/2)^2  # parabola open down/inward # plot(s,R)
Rb <- 4
baseline <- 1  # radial offset

# A length, s, parameterized Helix equation ( Kime-phase density model)
# x = R*cos(s), y = R*sin(t), z = s
x <- R * cos(per*s) + baseline/2
y <- R * sin(per*s) + baseline/2
z <- s

# A length, s, parameterized Helix equation (baseline helix)
# x = cos(s), y = sin(t), z = s
xb <- Rb*cos(per*sb)
yb <- Rb*sin(per*sb)

# pair start & end indices i <--> i+50
pair_id <- rep(c(1:length(x)), times=2)
  
df <- data.frame(   # concat the x-start and x-end X coordinates
  X = c(x,xb), Y = c(y,yb), Z = c(z,z), R=c(rep(R,resol),rep(Rb, resol)),
  phi=c(per*s,per*sb), color_col="lightgray", pair_id=pair_id)

plot_ly(x=x, y=y, z=z, type="scatter3d", line=list(width=10, color="navy blue"),
        mode='markers+lines', name="Phase Distribution") %>%
  # trace the main Z-axis
  add_trace(x=0, y=0, z=s, mode="lines",
              line = list(width = 3, color="navy blue"), name="core") %>%
              # , hoverinfo="none", showlegend=F) %>%
  # trace the baseline-helix
  add_trace(x=xb, y=yb, z=sb, mode="markers+lines", 
            marker=list(size = 5, color="gray"),
            line = list(width = 5, color="gray"), name="Phase Domain") %>%
  # trace the (radial) density heights
  add_trace(x=df$X, y=df$Y, z=df$Z, split=df$pair_id, name="Phase Density", 
            line = list(color = df$color_col), type = "scatter3d",
            hovertext = paste("phi=", round(df$phi,3), 
                              "<br>AxialPlane=", round(df$phi/(2*pi),0), 
                              "<br>PhaseDensity=", round(df$R, 3)),
            mode = "lines", showlegend = FALSE, inherit = FALSE) %>%
  layout(title="Diffraction Pattern Interpretation of Kime-Phases",
         scene = list(xaxis = list(title = 'k1'),
                      yaxis = list(title = 'k2'),
                      zaxis = list(title = 'Azimuthal Period')),
         legend = list(title="Functions", orientation = 'h'))
```


# Spherical Electron Orbit Densities 

In 4D Minkowski spacetime, the electron orbitals are temporally-dynamic covering the 3D spatial domain. In 3D Cartesian coordinates, the position/coordinates of points are described by $\pmb{r}=(x,y,z)\in \mathbb{R}^3$, whereas in spherical coordinates, the same points are described by the displacement magnitude (of the vector, $r$, from the coordinate space origin), and a pair of angles $(\phi, \theta)$ between the projection of the vector $r$ into the $(x,y)$ plane and the positive $x$-axis (*azimuthal angle* $\phi$), and between $r$ and the positive $z$-axis (*polar angle* $\theta$). Explicitly, the spherical coordinate representation of spatial locations involves $(r, \theta, \phi)\in \mathbb{R}^3$, where $r\in[0,\infty),\ \theta \in \left [-\frac{\pi}{2},\frac{\pi}{2}\right ], \ \phi\in [-\pi,\pi)$.

The bijective forward and inverse (*Descarte*) *Cartesian* to *spherical* coordinate transformations are: 

 - Spherical to Cartesian transform

$$x=r\sin\theta \cos\phi\\
y=r \sin\theta \sin\phi \\
z=r \cos\theta \ .$$

 - Cartesian to spherical transform

$$r=\sqrt{x^2+y^2+z^2} \\
\phi = \tan^{-1}\frac{y}{x} \\
\theta=\tan^{-1}\frac{\sqrt{x^2+y^2}}{z} \ .$$

## Angular momentum

In spherical coordinates, the momentum $p$ of a particle, such as an electron, has two orthogonal (independent) components; a *radial component* $p_r$ corresponding to the particle motion radially outward from the origin, and an *angular component* $L$, corresponding to rotational particle motion along the surface of a sphere of radius $r$ centered at the origin, e.g., atomic nucleus.
Hence, the particle momentum is $p=\left (p_r,\frac{L}{r}\right )$.
Since, the angular component $L$ is a vector of size $r$, $\frac{L}{r}$ is a vector in the same direction of length $1$, both tracking the particle motion along the surface of a sphere, according to the azimuthal ($\theta$) and polar ($\phi$) angles.

In Cartesian coordinates, the *angular momentum vector* is the [vector cross product](https://en.wikipedia.org/wiki/Cross_product)
$$L=\vec{r}\times p ,\ {\text{note that }} L\perp \vec{r} \ ,$$
and it's $x,y,z$ components are given in terms of the Cartesian coordinates of $r$ and $p$,

$$L_x=yp_z-zp_y\\
L_y=zp_x-xp_z\\
L_z=xp_y-yp_x \ .$$

For simplicity, with a slight abuse of notation, sometimes we may denote *scalar* and *vector* quantities with the same symbol, e.g., the vector $\vec{r}$ and it's size $|\vec{r}|=r$ may both be written as $r$.

When the potential energy $V$ depends only on $r$, the angular momentum is conserved, i.e., it's constant and has a trivial temporal derivative. Consider for instance
the [Coulomb potential](https://en.wikipedia.org/wiki/Coulomb%27s_law) $V(r)=-\frac{e^2}{4\pi \epsilon_o}\frac{1}{r}$.

Recall that the [Coulomb’s law](https://en.wikipedia.org/wiki/Coulomb%27s_law) implies that the electric field $E$ at point $r=r(x,y,z)$, due to a charge $q$ located at the origin,
is given by 
$E = k \frac{q}{r^2}r$. To show this angular momentum conservation, we start with the force 
acting on a classical electron $F=-\frac{k}{r^3}r$, and it's Cartesian components
$F_x=-\frac{k}{r^3}x$, $F_y=-\frac{k}{r^3}y$, and $F_z=-\frac{k}{r^3}z$.

[Newton's second law of motion](https://en.wikipedia.org/wiki/Newton%27s_laws_of_motion) is

  $$F=ma=m\frac{dv}{dt}=\frac{d(mv)}{dt}=\frac{dp}{dt}.$$

Let's compute the time-derivative of the first component of the angular momentum
vector,

$$\frac{dL_x}{dt}\equiv \frac{d (yp_z-zp_y)}{dt} =\\
\frac{dy}{dt}p_z + y\frac{dp_z}{dt}-\frac{dz}{dt}p_y-z\frac{dp_y}{dt}=
\frac{p_y}{m}p_z + yF_z- \frac{p_z}{m}p_y -zF_y=\\
\frac{p_y p_z}{m}-\frac{kyz}{r^3} - \frac{p_z p_y}{m}+ \frac{kzy}{r^3}=0 .$$

Similarly, $\frac{dL_y}{dt}=\frac{dL_z}{dt}=0$, and $L^2=L^2_x+L^2_y+L^2_z$, hence $\frac{dL}{dt}=0$, and therefore, 
the *energy* and the *angular momentum* are conserved. 
In quantum mechanics, the energy can have only certain allowed values. 
For instance, the [Bohr model](https://en.wikipedia.org/wiki/Bohr_model) predicts
that all possible energies for the hydrogen atom (electron layers $n=1,2,3,\cdots$) are in the form

$$E_n =- \frac{Z^2 e^4 m_e}{8\epsilon_o^2 h^2}\frac{1}{n^2} \ .$$
This Cartesian coordinate representation of the angular momentum $L=L(\theta,\phi)$ has an
analogous *spherical coordinate* representation in terms of the momentum 
$p=\left (p_r,\frac{L}{r}\right )$. In spherical coordinates, since $L=r\times p$, the classical energy is given by

$$E=\frac{p^2_r}{2m_e} + \frac{L^2}{2m_er^2} - \frac{Ze^2}{4\pi \epsilon_o r} =
\underbrace{\frac{p^2_r}{2m_e} - \frac{Ze^2}{4\pi \epsilon_o r}}_{\epsilon_r,\ radial\ energy} +
\underbrace{\frac{L^2}{2m_er^2}}_{angular\ energy} .$$

The energy $E$ depends on the nucleus of charge ($+Ze$) and $L^2$, which is conserved along with its three Cartesian components
$L_x,L_y,L_z$. In general, $L$ has three Cartesian components, however, in really is just a 2D parametric vector modeling the motion of an electron on a 2D manifold, a 2D sphere embedded in 
$\mathbb{R}^3$, parameterized by the angles $\theta,\phi$.

Cartesian coordinates in $\mathbb{R}^3$ require 3 components, but the 
sphere coordinates on $S^2$ have only 2 degrees of freedom, *polar* ($\theta$) and *azimuthal* ($\phi$). 
Hence, we can parameterize the *angular momentum* ($L$) and the *energy* ($E$)
using two parameters, e.g., the *magnitude* $L^2$ and the *angular momentum* 
$z$-component ($L_z$), which describes the particle motion projected in the $(x,y)$ 
plane parameterizing the motion in the azimuthal ($\phi$) direction.
As a result, the wavefunction describing the motion of the particle
will have a ring component 
$\frac{1}{\sqrt{2\pi}}e^{im\phi}, \ \forall\ m\in\mathbb{Z}$.
This limits the possibilities for the second energy component, the $z$-component,
only to $L_z=m\hbar$.

**Note**: [This video presents the basics of the raising and lowering eigenstates ](https://youtu.be/GLPSxU9I-Ao?list=PL65jGfVh1ilueHVVsuCxNXoxrLI3OZAPI&t=1226) for the [quantum ladder operator](https://en.wikipedia.org/wiki/Ladder_operator).

The  allowed values for the first component of the energy, the magnitude $L^2$,
are also discrete (energy levels) $L^2\to l(l+1)\hbar, \ l\in\{0, 1,2, \cdots, (n-1)\}$
for each energy level characterized by $n$. 
Notice that since $L_z\leq |L|$, $\forall\ l$, $m\leq l$ implies that
$m\in \{-l,-l+1,-l+2,\cdots, 0, \cdots,l-2,l-1,l\}$.

## Schrödinger Equation for a Hydrogen atom

In 1D, the *time-dependent* [Schrödinger equation](https://en.wikipedia.org/wiki/Schr%C3%B6dinger_equation)

$$i\hbar \frac {\partial }{\partial t}\Psi (x,t)=
\left[-\frac {\hbar ^{2}}{2m}\frac {\partial ^{2}}
{\partial x^{2}} + V(x,t)\right ]\Psi (x,t),$$

where $\Psi(x=space,t=time):\mathbb{R}\times \mathbb{R}^+ \to\mathbb{C}$ is a 
wavefunction,
$m$ is the mass of the particle, $V(x,t)$ is the potential that represents
the environment the particle is embedded in, and 
$$\hbar=\frac{h}{2\pi}=1.054571817\cdots\times 10^{-34} J\cdot s\equiv 	6.582119569\cdots\times 10^{-16} eV\cdot s$$ 
is the [reduced Planck 
constant](https://en.wikipedia.org/wiki/Planck_constant).

The above *time-dependent* Schrödinger equation illustrates that wavefunctions can form standing waves, *stationary states* that simplify the general solutions of the time-dependent Schrödinger equation for any state. 

Stationary states can also be described by the *time-independent* Schrödinger equation (TISE)

$${\hat {H}} |\Psi \rangle = E |\Psi \rangle ,$$
an *eigenvalue equation*, where $E$ is the energy of the system and the Hamiltonian operator ${\hat {H}}$ is time-independent. The total wavefunction $\Psi$ is still *time-dependent*, even though the TISE equation is time-agnostic, yet the eigen-solutions (eigenvalue=Energy, eigenfunction=wavefunction) retain their time-dependence. 

Let's consider a single-electron Coulomb system and express
the [Schrödinger equation](https://en.wikipedia.org/wiki/Schr%C3%B6dinger_equation)
in spherical coordinates, where we employed the [Laplacian to spherical coordinate conversion](https://mathworld.wolfram.com/LaplacesEquationSphericalCoordinates.html)

$$E |\psi \rangle \equiv E\psi(r,\theta,\phi)=\\
\underbrace{-\frac{\hbar^2}{2m_e r^2}
\left [\frac{\partial }{\partial r}r^2\frac{\partial }{\partial r}(\psi(r,\theta,\phi))
+ \frac{1}{\sin\theta} \frac{\partial}{\partial \theta} \sin\theta 
\frac{\partial}{\partial \theta} \psi (r,\theta,\phi) +
\frac{1}{\sin^2\theta} \frac{\partial^2}{\partial \phi^2} 
\psi(r,\theta ,\phi )\right]}_{angular\ energy} -
\underbrace{\frac{Ze^2}{4\pi\epsilon_o r} \psi(r,\theta,\phi)}_{\epsilon_r,\ radial\ energy}
\equiv\\ {\hat {H}} |\psi \rangle\ .$$

Solutions to this PDE yield the allowed values of the *angular momentum*, $L$, 
and the allowed *energies*, $E_n$.

The formulation of the Schrödinger equation is based on the classical energy

$$E=\frac{p^2_r}{2m_e} + \frac{L^2}{2m_er^2} - \frac{Ze^2}{4\pi \epsilon_o r} =
\underbrace{\frac{p^2_r}{2m_e} - 
\frac{Ze^2}{4\pi \epsilon_o r}}_{\epsilon_r,\ radial\ energy\ \perp (\theta,\phi)} +
\underbrace{\frac{L^2}{2m_er^2}}_{\epsilon_{\theta,\phi}\ angular\ energy,\ \perp r} .$$

In statistical data analysis, this independent component bifurcation of the energy$E$ resembles the classical
[variance decomposition of ANOVA](https://en.wikipedia.org/wiki/Analysis_of_variance#Algorithm)
into *variance between treatments* and an independent part,
*variance within treatments*, and the [bias and precision tradeoff](https://en.wikipedia.org/wiki/Bias%E2%80%93variance_tradeoff) in statistical learning.

This energy separability into *radial* and *angular* terms facilitates the decomposition of the wavefunction into a *product* of two terms

$$\psi(r,\theta,\phi)=
\underbrace{R(r)}_{radial} \cdot 
\underbrace{Y(\theta,\phi)}_{angular}\ , $$

where $r$ is the radius (in atomic units, $1$ Bohr radius 
$= 52.9\ pm = 52.9\times 10^{-12}$) from the atomic center, $0\leq \theta\leq\pi$
is the polar spherical coordinate, and $-\pi\leq \phi<\pi$ is the azimuthal spherical coordinate. 

The solutions $Y_{lm}(\theta,\phi),\ l,m\in\mathbb{Z}$ of the angular part,
$Y(\theta,\phi)$, are called *spherical harmonics* and yield the 
possible values of the *angular momentum* $L^2$ and the $z$-component $L_z$. 

Note that the [Euler formula](https://en.wikipedia.org/wiki/Euler%27s_formula) explicates the complex exponentials 

$$e^{i\phi}=\cos\phi+i\sin\phi \quad \& \quad e^{-i\phi}=\cos\phi-i\sin\phi. $$

Let's explore a few of the spherical harmonics for some $l,m$ parameters.

 - For $l=0$, there is just one possibility for $m$, $m=0$ corresponding to a unique spherical harmonic term (a constant)

$$Y_{00}(\theta,\phi)=\frac{1}{\sqrt{4\pi}}. $$

 - For $l=1$, there are 3 possible spherical harmonics corresponding to $m=-1,0,1$, 

$$Y_{11}(\theta,\phi)= -\sqrt{\frac{3}{8\pi}} \sin \theta\ e^{i\phi}. $$

$$Y_{1\ -1}(\theta,\phi)= -\sqrt{\frac{3}{8\pi}} \sin \theta\ e^{-i\phi}. $$

$$Y_{1\ 0}(\theta,\phi)= -\sqrt{\frac{3}{4\pi}} \cos \theta. $$

 - For $l=2$, there are 5 possibilities for $m$, $m\in\{-2,-1,0,1,2\}$ corresponding to 5 spherical harmonics,

$$Y_{2\ 2}(\theta,\phi)=\sqrt{ \frac{15}{32\pi}}\sin^2 \theta\ e^{2i\phi},$$
$$Y_{2\ -2}(\theta,\phi)=\sqrt{\frac{15}{32\pi}}\sin^2 \theta\ e^{-2i\phi},$$
$$Y_{2\ 1}(\theta,\phi)=-\sqrt{ \frac{15}{8\pi}}\sin \theta \cos \theta\ e^{i\phi},$$
$$Y_{2\ -1}(\theta,\phi)=\sqrt{ \frac{15}{8\pi}}\sin \theta \cos \theta\ e^{-i\phi},$$
$$Y_{2\ 0}(\theta,\phi)=\sqrt{ \frac{5}{16\pi}} (3\cos^2 \theta - 1).$$

Recall that 
$\psi(r,\theta,\phi)=\underbrace{R(r)}_{radial} \cdot \underbrace{Y(\theta,\phi)}_{angular}$ 
and the first factor in the Schrödinger equation, the *radial component* $R_{n\ l}(r)$,
depends on the parameters $n,l\in\mathbb{Z}$, whereas the second term, the *angular component*, $Y_{l\ m}(\theta,\phi)$, depends on $l,m\in\mathbb{Z}$.

$$E_n\ R_{nl}(r)=
-\frac{\hbar^2}{2m_e}\frac{1}{r^2}
\frac{d }{d r}r^2\frac{d }{d r}(R_{n\ l}(r))
+ \frac{l(l+1)\hbar^2}{2m_e r^2} R_{n\ l} (r)
- \frac{Z e^2}{4\pi\epsilon_o r} R_{n\ l} (r) .$$

The *angular functions* $Y_{l\ m}(\theta,\phi)$ are independent of the potential 
$V(r)$, whereas the *radial functions* $R_{n\ l}(r)$ directly
depend on the Coulomb potential. Hence, the solutions to the radial Schrödinger
equation, $R_{n\ l}(r)$, control the allowed energy levels. 
The boundary conditions that lead to the quantized energies are 
$\lim_{r\to 0}(rR_{n\ l}(r))=0$ and $\lim_{r\to \infty}(rR_{n\ l}(r))=0$. 

See the [University of Sheffields *Orbitron* website by Mark Winter](https://winter.group.shef.ac.uk/orbitron/atomic_orbitals/3d/3d_equations.html), which includes interactive content, images, and renditions representing atomic orbitals, radial and angular wavefunctions, animated plots of electron density and distribution functions.

The atomic orbitals are indexed by the quantum numbers and include 
$\{Level_1, Level_2, \cdots, Level_7\}$ and higher order orbitals
$\{g, h, i\}$ orbitals, hybrid orbitals, and molecular orbitals.

Let's denote the constant
$$a_o=\frac{4\pi\epsilon_o\hbar^2}{e^2m_e} =0.529177×10^{-10}, $$
and determine several radial functions part of the wavefunction solutions 
corresponding to the first few values of $n,l$

$$R_{1\ 0}(r)=2\left (\frac{Z}{a_o} \right)^{3/2} e^{\frac{-Zr}{a_o}},$$
$$R_{2\ 0}(r)=\frac{1}{2\sqrt{2}}\left (\frac{Z}{a_o} \right)^{3/2}
\left (2- \frac{Zr}{a_o}\right )e^{\frac{-Zr}{2a_o}},$$

$$R_{2\ 1}(r)=\frac{1}{2\sqrt{6}}\left (\frac{Z}{a_o} \right)^{3/2}\frac{Zr}{a_o} e^{\frac{-Zr}{2a_o}},$$

$$R_{3\ 0}(r)=
\frac{2}{81\sqrt{3}}\left (\frac{Z}{a_o} \right)^{3/2}
\left (27- 18\frac{Zr}{a_o} +2\left (\frac{Zr}{a_o} \right )^2\right )
e^{-\frac{Zr}{3a_o}}.$$

$$\underbrace{R_{3\ 2}(r)}_{3d_{z^2}}=
\frac{4}{81\sqrt{30}} \left (\frac{Z}{a_o}\right )^{3/2} 
\left (\frac{Zr}{a_o} \right)^{2} e^{-\frac{Zr}{3a_o}}.$$

Reconstructing the complete wavefunction solutions from the corresponding
radial and angular components yields

$$\psi_{n\ l\ m}(r,\theta,\phi)=\underbrace{R_{n\ l}(r)}_{radial} \cdot 
\underbrace{Y_{l\ m}(\theta,\phi)}_{angular}, \quad n,l,m\in\mathbb{Z} .$$

The 3 parameters $n,l,m\in\mathbb{Z}$ characterizing each state are called
the *quantum numbers* of the system and each parameter corresponds to a 
conserved quantity 

 - $n$ is the *principal* quantum number, 
 - $l$ is the *angular* quantum number, and 
 - $m$ is the *magnetic* quantum number.

The quantum interpretation of the wavefunction solution
$\psi_{n\ l\ m}(r,\theta,\phi)$ is as probability density amplitude
for finding the electron in a particular region of the 3D space.
The wavefunciton amplitudes are used to compute actual probabilities 
associated with measurements of the electron's position just like a CDF
can be computed via integration of the probability density function over
a specified domain. The likelihood of the electron being in a small voxel 
(volume element) $dV\subseteq \mathbb{R}^3$ around the point 
${\bf{r}}=(r,\theta,\phi)\in dV$ is

$$P(\text{electron in } dV\text{ about location } (r,\theta,\phi))\equiv
|\psi_{n\ l\ m}(r,\theta,\phi)|^2\ dV=
R^2_{n\ l}(r)|Y_{l\ m}(\theta,\phi)|^2\ dV. $$

In *Cartesian coordinates*, the small 3D volume element $dV$
represents a small box of dimensions $dx\cdot dy\cdot dz$, i.e.,
$dV=dx\ dy\ dz$.

In spherical coordinates, the volume element $dV$ corresponds 
the spherical volume of a ball around the point of interest

$$dV=r^2\sin\theta\ dr\ d\theta\ d\phi.$$

This can be verified by integrating the volume element $dV$ over a 3D ball and 
confirming the ball volume. Parameterize the 3D ball, $dV$, by 

$$Domain\equiv\left[0\leq \underbrace{r}_{radius}\leq R\right ] \times
\left[0\leq \underbrace{\theta}_{azimuthal}\leq \pi\right ] \times
\left[-\pi\leq \underbrace{\phi}_{angular}\leq \pi\right ].$$

By separability,
$$V_{ball}=\int_{Domain}{dV}= \int_{Domain} {r^2\sin\theta dr\ d\theta\ d\phi}=
\int_{0}^r {\int_{0}^{\pi} {\int_{-\pi}^{\pi} r^2\sin\theta dr\ d\theta\ d\phi}}=\\
\int_{0}^r { r^2\int_{0}^{\pi} { \sin\theta\int_{-\pi}^{\pi} dr\ d\theta\ d\phi}}=\\
\left (\frac{r^3}{3}\left |_{0}^R \right . \right )\times
\left (-\cos \theta \left |_{0}^{\pi} \right . \right ) \times 
\left (\phi \left |_{-\pi}^{\pi} \right .  \right )=
\frac{r^3}{3}\times 2\times 2\pi=\frac{4}{3}\pi R^3\ .$$


### Example 1: $n=2, l=1, m=0$ (Trivial)

Remember the constraints on the *quantum troika* $l=0,1,2,\cdots,n-1$ and 
$m\in\{-l,-l+1,\cdots,0,\cdots,l-1,l\}$. Consider the electron 
of a hydrogen atom $(Z=1)$ in a state with quantum numbers $n=2, l=1, m=0$. 
We'll compute the probability that the electron's position is farther from the nucleus than $a_o$, i.e., $Pr(r\geq a_o)$.

Start with the general wavefunction formulation
$$\psi_{n\ l\ m}(r,\theta,\phi)=\underbrace{R_{n\ l}(r)}_{radial} \cdot 
\underbrace{Y_{l\ m}(\theta,\phi)}_{angular}, \quad n,l,m\in\mathbb{Z} .$$

For the specific quantum triple $(n=2, l=1, m=0)$ we have
$$\psi_{n=2\ l=1\ m=0}(r,\theta,\phi) = \psi_{2\ 1\ 0}(r,\theta,\phi) =
\underbrace{R_{2\ 1}(r)}_{radial} \cdot 
\underbrace{Y_{1\ 0}(\theta,\phi)}_{angular}=\\
\frac{1}{2\sqrt{6}}\left (\frac{Z}{a_o} \right)^{3/2}\frac{Zr}{a_o} e^{\frac{-Zr}{2a_o}}
\times \left (-\sqrt{\frac{3}{4\pi}} \cos \theta \right ) 
\underbrace{=}_{Z=1}\\
-\frac{1}{4\sqrt{2\pi}} \frac{r}{a_o^{5/2}} e^{\frac{-r}{2a_o}} \cos \theta .$$

Integrating the amplitude of the wave function 
$|\psi_{2\ 1\ 0}(r,\theta,\phi)|^2$ against the volume element, $dV$
yields the probability $Pr(r\geq a_o)$, reflecting electron positioned
farther from the nucleus than $a_o$.

Next, derive the probability $Pr(r\geq a_o)$ using the wavefunction

$$\psi_{n=2\ l=1\ m=0}(r,\theta,\phi)=R_{n=2\ l=1}(r) \cdot Y_{l=1\ m=0}(\theta,\phi)=\\
\underbrace{\frac{1}{2\sqrt{6}}\left (\frac{Z}{a_o} \right)^{3/2}\frac{Zr}{a_o} e^{\frac{-Zr}{2a_o}}}_{R_{2\ 1}(r)} \cdot
\underbrace{\left (-\sqrt{\frac{3}{4\pi}} \cos \theta\right )}_{Y_{1\ 0}(\theta,\phi)}=\\
-\frac{e^{-\frac{r Z}{2 a_o}} r Z^2 \sqrt{\frac{Z}{a_o}} \cos(\theta)}{4 a_o^2 \sqrt{2\pi}} .$$

$$Pr(r\geq a_o)\equiv \int_{a_o}^{\infty} { \int_{0}^{\pi} { 
\int_{-\pi}^{\pi} |\psi_{2\ 1\ 0}(r,\theta,\phi)|^2\ r^3\sin \theta\ dr\ d\theta\ d\phi}}=\\
\int_{a_o}^{\infty} { \int_{0}^{\pi} { 
\int_{-\pi}^{\pi} \left(-\frac{e^{-\frac{r Z}{2 a}} r Z^2 \sqrt{\frac{Z}{a}} \cos(\theta)}{4 a^2 \sqrt{2\pi}}\right )^2 \ r^3\sin \theta\ dr\ d\theta\ d\phi}} =\\
\int_{a_o}^{\infty} { \int_{0}^{\pi} { 
\int_{-\pi}^{\pi} \left (\frac{rZ}{a_o}\right )^5 \frac{e^{-\frac{r Z}{a_o}} \cos^2\theta\ \sin\theta}{32 \pi}\ dr\ d\theta\ d\phi}} =
\cdots =\frac{1}{24} e^{-Z} (24 + Z (24 + Z (12 + Z (4 + Z)))), \\
\text{for}\ Re\left(\frac{Z}{a_o}\right )>0.$$

For $Z=1$,
$$Pr(r\geq a_o)=
    \begin{cases}
        \frac{65}{24\ e}\equiv 0.9963402 & \text{, if } a_o< 1\\
        1 & \text{, if } a_o\geq 1
    \end{cases}
.$$


### Example 2: $n=1, l=0, m=0$ (Non-Trivial)

Again, recall the constraints on the *quantum triple* $\forall\ n\ ,\ l=0,1,2,\cdots,n-1$ and 
$m\in\{-l,-l+1,\cdots,0,\cdots,l-1,l\}$. This is a non-trivial example for the electron 
of a hydrogen atom $(Z=1)$ in a state with quantum harmonic numbers $(n=1, l=0, m=0)$. 
Let's compute the same the probability that in this case, the electron's position is farther from the nucleus
than $a_o$, i.e., $Pr(r\geq a_o)$.

$$\psi_{n=1\ l=0\ m=0}(r,\theta,\phi)=\underbrace{R_{1\ 0}(r)}_{radial} \cdot 
\underbrace{Y_{0\ 0}(\theta,\phi)}_{angular} =$$


$$\underbrace{2\left (\frac{Z}{a_o} \right)^{3/2} e^{-\frac{Zr}{a_o}}}_
{(radial)\ R_{1\ 0}(r)} \times 
\underbrace{\left (\frac{1}{\sqrt{4\pi}}  \right )}_
{(angular)\ Y_{0\ 0}(\theta,\phi)}=
\frac{e^{-\frac{rZ}{a_o}}Z \sqrt{\frac{Z}{a_o}}}{a_o \sqrt{\pi}} .$$

Integrating the amplitude of the wave function 
$|\psi_{1\ 0\ 0}(r,\theta,\phi)|^2$ against the volume element, $dV$
yields the probability $Pr(r\geq a_o)$, reflecting electron positioned
farther from the nucleus than $a_o$.

$$Pr(r\geq a_o)\equiv \int_{a_o}^{\infty} { \int_{0}^{\pi} { 
\int_{-\pi}^{\pi} |\psi_{1\ 0\ 0}(r,\theta,\phi)|\ r^2\sin \theta\ dr\ d\theta\ d\phi}}=\\
\int_{a_o}^{\infty} { \int_{0}^{\pi} { 
\int_{-\pi}^{\pi} \left(\frac{e^{-\frac{-rZ}{a_o}}Z \sqrt{\frac{Z}{a_o}}}{a_o \sqrt{\pi}}
\right )^2 \ r^2\sin \theta\ dr\ d\theta\ d\phi}} =\\
\int_{a_o}^{\infty} { \int_{0}^{\pi} { 
\int_{-\pi}^{\pi} 
\frac{e^{-\frac{2r Z}{a}} r^2 Z^3 \sin\theta}{a^3 \pi}\ dr\ d\theta\ d\phi}}=
e^{2Z}(1+2Z(1+Z)), \ \text{for}\ Re\left(\frac{Z}{a_o}\right )>0.$$


For $Z=1$,
$$Pr(r\geq a_o)=
    \begin{cases}
        \frac{5}{e^2}\approx 0.6766764  & \text{if } a_o< 1\\
        1 & \text{if } a_o\geq 1
    \end{cases}
.$$

**Hence, in the case of the quantum triple $(n=1, l=0, m=0)$, given that $a_o< 1$, we have a relatively large probability, about $\frac{2}{3}$, that the 
particle is at least $1$ Bohr radius away from the hydrogen nucleus!**

Previously, in the case of a different quantum triple $n=2, l=1, m=0$, we had computed the probability $Pr(r\geq a_o)$ to be $0.9963402\approx \frac{3}{3}$, which is much larger than that for this quantum triple $(n=1, l=0, m=0)$.


## Radial Probability Distribution Functions

The following product, part of the above probability, is the 
*radial probability distribution function* 

$$P_{n\ l}(r)=r^2R^2_{n\ l}(r) . $$

We interpret $P_{n\ l}(r)dr$ as the *probability that the 
electron's position is in a radial shell of thickness $dr$ around a 2-sphere of 
radius $r$*.


*Notes*:

 - The radial probability distribution quantifies why the electron cannot 
 collapse into the nucleus, since $P_{n\ l}(r=0)\equiv 0$ in a measure theoretic sense corresponds to the likelihood of such collapse must be zero. Also, shrinking the radial shell rapidly
 reduces the probability of finding the electron in that shell, $\lim_{r\to 0}P_{n\ l}(r)= 0$.

 - Since each wavefunction corresponds to a probability distribution quantifying
 the likelihood an electron can be found for each energy, the number of allowed states corresponds to the electron properties and system behavior. For instance, 
 
   + For $n=1$, there is *only 1 energy* $E_1$ and 1 wavefunction $\psi_{1\ 0\ 0}$.
   + For $n=2$, there is *only 1 energy* $E_2$ with 4 possible states corresponding
   to the following possible harmonic triples $(n,l,m)$ and the four wavefunctions
   $\{\psi_{2\ 0\ 0}; \psi_{2\ 1\ -1}; \psi_{2\ 1\ 0}; \psi_{2\ 1\ 1} \}$. 
   Recall that $l=0,1,2,\cdots,n-1$ and $m\in\{-l,-l+1,\cdots,0,\cdots,l-1,l\}$.

$$l=0\quad, \quad m=0\\
l=1\quad, \quad m\in \{-1,0,1\} .$$

 - When more than $1$ wavefunctions correspond to a given energy level $E_n$, 
 that energy level is called *degenerate*. For the hydrogen atom in this example,
 the energy level, $E_2$ is 4-fold degenerate.


### Sharp ($s1$) $l=0$ orbitals

The wavefunctions $\psi_{n\ l\ m}(r,\theta,\phi)$ of the electron are called *orbitals*,
not because they represent the static 2D surface shells of orbits, but because they
quantify the probabilities of the electron occupying certain volume for a specific
allowable energy. The surface shapes of the orbitals are determined by the 
values of *angular momentum*. Let's examine the orbital shapes corresponding 
to the $l=0$ orbitals, called $s$ (*sharp*) orbitals. When $l=0 \implies m=0$
and the corresponding wavefunctions are 

$$\psi_{n\ l=0\ m=0}(r,\theta,\phi)=\underbrace{R_{n\ 0}(r)}_{radial} \cdot 
\underbrace{Y_{0\ 0}(\theta,\phi)}_{angular} =
R_{n\ 0}(r)\cdot \left (\frac{1}{4\pi}\right ) .$$

As $Y_{0\ 0}$ is a constant, 
$\psi_{n\ l=0\ m=0}(r,\theta,\phi)\equiv \psi_{n\ l=0\ m=0}(r)$, 
i.e., $\psi_{n\ l=0\ m=0}$ is independent of $(\theta, \phi)$. 
Therefore, all of these 4 orbitals must be *spherically symmetric*.


In particular, the electron probability density in any plane through the nucleus
(origin) would appear as concentric circles. Different orbital boundaries, 
also called *nodes*, correspond to trivial ($0$) probability density values. The contour surface boundaries enclose the largest proportions of places occupied by
electrons ($\sim 90\%$ of the electron probability), at different $1s$, $2s$, and $3s$
orbitals. Of course, the wave function may be either positive or negative, 
but its *amplitude* is always non-negative, $|\psi_{n\ l=0\ m=0}(r,\theta,\phi)|^2\geq 0$.

The density of points decreases exponentially as $r$ increases, 
because of the exponential dependence of the functions $R_{n\ 0}(r)$.

Also, $n\to\infty$ implies exponentials decay more slowly,
$e^{-r/a_0}$ for $n=1$, $e^{-r/2a_0}$, for $n=2$, and 
$e^{-r/3a_0}$, for $n=3$. The wavefunctions peak at $r=0$, suggesting
that the maximum amplitude to find the electron is on top of the nucleus (at the origin)!
However, since the *radial probability* density contains an extra $r^2$ factor
corresponding to the volume element $dV$, the probability 
 $\lim_{r\to 0}P_{n\ l}(r) =0$.

The most probable annulus for the hydrogen electron in the ground state. 
The following plot shows the dependence of 3 functions on the radius, $r$, 

 - The electron *probability density* $\psi^2$ smoothly decreasing with $r\uparrow$.  The density is greatest in the innermost shells of the state space, near the nucleus.
 - The *surface area* (spherical for hydrogen) of each shell, $S=4\pi r^2$, monotonically increases with $r\uparrow$. 
 - The distribution of dots in each spherical shell, reflecting the *total probability*  of finding the electron at a given value of $r$. With $r\uparrow$, the surface area of each shell increases faster than the electron probability density decreases. The electron *radial probability* has a peak corresponding to the  *(expected) most probable radius* for the $1s$ electron, $52.9\ pm$, the radius predicted by Bohr’s model of the hydrogen atom.


```{r electronDistributions, echo=T, warning=F, error=F, message=F}
# Bohr radius constant
a_o <- 0.529177 * (10^(-10))

# radial parts of the wavefunctions using pm picometer (1×10-¹² of a meter) distance
R_1_0 <- function (r, a=a_o*(10^(10)), Z=1) {
  return (2 * ((Z/a))^(3/2) * exp(-Z*r/a))
}

R_2_0 <- function (r, a=a_o*(10^(10)), Z=1) {
  return ( (1/(2*sqrt(2))) * (Z/a)^(3/2) * (2- (Z*r)/(a)) * exp(-Z*r/(2*a)) )
}

R_2_1 <- function (r, a=a_o*(10^(10)), Z=1) {
  return ( (1/(2*sqrt(6))) * (Z/a)^(3/2) * (Z*r/a) * exp(-Z*r/(2*a)))
}

R_3_0 <- function (r, a=a_o*(10^(10)), Z=1) {
  return ( (2/(81*sqrt(3))) * (Z/a)^(3/2) * (27- 18*(Z*r/a) +
                                                2*(Z*r/a)^2) * exp(-Z*r/(3*a)))
}

pdf_psi2 <- function (r, n, l) {
  #####  return ( R(r, n, l) * (1/(4*pi)) )
  if (n==1 & l==0) return ( (R_1_0(r) * (1/(4*pi)))^2 )
  else if (n==2 & l==0) return ( (R_2_0(r) * (1/(4*pi)))^2 )
  else if (n==2 & l==1) return ( (R_2_1(r) * (1/(4*pi)))^2 )
  else if (n==3 & l==0) return ( (R_3_0(r) * (1/(4*pi)))^2 )
  else { print('Invalid harmonic parameters (n,l) ...')
    return()
  }
} # Test: pdf_psi2(1, 3,0)

sphericalSurfaceArea <- function (r) {
  return (4*pi*r^2)
} # Test: sphericalSurfaceArea(1)

radialProbability <- function (r, n, l) {
  return (pdf_psi2(r, n, l) * r^2)
} # Test: radialProbability(1, 1, 0)

rad <- seq (from=10^(-5), to=2, length.out=1000)

df <- data.frame(cbind(radius=rad, probDensity=pdf_psi2(rad, 1,0), 
                       sphericalSurfaceArea=sphericalSurfaceArea(rad),
                       radialProbability=radialProbability(rad, 1,0)))
library(plotly)

df %>% plot_ly() %>%
  add_trace(x=~rad, y=~probDensity, type="scatter", mode="lines", 
            name="Radial Probability Density (Psi^2)") %>%
  add_trace(x=~rad, y=~(10^(-2) *sphericalSurfaceArea), type="scatter", mode="lines", 
            name="Spherical Surface Area (4 pi r^2)") %>%
  add_trace(x=~rad, y=~(10 * radialProbability), type="scatter", mode="lines", 
            name="Total Radial Probability (Psi^2 * r^2)") %>%
  layout(title="Radial Density, Surface Area & Total Probability",
           xaxis = list(title = 'Radius (r)'),
           yaxis = list(title ="Functional Values"),
           legend = list(title="Functions", orientation = 'h'))
```

Next, we can compare the *radial probability densities* of the first 3  radial wavefunctions $(n\in\{1,2,3\},\ l=0,\ m=0)$, 
$P_{n\ l}=R_{n\ l}(r)\cdot \left (\frac{1}{4\pi}\right )^2\cdot r^2$.
The electron *radial probability* is a function of distance from the nucleus ($r$) in all radial directions - the most probable radius increases with $n\uparrow$, however, compared to $1s$,
the $2s$ and $3s$ orbitals have regions of significant electron probability at smaller radii.
Note that the wavefunctions have radial *nodes* where the electron is unlikely to
appear; $R_{n\ 0}(r)$ has $n-1$ number of nodes - *shell boundaries.*

```{r electronProbabilities, echo=T, warning=F, error=F, message=F}
rad <- seq (from=10^(-5), to=15, length.out=1000)

df <- data.frame(cbind(radius=rad, RadProb_1_0=radialProbability(rad, 1,0), 
                       RadProb_2_0=radialProbability(rad, 2,0),
                       RadProb_3_0=radialProbability(rad, 3,0)))
df %>% plot_ly() %>%
  add_trace(x=~rad, y=~(RadProb_1_0), type="scatter", mode="lines", 
            name="Radial Prob 1s") %>%
  add_trace(x=~rad, y=~RadProb_2_0, type="scatter", mode="lines", 
            name="Radial Prob 2s") %>%
  add_trace(x=~rad, y=~RadProb_3_0, type="scatter", mode="lines", 
            name="Radial Prob 3s") %>%
  layout(title="Radial Probabilities (1s, 2s, 3s)",
           xaxis = list(title = 'Radius (r)'),
           yaxis = list(title ="Electron Probability (Psi^2 * r^2)"),
           legend = list(title="Functions", orientation = 'h'))
```


### Principal ($2p$) $l=1$ orbitals

The Principal $l=1$ orbital are described by the product of the *radial* and *angular* components

$$\psi_{n\ l\ m}(r,\theta,\phi)=\underbrace{R_{n\ 1}(r)}_{radial} \cdot 
\underbrace{Y_{1\ m}(\theta,\phi)}_{angular}.$$

The angular components are

$$Y_{1\ 0}(\theta,\phi)=\pm\left (\frac{3}{4\pi}\right )^{1/2} \cos\theta \ , \\
Y_{1\ \pm 1}(\theta,\phi)=\pm\left (\frac{3}{8\pi}\right )^{1/2} \sin\theta \ 
e^{\pm i\phi}\ .$$

These orbitals are not spherically symmetric, as the sharp $1S$ orbitals, however
they exhibit different types of symmetries.

Due to its synergy to the spherical coordinate transformation for $z$, $z=r\cos\theta$, the $m=0$ orbital is called $p_z$ *orbital*. This $p_z$ orbital
depends on $\theta$, via the $\cos\theta$ function, but is $\phi$ dependent. 

Below, we visualize the $2p$ in 3D. The *nodal plane* in the $p$
orbital corresponds to $\theta=\frac{\pi}{2}$, since $\forall \phi, \ \cos(\pi/2)=0$,
and the node boundary is a flat $(x,y)$ plane, a *nodal plane*.

The $\psi_{n\ 1\ 1}(r,\theta,\phi)$ and  $\psi_{n\ 1\ -1}(r,\theta,\phi)$
orbitals include *complex components*, since $\psi_{n\ 1\ \pm1}$ depends on
the term $e^{\pm i\phi}$. Even though these orbitals contain imaginary components, 
since $\psi_{n\ 1\ \pm1}$ are solutions of the Schrödinger 
equation with the same energy $E_n$, any linear combination of these two functions
will also be a solution of the Schrödinger equation at the same energy level. 
We can work with two combinations that both yield two real orbitals. For instance,
$p_x,p_y$,

$$\tilde{\psi}_{p_x} = \frac{1}{\sqrt{2}}(\psi_{n\ 1\ -1}(r,\theta,\phi)-
\psi_{n\ 1\ 1}(r,\theta,\phi))\ , \\
\tilde{\psi}_{p_y} = \frac{1}{\sqrt{2}}(\psi_{n\ 1\ -1}(r,\theta,\phi)+
\psi_{n\ 1\ 1}(r,\theta,\phi))\ .$$

The corresponding spherical harmonics are

$$Y_{p_x}(\theta,\phi) = \frac{1}{\sqrt{2}}(Y_{\ 1\ -1}(r,\theta,\phi)-
Y_{1\ 1}(r,\theta,\phi))=\left (\frac{3}{4\pi}\right )^{1/2}\sin\theta\ \cos\phi\ , \\
Y_{p_y}(\theta,\phi) = \frac{i}{\sqrt{2}}(Y_{\ 1\ -1}(r,\theta,\phi)+
Y_{1\ 1}(r,\theta,\phi))=\left (\frac{3}{4\pi}\right )^{1/2}\sin\theta\ \sin\phi\ .$$

Similarly to the $p_z$ orbital, the $p_x,\ p_y$ notation reflects the synergies with the spherical coordinate transformations 
$$x=r\sin\theta \ \cos\phi \\ y=r\sin\theta\ \sin\phi\ .$$

The $p_x,\ p_y$ orbitals have the same shape as the $p_z$ orbital, however
they are reoriented along the $x$-axis, for the $p_x$, and $y$-axis, for the $p_y$ orbitals.


### Diffuse ($3d$) $l=2$ orbitals

The $l=2$ *diffuse orbitals* also represent combinations of spherical harmonics 
of real orbitals. Corresponding linear combinations of the real orbitals
are denoted by $Y_{xy}, Y_{xz}, Y_{yz}, Y_{z^2}, Y_{x^2-y^2}$.
These five orbitals correspond to quantum numbers with $m\in \{-2,-1,0,1,2\}$. The notation signifies the *angular dependence* corresponding to the second-order
products $\{xy,\ xz,\ yz,\ z^2,\ x^2-y^2\}$.

The five $3d$ orbitals are shown in the 3D scene below. 
There are two *nodal surfaces*, i.e., *nodal planes*, in $4$ of the orbitals, 
$p_{xy},\ p_{xz},\ p_{yz},\ p_{x^2-y^2}$.

However, the last $3d_{z^2}$ orbital has a *nodal cone* boundary. 
The number of radial nodes is still $n-l-1$ and the overall number of 
nodes is unchanged. As $l\uparrow$, *radial nodes* are exchanged for 
*angular nodes.*

For any of the wavefunctions $\psi_{n\ l\ m}(r,\theta,\phi)$, averaging 
measurements of the distance between the electron and the origin (atomic nucleus)
multiple times yields the *expectation* $r$

$$\langle r \rangle=
\int_{0}^{\infty} { \int_{0}^{\pi} { 
\int_{-\pi}^{\pi} 
|\psi_{n\ l\ m}(r,\theta,\phi)|^2 r^3 \sin\theta\ dr\ d\theta\ d\phi}}=\\
\frac{n^2 a_o}{Z}\left (1+\frac{1}{2}\left ( 1-\frac{l(l+1)}{n^2}\right ) \right )\ .$$

## 3D orbitals

### Sharp ($s$) orbitals, $l=0$

```{r spherical_sharp_s_orbitals, echo=T, error=F, message=FALSE, warning=FALSE}
library(plotly)
# sweep or define (u,v) spherical coordinate parameter ranges
resol <- 50
theta <- seq(from = 0, to = pi, by = ((pi - 0)/(resol - 1)))
phi <- seq(from = -pi, to = pi, by = ((2*pi)/(resol - 1)))

#### l=0 Sharp orbitals (s) ####################################
R_n_0 <- function (r, n=1, a=a_o*(10^(10)), Z=1) {
  if (n==1) return ( R_1_0(r) )
  else if (n==2) return ( R_2_0(r) )
  else if (n==3) return ( R_3_0(r) )
  else { print('Invalid harmonic parameters (n,l) for R_n_0() ...')
    return()
  }
}

pdf_psi_n_0_0 <- function (r, n) {
  return ( R_n_0(r, n) * (1/(4*pi))^(1/2) )
} # Test: pdf_psi_n_0_0(1, 1)

shape_names <- c("1s", "2s", "3s")

# shape=="1s")  # Spheres
    # rendering (u,v) parametric surfaces requires x,y,z arguments to be 2D arrays
    # In out case, the three coordinates have to be 200*200 parameterized tensors/arrays
    r1 = 1.4                           # ~90% of RadialProbability  1s
    x1 = r1 * cos(phi) %o% sin(theta)   # x = r*cos(phi)*sin(psi)
    y1 = r1 * sin(phi) %o% sin(theta)   # y = r*sin(phi)*sin(psi)
    z1 = r1 * rep(1, resol) %o% cos(theta) # still need z to be 200*200 parameterized tensor/array

    
#shape=="2s") # Spheres
    # rendering (u,v) parametric surfaces requires x,y,z arguments to be 2D arrays
    # In out case, the three coordinates have to be 200*200 parameterized tensors/arrays
    r2_node1 = 1.3    # ~90% of RadialProbability  2s (node 1)
    x2_node1 = r2_node1 * cos(phi) %o% sin(theta)    # x = r*cos(phi)*sin(psi)
    y2_node1 = r2_node1 * sin(phi) %o% sin(theta)      # y = r*sin(phi)*sin(psi)
    z2_node1 = r2_node1*rep(1, resol) %o% cos(theta)  # z = r*cos(psi)
    
    r2_node2 = 3.1    # ~90% of RadialProbability  2s (node 2)
    x2_node2 = r2_node2 * cos(phi) %o% sin(theta)    # x = r*cos(phi)*sin(psi)
    y2_node2 = r2_node2 * sin(phi) %o% sin(theta)     # y = r*sin(phi)*sin(psi)
    z2_node2 = r2_node2*rep(1, resol) %o% cos(theta)  # z = r*cos(psi)
    
#shape=="3s")  # Spheres
    # rendering (u,v) parametric surfaces requires x,y,z arguments to be 2D arrays
    # In out case, the three coordinates have to be 200*200 parameterized tensors/arrays
    r3_node1 = 0.8    # ~90% of RadialProbability  2s (node 1)
    x3_node1 = r3_node1 * cos(phi) %o% sin(theta)      # x = r*cos(phi)*sin(psi)
    y3_node1 = r3_node1 * sin(phi) %o% sin(theta)      # y = r*sin(phi)*sin(psi)
    z3_node1 = r3_node1*rep(1, resol) %o% cos(theta)  # z = r*cos(psi)
    
    r3_node2 = 4    # ~90% of RadialProbability  2s (node 2)
    x3_node2 = r3_node2 * cos(phi) %o% sin(theta)     # x = r*cos(phi)*sin(psi)
    y3_node2 = r3_node2 * sin(phi) %o% sin(theta)      # y = r*sin(phi)*sin(psi)
    z3_node2 = r3_node2*rep(1, resol) %o% cos(theta)  # z = r*cos(psi)
    
    r3_node3 = 12    # ~90% of RadialProbability  2s (node 2)
    x3_node3 = r3_node3 * cos(phi) %o% sin(theta)     # x = r*cos(phi)*sin(psi)
    y3_node3 = r3_node3 * sin(phi) %o% sin(theta)      # y = r*sin(phi)*sin(psi)
    z3_node3 = r3_node3*rep(1, resol) %o% cos(theta)  # z = r*cos(psi)
    
# updatemenus component
updatemenus <- list(
  list(
    active = -1,
    type = 'buttons',
    buttons = list(
      list(
        label = shape_names[1],
        method = "update",
        args = list(list(visible = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE)),
                    list(title = shape_names[1]))),
      list(
        label = shape_names[2],
        method = "update",
        args = list(list(visible = c(FALSE, TRUE, TRUE, FALSE, FALSE, FALSE)),
                    list(title = shape_names[2]))),
      list(
        label = shape_names[3],
        method = "update",
        args = list(list(visible = c(FALSE, FALSE, FALSE, TRUE, TRUE, TRUE)),
                    list(title = shape_names[3])))
    )
  )
)

# p <- plot_ly(hoverinfo="none", legendshow=FALSE, showscale = FALSE) %>%
#   add_trace(x = ~x1, y = ~y1, z = ~z1, type = 'surface', opacity=0.3, visible=T,
#              contour=list(show=TRUE, color="#000", width=15, lwd=10,
#                           opacity=1.0, hoverinfo="none", legendshow=F)) %>%
#   add_trace(x = ~x2_node1, y = ~y2_node1, z = ~z2_node1, type='surface', opacity=0.6, visible=F,
#              contour=list(show=TRUE, color="#000", width=15, lwd=10,
#                           opacity=1.0, hoverinfo="none", legendshow=F)) %>%
#   add_trace(x = ~x2_node2, y = ~y2_node2, z = ~z2_node2, type='surface', opacity=0.3, visible=F,
#              contour=list(show=TRUE, color="#000", width=15, lwd=10,
#                           opacity=1.0, hoverinfo="none", legendshow=F)) %>%
#   add_trace(x = ~x3_node1, y = ~y3_node1, z = ~z3_node1, type = 'surface', opacity=0.9, visible=F,
#              contour=list(show=TRUE, color="#000", width=15, lwd=10,
#                           opacity=0.9, hoverinfo="none", legendshow=F)) %>%
#   add_trace(x = ~x3_node2, y = ~y3_node2, z = ~z3_node2, type = 'surface', opacity=0.4, visible=F,
#              contour=list(show=TRUE, color="#000", width=15, lwd=10,
#                           opacity=0.4, hoverinfo="none", legendshow=F)) %>%
#   add_trace(x = ~x3_node3, y = ~y3_node3, z = ~z3_node3, type = 'surface', opacity=0.2, visible=F,
#              contour=list(show=TRUE, color="#000", width=15, lwd=10,
#                           opacity=0.2, hoverinfo="none", legendshow=F)) %>%
#   layout(title = "Choose a Shape", showlegend = FALSE, updatemenus = updatemenus)
# p

plot_ly(hoverinfo="none",  showscale = FALSE) %>%
    add_trace(x = ~x1, y = ~y1, z = ~z1, type = 'surface', opacity=0.3, visible=T) %>%
    add_trace(x = ~x2_node1, y = ~y2_node1, z = ~z2_node1, type='surface', opacity=0.6, visible=F) %>%
    add_trace(x = ~x2_node2, y = ~y2_node2, z = ~z2_node2, type='surface', opacity=0.3, visible=F) %>%
    add_trace(x = ~x3_node1, y = ~y3_node1, z = ~z3_node1, type = 'surface', opacity=0.9, visible=F) %>%
    add_trace(x = ~x3_node2, y = ~y3_node2, z = ~z3_node2, type = 'surface', opacity=0.4, visible=F) %>%
    add_trace(x = ~x3_node3, y = ~y3_node3, z = ~z3_node3, type = 'surface', opacity=0.2, visible=F) %>%
    layout(title = "Choose a Shape", updatemenus = updatemenus)
```


### Principal orbitals ($p$), $l=1$ 

Recall the *quantum triple* $(n,l,m)$ includes $n$, the *principal quantum number*, $l$, the *angular quantum number*, and $m$, the *magnetic quantum number*.

### p1 orbitals: n=2, l=1, m=0
```{r spherical_principal_p_orbitals, echo=T, warning=F, error=F, message=F}
library(plotly)
Y_1_0 <- function (r, l=1, m=0, theta, phi) {
  return (  cos(theta) )
}
Y_1_1_mag <- function ( theta, phi) {
  return (  sin(theta)**2 )
}

R_2_1 <- function (r, a=a_o*(10^(10)), Z=1) {
  return ((r) * exp(-r))
}
pdf_psi_n_l_m <- function (r, n=2, l=1, m=0, theta, phi) {
  if (n==2 & l==1 & m==0) return ( r**2*R_2_1(r)**2  * Y_1_0(r, l=1, m=0, theta, phi)**2)
  else if (n==2 & l==1 & m==1) return ( r**2*R_2_1(r)**2  * Y_1_1_mag(theta, phi)**2)
  else { print('Invalid harmonic parameters (n,l) for pdf_psi_n_l_m() ...')
    return()
  }
}
meshgridn = function(L){
  out = list()
  n = sapply(L, length)
  w = 1:length(L)
  for(i in w){
    out[[i]] = rep(rep(L[[i]], rep(prod(n[w<i]), n[i])), prod(n[w>i]))
    # prod(NULL) == 1, so this works.
  }
  return(out)
}

gr <- meshgridn(list(x=seq(from=-5 , to=5, length.out=10), 
               y=seq(from=-5 , to=5, length.out=10),
               z=seq(from=-5, to=5, length.out=10)))



r = sqrt(gr[[1]]^2 + gr[[2]]^2 + gr[[3]]^2)
phi = atan2(gr[[2]],gr[[1]])
theta = acos(gr[[3]]/sqrt(gr[[1]]^2 + gr[[2]]^2+gr[[3]]^2))
# replace Nan with 0
theta[is.nan(theta)] <- 0
```

```{r p0_isosurface, echo=T, warning=F, error=F, message=F}
# n=2,l=1,m=0
p <-plot_ly(type='isosurface', x = gr[[1]], y = gr[[2]], z = gr[[3]],
  value = as.vector(100*abs(pdf_psi_n_l_m(r, n=2, l=1, m=0, theta, phi))), 
  isomin=6,  isomax=23.7, opacity=0.7, name="Psi_{nlm}(r,theta,phi)",
  colorscale='RdBu', surface = list(show = TRUE, count =2),
  slices = list(y = list(show = TRUE, locations = c(0)),
  x = list(show = TRUE, locations = c(0))))  %>%
  ### X-Y Plane (z=0)
  add_mesh(type = 'mesh3d', name="pz (m=0)",
             # Define all (4) zero-cells (points or vertices)
             # P_i(x,y,z), 0<=i<4
             x=c(-5, -5, 5, 5), y=c(-5, 5, -5, 5), z=c(0, 0, 0, 0),
             # Next define all triples (i,j,k) of vertices that form a 2-cell face.
             # All Tetrahedra have 4 faces
             i = c(0, 0), j = c(2, 1), k = c(3, 3),
             # Define the appearance of the 4 faces (2-cells)
             intensity = c(0.33, 0.33, 0.33, 0.33),
             color = "gray", opacity=0.6)%>% 
  hide_colorbar()
p
```

```{r p0_surface, echo=T, warning=F, error=F, message=F}
library(plot3D)
M <- mesh(seq(0, 2,length.out=2000), seq(0, 2*pi,length.out=50))
r <- exp(M$x)-1;  phi<- M$y
surface_para <- 0.1
costhe <- sqrt(surface_para/(r**4*exp(-2*r)))
sinthe <- sqrt(1-surface_para/(r**4*exp(-2*r)))
# patch the holes
sinthe[sinthe<0.1] <-0
x <- r*sinthe*cos(phi)
y <- r*sinthe*sin(phi)
z <- r*costhe
surface_para <- 0.15
costhe <- sqrt(surface_para/(r**4*exp(-2*r)))
sinthe <- sqrt(1-surface_para/(r**4*exp(-2*r)))
# patch the holes
sinthe[sinthe<0.1] <-0
x1 <- r*sinthe*cos(phi)
y1 <- r*sinthe*sin(phi)
z1 <- r*costhe
p <- plot_ly(x = x, y = y, z = z, type = "surface",colorscale = list(c(0,1),c("rgb(255,107,184)","rgb(128,0,64)")),showscale = FALSE, opacity=0.6) %>%
  layout(title = "Surface plot of the p0 orbital")%>% add_trace(x = x, y = y, z = -z, type="surface", colorscale = list(c(0,1),c("rgb(255,107,184)","rgb(128,0,64)")),showscale = FALSE, opacity=0.6) %>%add_trace(x = x1, y = y1, z = z1, type="surface",showscale = FALSE, opacity=0.6)%>%add_trace(x = x1, y = y1, z = -z1, type="surface",showscale = FALSE, opacity=0.6)

p
```

### p1 orbitals: n=2, l=1, m=1

```{r p1_isosurface, echo=T, warning=F, error=F, message=F}
r = sqrt(gr[[1]]^2 + gr[[2]]^2 + gr[[3]]^2)
phi = atan2(gr[[2]],gr[[1]])
theta = acos(gr[[3]]/sqrt(gr[[1]]^2 + gr[[2]]^2+gr[[3]]^2))
# replace Nan with 0
theta[is.nan(theta)] <- 0
p <-plot_ly(type='isosurface', x = gr[[1]], y = gr[[2]], z = gr[[3]],
  value = as.vector(100*abs(pdf_psi_n_l_m(r, n=2, l=1, m=1, theta, phi))), 
  isomin=6,  isomax=25, opacity=0.7, name="Psi_{nlm}(r,theta,phi)",
  colorscale='RdBu', surface = list(show = TRUE, count =1),
  slices = list(y = list(show = TRUE, locations = c(0)),
  x = list(show = TRUE, locations = c(0)),
  z = list(show = TRUE, locations = c(0))
  ))  %>%
  ### X-Y Plane (z=0)
  add_mesh(type = 'mesh3d', name="pz (m=0)",
             # Define all (4) zero-cells (points or vertices)
             # P_i(x,y,z), 0<=i<4
             x=c(-5, -5, 5, 5), y=c(-5, 5, -5, 5), z=c(0, 0, 0, 0),
             # Next define all triples (i,j,k) of vertices that form a 2-cell face.
             # All Tetrahedra have 4 faces
             i = c(0, 0), j = c(2, 1), k = c(3, 3),
             # Define the appearance of the 4 faces (2-cells)
             intensity = c(0.33, 0.33, 0.33, 0.33),
             color = "gray", opacity=0.6)%>% 
  hide_colorbar()
p
```

```{r p1_surface, echo=T, warning=F, error=F, message=F}
library(plot3D)
M <- mesh(seq(0, 2,length.out=2000), seq(0, 2*pi,length.out=50))
r <- exp(M$x)-1;  phi<- M$y
surface_para <- 0.1
sinthe <- sqrt(surface_para/(r**4*exp(-2*r)))
costhe <- sqrt(1-surface_para/(r**4*exp(-2*r)))
# patch the holes
# sinthe[sinthe>0.9] <-1
costhe[costhe<0.1] <-0
x <- r*sinthe*cos(phi)
y <- r*sinthe*sin(phi)
z <- r*costhe
surface_para <- 0.15
sinthe <- sqrt(surface_para/(r**4*exp(-2*r)))
costhe <- sqrt(1-surface_para/(r**4*exp(-2*r)))
# patch the holes
costhe[costhe<0.1] <-0
x1 <- r*sinthe*cos(phi)
y1 <- r*sinthe*sin(phi)
z1 <- r*costhe
p <- plot_ly(x = x, y = y, z = z, type = "surface",colorscale = list(c(0,1),c("rgb(255,107,184)","rgb(128,0,64)")),showscale = FALSE, opacity=0.6) %>%
  layout(title = "Surface plot of the p1 orbital")%>% add_trace(x = x, y = y, z = -z, type="surface", colorscale = list(c(0,1),c("rgb(255,107,184)","rgb(128,0,64)")),showscale = FALSE, opacity=0.6) %>%add_trace(x = x1, y = y1, z = z1, type="surface",showscale = FALSE, opacity=0.6)%>%add_trace(x = x1, y = y1, z = -z1, type="surface",showscale = FALSE, opacity=0.6)

p
```



### Diffuse orbitals ($d$), $l=2$ 

In this case, the *quantum triple* $(n,l,m)$ includes the *principal quantum number*,$n$, the *angular quantum number*, $l=2$, and the *magnetic quantum number*, $m\in\{-2,-1,0,1,2\}$.

### n=3,l=2,m=0
```{r spherical_diffuse_d_orbitals, echo=T, warning=F, error=F, message=F}
#### l=3 Diffuse orbitals (d) 
R_n_3 <- function (r, n=3, a=a_o*(10^(10)), Z=1) {
  if (n==1) { print('Invalid harmonic parameters (n,l) for R_n_0() ...')
    return()
  } else if (n==2) { print('Invalid harmonic parameters (n,l) for R_n_0() ...')
    return()
  }  else if (n==3) return ( R_3_3(r) )   # default
  else { print('Invalid harmonic parameters (n,l) for R_n_0() ...')
    return()
  }
}

# Radial Density: P_{n\ l}(r)=r^2R^2_{n\ l}(r)
P_n_l_r <- function (r, n=2, l=0) {
  return ( r^2 * (R_n_1(r))^2 )
}

# Y_1_0 <- function (r, l=1, m=0, theta, phi) {
#   return ( (3/(4*pi))^(1/2) * cos(theta) %o% rep(1, resol) )
# }
Y_1_0 <- function (r, l=1, m=0, theta, phi) {
  return ( rep((3/(4*pi))^(1/2), resol) %o% cos(theta) )
}

# radial parts of the wavefunctions using pm picometer (10^{-12} of a meter) distance
R_2_0 <- function (r, a=a_o*(10^(10)), Z=1) {
  return ( (1/(2*sqrt(2))) * (Z/a)^(3/2) * (2- (Z*r)/(a)) * exp(-Z*r/(2*a)) )
}

R_2_1 <- function (r, a=a_o*(10^(10)), Z=1) {
  return ( (1/(2*sqrt(6))) * (Z/a)^(3/2) * (Z*r/a) * exp(-Z*r/(2*a)))
}

R_3_0 <- function (r, a=a_o*(10^(10)), Z=1) {
  return ( (2/(81*sqrt(3))) * (Z/a)^(3/2) * (27- 18*(Z*r/a) +
                                                2*(Z*r/a)^2) * exp(-Z*r/(3*a)))
}

R_3_2 <- function (r, a=a_o*(10^(10)), Z=1) {
  rho=(2*Z*r)/3
  #return ( (1/(9*sqrt(30))) * (rho)^2 * Z^(3/2) * exp(-rho/2) )
  return ( (rho)^2 *exp(-rho/2) )
}

Y_2_0 <- function (r, a=a_o*(10^(10)), Z=1, theta=pi/2, phi=0) { # cone node
  #return  ( ((sqrt(5/4))) * (3*cos(theta)**2-1) * sqrt(1/(4*pi)) )
  return  (  (3*cos(theta)**2-1))
}  # Y_3_z2(r=10, theta=pi/3)

Y_2_1 <- function (r, a=a_o*(10^(10)), Z=1, theta=pi/2, phi=0) { # cone node
  return  (  -sin(theta)*cos(theta)) # ignore the phase phi factor
}  
Y_2_2 <- function (r, a=a_o*(10^(10)), Z=1, theta=pi/2, phi=0) { # cone node
  return  (  sin(theta)**2) # ignore the phase phi factor
}  

# \psi_{n\ l\ m}(r,\theta,\phi)=\underbrace{R_{n\ l}(r)}_{radial} \cdot 
# \underbrace{Y_{l\ m}(\theta,\phi)}_{angular}, \quad n,l,m\in\mathbb{Z}
pdf_psi_n_l_m <- function (r, n=3, l=2, m=0, theta, phi) {
  if (n==3 & l==2 & m==0) return ( r**2*R_3_2(r)**2  * Y_2_0(r, theta=theta, phi=phi)**2 )
  else if(n==3 & l==2 & m==1) return ( r**2*R_3_2(r)**2  * Y_2_1(r, theta=theta, phi=phi)**2 )
  else if(n==3 & l==2 & m==2) return ( r**2*R_3_2(r)**2  * Y_2_2(r, theta=theta, phi=phi)**2 )
  else { print('Invalid harmonic parameters (n,l) for pdf_psi_n_l_m() ...')
    return()
  }
}   # pdf_psi_n_l_m(r=10, theta=pi/3)
gr <- meshgridn(list(x=seq(from=-20 , to=20, length.out=20), 
               y=seq(from=-20 , to=20, length.out=20),
               z=seq(from=-20, to=20, length.out=20)))


r = sqrt(gr[[1]]^2 + gr[[2]]^2 + gr[[3]]^2)
phi = atan2(gr[[2]],gr[[1]])
theta = acos(gr[[3]]/sqrt(gr[[1]]^2 + gr[[2]]^2+gr[[3]]^2))
# replace Nan with 0
theta[is.nan(theta)] <- 0
# n=3,l=2,m=0
#value = as.vector(abs(pdf_psi_n_l_m(r, n=3, l=2, m=0, theta, phi))),

p <-plot_ly(type='isosurface', x = gr[[1]], y = gr[[2]], z = gr[[3]],
  value = as.vector(abs(pdf_psi_n_l_m(r, n=3, l=2, m=0, theta, phi))), 
  isomin=0.0,  isomax=955, opacity=0.7, name="Psi_{nlm}(r,theta,phi)",
  colorscale='RdBu', surface = list(show = TRUE, count =1),
  slices = list(y = list(show = TRUE, locations = c(0)),
  x = list(show = TRUE, locations = c(0))))  %>%
  ### X-Y Plane (z=0)
  add_mesh(type = 'mesh3d', name="pz (m=0)",
             # Define all (4) zero-cells (points or vertices)
             # P_i(x,y,z), 0<=i<4
             x=c(-5, -5, 5, 5), y=c(-5, 5, -5, 5), z=c(0, 0, 0, 0),
             # Next define all triples (i,j,k) of vertices that form a 2-cell face.
             # All Tetrahedra have 4 faces
             i = c(0, 0), j = c(2, 1), k = c(3, 3),
             # Define the appearance of the 4 faces (2-cells)
             intensity = c(0.33, 0.33, 0.33, 0.33),
             color = "gray", opacity=0.6)%>% 
  hide_colorbar()
p

```


```{r dz2_orbitals, echo=T, warning=F, error=F, message=F}
## Ignore scalar factors
## (r**2*exp(-r/3))**2*r**2*(3*cos-1)**2
# Changing the surface_para gives different probability isocontours
library(plot3D)
M <- mesh(seq(0, 100,length.out=2000), seq(0, 2*pi,length.out=50))
r <- M$x;  phi<- M$y
surface_para <- 600
costhe <- sqrt((sqrt(surface_para/(r**6*exp(-2/3*r)))+1)/3)
sinthe <- sqrt(1-costhe**2)
# patch the holes
sinthe[sinthe<0.1] <-0
x <- r*sinthe*cos(phi)
y <- r*sinthe*sin(phi)
z <- r*costhe
surface_para <- 600
costhe <- sqrt((-sqrt(surface_para/(r**6*exp(-2/3*r)))+1)/3)
sinthe <- sqrt(1-costhe**2)
# patch the holes
costhe[costhe<0.1] <-0
x1 <- r*sinthe*cos(phi)
y1 <- r*sinthe*sin(phi)
z1 <- r*costhe

p <- plot_ly(x = x, y = y, z = z, type = "surface",colorscale = list(c(0,1),c("rgb(255,107,184)","rgb(128,0,64)")),showscale = FALSE, opacity=0.6) %>%
  layout(title = "Surface plot of the dz^2 orbital")%>% add_trace(x = x, y = y, z = -z, type="surface", colorscale = list(c(0,1),c("rgb(255,107,184)","rgb(128,0,64)")),showscale = FALSE, opacity=0.6) %>%add_trace(x = x1, y = y1, z = z1, type="surface",showscale = FALSE, opacity=0.6)%>%add_trace(x = x1, y = y1, z = -z1, type="surface",showscale = FALSE, opacity=0.6)

p

```

### n=3,l=2,m=1
```{r dxz-orbi, echo=T, warning=F, error=F, message=F}

r = sqrt(gr[[1]]^2 + gr[[2]]^2 + gr[[3]]^2)
phi = atan2(gr[[2]],gr[[1]])
theta = acos(gr[[3]]/sqrt(gr[[1]]^2 + gr[[2]]^2+gr[[3]]^2))
# replace Nan with 0
theta[is.nan(theta)] <- 0
# n=3,l=2,m=0
#value = as.vector(abs(pdf_psi_n_l_m(r, n=3, l=2, m=0, theta, phi))),

p <-plot_ly(type='isosurface', x = gr[[1]], y = gr[[2]], z = gr[[3]],
  value = as.vector(abs(pdf_psi_n_l_m(r, n=3, l=2, m=1, theta, phi))), 
  isomin=0.0,  isomax=62.2, opacity=0.7, name="Psi_{nlm}(r,theta,phi)",
  colorscale='RdBu', surface = list(show = TRUE, count =1),
  slices = list(y = list(show = TRUE, locations = c(0)),x = list(show = TRUE, locations = c(0)),
  z = list(show = TRUE, locations = c(0))))  %>%
  ### X-Y Plane (z=0)
  add_mesh(type = 'mesh3d', name="pz (m=0)",
             # Define all (4) zero-cells (points or vertices)
             # P_i(x,y,z), 0<=i<4
             x=c(-20, -20, 20, 20), y=c(-20, 20, -20, 20), z=c(0, 0, 0, 0),
             # Next define all triples (i,j,k) of vertices that form a 2-cell face.
             # All Tetrahedra have 4 faces
             i = c(0, 0), j = c(2, 1), k = c(3, 3),
             # Define the appearance of the 4 faces (2-cells)
             intensity = c(0.33, 0.33, 0.33, 0.33),
             color = "gray", opacity=0.6)%>% 
  hide_colorbar()
p
```


```{r dxz-surf, echo=T, warning=F, error=F, message=F}
## Ignore scalar factors
## (r**2*exp(-r/3))**2*r**2*sin(2*theta)**2
library(plot3D)
M <- mesh(seq(0, 50,length.out=2000), seq(0, 2*pi,length.out=50))
r <- M$x;  phi<- M$y
surface_para <- 1
theta_d<- pi/2-asin((sqrt(surface_para/(r**6*exp(-2/3*r)))))
theta_d[is.nan(theta_d)] <- 0
theta_d <- theta_d
# patch the holes
# costhe[costhe<0.1] <-0
x1 <- r*sin(theta_d/2)*cos(phi)
y1 <- r*sin(theta_d/2)*sin(phi)
z1 <- r*cos(theta_d/2)
# clockwise rotation
x_tran <- cos(pi/4)*x1-sin(pi/4)*z1
z_tran <- sin(pi/4)*x1+cos(pi/4)*z1

p <- plot_ly() %>%add_trace(x = x_tran, y = y1, z = z_tran, type="surface",showscale = FALSE, opacity=0.6)%>%add_trace(x = -x_tran, y = -y1, z = -z_tran, type="surface",showscale = FALSE, opacity=0.6)
# counter clockwise
x_tran <- cos(-pi/4)*x1-sin(-pi/4)*z1
z_tran <- sin(-pi/4)*x1+cos(-pi/4)*z1
p <- p%>%add_trace(x = x_tran, y = y1, z = z_tran, type="surface",showscale = FALSE, opacity=0.6)%>%add_trace(x = -x_tran, y = -y1, z = -z_tran, type="surface",showscale = FALSE, opacity=0.6)%>%
  layout(title = "Surface plot of the dxz orbital")
p

```

The other orbitals can be obtained similarly through rotation of different axes.

### n=3,l=2,m=2

```{r dx2-y2-orbital, echo=T, warning=F, error=F, message=F}

r = sqrt(gr[[1]]^2 + gr[[2]]^2 + gr[[3]]^2)
phi = atan2(gr[[2]],gr[[1]])
theta = acos(gr[[3]]/sqrt(gr[[1]]^2 + gr[[2]]^2+gr[[3]]^2))
# replace Nan with 0
theta[is.nan(theta)] <- 0
# n=3,l=2,m=0
#value = as.vector(abs(pdf_psi_n_l_m(r, n=3, l=2, m=0, theta, phi))),

p <-plot_ly(type='isosurface', x = gr[[1]], y = gr[[2]], z = gr[[3]],
  value = as.vector(abs(pdf_psi_n_l_m(r, n=3, l=2, m=2, theta, phi))), 
  isomin=0.0,  isomax=253.2, opacity=0.7, name="Psi_{nlm}(r,theta,phi)",
  colorscale='RdBu', surface = list(show = TRUE, count =1),
  slices = list(y = list(show = TRUE, locations = c(0)),x = list(show = TRUE, locations = c(0)),
  z = list(show = TRUE, locations = c(0))))  %>%
  ### X-Y Plane (z=0)
  add_mesh(type = 'mesh3d', name="pz (m=0)",
             # Define all (4) zero-cells (points or vertices)
             # P_i(x,y,z), 0<=i<4
             x=c(-20, -20, 20, 20), y=c(-20, 20, -20, 20), z=c(0, 0, 0, 0),
             # Next define all triples (i,j,k) of vertices that form a 2-cell face.
             # All Tetrahedra have 4 faces
             i = c(0, 0), j = c(2, 1), k = c(3, 3),
             # Define the appearance of the 4 faces (2-cells)
             intensity = c(0.33, 0.33, 0.33, 0.33),
             color = "gray", opacity=0.6)%>% 
  hide_colorbar()
p
```
```{r dx2-y2-surf, echo=T, warning=F, error=F, message=F}
## Ignore scalar factors
## (r**2*exp(-r/3))**2*r**2*sin(theta)**2**2
library(mosaic)
M <- mesh(seq(pi/4, 3*pi/4,length.out=20), seq(0, 2*pi,length.out=50))
theta <- M$x;  phi<- M$y
surface_para <- 24
f_r <- surface_para/(sin(theta)**4)
x <- c()
y <- c()
z <- c()
for (i in 1:c(20*50)){
  res <- f_r[i]
  candidates <- findZeros(exp(-2/3*t)*t**6 -res~ t, xlim=c(0,30) )
  if (length(candidates$t)>0){
    for (j in candidates$t){
      if(j>0){
        x <- c(x, j*sin(theta[i])*cos(phi[i]))
    y <- c(y, j*sin(theta[i])*sin(phi[i]))
    z <- c(z, j*cos(theta[i]))
      }
    }
  }
}


```

```{r}
x <- c()
y <- c()
z <- c()
for (i in 1:2){
  res <- f_r[i]
  candidates <- findZeros(exp(-2/3*t)*t**6 -res~ t, xlim=c(0,30) )
  if (length(candidates$t)>0){
    for (j in candidates$t){
      if(j>0){
        x <- c(x, j*sin(theta[i])*cos(phi[i]))
        y <- c(y, j*sin(theta[i])*sin(phi[i]))
        z <- c(z, j*cos(theta[i]))
      }
    }
  }
}
```

```{r dx2-y2-scatter,warning=FALSE, message=FALSE}
# (x**2+y**2+z**2)*exp(-2/3*sqrt(x^2+y^2+z^2))*(x**2-y**2)
M <- mesh(seq(-30, 30,length.out=30), seq(-29, 31,length.out=30))
x_g <- M$x;  y_g<- M$y
surface_para <- 24
f_r <- surface_para/(sin(theta)**4)
x <- c()
y <- c()
z <- c()
for (i in 1:c(30*30)){
  x2my2 <- x_g[i]**2-y_g[i]**2
  x2py2 <- x_g[i]**2+y_g[i]**2
  candidates <- findZeros( (x2py2+t**2)*exp(-2/3*sqrt(x2py2+t**2))*x2my2**2-surface_para~ t, xlim=c(-30,30) )
  if (length(candidates)!=0 ){
    for (j in candidates$t){
        x <- c(x, x_g[i])
    y <- c(y, y_g[i])
    z <- c(z, j)
    }
  }
}
p <- plot_ly(x = x, y = y, z = z, type = "scatter3d",colorscale = list(c(0,1),c("rgb(255,107,184)","rgb(128,0,64)")),showscale = FALSE, opacity=0.6) %>%
  layout(title = "Scatter plot of the dx2-y2 orbital")
p
```







## Mathematica (Wolfram) script

```{r Mathematica_script_orbitals, eval=F, echo=T, warning=F, error=F, message=F}
In[1]:= R21[Z_, r_, a_] := (r (Z/a)^(5/2))/(2 Sqrt[6] E^((r Z)/(2 a)))

In[2]:= Expand[R21[Z, r, a]]

Out[2]= (E^(-((r Z)/(2 a))) r Z^2 Sqrt[Z/a])/(2 Sqrt[6] a^2)

In[3]:= R21[1, 1, 1]

Out[3]= 1/(2 Sqrt[6 E])

In[4]:= Y10[theta_, phi_] := -Sqrt[3/(4 Pi)] * Cos[theta]

In[5]:= Y10[Pi, Pi/3]

Out[5]= Sqrt[3/\[Pi]]/2

In[6]:= Psi210[Z_, r_, a_, theta_, phi_] :=
 R21[Z, r, a] * Y10[theta, phi]

In[7]:= Psi210[1, 1, 1, Pi, Pi/3]

Out[7]= 1/(4 Sqrt[2 E \[Pi]])

In[8]:= Expand[Psi210[Z, r, a, theta, phi]]

Out[8]= -((E^(-((r Z)/(2 a))) r Z^2 Sqrt[Z/a] Cos[theta])/(
 4 a^2 Sqrt[2 \[Pi]]))

In[9]:= MyIntegrand[Z_, r_, a_, theta_, phi_] :=
 Psi210[Z, r, a, theta, phi]^2 * r^2 * Sin[theta]

In[10]:= MyIntegrand[1, 1, 1, Pi/10, Pi/3]

Out[10]= ((5/8 + Sqrt[5]/8) (-1 + Sqrt[5]))/(128 E \[Pi])

In[30]:= Integrate[
 MyIntegrand[Z, r, a, theta, phi], {r, a, Infinity}, {theta, 0,
  Pi}, {phi, -Pi, Pi}]

Out[30]= ConditionalExpression[
 1/24 E^-Z (24 + Z (24 + Z (12 + Z (4 + Z)))), Re[Z/a] > 0]

In[31]:= Evaluate[%30, Z = 1]

Out[31]= Sequence[ConditionalExpression[65/(24 E), Re[1/a] > 0], 1]

In[12]:= Expand[MyIntegrand[Z, r, a, theta, phi]]

Out[12]= (E^(-((r Z)/a)) r^4 Z^5 Cos[theta]^2 Sin[
  theta])/(32 a^5 \[Pi])

In[13]:= Integrate[
 Psi210[Z, r, a, theta, phi]^2  , {r, 1, 10}, {theta, -Pi/2,
  Pi/2}, {phi, -Pi, Pi}]

Out[13]= (E^(-((10 Z)/
  a)) \[Pi] Z^2 (-2 a^2 - 20 a Z - 100 Z^2 +
   E^((9 Z)/a) (2 a^2 + 2 a Z + Z^2)))/(32 a^4)

In[14]:= Integrate[
 Psi210[1, r, a, theta, phi], {r, 1, 10}, {theta, -Pi/2,
  Pi/2}, {phi, -Pi, Pi}]

Out[14]= (1/a)^(3/2) E^(-5/a) (10 + 2 a - (1 + 2 a) E^((9/2)/a)) Sqrt[
 2 \[Pi]]

In[15]:= (* R10 *)
R10[Z_, r_, a_] := 2 ((Z/a)^(3/2))*( E^((-r Z)/(a)))

In[16]:= Expand[R10[Z, r, a]]

Out[16]= (2 E^(-((r Z)/a)) Z Sqrt[Z/a])/a

In[17]:= Y00[theta_, phi_] := Sqrt[1/(4 Pi)]

In[18]:= Expand[Y00[theta, phi]]

Out[18]= 1/(2 Sqrt[\[Pi]])

In[19]:= Expand[R10[Z, r, a] * Y00[theta, phi]]

Out[19]= (E^(-((r Z)/a)) Z Sqrt[Z/a])/(a Sqrt[\[Pi]])

In[21]:= MyIntegrand100[Z_, r_, a_, theta_,
  phi_] := (R10[Z, r, a] * Y00[theta, phi])^2 * r^2 * Sin[theta]

In[22]:= Expand[MyIntegrand100[Z, r, a, theta, phi]]

Out[22]= (E^(-((2 r Z)/a)) r^2 Z^3 Sin[theta])/(a^3 \[Pi])

In[29]:= Integrate[
 MyIntegrand100[Z, r, a, theta, phi], {r, a, Infinity}, {theta, 0,
  Pi}, {phi, -Pi, Pi}]

Out[29]= ConditionalExpression[E^(-2 Z) (1 + 2 Z (1 + Z)), Re[Z/a] > 0]

In[32]:= Evaluate[%29, Z = 1]

Out[32]= Sequence[ConditionalExpression[5/E^2, Re[1/a] > 0], 1]
```

# 3D Scene of Circular Kime-Phase Sampling

The following experiment shows a mock up of 5D spacekime where the three spatial dimensions are compressed into 1D (gray-colored vertical z-axis). The 3D scene shows the cross-product of 1D space and 2D kime. The random *circular kime-phase draws*, $\{\varphi_i\}_{i}^I \sim \Phi [-\pi,\pi)$ are shown for 3 different processes along with their corresponding pooled (averaged) estimates representing the *observed (measurable) data* reflecting the *observed sampling distribution* of the average phase for each of the 3 independent processes, 
$$\underbrace{\overline{\varphi}}_{measurable\\ quantity}=\frac{1}{I}\sum_i{\varphi_i}.$$
Recall that by [Central Limit Theorem](http://wiki.stat.ucla.edu/socr/index.php/SOCR_EduMaterials_Activities_GeneralCentralLimitTheorem), see [the SOCR CLT papper](https://doi.org/10.1080/10691898.2008.11889560), for all nice phase distributions, $\Phi [-\pi,\pi)$, $\overline{\varphi}\sim N\left (\mu_{\Phi}, \frac{\sigma_{\Phi}^2}{I}\right )$. As the phase distributions are symmetric, $\mu_{\overline{\varphi}}=\mathbb{E}(\overline{\varphi})=\mu_{\Phi}=0$, and 
$$\sigma_{\overline{\varphi}}= \frac{\sigma_{\Phi}^2}{I}\underbrace{\longrightarrow}_{I\to\infty}0.$$

**Caution**: In the simulation below, to emphasize the differences between different observable processes, we are *artificially offsetting (shifting) the random samples*, and their corresponding mean trajectories (observable patterns). These offsets are both in the vertical (*space*) dimension as well as in the radial (*phase*) space. These offsets are are completely artificial and just intended to enhance the interpretation of the kime-phase sampling process by avoiding significant overlaps of points and observable kime-dynamic trends.

```{r fig.height=10, fig.width=10, warning=FALSE, message=FALSE, error=FALSE}
library(animation)
epsilon <- 0.1
sampleSize <- 1000   # total number of phases to sample for 3 different processes (x, y, z)
sizePerTime <- 100   # number of phases to use for each fixed time (must divide sampleSize)
circleUniformPhi <- seq(from=-pi, to=pi, length.out=sizePerTime)

oopt = ani.options(interval = 0.2)
set.seed(1234)
# sample the the kime-phases for all 3 different processes and the r time points
x <- rvonmises(n=sampleSize, mu=circular(pi/5), kappa=3)
y <- rvonmises(n=sampleSize, mu=circular(-pi/3), kappa=5)
z <- rvonmises(n=sampleSize, mu=circular(0), kappa=10)
r <- seq(from=1, to=sampleSize/sizePerTime, length.out=10)

# Define a function that renormalizes the kime-phase to [-pi, pi)
pheRenormalize <- function (x) {
  out <- ifelse(as.numeric(x) <= pi, as.numeric(x)+pi, as.numeric(x)-pi)
  return (out)
}

# transform Von Mises samples from [0, 2*pi) to [-pi, pi)
x <- pheRenormalize(x)
y <- pheRenormalize(y)
z <- pheRenormalize(z)

# vectorize the samples
vectorX = as.vector(x)
vectorY = as.vector(y)
vectorZ = as.vector(z)
# Starting phases, set the first phase index=1
plotX = c(vectorX[1])
plotY = c(vectorY[1])
plotZ = c(vectorZ[1])

# # iterate over all time points (outer loop) and all phases (inner loop)
# for (t in 1:length(r)) {  # loop over time
#   for (i in 2:100) {      # loop over kime-phases
#     plotX = c(plotX,vectorX[i])
#     plotY = c(plotY,vectorY[i])
#     plotZ = c(plotZ,vectorZ[i])
# 
#     # Try to "stack=T the points ....
#     #r1 = sqrt((resx$x)^2+(resx$y)^2)/2;
#     #resx$x = r1*cos(resx$data)
#     #resx$x = r1*cos(resx$data)
# 
#     tempX = circular(plotX[[1]])
#     tempY = circular(plotY[[1]])
#     tempZ = circular(plotZ[[1]])
#     resx <- density(tempX, bw=25, xaxt='n', yaxt='n')
#     resy <- density(tempY, bw=25, xaxt='n', yaxt='n')
#     resz <- density(tempZ, bw=25, xaxt='n', yaxt='n')
#     res <- plot(resx, points.plot=TRUE, xlim=c(-1.1,1.5), ylim=c(-1.5, 1.5),
#                 main="Trivariate random sampling of\n kime-magnitudes (times) and kime-directions (phases)",
#                 offset=0.9, shrink=1.0, ticks=T, lwd=3, axes=F, stack=TRUE, bins=150)
#     lines(resy, points.plot=TRUE, col=2, points.col=2, plot.info=res, offset=1.0, shrink=1.45, lwd=3, stack=T)
#     lines(resz, points.plot=TRUE, col=3, points.col=3, plot.info=res, offset=1.1, shrink=1.2, lwd=3, stack=T)
#     segments(0, 0, r[i]*cos(vectorX[i]), r[i]*sin(vectorX[i]), lwd=2, col= 'black')
#     segments(0, 0, r[i]*cos(vectorY[i]), r[i]*sin(vectorY[i]), lwd=2, col= 'red')
#     segments(0, 0, r[i]*cos(vectorZ[i]), r[i]*sin(vectorZ[i]), lwd=2, col= 'green')
#     ani.pause()
#   }
# }
####################################################

# pl_list <- list()
pl_scene <- plot_ly(type='scatter3d', mode="markers")
plotX <- list() 
plotY <- list() 
plotZ <- list() 

plotX_df <- list()   # need separate dataframes to store all time foliations
plotY_df <- list()
plotZ_df <- list()

for (t in 1:length(r)) {  # loop over time
  # loop over kime-phases
  plotX[[t]] <- as.numeric(x[c(( (t-1)*length(r) + 1):((t-1)*length(r) + sizePerTime))])
  plotY[[t]] <- as.numeric(y[c(( (t-1)*length(r) + 1):((t-1)*length(r) + sizePerTime))])
  plotZ[[t]] <- as.numeric(z[c(( (t-1)*length(r) + 1):((t-1)*length(r) + sizePerTime))])
  
    # Try to "stack=T the points ....
    #r1 = sqrt((resx$x)^2+(resx$y)^2)/2;
    #resx$x = r1*cos(resx$data)
    #resx$x = r1*cos(resx$data)

  tempX = circular(unlist(plotX[[t]]))
  tempY = circular(unlist(plotY[[t]]))
  tempZ = circular(unlist(plotZ[[t]]))
  
  resx <- density(tempX, bw=25, xaxt='n', yaxt='n')
  resy <- density(tempY, bw=25, xaxt='n', yaxt='n')
  resz <- density(tempZ, bw=25, xaxt='n', yaxt='n')
  # res <- plot(resx, points.plot=TRUE, xlim=c(-1.1,1.5), ylim=c(-1.5, 1.5),
  # main="Trivariate random sampling of\n kime-magnitudes (times) and kime-directions (phases)",
  # offset=0.9, shrink=1.0, ticks=T, lwd=3, axes=F, stack=TRUE, bins=150)
  # pl_list[[t]] 
  unifPhi_df <- as.data.frame(cbind(t=t, circleUniformPhi=circleUniformPhi))
  plotX_df[[t]] <- as.data.frame(cbind(t=t, plotX=unlist(plotX[[t]])))
  plotY_df[[t]] <- as.data.frame(cbind(t=t, plotY=unlist(plotY[[t]])))
  plotZ_df[[t]] <- as.data.frame(cbind(t=t, plotZ=unlist(plotZ[[t]])))
  
  pl_scene <- pl_scene %>% add_trace(data=unifPhi_df, showlegend=FALSE,
                      x = ~((t-epsilon)*cos(circleUniformPhi)), 
                      y = ~((t-epsilon)*sin(circleUniformPhi)), z=0,
                      name=paste0("Time=",t), line=list(color='gray'),
                      mode = 'lines', opacity=0.3) %>%
    add_markers(data=plotX_df[[t]], x=~(t*cos(plotX)), y=~(t*sin(plotX)), z=0,
                      type='scatter3d', name=paste0("X: t=",t), 
                      marker=list(color='green'), showlegend=FALSE,
                      mode = 'markers', opacity=0.3) %>%
    add_markers(data=plotY_df[[t]], x=~((t+epsilon)*cos(plotY)),
                    y=~((t+epsilon)*sin(plotY)), z=0-epsilon, showlegend=FALSE,
                    type='scatter3d', name=paste0("Y: t=",t), 
                    marker=list(color='blue'),
                    mode = 'markers', opacity=0.3) %>%
    add_markers(data=plotZ_df[[t]], x=~((t+2*epsilon)*cos(plotZ)),
                y=~((t+2*epsilon)*sin(plotZ)), z=0+epsilon, showlegend=FALSE,
                type='scatter3d', name=paste0("Z: t=",t), 
                marker=list(color='red'),
                mode = 'markers', opacity=0.3)
} 

means_df <- as.data.frame(cbind(t = c(1:length(r)),
                                plotX_means=unlist(lapply(plotX, mean)),
                                plotY_means=unlist(lapply(plotY, mean)),
                                plotZ_means=unlist(lapply(plotZ, mean))))
pl_scene <- pl_scene %>% 
  # add averaged (donoised) phase trajectories
  add_trace(data=means_df, x=~(t*cos(plotX_means)), 
        y=~(t*sin(plotX_means)), z=0,
        type='scatter3d', showlegend=FALSE, mode='lines', name="Expected Obs X", 
        line=list(color='green', width=15), opacity=0.8) %>%
  add_trace(data=means_df, x=~(t*cos(plotY_means)), 
        y=~(t*sin(plotY_means)), z=0-epsilon,
        type='scatter3d', showlegend=FALSE, mode='lines', name="Expected Obs X", 
        line=list(color='blue', width=15), opacity=0.8) %>%
  add_trace(data=means_df, x=~(t*cos(plotZ_means)), 
        y=~(t*sin(plotZ_means)), z=0+epsilon,
        type='scatter3d', showlegend=FALSE, mode='lines', name="Expected Obs X", 
        line=list(color='red', width=15), opacity=0.8) %>%
  add_trace(x=0, y=0, z=c(-2,2), name="Space", showlegend=FALSE,
             line=list(color='gray', width=15), opacity=0.8) %>%
  layout(title="Pseudo Spacekime (1D Space, 2D Kime) Kime-Phase Sampling and Foliation",
          scene = list(xaxis=list(title="Kappa1"), yaxis=list(title="Kappa2"),
                        zaxis=list(title="Space"))) %>% hide_colorbar()
pl_scene

# pl_list %>%
#   subplot(nrows = length(pl_list)) %>%
#      layout(title="Integral Approximations by Riemann Sums", legend = list(orientation = 'h'))
```


<!--html_preserve-->
<div>
    	<footer><center>
			<a href="https://www.socr.umich.edu/">SOCR Resource</a>
				Visitor number <img class="statcounter" src="https://c.statcounter.com/5714596/0/038e9ac4/0/" alt="Web Analytics" align="middle" border="0">
				<script type="text/javascript">
					var d = new Date();
					document.write(" | " + d.getFullYear() + " | ");
				</script> 
				<a href="https://socr.umich.edu/img/SOCR_Email.png"><img alt="SOCR Email"
	 			title="SOCR Email" src="https://socr.umich.edu/img/SOCR_Email.png"
	 			style="border: 0px solid ;"></a>
	 		 </center>
	 	</footer>

	<!-- Start of StatCounter Code -->
		<script type="text/javascript">
			var sc_project=5714596; 
			var sc_invisible=1; 
			var sc_partition=71; 
			var sc_click_stat=1; 
			var sc_security="038e9ac4"; 
		</script>
		
		<script type="text/javascript" src="https://www.statcounter.com/counter/counter.js"></script>
	<!-- End of StatCounter Code -->
	
	<!-- GoogleAnalytics -->
		<script src="https://www.google-analytics.com/urchin.js" type="text/javascript"> </script>
		<script type="text/javascript"> _uacct = "UA-676559-1"; urchinTracker(); </script>
	<!-- End of GoogleAnalytics Code -->
</div>
<!--/html_preserve-->