function alpha=alphaPacejka(Fy,Fzi_N,Fzo_N,parvec)

%alphaPacejka
%
%   X = alphaPacejka(Fy,Fzi,Fzo,a) gives the slip angle (in degrees)
%   necessary to produce the side force Fy when the normal force on 
%   the inside tire is Fzi and the normal force on the outside tire
%   is Fzo for the parameter vector a.  All forces are in Newtons
%   and the parameter vector a can take two forms.  The parameter 
%   vector takes the form:
%
%        a = [a1 a2 a3 a4 a5 a6 a7 a8]
%
%   where the coefficients are as given in SAE Paper #     .
%
%   In the event that the side force cannot be physically achieved
%   with the normal forces and parameters given, the function returns
%   the value X = 1e6.
%
%   The magic formula tire model is given by:
%        Fy =  -D*sin(C*atan(B*((1-E)*x+(E/B)*atan(B*x))));
%
%   This function follows the convention that a positive slip angle
%   produces a negative side force and a negative slip angle produces
%   a positive side force.  

% Version 1.1
% Chris Gerdes 3/30/2002
% Function for use with ME106/227


% First, convert the parameter vector into something useful.
% If it has 8 elements, interpret these as coefficients a1-a8.

if (length(parvec) == 8)

% Convert the normal force in Newtons into the kN used in the paper.

    Fzi=Fzi_N/1000;
    Fzo=Fzo_N/1000;
    
    a1=parvec(1);
    a2=parvec(2);
    a3=parvec(3);
    a4=parvec(4); 
    a5=parvec(5);
    a6=parvec(6);
    a7=parvec(7);
    a8=parvec(8);

% Now calculate the values B, C, D and E for each tire.
% If the normal load is zero, set D = 0 and set the other
% parameters to some dummy values.  This will ensure that
% the tire produces zero lateral force with zero normal
% force.

    if (Fzi > 0),
        Ci=1.3;
        Di=a1*Fzi*Fzi+a2*Fzi;
        Bi=a3*sin(a4*atan(a5*Fzi))/(Ci*Di);
        Ei=a6*Fzi*Fzi+a7*Fzi+a8;
    else
        Bi=1;
        Ci=1.3;
        Di=0;
        Ei=1;
    end
    
    if (Fzo > 0),
        Co=1.3;
        Do=a1*Fzo*Fzo+a2*Fzo;
        Bo=a3*sin(a4*atan(a5*Fzo))/(Co*Do);
        Eo=a6*Fzo*Fzo+a7*Fzo+a8;
    else
        Bo=1;
        Co=1.3;
        Do=0;
        Eo=1;
    end
    
% If the parameter vector doesn't have 8 elements,
% then we have a problem...

else
    error('Incorrect parameter vector.  See help for syntax.');
end
    
% Next see what maximum side force the tires can
% produce for this normal force.  Actually, this
% finds the minimum side force which gives us the
% negative of the maximum side force.  

functofindmin=inline('Fycomb(x,P1,P2,P3,P4,P5,P6,P7,P8)',8);

[alphamin,Fymin]=fminsearch(functofindmin,0,[],Bi,Ci,Di,Ei,Bo,Co,Do,Eo);

% If the value of Fy sent to the function can be
% physically achieved, use fzero to find it.  If 
% it cannot be, return 1e6 as the value of the
% slip angle.

functomin=inline('(P1 - Fycomb(x,P2,P3,P4,P5,P6,P7,P8,P9))',9);

if (abs(Fy) < abs(Fymin))
   [alpha,Fyval]=fzero(functomin,10,[],Fy,Bi,Ci,Di,Ei,Bo,Co,Do,Eo);
else
   alpha = 1e6;
end;