MODEL UI2PQ3
COMMENT-------------------------------------------------------------------------
Calculates the input active and reactive power of time domain voltage/current 
signals by a user selectable algorithm. Fixed frequency.
Output in 3-phase VA quantities. Downsampling approach ala COMTRADE component.
Algorithm=0: 8-point RADIX-2 FFT algorithm.
Algorithm=1: Recursive DFT algorithm. HKH Apr. 2021.
Algorithm=2: alpha-beta transformation by Peņaloza
----------------------------------------------------------------------ENDCOMMENT
INPUT 
  V[1..3], I[1..3]     --input signals to be transformed to phasors and power
DATA  
  FREQ {DFLT:50}       --power frequency
  SampleFreq {dflt:400}--sample frequency in this model in Hz
  ScaleV {dflt:1}      --multiply output by this number 
  ScaleI {dflt:1}      --multiply output by this number
  Algorithm {dflt:0}   --0: FFT radix2-8, 1: DFT recursive
OUTPUT 
  P, Q     --Active and reactive power output signals (rms)
VAR    
  P,Q,reV[1..3],imV[1..3],reI[1..3],imI[1..3],OMEGA,D,alpha,x1,x3,x5,x7,
  Ti,dt,DTi, xre,xim,n,delta_T[0..7],
  P3A,Q3A,P3I,Q3I,DP3,DQ3,VAT,IAT,VBT,IBT,Period
CONST SQR3 {VAL:1.732050807569}
HISTORY
  V[1..3] {DFLT:0}
  I[1..3] {DFLT:0}
  P3A {DFLT: 0.0}, Q3A {DFLT: 0.0}
  P3I {DFLT: 0.0}, Q3I {DFLT: 0.0}
DELAY CELLS (V[1..3]) : recip(FREQ*timestep)+2
DELAY CELLS (I[1..3]) : recip(FREQ*timestep)+2
DELAY CELLS (P3I, Q3I): recip(FREQ*timestep)+2

INIT 
  OMEGA:= 2*PI*FREQ
  reV[1..3] :=0
  imV[1..3] :=0
  reI[1..3] :=0
  imI[1..3] :=0
  alpha:=1/sqrt(2)
  for n:= 0 to 7 do
    delta_T[n] := n/(FREQ*8)
  endfor
  Period:=recip(FREQ)
  dt:=recip(SampleFreq)
  Ti:=0
ENDINIT
EXEC
  DTi:=T-Ti
  if DTi>=0
  then
    if Algorithm=2
    then
      VAT := (2.0*V[1] - V[2] - V[3])/SQR3
      IAT := (2.0*I[1] - I[2] - I[3])/SQR3
      VBT := V[2] - V[3]
      IBT := I[2] - I[3]
      P3I := (VAT*IAT + VBT*IBT)/2.0*ScaleV*ScaleI
      Q3I := (VBT*IAT - VAT*IBT)/2.0*ScaleV*ScaleI
      P3A := P3A + (delay(P3I, DTi, 2) - delay(P3I, Period+DTi, 2))*dt*FREQ
      Q3A := Q3A + (delay(Q3I, DTi, 2) - delay(Q3I, Period+DTi, 2))*dt*FREQ
      P := P3I*(Ti<=Period) + P3A*(Ti>Period)
      Q := Q3I*(Ti<=Period) + Q3A*(Ti>Period)
    else
      P:=0
      Q:=0
      for n:=1 to 3 do
        if Algorithm=0
        then
          x1 := delay(V[n],delta_T[0]+DTi,2) - delay(V[n],delta_T[4]+DTi,2)
          x3 := delay(V[n],delta_T[2]+DTi,2) - delay(V[n],delta_T[6]+DTi,2)
          x5 := delay(V[n],delta_T[1]+DTi,2) - delay(V[n],delta_T[5]+DTi,2)
          x7 := delay(V[n],delta_T[3]+DTi,2) - delay(V[n],delta_T[7]+DTi,2)
          xRe := x1 + (x5 - x7)*alpha
          xIm := x3 + (x5 + x7)*alpha
          ReV[n] := (xRe*cos(OMEGA*Ti)+xIm*sin(OMEGA*Ti))/4
          ImV[n] := (xIm*cos(OMEGA*Ti)-xRe*sin(OMEGA*Ti))/4
          x1 := delay(I[n],delta_T[0]+DTi,2) - delay(I[n],delta_T[4]+DTi,2)
          x3 := delay(I[n],delta_T[2]+DTi,2) - delay(I[n],delta_T[6]+DTi,2)
          x5 := delay(I[n],delta_T[1]+DTi,2) - delay(I[n],delta_T[5]+DTi,2)
          x7 := delay(I[n],delta_T[3]+DTi,2) - delay(I[n],delta_T[7]+DTi,2)
          xRe := x1 + (x5 - x7)*alpha
          xIm := x3 + (x5 + x7)*alpha
          ReI[n] := (xRe*cos(OMEGA*Ti)+xIm*sin(OMEGA*Ti))/4
          ImI[n] := (xIm*cos(OMEGA*Ti)-xRe*sin(OMEGA*Ti))/4
        else
          D:=2*FREQ*dt*(delay(V[n],DTi,2)-delay(V[n],Period+DTi,2)) 
          ReV[n]:=ReV[n]+D*cos(OMEGA*Ti)
          ImV[n]:=ImV[n]-D*sin(OMEGA*Ti)
          D:=2*FREQ*dt*(delay(I[n],DTi,2)-delay(I[n],Period+DTi,2))
          ReI[n]:=ReI[n]+D*cos(OMEGA*Ti)
          ImI[n]:=ImI[n]-D*sin(OMEGA*Ti)
        endif
        P:=P+(reV[n]*reI[n]+imV[n]*imI[n])*ScaleV*ScaleI/2
        Q:=Q+(imV[n]*reI[n]-reV[n]*imI[n])*ScaleV*ScaleI/2
      endfor
    endif
    Ti:=Ti+dt
  endif  
ENDEXEC
ENDMODEL