function [ecc_dim,attitude_angle] = JBe(R,L,C,Ns,load,mu_in)

% R è il raggio dell'albero in [m]
% L è la lunghezza del bearing [m]
% C è la clearance [m]
% Ns numero di giri al secondo
% load è il carico [N]
% mu_in è la viscosità in [Pa*s]

Epsilon=[0.05:0.05:1];
output=zeros(length(Epsilon),9);

for n=1:length(Epsilon)
epsilon=Epsilon(n);
%   epsilon=0.1;

% R=0.051;        % Raggio dell'albero
% L=0.0765;
D=2*R;
lambda=L/D; % Aspect Ratio
% C=6.375*10^-5;
% Ns=25;
% load=11000;
inlet_spec=2;

if inlet_spec==1 %inlet hole

    Dh=0;
    
elseif inlet_spec==2   %axial groove unloaded section

    Lg=0.067;
    wg=4.76*10^-3;

elseif inlet_spec==3   %axial groove at maximum film thickness

    Lg=0.067;
    wg=4.76*10^-3;

elseif inlet_spec==4   %Two axial groove at +-90°

    Lg=0.067;
    wg=4.76*10^-3;

end

% mu_in=0.03;         % lubricant viscosity
rho=0.0315;         % lubricant density
beta=0.0414;
Cp=0.48;
PS=2.76*10^5;       % lubricant supply pressure
Tin=36.8;           % lubricant inlet temperature

N=70;               % numero di divisioni lungo theta
M=70;               % numero di divisioni lungo y

p_in=0;
p_out=0;

p_in_adim=p_in/(mu_in*Ns)*((C/R)^2);
p_out_adim=p_out/(mu_in*Ns)*((C/R)^2);

vp_in=p_in_adim*ones(M,1);
vp_out=p_out_adim*ones(M,1);


delta_theta=2*pi/(N-1);    % intervallo angolare
delta_y=1/(M-1);           % intervallo lungo l'asse

theta=zeros(N,1);

for i=1:N
    theta(i)=delta_theta*i-delta_theta;
end


Y=zeros(M,1);

for k=0:length(Y)-1
Y(k+1)=(k*delta_y);
end

P=zeros(N,M);   % allocazione della matrice delle pressioni
omega=1.95;
tol=10^(-3);
k=0;
err=1;
norma=0;
sump=0;

h=zeros(N,1);
hd=zeros(N,1);

for i=1:length(theta) 

    h(i)=1+epsilon*cos(theta(i));
    hd(i)=-epsilon*sin(theta(i));

end

P(1,:)=vp_in;
P(N,:)=vp_out;

while err>tol
 for i=2:N-1
  for j=2:M-1
   
     
 k=k+1; %contatore
     
     a=2*((lambda^-2)*0.25/(delta_y^2)+(1/(delta_theta)^2));
     b=1/(delta_theta^2)+3*(hd(i))/((h(i))*(2*delta_theta));
     c=1/(delta_theta^2)-3*(hd(i))/((h(i))*(2*delta_theta));
     d=(lambda^-2)*0.25/(delta_y^2);
     e=-(12*pi*hd(i))/((h(i))^3);
    

    P_prec=P(i,j);     
    P(i,j)=P(i+1,j)*(b/a)+P(i-1,j)*(c/a)+P(i,j-1)*(d/a)+P(i,j+1)*(d/a)+(e/a); 
    P(i,j)=P_prec*(1-omega)+omega*P(i,j);
 
    if P(i,j)< 0
       
            P(i,j)=0;
    end
    
   norma=norma+abs(P(i,j)-P_prec);
   sump=sump+P(i,j); 
    
   end
 end

  err=norma/sump;
       
end
  
Theta_grad=theta*360/(2*pi);     % espressione della coordinata angolare in gradi
 
