function [ transition  ] = gillespie(rmatrix,input)
%gillespie step: calculates a new state and a step time given 
%a rate matrix and an initial state as input
    
%calculate total rate
tot_states = size(rmatrix,1);
rate =0;
for i=1:tot_states
    rate = rate + rmatrix(input,i);
end
rate

%Calculate the time until next reaction
deltat = exprnd(1/rate);

%Choose reaction
reactvar = rand*rate;
temp=0;
for i=1:tot_states
    temp = temp+ rmatrix(input,i);
    if temp > reactvar
        reactid=i;
        break
    end
end

%Return new state and time taken
transition = [reactid,deltat];

