%%%%valori iniziali delle coperture
%copertura iniziale di prato semplice(tra 0 e 1)
global init_l = 0.148;
%copertura iniziale di prato incolto(tra 0 e 1)
global init_f = 0.147;
%copertura iniziale di prato complesso(tra 0 e 1)
global init_m = 0.303;
%copertura iniziale di arbusti(tra 0 e 1)
global init_s = 0.123;
%copertura iniziale di alberi(tra 0 e 1)
global init_t = 0.402;
%copertura iniziale di sottobosco(tra 0 e 1)
global init_u = (1 - init_l - init_f - init_m);

%%%%variabili di controllo
%altitudine (m)
global A = 600;
%Carico Animale Globale (UBA giorni/ha/anno)
global  GSD =0;

%%%%parametri
%Valore Pastorale medio incolto
global PVF = 10;
%Valore Pastorale medio pascolo semplice
global PVL = 20;
%Valore Pastorale medio pascolo complesso
global PVM = 40;
%Valore Pastorale medio sottobosco
global PVU = 5;
%Valore Pastorale globale
global GPV = 21;

%%%%equazioni
%valore pastorale locale
global LPV = PVF*init_f+PVL*init_l+PVM*init_m+PVU*init_u;
%carico animale locale(UBA giorno/ha/anno)
global init_lsd = GSD*LPV/GPV ;
%parametri
global d = -0.013305; 
global b = -0.005;
global c = 1081.093;
%effetti dell'altitudine(tra 0 e 1)
global AE = exp(b*(A-c))/(exp(b*(A-c))+1);
%capacità portante globale(UBA giorno/ha/anno)
global GCC = 115000*GPV/18/A;
%tasso di utilizzazione globale(tra 0 e 1)
global GU =GSD/GCC;


%%%%%%%% parametri temporali di simulazione
%t inizio simulazione
%global t_initial = 0;
%t fine simulazione
%global t_final = 300;
	
%global monitor_points = 200;
%global delta_t = (t_final - t_initial) / monitor_points;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = fun (t,x,xd)

global d
global GSD 
global PVF
global PVL
global PVM
global PVU
global GPV
global b 
global A 
global c 
global AE 
global SCr
global SCs
global SCp
global TCr
global TCs
global TCp
global GCC
global GU 
global init_t
global t_initial 
global t_final 
global monitor_points 
global delta_t 

%x1
global init_l   
%x2
global init_f 
%x3
global init_m 
%x4
global init_u 
%x5
global LPV 
%x6
global init_lsd 
%x7
global init_s 

%%%%%% L(t) = L(t−dt)+(−LtoM−LtoF)×dt Prato Semplice
%%%%%% L(t)= -((0.06*L*GP)-(0,05*M*(1-GP))) - ((0,5*L*(1-GP)*(O.3+S)*AE)-(0.1*F*GP)) 
f(1) = -((0.06*x(1)*(1-exp(d*x(5))))-(0.05*x(3)*(exp(d*x(5)))))-( (0.5*x(1)*(exp(d*x(5)))*(0.3+x(6)))*AE-(0.1*x(2)*(1-exp(d*x(5)))));

%%%%%% F(t) = F(t−dt)+(LtoF+MtoF−FtoU)×dt  Incolto
%%%%%% F(t)= ((0,5*L*(1-GP)*(O.3+S)*AE)-(0.1*F*GP)) + (O.8*M*(1-GP)*(O.3+S)-(0.4*F*(GP)^2) - (F*(T-U)) 
f(2) =(((0.5*x(1)*(exp(d*x(5)))*(0.3+x(6)))*AE-(0.1*x(2)*(1-exp(d*x(5))))) + (0.8*x(3)*(exp(d*x(5)))*(0.3+x(6))-(0.4*x(2)*((1-exp(d*x(5))^2)))) - (x(2)*(x(6)-x(4))));

%%%%%% M(t) = M(t−dt)+(LtoM−MtoF)×dt Prato complesso
%%%%%% M(t)=((0.06*L*GP)-(0,05*M*(1-GP))) - (O.8*M*(1-GP)*(O.3+S)-(0.4*F*(GP)^2)
f(3) = ((0.06*x(1)*(1-exp(d*x(5))))-(0.05*x(3)*(+exp(d*x(5))))) - (0.8*x(3)*(exp(d*x(5)))*(0.3+x(6))-(0.4*x(2)*((1-exp(d*x(5))^2))));

%%%%%% U(t) = U(t−dt)+(FtoU)×dt
%%%%%% U(t) = (F*(T-U))
f(4) = (x(2)*(x(6)-x(4)));

%%%%% if per calcolare LSD
LCC =LPV*115000/18/A;
LU = x(5)/LCC;
if (LU < 1)
  dLSD = x(5) * (GU-LU);
else
  dLSD = x(5) *(1-LU);
endif

%%%%%% LSD(t) =LSD(t−dt)+(dLSD)×dt capacità di carico locale
f(5) = dLSD ;



%%%%% S(t)=S(t−dt)+(dSi−dSo)×dt        Arbusti

dSo = (0.1 + 1 - exp(d*x(5))) * (0.1 + x(7)) * x(6) * 0.2;
%
%%%%% funzione DELAY sulla variabile S(t)
%%%%% dSi = 0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE -----definizione di Stella
%
dSi = 0.2*((xd(2,1)+0.5*xd(4,1))*(exp(d*xd(5,1))))*(1-x(6))*(1-x(7))*AE;


f(6) = dSi - dSo;



%%%%% T(t) = T(t−dt)+(dTi−dTo)×dt   Alberi

dTo = x(7)*0.005;
%
%%%%% funzione DELAY sulla variabile T(t)
%%%%% dTi =DELAY(S×(1−GP),10)×(1−T)×0.1×AE  -----definizione di Stella
%
dTi = ((xd(6,2)*(exp(d*xd(5,2))))) * (1-x(7)) * 0.1 * AE;


f(7) = dTi - dTo ;

endfunction


t=[0;300];


res = ode45d (@fun, t, [init_l;init_f;init_m;init_u;init_lsd;init_s;init_t], [5,10],ones(7,2));


plot (res.x, res.y);
xlabel('Tempo (anni)')
ylabel('Variabili')
title('Patumod')
legend ("L", "F", "M", "U", "LSD", "S", "T")
axis([0,200, 0, 1])
replot
print('figuregsd0a600.ps', '-deps');