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 first is for
%   tire force calculation assuming zero camber angle.  In this case
%
%        a = [a1 a2 a3 a4 a5 a6 a7 a8]
%
%   where the coefficients are as given in SAE Paper #870421.  The 
%   other form is
%
%        a = [a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 cain caout]
%
%   where the coefficients are again taking from the referenced
%   paper and the angles cain and caout represent the inner and outer
%   tire inclination angles in degrees.
%
%   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.2
% Chris Gerdes 
% Dynamic Design Lab
% Stanford University
% 4/6/2002
% Function for use with ME106/227
% This revised version also includes camber angle change.

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

if ((length(parvec) == 8)|(length(parvec)==14))

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

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

% Get the remaining values if camber coefficients and angles
% were also passed to the function.

    if (length(parvec) == 14)
       a9=parvec(9);
       a10=parvec(10);
       a11=parvec(11);
       a12=parvec(12);
       camberi=parvec(13);  % Camber angle of inside tire
       cambero=parvec(14);  % Camber angle of outside tire
    end
    
% Now calculate the values B, C, D, E, DSh and DSv for each 
% tire.

    Ci=1.3;
    Di=a1*Fzi*Fzi+a2*Fzi;
    if (Di ~= 0) 
       Bi=a3*sin(a4*atan(a5*Fzi))/(Ci*Di);
    else
       Bi = 1;
    end
    Ei=a6*Fzi*Fzi+a7*Fzi+a8;
    
    Co=1.3;
    Do=a1*Fzo*Fzo+a2*Fzo;
    if (Do ~= 0) 
       Bo=a3*sin(a4*atan(a5*Fzo))/(Co*Do);
    else
       Bo = 1;
    end
    Eo=a6*Fzo*Fzo+a7*Fzo+a8;

    if (length(parvec) == 14)
       Bi=Bi-a12*abs(camberi)*Bi;
       Bo=Bo-a12*abs(cambero)*Bo;
       DShi=a9*camberi;
       DSho=a9*cambero;
       DSvi=(a10*Fzi*Fzi + a11*Fzi)*camberi;
       DSvo=(a10*Fzo*Fzo + a11*Fzo)*cambero;
       
       if (Fzi <= 0)
          DSvi = 0;
       end
       if (Fzo <= DSvo)
          DSvo = 0;
       end
       
    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.  If camber angle
% is included we need to pay attention to the sign
% of the side force since the curve is no longer
% symmetric and the max and min forces are distinct.  

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

   [alphamin,Fymax]=fminsearch(functofindmin,0,[],Bi,Ci,Di,Ei,Bo,Co,Do,Eo);
else
   if (Fy <= 0) 
      functofindmin=inline('Fycomb(x,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12)',12);
   else
      functofindmin=inline('-Fycomb(x,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12)',12); 
   end
   
   [alphamax,Fymax]=fminsearch(functofindmin,0,[],Bi,Ci,Di,Ei,DShi,DSvi,Bo,Co,Do,Eo,DSho,DSvo);
end;

% 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.

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

   if (abs(Fy) < abs(Fymax))
      [alpha,Fyval]=fzero(functomin,alphamax*sign(Fy),[],Fy,Bi,Ci,Di,Ei,Bo,Co,Do,Eo);
   else
      alpha = -sign(Fy)*1e6;
   end;
else
   functomin=inline('(P1 - Fycomb(x,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13))',13);

   if (abs(Fy) < abs(Fymax))
      [alpha,Fyval]=fzero(functomin,alphamax*sign(Fy),[],Fy,Bi,Ci,Di,Ei,DShi,DSvi,Bo,Co,Do,Eo,DSho,DSvo);
   else
      alpha = -sign(Fy)*1e6;
   end;
end;