clear all

N0=2;
h=0.1;

%kameno doba
t1_min=0;
t1_max=100;
K1=5; 

% poljoprivredna revolucija
t2_min=t1_max; % 10.000
t2_max=t1_max+9000;
K2=450;

% kuga
t3_min=t2_max; %1300
t3_max=t2_max+100;
K3=K2*0.7;

t4_min=t3_max; % 1400.
t4_max=t3_max+300;
K4=600;

% industrijska revolucija 1700.
t5_min=t4_max;
t5_max=t4_max+250;
K5=3000;

% zelena revolucija (pesticidi, genetske modifikacije...) 1950.
t6_min=t5_max;
t6_max=t5_max+70;
K6=7000;

%N1(1)=N0;
i=1;
j=1;
k=1;
l=1;
s=1;
q=1;
for t=t1_min:h:t6_max
    if t<=t1_max   % Kameno doba
        %N0=N(end);
        r=0.04;
        N1(i)=K1/(1+(K1/N0-1)*exp(-r*t)); 
        i=i+1;               
    else if t<=t2_max   % 8.000pne -1300
            N0=N1(end);
            r=0.04;
            N2(j)=K2/(1+(K2/N0-1)*exp(-r*(t-t2_min))); % 
            j=j+1;    
        else if t<=t3_max  % 1300-1400.
                r=0.04;
                N0=N2(end);
                N3(k)=K3/(1+(K3/N0-1)*exp(-r*(t-t3_min))); % 
                k=k+1;    
            else if t<=t4_max % 1400-1700.
                    N0=N3(end);
                    r=0.04;
                    N4(l)=K4/(1+(K4/N0-1)*exp(-r*(t-t4_min))); % 
                    l=l+1;   
                else if t<=t5_max % 1700-1950.                        
                        N0=N4(end);
                        r=0.04;
                        N5(s)=K5/(1+(K5/N0-1)*exp(-r*(t-t4_min))); % 
                        s=s+1;   
                    else if t<=t6_max %1950-2020
                            N0=N5(end);
                            r=0.04;
                            N6(q)=K6/(1+(K6/N0-1)*exp(-r*(t-t4_min))); % 
                            q=q+1;  
                        end
                    end
                end
            end
        end
    end
end


%hold on
plot(t1_min:h:t1_max,N1)
legend('Kameno doba')
pause
plot(t2_min+h:h:t2_max,N2)
legend('cc 8.000pne -1300.')
pause
plot(t3_min+h:h:t3_max,N3)
legend('1300-1400.')
pause
plot(t4_min+h:h:t4_max,N4)
legend('1400-1700.');
pause
plot(t5_min+h:h:t5_max,N5)
legend('1700-1950.')
pause
plot(t6_min+h:h:t6_max,N6)
legend('1950-2020.')
pause
plot(t1_min:h:t6_max, [N1,N2,N3,N4,N5,N6])
axis equal

