function [observed, totalpop] = gensalmondata();

run('marshmat') %make sure the mat file is up to date
load('marshmat')

burn = 20;
simlen = 70;
agestruc = 0; %age struc == 2 means put in perturbation (0 reproduction one year)
densdep = 1;  %means add density dependence
matsize=length(stableagestruc);
census = ones(matsize,simlen+burn);
census(:,1) = 100000*stableagestruc; %stable and pulse age strucs defined in the marshmat.m, snakefall.m and upcolsteelhead.m files

totalpops=zeros(1,simlen);

for i = 2:(burn+simlen),  
   if(agestruc==2 & i == (burn+2)) census(1,i-1)=0; end %perturbation means no reproduction
	e1=mymvnrnd(length(evarstream), 0, evarstream, 0.2);
   	e2=mymvnrnd(length(evarocean), 0, evarocean, 0.8);
   e=[e1 e2];
   stochmat = stochplace;
   stochmat(stochmat>0)=e(stochmat(stochmat>0)); %this assigns the e variable to the places denoted by stochplace
   if(densdep & i>11) 
   		aveeggdens = mean(census(end,(i-11):(i-2)));
   		prevreleggdens = min(census(end,i-1)/aveeggdens, 3.33);
   		prevreleggdens = max(prevreleggdens,0.333);
   		stochmat(:,end)=stochmat(:,end)*((prevreleggdens-1)*(-0.782/3)+1);
   end
   thisyearmat = themat .* stochmat;
   tmp=thisyearmat(2:matsize,:);
   tmp(tmp>1)=1;
   thisyearmat(2:matsize,:)=tmp;
   
   census(:,i) = thisyearmat*census(:,i-1);
end % i years

%totalpop is the number of eventual spawners alive at time t
totalpop = sum(census(:,(burn+1):end).*(weight'*ones(1,simlen)),1);
observed = sum(census(obsstage,(burn+1):end),1);

sampleerrorsigma = .1*ones(1,simlen); bias=0; 
samperror = normrnd(bias,sampleerrorsigma);
samperror = exp(samperror); % make it lognormal

%observed = observed.*samperror; %add sampling error to spawners count