function [] = SPPplus_numerical_v1(time_wave,resp_wave,L,k,num_mode,...
    print_state,filename)
% This Function applies the SPP Analysis by means of a numerical
% differentiation. The results are then printed/displayed/saved by other
% functions as requested
    %Inputs: time_wave = Lx1 vector of time at each measurement point
            %resp_wave = Lx3 matrix of the strain, rate and stress data,
                %with each row representing a measuring point
            %L = number of measurement points in the extracted data
            %k = step size for numerical differentiation
            %num_mode = method# to use for numerical differentiation,
                %specified by the user
            %print_state = vector identifying which optional outputs have
                %been requested by the user
            %filename = name of file from which the data originated

%Calculate the frequency of the response
omega = max(abs(resp_wave(:,2)))/max(abs(resp_wave(:,1)));

save_resp_wave = resp_wave;
resp_wave(:,2) = resp_wave(:,2)/omega;

%Calculate the average timestep size
dt = sum(time_wave(2:L)-time_wave(1:(L-1)))/(L-1);

%Initialize result vectors
N_vect = zeros(L,3);
B_vect = zeros(L,3);
rd = zeros(L,3);
mag_rd = zeros(L,1);
rdd = zeros(L,3);
rddd = zeros(L,3);
rd_cross_rdd = zeros(L,3);
mag_rd_cross_rdd = zeros(L,1);
rd_cross_rd_cross_rdd = zeros(L,3);
Gp_t = zeros(L,1);
Gpp_t = zeros(L,1);
Gp_t_dot = zeros(L,1);
Gpp_t_dot = zeros(L,1);

