Night-Day Cycles in S-Type Planetary Systems
Although this question would probably be more suited for the worldbuilding.se, you're asking for scientific references, so I guess this is ok to ask here.

What you describe is a hierarchical system, akin to Sun-Jupiter-Anymoon. I assume the planet orbits star $S_1$ directly, and $S_2$ is the larger star, and the planet will be abbreviated with $P$.  

I don't know any references, but just from geometric considerations I'll state a formula here that should work in a quasi-2D setting (what I mean by this, later).
The surface illumination function $I(t)$ should be a simple sum of the two-star contributions $$I(t)=I_1\,\sin(\omega_0 t) + I_2(t)\, \sin(\omega_1\,t), $$ where $\omega_0$ is the spin rotation frequency of P around itself. The amplitudes of the sinii are $I_1$, the constant amplitude related to the constant distance of $S_1$ to $P$, and $I_2(t)$ being the irradiation amplitude coming $S_2$. As example, $I_1$ could be a number like the familiar $1360 \, kW/m^2$ for the sun.
I wrote $I_2(t)$ as a function of time, because distance from the planet to $S_2$ varies over time, as it orbits $S_1$. Therefore, as simplest approximation, I would assume that this is also a sinus, hence $$I_2(t)=I_{2,flct}\; \sin(\omega_1\,t)+I_{2,avg},$$ with $I_{2,avg}$ being the average irradiation at $S_1$'s orbit around $S_2$, and $I_{2,flct}$ is given by the minimum and maximum distance of the planet towards $S_2$, and finally $\omega_1$ being the orbital frequency of P around $S_1$.

Once you have some numbers for the constants, you can use python, or any other routine, in order to investigate the behaviour of your day/night cycles. Here it is important that your routine can replace the negative contributions in $I(t)$ with zero values, otherwise the negative values of $I_1$ could 'eat' the positive values of $I_2(t)$ and instead of perpetual day you'd get a perpetual night from the formula, when P would be directly between $S_1$ and $S_2$. In python you would do this via some command like "I_updated = [0 if i &lt; 0 else i for i in I_raw]", where "I_raw" would be the values from the individual $S_1$ or $S_2$ contributions, before summing them up to $I(t)$.

If you want to know 'which sun is up' at any given moment, you can plot the sinii separately. 

This quickly scribbled formula ignores possible 3D-effects on the illumination, such as a possible Kozai resonance in your system, but also misrepresents 'shadowing' of $S_2$ by $S_1$. Therefore the formula represents a very slightly inclined orbital plane w.r.t $S_1$.

Please observe that your system, as you want to build it, might not be stable for too long:
For the hierarchical structure of your orbits to exist, the stars must have a mass ratio of at least $m_2/m_1\approx 23$. If $m_1$ is a sun-like G0 star, then $m_2$ will be an O9, which has a lifetime of only 4 million years on the main sequence.
However, if you assign $m_2$ to be the G0 star, then $m_1$ could be a M2 dwarf, and that should work out fine. Also remember, that in your chosen configuration, there's always one part of the planets orbit around $S_1$, namely when it's exactly between $S_1$ and $S_2$, where the whole planet experiences perpetual day, and depending on the choice of stellar type, might get quite cooked.

For more on three-body motion, I'd recommend diving into Murray &amp; Dermott: Solar system dynamics.