function x_trim = trimConfig(plane)
v0 = plane.v0;
m = plane.m;
M = plane.M;
rho = plane.rho;
g = 9.81;

Ixx = plane.Jxx;    % moment setrvacnosti pro osu x
Izz = plane.Jyy;    % moment setrvacnosti pro osu y
Iyy = plane.Jzz;    % moment setrvacnosti pro osu z
Ixz = plane.Jxy;    % moment setrvacnosti pro osu z
cfgmatfile = 'MyAircraft';


% stredni aerodynicka tetiva
c = 1.91827; % m
% rozpeti
b = 19.0789; % m
% plocha kridla
S = 34.88; % m^2


%% Aerodynamicke charakteristiky
% vztlakovy koeficient
    % pocatecni hodnota vztlakoveho soucinitele
    x_cy0       = [0,0.4,0.6,0.7,0.75,0.8,0.85,0.9]; 
    y_cy0       = [0.09,0.09,0.09,0.105,0.13,0.12,0.07,0.02];
    CL0 = interp1(x_cy0,y_cy0,M);
    
    % derivace soucinitele vztlaku podle uhlu nabehu
    if M <= 0.3
        CLa = 4.4418;
    else
        CLa = - 62.716*M^5 + 122.09*M^4 - 85.797*M^3 + 28.153*M^2 ...
                   - 3.8433*M + 4.541;
    end;

    % derivace soucinitele vztlaku podle vychylky vyskovky
    CLde = 0.305;

    % derivace soucinitele vztlaku podle bezrozmerne zmeny uhlu nabehu
    if M<0.3
        Mp = 0.3;
        CLalphadot = -49.167 + 558.542*Mp - 2284.051*Mp^2 ... 
        + 4412.037*Mp^3 - 4067.130*Mp^4 + 1435.185*Mp^5;
    elseif M>0.75
        Mp = 0.75;
        CLalphadot = -49.167 + 558.542*Mp - 2284.051*Mp^2 ... 
        + 4412.037*Mp^3 - 4067.130*Mp^4 + 1435.185*Mp^5;
    else
        CLalphadot = -49.167 + 558.542*M - 2284.051*M^2 ... 
        + 4412.037*M^3 - 4067.130*M^4 + 1435.185*M^5;
    end

    % derivace soucinitele vztlaku podle bezrozmerne rychlosti kloneni
    CLq = 5.7414 - 9.3115*M + 20.8*M^2 - 13.465*M^3;


% odporovy soucinitel
    % pocatecni hodnota odporoveho soucinitele
    CD0 = 0.08;
    
    % derivace soucinitele odporu podle uhlu nabehu
    if M<0.3
        Mp = 0.3;
        CDa = -1*(10.081 - 23.378*Mp + 13.048*Mp^2);
    else
        CDa = -1*(10.081 - 23.378*M + 13.048*M^2);
    end

    % derivace soucinitele odporu podle vychylky vyskovky
    if M<0.3
        Mp = 0.3;
        CDde =( 0.1857*Mp^2 - 0.1897*Mp + 0.0649);
    elseif M>0.7
        Mp = 0.7;
        CDde =( 0.1857*Mp^2 - 0.1897*Mp + 0.0649);
    else
        CDde =( 0.1857*M^2 - 0.1897*M + 0.0649);
    end


% bocny soucinitel
    cyPom = m*9.81/(0.5*rho*v0^2*S);
    
    % derivace soucinitele bocne sily podle uhlu vyboceni
    CYbeta = (-0.946 - 0.0333*cyPom + 0.0032*M - 0.16*M^2);
    
    % derivace soucinitele bocne sily podle vychylky kridelek
    CYda = -0.02;

    % derivace soucinitele bocne sily podle vychylky smerovky
    CYdr = (-0.1845 + 0.0325*cyPom -0.0093*M -0.0228*M^2 ...
        -0.012*M^3);

    %  derivace soucinitele bocne sily podle kloneni
    CYp = 0;
    
    %  derivace soucinitele bocne sily podle zatacive rychlosti
    CYr = 0;

    
% soucinitel klopiveho momentu
    y_mz0 = [0.025,0.025,0.03,0.025,0.02,0.05,-0.02,-0.03];
    % pocatecni hodnota klopiveho soucinitele
    Cm0 = interp1(x_cy0,y_mz0,M);

    % derivace soucinitele klopiveho momentu podle uhlu nabehu
    if M <= 0.3
        Cma = -0.4899;
    else
        Cma = 5.9964*M^4 - 10.271*M^3 + 6.7921*M^2 - 2.0056*M -0.2708;
    end    

    % derivace soucinitele klopiveho momentu podle vychylky vyskovky
    Cmde = -0.88;

    % derivace soucinitele klopiveho momentu podle zmeny uhlu nabehu
    if M <= 0.3
        Cmalphadot = -1.38271969;
    else
        Cmalphadot = 502.427*M^5 - 1303.11*M^4 + 1330.2481*M^3 ...
            - 668.0225*M^2 + 164.08*M - 17.0671;
    end    
    
    % derivace soucinitele klopiveho momentu podle klopeni
    if M <= 0.3
        Cmq = -4.9934668;
    else
        Cmq = 2.5926*M^3 - 8.7063*M^2 + 4.8590*M - 5.7376;
    end    

    