% surf(Theta_grad,Y,P')           %generazione superficie di pressione
% 
% plot(Theta_grad,P(:,M/2))

Pmax=max(P(:,M/2));
[Pmax1,I]=max(P(:,M/2),[],'linear');
theta_max=theta(I)*360/(2*pi);
[Pmin1,J]=min(P(2:N,M/2),[],'linear');
theta_cav=theta(J)*360/(2*pi);

for i=1:length(theta)
    v(i)=cos(theta(i));
    w(i)=sin(theta(i));
end

for i=1:M
    z(:,i)=v'.*P(:,i);
    t(:,i)=w'.*P(:,i);
end

i1=trapz(Y,trapz(theta,z));
i2=-trapz(Y,trapz(theta,t));

SOM=((i1^2)+(i2^2))^(-0.5)*2; %numero di sommerfeld

phi=atan(i2/i1)*360/(2*pi);

%  P_dim=P*mu*Ns*(R/C)^2;

pp=(h.^3).*(P(:,M-1)/(delta_y));
x=trapz(theta,pp);

QL=2/(12*pi)*(R^2/L^2)*x;     % Portata di leakage adimensionale

coeff=(-1/(pi*12)*h(2)^3)/delta_theta;
oo=coeff*(P(2,:)-P(2-1,:))+h(2);
xx=trapz(Y,oo);
QI=xx;                       % Portata di inlet adimensionale


%  psegn=(P*mu*Ns*R^3*L)/(C^2);

  psegn=P;

W1=trapz(Y,trapz(theta,psegn));
W=((i1^2)+(i2^2))^(0.5);
ww=(W/W1)^-1;

%  SOM=mu*Ns*D*L*(R/C)/W;

nn=zeros(N,M);

for j=1:M
 for i=1:N-1
     
    nn(i,j)=D*L*mu_in*Ns*(R/C)*pi*(((P(i+1,j)-P(i,j))/delta_theta)*h(i)/(4*pi)+1/h(i));
          
 end
end

for j=1:M
nn(N,j)=D*L*mu_in*Ns*(R/C)*pi/h(N);
end

f_mu=trapz(Y,trapz(theta,nn));

f=f_mu/(W*mu_in*Ns*R^3*L/(C^2));
factor=R/C*f;

output(n,1)=epsilon;
output(n,2)=SOM;
output(n,3)=QL;
output(n,4)=QI;
output(n,5)=factor;
output(n,6)=Pmax;
output(n,7)=theta_max;
output(n,8)=phi;
output(n,9)=theta_cav;

end

% toc
% tic
% plot(output(:,1),output(:,2))

% surf(Theta_grad,Y,P')
% plot(Theta_grad,P(:,M/2))

% Tout_i=1;
% mu_i=mu_in;
% Tout=0;
kk=1;
% diff=10;
Qm=0;
Qp=0;
S_exponential=0;

% if inlet_spec==1
% 
% flow_coeff=(Qm^S_exponential)*Qp^(1-S_exponential);
% 
%    elseif inlet_spec==2 
% 
%        m=0.27*(leakage_flow/Qp)^0.27;
%        flow_coeff=Qm*(Lg/L)^m;
% 
% end


% while  abs(diff)>1
% 
%     
% kk=kk+1; %contatore ciclo 



mu_i=mu_in; %%%%
S=(mu_i*Ns*L*D*(R/C)^2)/load;


        
for i=1:length(Epsilon)-1

    if S<output(i,2)
    if S>output(i+1,2)

        delta_S=output(i,2)-output(i+1,2);
        x=(S-output(i+1,2))/delta_S;
    
        ecc=(output(i,1)*x+output(i+1,1)*(1-x));
        attitude_angle=(output(i,8)*x+output(i+1,8)*(1-x));
        cavitation_angle=(output(i,9)*(x)+output(i+1,9)*(1-x));
        Adim_maximum_Pressure=(output(i,6)*(x)+output(i+1,6)*(1-x));
        friction_factor=(output(i,5)*(x)+output(i+1,5)*(1-x));
        adim_leakageflow=(output(i,3)*(x)+output(i+1,3)*(1-x));
        adim_flow=(output(i,4)*(x)+output(i+1,4)*(1-x));
    
        break

    end
    end
end

friction_coefficient=friction_factor/(R/C);
Friction_force=friction_coefficient*load;
Power_loss=Friction_force*(2*pi*R*Ns);
leakage_flow=adim_leakageflow*(pi/2)*Ns*D*L*C;
% Qp=(f_g*((h_g/C)^3)*PS*C^3)/mu_i;
% Qm=leakage_flow+Qp-0.3*(leakage_flow*Qp)^0.5;

if inlet_spec==1 %oil hole in the unloaded section of the bearing

f_g=0.675*(Dh/L+0.4)^1.75;
Qp=(f_g*((1+ecc*cos(phi))^3)*PS*C^3)/mu_i;
Qm=leakage_flow+Qp-0.3*(leakage_flow*Qp)^0.5;    
S_exponential=0.6;
flow_coeff=(Qm^S_exponential)*Qp^(1-S_exponential);

elseif inlet_spec==2 % axial groove positioned at the maximum film thickness

f_g=(1.25-0.25*(Lg/L))/(3*(L/Lg-1)^(1/3))+(wg/D)/(3*(L/D)*(1-Lg/L));
Qp=(f_g*((1+ecc*cos(phi))^3)*PS*C^3)/mu_i;
Qm=leakage_flow+Qp-0.3*(leakage_flow*Qp)^0.5; 
S_exponential=0.4+0.75*(Lg/L)^0.7;
flow_coeff=(Qm^S_exponential)*Qp^(1-S_exponential);

elseif inlet_spec==3

f_g=(1.25-0.25*(Lg/L))/(3*(L/Lg-1)^(1/3))+(wg/D)/(3*(L/D)*(1-Lg/L));
Qp=(f_g*((1+ecc)^3)*PS*C^3)/mu_i;
Qm=leakage_flow+Qp-0.3*(leakage_flow*Qp)^0.5;
m=0.27*(leakage_flow/Qp)^0.27;
flow_coeff=Qm*(Lg/L)^m;


elseif inlet_spec==4 

f_g=(1.25-0.25*(Lg/L))/(3*(L/Lg-1)^(1/3))+(wg/D)/(3*(L/D)*(1-Lg/L));
Qp=(f_g*((1+ecc*cos(phi))^3+(1-ecc*sin(phi))^3)*PS*C^3)/mu_i;
Qm=leakage_flow+Qp-0.3*(leakage_flow*Qp)^0.5
m=0.27*(leakage_flow/Qp)^0.27;
flow_coeff=Qm*(Lg/L)^m;

end

leakage_flow_total=flow_coeff;
temperature_rise=Power_loss/((1.72*10^6)*4.186*leakage_flow_total);
Maximum_pressure=(Adim_maximum_Pressure*mu_i*Ns*(R/C)^2)+PS;

% Tout=Tout_i;
Tout_i=Tin+temperature_rise;
% mu=mu_in;
% mu_i=mu_in*exp(-beta*(Tout_i-Tin));
% 
% diff=Tout_i-Tout;


% risultati

risultati(1,kk)=ecc;
risultati(2,kk)=attitude_angle+180;
risultati(3,kk)=cavitation_angle;
risultati(4,kk)=Maximum_pressure;
risultati(5,kk)=Friction_force;
risultati(6,kk)=Power_loss;
risultati(7,kk)=leakage_flow_total;
risultati(8,kk)=temperature_rise;
risultati(9,kk)=Tout_i;
risultati(10,kk)=S;

ecc_dim=ecc*C;

%     if kk>=1000
% 
% return
%     end
% end

% toc

end