function X = ElipseContacto(i,fn,s,q) %s es de 4 comps en cada i (1=1 trasera)
global MBody
sj = s(1:2);
sk = s(3:4);
Ej = MBody.ContactoElastico(i).EI;
Nuj = MBody.ContactoElastico(i).NuI;
Ek = MBody.ContactoElastico(i).EJ;
Nuk = MBody.ContactoElastico(i).NuJ;
mu = MBody.ContactoElastico(i).mu;
[Rj,T1j,T2j,DDRs1j, DDRs2j, DDRs1s2j, Curvj] = feval(MBody.ContactoElastico(i).GeomI,q,sj,2);
[Rk,T1k,T2k,DDRs1k, DDRs2k, DDRs1s2k, Curvk] = feval(MBody.ContactoElastico(i).GeomJ,q,sk,2);
Am=-1.0864;
Bm=-0.1065;
Cm=1.35;
Dm=1.0579;
An=-0.7744;
Bn=0.2567;
Cn=0.2;
Dn=-0.281;
K1=(1-Nuj^2)/(pi*Ej);
K2=(1-Nuk^2)/(pi*Ek);
K3=0.5*(Curvk(1)+Curvk(2)+Curvj(1)+Curvj(2));
t2k=T2k/norm(T2k);
t1j=T1j/norm(T1j);
t2j=T2j/norm(T2j);
n=cross(t1j,t2j);
Psi=acos(transpose(t2k)*t2j);
K4=0.5*sqrt((Curvk(1)-Curvk(2))^2+(Curvj(1)-Curvj(2))^2+2*(Curvk(1)-Curvk(2))*(Curvj(1)-Curvj(2))*cos(2*Psi));
theta=acos(K4/K3);
mthet=Am*tan(theta-(pi/2))+(Bm/(theta^(Cm)))+Dm;
nthet=(1/(An*tan(theta-(pi/2))+1))+Bn*(theta^(Cn))+Dn*sin(theta);
aelip=norm(mthet*((3*pi*fn*(K1+K2)/(4*K3))^(1/3)));
belip=norm(nthet*((3*pi*fn*(K1+K2)/(4*K3))^(1/3)));
X = [aelip belip]; %devuelve las dimensiones de la elipse contacto
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
