clear all
clc


% INDUCTION MACHINE
Pn = 3600;                      % nominal stator power [W]
fn = 50;                        % nominal frequency [Hz]
Temn = 37;                      % nominal torque [Nm] (calculated, formula in Notes)
Un = 380;                       % nominal stator voltage [V];
pp = 3;                         % pole pairs [-]
ns = 60*fn/pp;                  % synchronous rotational speed [rpm]
J = 0.2286;                     % estimate of whole system inertia [kg*m^2]
smax = 1/2;                     % maximum slip [-]
D = 1.8e-6;                     % friction factor [-]
%nn = 935;                      % nominal rotational speed from machine label [rpm], actual value from simulation is 860

Rs = 1.688;                     % stator resistance [ohm]
Rr = 3.5;                       % rotor resistance [ohm]
Ls_sigma = 0.01196;             % leakage stator inductance [H]
Lr_sigma = 0.01298;             % leakage rotor inductance [H]
Lm = 0.1735;                    % magnetizing inductance [H]
Ls = Lm + Ls_sigma;             % stator inductance [H]
Lr = Lm + Lr_sigma;             % rotor inductance [H]
sigma = 1 - Lm^2/(Ls*Lr);       % leakage coefficient[-]
Psi_stator = sqrt(2/3)*Un/(2*pi*fn);    % cca stator flux [Wb]


% SIMULATION PARAMETERS
fsw = 5e3;                      % frequency of the triangular carrier signal [Hz]
Ts = 1e-6;                      % simulation step [sec]


% WIND TURBINE MODEL
radius = 1.70;                  % radius of the wind turbine [m]
N = ns/325;                     % gearbox ratio [-]
ro = 1.229;                     % avg air density, sea level at 15°C [kg/m3]


% MPPT
Cp_max = 0.44;
lambda_opt = 7.2;
k_opt = (0.5*ro*pi*(radius^5)*Cp_max)/((N^3)*(lambda_opt^3));


% ROTOR AND GRID SIDE CONVERTERS PARAMETERS (RSC AND GSC)
U_DC = 2*sqrt(2/3)*Un;          % DC bus voltage [V]
C_DC = 4.7e-3;                  % DC bus capacitance
L_grid = 62.5e-3;               % grid side L filter's inductance


% RSC CONTROL
kp_id = 150;
ki_id = 10000;

kp_iq = kp_id;
ki_iq = ki_id;

kp_n = 50;
ki_n = 0;


% GSC CONTROL
K_Pg = 1/(1.5*Un*sqrt(2/3));
K_Qg = -K_Pg;

kp_idg = 50;
ki_idg = 9000;

kp_iqg = kp_idg;
ki_iqg = ki_idg;

kp_v = -950;
ki_v = -285000;


% POWER AND TORQUE COEFFICIENTS
beta = 0;                       % blade pitch angle [?]
index = 1;                      % [-]

    for lambda = 0.1:0.01:11.8
        lambdai(index) = (1./((1./(lambda-0.02.*beta)+(0.003./(beta^3+1)))));
        Cp(index) = 0.73.*(151./lambdai(index)-0.58.*beta-0.002.*beta^2.14-13.2).*(exp(-18.4./lambdai(index)));
        Ct(index) = Cp(index)/lambda;
        index = index+1;
    end
    tab_lambda = [0.1:0.01:11.8];