Is $g(\eta)$, as defined in the question, a Dirac delta function? Section 1: Maths Question (TL;DR version)
In the course of trying to solve a physics problem (ref. Section 2), I encountered a mathematical question. To make my post brief, I'll write only the maths question here that needs to be addressed:
\begin{align} &f(r,\eta)= -\frac{r-R\eta}{(r^2+R^2-2rR\eta)^{3/2}} &\text{where, }0\leq r \leq \infty \text{ & }-1 \leq \eta \leq 1\end{align}
When one plots $f$ as a function of $r$ for various values of $\eta$, one observes that $f$ is continuous at $r=R$ for all values of $\eta$ except $\eta=1$. In the case of $\eta=1$, $f$ diverges to $+\infty$ and $-\infty$ on the left and right sides of $r=R$ respectively $\left(\because f(r,1)=-\frac{r-R}{|r-R|^3}\right)$.
This implies the following,
\begin{align}g(\eta) \equiv \lim_{r \to R+}f(r,\eta)-\lim_{r \to R-}f(r,\eta) \; &\text{is zero for }\eta \neq 1 \\ & \text{ blows up for }\eta=1 \end{align}
This is similar to how a Dirac delta function behaves (blows up at one point and zero everywhere else). A stronger motivation for why I believe it might be a Dirac delta function is given in the next section.
Question: Is $g(\eta)$ as defined above, a Dirac delta function in $\eta$ (up to some scale factor)?

