function [] =...
split_step_fourier_method_test(...
parameters,...
presentation_type,...
mode_3D,...
fiber_length,...
baudrate,...
impulse_width,...
impulse_spacing,...
bit_sequence_str,...
shape,...
attenuation,...
dispersion_coefficient,...
nonlinearity_coefficient,...
amp_scaling,...
width_scaling,...
static_mode,...
filename)

%% parameters


lmd0=1500*10^-9; %carrier wavelength in m
n0=1.5;
c=299792458;
eps0=8.854187817*10^-12;
mu0=1.256637061*10^-6;
k0=2*pi/lmd0/n0;
wave_impedance=mu0*c/n0;
effective_mode_area=50*(10^-6)^2;
b0=k0*n0;
b1=n0/c;


switch parameters

case 'self_defined'


	z=fiber_length*1000; %distance from signal origin [km] -> [m]
	dz=1;   %step size in m
	D=dispersion_coefficient*10^-12/10^-9/10^3;  %dispersion coefficient [ps/(nm km)] -> [s/m^2]
	n2=nonlinearity_coefficient*10^-23;%0.0024*10^-18*c*eps0;
	b2=-D*lmd0^2/2/pi/c; %group velocity dispersion parameter beta2 in s^2/m
	g=k0*n2;    %nonlinearity coefficient gamma
	a=attenuation/10000/log10(exp(1))/2;  
	impulse_width=impulse_width*10^-12; %impulse_width [ps] -> [s]
	impulse_spacing=impulse_spacing*10^-12;
	baudrate=baudrate*10^9; %baudrate [Gbd] -> [Bd]
	oversampling_factor=8;

	bit_sequence=zeros(0,0);
        bit_start=0;
        for l=1:length(bit_sequence_str)
    	   if bit_sequence_str(l)=='['
              bit_start=l;
           elseif bit_sequence_str(l)==']'
              bit_sequence=[bit_sequence,str2num(bit_sequence_str(bit_start+1:l-1))];
           end    % bit_sequence(l-1)=str2num(bit_sequence_str(l));
	end
        absorbing_boundaries=true;
	dz_variable=true;
	zeta=zeros(1,length(bit_sequence));

	if strcmp(shape,'raised_cos')
		impulse_width=impulse_width/2;
		impulse_spacing=impulse_spacing+impulse_width;
	elseif strcmp(shape,'rect')
		impulse_width=impulse_width/4;
		impulse_spacing=impulse_spacing+3*impulse_width;
	end


case 'one_soliton'

        z=100*1000; %distance from signal origin [km] -> [m]
        dz=1;   %step size in m
        D=17*10^-12/10^-9/10^3;  %dispersion coefficient [ps/(nm km)] -> [s/m^2]
        n2=5.8*10^-23;%0.0024*10^-18*c*eps0;
        b2=-D*lmd0^2/2/pi/c; %group velocity dispersion parameter beta2 in s^2/m
        g=k0*n2;    %nonlinearity coefficient gamma
        a=0;
        impulse_width=100*10^-12; %impulse_width [ps] -> [s]
        impulse_spacing=0*10^-12;
        baudrate=0; %baudrate [Gbd] -> [Bd]
        shape='sech';
        bit_sequence=[1];
        amp_scaling=1;
        width_scaling=1;
        static_mode='3D';
        oversampling_factor=16;
        absorbing_boundaries=true;
        dz_variable=true;
	zeta=zeros(1,length(bit_sequence));



case 'two_soliton'

        z=20*1000; %distance from signal origin [km] -> [m]
        dz=1;   %step size in m
        D=17*10^-12/10^-9/10^3;  %dispersion coefficient [ps/(nm km)] -> [s/m^2]
        n2=5.8*10^-23;%0.0024*10^-18*c*eps0;
        b2=-D*lmd0^2/2/pi/c; %group velocity dispersion parameter beta2 in s^2/m
        g=k0*n2;    %nonlinearity coefficient gamma
        a=0;
        impulse_width=100*10^-12; %impulse_width [ps] -> [s]
        impulse_spacing=0*10^-12;
        baudrate=0; %baudrate [Gbd] -> [Bd]
        shape='sech';
        bit_sequence=[1];
        amp_scaling=2;
        width_scaling=1;
        static_mode='3D';
        oversampling_factor=16;
        absorbing_boundaries=true;
	dz_variable=true;
	zeta=zeros(1,length(bit_sequence));


