function  F=Get_Fastsim(C,Cp,ab,N)
%FASTSIM求解蠕滑力程序
f=0.25; %許用粘著系數(shù)
a=ab(1);b=ab(2);%接觸斑半徑
fx=Cp(1);fy=Cp(2);fin=Cp(3);  %蠕滑力
C11=C(1);C22=C(2);C23=C(3);%kalker系數(shù)
if N<=0
    Fx=0;Fy=0;Mz=0;MU=0;
else         %法向力小于0，蠕滑力也為0
Nx=10;Ny=10;%接觸斑被分割的數(shù)目
v=0.25;
E=2.0*10^11;
G=E/2/(1+v);
Tz=pi/2;
L1=8*a/3/C11/G;L2a=8*a/3/C22/G;L2b=pi*a*sqrt(a/b)/4/C23/G;
z0=2*N/pi/a/b; 
dy=2/Ny;
Tx=0;Ty=0;Mz=0;
n1=a*fx/f/z0/L1;
n2=a*fy/f/z0/L2a;
fin1=a*b*fin/f/z0/L2b;
fin2=a*a*fin/f/z0/L2b;
for i=1:Ny
    Y=(-1+i*dy-dy/2);
    SX=n1-Y*fin1;%接觸位移失量
    px=0;py=0;%p0的初值為0
    A2=1-Y^2;
    X0=sqrt(A2);
    A=X0;
    h=2*A/Nx;%分成10塊
    AR=h*dy;%每一小塊的面積
    X=A-h/2;%
    for j=1:Nx
        SY=n2+fin2*(X);%???
        px=px-SX*h;
        py=py-SY*h;
        Z=A2-X^2;
        p=sqrt(px^2+py^2)/Z;%相差一個f？
        if p>1
            px=px/p;
            py=py/p;
        end
        Tx=Tx+px*AR;
        Ty=Ty+py*AR;
        %Mz=Mz-(px*f*z0)*(Y*b)+(py*f*z0)*(X*a);
        Mz=Mz-((px*AR)/Tz*f*N)*(Y*b)+((py*AR)/Tz*f*N)*(X*a);  %蠕滑力矩
        X0=X;
        X=X-h;
    end
end
Fx=Tx/Tz*f*N;
Fy=Ty/Tz*f*N;
MU=sqrt(Tx^2+Ty^2)/Tz*f;
if MU>=f
    MU=f;
end
%Mz=Mz/Tz;
end
F=[Fx Fy Mz MU];