Section 2: Physics Problem
The physics problem setup is a general spherical surface charge distribution $\sigma(\theta,\phi)$ of radius $R$.
It is known that the component of the electric field, $\mathbf{E}=-\nabla\Phi$, that is normal to the spherical surface is discontinuous. i.e.,
$$\lim_{r \to R+}\partial_r \Phi(r,\theta,\phi)-\lim_{r \to R-}\partial_r\Phi(r,\theta,\phi)=-\frac{\sigma(\theta,\phi)}{\epsilon_0} \tag{1; eq. 2.31 in [1]}$$
The above result is commonly proved by applying Gauss' law to an infinitesimal Gaussian "pill-box" covering the region of interest.
However, I wish to prove the above result (eq. 1) by only using the following Green's function solution for the electric potential (eq. 2).
\begin{align}&\Phi(\mathbf{r}) =\frac{1}{4\pi \epsilon_0}\int \frac{\rho(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}d^3\mathbf{r}' &\rho(\mathbf{r})=\sigma(\theta,\phi)\delta(r-R) \tag{2}\\ 
\Rightarrow \;&\Phi(\mathbf{r})=\frac{1}{4\pi \epsilon_0}\int\frac{\sigma(\theta',\phi')}{|r \hat{r}-R\hat{r}'|}R^2\sin\theta' d\theta' d\phi' &\text{where, }\hat{r}=\hat{r}(\theta,\phi)  \text{ & }\hat{r}'=\hat{r}(\theta\,',\phi') \end{align}
Using $|r \hat{r}-R\hat{r}'|=\sqrt{r^2+R^2-2rR\hat{r}\cdot\hat{r}'}$, we have,
\begin{align}\partial_r \Phi=-\frac{1}{4\pi \epsilon_0}\int\sigma(\theta',\phi')\frac{r-R\hat{r}\cdot\hat{r}'}{|r \hat{r}-R\hat{r}'|^3}R^2\sin\theta' d\theta' d\phi'  \end{align}
I encountered a question in the course of my attempt to prove eq. 1. I'll describe it below.
\begin{align}&\partial_r \Phi=-\frac{1}{4\pi \epsilon_0}\int\sigma(\theta',\phi')\frac{r-R\eta}{(r^2+R^2-2rR\eta)^{3/2}}R^2\sin\theta' d\theta' d\phi'  &\text{where, }\eta \equiv \hat{r}\cdot\hat{r}' \end{align}
$$\lim_{r \to R+}\partial_r \Phi-\lim_{r \to R-}\partial_r\Phi =\frac{1}{4\pi \epsilon_0}\int\sigma(\theta',\phi')(\lim_{r \to R+}f-\lim_{r \to R-}f)R^2\sin\theta' d\theta' d\phi'\tag{3}$$
$$\text{where, }f(r,\eta)\equiv -\frac{r-R\eta}{(r^2+R^2-2rR\eta)^{3/2}} $$
When one plots this function $f$ online as a function of $r$ for various values of $\eta$, one observes that $f$ is continuous at $r=R$ for all values of $\eta$ ($\eta \in [-1,1]$) except $\eta=1$. For $\eta=1$, the function $f$ diverges to $+ \infty$ and $- \infty$ on the left and right  sides of $r=R$ respectively $\left(\because f(r,1)=-\frac{r-R}{|r-R|^3}\right)$.
This implies the following,
\begin{align}g(\eta) \equiv \lim_{r \to R+}f(r,\eta)-\lim_{r \to R-}f(r,\eta) \; &\text{is zero for }\eta \neq 1 \tag{4}\\ & \text{ blows up for }\eta=1 \text{ ($\eta=1$ $\Leftrightarrow$ $\theta'=\theta$ and $\phi'=\phi$)}\end{align}
This looks promising because the above behavior is similar to a Dirac delta function (blows up at one point and zero everywhere else). The discontinuity in the electric field at $(\theta,\phi)$ is only "aware" of the value of the surface charge density $\sigma$ at $(\theta,\phi)$ (ref. eq. 1) and hence, I believe I need a Dirac delta function in the integral in eq. 3 to get the $\sigma$ out of the integral.
Question: Is $g(\eta)$ as defined in eq. 4, a Dirac delta function (up to some scale factor $\#$)? That is,
$$\text{Is }g= (\#)\; \delta(\theta'-\theta)\delta(\phi'-\phi)?$$
I'd really appreciate any insight that addresses my problem.
References
$[1]$ Griffiths, Introduction to Electrodynamics (3rd ed.)
 A: This is an interesting question, and the answer to it is definitely YES. Let us deal with the problem using expansions in spherical harmonics.
First note that the expansion of the Laplacian's Green's function in terms of Legendre polynomials is known
$$\frac{1}{|\mathbf{r}-\mathbf{r}'|}=\begin{Bmatrix}\frac{1}{r}\sum_{n=0}^{\infty}P_n(\hat{r}\cdot\hat{r}')\Big(\frac{R}{r}\Big)^{n}~~~~,r> R\\ \frac{1}{R}\sum_{n=0}^{\infty}P_n(\hat{r}\cdot\hat{r}')\Big(\frac{r}{R}\Big)^{n}~~~~, r<R\end{Bmatrix}$$
Then we compute the derivative of this function with respect to the radial coordinate $r$ which will help us compute the electric field:
$$f(r,\eta)=\frac{\partial}{\partial r}\frac{1}{|\mathbf{r}-\mathbf{r}'|}=\begin{Bmatrix}-\frac{1}{R^2}\sum_{n=0}^{\infty}(n+1)P_n(\hat{r}\cdot\hat{r}')\Big(\frac{R}{r}\Big)^{n+2}&,~r> R\\ \frac{1}{R^2}\sum_{n=0}^{\infty}nP_n(\hat{r}\cdot\hat{r}')\Big(\frac{r}{R}\Big)^{n-1}&,~r<R\end{Bmatrix}$$
Finally we deal with the quantity claimed to be a delta-function:
$$g(\eta)=\frac{\partial}{\partial r}\frac{1}{|\mathbf{r}-\mathbf{r}'|}\Bigg|_{r\to R^+}-\frac{\partial}{\partial r}\frac{1}{|\mathbf{r}-\mathbf{r}'|}\Bigg|_{r\to R^-}=-\frac{1}{R^2}\sum_{n=0}^{\infty}(2n+1)P_n(\hat{r}\cdot\hat{r}')$$
Using the addition theorem for spherical harmonics, which states that
$$P_{n}(\hat{\mathbf{r}}\cdot\hat{\mathbf{r}}')=\frac{4\pi}{2n+1}\sum_{m=-n}^{n}Y_{nm}(\hat{\mathbf{r}})(Y_{nm})^*(\hat{\mathbf{r}}')$$
and the expansion of the delta function in spherical harmonics:
$$\delta(\theta-\theta')\delta(\phi-\phi')=\sin\theta'\sum_{lm}Y_{lm}(\theta, \phi)Y^*{}_{lm}(\theta', \phi')$$
we have proven that
$$\frac{\partial}{\partial r}\frac{1}{|\mathbf{r}-\mathbf{r}'|}\Bigg|_{r\to R^+}-\frac{\partial}{\partial r}\frac{1}{|\mathbf{r}-\mathbf{r}'|}\Bigg|_{r\to R^-}=-\frac{4\pi}{R^2}\frac{\delta(\theta-\theta')\delta(\phi-\phi')}{\sin\theta'}$$
and thus 
$$\begin{align}E_r(r\to R^{+},\theta, \phi)-E_r(r\to R^{-},\theta, \phi)&=-\frac{1}{4\pi\epsilon_0}\int R^2\sin\theta'd\theta'd\phi'\sigma(\theta', \phi')\Big(-\frac{4\pi}{R^2}\frac{\delta(\theta-\theta')\delta(\phi-\phi')}{\sin\theta'}\Big)\\&=\frac{\sigma(\theta,\phi)}{\epsilon_0}\end{align}$$
and our faith to mathematical consistency of physical theories has been successfully restored.
EDIT: Upon @mrc ntn's prompt,and for completeness, I want to close this by showing that $g(\eta)\propto\delta(\eta-1)$.
Without repeating any of the calculations above because they are identical (substituting $\eta=\hat{r}\cdot\hat{r}'$, it is true that
$$g(\eta)=-\frac{1}{R^2}\sum_{n=0}^{\infty}(2n+1)P_n(\eta)$$
But $P_n(1)=1$ and due to the completeness relation of the Legendre polynomials given here we conclude that
$$g(\eta)=-\frac{1}{R^2}\sum_{n=0}^{\infty}(2n+1)P_n(\eta)P_n(1)=-\frac{2}{R^2}\delta(\eta-1)$$
A: [Note: I'm still not happy with this.  However, I have to take a break and come back to reconsider it.]
Partial (and tentative) answer.  I think the answer may be yes.  I think that we can assume $R = 1$ without loss of generality.  (OP will please correct me if this is not the case.)
In that case, we have
$$
f(r, \eta) = \frac{r-\eta}{(r^2+1-2r\eta)^{3/2}}
$$
and then
$$
g(\eta) = \lim_{r\to1^+} f(r, \eta) - \lim_{r\to1^-} f(r, \eta)
$$
As I indicate in the comments, I would consider this a Dirac delta function (up to scaling) if we can reasonably interpret
$$
\int_{\eta = -1}^1 g(\eta) h(\eta) \, d\eta = k h(1)
$$
for arbitrary $h(\cdot)$ defined on $[-1, 1]$, where $k$ is a constant that doesn't depend on $h(\cdot)$.  As a plausibility argument, let $h(\eta) = 1$ identically and integrate
$$
F(r, \eta) \stackrel{\text{def}}{=} \int f(r, \eta) \, d\eta
    = \frac{1-r\eta}{r^2\sqrt{r^2-2r\eta+1}} + C
$$
for $r \geq 0, -1 \leq \eta \leq 1$ and $(r, \eta) \not= (1, 1)$.  We now write
$$
G(\zeta) \stackrel{\text{def}}{=} \int_{\eta=-1}^\zeta g(\eta) \, d\eta
$$
where we might reasonably interpret
$$
G(\zeta)
    = \left[\lim_{r\to1^+} F(r, \zeta) - \lim_{r\to1^-} F(r, \zeta)\right]
    - \left[\lim_{r\to1^+} F(r, -1) - \lim_{r\to1^-} F(r, -1)\right]
$$
For $\zeta < 1$, we find that the two-sided limit exists and is equal to
$$
\lim_{r\to1} F(r, \zeta) = \sqrt{\frac{1-\zeta}{2}}
$$
and so $G(\zeta) = 0$ (giving us a sense that $g(\eta)$ doesn't "pick off" any values for $\eta < 1$), but when $\zeta = 1$, we have
$$
\lim_{r\to1^+} F(r, 1) = -1
$$
but
$$
\lim_{r\to1^-} F(r, 1) = 1
$$
yielding $G(1) = -2$.  So it looks to me like $g(\eta)$ is $-2$ times a Dirac delta function.  I'm not prepared at the moment to follow this to its logical conclusion for general $h(\cdot)$, but perhaps this will give you an idea of how to proceed.
A: The short answer to your quetion is: yes.  Your function $f(r,\eta)$ is in fact the Dirac delta measure when $\eta = 1$.  The broader explanation for this is given by the divergence of the expression $\frac{\textbf{r-r'}}{||\textbf{r-r'}||^3}$.  Observe that in computing the divergence, we expect that $\textbf{r} \neq \textbf{r'}$ (otherwise the quantity is undefined):
\begin{eqnarray}
\nabla \cdot \frac{\textbf{r - r'}}{||\textbf{r-r'}||^3} & = & \nabla \cdot \frac{(x-x', y-y', z-z')}{\left ((x-x')^2 + (y-y')^2 + (z-z')^2\right )^{3/2}} \\
& = & \sum_{k=1}^3\frac{||r-r'||^3 - 3(x_k-x_k')^2||r-r'||}{||r-r'||^6} \\
& = & \frac{3||\textbf{r-r'}||^3 - 3||\textbf{r-r'}||^3}{||\textbf{r-r'}||^6}\\
& = & 0.
\end{eqnarray}
However, taking $\textbf{r}$ to be the variable of integration we can verify via the divergence theorem that when we integrate over a sphere of any radius centered at $\textbf{r'}$ that we get a nonzero quantity.  Consider the set of $\textbf{r} \in \mathbb{R}^3$ such that $||\textbf{r-r'}|| = \epsilon$.  Then we have by the divergence theorem
\begin{eqnarray*}
\iiint \nabla \cdot \frac{\textbf{r-r'}}{||\textbf{r-r'}||^3}d\textbf{r}^3 &  = & \iint \frac{\textbf{r-r'}}{||\textbf{r-r'}||^3}\cdot d\textbf{S} \\
& = & \int_0^\pi\int_0^{2\pi} \epsilon^2\sin\theta \frac{\textbf{r-r'}}{||\textbf{r-r'}||^3}\cdot \frac{\textbf{r-r'}}{||\textbf{r-r'}||}d\phi d\theta \\
& = & 4\pi\epsilon^2\frac{||\textbf{r-r'}||^2}{||\textbf{r-r'}||^4} \\
& = & 4\pi.
\end{eqnarray*}
Despite the divergence being almost everywhere zero, we can verify that it has a nonzero integral along a sphere of arbitrary radius around the point $\textbf{r'}$.  We therefore conclude that 
$$
\nabla \cdot \frac{\textbf{r - r'}}{||\textbf{r-r'}||^3} \;\; =\;\; 4\pi \delta\left (\textbf{r-r'}\right ).
$$
Note that this matches the physical intuition we would expect if a charge distribution was placed on a spherical shell.  The result also matches mathematical intuition as well: The reason why we don't run across this issue as much with point charges is because Riemann integration (consistent with Lebesgue integration for smooth functions) is not sensitive to point discontinuities (i.e. the endpoints of integration do not contribute to the final integral in the sense that a single-variable integral over the interval $(a,b)$ is no different than that over $[a,b), \; (a,b]$, or $[a,b]$).  However when a distribution to occupies a larger portion of physical space, then we need to alter the measure of integration to reflect this.  The Dirac delta arises precisely because the charge distribution takes up a region of nonzero measure in the 2-dimensional $\theta$-$\phi$ space, but a spherical shell has zero measure in 3-dimensional $r$-$\theta$-$\phi$ space.
A: It is not. Define $g(\eta , a)= f( R + a, \eta)  -  f(R - a, \eta)$ for $a>0$ and small. Expand in $a$, 
$$
g(\eta , a) = \frac{( 3 \eta -1) a}{ 2 \sqrt{2 - 2 \eta } (\eta-1) R^3} +O(a^3)
$$ 
To have a Dirac delta you need something that peaks at some point while shrinking for $a \rightarrow 0^+$ (i.e. loosely speaking  $\delta(\eta) \approx a^{-1}e^{-\eta^2/a^2}$ apart from a normalization factor). 
In your case $g(\eta , a)$ does not shrink , it just goes to zero. 
So it is not a Delta, it is (from the philosophical point of view) just something like $a/f(\eta)$ where $f(\eta)$ has some zero for some $\eta$ in the domain. 
Edit: try to plot the full $g(\eta,a)$, fixing a certain value of $R$ (or its linear expansion in $a$, they are clearly the same function for small $a$). Make several plots for $a=0.1, 0.01..$. Maybe not what a "true mathematician" wants to do, but you will convince yourself that there is no way it can provide a good representation of the Delta. You can also do some numerical experiment integrating a test function for smaller and smaller values of $a$. Still the same conclusion.
