%dosims - 
%this script conducts a simple metapopulation simulation.
%The parameters are assigned and simple graphs drawn from 
%the data structure (i.e. the output from the simulation).
%The simulations are performed by the function "fsim"

%Graham Medley, February 2013; University of Warwick

%scrub all previous work
clear
fprintf ( '\nRUNNING METAPOPULATION SIMULATION...\n' );

%Size of the metapopulation
Npops = 10;
fprintf ( 'Number of populations: %i\n', Npops );

%Birth and death rates - these are equal for each population
mu0 = ones ( 1, Npops ) * 0.5; %2 years life expectancy optimum
mu1 = ones ( 1, Npops ) * 0.05; %1 year life expectancy at N=10
gamma0 = ones ( 1, Npops ) * 1.5 / 2; %birth rate optimum
gamma1 = ones ( 1, Npops ) * 0; %density independent birth rate
K = 10; %carrying capacity of each patch

%Initial population sizes
IC = [ 0 0 0 0 0 10 10 10 10 10 ];
fprintf ( 'Initial sizes: ' );
fprintf ( '%i ', IC )
fprintf ( '\n' )

%Maximum time over which the simulation is to be performed
MaxTime = 50;

%Number of repetitions - if you change this then the multiplot graphs below will
%have to be changed.
NReps = 20;
fprintf ( 'Maximum time: %i, and %i repetitions\n', MaxTime, NReps );

D = fsim ( Npops, mu0, mu1, gamma0, gamma1, IC, NReps, MaxTime );

%D is a multidimensional array, with dimensions of (repetitions, time, populations)
%the functions below calculate statistics relating to the simulations in D
%and produce some graphs

plotSummary ( D );

plotRepetitions ( D );

[ TE ] = calcExtinctions ( D ); %calculate the mean times to extinction; TE==0 is not extinct

fprintf ( 'The proportion of extinct metapops was %6.3f\n', sum ( TE>0 ) / NReps );
fprintf ( 'with times to extinction\n' );
fprintf ( '%5.2f ', TE(TE>0)' ); %times to extinction
fprintf ( '\n' );