//\begin{verbatim}
// exec anchovy.sce 

// ANCHOVY
// parameters for the dynamical model

sex_ratio = 0.5 ; 
mature=sex_ratio*[1 1 1]; // proportion of matures at ages
weight=10^(-3)*[16, 28, 36]; // mean weights at ages (kg)
M=1.2; // natural mortality
F=[0.4 0.4 0.4]; // exploitation pattern
pi=0 ; // no plus-group
A=sum(ones(F)); // maximum age

// STOCK-RECRUITMENT RELATIONSHIPS

RR= 10^6 * [14016 7109 3964 696] ; 
// R_mean R_gm R_min 2002 (ICES) R_min 2004 (ICES)

// // CONSTANT STOCK-RECRUITMENT

function y=mini_constant(x)
    y= mini(RR) ;
endfunction

function y=mean_constant(x)
    y= RR(1) ;
endfunction


// // RICKER STOCK-RECRUITMENT 

a=0.79*10^6;
b=1.8*10^(-5);
// Ricker coefficients for tons

function y=Ricker(x)
  xx=10^{-3}*x // xx measured in tons
    y= a *( xx .* exp(-b* xx ) );
endfunction


// // LINEAR STOCK-RECRUITMENT

r= (500* 10^3) *21* 0.5 * 10^{-5} ;
function y=linear(x)
    y= r * x ;
endfunction


function y=SSB(N)
  // Spawning biomass
  y=  (mature.*weight) * N ;
endfunction

function Ndot=dynamics(N,lambda,phi)
  // Population dynamics
  mat=diag( exp(-M - lambda * F(1:($-1)) ) ,-1 ) +...
  diag( [zeros(F(1:($-1))) pi*exp(-M - lambda * F($)) ]);
	    // sub diagonal terms // diagonal terms
  Ndot= mat*N + [phi(SSB(N)) ; zeros(N(2:$)) ] ;
endfunction

// initial values

N1999=10^6*[4195 2079 217]';
N2000=10^6*[7035 1033 381]';
N2001=10^6*[6575 1632 163]';
N2002=10^6*[1406 1535 262]';
N2003=10^6*[1192 333 255]';
N2004=10^6*[2590 254 43]';
//  N0=10^6*[1379 1506 256]';

// ICES values

B_pa = 10^6 * 33 ; // kg
F_pa=1;
B_lim=21000 * 10^3 ; // kg
 // F_lim=Inf;
 
// TRAJECTORIES

T=10; // horizon in years
multiplier=0.5;

stock_recruitment=list();
stock_recruitment(1)=mini_constant;
stock_recruitment(2)=mean_constant;
stock_recruitment(3)=Ricker;
stock_recruitment(4)=linear;

for i=1:4 do

  phi=stock_recruitment(i) ; 
  // selecting a stock-recruitment relationship
  traj=[N1999];
  for t=0:(T-1)
	traj=[traj, dynamics(traj(:,$),multiplier,phi)] ; 
  end
// 
  total_biomass= weight * traj ; 
// 
  xset("window",i) ;  xbasc(i);
  plot2d(0:T,total_biomass,rect=[0,0,T,4*10^8])
  xtitle('Projections total biomass...
		 with catch pressure u(t)='...
		 +string(multiplier),'time t (years)',...
		 'total biomass B(t) (kg)')
  
  xset("window",10+i) ; xbasc(10+i);
  plot2d(0:T,traj',rect=[0,0,T,2*10^10]);
// drawing diamonds, crosses, etc. to identify the curves
  plot2d(0:T,traj',style=-[1,2,3]);
  legends(['ages 0--1';'ages 1--2';'ages 2--'],-[1,2,3],'ur')
  xtitle('Projections for the 3 age-classes...
		 with catch pressure u(t)='...
		 +string(multiplier),'time t (years)','abundances N_a(t)')
end

//\end{verbatim}