function D = fsim ( Npops, mu0, mu1, gamma0, gamma1, IC, NReps, MaxTime )
%function to perform metapopulation simulations

%Graham Medley, February 2013; University of Warwick

%data structure for the simulation output. Initially this is a 3-D array
%consisting of output at the start of each time interval (and the end of
%the last one) for each repetition and each population
D = zeros ( NReps, MaxTime+1, Npops );

for irep = 1 : NReps
    
    pops = IC; %pops is set to the initial condition
    time = 0; %time is set to zero
    outcount = 1; %this is a counter to record the time interval
    D ( irep, outcount, :  ) = pops; %record the initial condition as the first output
    outtime = 1; %this is the time of the next output
    
    %perform a translocation ( at the start )
    pops = translocate ( time, pops );
    
    %continue until MaxTime or metapopulation extinction or one population
    %gets silly
    while ( time < MaxTime ) && any ( pops > 0 ) && all ( pops < 100 )
        
        BirthRates = max ( 0, gamma0 -  gamma1 .* pops );
        DeathRates = mu0 + mu1 .* pops;
        BR = pops .* BirthRates; %calculate the population level birth rates
        DR = pops .* DeathRates; %calculate the population level death rates
        
        EventRates = cumsum ( [ sum(BR) sum(DR) ] );
        TotalRate = EventRates ( end ); %the total rate of events
        
        randnums = rand ( 1, 3 ); %get 3 random numbers from a uniform distribution btn. 0 and 1
        tstep = - log( randnums(1) ) / TotalRate; %calculate the time to the next event using the 1st random number
        time = time + tstep; %time moves ever on...
        oldpops = pops; %remember the population in its current state
        
        %this section chooses the event and the population, and changes the
        %population accordingly
        %use the second random number to decide btn. events
        eventindex = find ( EventRates > randnums ( 2 ) * TotalRate );
        switch eventindex ( 1 )
            case 1
                PopRate = cumsum ( BR );
                index = find ( PopRate > ( randnums ( 3 )*PopRate ( end ) ) );
                pops ( index ( 1 ) ) = pops ( index ( 1 ) ) + 1; %add one to the population
            case 2
                PopRate = cumsum ( DR );
                index = find ( PopRate > ( randnums ( 3 )*PopRate ( end ) ) );
                pops ( index ( 1 ) ) = pops ( index ( 1 ) ) - 1; %remove one from the population
            otherwise
                %this indicates an error - report it and stop simulating
                disp ( 'error in fsim' );
                return
        end
        
        while ( time > outtime ) && ( outtime <= MaxTime )
            outcount = outcount + 1; %record the next time point
            D ( irep, outcount, : ) = oldpops; %population as it is at outtime
            outtime = outtime + 1; %next time point for recording
            pops = translocate ( time, pops ); %do the translocation annually
        end
        
    end %while
    
end %for NReps