// Variables:
//      t: Time, passed from PSIM by value
//   delt: Time step, passed from PSIM by value
//     in: input array, passed from PSIM by reference
//    out: output array, sent back to PSIM (Note: the values of out[*] can
//         be modified in PSIM)

// The maximum length of the input and output array "in" and "out" is 30.

// Las entradas del sistema son
// in[0]= fr; frecuencia de la red
// in[1]= Vref; tension de referencia sin normalizar
// in[2]= fs; frecuencia de conmutación: La frecuencia de muestreo es 2xfs
// in[3]= Vcc; tension de la barra CC
// in[4]= vc; tension de salida

#include<math.h>

//#define	K11		  0.440787706646521		// Gan. RE iLs (para fm=20KHz).
//#define	K12		 25.848651939217628		// Gan. RE Vc (para fm=20KHz).
//#define	K13		  1.381029356307010		// Gan. RE ukm1 (para fm=20KHz).
//#define	K1			  0.798890354811086		// Gan. Servo (para fm=20KHz).

// Aqui se implementa un controlador proyectado en matlab, que consta de un compensador convencional
// proporcional derivativo y un modelo interno, con N =100 y filtro Q(z)=0.25 z^-1+0.5+0.25 z.
// La frecuencia de muestreo se proyecto para 5 kHz.
void __declspec(dllexport) simuser (double t, double delt, double *in, double *out)

{
   static double Vbase=311;
   static double Ibase=64.28;
   static double vc=0; // tension de fase en bornes del capacitor de filtro
   static double iL=0;   // Corriente en el inductor.
   static double uk=0, ukm1=0, upd=0;
   static double uc=0;
   static double e=0;
   static double vk=0, vkm1=0;
   static double ekm1=0, ekm2=0;
   static double r, rd;
   static double rampa;
   static double fr=in[0];
   static double Tr=1/fr;
   static double K11= 0.880504424685516;
   static double K12=-0.365608376983481;
   static double K13= 1.636148773744687;
   static double K1= 0.522979790723460;
   static double Ksv=1.816513786671996;
   static double k1=0.02, k2=-0.135;
   static double uklim=1.0;
   static double Tiempo=6;						// N° de periodos
   static int cont=0, Ncont, calculo_u, muestreo, direccion = 1;
   static int PWM1, PWM2;
   static int ComparReg1, ComparReg2;
   static double Vcc=in[3];
   static double fn=Vbase/Vcc;

   /* Periodo de muestreo */
	static double Ts = 1/in[2];

// Se calcula el numero de conteos en un periodo
   Ncont=Ts/delt;


   if (calculo_u==1)
      {
      calculo_u=0;

      r = in[1]/Vbase; // poner una senoidal con la amplitud deseada para la tension de referencia
      rampa = t;
      // Variación inicial de Vref (generación de rampa):
      if (t < 0.12)
           {rd = (1/(Tiempo*Tr))*rampa*r;} // Calcula incremento de referencia.
      else
      	  {rd = r;}
 	 /*********************************************************************************/
    /*                         Se calcula el error de tension                        */
    /*********************************************************************************/
    	//rd = 1;
      e = rd - vc;
    /*********************************************************************************/
    /*                   Calculo de la estrategia de control                         */
    /*********************************************************************************/
      upd = k1*ekm1 + k2*ekm2 + 1.0*rd;

	   vk = vkm1 + e;
      //vk = (1-Ksv*K1)*vkm1 + Ksv*ukm1 + Ksv*K11*iL + Ksv*K12*vc + Ksv*K13*ukm1 + e;

      uk = K1*vk - K11*iL - K12*vc - K13*ukm1;

      //if(uk < 0) {uk = 0;}					// SI uk es negativa, se limita a uk = 0.

      if(uk > uklim)                           // Si uk supera un valor dado,
      	{uk = uklim;                          // limita acción de control.
         vk = (uk + K11*iL + K12*vc + K13*ukm1)/K1;}  // Recalcula vkm1.
      if(uk < -uklim)                           // Si uk supera un valor dado,
      	{uk = -uklim;                          // limita acción de control.
         vk = (uk + K11*iL + K12*vc + K13*ukm1)/K1;}  // Recalcula vkm1.

      vkm1 = vk;				 		       // Obtiene servo anterior.
      ukm1 = uk;				 		       // Obtiene acción de control anterior.
      ekm2 = ekm1;
      ekm1 = e;

      //uk = upd;
      //uk = rd; //Para lazo abierto

    uc=uk*(Ncont/4)*fn;
    } //termina el calculo de leyes de control, muestreo, determinación del sector, y tiempo de aplicación de las llaves


   //**********************************************************//
	// Incrementamos o decrementamos el contador segun corresponda
   cont=cont+direccion;
   if (cont < 0)
      {
		direccion=1;
      muestreo=0;
		}
   if(cont == Ncont/2)
		{
		direccion=-1;
      muestreo=1;
		}

   if (cont==Ncont/2)  //((cont==0) || (cont==Ncont/2))
  		{
      ComparReg1= uc+Ncont/4;
      ComparReg2= -uc+Ncont/4;
      }

   // Maquina de estado que define el patron PWM
    if (ComparReg1 > cont)
        PWM1=1;
    else
        PWM1=0;

    if (ComparReg2 > cont)
        PWM2=1;
    else
        PWM2=0;

    if(muestreo==1)
    	{
       muestreo=0;
       // Tension de la salida del filtro
        vc = in[4]/Vbase;
        iL = in[5]/Ibase;
        calculo_u=1;
      }
// Se actualizan los valores de las salidas
	out[0] = PWM1; //brazo 1
   out[1] = PWM2; //brazo 2
	out[2] = uk;
   out[3] = rd;
   out[4] = vk;
   out[5] = e;
}