import com.comsol.model.*
import com.comsol.model.util.*

set(0,'DefaultAxesFontSize',15.0);
set(0,'DefaultTextFontSize',15.0);

% Now connect Comsol to Matlab server

model=mphload('Nagumo_1D_Traveling1Front.mph');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         BILDER                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% (a) Data Sets:
% 'dset1' <- Nonfrozen Equation/Solution 1
% 'dset2' <- Nonfrozen Equation/Solution Store (Nonfrozen Solution)
% 'dset3' <- Frozen Equation/Solution 3
% 'dset4' <- Frozen Equation/Solution Store (Frozen Solution)
% 'dset5' <- Frozen Equation/Solution Store (Eigenvalues and Eigenfunctions):
% (b) Dependent Variables:
% 'u'      <- dset2
% 'v'      <- dset4
% 'g1'     <- dset4
% 'mu1'    <- dset4
% 'w'      <- dset5
% 'lambda' <- dset5

% a) Nonfrozen solution
u=mpheval(model,'u','dataset','dset2');
solinfo1=mphsolinfo(model,'dataset','dset2');
time1=solinfo1.solvals;
space1=u.p;

% b) Frozen solution
v=mpheval(model,'v','dataset','dset4');
mu1=mphevalpoint(model,'mu1','selection',1,'dataset','dset4');
mu1t=mphevalpoint(model,'mu1t','selection',1,'dataset','dset4');
g1=mphevalpoint(model,'g1','selection',1,'dataset','dset4');
l2vt=mphevalpoint(model,'sqrt(sqr_vt)','selection',1,'dataset','dset4');
solinfo2=mphsolinfo(model,'dataset','dset4');
time2=solinfo2.solvals;
space2=v.p;

% c) Solution of eigenvalue problem
w=mpheval(model,'v','dataset','dset5');
solinfo3=mphsolinfo(model,'dataset','dset5');
lambda=solinfo3.solvals;
space3=w.p;

%  ------------------------------
% | a) Nonfrozen Nagumo equation |
%  ------------------------------

% Figure 01: (Traveling Front)
figure
hold on
plot(u.p,u.d1(1,:),'b','Linewidth',3);
plot(u.p,u.d1(401,:),'g','Linewidth',3);
plot(u.p,u.d1(601,:),'r','Linewidth',3);
plot(u.p,u.d1(801,:),'c','Linewidth',3);
hold off
xlabel('x');
ylabel('u(x,t)');
legend('t=0','t=40','t=60','t=80','Location','NorthEast')
axis tight
%File -> Export setup -> Expand axes to fill figure
print('-dpng','-r100','Nagumo_1D_Traveling1Front_Nonfrozen_TravelingFront')
close(gcf)

% Figure 02: (Space time plot of u)
figure
surf(u.p,time1,u.d1,'EdgeColor','none');
colormap jet
shading interp
%set(gcf,'renderer','zbuffer');
axis([space1(1) space1(end) time1(1) time1(end) 0 1]);
xlabel('x'); ylabel('t'); zlabel('u(x,t)');
%File -> Export setup -> Expand axes to fill figure
print('-dpng','-r100','Nagumo_1D_Traveling1Front_Nonfrozen_SpaceTime_3DViewA')
view([0,90]);
print('-dpng','-r100','Nagumo_1D_Traveling1Front_Nonfrozen_SpaceTime')
close(gcf)

%  ---------------------------
% | b) Frozen Nagumo equation |
%  ---------------------------

% Figure 03: (Profile)
figure
plot(space2,v.d1(end,:),'b','Linewidth',3);
axis([space2(1) space2(end) 0 1]);
xlabel('x'); ylabel(['v(x,',num2str(time1(end)),')']);
%File -> Export setup -> Expand axes to fill figure
print('-dpng','-r100','Nagumo_1D_Traveling1Front_Frozen_Profile')
close(gcf)

% Figure 04: (Space time plot of profile)
figure
surf(v.p,time2,v.d1,'EdgeColor','none');
colormap jet
shading interp
%set(gcf,'renderer','zbuffer');
axis([space2(1) space2(end) time2(1) time2(end) 0 1]);
xlabel('x'); ylabel('t'); zlabel('v(x,t)');
%File -> Export setup -> Expand axes to fill figure
print('-dpng','-r100','Nagumo_1D_Traveling1Front_Frozen_SpaceTime_3DView')
view([0,90]);
print('-dpng','-r100','Nagumo_1D_Traveling1Front_Frozen_SpaceTime')
close(gcf)

