clear all
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GLOBAL INPUTS TO BE CROSS_USED IN OTHER FUNCTIONS
global pcav 
global rhocav 
global Z
global w
global omega
global TC
global LC
global DR
global pa
global pD
global pm
global BM
global mu 
global BC
global GArea
global Gangle

tic 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%  USER INPUT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Nominal gear radial and axial clearance, bearing block clearance
TC = 0.05; %define tooth tip clearance with casing [mm]
LC = 0.01; %define side tooth clearance with casing [mm]
BC = 0.0; %define bearing block clearance with casing sides [mm]

DR = 17.78 * 0.8; %define gear drainage radius [mm]

%Backend groove specs: groove are and angular extension
GArea = 6.93; %define groove area [mm^2]
Gex = [12, 25]; %define groove extension limits 
                %from horizontal [deg] [(+)if going down]

%Hydrodynamic bearings specs (see dedicated model)
AR = 0.90;  %define bearing Aspect Ratio [mm]
BearR = DR; %define bearing internal radius [mm]
Clear = 0.04; %define bearing internal clearance [mm]

%Working conditions: pressure and speed
RPM = 1000; % speed [rpm]
pa = 0e5; %[Pa]
pm = 10e5; %[Pa]

%Fluid properties
BM = 1.635e9; %Bulk Modulus [Pa]
pcav = 0e5; %Cavitation [Pa]
rhocav = 860; %[kg/m3] 
rho = rhocav; %define fluid density
mu = 36.5e-3; %define dynamic viscosity [Pa*s]

%Settings for iterative method and convergence condition
max_count = 20; %max iterations for predicting Static Equilibrium Position (SEP)
                %(max_count =10÷15 should be enough to reach desired
                %tolerances)
toll_e = 1e-03; %relative tolerance on eccentricity module, e
toll_G = 1e-03; %relative tolerance on eccentricity direction, Gamma
%NOTE: Keep tolerances around 1e-03÷1e-06 for best tradeoff of accuracy and
%computational time!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pD = pa; %drainage is connected to suction, same pressure [Pa]
Gangle = Gex + 180; %backend groove angular coord.
                    %in driver gear reference system 
omega = RPM * 2*pi/60; %from [rpm] to [rad/sec]
Lbear = AR*2*BearR; %bearing length given radius and aspect radius [mm]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% READ INPUT from .xlsx FILE %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Read data from Excel
filename = 'Gprofile.xlsx';
M = readmatrix(filename,'Sheet', 1, 'Range','A2:E2');
Z = M(1); %teeth number
w = M(2); %gear thickness [mm]
alpha_w = M(3); %work pressure angle [deg]
Rbase = M(5); %base radius [mm]
M = readmatrix(filename,'Sheet', 1, 'Range','A4:F4');
TH = M(1); %tooth height [mm]
TW_r = M(2); %tooth width at root radius [mm]
TW_t = M(3); %tooth width at tip radius [mm]
VV = M(5); %inter teeth vane volume [mm3]
MT = M(6); % mean tooth thickness [mm] (for tooth side leakages)
M = readmatrix(filename,'Sheet', 2);
TR = M(:, 1); TW = M(:, 2); %tooth profile radii and thickness vectors [mm]

