####### Written by Doug Millner of NERX Power Consultants LLC.  
#######  You are allowed to use freely provided you never remove the above credits 
######   for the code. 

##The code is written to only take in phase values.  If you have sequence values
##convert them with the A matrix and hand punch it in.   To use code, put in the values for the 
##magnitude and angle of the phasors.  


############ Set GenerateGif to 1 if you wish to not display a figure and dump out gif files to animate. 
############ Set GenerateGif to 0 to no generate gifs and display plots. 
GenerateGif = 0;
##########  ENTER PHASOR MAGNITUDE AND ANGLES HERE

%magnitude of phase phasors

VanMag = 1
VbnMag = 1
VcnMag = 1
%angle of phase to neutral phasors
VanAngle = 0; 
VbnAngle = 240;
VcnAngle = 120;
######################################################################3

#degrees per frame
FrameDegree = 2;
## total number of frames
totalframes = 180

#set to 1 for axis for phasors, 0 for no.
AxisForPhasors = 1


###### end of inputs

#######################################################################################################




%uncomment next line to when making gif to prevent a plot from being made during the creation of gifs.

set(0, 'defaultfigurevisible', 'on');
if GenerateGif == 1
  
set(0, 'defaultfigurevisible', 'off');


endif

%degrees per frame. 


a = exp(i*(2/3)*pi)

A = [1 1 1; 1 a^2 a; 1 a a^2]

%phasors
Van = VanMag*(real(cos(VanAngle*pi/180)) + i*(sin(VanAngle*pi/180)));
Vbn = VbnMag*(real(cos(VbnAngle*pi/180)) + i*(sin(VbnAngle*pi/180)));
Vcn = VcnMag*(real(cos(VcnAngle*pi/180)) + i*(sin(VcnAngle*pi/180)));

VAN = Van
VBN = Vbn
VCN = Vcn


Sequence = (A^-1)*[Van;Vbn;Vcn]
Va0 = Sequence(1,1)
Va1 = Sequence(2,1)
Va2 = Sequence(3,1)




m = 2 * max([real([Van Vbn Vcn]) imag([Van Vbn Vcn])]);
scale = max([Va0 Va1 Va2])

% Calculate a suitable axis limit based on phasors


Va0scale = (m/m)*(Va0/scale)
Va1scale = (m/m)*(Va1/scale)
Va2scale = (m/m)*(Va2/scale) 



 figure 
 hold on

   plot1 = plot([-4*m 4*m],[0 0])
  plot2 = plot([0 0],[-4*m 4*m])

  if AxisForPhasors == 1
  plot([-1.5*m -1.5*m],[-2*m 0*m])
  plot([1.5*m 1.5*m],[0*m -2*m])
 endif
 
 text1 = text (-1.5*m, m*.5, "Positive Phasors");
text2 = text(0.1*m, -.6*m, "Phase Phasors")
text3 = text (0.1*m, .5*m, "Zero Phasors")
text4 = text(1.9*m, m*.5, "Negative Phasors")
text5 = text(-1.5*m, -1.7*m, "Positive Values")
text6 = text(0.1*m, -1.7*m, "Zero Values")
text7 = text(1.9*m, -1.7*m, "Negative Values")
text8 = text(1.9*m, m*-.6, "Phase Values")
text9 = text(-2*m,2.1*m, "Phase A Blue, Phase B Green, Phase C Red, Only Phase Values to Scale")

set(text1, 'fontsize', 20)
set(text2, 'fontsize', 20)
set(text3, 'fontsize', 20)
set(text4, 'fontsize', 20)
set(text5, 'fontsize', 20)
set(text6, 'fontsize', 20)
set(text7, 'fontsize', 20)
set(text8, 'fontsize', 20)
set(text9, 'fontsize', 8)







%this is how many frames are generated
for t=[0:totalframes]
 
 
 
 


Sequence = (A^-1)*[Van;Vbn;Vcn]
Va0 = Sequence(1,1)
Va1 = Sequence(2,1)
Va2 = Sequence(3,1)
 
 
 
Va0scale = (m/m)*(Va0/scale)
Va1scale = (m/m)*(Va1/scale)
Va2scale = (m/m)*(Va2/scale) 


 
 
 
 
%place any varying function of the phasors for the animations here and frame to totalframes
%Van = Van+(VCN*(-1)/totalframes)
%Vbn = Vbn+(VCN*(-1)/totalframes)
%Vcn = Vcn+(VCN*(-1)/totalframes)
 
 
 
 
 
 
 
 
 
 
 
 
 
 
  rotationconstant = exp(i*(FrameDegree*pi/180));
  Vanrotated = Van*rotationconstant^t;
  Vbnrotated =  Vbn*rotationconstant^t;
  Vcnrotated = Vcn*rotationconstant^t;
  
  Va0rotated = Va0scale*rotationconstant^t;
  Va1rotated = Va1scale*rotationconstant^t;
  Va2rotated = Va2scale*rotationconstant^t;
  
  
  
  
  
  
  

  
  
  
  
  
 axis([-2*m,2*m,-2*m,2*m], 'square')
 
  %VanPlot.Color = 'red';
  %VbnPlot.Color = 'blue';
  %VcnPlot.Color =  'green';