case '0_9_soliton'

        z=100*1000; %distance from signal origin [km] -> [m]
        dz=1;   %step size in m
        D=17*10^-12/10^-9/10^3;  %dispersion coefficient [ps/(nm km)] -> [s/m^2]
        n2=5.8*10^-23;%0.0024*10^-18*c*eps0;
        b2=-D*lmd0^2/2/pi/c; %group velocity dispersion parameter beta2 in s^2/m
        g=k0*n2;    %nonlinearity coefficient gamma
        a=0;
        impulse_width=50*10^-12; %impulse_width [ps] -> [s]
        impulse_spacing=200*10^-12;
        baudrate=0; %baudrate [Gbd] -> [Bd]
        shape='sech';
        bit_sequence=[1];
        amp_scaling=0.8;
        width_scaling=1;
        static_mode='3D';
        oversampling_factor=8;
        absorbing_boundaries=true;
        dz_variable=true;
	zeta=zeros(1,length(bit_sequence));


case '1_1_soliton'

        z=50*1000; %distance from signal origin [km] -> [m]
        dz=1;   %step size in m
        D=17*10^-12/10^-9/10^3;  %dispersion coefficient [ps/(nm km)] -> [s/m^2]
        n2=5.8*10^-23;%0.0024*10^-18*c*eps0;
        b2=-D*lmd0^2/2/pi/c; %group velocity dispersion parameter beta2 in s^2/m
        g=k0*n2;    %nonlinearity coefficient gamma
        a=0;
        impulse_width=50*10^-12; %impulse_width [ps] -> [s]
        impulse_spacing=100*10^-12;
        baudrate=0; %baudrate [Gbd] -> [Bd]
        shape='sech';
        bit_sequence=[1];
        amp_scaling=1.2;
        width_scaling=1;
        static_mode='3D';
        oversampling_factor=8;
        absorbing_boundaries=true;
        dz_variable=true;
	zeta=zeros(1,length(bit_sequence));



case 'breather_phi_0'

        z=100*1000; %distance from signal origin [km] -> [m]
        dz=1;   %step size in m
        D=17*10^-12/10^-9/10^3;  %dispersion coefficient [ps/(nm km)] -> [s/m^2]
        n2=5.8*10^-23;%0.0024*10^-18*c*eps0;
        b2=-D*lmd0^2/2/pi/c; %group velocity dispersion parameter beta2 in s^2/m
        g=k0*n2;    %nonlinearity coefficient gamma
        a=0;
        impulse_width=40*10^-12; %impulse_width [ps] -> [s]
        impulse_spacing=-20*10^-12;
        baudrate=0; %baudrate [Gbd] -> [Bd]
        shape='sech';
        bit_sequence=[1 1];
        amp_scaling=1;
        width_scaling=1;
        static_mode='3D';
        oversampling_factor=8;
        absorbing_boundaries=true;
	dz_variable=false;
	zeta=zeros(1,length(bit_sequence));


case 'breather_phi_pi'

        z=80*1000; %distance from signal origin [km] -> [m]
        dz=1;   %step size in m
        D=17*10^-12/10^-9/10^3;  %dispersion coefficient [ps/(nm km)] -> [s/m^2]
        n2=5.8*10^-23;%0.0024*10^-18*c*eps0;
        b2=-D*lmd0^2/2/pi/c; %group velocity dispersion parameter beta2 in s^2/m
        g=k0*n2;    %nonlinearity coefficient gamma
        a=0;
        impulse_width=40*10^-12; %impulse_width [ps] -> [s]
        impulse_spacing=-20*10^-12;
        baudrate=0; %baudrate [Gbd] -> [Bd]
        shape='sech';
        bit_sequence=[0 1 exp(j*pi) 0];
        amp_scaling=1;
        width_scaling=1;
        static_mode='3D';
        oversampling_factor=8;
        absorbing_boundaries=true;
	dz_variable=false;
	zeta=zeros(1,length(bit_sequence));



