%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "MINT"
% Written for GNU Octave
%
% INSTITUTE OF TELECOMMUNICATIONS
% University of Stuttgart
% www.inue-uni-stuttgart.de
% author: Maximilian Arnold
% date: 28.6.2017

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [] = Interference_Map(Pos_X_BS1,Reflection,Surface_Map,webdemo_en,svg_en,filename)
%% Linewidth Arial
set(0, "defaultlinelinewidth", 2);

if (webdemo_en)
    set(0, "defaultaxesfontsize", 8);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
  elseif (svg_en)
    set(0, "defaultaxesfontsize", 6);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 6);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 1;
  else
    set(0, "defaultaxesfontsize", 10);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
  endif

%% Definition of Domain
% Positions
Position_X_1 = 20;
Position_Y_1 = 40;
Position_X_2 = Pos_X_BS1;
Position_Y_2 = 50;
Pos_X_User1  = 54;
Pos_Y_User1  = 55;

% Space
Lx = 12;
Ly = 10;
dx = 0.1;
dy = dx;
nx = fix(Lx/dx);
ny = fix(Ly/dy);
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);

% Time
T = 40;

% Field variable
wn = zeros(nx,ny);
wnm1 = wn; % w at time n-1
wnp1 = wn; % w at time n+1

% Set Settings for Wall Reflection
if(Reflection==0 & Surface_Map ~=2)
    % Switch User and Plot Wall
    Pos_X_User1 = 100;
    Pos_Y_User1 = 50;
end

% Parameters
CFL = 0.2; % CFL = c.dt/dx
c = 1;
dt = CFL*dx/c;
f  = 1;

%% Go over time propagation
t=0;
index =0;
counter =1;

