
%%%%%%%%%%%%%%% How do we get traveling waves? %%%%%%%%%%%%%

%%%%%%%%%%%% Choose your spatial C profile %%%%%%%%%%%%%%%%%

Profile  = 3;     %%%%%%% Choose between 1, 2 or 3

%%%% Input parameters (keep one fixed and vary the other) %%%%%%

gamma  = 0.5; % Slope of the ramp
DeltaX = 40; % Width of the ramp


%%%%%%%%%%%%% Do not change anything below %%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%% Simulation parameters %%%%%%%%%%%%%%%%%

T = 40;
Tau = 2.0;
H = 5;
Co = 2;

%%%%% Time steps and resolution %%%%%%%%%
N = 15000;
M = 100;
dt = 0.01;





%%%%%%%%%%%% Setting up the Ramp function %%%%%%%%%%

time = (1:N+1)'*dt;

C = gamma*(heaviside(time-10)-heaviside(time-T)).*(time-10);

%%%%%%% Initial conditions %%%%%%%%%%

w = zeros(Tau/dt+1,M);

%%%%%%%%%%%%% Defining the C in the tailbud %%%%%%%%%%%

v  = 0.01; %%%%%%% Traveling speed of front

Ctb = Co*(1-heaviside((1:M)-M+v*(1:N+1)'));
Cpsm = (Co + gamma*((1:M)-M+v*(1:N+1)')).*(heaviside((1:M)-M+v*(1:N+1)') -  heaviside((1:M)-M-DeltaX+v*(1:N+1)'));



for i = Tau/dt+1 : N
    
    if Profile == 1
        
        C(i,1:M) = Ctb(i,1:M);
        
    else if Profile == 2
            
            C(i,1:M) = Cpsm(i,1:M);
            
        else if Profile == 3
                
                C(i,1:M) = Ctb(i,1:M) + Cpsm(i,1:M);
                
            end
        end
    end
    
    w(i+1,1:M) = (-w(i,1:M)+C(i,1:M).*(1+w(i-Tau/dt,1:M).^H).^-1)*dt + w(i,1:M);
    
end




video = VideoWriter('Traveling_test.avi');
video.FrameRate = 20;
open(video);



figure
for j = 1:25:N
    % Example of plot
    
subplot(2,1,1)
plot(w(j,:),'-b','LineWidth',2 )
set(gca,'fontsize', 16)
ylabel('w')
xlim([1 100])
ylim([0 max(max(w))+.5])
subplot(2,1,2)
plot(C(j,:),'-r','LineWidth',2)
set(gca,'fontsize', 16)
xlabel('Cell Number')
ylabel('C')

xlim([1 100])
ylim([0 max(max(C))])
    
    % Store the frame
    frame = getframe(gcf);
    writeVideo(video,frame);
   
end



close(video)
close(gcf)



video = VideoWriter('Traveling_test2.avi');
video.FrameRate = 20;
open(video);


 mymap_green = horzcat(zeros(1001,1),(0:.001:1)',zeros(1001,1));

figure
for j = 1:25:N
    % Example of plot
    
imagesc(repmat(w(j,:),30,1)')
set(gca,'YDir','Normal')
set(gca,'xtick',[])
set(gca,'xticklabel',[])
set(gca,'ytick',[])
set(gca,'yticklabel',[])
set(gcf,'units','points','position',[10,10,250,800])
caxis([0 max(max(w))/2+.5])
colormap(mymap_green);
    
    % Store the frame
    frame = getframe(gcf);
    writeVideo(video,frame);
   
end



close(video)
close(gcf)

 