case 'kuznetsov_ma_breather'


        z=500*1000; %distance from signal origin [km] -> [m]  
        D=17*10^-12/10^-9/10^3;  %dispersion coefficient [ps/(nm km)] -> [s/m^2]
        n2=5.8*10^-23;%0.0024*10^-18*c*eps0;
        b2=-D*lmd0^2/2/pi/c; %group velocity dispersion parameter beta2 in s^2/m
        g=k0*n2;    %nonlinearity coefficient gamma
        a=0;
        impulse_width=200*10^-12; %impulse_width [ps] -> [s]
        impulse_spacing=200*10^-12;
        baudrate=0; %baudrate [Gbd] -> [Bd]
        b0=k0*n0;
        shape='breather';
        af=3;%2
        func_sc=1;%1.2;
        bit_sequence=[1];
        amp_scaling=1;
        width_scaling=1;
        static_mode='3D';
        oversampling_factor=32;
        absorbing_boundaries=false;
	dz_variable=true;
	zeta=zeros(1,length(bit_sequence));


case 'peregrine'


        z=1000*1000; %distance from signal origin [km] -> [m]
        D=17*10^-12/10^-9/10^3;  %dispersion coefficient [ps/(nm km)] -> [s/m^2]
        n2=5.8*10^-23;%0.0024*10^-18*c*eps0;
        b2=-D*lmd0^2/2/pi/c; %group velocity dispersion parameter beta2 in s^2/m
        g=k0*n2;    %nonlinearity coefficient gamma
        a=0;
        impulse_width=1800*10^-12; %impulse_width [ps] -> [s]
        impulse_spacing=0;
        baudrate=0; %baudrate [Gbd] -> [Bd]
        shape='breather';
        af=0.51;%0.6
        func_sc=1;
        bit_sequence=[1];
        amp_scaling=1;
        width_scaling=1;
        static_mode='3D';
        oversampling_factor=64;
        absorbing_boundaries=false;
	dz_variable=true;
	zeta=zeros(1,length(bit_sequence));



case 'collision'

        z=150*1000; %distance from signal origin [km] -> [m]
        D=17*10^-12/10^-9/10^3;  %dispersion coefficient [ps/(nm km)] -> [s/m^2]
        n2=5.8*10^-23;%0.0024*10^-18*c*eps0;
        b2=-D*lmd0^2/2/pi/c; %group velocity dispersion parameter beta2 in s^2/m
        g=k0*n2;    %nonlinearity coefficient gamma
        a=0;
        impulse_width=40*10^-12; %impulse_width [ps] -> [s]
        impulse_spacing=20*10^-12;
        baudrate=0; %baudrate [Gbd] -> [Bd]
        shape='sech';
        bit_sequence=[1 1 0];
        amp_scaling=1;
        width_scaling=1;
        static_mode='3D';
        oversampling_factor=16;
        absorbing_boundaries=true;
        dz_variable=true;
        zeta=[1,-0.1,0];



end

%%


switch presentation_type



    case 'static'




%% time axis


fs=oversampling_factor*200/(impulse_width*pi^2);
%t1=impulse_width/2+impulse_spacing/2;
t1=impulse_width/2+(abs(impulse_spacing/2)+impulse_spacing/2)/2;
%nb_points=ceil(t1*fs);
nb_points_imp=ceil(impulse_width/2*fs);
nb_points_spa=ceil(impulse_spacing/2*fs);
%nb_points=nb_points_imp+nb_points_spa;
nb_points=nb_points_imp+(abs(nb_points_spa)+nb_points_spa)/2;
tx_signal_time =linspace(-t1,t1,2*nb_points-1);% linspace(0,t1,nb_points);
%tx_signal_time=[-tx_signal_time(end:-1:2),tx_signal_time];
delta_t=tx_signal_time(2)-tx_signal_time(1);



