function [x,y,z,alpha,beta,gamma] = fcn(theta1,theta2,theta3,dtheta1,dtheta2,dtheta3)

La=0.069;
Lc=0.33	;
Ra=0.169;
Rc=0.08;
b=0.049;

%Centre des spheres

    sphere1=[Ra - Rc + La*cos(theta1) ; 0 ; La*sin(theta1)];
    sphere2=[Rc/2 - Ra/2 - (La*cos(theta2))/2 ; (3^(1/2)*Rc)/2 - (3^(1/2)*(Ra + La*cos(theta2)))/2 ; La*sin(theta2)];
    sphere3=[Rc/2 - Ra/2 - (La*cos(theta3))/2 ; (3^(1/2)*(Ra + La*cos(theta3)))/2 - (3^(1/2)*Rc)/2 ; La*sin(theta3)];
                                 
%Determination de l'intersection des 3 spheres de centres sphere1,2,3 et de
%rayon Lc 

    dists1s2 = sphere2-sphere1; 
    e_x = dists1s2/norm(dists1s2); %u1
    dists3s1 =sphere3-sphere1; 
    i = dot(e_x,dists3s1); %scalaire dist et vecteur unitaire
    temp3 = dists3s1 - i*e_x; 
    e_y =temp3/norm(temp3); %u2
    e_z = cross(e_x,e_y); %u3
    d = norm(sphere2-sphere1); 
    j= dot(e_y,dists3s1); 
    x = d / 2; 
    y = (-2*i*x + i*i + j*j) / (2*j); 
    temp4 = Lc*Lc - x*x - y*y; 
    z =sqrt(temp4); 
    p_1 = +sphere1+x*e_x + y*e_y -z*e_z; 
    x=p_1(1);
    y=p_1(2);
    z=p_1(3);
    
    alpha=0;
    beta=0;
    gamma=0;
end