The reason you're having trouble is because you don't know what reactive power really is. It is the quantity of real power that flows from source to load that has an equivalent quantity of real power that flows from load to source, thus producing no work. It's just the quantity power circulating back and forth doing no work and not getting dissipated. But it is all just real power.
See illustrated example here:
How does power factor show itself in this data from sensor readings?
Knowing that, you should realize you no longer need to know anything but V and I. No RMS, no \$\sqrt 3\$, no wye vs. delta, no balanced vs unbalanced, no lin-to-neutral vs lin-to-line. You just do everything numerically:
Each instantaneous voltage and instantaneous current sample taken at the same time are multiplied together to produce an instantaneous power point and a power waveform is constructed.
The time average of the power waveform is the real power (that is, the net power, since it summed up all positive and negative power points) that is flowing from source-to-load, or load-to-source, depending in the sign of the result. This is performed on-the-fly by numerically integrating the power and dividing by the integration interval.
Reactive power is the smaller value between the time average of all positive power points and time average of negative power points. Why? The smaller value gets swamped out by the larger value cancelling out which is reactive power.
Apparent power is a bit meaningless when you have all this data so
if you want it you jut construct it using Pythagoras on the real and
reactive powers calculated. Or you can use numerically integration
to calculate RMS of the voltage and current waveforms and multiply
the values together.
$$\text{RMS}=\sqrt{\frac{1}{T}\int_0^T{x^2 dt}}$$ (If you did not
know, this is how you calculate RMS for anything. This is what you should have been using to get RMS values from real world measurements of arbitrary signals since only uses measurements and makes no assumptions)
You can then calculate PF.
You don't need either apparent power or PF if you actually know the real power running in both directions. They are just abstractions that exist to model the values you get when you take voltage and current measurements independently thus losing all timing/phase information.
As for the integration interval, you can try to make it one period if you know it or can detect it, but if you don't/can't then you can make it a time much longer than one period (a hundred or a thousand times). The integral sum of all the whole cycles makes the contribution of the last cycle infinitesmal so it does not matter where the last cycle gets truncated due to timing imprecision.
Regardless of whether you have a wye or delta, if you have instrumentation measuring the phase voltage (not necessarily the line-to-line voltage specially for the wye) and if you are measuring the phase currents (not necessarily the same as line currents specially for the delta) then you can simply use the method above as is to calculate the instantaneous power of each phase and sum them. That will give you the real power of the three-phase system.
However, if you want to treat the load is unknown things get a bit trickier since now you only have measurements taken from outside the load: line-to-line voltages and line currents. That means that the phase voltage is obscured for wyes and the phase currents are obscured for deltas. It is not immediately obvious how to separate line currents into phase currents for a delta, nor is it obvious how to separate line-to-line voltages into phase voltages for a wye. Things like nodal analysis don't work and give you infinite solutions.
In this case, use Blondel's thereom:
\$p(t) = v_{ab}(t)i_a(t) + v_{cb}(t)i_c(t)\$
This is basically a clever way to use superposition since it expands out to be:
\$p(t) = v_a(t)v_a(t) + v_b(t)v_b(t) + v_c(t)v_c(t)\$
You can verify yourself by modelling it as a wye and using:
\$i_a(t)+i_b(t)+i_c(t) = 0\$
\$v_{\phi}(t) = v_{\phi} - v_{{\phi}n}\$, where '\$\phi\$' is phase A, B, or C, and 'n' is the neutral
\$v_{x}(t) = v_{y} - v_{x-y}\$, where 'x' and 'y' are Phases A, B, or C
It's about five simple lines long and if you think about it as a wye, it basically calculates the power of the phases A and B (where currents are being measured through normally). However, in Phase C, where the measured currents merge, the power here is calculated through superposition of powers via superposition of the currents without actually needing to separate the line-to-line voltages of the wye into distinct phase voltages. However because of this, Blondel's Thereom only tells you overall real power and not the power in individual phases.
Beyond that, for power such as in individual phases or power factors, etc you need PQ Theory which I have only just become aware of with complications trying to answer this question and it is a lot more complicated and I am only just reading on it now.
Though, this equation does exist for imaginary power outside of PQ Theory but I only became aware of it reading about PQ theory:
\$q(t) = \frac{1}{\sqrt(3)}[v_{ab}(t)i_c(t) + v_{bc}(t)i_a(t) +v_{ca}(t)i_b(t)]\$
