function [Uage,Fage] = ...
            HPVprogression(U0,age,foi,p0,p1,ps,h0,h1,w,sage,satt,sens)
% function  [Uage,Fage,age] =
%            HPVprogression(U0,age,foi,p0,p1,ps,h0,h1,w,sage,satt,sens)
% HPV progression model, combines natural history with screening.
% The states in U are ordered as X = (S,R,X0,X1,X2,X3,XA,DCA)
% 
% INPUT:   size    
%     U0   (xm,1)  Initial states for state pdf U. 
%     age  (1,am)  Ages for the population (age(j+1)-age(j) will be dt!)
%     foi (1,am)   foi(j) is the force of infection at age age(j)
%     p0  (1,1)    progression rate from X0
%     p1  (1,1)    progression rate from X1-X3
%     ps  (1,1)    symptom rate from XA to DCA (detected cancer)
%     h0  (1,1)    clearance rate from X0 to R
%     h1  (1,1)    clearance rate from X1-X3 to R
%     w   (1,1)    waning natural immunity rate from R to S
%     sage (1,sm)  screening ages
%     satt (1,sm)  satt(j1)=screening attendance probability at sage(j)
%     sens (xm,1)  sens(j) = sensitivity for x(j)
% OUTPUT: size     
%     Uage (xm,am) Uage(:,j) is pdf for population at age(j)
%     Fage (xm,sm) Fage(:,j) are the findings at screening age sage(j)

% Some sizes
xm = length(U0);    % Number of X states
am = length(age);  % Number of age ticks
sm = length(sage); % Number of screening rounds
dt = age(2)-age(1); % Define dt here to avoid contradictions

Uage = zeros(xm,am);
Fage = zeros(xm,sm);
snow = 1; % count for screening rounds (next one to perform)

Unow = U0;
Uage(:,1) = Unow;

for j1 = 1:(am-1) % update time am-1 times to get to age(am)
   
  % If screening round, perform screening:
  if snow<=sm   % Screening rounds still to come
    if age(j1)>=sage(snow)  % We achieved the next screening age
       [Unow,Fnow] = PerfScr(Unow,satt(snow),sens); % Perform screening
       Fage(:,snow) = Fnow;
       snow = snow+1;
    end
  end
  
  % Natural history step to the next time/age tick:
  Unow = NatHistStep(Unow,dt,foi(j1),p0,p1,ps,h0,h1,w);

  % Save U after step:
  Uage(:,j1+1) = Unow;
    
end
