Issue when converting to radians for trig functions and using lsqcurvefit ?

Hello all, I am having issues generating coefficients for a modeling formula when utilizing the lsqcurvefit function. The modeling formula is a 4-term of the Pacejka Magic Formula, therefore it utilizes sin and atan functions. The problem I am having is that when I convert the angular data to radians to solve with lsqcurvefit, then the generated coefficients are wildly off from expected. When I keep the angular data in degrees then the coefficients are close to what they should be, except for varying parameters then the model will behave unexpectedly at the origin (this may be due to the data and not the model). Take in mind that I am using sin and atan, therefore the input should be in radians !

**The angular data x1 is a 273x1 vector consisting of double type values from -11>x>11 degrees

I don't know if I'm missing something with my code as to why the trig functions & lsqcurvefit are not acting as expected. This is why I'm posting here first and not on the TTC forum.

function fy = Pacejka4_Model(P,X)
% 4 Term Pacejka Tire Model
% Solve coefficients P using Least Squares Method
x1 = X(:,1); %Slip
x1 = deg2rad(x1);
x2 = X(:,2); %FZ
D1 = P(1);
D2 = P(2);
B = P(3);
C = P(4);
D = (D1 + D2/1000.*x2).*x2; %peak value normalized
fy = D.*sin(C.*atan(B.*x1));

   inxx = 1:length(saSingle1);
xdata1 = [saSingle1(inxx),fzSingle1(inxx)];
ydata1 = fySingle1(inxx);
xdata2 = [saSingle2(inxx),fzSingle2(inxx)];
ydata2 = fySingle2(inxx);
xdata3 = [saSingle3(inxx),fzSingle3(inxx)];
ydata3 = fySingle3(inxx);
xdata4 = [saSingle4(inxx),fzSingle4(inxx)];
ydata4 = fySingle4(inxx);
xdata5 = [saSingle5(inxx),fzSingle5(inxx)];
ydata5 = fySingle5(inxx);
f0 = [1.33,.05,.1335,1.00];
options = optimset('MaxFunEvals',1000000000,...
    'TolFun',.00000001,'Display','off');
%Use the least squares method to obtain coeff for model
[FY41,f_resnorm,f_residual] = lsqcurvefit('Pacejka4_Model',f0,...
    xdata1,ydata1,[],[],options);
[FY42,f_resnorm,f_residual] = lsqcurvefit('Pacejka4_Model',f0,...
    xdata2,ydata2,[],[],options);
[FY43,f_resnorm,f_residual] = lsqcurvefit('Pacejka4_Model',f0,...
    xdata3,ydata3,[],[],options);
[FY44,f_resnorm,f_residual] = lsqcurvefit('Pacejka4_Model',f0,...
    xdata4,ydata4,[],[],options);
[FY45,f_resnorm,f_residual] = lsqcurvefit('Pacejka4_Model',f0,...
    xdata5,ydata5,[],[],options);
%Generate FY values from the model
fyModel1 = Pacejka4_Model(FY41,xdata1);
fyModel2 = Pacejka4_Model(FY42,xdata2);
fyModel3 = Pacejka4_Model(FY43,xdata3);
fyModel4 = Pacejka4_Model(FY44,xdata4);
fyModel5 = Pacejka4_Model(FY45,xdata5);

[dpb]

If your data are in degrees, why not simplify and use the XXXD() forms of the trig functions and eliminate the hassle step?