% Population dynamics of a group of neanderthals

%close all; 
clear all;

%% maximum age:
Nages=80;
% number of years to run the model:
Tend=5000;

%% birth rage per woman per year in fertile age group:
births_per_year=1.0/3.0;

%% fertility age range:
fertility_s=19;
fertility_e=28;

% different possible scenarios of mortality rate, see below:
scenario=2;

%% this is used to modify the survivorship curve taken from a book:
survivorship_factor=1.2;

% initizlize population distribution variable:
P=zeros(Nages,Tend);
% initial conditions:
load Output/restart.dat;
P(:,1)=restart*0.25;

%% frequency of plotting:
Tplot=50;


%% some useful counters:
Pmax=max(P(:,1));
Ptotal=zeros(Tend,1);
Ptotal(1)=floor(sum(P(:,1)));

% time stepping:
for t=2:Tend

  %% plot current age distribution:
  
  Ptotal(t)=floor(sum(P(:,t-1)));
  if t==2 || floor((t-1)/Tplot)*Tplot==t-1
    figure(1);
    clf;
    bar(P(:,t-1));
    Pmax_plot=max(Pmax,max(P(:,t-1)));
    axis([1 Nages 0 Pmax_plot])
    h1=xlabel('age');
    h2=ylabel('population');
    h3=title(sprintf('population distribution at year=%g, total=%g',t-1,Ptotal(t)));
    set([gca,h1,h2,h3],'fontsize',16);
    pause
  end

  %% multiply by 0.5 to get the number of females:
  P(1,t)=floor(births_per_year*0.5*sum(P(fertility_s:fertility_e,t-1)));

  %% mortality:

  %% survivorship:
  
  if scenario==1 || scenario==2
    %% from some book that ofer brought:
    survivorship=[2  5 .61;
                  6  10 .58;
                  11 15 .55;
                  16 20 .53;
                  21 25 .50;
                  26 30 .46;
                  31 35 .42;
                  36 40 .39;
                  41 45 .35;
                  46 50 .31;
                  51 55 .27;
                  56 60 .23;
                  61 65 .17;
                  66 70 .12;
                  71 75 .07;
                  76 80 .03];
  end
  
  if scenario==2
%     figure
%     plot(survivorship(:,3));

    survivorship(:,3)=survivorship(:,3) ...
    *survivorship_factor .* (1-tanh((survivorship(:,1)-35)/7))/2.0;
%     figure
%     plot(survivorship(:,3));
%     pause
  
  end
  if max(survivorship(:,3))>0.99
    fprintf(1,'*** Error: survivorship>1!\n')
    survivorship(:,3)
  end
  
  num_age_groups=length(survivorship(:,1));
  for age_group=1:num_age_groups
  
    age_s=survivorship(age_group,1);
    age_e=survivorship(age_group,2);
    if age_group==1
      survival_rate=survivorship(age_group,3);
    else
      survival_rate=survivorship(age_group,3)/survivorship(age_group-1,3);
    end
    survival_rate_per_year=survival_rate^(1/(age_e-age_s+1));
    for age=age_s:age_e
      P(age,t)=P(age-1,t-1)*survival_rate_per_year;
    end

  end
  
end

figure(2)
clf
% plot survival rate as function of age:
subplot(2,1,1)
survival=diag(P/P(1,1));
h=plot(survival);
set(h,'linewidth',2);
h1=xlabel('Age');
h2=ylabel('Fraction');
h3=title('Surviving fraction');
set([gca,h1,h2,h3],'fontsize',16);

% plot total population as function of time:
subplot(2,1,2)
h=plot(Ptotal);
set(h,'linewidth',2);
h1=xlabel('time (years)');
h2=ylabel('population');
h3=title('Total population');
set([gca,h1,h2,h3],'fontsize',16);

filename=sprintf(...
  'Output/Neanderthal_birth_every_%gyrs_Fert-range=%g-%g_surv_fact=%g.jpg' ...
  ,1/births_per_year,fertility_s,fertility_e,survivorship_factor);
saveas(gcf,filename)

%% to save current final state for restarting, uncomment these:
fid=fopen('Output/restart-data.dat','w');
fprintf(fid,'%g\n',P(:,t));
fclose all;