%%

signal_data_z0_array=zeros(0,0);

%% create tx data


switch shape

    case 'sech'

    for m=1:length(bit_sequence)
       
            if nb_points_spa<0 && m>1
                 signal_data_z0_array=...
                 [signal_data_z0_array(1:end+nb_points_spa*2),...
                 signal_data_z0_array(end+nb_points_spa*2+1:end)+abs(bit_sequence(m))*15*sqrt(abs(b2))/g^(1/2)/impulse_width*sech((1/width_scaling)*15/impulse_width*tx_signal_time(1:-nb_points_spa*2)).*exp(-1i*2*g^(1/6)*zeta(m)/sqrt(abs(b2))*tx_signal_time(1:-nb_points_spa*2))*exp(1i*angle(bit_sequence(m))),...
                 abs(bit_sequence(m))*15*sqrt(abs(b2))/g^(1/2)/impulse_width*sech((1/width_scaling)*15/impulse_width*tx_signal_time(-nb_points_spa*2+1:end)).*exp(-1i*2*g^(1/6)*zeta(m)/sqrt(abs(b2))*tx_signal_time(-nb_points_spa*2+1:end))*exp(1i*angle(bit_sequence(m)))];
                 signal_data_z0_array=[signal_data_z0_array(1:end-((2*nb_points-1)+nb_points_spa)),signal_data_z0_array(end-((2*nb_points-1)+nb_points_spa)+2:end)];
            else
                 signal_data_z0_array=...
                 [signal_data_z0_array(1:end-1),...
                 abs(bit_sequence(m))*15*sqrt(abs(b2))/g^(1/2)/impulse_width*sech((1/width_scaling)*15/impulse_width*tx_signal_time).*exp(-1i*2*g^(1/6)*zeta(m)/sqrt(abs(b2))*tx_signal_time)*exp(1i*angle(bit_sequence(m)))];
            end
        
    end


    case 'gauss'

    for m=1:length(bit_sequence)
        if nb_points_spa<0 && m>1
            signal_data_z0_array=...
            [signal_data_z0_array(1:end+nb_points_spa*2),...
            signal_data_z0_array(end+nb_points_spa*2+1:end)+abs(bit_sequence(m))*15*sqrt(abs(b2))/g^(1/2)/impulse_width*exp(-0.3*(15/impulse_width)^2*tx_signal_time(1:-nb_points_spa*2).^2)*exp(1i*angle(bit_sequence(m))),...
            abs(bit_sequence(m))*15*sqrt(abs(b2))/g^(1/2)/impulse_width*...
            exp(-0.3*((1/width_scaling)*15/impulse_width)^2*tx_signal_time(-nb_points_spa*2+1:end).^2)*exp(1i*angle(bit_sequence(m)))];
            signal_data_z0_array=[signal_data_z0_array(1:end-((2*nb_points-1)+nb_points_spa)),signal_data_z0_array(end-((2*nb_points-1)+nb_points_spa)+2:end)];
        else
            signal_data_z0_array=...
            [signal_data_z0_array(1:end-1),...
            abs(bit_sequence(m))*15*sqrt(abs(b2))/g^(1/2)/impulse_width*exp(-0.3*(15/impulse_width)^2*tx_signal_time.^2)*exp(1i*angle(bit_sequence(m)))];
        end
        
    end


    case 'rect'

    for m=1:length(bit_sequence)
        if nb_points_spa<0 && m>1
            win=zeros(1,length(tx_signal_time));
            win(nb_points-(nb_points_imp-1):nb_points+(nb_points_imp-1))=...
                15*sqrt(abs(b2))/g^(1/2)/(4*impulse_width);
            signal_data_z0_array=...
            [signal_data_z0_array(1:end+nb_points_spa*2),...
            signal_data_z0_array(end+nb_points_spa*2+1:end)+abs(bit_sequence(m))*win(1:-nb_points_spa*2)*exp(1i*angle(bit_sequence(m))),...
            abs(bit_sequence(m))*win(-nb_points_spa*2+1:end)*exp(1i*angle(bit_sequence(m)))];
            signal_data_z0_array=[signal_data_z0_array(1:end-((2*nb_points-1)+nb_points_spa)),signal_data_z0_array(end-((2*nb_points-1)+nb_points_spa)+2:end)];
        else
            win=zeros(1,length(tx_signal_time));
            win(nb_points-(nb_points_imp-1):nb_points+(nb_points_imp-1))=...
                15*sqrt(abs(b2))/g^(1/2)/(4*impulse_width);            
            signal_data_z0_array=...
            [signal_data_z0_array(1:end-1),abs(bit_sequence(m))*win*exp(1i*angle(bit_sequence(m)))];
        end
    end


    case 'raised_cos'

    for m=1:length(bit_sequence)
        if nb_points_spa<0 && m>1
            win=zeros(1,length(tx_signal_time));
            win(nb_points-(nb_points_imp-1):nb_points+(nb_points_imp-1))=...
                15*sqrt(abs(b2))/g^(1/2)/(2*impulse_width)*...
                cos((pi/2)/(impulse_width/2).*...
                tx_signal_time(nb_points-(nb_points_imp-1):nb_points+(nb_points_imp-1))).^2;
            signal_data_z0_array=...
            [signal_data_z0_array(1:end+nb_points_spa*2),...
            signal_data_z0_array(end+nb_points_spa*2+1:end)+abs(bit_sequence(m))*win(1:-nb_points_spa*2)*exp(1i*angle(bit_sequence(m))),...
            abs(bit_sequence(m))*win(-nb_points_spa*2+1:end)*exp(1i*angle(bit_sequence(m)))];
            signal_data_z0_array=[signal_data_z0_array(1:end-((2*nb_points-1)+nb_points_spa)),signal_data_z0_array(end-((2*nb_points-1)+nb_points_spa)+2:end)];
       else
            win=zeros(1,length(tx_signal_time));
            win(nb_points-(nb_points_imp-1):nb_points+(nb_points_imp-1))=...
                15*sqrt(abs(b2))/g^(1/2)/(2*impulse_width)*...
                cos((pi/2)/(impulse_width/2).*...
                tx_signal_time(nb_points-(nb_points_imp-1):nb_points+(nb_points_imp-1))).^2;
            signal_data_z0_array=...
            [signal_data_z0_array(1:end-1),abs(bit_sequence(m))*win*exp(1i*angle(bit_sequence(m)))];
        end
    end

    
    case 'breather'

    bf=sqrt(8*af*(1-2*af));
    z0=-500000; 

    %func_sc=(sqrt(abs(b2))/((g^(1/6))*sqrt(2*af-1)*impulse_width))*...
            %acosh((4*af-1)/sqrt(2*af));
    for m=1:length(bit_sequence)
        if bit_sequence(m)~=0
           signal_data_z0_array=...
           [signal_data_z0_array(1:end-1),...
           abs(bit_sequence(m))/(g^(1/3))*func_sc*exp(-1i*g^(1/3)*z0)*...
           (ones(1,length(tx_signal_time))+...
           (((2-4*af)*cosh(bf*g^(1/3)*z0)+1i*bf*sinh(-bf*g^(1/3)*z0))./((sqrt(2*af)*cos(func_sc*2*sqrt(1-2*af)*(g^(1/6))/sqrt(abs(b2)).*tx_signal_time))-...
           cosh(bf*g^(1/3)*z0)*ones(1,length(tx_signal_time)))))*exp(1i*angle(bit_sequence(m)))];
        else
           signal_data_z0_array=...
           [signal_data_z0_array(1:end-1),0.*tx_signal_time];
        end
    end


