pkg load control 
pkg load io  
pkg load signal 

%eig_A=eig(sys.A);
%plot(real(eig_A),imag(eig_A),"rx")
%grid on



%% se define la variable para hacer el barrido
%cycle_vector = logspace(log10(0.5),log10(2),100);
cycle_vector = linspace(0,15,16);
color_plot = colormap(jet(length(cycle_vector))); %los colores van variando del azul al rojo


%all = axes('position',[.1,.1,.8,.8])
length(cycle_vector);
%% cycle (secalcula las nuevas matrices para cada valor del barrido)
for n_cycle = 1:length(cycle_vector)
  %hay que recalcular la matrz A para el barridod e la variable puede que haya que Volver a linealizar
  kv=cycle_vector(n_cycle)
  %kv=15;
  Tv=0.05;%1;%Constante de tiempo del PI de la componente d 
  kq=0.4;%Constante proporcional del PI de la componente q 
  Tq=0.02;%Constante de tiempo del PI de la componente q 
  kdroop_AC=0;%-3;%Pendiente droop Q/U_AC 
  droop_FSM=0.05;%Pendiente droop P/f 
  FSM_slope=droop_FSM/0.002;%Pendiente droop P/f en p.u.

  %Inner Loop
  ki=5000;%1% Constante integral de los controladores proporcional resonantes 
  kp=1.5;% Constante proporcional de los controladores proporcional resonantes 
  w=2*pi*50;
  wg0=2*pi*50.01;
  w_r=2*pi*50;
  wr=2*pi*50;
  a23 = 2*(wg0^2+wr^2);
  da23 = 4*wg0;
  a03 = (wr^2-wg0^2)^2;
  da03 = -4*wg0*(wr^2-wg0^2);
  b33 = ki;
  db33 = 0;
  b13 = ki*(wg0^2+wr^2);
  db13 = ki*2*wg0;
  c23 = -ki*wg0;
  dc23 = -ki;
  c03 = ki*wg0*(wr^2-wg0^2);
  dc03 = ki*(wr^2-3*wg0^2);

  %2-level VSC&Normalization block
  Vdref_perm=1.002;%1.003;
  Vqref_perm=0.085;%0.077;
  Vdc_diff_perm=0.5;%0.492;
  kdiff=(0.7)*Vdc_diff_perm;%Según DIgilent debería ser multiplicada por 1.4 en vez de 0.7 
  K_0=sqrt(3)/(2*sqrt(2));%Factor de modulación sinusoidal del VSC
  Ub_DC=640/2;%Tensión base en DC del VSC-Según DIgilent debería ser 640 en vez de 640/2 
  Ub_AC=275;%Tensión base en AC del VSC
  k0=(K_0)*Ub_DC/Ub_AC;
  Vdc_perm=1.001;%1.017;
  Pmd_perm=0.702;%0.702;%0.69;
  Pmq_perm=0;%0.053;
  Kvsc=k0*Vdc_perm;

  x1d_0= 0;
  x1q_0= 0;
  x2d_0=0;
  x2q_0=0;
  x3d_0=0;
  x3q_0=0;
  x4d_0=0;
  x4q_0=0;

%Función enable: para marcar la activación de la rama FSM
  enable=0;

%Definición de matrices A,B,C y D

  %Outer loop:Vdc-Q

  A_ol=[0 0;0 0];
  B_ol=[-kv/Tv kv/Tv 0 0 0 0 -kv*FSM_slope*enable/Tv kv*enable/Tv ...
  -kv*enable/Tv kv*enable/Tv;0 0 kq/Tq -kq/Tq -(kq/Tq)*kdroop_AC ...
  (kq/Tq)*kdroop_AC 0 0 0 0];
  C_ol=[1 0;0 1];
  D_ol=[-kv kv 0 0 0 0 -kv*FSM_slope*enable kv*enable -kv*enable kv*enable;0 ...
  0 kq -kq -kq*kdroop_AC kq*kdroop_AC 0 0 0 0];


  %Inner loop:PR current controller

  A_il=[0 1 0 0 0 0 0 0; 0 0 1 0 0 0 0 0;0 0 0 1 0 0 0 0; -(w^2-wg0^2) 0 ...
  -4*w^2 0 0 0 0 0;0 0 0 0 0 1 0 0;0 0 0 0 0 0 1 0;0 0 0 0 0 0 0 1;0 0 0 0 ...
  -(w^2-wg0^2) 0 -4*w^2 0];
  B_il=[0 0 0 0 0 0 0;0 0 0 0 0 0 0;0 0 0 0 0 0 0;0 0 -1 0 1 0 -4*w*x3d_0;0 ...
  0 0 0 0 0 0;0 0 0 0 0 0 0;0 0 0 0 0 0 0;0 0 0 -1 0 1 -4*w*x3q_0];
  C_il=[0 2*ki*w^2 0 ki 0 0 ki*w 0; 0 0 -ki*w 0 0 2*ki*w^2 0 ki];
  D_il=[1 0 -kp 0 kp 0 ki*(2*w*x2d_0+x3q_0+2*w^2*x1q_0); 0 1 0 -kp 0 kp ...
  ki*(2*w*x2q_0-2*w^2*x1d_0-x3d_0)];

  %Normalization block 
  D_nb=[kdiff 0 (0.7)*Vdref_perm;0 kdiff (0.7)*Vqref_perm];%Según DIgilent
#debería ser multiplicada por 1.4 en vez de 0.7 

  %2-level VSC:controlled voltage source-average model
  D_vsc=[k0*Vdc_perm 0 k0*Pmd_perm;0 k0*Vdc_perm k0*Pmq_perm];

  %Outer loop
  stname_ol = {'xv','xq'};

  inname_ol={'vdcref','vdc','qacref','qac','vacref','vac','fmeas','pdcref','paux','FSMcoordpos'};
  outname_ol={'idref','iqref'};
  sys_ol = ss(A_ol,B_ol,C_ol,D_ol,'stname',stname_ol,'inname',inname_ol,'outname',outname_ol);

  %Inner loop
  stname_il = {'xd1','xd2','xd3','xd4','xq1','xq2','xq3','xq4'};
  inname_il={'Vd','Vq','Id','Iq','idref','iqref','omegagiro'};
  outname_il={'Vdref','Vqref'};
  sys_il = ss(A_il,B_il,C_il,D_il,'stname',stname_il,'inname',inname_il,'outname',outname_il);

  %Normalization block 
  inname_nb={'Vdref','Vqref','Vdcdiff'};
  %inname_nb={'Vdref','Vqref'};
  outname_nb={'Pmd','Pmq'};
  sys_nb = ss (D_nb,'inname',inname_nb,'outname',outname_nb);

  %2-ln_cycleevel VSC
  inname_vsc={'Pmd','Pmq','vdc'};
  outname_vsc={'Vd_AC','Vq_AC'};
  sys_vsc = ss (D_vsc,'inname',inname_vsc,'outname',outname_vsc);

  %Assembled system1

  inname_group={'vdcref','vdc','qacref','qac','vacref','vac','fmeas','pdcref','paux','FSMcoordpos','Vd','Vq','Id','Iq','omegagiro','Vdcdiff'};

%inname_group={'vdcref','vdc','qacref','qac','vacref','vac','fmeas','pdcref','paux','FSMcoordpos','Vd','Vq','Id','Iq'};

  outname_group={'idref','iqref','Vdref','Vqref','Pmd','Pmq','Vd_AC','Vq_AC'};
  sys=connect(sys_ol,sys_il,sys_nb,sys_vsc,inname_group,outname_group);
  Acl=sys.A + sys.B*ones(16,10);

 endfor