function simulate_array(hl,redeploy);

% simulate the float array
 

load transit_matrix lonmat latmat dT p
Lon=unique(lonmat); Lat=unique(latmat);

%hl=5*365; % half-life of "quit" floats (days); in practice,
        % the half-life will be shorter because of "death zones"

         
   
t=0:dT:5*12*dT;

% Ndead: fraction of drifters in box that don't survive next dT
% days ... due to dying, or box never sampled.
Ndead=1-reshape(full(sum(p,2)'),[length(Lon) length(Lat)])';
% Lots of bins have 0 deaths, so median isn't useful.
%  Global mean away from death zones (Ndead>=.05) is 0.0023.
%  This suggests that the half-life of "quit" floats is -dT/0.0023*log(0.5)=9051 days. 
%  That's not very realistic ... actual median lifetime is about 1262 days
%  for all floats in YoMaHa data set.
% For the simulations, impose a floor to give half-life hl.
% Let Ndead stay large where many drifters are dying, e.g. due
% to many run-agrounds. 
Ndead(find(Ndead<.2))=-dT/hl*log(.5); 


load eval_grid long* latg*

% only have floats between 5S-5N
%q=find(abs(latg)>5);long(q)=[]; latg(q)=[]; 
%q=find(abs(latg2)>5); long2(q)=[]; latg2(q)=[];

if (1)  % start with grid at long/latg resolution 
    lon=long2; lat=latg2;
else  % start with prerun simulation
    load sim_450 lon lat  
    lon=lon(:,size(lon,2)); lat=lat(:,size(lat,2));
    q=find(isfinite(lon)); lon=lon(q); lat=lat(q);
end;
 

% don't deploy (or even evaluate array) in death zones where odds
% of dying are >=5%. 
ndead=interp2([Lon;Lon+360],Lat,[Ndead,Ndead],long,latg,'linear');
q=find(ndead<0.05);
long=long(q); latg=latg(q);



lon2=NaN*ones(6000,length(t));
lat2=lon2;  
lon2(1:length(lon),1)=lon;
lat2(1:length(lat),1)=lat;
id(1:length(lon),1)=(1:length(lon))';
lon=lon2; clear lon2;
lat=lat2; clear lat2;

lon=360-139*ones(1000,1);
lat=-1*ones(1000,1);
id(1:length(lon),1)=(1:length(lon))';

N=zeros(size(t));
Nd=N; Nn=N;
N(1)=length(find(isfinite(lon(:,1))));  

dep_lon=[]; dep_lat=[]; dep_t=[];

clf;
plot(lon(:,1),lat(:,1),'.k');
axis([100 300 -20 20]); overlay;
pause(.1);
for j=2:length(t)
    ndead=interp2(Lon,Lat,Ndead,lon(:,j-1),lat(:,j-1),'linear');
    % if lon=NaN (not a live float), then ndead=NaN
    diceroll=rand(size(ndead))/1;
    for i=1:size(lon,1)
        if (ndead(i)<diceroll(i)) % it LIVES!!
            ind=( (lonmat-lon(i,j-1)).^2+(latmat-lat(i,j-1)).^2);
            ind=find(ind==min(ind),1);
            if (sum(p(ind,:)))
                r=rand(1,1)*sum(p(ind,:));
                ind2=find(cumsum(p(ind,:))>=r, 1 );
                lon(i,j)=lonmat(ind2); 
                lat(i,j)=latmat(ind2);
                id(i,j)=id(i,j-1);
            else
                lon(i,j)=lon(i,j-1); 
                lat(i,j)=lat(i,j-1);
                id(i,j)=id(i,j-1);
            end;
        else  % float died
            lon(i,j)=NaN; lat(i,j)=NaN;  id(i,j)=NaN;
            if (isfinite(ndead(i)))
                Nd(j)=Nd(j)+1;
            end;
        end; 
    end; 
    
    if (redeploy)  % add drifters as needed
        % identify distance between eval_grid and closest drifter
        del=NaN*ones(size(long));
        lon2=lon(:,j); lat2=lat(:,j);
        lon2=[lon2-360;lon2;lon2+360]; lat2=[lat2;lat2;lat2]; 
        for k=1:length(long);  % matricized version takes 1.6 times longer
            del(k)=min(sqrt( (lon2-long(k)).^2 + (lat2-latg(k)).^2));
        end;   
        %del=min(sqrt( (lon2*ones(1,length(long))-ones(length(lon2),1)*long').^2 ...
        %    +(lat2*ones(1,length(latg))-ones(length(lat2),1)*latg').^2 ))'; 
        
        % add floats where spacing is greater than 2.1 degrees 
        disp('adding floats where spacing is > 2.1 degrees')
        while (max(del)>2.1)
            q=find(del==max(del),1);  % deploy where need is greatest
            qq=find(isnan(id(:,j)));  % unused index to assign new float
            for k=1:length(q)
                lon(qq(k),j)=long(q(k));
                lat(qq(k),j)=latg(q(k));
                id(qq(k),j)=max(max(id))+1;
                Nn(j)=Nn(j)+1;
                dep_lon=[dep_lon;long(q(k))];
                dep_lat=[dep_lat;latg(q(k))];
                dep_t=[dep_t;t(j)]; 
            end;
            del2=NaN*ones(size(long));
            lon2=lon(:,j); lat2=lat(:,j);
            lon2=[lon2-360;lon2;lon2+360]; lat2=[lat2;lat2;lat2];
            for k=1:length(long);
                del2(k)=min(sqrt( (lon2-long(k)).^2 + (lat2-latg(k)).^2));
            end;
            del=del2;
        end;
         
    end;
    
    % count number of floats
    N(j)=length(find(isfinite(lon(:,j))));
    
    plot(lon(:,j),lat(:,j),'.k'); 
    title(['Day ',num2str(t(j)),': ',num2str(N(j)),' floats']); 
    axis([0 360 -90 90]); overlay;
    pause(.1);
    
end;

figure(1);  if (0)
for j=1:length(t);
    plot(lon(:,j),lat(:,j),'.k'); 
    title(['Day ',num2str(t(j)),': ',num2str(N(j)),' floats']); 
    axis([0 360 -90 90])
    overlay;
    pause(1);
end;        end;

figure(2); 
subplot(211)
plot(t,N,'k','linewidth',2);  % note that N=N(1)+cumsum(Nn-Nd)
title('Number of floats in array','fontsize',12);
subplot(212)
plot(t,Nd,'r',t,Nn,'b','linewidth',2); 
title('Number deployed (blue) and died (red)','fontsize',12);

figure(3);
wysiwyg 
q=find(dep_lon<30); dep_lon(q)=dep_lon(q)+360;
plot(dep_lon+.5*randn(size(dep_lon)), ...
    dep_lat+.5*randn(size(dep_lat)),'.b','markersize',6); 
dep_lon(q)=dep_lon(q)-360;
overlay;
mapaxis([30 30+360 -90 90])
title('Deployment locations','fontsize',12);


save sim lon lat t id N* hl dep_* Lon Lat long latg

return




% to create shaded plot of long, latg (eval) region: 
Lon2=unique(long); Lat2=unique(latg);
A=NaN*ones(length(Lat2),length(Lon2));
for i=1:length(Lat2);
    for j=1:length(Lon2);
        q=find(latg==Lat2(i) & long==Lon2(j));
        A(i,j)=length(q);
    end;
end;


Lon2=unique(long); Lat2=unique(latg);
Ndep=zeros(length(Lat2),length(Lon2));
for i=1:length(Lat2);
    q=find(dep_lat>Lat2(i)-2.5 & dep_lat<=Lat2(i)+2.5);
    for j=1:length(Lon2);
        q2=find(dep_lon(q)>Lon2(j)-2.5 & dep_lon(q)<=Lon2(j)+2.5);
        Ndep(i,j)=length(q2);
    end;
end;
% rescale to be number deployed per year  
Ndep=Ndep/max(t)*365.24;

