// \begin{verbatim}
// PARAMETERS // 
 
// initial time
t_0=1990; 
// final Time
t_F=2100; 
// Time step
delta_t=1; 

mean_taux_Q=0.03;						
sig_taux_Q=0.03; 	
// mean and standard deviation of growth rate 
alphaa=0.64;											
// dynamics parameter: atmospheric retention 
// (uncertain +- 0.15)	
sigma=0.519;
deltaa=1/120;

// concentration target (ppm)
M_sup=550;
M_inf=274;
// Initial conditions  
t=t_0;
M0=354; //in (ppm)
M_bau=M0; M_g=M0;
Q0 = 20.9; // in (T US$)
	
function E=emissions(Q,a)   						
	E = sigma * Q * (1-a);								
	endfunction
function Mnext=dynamics(t,M,Q,a)   						
	E = emissions(Q,a) ;								
	Mnext =M  + alphaa* E -deltaa*(M-M_inf);			
	 endfunction
// Distinct abatment policies
	a = 1*ones(1,t_F-t_0+1);   
	// Strong mitigation 
	a = 0*ones(1,t_F-t_0+1);  
	// No mitigation (BAU)
	a = 0.9*ones(1,t_F-t_0+1);  
	// medium mitigation  	
//a = 1*rand(1,t_F-t_0+1);  
// random mitigation  

xbasc(1);xbasc(2);xbasc(4)

// System Dynamics 
N_simu=3;
for ii=1:N_simu
  t=t_0;
  M=M0; M_bau=M; M_g=M;
  Q = Q0;
  //Initialisation
  L_t=[t_0]; 
  L_M=[M0]; L_bau=[M0]; L_g=[M0];
  L_Q=[];
  L_E=[]; L_Ebau=[];L_Eg=[]; 

  for (t=t_0:delta_t:t_F)	
	g_Q=sig_taux_Q*2*(rand(1,1)-0.5)+mean_taux_Q; 
	// random growth rate 
	E = emissions(Q,a(t-t_0+1));		
	M=dynamics(t,M,Q,a(t-t_0+1));
	// Mitigation concentration
    E_bau= emissions(Q,0);		
	M_bau=dynamics(t,M_bau,Q,0);
	//  Business as usual (BAU)
	E_g = emissions(Q,1);		
	M_g=dynamics(t,M_g,Q,1);			
	// total abatement
	Q=(1+g_Q)*Q;	
	L_Q=[L_Q Q];
	L_t=[L_t t+1];
	L_M=[L_M M];
	L_E=[L_E E];
	L_bau=[L_bau M_bau];
	L_Ebau=[L_Ebau E_bau];
	L_g=[L_g M_g];
	L_Eg=[L_Eg E_g];
  end,

// Results printing

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

xset("window",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], rect=[t_0,0,t_F,2000]) ;
legends(["Mitigation";"BAU";"green";"threshold"],...
      -[4,5,3,-1],'ul');
xtitle('Concentration CO2','t','M(t) (ppm)');

xset("window",4);
plot2d(L_t(abcisse),L_Q(abcisse),style=-[8]);
xtitle('Economie: Production Q(t)','t','Q(t) (T US$)'); 

end
// \end{verbatim}