This Mathematica code for the inverse kinematic resolution of a 6 dof robot was thought to work with atan2 which range is between -pi to piin Matlab but in Mathematica ArcTan works between -pi/2 to pi/2. Can you modify the code in order to correctly handle this?
Clear[InverseKinematics6DOF]
InverseKinematics6DOF[pe_,quat_, params_]:=Module[
{
(*Extract position components*) 
px,py,pz,
(*Extract quaternion components*) 
(*q0,q1,q2,q3,*)
qw,qx,qy,qz,
(*Compute rotation matrix from quaternion*) 
R,
(*End-effector orientation vectors*) 
ne,se,ae,
(*Wrist position*) 
pw,
(*Joint angles*) 
q1Sol,q2Sol,q3Sol,q4Sol,q5Sol,q6Sol,
(*Intermediate calculations*) 
L2,L3,L4,A2,Pw,
R03,
T01,
phi,beta,gamma,eta,r,a1=a1/.params,a2=a2/.params,a3=a3/.params,d1=d1/.params,d4=d4/.params,d6=d6/.params,
(*List of joint angle solutions*) 
jointAngles123,q456Sol,allSolutions
},{px,py,pz}=pe;
{qw,qx,qy,qz}=quat;
R03=Evaluate[(A0iall[[4]][[1;;3,1;;3]])/.params/.q1->q1temp/.q2->q2temp/.q3->q3temp];
T01=Evaluate[A0iall[[2]]/.params/.q1->q1temp];
(*Define your robot's DH parameters here*) 

(*Step 2:Compute the Rotation Matrix from the Quaternion*) 
R={{1-2 qy^2-2 qz^2,2 qx qy-2 qz qw,2 qx qz+2 qy qw},{2 qx qy+2 qz qw,1-2 qx^2-2 qz^2,2 qy qz-2 qx qw},{2 qx qz-2 qy qw,2 qy qz+2 qx qw,1-2 qx^2-2 qy^2}};
(*R={{1-2 (q2^2+q3^2),2 (q1 q2-q0 q3),2 (q1 q3+q0 q2)},{2 (q1 q2+q0 q3),1-2 (q1^2+q3^2),2 (q2 q3-q0 q1)},{2 (q1 q3-q0 q2),2 (q2 q3+q0 q1),1-2 (q1^2+q2^2)}};*)
(*Step 3:Compute the End-Effector Orientation Vectors*)
 ne=R[[All,1]];
se=R[[All,2]];
ae=R[[All,3]];
(*Step 4:Compute the Wrist Position*) 
pw={px,py,pz}-d6*ae;
Pw=Append[pw,1];
(*Step 5:Compute q1*) 

q1Sol=ArcTan[pw[[1]],pw[[2]]];
q1Sol={q1Sol,q1Sol+Pi};
(*Two possible solutions*) 
(*Step 6:Compute q2 and q3 for each q1*)
 L2=a2;
L3=Sqrt[d4^2+a3^2];
A2=a3;
L4=Sqrt[A2^2+d4^2];
phi=ArcCos[(A2^2+L4^2-L3^2)/(2 A2 L4)];
(*Initialize list of joint angles*) 
jointAngles123={};
Do[
(*Transform pw to frame 1*) 
T01=T01/.q1temp->q1p;
pw1=Inverse[T01].Pw;
r=Sqrt[pw1[[1]]^2+pw1[[2]]^2];
beta=ArcTan[pw1[[1]],-pw1[[2]]];
gamma=ArcCos[(L2^2+r^2-L3^2)/(2 L2 r)];
Print["beta:"];
Print[beta];
Print["gamma:"];
Print[gamma];
q2ElbowUp=(Pi/2)-beta-gamma;
q2ElbowDown=(Pi/2)-beta+gamma;
eta=ArcCos[(L2^2+L4^2-r^2)/(2 L2 L4)];
q3ElbowUp=Pi-phi-eta;
q3ElbowDown=Pi-phi+eta;
(*Append both elbow up and elbow down solutions*)
 AppendTo[jointAngles123,{q1p,q2ElbowUp,q3ElbowUp}];
AppendTo[jointAngles123,{q1p,q2ElbowDown,q3ElbowDown}];,{q1p,q1Sol}];
(*Step 7:Compute q4,q5,q6 for each set of q1,q2,q3*) 
q456Sol=Table[Module[{R01,R12,R23,R36,q4Up,q5Up,q6Up,q4Down,q5Down,q6Down,a33},{q1sol,q2sol,q3sol}=joint;
(*Compute the rotation matrix R0_3*)

R03=R03/.q1temp->q1sol/.q2temp->q2sol/.q3temp->q3sol;
R36=Transpose[R03].R;
a33=R36[[3,3]];
(*ifPositiveQ5=ArcTan[Sqrt[R36[[1,3]]^2+R36[[2,3]]^2],a33];
ifNegativeQ5=ArcTan[-Sqrt[R36[[1,3]]^2+R36[[2,3]]^2],a33];*)
ifPositiveQ5=ArcTan[a33, Sqrt[R36[[1,3]]^2+R36[[2,3]]^2]];
ifNegativeQ5=ArcTan[a33, -Sqrt[R36[[1,3]]^2+R36[[2,3]]^2]];
(*Wrist up solution*)
q5Up=ifPositiveQ5;

q4Up=ArcTan[R36[[1,3]],R36[[2,3]]];
q6Up=ArcTan[-R36[[3,1]],R36[[3,2]]];
(*Wrist down solution*) 
q5Down=ifNegativeQ5;

q4Down=ArcTan[-R36[[1,3]],-R36[[2,3]]];
q6Down=ArcTan[R36[[3,1]],-R36[[3,2]]];
{{q4Up,q5Up,q6Up},{q4Down,q5Down,q6Down}}],{joint,jointAngles123}];
(*Step 8:Combine all solutions*)
 allSolutions=Flatten[Table[{Join[jointAngles123[[i]],q456Sol[[i,1]]],Join[jointAngles123[[i]],q456Sol[[i,2]]]},{i,Length[jointAngles123]}],1];
(*Return all feasible joint angle configurations*) 
allSolutions];