function [x,y]=trajectory(A,O,PD,PHI0,N,draw_plot)

#-- Robot parameters
r=30;    #-- Wheels radius (mm)
l=85/2;  #-- Distance between the com and the wheel (mm)

#-- Oscillator's parameters
A= A   * pi/180;
O= O    * pi/180;
T= 30;
PD= PD  * pi/180;
PHI0= PHI0  * pi/180;

tics = T*N+1;

n=[0:tics];

S1 = sin(2*pi*n/T + PHI0);
S2 = sin(2*pi*n/T + PD + PHI0);

C1 = cos(2*pi*n/T + PHI0);
C2 = cos(2*pi*n/T + PD + PHI0);


#-- These equations are correct! Validated!!
teta = r/(2*l) * (A*S1 - A*S2 + 2*O) + pi/2;  #-- correction + 0.06
dx = r*pi*A/T * (C1 + C2) .* cos(teta);
dy = r*pi*A/T * (C1 + C2) .* sin(teta);

#-- Integrate dx and dy to obtain the x and y position
x=zeros(1,tics);
y=zeros(1,tics);
for t=[2:tics]
  x(t) = x(t-1)+dx(t-1);
  y(t) = y(t-1)+dy(t-1);
endfor;

if (draw_plot==true)
  plot(x,y);
  axis(["square"])
  axis([0,500,-250,250])
  xlabel ("x (mm)");
  ylabel ("y (mm)");
  title("Robot trayectory (A=90, PD=90, PHI0=0)");
  grid on;
endif;