end
%%



signal_data_z0=amp_scaling*signal_data_z0_array;


%tx_signal_time=linspace(tx_signal_time(1),...
%                        t1+(length(bit_sequence)-1)*2*tx_signal_time(end-1),...
%                        2*length(bit_sequence)*nb_points-length(bit_sequence)-(length(bit_sequence)-1));

tx_signal_time=-t1:delta_t:-t1+(length(signal_data_z0_array)-1)*delta_t;



tx_signal_data=signal_data_z0;
rx_signal_time=tx_signal_time;


peak_power=(max(abs(signal_data_z0))^2*effective_mode_area)/(2*wave_impedance);
%mean_power=effective_mode_area/(2*wave_impedance*(impulse_width+impulse_spacing))*delta_t*sum(signal_data_z0.^2)/length(bit_sequence);
mean_power=effective_mode_area/(2*wave_impedance*(length(signal_data_z0)*delta_t))*delta_t*sum(abs(signal_data_z0).^2);
%impulse_energy=mean_power*(impulse_width+impulse_spacing);
impulse_energy=sum(abs(signal_data_z0).^2)*delta_t*effective_mode_area/(2*wave_impedance);





f=linspace(-fs/2,fs/2,length(tx_signal_data));
fc=c/lmd0;



if length(tx_signal_time)<3000
	data_memory=zeros(1000,1+length(tx_signal_time)+1);
        data_memory(1,2:end-1)=signal_data_z0;
        data_memory(1,end)=impulse_energy;