% soucinitel kloniveho momentu
    % derivace soucinitele klopiveho momentu podle uhlu vyboceni
    if M<0.45
        Clbeta = -0.056 - 0.0576*cyPom - 0.0088*M;
    else
        Clbeta = -0.06 - 0.0576*cyPom - 0.0285*(M - 0.45);
    end
    
    % derivace soucinitele klopiveho momentu podle vychylky kridelek
    if M<=0.6
        Clda = (-0.156 + 0.04*(M -0.2));
    else
        Clda = -1.1039 + 4.6696*M - 7.5941*M^2 ...
            + 4.1483*M^3;
    end
    
    % derivace soucinitele klopiveho momentu podle vychylky smerovky
    Cldr = (-0.03 + 0.008*cyPom)*1;
    
    % derivace soucinitele klopiveho momentu podle kloneni
    Clp = -0.5365 - 1.5267*M + 10.781*M^2 ...
        -20.666*M^3 + 12.4242*M^4;
    
    % derivace soucinitele klopiveho momentu podle zatacive rychlosti
    Clr = (-0.075 - 0.215*cyPom -0.016*M);

% soucinitel zataciveho momentu
    % derivace soucinitele zataciveho momentu podle uhlu vyboceni
    Cnbeta = (-0.173 + 0.01667*cyPom + 0.0092*M -0.1*M^2 +...
        0.0208*M^3);
    
    % derivace soucinitele klopiveho momentu podle vychylky kridelek
    Cnda = -0.0015;

    % derivace soucinitele klopiveho momentu podle vychylky smerovky
    Cndr = (0.015*cyPom - 0.095 -0.0056*M -0.0089*M^2 ...
        - 0.0068*M^3);
    
    % derivace soucinitele klopiveho momentu podle kloneni
    if M<0.7
        Cnp = -0.07435 + 0.045*cyPom - 0.0092*M;
    else
        Cnp = -0.08079 + 0.044*cyPom - 0.0333*(M-0.7);
    end

    % derivace soucinitele klopiveho momentu podle zatacive rychlosti
    Cnr = -0.2291 - 0.0065*cyPom -0.0236*cyPom^2 + 0.0069*cyPom^3 ...
        - 0.0219*M + 0.0097*M^2 - 0.0447*M^3;


%% trimovani letounu

CG = 0.245 ; 
rAC = plane.x_T;
Jgross = [Ixx Iyy Izz Ixz];
qb = 1/2*rho*v0^2;

vyber = 0;
x_trim = [ 0 0 0 0];
x_T0 = 0.2450;
x_T = plane.x_T;
if vyber == 1
    if M<0.4
        CD0p = 0.0223;
    elseif M>0.92 
        CD0p = 0.1032;
    else
        CD0p = 7.861*M^4 - 18.378*M^3 + 15.832*M^2 - 5.9269*M + 0.8349;
    end
    qb = 1/2*rho*v0^2;
    fnc = @(x)((Cma*x(1) + Cmde*x(2) + Cm0 + (x_T - x_T0)*(CL0 + CLa*x(1) + CLde*x(2)))^2 +...
        (CLa*x(1) + CLde*x(2) - m*g/(qb*S)*cos(x(3)) + CL0)^2 +...
        (CD0p + 0.00103*(x(1)*180/pi)^2 - 0.00308*x(1)*180/pi +... 
        6.2863E-07*(x(2)*180/pi)^3 + 4.1831E-05*(x(2)*180/pi)^2 -...
        3.1917E-04*(x(2)*180/pi) + 1.6834E-04 +...
        m*g/(qb*S)*sin(x(3)) - (x(4))/(qb*S))^2 + ...
        (x(1) - x(3))^2);
    [x,fval,exitflag,output] = fminsearch(fnc, ...
                [0, 0, 0, 8000], optimset('TolX',1e-15));
    x_trim = [x(1) x(2) x(3) x(4)];
    disp('-------------------------------------')
    disp('----- ---- Trim function ----- -----')
    disp(sprintf('AoA           [°]: %5.3d',x_trim(1)*180/pi))
    disp(sprintf('Elevator      [°]: %5.3d',x_trim(2)*180/pi))
    disp(sprintf('Pitch angle   [°]: %5.3d',x_trim(3)*180/pi))
    disp(sprintf('Engine thrust [N]: %6.1f',x_trim(4)))
    disp('-------------------------------------')
end

% ulozeni velicin do mat souboru
save(cfgmatfile);
