%solve the pressure of bearing
%d1=input('please input the value of diameter of bearing')
%d2=input('please input the value of diameter of shaft')
%e=input('please input the value of the rate of off-center')
%l=input('please input the value of shaft length')
%c=(d1-d2)/2
%t=d2/l
m=25
n=25
error=1.0e-2
T=[0.1,0.2,0.4,0.6,0.8,1.0,1.5,2.0]
W=[0.1,0.2,0.3,0.4,0.6,0.7,0.8,0.9,0.95,0.99]
p=zeros([m+1,n+1])
q=(2*pi)/m
v=2/n
k=2
F=zeros([8,10])
for jj=1:8
  for ii=1:10 
  t=T(jj)
  w=W(ii)
  while k>0
    p1=p;
     for j=2:n
       for i=2:m
        a=(1+w*cos((pi/180)*360*(2*i+1)/(2*m)))^3
        b=(1+w*cos((pi/180)*360*(2*i-1)/(2*m)))^3
        c=((t*q)/v)^2*(1+w*cos((pi/180)*360*(2*j-1)/(2*m)))^3
        d=((t*q)/v)^2*(1+w*cos((pi/180)*360*(2*j+1)/(2*m)))^3
        e1=a+b+c+d
        f=q*((1+w*cos((pi/180)*(360*(2*i+1))/(2*m)))-(1+w*cos((pi/180)*(360*(2*i-1))/(2*m))))
        p(i,j)=(a*p(i+1,j)+b*p(i-1,j)+c*p(i,j+1)+d*p(i,j-1)-f)/e1
           if p(i,j)<=0 
             p(i,j)=0
             break;
           end
      end
     end
    sum=0;
    sum1=0;
    sum2=0;
    for s=1:m+1
      for t=1:n+1
         sum1=sum1+abs((p(s,t)-p1(s,t))) %收斂準則的分子
         sum2=sum2+abs(p(s,t))     %收斂準則的分母
      end
    end
  sum=sum1/sum2;
    if sum<=error
      break;
    end
  end
  F1=0
  F2=0
  for j=1:n+1
    for i=1:m+1
       F1=F1+p(i,j)*cos(i*q)*q*v
       F2=F2+p(i,j)*sin(i*q)*q*v
    end
  end
  F(jj,ii)=sqrt(F1^2+F2^2)
  end
  end
plot(W,F)