clear all
close all
dt=.1;
nfish=30; %% number of fish
fishdensity=12; %% starting area per fish
maxv=.9; %% starting velocity
maxturn=pi/4*dt;%% maximum random turn from previous direction by followers
lmaxturn=pi/2*dt; %% ditto by leader
startangle=2*pi; %% can start facing anywhere 0-this angle
maxa=.1*maxv; %% max acceleration
fish=zeros(nfish,4); %% components: direction velocity x y
timetotal=120;
idist=10; %% max distance at that followers check for influences

for i=1:nfish;
    fish(i,1)=rand()*startangle;
    fish(i,2)=.8*maxv+.2*rand()*maxv;
    fish(i,3)=rand()*sqrt(nfish*fishdensity)-sqrt(fishdensity*nfish)/2;
    fish(i,4)=rand()*sqrt(nfish*fishdensity)-sqrt(fishdensity*nfish)/2;
end


figure(1)
filename='fishschool.gif';


for it=1:timetotal/dt;
    
    for j=2:nfish %% slight randomization of direction and velocity
        fish(j,2)=abs(fish(j,2)-maxa+rand()*2*maxa);
        fish(j,1)=fish(j,1)-maxturn+rand()*2*maxturn; 
        if j==nfish %% leader changed
            fish(1,1)=(fish(1,1)-lmaxturn+rand()*2*lmaxturn);
            fish(1,2)=abs(fish(1,2)-maxa+rand()*2*maxa);
        end
        fish(1,1)=fish(1,1);
    end
    for j=1:nfish %% If velocity is too low move it up
        if fish(j,2)<.6*maxv
            fish(j,2)=.6*maxv;
        end
    end
    for j=2:nfish;
        velocity=0;
        angles=0;
        nearfish=0;
        fdistances=zeros(2,nfish);
        cfishloc=[fish(j,3) fish(j,4)];
        for n=1:nfish;
            nfishloc=[fish(n,3) fish(n,4)];
            if abs(sum(cfishloc-nfishloc))<2*idist %% filter out very far away fish (speeds this up)
                fdistances(n)=norm(cfishloc-nfishloc); %% calculates distance away
                
            else
                fdistances(n)=2*idist;
            end
        end
        for n=1:nfish
            if fdistances(n)<=idist %% examines nearby fish in influence zone
                angles=angles+fish(n,1); 
                velocity=velocity+fish(n,2);
                nearfish=nearfish+1;
            end
        end
        fish(j,1)=((angles)/nearfish); %% finds average direction/velocity
        fish(j,2)=(velocity/nearfish);
    end
    for j=1:nfish
        fish(j,3)=fish(j,3)+fish(j,2)*cos(fish(j,1))*dt; %% calculates new position
        fish(j,4)=fish(j,4)+fish(j,2)*sin(fish(j,1))*dt;
    end
    
    if mod(it,9)==1 %% Plots and makes GIF
    plot(fish(2:nfish,3),fish(2:nfish,4),'k.', ...
    fish(1,3),fish(1,4),'r.','markersize',10);
    xlim([-120 120]);
    ylim([-120 120]);
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            f=12;
    frame = getframe(1);
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,256);
    if it == 1;
    imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
    else
    imwrite(imind,cm,filename,'gif','WriteMode','append');
    end
    saveas(gcf,['wave_over_ridge',sprintf('%03d',f),'.png']) 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
end
        
        
        
        
            