% Figure 05: (Velocities)
figure;
plot(time2(1:501),mu1(1:501),'b-','LineWidth',3);
xlabel('t'); ylabel('\mu(t)');
axis([time2(1) time2(501) -0.4 -0.1]);
%File -> Export setup -> Expand axes to fill figure
print('-dpng','-r100','Nagumo_1D_Traveling1Front_Frozen_Velocity')
close(gcf)

% Figure 06: (Position)
figure;
plot(g1,time2,'b-','LineWidth',3);
xlabel('\gamma(t)'); ylabel('t');
axis([space2(1) space2(end) time2(1) time2(end)]);
%File -> Export setup -> Expand axes to fill figure
print('-dpng','-r100','Nagumo_1D_Traveling1Front_Frozen_Velocity')
close(gcf)

% Figure 07: (L^2-norm of vt and absolute value of mu1t)
figure;
hold on
plot(time2,l2vt,'b-','LineWidth',3);
plot(time2,abs(mu1t),'r-','LineWidth',3);
hold off
set(gca, 'yscale', 'log');
xlabel('t');
legend('||v_t||_{L^2}','|\mu_{1,t}|','Location','NorthEast');
axis([time2(1) time2(end) 10^(-16) 2]);
%File -> Export setup -> Expand axes to fill figure
print('-dpng','-r100','Nagumo_1D_Traveling1Front_Frozen_Convergence')
close(gcf)

%  -------------------------------------------
% | c) Eigenvalue Problem for Nagumo equation |
%  -------------------------------------------

% Figure 08: (Eigenvalues and spectrum)
etol=0.01;
sb=-0.25;  %spectral bound of Df(v_+/-)
figure;
hold on;
lambdaN=-lambda;
lr=real(lambdaN);
li=imag(lambdaN);
plot([0 0],[-1 1],'k','Linewidth',2);
plot([-2 0.5],[0 0],'k','Linewidth',2);

lam_pt_sym=[]; lam_pt=[]; lam_ess=[]; 

for i=1:length(lr)
  if abs(lr(i))<=etol
   lam_pt_sym=[lam_pt_sym,lambdaN(i)];
  elseif lr(i)<=(sb+etol)
   lam_ess=[lam_ess,lambdaN(i)];
  else
   lam_pt=[lam_pt,lambdaN(i)];
  end
end

% a) point spectrum (on the imaginary axis)
plot(real(lam_pt_sym),imag(lam_pt_sym),'marker','o','Markersize',15,'linewidth',3,'linestyle','none','color','b'); 
% b) point spectrum (not on the imaginary axis)
plot(real(lam_pt),imag(lam_pt),'marker','+','Markersize',15,'linewidth',3,'linestyle','none','color','b'); 
% c) essential spectrum
plot(real(lam_ess),imag(lam_ess),'marker','.','Markersize',20,'linewidth',3,'linestyle','none','color','r'); 
hold off;
xlabel('Re \lambda'); ylabel('Im \lambda');
axis([-2 0.5 -1 1]);
%File -> Export setup -> Expand axes to fill figure
print('-dpng','-r100','Nagumo_1D_Traveling1Front_Frozen_Spectrum')
close(gcf)

% Figure 10: (Eigenfunctions)
figure
solnum=1;
plot(space,w.d1(solnum,:),'b','Linewidth',2);
xlabel('x');
ylabel('w(x)');
title(num2str(lambda(solnum)));

%  -----------------------
% | d) Movies (Nonfrozen) |
%  -----------------------
set(0,'DefaultAxesFontSize',15.0);
set(0,'DefaultTextFontSize',15.0);
       
colors={'b','r','g','k'};
linestyles={'-','-','-','-'};

w1 = 1024; w2 = 768;
figure('Renderer','zbuffer','Nextplot','ReplaceChildren', ...
           'PaperPositionMode','auto','Position',[100 100 w1 w2]);
sec=10;
fps=10;
frames=sec*fps; %frames
ats=length(time1); %all time steps
rts=floor(ats/frames); %reduced time steps

writerObj = VideoWriter('Nagumo_1D_Traveling1Front.avi');
writerObj.FrameRate = fps;
writerObj.Quality = 99;
open(writerObj); 
set(gca,'nextplot','replacechildren');
set(gcf,'Renderer','zbuffer')
 
