%%%-------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;


%%%%-------relazioni
%--valore pastorale locale
global init_LPV = 1*PVF*init_F+PVL*init_L+PVM*init_M+1*PVU*init_U;
%--carico animale locale (UBA giorno/ha/anno)
global init_LSD = GSD*init_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;
%--capacità portante locale (UBA giorno/ha/anno)
global init_LCC = init_LPV*115000/18/A;
%--tasso di utilizzazione locale (tra 0 e 1)
global init_LU = init_LSD/init_LCC

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

global GSD 
global A 
global d
global b 
global c 
global PVF
global PVL
global PVM
global PVU
global GPV
global AE 
%global SCr
%global SCs
%global SCp
%global TCr
%global TCs
%global TCp
global GCC
global GU 
global LPV 
global LCC
global LU
global init_LU

%x1
global init_L
%x2
global init_F 
%x3
global init_M 
%x4
global init_U 
%x5
global init_LSD
%x6
global init_S 
%x7
global init_T





%%%%%% 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(7)-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(7)-x(4)));



%%%%% if per calcolare LSD


if (t=0)
	if (init_LU < 1)
	  dLSD = x(5) * (GU-LU);
	else
	  dLSD = x(5) *(1-LU);
	endif
else 
	if (LU < 1)
	  dLSD = x(5) * (GU-LU);
	else
	  dLSD = x(5) *(1-LU);
	endif 
endif

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


%%%%% S(t)=S(t−dt)+(dSi−dSo)×dt        Arbusti
%%%%% S(t) = (0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE)-((0.1+GP)*(0.1+T)*S*0.2)---definizione del delay secondo Stella
f(6) =(0.2*((xd(2,1)+0.5*xd(4,1))*(exp(d*xd(5,1))))*(1-x(6))*(1-x(7))*AE) - ((0.1 +(1-exp(d*x(5))))*(0.1+x(7))*x(6)* 0.2);


%%%%% T(t) = T(t−dt)+(dTi−dTo)×dt   Alberi
%%%%% T(t) = (DELAY(S×(1−GP),10)×(1−T)×0.1×AE) - T*0.005 ---definizione del delay secondo Stella
f(7) = ((xd(6,2)*(exp(d*xd(5,2))))*(1-x(7))*0.1*AE) - (x(7)*0.005);


LPV = 1*PVF*x(2)+PVL*x(1)+PVM*x(3)+1*PVU*x(4)

LCC = LPV*115000/18/A

LU = x(5)/LCC


endfunction


t=[0:200];


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,100, 0, 1])
replot
print('OctDelaygsd0a600.ps', '-deps');