//\begin{verbatim}
clear


// PARAMETERS // 
 
//  initial time
t_0=1990; 
// Final Time
t_F=2100; 
// Time step
delta_t=1; 

taux_Q=0.01;										    
// economic growth rate
alphaa=0.64;											
// marginal ratio marginal atmospheric retention
// (uncertain +- 0.15)	
sigma=0.519;
absortion=1/120 ;
// concentration target (ppm)
M_sup=550;
// Initial conditions  
t=t_0;
M=354; //in (ppm)
M_bau=M; M_g=M;
Q = 20.9; // in (T US$)
E = sigma * Q ;		



// Distinct abatment policies
	u = 1*ones(1,t_F-t_0+1);   // Strong mitigation 
	u = 0*ones(1,t_F-t_0+1);   // No mitigation (BAU)
	u = 0.6*ones(1,t_F-t_0+1);  // medium mitigation  	
//u = 1*rand(1,t_F-t_0+1);  // random mitigation  

// Initialisation (empty lists)
L_t=[]; L_M=[]; L_bau=[]; L_E=[];
L_Eg=[];L_Eb=[]; L_Q=[];L_g=[];
// System Dynamics 
for (t=t_0:delta_t:t_F)
		
	L_t=[L_t t];
   L_M=[L_M M];
 L_Q=[L_Q Q];
	L_bau=[L_bau M_bau];
	L_g=[L_g M_g];

	Q=(1+taux_Q)*Q;							
	E = sigma * Q * (1-u(t-t_0+1));				 
	L_E=[L_E E];
	// Emissions CO2
	M = M* (1-absortion) + alphaa* E;						
	// dynamics concentration CO2
	
	E_bau = sigma * Q ;								
	L_Eb=[L_Eb E_bau];
	// Emissions Business as usual (BAU)
	M_bau = M_bau* (1-absortion) + alphaa* E_bau;			
	// dynamics BAU
	
	E_g = 0;			 
	L_Eg=[L_Eg E_g];					
	// Green: no emissions
	M_g = M_g* (1-absortion) + alphaa* E_g;							
	// dynamics without pollution
end, 
 
// Results printing

long=prod(size(L_t)); 
step=floor(long/20); 
abcisse=1:step:long;
xset("window",1);xbasc(1)
plot2d(L_t(abcisse),[L_E(abcisse)' L_Eb(abcisse)' ...
    L_Eg(abcisse)'],style=-[4,5,3]) ;
legends(["viable";"BAU";"green"],-[4,5,3],'ul');
xtitle('Emissions E(t)','t','E(t) (GtC)');

xset("window",2);xbasc(2)
plot2d(L_t(abcisse),[L_M(abcisse)' L_bau(abcisse)' ...
					L_g(abcisse)' ones(L_t(abcisse))'*M_sup],...
	   style=-[4,5,3,-1]) ;
legends(["viable";"BAU";"green";"threshold"],...
      -[4,5,3,-1],'ul');
xtitle('Concentration CO2','t','M(t) (ppm)');

xset("window",4); xbasc(4)
plot2d(L_t(abcisse),L_Q(abcisse)); 
xtitle('Economie: Production Q(t)','t','Q(t) (T US$)'); 

//\end{verbatim}