%function [rOpt,state] = lin_step(state,aOpt,aStepSize)

function [rOpt_x,rOpt_y,z_prop,lin_z_test,...
    corr_length,n_plates_done,missing_dz,n_step,test_plates,...
    test_plate_numbers,brf,common_beta_x,common_beta_y,alpha_lin_x,alpha_lin_y]...
    = lin_step(...
    opt_x,opt_y,z_prop,lin_z_test,aStepSize,corr_length,...
    n_plates_done,missing_dz,n_step,test_plates,test_plate_numbers,...
    brf,common_beta_x,common_beta_y,alpha_lin_x,alpha_lin_y)

%%%%%% 1) Update and Check Distances etc. %%%%%%

% update propgated distance z_prop
z_prop = z_prop + aStepSize;

% calculate the number of plates needed for the so far propagated fiber length
n_plates = ceil(z_prop/corr_length);

% subtract the number of plates which were already processed
n_plates_left = n_plates - n_plates_done;

% compute last plate size ( if it fits, it should be 0)
if missing_dz > aStepSize
    
    last_plate = aStepSize;
    missing_dz = missing_dz-aStepSize;
    plate_sizes = last_plate;
    plate_numbers = n_plates;
    
else
    
    last_plate = aStepSize - missing_dz - (n_plates_left-1)*corr_length;
    
    if missing_dz == 0
        missing_dz = [];
    end
    
    %build vector of plate lengths with missing plate part from prev.
    %iterartion , then some normal plates and finally a fraction of a plate
    %to fit into the step length
    plate_sizes = [missing_dz corr_length*ones(1,n_plates_left-1) last_plate];
    
    if n_plates_done == 0
        plate_numbers =[(n_plates_done+1):(n_plates-1) n_plates];
    else
        plate_numbers = [n_plates_done (n_plates_done+1):(n_plates-1) n_plates];  % not wrking yet
    end
    
    %remember for next step 
    missing_dz = corr_length - last_plate;
    
end

plate_steps = repmat(n_step,1,length(plate_sizes));

%
%figure;stem(plate_sizes);
test_plates = [test_plates,plate_sizes];

test_plate_numbers = [test_plate_numbers, plate_numbers];

%%%%%% 2) Apply Waveplate Model %%%%%%

% transfer optical envelope to frequency domain for effective convolution with transfer function h
opt_x=fft(opt_x);
opt_y=fft(opt_y);

% db1 = gpuArray(brf.db1);
% db0 = gpuArray(brf.db0);
% common_beta_x = gpuArray(common_beta_x);
% common_beta_y = gpuArray(common_beta_y);

db1 = (brf.db1);
db0 = (brf.db0);
common_beta_x = (common_beta_x);
common_beta_y = (common_beta_y);

% process every waveplate with given sizes in plate_sizes
for n=1:length(plate_sizes)
    dz = plate_sizes(n);
    
    %         figure(87);subplot(2,1,1);plot(real(x(900:1150)));subplot(2,1,2);plot(real(y(900:1150)));
    %         MOV1=[MOV1 getframe(87)];
    
    % extract rotation matrix from pre calculated matrices
    matR = brf.matR{plate_numbers(n)};
    
    % transform to eigenvalue of of fiber segment
    tOpt.X = conj(matR(1,1))*opt_x + conj(matR(2,1))*opt_y;
    tOpt.Y = conj(matR(1,2))*opt_x + conj(matR(2,2))*opt_y;
    
    % calculate statistical delta beta for pmd
    delta_beta = 0.5*(db1+db0(n))/corr_length;
    % build transfer function with delta beta
    %common.beta = beta1+beta2*omega^2
    
    %accumulate delta beta for log...
    brf.simdgd = brf.simdgd + (db1(length(db1)/2+1)+db0(n))/corr_length;
    
    h.X = exp(-1j*(common_beta_x-delta_beta)*dz);
    h.Y = exp(-1j*(common_beta_y+delta_beta)*dz);
    % delta_beta has to be added to the transfer function
    
    % process with transfer function
    tOpt.X = h.X.*tOpt.X ;
    tOpt.Y = h.Y.*tOpt.Y ;
    
    % rotate back
    opt_x = matR(1,1)*tOpt.X + matR(1,2)*tOpt.Y;
    opt_y = matR(2,1)*tOpt.X + matR(2,2)*tOpt.Y;
    
end

lin_z_test = lin_z_test + sum(plate_sizes,2);

%update the number of processed plates so far
n_plates_done = n_plates_done + n_plates_left;

% attanuate the signal each linear state with alpha
% ( 0.2dB = 4.6052e-05 )
rOpt_x=ifft(exp(-alpha_lin_x*aStepSize/2).*opt_x); % /2 not sure why (have to find it in formulas) 
rOpt_y=ifft(exp(-alpha_lin_y*aStepSize/2).*opt_y); % but not relevant for now

end