%Read data from Excel
filename = 'GPgeometry.xlsx'; 
M =  readmatrix(filename,'Sheet', 1, 'Range', 'A2');
TETA = M(:, 1); %sampling angular positions vector [deg]
TV = M(:, 2); %trapped volume vector [mm3]
dVol = M(:, 3); %trapped dvolume derivative
UDA = M(:, 4); %Upper Discharge Area (connection area w/ discharge relief groove [mm2]
LDA = M(:, 5); %Lower Discharge Area (connection area w/ suction relief groove [mm2]
NoSV = M(:, 6); %Number of (contemporary) Sealed Volumes vector 
thcontact1 = M(:, 7);  %angular coord. of meshing contact #1 [deg] (driving gear ref. sys)
thcontact2 = M(:, 8);  %angular coord. of meshing contact #2 [deg] (driving gear ref. sys)
rcontact1 = M(:, 9);  %radial coord. of meshing contact #1 [deg] (driving gear ref. sys)
rcontact2 = M(:, 10); %radial coord. of meshing contact #2 [deg] (driving gear ref. sys)
thseal1 = M(:, 11);  %angular coord. of last tooth in contact w/ casing [deg]
thseal2 = M(:, 12);  %angular coord. of last tooth in contact w/ casing [deg]
r2contact1 = M(:, 13); %radial coord. of meshing contact #1 [deg] (driven gear ref. sys)
r2contact2 = M(:, 14); %radial coord. of meshing contact #2 [deg] (driven gear ref. sys)

%Cropping trapped volume vectors when it does not exist
dVol(isnan(dVol)) = []; 
for ii = 2 : length(thcontact2)
    if thcontact2(ii-1) == 10000
        rcontact2(ii) = 10000; thcontact2(ii) = 10000;  r2contact2(ii) = 10000;
    end
end
dVol(isnan(dVol)) = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Interpolate trapped volume, UDA, LDA data from Excel vectors to make
%continuous functions. 
%Trapped Volume is a second degree polynomial,
%derivative is a first degree polinomial (polyfit & polyval)
cdvol = polyfit( TETA(1 : length(dVol)), dVol, 1);
cVol = polyfit( TETA(1 : length(dVol)), TV(1: length(dVol)), 2);
dVdth = @(t) polyval(cdvol, t);
VolTr = @(t) polyval(cVol, t);
%UDA, LDA interpolation w/ makima interpolating MATLAB function (the best
%to interpolate data w/ very smooth curves)
UDA = UDA(1 : length(dVol)); LDA = LDA(1 : length(dVol));
th = TETA(1 : length(dVol));
Adel = @(t) makima(th, UDA, t); 
Asuc = @(t) makima(th, LDA, t);

% De-comment if you want to plot UDA, LDA, Trapped volume and its derivative
%   figure(10)
%   plot(th, UDA, '-r', th, LDA, '-b', th, Adel(th), '*r', th, Asuc(th), '.b' )   
%   figure(20)
%   plot(th, TV(1:length(th)), th, dVol , th, VolTr(th), '.',  th, dVdth(th), '.' )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%TRAPPED VOLUME - PRESSURE EVOLUTION %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%In this section Trapped Volume ODE is set and solved.

%Define main leakage flows in trapped volume: delivery relief groove
%(Qtin), suction relief groove (Qtout), drainage (Qtdr),
%tooth sides at delivery (Qtlin), tooth sides at suction (Qtlout). 
% All are in [mm3/sec]!!!

Qtin = @(t, p) (0.65*Adel(t).*(1e3*sqrt(2*abs(pm - p)./density(pm,BM,pcav,rhocav)))).*(pm>=p) + ...
               (-0.65*Adel(t).*(1e3*sqrt(2*abs(pm - p)./density(p,BM,pcav,rhocav)))).*(pm<p);
                                                   %[mm3/s] !!!!!!!!!!!!!!
Qtout = @(t, p) (0.65*Asuc(t).*(1e3*sqrt(2*abs(p-pa)./density(p,BM,pcav,rhocav)))).*(p>=pa) + ... 
                (-0.65*Asuc(t).*(1e3*sqrt(2*abs(p-pa)./density(pa,BM,pcav,rhocav)))).*(p<pa);    
                                                    %[mm3/s] !!!!!!!!!!!!!
Qtdr = @(t, p) (.5*TR(1)*2*pi/Z)*LC.^3/(12*mu*(TR(1)-DR))*(p-pD);
                                                   %[mm3/s] !!!!!!!!!!!!!!
Qtlin = @(t, p) TH*LC.^3/(12*mu*MT)*(pm - p) + .5*(TH*LC*omega*Rp);
                                                   %[mm3/s] !!!!!!!!!!!!!!
Qtlout = @(t, p)  TH*LC.^3/(12*mu*MT)*(p - pa) +.5*(TH*LC*omega*Rp);
                                                   %[mm3/s] !!!!!!!!!!!!!!

%Define Trapped Volume ODE
funTrap = @(t, p) pi/180 * (BM./VolTr(t)).*((2*Qtin(t, p) - 2*Qtout(t, p) - ...
                            4*Qtdr(t, p) + 2*Qtlin(t, p) - 2*Qtlout(t, p))./omega - dVdth(t));

options = odeset('Events', @pressureEventsFcn);

   figure(100)
   %[tt, pt, tend, pte, ie] = ode15s(funTrap, [0, th(end)], pm , options);
   [tt, pt] = ode15s(funTrap, [0, th(end)], pm );
   plot(tt, pt)
   xlabel('\theta [deg]')
   ylabel('p_{trap}/p_m')
   hold on
   
   if tt(end) < th(end)
         plot( [tt(end), th(end)], [pcav, pcav], '-r')
   end
   
   hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% SEALED VOLUMES - PRESSURE EVOLUTION %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%In this section Sealed Volume ODE set is solved for a driver gear STATIC
%EQUILIBRIUM POSITION guess until convergence tolerance is achieved.
%Main leakages and ODE set will be defined in called function 'odesealfcn.m'.

%INITIAL GUESS ON BEARING STATIC EQUILIBRIUM POSITION (SEP)
xbear = 0.000; ybear = -0.8*TC;  %initial guess on SEP [mm]

%From (X,Y) coord. to (e,Gamma) coord.
BEAR = [xbear, ybear, sqrt(xbear^2+ybear^2), rad2deg(mod(atan2(ybear, xbear), 2*pi))];

%Write SEP initial guess on Excel FILE
filename = 'BearingCenterCoord.xlsx';
string = ["xbear" "ybear" "e" "GAMMA"]; 
writematrix(string, filename, 'Sheet', 1, 'Range', 'A1')
writematrix(BEAR, filename, 'Sheet', 1, 'Range', 'A2')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Two possible schemes for NoSV: either N->N+1->N or N->N-1->N.
%following code lines are for detecting the type of scheme. 
%When a change in NoSV values occurs the corresponding TETA index
%is registered in 'index' vector.
%Scheme type is in vector 'nsv' (length = 3) 
cont = 0; index = -ones(10, 1); nsv = -ones(10, 1);
for ii = 1 : length(NoSV) - 1
    if NoSV(ii) ~=  NoSV(ii+1)
        cont = cont + 1;
        index(cont) = ii + 1;
        nsv(cont) = NoSV(ii);
    end 
    if ii == length(NoSV) - 1
        cont = cont + 1;
        index(cont) = ii + 1;
        nsv(cont) = NoSV(ii);
    end
end 
index(index==-1) = []; nsv(nsv==-1) = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Some initialized values (please do not modify)
tseal = []; pseal = [];
iter = 1;
rel_err_e = 10000; %relative error on eccentricity magnitude
rel_err_G = 10000; %relative error on eccentricity direction angle

%Set how many times to re-run ode solver over angular pitch angle. You may 
%change this (however, repeat = 5 lets curve converge to regime w/out 
%useless overcycling)
repeat = 5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Start iterative process to evaluate SEP
while iter <= max_count && (rel_err_e >= toll_e || rel_err_G >= toll_G)

%following 'for cycle jj' is to determine Sealed Volumes pressure evolution @ regime.
%tseal will be a vector containing the angular positions of integration
%(ode15s is a variable integration step solver).
%pseal will be a matrix: each row is a vector containing the integrated pressure
%in the corresponding sealed volume for a whole angular pitch (e.g. row 1 
%contains the pressure inside sealed volume 1, row 2 sealed volume 2, and
%so on...).
%tseal and pseal are updated for each repetition of the 'for cycle jj'.
%Internal 'if' splits between scheme N->N+1->N and N->N-1->N.
for jj = 1 : repeat 
    if nsv(1) ~= max(nsv)    % if the scheme is N -> N+1 -> N  (ex. 6 -> 7 -> 6)
        for ii = 1 : length(index)
            if ii == 1 && jj == 1  %very first integration step                  
               thseal = [TETA(1), TETA(index(1))]; %ends of integration 
                                                   %(until transition N -> N+1)
               CI =  pa * ones(nsv(ii), 1);  %initial condition setting: pressure in all
                                             %N volumes is set to pa
               %ODE set solving (see 'odesealfcn.m' inputs)
               [ts, ps] = ode15s(@(t, p) odesealfcn(t,p, nsv(ii), index), [thseal(1), thseal(2)], CI);
               tseal = ts'; 
               pseal = [pa*ones(1, size(ps', 2)); ps']; %add first line for volume #N+1
                                                        
            elseif ii == 1 && jj > 1 %first integration step (repeated)
               thseal = [TETA(1), TETA(index(1))]; %ends of integration 
                                                   %(until transition N -> N+1)
               CI = pseal(1:end-1, end); %initial condition setting: pressure in all
                                         %N volumes is set to end values of
                                         %previous complete integration
               %ODE set solving (see 'odesealfcn.m' inputs)
               [ts, ps] = ode15s(@(t, p) odesealfcn(t,p, nsv(ii), index), [thseal(1), thseal(2)], CI);
               tseal = ts';  
               pseal = [pa*ones(1, size(ps', 2)); ps']; %add first line for volume #N+1
            elseif ii == 2 %second integration step
               thseal = [TETA(index(ii-1)), TETA(index(ii))] ; %ends of integration 
                                                               %(from transition N -> N+1 to transition N -> N-1)
               CI = pseal(:, end); %initial condition setting: pressure in all
                                   %N+1 volumes is set to end values of
                                   %first integration step
               [ts, ps] = ode15s(@(t, p) odesealfcn(t,p, nsv(ii), index), [thseal(1), thseal(2)], CI);
               tseal = [tseal, ts']; pseal = [pseal, ps']; %results from integration steps 1 and 2 are seamed together
            elseif ii == 3  %third integration step
               thseal = [TETA(index(ii-1)), TETA(index(ii))] ; %ends of integration 
                                                               %(from transition N -> N-1 to end of angular pitch)
               CI = pseal(1:end-1, end); %initial condition setting: pressure in all
                                         %N volumes is set to first N end values of
                                         %second integration step
               [ts, ps] = ode15s(@(t, p) odesealfcn(t,p, nsv(ii), index), [thseal(1), thseal(2)], CI);
               PS = [ps'; pm*ones(1, size(ps', 2)) ];
               tseal = [tseal, ts']; pseal = [pseal, PS]; %results from integration steps 1, 2 and 3 are seamed together
            end
        end

    else                 %if the scheme is N -> N-1 -> N   (ex. 7 -> 6 -> 7 )
                         %Procedure is similar to previous, but initial
                         %conditionss vectors are set in a slightly
                         %different way to match ODE system dimensions
       for ii = 1 : length(index)
            if ii == 1 && jj == 1  %very first integration step
               thseal = [TETA(1), TETA(index(1))];
               CI =  pa * ones(nsv(ii), 1);  
               [ts, ps] = ode15s(@(t, p) odesealfcn(t,p, nsv(ii), index), [thseal(1), thseal(2)], CI);
               tseal = ts'; 
               pseal = ps';
            elseif ii == 1 && jj > 1 %first integration step (repeated)
               thseal = [TETA(1), TETA(index(1))];
               CI = pseal(:, end);
               [ts, ps] = ode15s(@(t, p) odesealfcn(t,p, nsv(ii), index), [thseal(1), thseal(2)], CI);
               tseal = ts';  
               pseal = ps';  
            elseif ii == 2  %second integration step
               thseal = [TETA(index(ii-1)), TETA(index(ii))] ; 
               CI = pseal(1:end-1, end);
               [ts, ps] = ode15s(@(t, p) odesealfcn(t,p, nsv(ii), index), [thseal(1), thseal(2)], CI);
               tseal = [tseal, ts']; 
               PS = [ps'; pm*ones(1, size(ps', 2))];
               pseal = [pseal, PS];
            elseif ii == 3  %third integration step
               thseal = [TETA(index(ii-1)), TETA(index(ii))] ;
               CI = [pa; pseal(1:end-1, end)];
               [ts, ps] = ode15s(@(t, p) odesealfcn(t,p, nsv(ii), index), [thseal(1), thseal(2)], CI);
               tseal = [tseal, ts']; 
               pseal = [pseal, ps'];
            end
        end  
    end
%De-comment following lines to plot provisional Sealed Volumes pressure evolution
%along with SEP process run (not advised, better to refer to final pressure
%history plot when 'while cycle' iterations are ended)
%     figure(200 + iter)
%     for kk = 1 : size(pseal, 1)
%         plot(tseal + (jj-1)*360/Z, pseal(kk, :), 'LineWidth', 1.5);
%         hold on
%     end
%     plot(TETA(index(1))*[1, 1] + (jj-1)*360/Z, [0, 1]*pm, '--k')
%     plot(TETA(index(2))*[1, 1] + (jj-1)*360/Z, [0, 1]*pm, '--k')
%     plot(TETA(index(3))*[1, 1] + (jj-1)*360/Z, [0, 1]*pm, '--k')
%     hold on
end

hold off

%SAMPLING OF RESULTS
%following lines are for sampling pseal pressure matrix and ptrap vector
%at fixed angular steps. Angular step must be the same used for geometrical
%sampling in 'geometry.m'.
SampleDeg = TETA(2) - TETA(1); %keep same as geometry module sampling angular step [deg]
t_trap_samp = 0:SampleDeg:th(end);
p_trap_samp = makima(tt, pt, t_trap_samp); %Trapped volume pressure in sampling
                                           %angular positions
t_seal_samp = 0:SampleDeg:tseal(end);
p_seal_samp = zeros(size(pseal, 1), length(t_seal_samp));
%the following lines are to delete eventual duplicate values at ode15s
%integration boundaries ends
for ii = 2 : length(tseal) - 2  
    if tseal(ii) == tseal(ii - 1),  tseal(ii) = []; pseal(:, ii) = [];  end
end

for ii = 1 : size(pseal, 1)
 p_seal_samp(ii, :) = makima(tseal, pseal(ii, :), t_seal_samp); %Sealed volumes pressure in
                                                                %sampling angular positions
end

%PRESSURE MEAN FORCE + CONTACT MEAN FORCE
%following lines are for evaluating resultant pressure force + contact
%force. Resultant force is the input of the hydrodynamic bearing function
%'JBE.m'.

%following lines are for building field vectors around driver gear. 
%thsealvec will contain the angular portions associated with every sealed
%volume (counter clockwise order)
%psealvec will contain corresponding pressure values (counter clockwise order)
%Evaluation is split depending on the scheme (N -> N+1 -> N or N -> N-1-> N).
for ii = 1 : length(TETA) 
 thsealvec = (thseal2(ii)-thseal1(ii))/NoSV(ii) * ones(1, NoSV(ii)); 
 if nsv(1) ~= max(nsv)    %if scheme is N -> N+1 -> N  (ex. 6 -> 7 -> 6 )
         if TETA(ii) < TETA(index(1)) %(N contemporary sealed volumes)
                 psealvec = p_seal_samp(end:-1:2, ii);
         elseif TETA(ii) >= TETA(index(1)) && TETA(ii) < TETA(index(2))
                                      %(N+1 contemporary sealed volumes)
                 psealvec = p_seal_samp(end:-1:1, ii);
         else                         %(N contemporary sealed volumes)
                 psealvec = p_seal_samp(end-1:-1:1, ii);
         end
  
 else                    %if scheme is N -> N-1 -> N  (ex. 6 -> 5 -> 6 )
         if TETA(ii) < TETA(index(1))
                 psealvec = p_seal_samp(end:-1:1, ii);
         elseif TETA(ii) >= TETA(index(1)) && TETA(ii) < TETA(index(2))
                 psealvec = p_seal_samp(end-1:-1:1, ii); 
         elseif ii == index(2)
                 psealvec = [p_seal_samp(end-1:-1:1, ii); 0]; 
         else
                 psealvec = p_seal_samp(end:-1:1, ii);
         end
 end  

 %Build complete field vectors, including every distinct control volume
 %(counterclockwise order) : delivery, sealed volumes, suction, trapped
 %volume (when existing).
 %Angular portions start from teta=0 in XY reference system and go
 %counterclockwise.
if thcontact2(ii) < 10000  && ii <= length(p_trap_samp) %trapped volume exists
    thfield = pi/180 * [thcontact1(ii), thseal1(ii)-thcontact1(ii),  thsealvec, ...
                        thcontact2(ii)-thseal2(ii), 360-thcontact2(ii)]; 
                            %vector contains all angular portions
    pfield = [p_trap_samp(ii), pm, psealvec', pa, p_trap_samp(ii)];
                           %vector contains all corresponding pressure values
else 
   if thcontact1(ii) < 90  %only one contact, on delivery side
     thfield = pi/180 * [thcontact1(ii), thseal1(ii)-thcontact1(ii),  thsealvec, ...
                           360-thseal2(ii)];  %vector contains all angular portions
     pfield = [pa, pm, psealvec', pa];  %vector contains all corresponding pressure values

   else                   %only one contact, on suction side
     thfield = pi/180 * [thseal1(ii), thsealvec, ...
                           thcontact1(ii)-thseal2(ii), 360-thcontact1(ii)]; 
                             %vector contains all angular portions
     pfield = [pm, psealvec', pa, pm];
                             %vector contains all corresponding pressure values
   end
end
   %PRESSURE FORCE on driver gear evaluation (see function 'pressureForce.m' )
   [fxp(ii), fyp(ii)] = pressureForce(pfield, thfield, w, TR(end));

   %PRESSURE TORQUE [Nm] on driven gear evaluation (see function 'meshingTorque.m')
   if thcontact2(ii) < 10000 %if there are 2 contact points and trapped volume exists
   Torque(ii) = meshingTorque(p_trap_samp(ii), r2contact1(ii), r2contact2(ii), TR(end), w);
   else %if ther is only 1 contact point and trapped volume does not exist
       Torque(ii) = meshingTorque(NaN, r2contact1(ii), r2contact2(ii), TR(end), w); 
   end

end

fp = [mean(fxp), mean(fyp)]; %Mean pressure force  in XY components vector [N]
fmesh= mean(Torque)/(0.001*Rbase); %Mean contact force magnitude [N]
                                   %NOTE: base circle radius = d_prim/2 *cos(alpha_w)
fmg = fmesh*[-sind(alpha_w), cosd(alpha_w)];   %Mean contact force in XY components vector [N]
fgear = fp + fmg; %Total force on gear [N]
load = .5*norm(fgear, 2);  %load magnitude on bearing [N]
                           %(load is divided between two JBs)
load_angle = rad2deg(mod(atan2(fgear(2), fgear(1)), 2*pi)); %load direction on bearing [deg]

%Current bearing position is stored to be later compared with new bearing
%position (i.e bearing function ouput with updated load)
e_old = sqrt(BEAR(iter, 1)^2 + BEAR(iter, 2)^2 );
G_old =  rad2deg(mod(atan2(BEAR(iter, 2), BEAR(iter, 1)), 2*pi));

%Bearing function
%INPUT: bearing geometric dimensions [m] (not [mm]!!!), speed [rev/sec], load, dynamic viscosity
%OUTPU: eccentricity and attitude angle (direction w/ respect to load direction)
[ecc,attitude_angle] = JBe(BearR*0.001, Lbear*0.001, Clear*0.001, RPM/60, load, mu);
ecc = ecc * 1000; %[mm]
GAMMA = load_angle - attitude_angle; %direction in XY reference system
xbear = ecc*cosd(GAMMA);
ybear = ecc*sind(GAMMA);

%Write updated provisional SEP in Excel FILE
BEAR = [BEAR; xbear, ybear, ecc, GAMMA];
writematrix(BEAR, filename, 'Sheet', 1, 'Range', 'A2')

%evaluate relative distance between current and previous SEP guess
rel_err_e = abs(ecc - e_old)/abs(e_old);
rel_err_G = abs(GAMMA - G_old)/abs(G_old);

iter = iter + 1;

end

%SEP iterative process is now ended !!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% VANE PRESSURE HISTORY PLOT%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Vane pressure history plot is built by seaming together sealed volumes
%pressures, delivery pressure, trapped volume pressure, suction pressure
figure(300)
for ii = 1 : size(p_seal_samp, 1)
    if nsv(1) ~= max(nsv)
        plot(t_seal_samp + (ii-1)*360/Z - TETA(index(1)), p_seal_samp(ii, :),':k', 'linewidth', 1.5)
        xlim([0, 360])
        hold on
    else 
       vv1 = find(t_seal_samp > TETA(index(2)));
       if ii > 1
           p_seal_samp1 = [p_seal_samp(ii, vv1), p_seal_samp(ii, 1:vv1(1)-2)];
           t_seal_samp1 =  t_seal_samp(t_seal_samp ~= TETA(index(2)));
       else
           p_seal_samp1 = [0, p_seal_samp(ii, vv1), p_seal_samp(ii, 1:vv1(1)-2)]; 
           t_seal_samp1 = t_seal_samp;
       end
        plot(t_seal_samp1 + (ii-1)*360/Z, p_seal_samp1,':k', 'linewidth', 1.5)
        xlim([0, 360])
        hold on
    end

    if ii == size(p_seal_samp, 1) 
        tman = t_seal_samp(end)+(ii-1)*360/Z:SampleDeg:(max(thseal2)-thcontact1(1) - 10);
        tasp = t_trap_samp(end)+tman(end):SampleDeg:360;
        plot(tman, pm*ones(length(tman), 1) , ':k', 'linewidth', 1.5)
        plot(t_trap_samp + tman(end), p_trap_samp, ':r', 'linewidth', 1.5)
        plot(tasp, [p_trap_samp(end); pa*ones(length(tasp)-1, 1)], ':k', 'linewidth', 1.5)
    end
end

%Plot pressure force on driver gear over a complete angular pitch
figure(400)
plot(TETA, fxp, '-r', TETA, fyp, '-b')
ylabel('Force [N]')
xlabel('theta [deg]')
legend('fxp', 'fyp')
%Plot pressure torque on driven gear over a complete angular pitch
figure(500)
plot(TETA, Torque)
ylabel('Torque [Nm]')
xlabel('theta [deg]')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% Volumetric efficiency evaluation %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The following lines are for evaluating volumetric efficiency

Cil = pi/4*w*((2*TR(end))^2 - (2*TR(1))^2); %theoretic pump displacement [mm^3/rev]

%volumetric efficiency evaluation (see function 'etaVol.m'). Sampled
%pressure values for both sealed volumes and trapped volume are used here.
%NOTE: if the pump does not have backend grooves Q_back is always a null vector).
%NOTE : fluid density is set at delivery pressure.
[nVol, Q_tip, Q_lat, Q_drain, Q_mesh, Q_relief, Q_back, Q_leak] = etaVol(t_seal_samp, p_seal_samp, p_trap_samp, Cil, density(pm,BM,pcav,rhocav));

%Plot mean leakages (tooth tip, tooth side, delivery relief groove, backend
%groove, meshing area, drainage) in a pie chart [liter/sec].
Qpie = [mean(Q_tip)/mean(Q_leak), mean(Q_lat)/mean(Q_leak),...
       mean(Q_relief)/mean(Q_leak), mean(Q_back)/mean(Q_leak),...
       mean(Q_drain+Q_mesh)/mean(Q_leak), ];
figure(600)
labels = {'Tooth tip','Tooth sides', 'Relief groove', 'Backend groove', 'Meshing area, Drainage'};
pie(Qpie,labels)
title(['Q_{loss} = ', num2str(mean(Q_leak)*1e-6,'%.4f'), ' [l/s]'])

toc
beep  %sound notification when script run ends

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% ! After the run is done, please delete the Excel FILE %%%%%%%%%%%%%
%%%%%%% 'BearingCenterCoord.xlsx before doing a new run !!!!! %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