% Plot phasors one at a time

  VanPlot = quiver(0,0,real(Vanrotated),imag(Vanrotated));
  VbnPlot = quiver(0,0,real(Vbnrotated),imag(Vbnrotated));
  VcnPlot = quiver(0,0,real(Vcnrotated),imag(Vcnrotated));
  
  Va0Plot = quiver(0, 1*m,real(Va0rotated), imag(Va0rotated));
  Vb0Plot = quiver(0, 1*m,real(Va0rotated), imag(Va0rotated));
  Vc0Plot = quiver(0, 1*m,real(Va0rotated), imag(Va0rotated));
  
  
  Va1Plot = quiver(-1.5*m, 1*m,real(Va1rotated), imag(Va1rotated));
  Vb1Plot = quiver(-1.5*m, 1*m,real(Va1rotated*a^2), imag(Va1rotated*a^2));
  Vc1Plot = quiver(-1.5*m, 1*m,real(Va1rotated*a), imag(Va1rotated*a));
  
  Va2Plot = quiver(1.5*m, 1*m,real(Va2rotated), imag(Va2rotated));
  Vb2Plot = quiver(1.5*m, 1*m,real(Va2rotated*a), imag(Va2rotated*a));
  Vc2Plot = quiver(1.5*m, 1*m,real(Va2rotated*a^2), imag(Va2rotated*a^2));
  axis1 = plot(-1.5*m -1.5*m,-2*m,-2*m)
  
  
  
  Steps = 360/FrameDegree;
  
  Xstep = (m)/Steps
  
  
  x = -0.5*m:Xstep:0.5*m;
  
  
  
  Va0PlotWave = plot(x,-1*m+real(Va0rotated*(rotationconstant).^(x/Xstep)));
  
  Va1PlotWave = plot(x-1.5*m,-1*m+real(Va1rotated*(rotationconstant).^(x/Xstep)));
  Vb1PlotWave = plot(x-1.5*m,-1*m+real(Va1rotated*a^2*(rotationconstant).^(x/Xstep)));
  Vc1PlotWave = plot(x-1.5*m,-1*m+real(Va1rotated*a*(rotationconstant).^(x/Xstep)));
  
  Va2PlotWave = plot(x+1.5*m,-1*m+real(Va2rotated*(rotationconstant).^(x/Xstep)));
  Vb2PlotWave = plot(x+1.5*m,-1*m+real(Va2rotated*a*(rotationconstant).^(x/Xstep)));
  Vc2PlotWave = plot(x+1.5*m,-1*m+real(Va2rotated*a^2*(rotationconstant).^(x/Xstep)));
  
VanWave =   plot(x+1.5*m,real(Vanrotated*(rotationconstant).^(x/Xstep)));
VbnWave =   plot(x+1.5*m,real(Vbnrotated*(rotationconstant).^(x/Xstep)));
VcnWave =   plot(x+1.5*m,real(Vcnrotated*(rotationconstant).^(x/Xstep)));  
  
  
set(VanPlot, 'Color', 'blue','linewidth',5);
set(Va1Plot, 'Color', 'blue','linewidth',5);
set(Va0Plot, 'Color', 'blue','linewidth',5);
set(Va2Plot, 'Color', 'blue','linewidth',5);


set(VbnPlot,'Color', 'green','linewidth',5);
set(Vb0Plot,'Color', 'green','linewidth',5);
set(Vb1Plot,'Color', 'green','linewidth',5);
set(Vb2Plot,'Color', 'green','linewidth',5);


set(VcnPlot, 'Color', 'red','linewidth',5);
set(Vc1Plot, 'Color', 'red','linewidth',5);
set(Vc0Plot, 'Color', 'red','linewidth',5);
set(Vc2Plot, 'Color', 'red','linewidth',5);

set(Va0PlotWave,'Color','red','linewidth',5);

set(Va1PlotWave, 'Color', 'blue','linewidth',5);
set(Vb1PlotWave, 'Color', 'green','linewidth',5);
set(Vc1PlotWave, 'Color', 'red','linewidth',5);

set(VanWave, 'Color', 'blue','linewidth',5);
set(VbnWave, 'Color', 'green','linewidth',5);
set(VcnWave, 'Color', 'red','linewidth',5);


set(Va2PlotWave, 'Color', 'blue','linewidth',5);
set(Vb2PlotWave, 'Color', 'green','linewidth',5);
set(Vc2PlotWave, 'Color', 'red','linewidth',5);

if GenerateGif == 1
filename=sprintf('output/%05d.gif',t);
print(filename);

endif

########UNCOMMENTT THE NEXT TWO LINES TO MAKE THE PROGRAM WRITE GIF FILES FOR EACH FRAME
##

pause(0.01);
  delete(VanPlot);
  delete(VbnPlot);
  delete(VcnPlot);
  delete(Va1Plot);
  delete(Va0Plot);
  delete(Va2Plot);
  delete(Vb1Plot);
  delete(Vb0Plot);
  delete(Vb2Plot);
  delete(Vc1Plot);
  delete(Vc0Plot);
  delete(Vc2Plot);
  delete(Va0PlotWave);
  delete(Va1PlotWave)
  delete(Vb1PlotWave)
  delete(Vc1PlotWave)
  delete(Va2PlotWave)
  delete(Vb2PlotWave)
  delete(Vc2PlotWave)
  delete(VanWave)
  delete(VbnWave)
  delete(VcnWave)
  

close figure





endfor