classdef Fiber
    % Fiber: Simulate optical fiber signal propagation using
    % the split-step Fourier method (SSFM).

    properties
        % Simulation sampling frequency [Hz]
        fsimu
        % Fiber length [m]
        fiber_length
        % Attenuation coefficient [dB/km]
        alpha
        % Dispersion parameter D [s/m^2]
        D
        % Dispersion slope [s/m^3]
        Dslope
        % Reference wavelength [m]
        lambda0
        % Nonlinear coefficient [1/(W·m)]
        gamma
        % Maximum allowed nonlinear phase change per step [rad]
        dphimax

        % Derived parameters:
        % Second-order dispersion coefficient β2 [s^2/m]
        b2
        % Third-order dispersion coefficient β3 [s^3/m]
        b3
        % Linear attenuation constant [1/m]
        alpha_lin
        % Frequency-domain linear operator per unit length
        linstep
    end

    methods
        function obj = Fiber(options)
            % Constructor: initialize fiber physical and simulation parameters.
            arguments
                options.fsimu            % Sampling frequency [Hz]
                options.fiber_length = 0 % Fiber length [km]
                options.alpha = 0.2      % Attenuation [dB/km]
                options.D = 17           % Dispersion D parameter [ps/(nm·km)]
                options.Dslope = 0.06    % Dispersion slope [ps/(nm^2·km)]
                options.lambda0 = 1550   % Reference wavelength [nm]
                options.gamma = 0        % Nonlinear coefficient [1/(W·km)]
                options.dphimax = 5e-3   % Max phase step [rad]
            end

            % Assign inputs to object properties, converting units to SI.
            obj.fsimu        = options.fsimu;
            obj.fiber_length = options.fiber_length * 1e3;   % km → m
            obj.alpha        = options.alpha;
            obj.D            = options.D * 1e-6;            % ps/(nm·km) → s/(m·m)
            obj.Dslope       = options.Dslope * 1e3;        % ps/(nm^2·km) → s/m^2
            obj.lambda0      = options.lambda0 * 1e-9;      % nm → m
            obj.gamma        = options.gamma;               % 1/(W·km) (assumed → 1/(W·m) internally)
            obj.dphimax      = options.dphimax;
        end

        function signalclass_out = process(obj, signalclass_in)
            % process: Apply fiber propagation to input signal class.
            %   Calls the internal SSFM routine and logs the operation.

            % Run internal split-step propagation
            signalclass = obj.process_(signalclass_in);

            % Add entry to logbook for tracking
            signalclass = signalclass.logbookentry('Fiber ');

            % Return processed signal class
            signalclass_out = signalclass;
        end

        function opt_sig = process_(obj, opt_sig)
            % process_: Internal routine for one-step fiber propagation.
            %   Computes linear and (optionally) nonlinear effects.

            % Extract time-domain field and wavelength
            signal        = opt_sig.signal;
            lambda_signal = opt_sig.lambda;

            % Update dispersion parameter D for current wavelength
            obj.D = obj.D + (lambda_signal - obj.lambda0) * obj.Dslope;

            % Compute dispersion coefficients (β2, β3)
            obj.b2 = -obj.D * lambda_signal^2 / (2 * pi * Constant.LightSpeed);
            obj.b3 = (lambda_signal^2 / (2 * pi * Constant.LightSpeed))^2 * obj.Dslope;

            % Convert attenuation from dB/km to linear 1/m
            obj.alpha_lin = obj.alpha / 10 * log(10) / 1e3;

            % Build frequency axis for FFT operations
            N = length(signal);
            faxis = linspace(-obj.fsimu/2, obj.fsimu/2, N+1);
            faxis = faxis(1:end-1);          % drop redundant endpoint
            faxis = ifftshift(faxis)';       % center zero frequency

            % Define linear operator per meter in frequency domain
            obj.linstep = -obj.alpha_lin/2 ...                 % half-step loss
                          - 1j*2*pi^2*obj.b2 .* faxis.^2 ... % second-order dispersion
                          - 1j*(4/3)*pi^3*obj.b3 .* faxis.^3; % third-order dispersion

            % Choose linear-only or nonlinear SSFM based on gamma
            if obj.gamma ~= 0
                % Full adaptive SSFM with nonlinear Schrödinger solver
                opt_out = obj.NLSE(signal);
            else
                % Single-step linear propagation in frequency domain
                opt_out = ifft( fft(signal) .* exp(obj.linstep * obj.fiber_length) );
            end

            % Update output field in signal class
            opt_sig.signal = opt_out;
        end

        function yout = NLSE(obj, xin)
            % NLSE: Solve nonlinear Schrödinger equation by adaptive split-step
            %   xin: input time-domain field
            %   Returns yout: output time-domain field after propagation

            % Initial estimate of effective length per step
            maxPow = obj.gamma * max(abs(xin).^2);
            Leff   = obj.dphimax / maxPow;
            dz     = Leff;
            z_prop = 0;

            % Apply initial half-step linear operator
            yout = fft(xin) .* exp(obj.linstep * dz/2);

            % Loop until full fiber length is reached
            while true
                % Inverse FFT to time domain for nonlinear phase shift
                yout = ifft(yout);

                % Nonlinear phase rotation per segment
                power = abs(yout).^2;
                Hnl   = exp(-1j * obj.gamma * power * dz);
                yout  = yout .* Hnl;

                % Accumulate propagation distance
                z_prop = z_prop + dz;

                % Compute adaptive step for next interval
                maxPow = obj.gamma * max(abs(yout).^2);
                dz_new = obj.dphimax / maxPow;

                % If remaining length shorter than new step, finish loop
                if z_prop + dz_new > obj.fiber_length
                    dz_new = obj.fiber_length - z_prop;
                    break;
                end

                % Half-step linear operator bridging segments
                yout = fft(yout) .* exp(obj.linstep * ((dz/2) + (dz_new/2)));
                dz   = dz_new;
            end

            % Final propagation segment: combine half-steps and nonlinear
            yout = fft(yout) .* exp(obj.linstep * ((dz/2) + (dz_new/2)));
            yout = ifft(yout);

            % Last nonlinear phase shift
            power = abs(yout).^2;
            Hnl   = exp(-1j * obj.gamma * power * dz_new);
            yout  = yout .* Hnl;

            % Final half-step linear operator to complete SSFM
            yout = fft(yout) .* exp(obj.linstep * (dz_new/2));
            yout = ifft(yout);
        end
    end
end