NIntegrate issue without any singularity

I'm evaluating an integral using NIntegrate, I checked the expression and there shouldn't be any singularity at least on paper. The parameter zs can take on values in [0,9.9] and zh=10. I took zs=9 as an example, but when I evaluate it gave me an error about singularity/WorkingPrecision. I have no idea why this is happening since I already tried to adjust WorkingPrecison and other ways. Any help?

d = 3;
ag = 8;
pg = 8;
wp = 30;
f[z_, zh_] := 1 - (z/zh)^(d + 1);
Sinteg[z_?NumericQ, zs_?NumericQ, zh_?NumericQ] := Module[{zr, zsr, zhr}, {zr, zsr, zhr} = Rationalize[{z, zs, zh}, 0]; SetPrecision[(Sqrt[1/(f[zr, zhr] (1 - (zr/zsr)^(2 d) (1 + ((2 d f[zsr, zhr])/(d + 1)) (zhr/zsr)^(d + 1) (1 - (f[zsr, zhr]/f[zr, zhr])))))]), wp]]
SFull[zs_?NumericQ, zh_?NumericQ] := Module[{zsr, zhr}, {zsr, zhr} = Rationalize[{zs, zh}, 0]; NIntegrate[SetPrecision[Sinteg[z, zsr, zhr], wp], {z, 0, zsr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 300]]

SFull[9,10]
(*During evaluation of In[214]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
During evaluation of In[214]:= NIntegrate::zeroregion: Integration region {{9.00000000000000000000000000000,9.0000000000000000000000000000000000000000000000000000000000000000000000000000000}} cannot be further subdivided at the specified working precision. NIntegrate assumes zero integral there and on any further indivisible regions.
330.387118746134454732010635945*)
$$$$

• Your integrand diverges as $1/(z_s-z)$ as $z\to z_s$, and so the integral is unbounded. More precisely, for $z\to z_s$ the integrand is $\frac{z_h^2 z_s}{\sqrt{6} (z_s-z) \sqrt{5 z_h^4-z_s^4}} + O(1)$. May 30, 2022 at 10:48
• @Roman The original expression is under a square root sign so it should be finite right under an integral? As $z \rightarrow z_s$ the denominator approaches $\sqrt{0}$. May 30, 2022 at 10:59
• @Roman Did you use Mathematica to expand out the integrated as $z \rightarrow zs$? If yes, how did you do it? May 30, 2022 at 15:29
• See answer below, too long for a comment. May 30, 2022 at 15:46

As per comment, the integrand diverges for $$z\to z_s$$:

d = 3;
f[z_, zh_] = 1 - (z/zh)^(d + 1);
Sinteg[z_, zs_, zh_] = Assuming[0 < z < zs < zh,
Sqrt[1/(f[z, zh] (1 - (z/zs)^(2 d) (1 + ((2 d f[zs, zh])/(d + 1))
(zh/zs)^(d + 1) (1 - (f[zs, zh]/f[z, zh])))))] // FullSimplify]


$$\frac{\sqrt{2} z_h^2 z_s^5}{\sqrt{3 z^{10} z_h^4+3 z^6 z_s^8+2 z_s^{10} \left(z_h^4-z^4\right)-z^6 z_s^4 \left(z^4+5 z_h^4\right)}}$$

Assuming[0 < z < zs < zh,
Series[Sinteg[z, zs, zh], {z, zs, 0}] // FullSimplify]
`

$$-\frac{z_h^2 z_s}{(z-z_s) \sqrt{30 z_h^4-6 z_s^4}}+\frac{65 z_h^6-17 z_h^2 z_s^4}{6 \sqrt{6} \left(5 z_h^4-z_s^4\right)^{3/2}}+O(z-z_s)$$