coding for a finite iteration

% Specify the variable %
Temp_in = 23 ;      % Sump Temperature (°C)
d       = 0.150;    % Shaft Journal Diameter (Meter, m) 
r       = (d/2);    % Shaft Journal Radius (Meter, m)
N       = (166.69/60); % Shaft Speed (revolution per second,rps)
W       = 9000;     % Radial Load (Newton, N)
L       = 0.300 ;   % Bearing Length (Meter, m)
c       = 0.00065 ;  % Clearance (Meter, m)

%-------------------------------------------------------------------------%
% Calculate the variable %
    
C_ratio     = r/c;        % Clearance ratio
P           = W/(2*r*L);  % Nominal Pressure (Pascal)
Slenderness = L/d  ;      % Bearing Aspect Ratio

%-------------------------------------------------------------------------%
% Calculation to find average temperature and viscosity %

Tf = Temp_in + 3;               % Assumed Average Temperature (°C)
V  = 0.0018*exp(-0.026*Tf);     % Viscosity (Pa.s)
S  = ((C_ratio)^2)*((V*N)/P);   % Sommerfeld Numbert


if Slenderness      == 1
     delta_T        = (0.349109 + (6.00940*S)+(0.047467*(S^2)))*(P*(10^-6))/0.120;     % Temperature Rise(°C) at l/d=1
 elseif Slenderness == 1/2
     delta_T        = (0.394552 + (6.392527*S)-(0.036013*(S^2)))*(P*(10^-6))/0.120;    % Temperature Rise(°C) at l/d=1/2
 elseif Slenderness == 1/4
     delta_T        = 0.933828 + (6.437512*S)- (0.011048*(S^2))*(P*(10^-6))/0.120;     % Temperature Rise(°C) at l/d=1/4
else 
     delta_T        = ((33405*(S)^6)-(51557*(S)^5)+(30669*(S)^4)-(8871.9*(S)^3)+(1297.4*(S)^2)-(85.088*(S))+2.6111)*(P*(10^-6))/0.120; % Temperature Rise(°C) at l/d=Infinity
end

T_average  = Temp_in + (delta_T/2); % Average temperature (°C)


 while abs (Tf-T_average)>0.001
        
        delta = (Tf-T_average)/2;
        Tf    = Tf-delta;
    
        V  = 0.0018*exp(-0.026*Tf);   % Viscosity (Pa.s)
        S  = (C_ratio)^2*((V*N)/P);   % Sommerfeld Number
        
        if Slenderness     == 1
            delta_T        = (0.349109 + (6.00940*S)+(0.047467*(S^2)))*(P*(10^-6))/0.120;     % Temperature Rise(°C) at l/d=1
        elseif Slenderness == 1/2
            delta_T        = (0.394552 + (6.392527*S)-(0.036013*(S^2)))*(P*(10^-6))/0.120;    % Temperature Rise(°C) at l/d=1/2
        elseif Slenderness == 1/4
            delta_T        = (0.933828 + (6.437512*S)- (0.011048*(S^2)))*(P*(10^-6))/0.120;    % Temperature Rise(°C) at l/d=1/4
        elseif Slenderness > 1
            delta_T        = ((33405*(S)^6)-(51557*(S)^5)+(30669*(S)^4)-(8871.9*(S)^3)+(1297.4*(S)^2)-(85.088*(S))+2.6111)*(P*(10^-6))/0.120; % Temperature Rise(°C) at l/d=Infinity
        end
    
         T_average  = Temp_in + (delta_T/2); % Average temperature (°C)end
    
     
    disp ('Table for iteration');disp ('         Tf          V            S          deltaT       T average ');
    format shortG
    disp ([Tf', V',S',delta_T', T_average',]);
  
 end

I want to get the value of clearance when Tf-Taverage=0 for only when the iteration stop at the 36th iteration. Can somebody help me with this?Thank you.

[Walter Roberson]

What leads you to think that Tf-Taverage == 0 will ever hold? You are using floating point calculations and there will be round-off issues. It is entirely possible that evaluating at two adjacent representable numbers will lead to values that are negative for one of the two choices, and positive for the other choice, without it being possible to find an exact 0.

For example, if you were trying to solve x^2-2 == 0 then

>> (sqrt(2)*(1-eps))^2-2
ans =
     -4.44089209850063e-16
>> (sqrt(2)*(1))^2-2
ans =
      4.44089209850063e-16
>> (sqrt(2)*(1+eps))^2-2
ans =
      8.88178419700125e-16

No exact solution