function dx = Contact_friction(t,x,w)

% Parameters
mb1 = 4;m1= 32.1;
c1 = 1050; c2 = 2100; k = 2.5*10^7; 
e = 0.05*1e-3;
g  = 9.81;
c = 0.11*1e-3;

R =25*1e-3;
L = 12*1e-3;
mu = 0.018;
P = 16.05;
dentak1=0.1*k;
dentak2=0.1*k;

beta=pi/2;
A=0;
delta1 = mu*w*R*L/P*(R/c)^2*(L/2/R)^2;

%油膜力摩擦系数
% Oil force
% [fx,fy] = oil_force(xb1,yb1,dxb1,dyb1)
[fx1,fy1] = oil_force(x(1),x(3),x(2),x(4));

phi=t+beta;
liewen=((1+cos(phi))/2)^A;

% 裂纹刚度
kxx=k-liewen*(dentak1*cos(phi)*cos(phi)+dentak2*sin(phi)*sin(phi));
kyx=-liewen*((dentak1-dentak2)*sin(phi)*cos(phi));
kxy=-liewen*((dentak1-dentak2)*sin(phi)*cos(phi));
kyy=k-liewen*(dentak1*cos(phi)*cos(phi)+dentak2*sin(phi)*sin(phi));

%方程式
dx(1)=x(2);
fs1 = -c1*x(2)/mb1/w+0.5*kxx*(x(1)-x(5))/mb1/w^2+0.5*kxy*(x(3)-x(7))/mb1/w^2+delta1*P*fx1/mb1/c/w^2;
dx(3)=x(4);
fs2 = -c1*x(4)/mb1/w+0.5*kyx*(x(1)-x(5))/mb1/w^2+0.5*kyx*(x(3)-x(7))/mb1/w^2+delta1*P*fy1/mb1/c/w^2-g/c/w^2;
dx(5)=x(6);
fs3 = -c2*x(6)/m1/w-kxx*(x(1)-x(5))/m1/w^2-kxy*(x(3)-x(7))/m1/w^2+e*cos(t-beta)/c;
dx(7)=x(8);
fs4 = -c2*x(8)/m1/w-kyx*(x(1)-x(5))/m1/w^2-kyy*(x(3)-x(7))/m1/w^2+e*sin(t-beta)/c-g/c/w^2;

dx = [x(2);
     fs1;
    x(4);
     fs2;
    x(6);
     fs3;
    x(8);
     fs4];
    
end
    


function [fx1,fy1] = oil_force(xb1,yb1,dxb1,dyb1)

CC1 = -sqrt((xb1-2*dyb1)^2+(yb1+2*dxb1)^2)/(1-xb1^2-yb1^2);

alpha= atan((yb1+2*dxb1)/(xb1-2*dyb1))-pi/2*sign((yb1+2*dxb1)/(xb1-2*dyb1))-pi/2*sign(yb1+2*dxb1);
Gxya= 2*(pi/2+atan((yb1*cos(alpha)-xb1*sin(alpha))/(sqrt(1-xb1^2-yb1^2))))/(sqrt(1-xb1^2-yb1^2));
V     = (2+(yb1*cos(alpha)-xb1*sin(alpha))*Gxya)/(1-xb1^2-yb1^2);
S     = (xb1*cos(alpha)+yb1*sin(alpha))/(1-(xb1*cos(alpha)+yb1*sin(alpha))^2);

fx1 = CC1*(3*xb1*V-sin(alpha)*Gxya-2*cos(alpha)*S);
fy1 = CC1*(3*yb1*V+cos(alpha)*Gxya-2*sin(alpha)*S);


clear;clc;
w= 200;
T = 2*pi;
x0 = ones(8,1)*0.1;
% x0 = zeros(8,1)*0.1;

[t,x]=ode45(@Contact_friction,[0,T*250],x0,[],w);
x0 = x(end,:)     ;
w = linspace(200,2500,50);
for h = 1:length(w)
    [t,x]=ode45(@Contact_friction,[0:T/500:T*250],x0,[],w(h));
    plot(w(h),x(200*500:500:end,5),'k.');hold on;
   
%     xrms(h) =
x0=x(end,:)   ;
    h
end
    
set(gcf,'PaperPositionMode','manual');
set(gcf,'PaperUnits','points');
xx=get(gcf,'position');
set(gcf,'PaperPosition',[0,0,xx(3)/1,xx(4)/1.5]);
print(gcf,'-dtiff','-r600',['E:\HUNDUN'])
print(gcf,'-deps','-r600',['E:\HUNDUN'])