else 
	data_memory=zeros(1000,3002);
        data_memory(1,2:end-1)=signal_data_z0(floor(linspace(1,length(signal_data_z0),3000)));
        data_memory(1,end)=impulse_energy;
end
%data_memory=zeros(10*z/1000,length(tx_signal_time));




%% iterative calculation of pulse propagation


dz_sum=0;
t_elapsed=0;


absorption_function=exp(-0.02*(sech(2*15/(impulse_width+impulse_spacing)*(tx_signal_time-tx_signal_time(end))).^2+sech(2*15/(impulse_width+impulse_spacing)*(tx_signal_time-tx_signal_time(1))).^2));


%absorption_function=exp(-0.02*(sech(2*((1/width_scaling)*15/impulse_width)*(tx_signal_time-tx_signal_time(end))).^2+sech(2*((1/width_scaling)*15/impulse_width)*(tx_signal_time-tx_signal_time(1))).^2));

limit_exceeded=false;

while dz_sum<z

tic;

    if dz_variable

    	dz=max(z*t_elapsed/240,min(pi/(500*g*max(abs(signal_data_z0))^2),z/1000));
        
        if dz==z*t_elapsed/240
           limit_exceeded=true;
        end

    end

    H=exp(-1i*b2*4*pi^2*f.^2*dz/2);


    rx_signal_data_kerr=signal_data_z0.*exp(-1i*g*abs(signal_data_z0).^2*dz); % consider kerr effect

    rx_signal_data_kerr_att=rx_signal_data_kerr.*exp(-a*dz); % consider attenuation
    rx_signal_data_kerr_att_freq=fftshift(fft(rx_signal_data_kerr_att));

    rx_signal_data_kerr_att_disp_freq=rx_signal_data_kerr_att_freq.*H; % consider chromatic dispersion
    rx_signal_data_kerr_att_disp=ifft(ifftshift(rx_signal_data_kerr_att_disp_freq));

    signal_data_z0=rx_signal_data_kerr_att_disp;

    
    if absorbing_boundaries

	signal_data_z0=signal_data_z0.*(absorption_function.^dz);


    end

    impulse_energy=sum(abs(signal_data_z0).^2)*delta_t*effective_mode_area/(2*wave_impedance);

    dz_sum=dz_sum+dz;


	if length(tx_signal_time)<3000
		 data_memory(ceil(dz_sum/z*1000),:)=[dz_sum,signal_data_z0,impulse_energy];
	else
		data_memory(ceil(dz_sum/z*1000),:)=[dz_sum,signal_data_z0(floor(linspace(1,length(signal_data_z0),3000))),impulse_energy];
	end