if(Surface_Map~=2)
    while(t< T)
        t = t+dt;
        wnm1= wn; % Save current and previous arrays
        wn=wnp1;  % Save current and previous arrays
        
        % Source
        if (Reflection==1)
            wn(Position_X_1,Position_Y_1) = dt^1.5*20*sin(f*pi*t);
            wn(Position_X_2,Position_Y_2) = dt^1.5*20*sin(f*pi*t);
        else
            wn(Position_X_2,50) = dt^1.5*20*sin(f*pi*t);
        end
        
        % Calculate new Interference_Map
        wnp1 = 2*wn - wnm1 + CFL^2 * (circshift(wn,[1,0])+circshift(wn,[0,1])-4*wn + circshift(wn,[-1,0]) + circshift(wn,[0,-1]));
        
        % Absorbing Boundary Condition besides the right Wall
        wnp1(1,:)=wn(2,:) + ((CFL-1)/(CFL+1))*(wnp1(2,:)-wn(1,:));
        wnp1(:,1) = wn(:,2) + ((CFL-1)/(CFL+1))*(wnp1(:,2)-wn(:,1));
        wnp1(:,end) = wn(:,end-1) + ((CFL-1)/(CFL+1)) * (wnp1(:,end-1)-wn(:,end));
        if (Reflection ==1)
            % Absorbing Boundary for the right Wall
            wnp1(end,:) = wn(end-1,:) + ((CFL-1)/(CFL+1))*(wnp1(end-1,:)-wn(end,:));
        end
        
        % Save this plot each n-th Point
        if(mod(counter,12) ==1)
            % Create New Filname
            newFilename = [filename '_' num2str(index)];
            index = index+1;
            % Create Plot
            surf(abs(10*wn')); hold on;
            if(Reflection==0)                
                plot([120, 120],[0 , 1000000],'linewidth',10,'color','b');
                text (105,80, "Wand",'color','b');
            end
            shading interp;
            if(Surface_Map)
                grid off;
            else
                view([0 90]);
            end
            setGraphSettings();
            
            text (Pos_X_User1+2,Pos_Y_User1+3, "Nutzer",'color','r');
            plot(Pos_X_User1,Pos_Y_User1,'marker','x','markersize',5,'color','r');
            title(['Interferenzmuster [t = ' num2str(t) ' ns ,  ' '\lambda' '= 20 cm, Betragsquadrat]']);
            print(newFilename, '-dpng','-S640,480');
            hold off;
        end
        counter = counter+1
    end
else
    % Field variable
    wn = zeros(nx,ny);
    wnm1 = wn; % w at time n-1
    wnp1 = wn; % w at time n+1
    t=0;
    counter =1;
    while(t< T)
        t = t+dt;
        wnm1= wn; % Save current and previous arrays
        wn=wnp1;  % Save current and previous arrays
        
        % Source
        wn(Position_X_1,Position_Y_1) = sin(f*pi*t);
        
        % Calculate new Interference_Map
        wnp1 = 2*wn - wnm1 + CFL^2 * (circshift(wn,[1,0])+circshift(wn,[0,1])-4*wn + circshift(wn,[-1,0]) + circshift(wn,[0,-1]));
        
        % Absorbing Boundary Condition
        wnp1(1,:)=wn(2,:) + ((CFL-1)/(CFL+1))*(wnp1(2,:)-wn(1,:));
        wnp1(:,1) = wn(:,2) + ((CFL-1)/(CFL+1))*(wnp1(:,2)-wn(:,1));
        wnp1(:,end) = wn(:,end-1) + ((CFL-1)/(CFL+1)) * (wnp1(:,end-1)-wn(:,end));
        wnp1(end,:) = wn(end-1,:) + ((CFL-1)/(CFL+1))*(wnp1(end-1,:)-wn(end,:));
        
        % Save Away User Data
        Source_1(counter,1) = t;
        Source_1(counter,2) = wnp1(Pos_X_User1,Pos_Y_User1);
        counter = counter +1;
    end
    
    % Field variable
    wn = zeros(nx,ny);
    wnm1 = wn; % w at time n-1
    wnp1 = wn; % w at time n+1
    t=0;
    counter =1;
    while(t< T)
        t = t+dt;
        wnm1= wn; % Save current and previous arrays
        wn=wnp1;  % Save current and previous arrays
        
        % Source
        wn(Position_X_2,Position_Y_2) = sin(f*pi*t);
        
        % Calculate new Interference_Map
        wnp1 = 2*wn - wnm1 + CFL^2 * (circshift(wn,[1,0])+circshift(wn,[0,1])-4*wn + circshift(wn,[-1,0]) + circshift(wn,[0,-1]));
        
        % Absorbing Boundary Condition
        wnp1(1,:)=wn(2,:) + ((CFL-1)/(CFL+1))*(wnp1(2,:)-wn(1,:));
        wnp1(:,1) = wn(:,2) + ((CFL-1)/(CFL+1))*(wnp1(:,2)-wn(:,1));
        wnp1(:,end) = wn(:,end-1) + ((CFL-1)/(CFL+1)) * (wnp1(:,end-1)-wn(:,end));
        wnp1(end,:) = wn(end-1,:) + ((CFL-1)/(CFL+1))*(wnp1(end-1,:)-wn(end,:));
        
        % Save Away User Data
        Source_2(counter,1) = t;
        Source_2(counter,2) = wnp1(Pos_X_User1,Pos_Y_User1);
        counter = counter +1;
    end
    % Scale to 1
    if(Pos_X_BS1>72)
    Source_2(find(Source_2(:,1)<3.8& Source_2(:,1)>0),2) = 0.9*Source_2(find(Source_2(:,1)<3.8 & Source_2(:,1)>0),2);
      Source_1(find(Source_1(:,1)<4.9 & Source_1(:,1)>3.35),2) = 0.9*Source_1(find(Source_1(:,1)<4.9 & Source_1(:,1)>3.35),2);
    else
    Source_1(find(Source_1(:,1)<5.05 & Source_1(:,1)>3.35),2) = 0.9*Source_1(find(Source_1(:,1)<5.05 & Source_1(:,1)>3.35),2);
    Source_2(find(Source_2(:,1)<3& Source_2(:,1)>0),2) = 0.92*Source_2(find(Source_2(:,1)<3 & Source_2(:,1)>0),2);
    end
    Scale = max([Source_1(:,2)+Source_2(:,2);Source_1(:,2);Source_2(:,2)]);
    for i=0:60
       newFilename = [filename '_' num2str(i)];
       plot(Source_2(:,1),Source_2(:,2)./Scale,'b','linewidth',1.5);hold on;
       plot(Source_1(:,1),Source_1(:,2)./Scale,'g','linewidth',1.5);
       plot(Source_2(:,1),(Source_1(:,2) +Source_2(:,2))./Scale,'r','linewidth',1.5);hold on
       xlim([0 30]);
       ylim([-1.35 1.35]);
       xlabel('t in [ps]');
       ylabel('Amplitude');
       legend('Direkter Pfad','Reflektierter Pfad','Nutzsignal');
       str = sprintf('Zeitsignal an der Stelle des Nutzers');
       title(str);
       hold off;
       print(newFilename, '-dpng','-S640,480');
    end
end
end

function [] = setGraphSettings()
caxis([0 0.3]);
xlim([0 120]);
ylim([0 100]);
set (gca, 'xticklabel', {'0',  '20', '40', '60', '80', '100',  '120'});
set (gca, 'yticklabel', {'0', '20', '40', '60',  '80', '100'});
xlabel ("x [cm]");
ylabel ("y [cm]");
grid off;
axis equal;
% still in testing mode: new color map, good for printing (white background), colorful like jet
nb_cols=256;
newmap= jet(nb_cols);
% generate smooth transition of jet() from "black to white" brightness
for (k=1:nb_cols)
    newmap(k,1:3) = (newmap(k,1:3)*(k-1) + [nb_cols-k, nb_cols-k, nb_cols-k])/(nb_cols-1);
end
jet_light_inue=newmap;
colormap( jet_light_inue );
end
