function [F_y,M_z]=extension(alpha,phi_t)
% phi_t=6;
% Magic Formula
F_z=3000;
F_z0=3000;
R_0=0.3;
R_e=0.3;
p_Cx1=1.65;
p_Kx1=12;
p_Cy1=1.3;
p_Ey1=-1;
p_Ky1=10;
p_Vy1=0;
p_Dx1=1;
p_Kx2=10;
p_Dy1=1;
p_Ey2=0;
p_Ky2=1.5;
p_Vy2=0;
p_Dx2=0;
p_Kx3=-0.6;
p_Dy2=0;
p_Ey3=0;
p_Ky3=0;
p_Vy3=0.15;
p_Ex1=-0.5;
p_Hx1=0;
p_Dy3=0;
p_Ey4=0;
p_Hy1=0;
p_Vy4=0;
p_Ex2=0;
p_Hx2=0;
p_Hy2=0;
p_Ex3=0;
p_Vx1=0;
p_Hy3=0.25;
p_Ex4=0;
p_Vx2=0;
q_Bz1=6;
q_Bz2=-4;
q_Bz3=0.6;
q_Bz4=0;
q_Bz5=0;
q_Bz9=0;
q_Bz10=0.7;
q_Cz1=1.05;
q_Dz1=0.12;
q_Dz2=-0.03;
q_Dz3=0;
q_Dz4=0;
q_Dz6=0;
q_Dz7=0;
q_Dz8=0.6;
q_Dz9=0.2;
q_Ez1=-10;
q_Ez2=0;
q_Ez3=0;
q_Ez4=0;
q_Ez5=0;
q_Hz1=0;
q_Hz2=0;
q_Hz3=0;
q_Hz4=0;
r_Bx1=5;
r_Bx2=8;
r_Cx1=1;
r_Hx1=0;
r_By1=7;
r_By2=2.5;
r_By3=0;
r_Cy1=1;
r_Hy1=0.02;
r_Vy1=0;
r_Vy2=0;
r_Vy3=-0.2;
r_Vy4=14;
r_Vy5=1.9;
r_Vy6=10;
s_sz1=0;
s_sz2=-0.1;
s_sz3=-1.0;
s_sz4=0;
p_Dxphi1=0.4;
p_Dxphi2=0;
p_Dxphi3=0;
p_Kyphi1=1;
p_Dyphi1=0.4;
p_Dyphi2=0;
p_Dyphi3=0;
p_Dyphi4=0;
p_Hyphi1=1;
p_Hyphi2=0.15;
p_Hyphi3=0;
p_Hyphi4=-4;
p_epsilongammaphi1=0;
p_epsilongammaphi2=0;
q_Dtphi1=10;
q_Crphi1=0.2;
q_Crphi2=0.1;
q_Brphi1=0.1;
q_Drphi1=1;
q_Drphi2=-1.5;
A_mu=10;
%---------------------
lambda_Fz0=1;
lambda_mux=1;
lambda_muy=1;
lambda_muV=0;
lambda_Kxkappa=1;
lambda_Kya=1;
lambda_Cx=1;
lambda_Cy=1;
lambda_Ex=1;
lambda_Ey=1;
lambda_Hx=1;
lambda_Hy=1;
lambda_Vx=1;
lambda_Vy=1;
lambda_Kygamma=1;
lambda_Kzgamma=1;
lambda_t=1;
lambda_Mr=1;
lambda_xa=1;
lambda_yk=1;
lambda_Vyk=1;
lambda_s=1;
lambda_Cz=1;
lambda_Mx=1;
lambda_My=1;
r_Ex1=0.3;
r_Ex2=0.3;
% alpha=(-20:0.1:20)*pi/180;
V_cx=30;
epsilon_x=0.001;
V_0=sqrt(9.8.*0.3);
kappa=0;
% zeta_i=1;
F_z0_apostrophe=3000;
q_Bz6=1;
cos_alpha_apostrophe=cos(alpha);
lambda_Ky=1;
gamma=8.*pi./180;
epsilon_gamma=0.001;
epsilon_K=0.001;
lambda_Mphi=1;
epsilon_r=0.001;
%---------------------
%longgitudinal force
F_z0_star=lambda_Fz0.*F_z0;%1
df_z=(F_z-F_z0_star)./F_z0_star;%2
alpha_star=tan(alpha);%3 unknown
gamma_star=sin(gamma);%4 unknown
% kappa=-V_sx./abs(V_cx);%5 unknown
% cos(alpha_star)=V_cx./(V_c+epsilon_V);%6,6a unknown
zeta_i=1;
lambda_mux_star=lambda_mux;%7 unknown
lambda_muy_star=lambda_muy;% unknown
% lambda_mux_star=lambda_mux./(1+lambda_muV.*V_s./V_0);%7 unknown
% lambda_muy_star=lambda_muy./(1+lambda_muV.*V_s./V_0);% unknown
lambda_mux_apostrophe=A_mu.*lambda_mux_star./(1+(A_mu-1).*lambda_mux_star);%8
lambda_muy_apostrophe=A_mu.*lambda_muy_star./(1+(A_mu-1).*lambda_muy_star);
C_x=p_Cx1.*lambda_Cx;%11 >0
mu_x=(p_Dx1+p_Dx2.*df_z).*lambda_mux_star;%13 >0
D_x=mu_x.*F_z.*zeta_i;%12 >0
S_Hx=(p_Hx1+p_Hx2.*df_z).*lambda_Hx;%17
K_xk=F_z.*(p_Kx1+p_Kx2.*df_z).*exp(p_Kx3.*df_z).*lambda_Kxkappa;%15 =C_Fk
kappa_x=kappa+S_Hx;%10
E_x=(p_Ex1+p_Ex2.*df_z+p_Ex3.*df_z.^2).*(1-p_Ex4.*sign(kappa_x)).*lambda_Ex;%14 <=1
B_x=K_xk./(C_x.*D_x+epsilon_x);%16
S_Vx=F_z.*(p_Vx1+p_Vx2.*df_z).*lambda_Vx.*lambda_mux_apostrophe.*zeta_i;%18
F_x0=D_x.*sin(C_x.*atan(B_x.*kappa_x-E_x.*(B_x.*kappa_x-atan(B_x.*kappa_x))))+S_Vx;%9
% S_Vx=F_z.*(p_Vx1+p_Vx2.*df_z).*(abs(V_cx)./(epsilon_Vx+abs(V_cx))).*lambda_Vx.*lambda_mux_apostrophe.*zeta_i;%18
%laterral force
C_y=p_Cy1.*lambda_Cy;%21 >0
mu_y=(p_Dy1+p_Dy2.*df_z).*(1-p_Dy3.*gamma_star.^2).*lambda_muy_star;%23 >0
D_y=mu_y.*F_z.*zeta_i;%22
K_ya0=p_Ky1.*F_z0_apostrophe.*sin(2.*atan(F_z./(p_Ky2.*F_z0_apostrophe))).*lambda_Kya;%25
zeta_3=cos(atan(p_Kyphi1.*R_0.^2.*phi_t.^2));%79
K_ya=K_ya0.*(1-p_Ky3.*gamma_star.^2).*zeta_3;%26
B_y=K_ya./(C_y.*D_y);%27
zeta_0=0;%83
phi=phi_t-sin(gamma)./0.3;
B_yphi=p_Dyphi1.*(1+p_Dyphi2.*df_z).*cos(atan(p_Dyphi3.*tan(alpha)));%78
zeta_2=cos(atan(B_yphi.*(R_0.*abs(phi_t)+p_Dyphi4.*sqrt(R_0.*abs(phi_t)))));%77
D_Hyphi=(p_Hyphi2+p_Hyphi3.*df_z).*sign(V_cx);%86
C_Hyphi=p_Hyphi1;%85
E_Hyphi=p_Hyphi4;%87
K_ygamma0=(p_Hy3.*K_ya0+F_z.*(p_Vy3+p_Vy4.*df_z)).*lambda_Kygamma;%30
K_yRphi0=K_ygamma0./(1-epsilon_gamma);%89
B_Hyphi=K_yRphi0./(C_Hyphi.*D_Hyphi.*(K_ya0+epsilon_K));%88
S_Hyphi=D_Hyphi.*sin(C_Hyphi.*atan(B_Hyphi.*R_0.*phi-E_Hyphi.*(B_Hyphi.*R_0.*phi-atan(B_Hyphi.*R_0.*phi))));%80
S_Vygamma=F_z.*(p_Vy3+p_Vy4.*df_z).*gamma_star.*zeta_2.*lambda_Kygamma.*lambda_muy_apostrophe;%82
K_ya_apostrophe=K_ya+0.001;%39
S_Hy=(p_Hy1+p_Hy2.*df_z).*lambda_Hy+S_Hyphi-S_Vygamma./K_ya_apostrophe;%81
zeta_4=1+S_Hyphi-S_Vygamma./K_ya_apostrophe;%84
S_Hy=(p_Hy1+p_Hy2.*df_z).*lambda_Hy+p_Hy3.*gamma_star.*lambda_Kygamma.*zeta_0+zeta_4-1;%28
alpha_y=alpha_star+S_Hy;%20
E_y=(p_Ey1+p_Ey2.*df_z).*(1-(p_Ey3+p_Ey4.*gamma_star.^2).*sign(alpha_y)).*lambda_Ey;%24 <=1
S_Vy=F_z.*((p_Vy1+p_Vy2.*df_z).*lambda_Vy+(p_Vy3+p_Vy4.*df_z).*gamma_star.*lambda_Kygamma).*lambda_muy_apostrophe.*zeta_2;%29
F_y0=D_y.*sin(C_y.*atan(B_y.*alpha_y-E_y.*(B_y.*alpha_y-atan(B_y.*alpha_y))))+S_Vy;%19
% M_z0=M_z0_apostrophe+M_zr0;%31
% M_z0_apostrophe=-t_0.*F_y0;%32
% t_0=D_t.*cos(C_t.*atan(B_t.*alpha_t-E_t.*(B_t.*alpha_t-atan(B_t.*alpha_t)))).*cos_alpha_apostrophe;%33
S_Ht=q_Hz1+q_Hz2.*df_z+(q_Hz3+q_Hz4.*df_z).*gamma_star;%35
alpha_t=alpha_star+S_Ht;%34
% M_zr0=D_r.*cos(C_r.*atan(B_r.*alpha_r));%36
S_Hf=S_Hy+S_Vy./K_ya_apostrophe;%38
alpha_r=alpha_star+S_Hf;%37
B_t=(q_Bz1+q_Bz2.*df_z+q_Bz3.*df_z.^2).*(1+q_Bz5.*abs(gamma_star)+q_Bz6.*gamma_star.^2).*lambda_Kya./lambda_muy_star;%40 (>0)
C_t=q_Cz1;%41 (>0)
D_t0=F_z.*(R_0./F_z0_apostrophe).*(q_Dz1+q_Dz2.*df_z).*lambda_t.*sign(V_cx);%42
zeta_5=cos(atan(q_Dtphi1.*R_0.*phi));%91
D_t=D_t0.*(1+q_Dz3.*abs(gamma_star)+q_Dz4.*gamma_star.^2).*zeta_5;%43
E_t=(q_Ez1+q_Ez2.*df_z+q_Ez3.*df_z.^2).*(1+(q_Ez4+q_Ez5.*gamma_star).*(2./pi).*atan(B_t.*C_t.*alpha_t));%44
zeta_6=cos(atan(q_Brphi1.*R_0.*phi));% 102
epsilon_gamma=p_epsilongammaphi1.*(1+p_epsilongammaphi2.*df_z);%90
M_zphiinfinity=q_Crphi1.*mu_y.*R_0.*F_z.*sqrt(F_z./F_z0_apostrophe).*lambda_Mphi;%95
C_Drphi=q_Drphi1;%96
E_Drphi=q_Drphi2;%97
K_zgammar0=F_z.*R_0.*(q_Dz8+q_Dz9.*df_z).*lambda_Kzgamma;%99
D_Drphi=M_zphiinfinity./sin(0.5.*pi.*C_Drphi);%94
B_Drphi=K_zgammar0./((C_Drphi.*D_Drphi+epsilon_r).*(1-epsilon_gamma));%98
K_zgamma0=K_zgammar0-D_t0.*K_ygamma0;% (=C_Mgamma) 100
K_zRphi0=K_zgamma0./(1-epsilon_gamma);% (C_Mphi./R_0) 101
C_yk=r_Cy1;%63
E_yk=0.3+0.3.*df_z;%64
S_Hyk=r_Hy1;%65
B_yk=r_By1.*cos(atan(r_By2.*(alpha_star-r_By3))).*lambda_yk;%62
D_Vyk=mu_y.*F_z.*(r_Vy1+r_Vy2.*df_z+r_Vy3.*gamma_star).*cos(atan(r_Vy4.*alpha_star));%67
S_Vyk=D_Vyk.*sin(r_Vy5.*atan(r_Vy6.*kappa));%66
kappa_s=kappa+S_Hyk;%61
G_yk0=cos(C_yk.*atan(B_yk.*S_Hyk-E_yk.*(B_yk.*S_Hyk-atan(B_yk.*S_Hyk))));%60
G_yk=cos(C_yk.*atan(B_yk.*kappa_s-E_yk.*(B_yk.*kappa_s-atan(B_yk.*kappa_s))))./G_yk0;%59
F_y=G_yk.*F_y0+S_Vyk;%58
M_zphi90=M_zphiinfinity.*(2./pi).*atan(q_Crphi2.*R_0.*abs(phi)).*G_yk;%103
D_rphi=D_Drphi.*sin(C_Drphi.*atan(B_Drphi.*R_0.*phi-E_Drphi.*(B_Drphi.*R_0.*phi-atan(B_Drphi.*R_0.*phi))));%93
zeta_8=1+D_rphi;%92
B_xphi=p_Dxphi1.*(1+p_Dxphi2.*df_z).*cos(atan(p_Dxphi3.*kappa));%106
zeta_7=(2./pi).*acos(M_zphi90./(abs(D_rphi)+epsilon_r));%104
zeta_1=cos(atan(B_xphi.*R_0.*phi));%105
B_r=(q_Bz9.*lambda_Ky./lambda_muy_star+q_Bz10.*B_y.*C_y).*zeta_6;%45
C_r=zeta_7;%46
D_r=F_z.*R_0.*((q_Dz6+q_Dz7.*df_z).*lambda_Mr.*zeta_2+(q_Dz8+q_Dz9.*df_z).*gamma_star.*lambda_Kzgamma.*zeta_0).*cos_alpha_apostrophe.*lambda_muy_star.*sign(V_cx)+zeta_8-1;%47
% K_za0=D_t0.*K_ya0;%48
% K_zgamma0=F_z.*R_0.*(q_Dz8+q_Dz9.*df_z).*lambda_Kzgamma-D_t0.*K_ygamma0;%49
C_xa=r_Cx1;%55
S_Hxa=r_Hx1;%57
E_xa=r_Ex1+r_Ex2.*df_z;%56
B_xa=r_Bx1.*cos(atan(r_Bx2.*kappa)).*lambda_xa;%54
alpha_s=alpha_star+S_Hxa;%53
G_xa0=cos(C_xa.*atan(B_xa.*S_Hxa-E_xa.*(B_xa.*S_Hxa-atan(B_xa.*S_Hxa))));%52
G_xa=cos(C_xa.*atan(B_xa.*alpha_s-E_xa.*(B_xa.*alpha_s-atan(B_xa.*