t_elapsed=toc;

   
end


rx_signal_data=signal_data_z0;
rx_signal_data_hf=rx_signal_data.*exp(1i*(2*pi*c/lmd0*(tx_signal_time+b1*z)-b0*z));


switch static_mode

   case '1D'

   %% static plot for z=fiber_length

     g=figure(1);
     set(g,'visible','off');
     set(0, 'defaultaxesfontsize', 8);
     set(0, 'defaultaxesfontname', "Helvetica");
     set(0, 'defaulttextfontsize', 8);
     set(0, 'defaulttextfontname', "Helvetica");
     exp_max=floor(max(max(log10(abs(data_memory(:,2:end-1))))));    
     plot(tx_signal_time*10^12,tx_signal_data/(10^exp_max),'Color','blue','LineWidth',2);
     hold on;
     %plot(rx_signal_time*10^12,abs(rx_signal_data)/(10^exp_max),'Color','red');
     axis([min(rx_signal_time*10^12)-0.5*impulse_width/2*10^12,...
         max(rx_signal_time*10^12)+0.5*impulse_width/2*10^12,...
         0,...
         1.1*max(max(abs(tx_signal_data)/(10^exp_max)),max(abs(rx_signal_data)/(10^exp_max)))]);
     xlabel('time \tau [ps]');
     ylabel(strcat('|A| / ','10^',num2str(exp_max),' [V/m]'));
     %legend('Tx signal','Rx signal');
     title(['P_{max,in} = ',num2str(peak_power*1000),' mW   -   ',...
           'P_{avg} = ',num2str(mean_power*1000),' mW   -   ',...
           'W = ', num2str(impulse_energy*10^12),' pJ']);

print(filename, '-dpng',"-S930,350");



   case '1D_dynamic'

     g=figure(1);
     set(g,'visible','off');
     set(0, 'defaultaxesfontsize', 8);
     set(0, 'defaultaxesfontname', "Helvetica");
     set(0, 'defaulttextfontsize', 8);
     set(0, 'defaulttextfontname', "Helvetica");
     exp_max=floor(max(max(log10(abs(data_memory(:,2:end-1))))));
     plot(tx_signal_time*10^12,abs(tx_signal_data)/(10^exp_max),'Color','blue','LineWidth',2);
     axis([min(tx_signal_time*10^12),max(tx_signal_time*10^12),0,max(max([abs(data_memory(:,2:end-1))/(10^exp_max);tx_signal_data/(10^exp_max)]))]);
     xlabel('time \tau [ps]');
     ylabel(strcat('|A| / ','10^',num2str(exp_max),' [V/m]'));
     hold on;


     for l=1:20:length(data_memory(:,1))

          h0=plot(rx_signal_time(floor(linspace(1,length(rx_signal_time),length(data_memory(1,1:end-1))-1)))*10^12,abs(data_memory(l,2:end-1))/(10^exp_max),'Color','red','LineWidth',2);
          title(['z=',num2str(data_memory(l,1)/1000),'km','     ',...
                 'E_0=',num2str(data_memory(1,end)*10^12),'pJ','     ',...
                 'E=',num2str(data_memory(l,end)*10^12),'pJ']);
          %legend('Tx signal','Rx signal');


          new_filename=[filename '_'  num2str((l-1)/20)];
          print(new_filename, '-dpng',"-S640,480");

          delete(h0);


     end


   case '3D'


