pkg load signal;

%Transistor Parameters
Ron = 100E-3;
tdon = 18E-9;
tr = 13E-9;
tdoff = 59E-9;
tf = 11E-9;
Crss = 2.5E-12; %Gate-Drain Capacity
Vth = 5.5;  %Conduction Voltage Threshold
Vgg = 8;  %Trigger Voltage
Gm = 25;   %Transconductance
Rg = 10;   %Gate Resistor

%Grid and Inverter Parameters
Vg = 230;           %Grid voltage
f_net = 50;         %Grid Frecuency
T_net = 1/f_net;
w_net = 2*pi*f_net;
Vdc = 400;        %Tensión bus de continua
P = 2300;         %Inverter Power
I = P/Vg;         %RMS current
Ip = I*sqrt(2);   %Peak current

%Commutation Parameters
f_sw = 20000; %Switching Frecuency
T_sw = 1/f_sw;
w_sw = 2*pi*f_sw;
modulation_index = 0.6;

%Current, carrier, modulator waves generation
delta = T_sw/20000;      %Time interval between test points

t = 0:delta:T_net;    %Discrete time domain generation
It = Ip*sin(w_net*t);
carrier = sawtooth(w_sw*t,1/2);
modula = modulation_index*sin(w_net*t);
trigger = modula>carrier;

%Deltas quantity transistor parameters
deltas_tdon = uint32(tdon/delta);
deltas_tr = uint32(tr/delta);
deltas_tdoff = uint32(tdoff/delta);
deltas_tf = uint32(tf/delta);

sw_on = false;
Eon = 0;
ERon = 0;
Eoff = 0;

tic
i = 1;
do
  if((trigger(i) == true) && (sw_on == false))
  i = (i + deltas_tdon + deltas_tr);
  %Current Rise
  i_t = abs(It(i));
  Eon = Eon + i_t*Vdc*tr/2;
  %Voltage Fall
  deltaVds = i_t*Ron - Vdc;
  deltaVdsdt = (Vgg - ((i_t/Gm) + Vth))/(-Crss*Rg);
  tfv = (deltaVds/deltaVdsdt);
  Eon = Eon + i_t*Vdc*tfv/2;
  deltas_tfv = uint32(tfv/delta);
  i = (i + deltas_tfv);
  sw_on = true;

  elseif((trigger(i) == true) && (sw_on == true))
  %Conduction Losses
  ERon = ERon + (It(i)^2)*Ron*delta;
  i = (i + 1);

  elseif((trigger(i) == false) && (sw_on == true))
  i = (i + deltas_tdoff);
  i_t = abs(It(i));
  Eoff = Eoff + (i_t^2)*Ron*tdoff;
  %Voltage Rise
  deltaVds = i_t*Ron - Vdc;
  deltaVdsdt = (-Rg*i_t)/(Rg*Crss);
  trv = deltaVds/deltaVdsdt;
  Eoff = Eoff + i_t*Vdc*trv/2;
  deltas_trv = uint32(trv/delta);
  i = (i + deltas_trv);
  %Current Fall
  Eoff = Eoff + abs(It(i)*Vdc*tf/2);
  i = (i + deltas_tf);
  sw_on = false;

  else
  i = (i + 1); %No case happened
  endif

until(i>=length(t))
toc

Won = 4*Eon*f_net
WRon = 4*ERon*f_net
Woff = 4*Eoff*f_net
Wtotal = Won+WRon+Woff