for j = 1:(frames+1)
    % Plot 1
    plot(u.p,u.d1(1+(j-1)*rts,:),'b','Linewidth',3);
    axis([-50 50 -0.05 1.05])
    xlabel('x');
    ylabel('u(x,t)');
    
    % Plot 2: Space-Time
    %if (j>=2)
    %surf(u.p,time1(1:1+(j-1)*rts),u.d1(1:1+(j-1)*rts,:),'EdgeColor','none');
    %end
    %colormap jet
    %shading interp
    %colorbar
    %caxis([-0.05 1.05])
    %axis([space1(1) space1(end) time1(1) time1(end) -0.05 1.05]);
    %xlabel('x'); ylabel('t'); zlabel('u(x,t)');
    %view([0,90])
    
    title(['\fontsize{15}Nagumo Equation, Traveling Front, t=',num2str(time1(1+(j-1)*rts),'%5.1f')]);
    frame = getframe(gcf);
    writeVideo(writerObj,frame);
end
close(writerObj);
close(gcf);

%  --------------------
% | e) Movies (Frozen) |
%  --------------------
set(0,'DefaultAxesFontSize',15.0);
set(0,'DefaultTextFontSize',15.0);
       
colors={'b','r','g','k'};
linestyles={'-','-','-','-'};

w1 = 1024; w2 = 768;
hf=figure('Renderer','zbuffer','Nextplot','ReplaceChildren', ...
           'PaperPositionMode','auto','Position',[100 100 w1 w2]);
hwf = subplot(2,2,1);
hlf = subplot(2,2,2);
hcf = subplot(2,2,3);
hpf = subplot(2,2,4);
pos1=get(hwf,'Position');
pos2=get(hlf,'Position');
pos3=get(hcf,'Position');
pos4=get(hpf,'Position');
shift=[0 -0.02 0 0];
set(hwf,'Position',pos1+shift);
set(hlf,'Position',pos2+shift);
set(hcf,'Position',pos3+shift);
set(hpf,'Position',pos4+shift);
sec=10;
fps=10;
frames=sec*fps; %frames
ats=length(time2); %all time steps
rts=floor(ats/frames); %reduced time steps

writerObj = VideoWriter('Nagumo_1D_Frozen1Front.avi');
writerObj.FrameRate = fps;
writerObj.Quality = 99;
open(writerObj); 
set(gca,'nextplot','replacechildren');
set(gcf,'Renderer','zbuffer')
 
for j = 1:(frames+1)
    % Profile
    axes(hwf);
    plot(v.p,v.d1(1+(j-1)*rts,:),'b','Linewidth',3);
    title('Profile');
    axis([-50 50 -0.05 1.05])
    xlabel('x');
    ylabel('v(x,t)');
    axis square
    % Velocity
    axes(hlf);
    plot(time2(1:1+(j-1)*rts),mu1(1:1+(j-1)*rts),'r-','LineWidth',3);
    title('Velocity');
    axis([time2(1) time2(end) -0.4 -0.1]);
    xlabel('t');
    ylabel('\mu(t)');
    axis square
    % Position
    axes(hpf);
    plot(g1(1:1+(j-1)*rts),time2(1:1+(j-1)*rts),'k-','LineWidth',3);
    title('Position');
    axis([space2(1) space2(end) time2(1) time2(end)]);
    xlabel('\gamma(t)');
    ylabel('t');
    axis square
    % Convergence (of Profile and Velocity)
    axes(hcf);
    hold on
    plot(time2(1:1+(j-1)*rts),l2vt(1:1+(j-1)*rts),'b-','LineWidth',3);
    plot(time2(1:1+(j-1)*rts),abs(mu1t(1:1+(j-1)*rts)),'r-','LineWidth',3);
    hold off
    title('Convergence');
    set(gca, 'yscale', 'log');
    xlabel('t');
    legend('||v_t||_{L^2}','|\mu_{t}|','Location','NorthEast');
    axis([time2(1) time2(end) 10^(-16) 2]);
    box on
    axis square
    
    % Titel                                                                                    ');
    [ax,h]=subtitle(['Nagumo Equation, Frozen 1-Front, t=',num2str(time1(1+(j-1)*rts),'%5.1f')]);
    
    frame = getframe(gcf);
    writeVideo(writerObj,frame);
    
    axes(hwf); cla;
    axes(hlf); cla;
    axes(hpf); cla;
    axes(hcf); cla;
    delete(h);
end
close(writerObj);
close(gcf);

%%% USE SUBTITLE FUNCTION:
%%%

function [ax,h]=subtitle(text)
%
%Centers a title over a group of subplots.
%Returns a handle to the title and the handle to an axis.
% [ax,h]=subtitle(text)
%           returns handles to both the axis and the title.
% ax=subtitle(text)
%           returns a handle to the axis only.

ax=axes('Units','normalized','Position',[.075 .1 .85 .85],'Visible','off');
set(get(ax,'Title'),'Visible','on')
title(text,'fontsize',16.5);
if (nargout < 2)
    return
end
h=get(ax,'Title');
end