switch mode_3D

     case 'surface'

     g=figure(1);
     set(g,'visible','off');
     set(0, 'defaultaxesfontsize', 8);
     set(0, 'defaultaxesfontname', "Helvetica");
     set(0, 'defaulttextfontsize', 8);
     set(0, 'defaulttextfontname', "Helvetica");
     exp_max=floor(max(max(log10(abs(data_memory(:,2:end))))));
     surf(rx_signal_time(floor(linspace(1,length(rx_signal_time),length(data_memory(1,:))-1)))*10^12,...
         data_memory(:,1)/1000,...                                                              %(1:length(data_memory(:,1)))*z/1000/length(data_memory(:,1)),
         abs(data_memory(:,2:end))/(10^exp_max),'edgecolor','none');
     set(gca,'YDir','rev');
     xlabel('time \tau [ps]');
     ylabel('distance z [km]');
     zlabel(strcat('|A| / ','10^',num2str(exp_max),' [V/m]'));
     bit_sequence_str='';
     for k=1:length(bit_sequence)
        bit_sequence_str=[bit_sequence_str,'['];
        if isreal(bit_sequence(k))
           bit_sequence_str=[bit_sequence_str,num2str(bit_sequence(k)),']'];
        else
           if abs(bit_sequence(k))~=1 
              bit_sequence_str=[bit_sequence_str,num2str(abs(bit_sequence(k))),'*'];
           end
           if angle(bit_sequence(k))/pi~=1
              bit_sequence_str=[bit_sequence_str,' exp(j*',num2str(angle(bit_sequence(k))/pi),'*pi)]'];
           else
              bit_sequence_str=[bit_sequence_str,' exp(j*pi)]'];
           end           
        end
     end

     if (~strcmp(parameters,'self_defined'))
     	title_str=['T=',num2str(impulse_width*10^12),' ps,  ',...
               '\Delta T=',num2str(impulse_spacing*10^12),' ps,   ',...
               'pulse sequence=',bit_sequence_str,',   ',...
               'pulse shape = ',shape;...	       
               'L=',num2str(z/1000),' km,   ',...
               '\alpha=',num2str(attenuation),' dB/km,   ',...
               'D=',num2str(dispersion_coefficient),' ps/nm/km,   ',...
               'A*/A=',num2str(amp_scaling),',   ',...
               'T*/T =',num2str(width_scaling)];
     
     else
     	title_str='';
     end
     axis_vals=axis;
     title(title_str);
     text(axis_vals(1),axis_vals(3),-0.7*axis_vals(6),['P_{max,in} = ',num2str(peak_power*1000),' mW   -   ',...
           'P_{avg} = ',num2str(mean_power*1000),' mW   -   ',...
           'W = ', num2str(impulse_energy*10^12),' pJ']);
    
    
     case 'heatmap'

     g=figure(1);
     set(g,'visible','off');
     set(0, 'defaultaxesfontsize', 8);
     set(0, 'defaultaxesfontname', "Helvetica");
     set(0, 'defaulttextfontsize', 8);
     set(0, 'defaulttextfontname', "Helvetica");
     imagesc(rx_signal_time(floor(linspace(1,length(rx_signal_time),length(data_memory(1,:))-1)))*10^12,...
             data_memory(:,1)/1000,...                                                              %(1:length(data_memory(:,1)))*z/1000/length(data_memory(:,1)),...
             abs(data_memory(:,2:end)));
     xlabel('time \tau [ps]');
     ylabel('distance z [km]');
     title(['P_{max,in} = ',num2str(peak_power*1000),' mW   -   ',...
           'P_{avg} = ',num2str(mean_power*1000),' mW   -   ',...
           'W = ', num2str(impulse_energy*10^12),' pJ']);



end;

print(filename, '-dpng',"-S640,480");


end



     case 'chart'


	peak_power=linspace(10^-6,10^2,100000);

	width=15*sqrt(abs(b2))/g^(1/2)./sqrt(peak_power*2*wave_impedance/effective_mode_area);

	g=figure(1);
        set(g,'visible','off');
        set(0, 'defaultaxesfontsize', 8);
        set(0, 'defaultaxesfontname', "Helvetica");
        set(0, 'defaulttextfontsize', 8);
        set(0, 'defaulttextfontname', "Helvetica");
	loglog(peak_power*10^3,width*10^12,'LineWidth',2);
	ylabel('pulse duration T [ps]');
	xlabel('peak power P_{max} [mW]');
        grid on; 


print(filename, '-dpng',"-S930,350");



end

end