%Perform Numerical differentiation
for p=1:L
    if num_mode == 0
        %use "standard" differentiation on the response data (does not make
            %assumtions abot the form of the data, uses forward/backward
            %difference at ends and centered derivative elsewhere)
        if p<=2*k
            rd(p,:) = (-resp_wave(p+2*k,:)+4*resp_wave(p+k,:)-3*resp_wave(p,:))/(2*k*dt);
            rdd(p,:) = (-resp_wave(p+3*k,:)+4*resp_wave(p+2*k,:)-5*resp_wave(p+k,:)+2*resp_wave(p,:))/((k*dt)^2);
            rddd(p,:) = (-3*resp_wave(p+4*k,:)+14*resp_wave(p+3*k,:)-24*resp_wave(p+2*k,:)+18*resp_wave(p+k,:)-5*resp_wave(p,:))/(2*(k*dt)^3);
        elseif p>=(L-2*k)
            rd(p,:) = (resp_wave(p-2*k,:)-4*resp_wave(p-k,:)+3*resp_wave(p,:))/(2*k*dt);
            rdd(p,:) = (-resp_wave(p-3*k,:)+4*resp_wave(p-2*k,:)-5*resp_wave(p-k,:)+2*resp_wave(p,:))/((k*dt)^2);
            rddd(p,:) = (3*resp_wave(p-4*k,:)-14*resp_wave(p-3*k,:)+24*resp_wave(p-2*k,:)-18*resp_wave(p-k,:)+5*resp_wave(p,:))/(2*(k*dt)^3);
        else
            rd(p,:) = (-resp_wave(p+2*k,:)+8*resp_wave(p+k,:)-8*resp_wave(p-k,:)+resp_wave(p-2*k,:))/(12*k*dt);
            rdd(p,:) = (-resp_wave(p+2*k,:)+16*resp_wave(p+k,:)-30*resp_wave(p,:)+16*resp_wave(p-k,:)-resp_wave(p-2*k,:))/(12*(k*dt)^2);
            rddd(p,:) = (resp_wave(p+2*k,:)-2*resp_wave(p+k,:)+2*resp_wave(p-k,:)-resp_wave(p-2*k,:))/(2*(k*dt)^3);
        end
    end
    
    if num_mode == 1
        % use "looped" differentiation on the response data (assumes steady
            %state oscilatory, uses centered difference everywhere by
            %connecting the data in a loop)
        p2 = p+2*k;
        p1 = p+k;
        pn1 = p-k;
        pn2 = p-2*k;
        if p2>L
            p2 = p2-L;
        end
        if p1>L
            p1 = p1-L;
        end
        if pn1<1
            pn1 = pn1+L;
        end
        if pn2<1
            pn2 = pn2+L;
        end
        rd(p,:) = (-resp_wave(p2,:)+8*resp_wave(p1,:)-8*resp_wave(pn1,:)+resp_wave(pn2,:))/(12*k*dt);
        rdd(p,:) = (-resp_wave(p2,:)+16*resp_wave(p1,:)-30*resp_wave(p,:)+16*resp_wave(pn1,:)-resp_wave(pn2,:))/(12*(k*dt)^2);
        rddd(p,:) = (resp_wave(p2,:)-2*resp_wave(p1,:)+2*resp_wave(pn1,:)-resp_wave(pn2,:))/(2*(k*dt)^3);
    end
    
    %Calculate SPP metrics by linear algebra using derivatives
    rd_cross_rdd(p,:) = [rd(p,2)*rdd(p,3)-rd(p,3)*rdd(p,2);...
        rd(p,3)*rdd(p,1)-rd(p,1)*rdd(p,3);rd(p,1)*rdd(p,2)-rd(p,2)*rdd(p,1)];
    rd_cross_rd_cross_rdd(p,:) = [rd(p,2)*rd_cross_rdd(p,3)-rd(p,3)*rd_cross_rdd(p,2);...
        rd(p,3)*rd_cross_rdd(p,1)-rd(p,1)*rd_cross_rdd(p,3);rd(p,1)*rd_cross_rdd(p,2)-rd(p,2)*rd_cross_rdd(p,1)];
    mag_rd(p) = sqrt(rd(p,1)^2+rd(p,2)^2+rd(p,3)^2);
    mag_rd_cross_rdd(p) = sqrt(rd_cross_rdd(p,1)^2+rd_cross_rdd(p,2)^2+rd_cross_rdd(p,3)^2);
    N_vect(p,:) = -rd_cross_rd_cross_rdd(p,:)/(mag_rd(p)*mag_rd_cross_rdd(p));
    B_vect(p,:) = rd_cross_rdd(p,:)/mag_rd_cross_rdd(p);
    Gp_t(p) = -rd_cross_rdd(p,1)/rd_cross_rdd(p,3);
    Gpp_t(p) = -rd_cross_rdd(p,2)/rd_cross_rdd(p,3);
    Gp_t_dot(p) = -rd(p,2)*(rddd(p,1)*rd_cross_rdd(p,1)+rddd(p,2)*rd_cross_rdd(p,2)+rddd(p,3)*rd_cross_rdd(p,3))/(rd_cross_rdd(p,3))^2;
    Gpp_t_dot(p) = rd(p,1)*(rddd(p,1)*rd_cross_rdd(p,1)+rddd(p,2)*rd_cross_rdd(p,2)+rddd(p,3)*rd_cross_rdd(p,3))/(rd_cross_rdd(p,3))^2;
end

%Determine the speed of the moduli
G_speed = sqrt((Gp_t_dot).^2+(Gpp_t_dot).^2);

%=========================Print SPP Results================================

SPPplus_print_v1(filename,k,omega,save_resp_wave,Gp_t,Gpp_t,G_speed,Gp_t_dot,...
    Gpp_t_dot,print_state(1),rd,N_vect,B_vect,'NUMERICAL',...
    'Data calculated via Numerical Differentiation','Step size of ')

%=========================Display SPP Results==============================

getfigs = [1,0];%specifies to produce Standard SPP figure
            %(ex. SPP_LAOS_NUMERICAL_PLOT_sample_file1.jpg)

SPPplus_figures_v1(resp_wave,0,Gp_t,Gpp_t,G_speed,Gp_t_dot,Gpp_t_dot,...
    time_wave,0,0,0,0,print_state(2),'NUMERICAL',filename,getfigs)

end

