function [] = SPPplus_fourier_v1(time_wave,resp_wave,L,omega,M,p,...
    print_state,filename)
% This Function applies the SPP Analysis by means of a fourier series. 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
            %omega = frequency of oscilation (rad/s)
            %M = number of harmonics for stress
            %p = number of cycles
            %print_state = vector identifying which optional outputs have
                %been requested by the user
            %filename = name of file from which the data originated

%=============Calculate SPP Metrics via Fourier Series=====================

format long G;

%Determine maximum number of calculatable harmonics
W = round(L/(2*p));

%Perform fourier transform of data
fft_resp = fft(resp_wave,L);
sign_fft_resp = fft_resp/L;
f_domain = linspace(0,W,W+1);
FWave = fft_resp(1:L/2,:);

%Determine indeces of relavent harmonics
k = zeros(1,W+1);
for n=1:W+1
    k(n)=p*(n-1)+1;
end

%Normalize size of fourier harmonics
sign_fft_resp_s = sign_fft_resp(:,3);
ft_amp_sign_fft_resp_s=2*abs(sign_fft_resp_s(k));
ft_amp_sign_fft_resp_s=ft_amp_sign_fft_resp_s./ft_amp_sign_fft_resp_s(2);

%separate real and imaginary components of the fourier response
FWreal = real(FWave);
FWimag = imag(FWave);

%Calculate Fourier harmonics of strain, rate, and stress
An1_n=2*FWreal(:,1)/L;
An1_n(1)=0; %becasue stress has no offset term.       %F1(1)=F1(1)/2;
Bn1_n=-2*FWimag(:,1)/L;
An1_r=2*FWreal(:,2)/L;
An1_r(1)=0; %becasue stress has no offset term.       %F1(1)=F1(1)/2;
Bn1_r=-2*FWimag(:,2)/L;
An1=2*FWreal(:,3)/L;
An1(1)=0; %becasue stress has no offset term.       %F1(1)=F1(1)/2;
Bn1=-2*FWimag(:,3)/L;

%Calculate offset factor
Delta=atan(An1_n(p+1)/Bn1_n(p+1));

%Compensate for ofset in harmonics
J=length(An1_n);
for nn=1:J
    An_n(nn)=An1_n(nn)*cos(Delta/(p))-Bn1_n(nn)*sin(Delta/(p));
    Bn_n(nn)=Bn1_n(nn)*cos(Delta/(p))+An1_n(nn)*sin(Delta/(p));
end
J=length(An1_r);
for nn=1:J
    An_r(nn)=An1_r(nn)*cos(Delta/(p))-Bn1_r(nn)*sin(Delta/(p));
    Bn_r(nn)=Bn1_r(nn)*cos(Delta/(p))+An1_r(nn)*sin(Delta/(p));
end
J=length(An1);
for nn=1:J
    An(nn)=An1(nn)*cos(Delta/(p))-Bn1(nn)*sin(Delta/(p));
    Bn(nn)=Bn1(nn)*cos(Delta/(p))+An1(nn)*sin(Delta/(p));
end

%Set Result Vectors
Recon_Wave = [ones(L,1)*An_n(1),ones(L,1)*An_r(1),ones(L,1)*An(1)];
Recon_Wave_dot = zeros(L,3);
Recon_Wave_ddot = 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);

%Find fourier series for strain
for n=1:2:1
    Recon_Wave(:,1) = Recon_Wave(:,1) + An_n(p*n+1)*cos(n*omega*time_wave)+ Bn_n(p*n+1)*sin(n*omega*time_wave);
    Recon_Wave_dot(:,1) = Recon_Wave_dot(:,1) - n*omega*An_n(p*n+1)*sin(n*omega*time_wave) + n*omega*Bn_n(p*n+1)*cos(n*omega*time_wave);
    Recon_Wave_ddot(:,1) = Recon_Wave_ddot(:,1) - n^2*omega^2*An_n(p*n+1)*cos(n*omega*time_wave) - n^2*omega^2*Bn_n(p*n+1)*sin(n*omega*time_wave);
end

%Find fourier series for rate
for n=1:2:1
    Recon_Wave(:,2) = Recon_Wave(:,2) + An_r(p*n+1)*cos(n*omega*time_wave)+ Bn_r(p*n+1)*sin(n*omega*time_wave);
    Recon_Wave_dot(:,2) = Recon_Wave_dot(:,2) - n*omega*An_r(p*n+1)*sin(n*omega*time_wave) + n*omega*Bn_r(p*n+1)*cos(n*omega*time_wave);
    Recon_Wave_ddot(:,2) = Recon_Wave_ddot(:,2) - n^2*omega^2*An_r(p*n+1)*cos(n*omega*time_wave) - n^2*omega^2*Bn_r(p*n+1)*sin(n*omega*time_wave);
end

%Find harmonics for Gp_t and Gpp_t
amp_strain=max(abs(Recon_Wave(:,1)));
if (resp_wave(3,3)<0)
    Gp_n=-Bn/amp_strain;
    Gpp_n=-An/amp_strain;
else
    Gp_n=Bn/amp_strain;
    Gpp_n=An/amp_strain;
end

%Find fourier series for results
for n=1:2:M
    %Stress
    Recon_Wave(:,3) = Recon_Wave(:,3) + An(p*n+1)*cos(n*omega*time_wave)+ Bn(p*n+1)*sin(n*omega*time_wave);
    Recon_Wave_dot(:,3) = Recon_Wave_dot(:,3) - n*omega*An(p*n+1)*sin(n*omega*time_wave) + n*omega*Bn(p*n+1)*cos(n*omega*time_wave);
    Recon_Wave_ddot(:,3) = Recon_Wave_ddot(:,3) - n^2*omega^2*An(p*n+1)*cos(n*omega*time_wave) - n^2*omega^2*Bn(p*n+1)*sin(n*omega*time_wave);
    
    %Dynamic moduli
    Gp_t=Gp_t+cos(omega*time_wave).*(n*(Gp_n(p*n+1)*cos(n*omega*time_wave)-Gpp_n(p*n+1)*sin(n*omega*time_wave)))...
        +sin(omega*time_wave).*(n^2*(Gp_n(p*n+1)*sin(n*omega*time_wave)+Gpp_n(p*n+1)*cos(n*omega*time_wave)));
    Gpp_t=Gpp_t-sin(omega*time_wave).*(n*(Gp_n(p*n+1)*cos(n*omega*time_wave)-Gpp_n(p*n+1)*sin(n*omega*time_wave)))...
        +cos(omega*time_wave).*(n^2*(Gp_n(p*n+1)*sin(n*omega*time_wave)+Gpp_n(p*n+1)*cos(n*omega*time_wave)));
    
    %Derivatives of dynamic moduli
    Gp_t_dot=Gp_t_dot-omega*n*sin(omega*time_wave).*(Gp_n(p*n+1)*cos(n*omega*time_wave)-Gpp_n(p*n+1)*sin(n*omega*time_wave))...
        +omega*n^3*sin(omega*time_wave).*(Gp_n(p*n+1)*cos(n*omega*time_wave)-Gpp_n(p*n+1)*sin(n*omega*time_wave));
    Gpp_t_dot=Gpp_t_dot-omega*n*cos(omega*time_wave).*(Gp_n(p*n+1)*cos(n*omega*time_wave)-Gpp_n(p*n+1)*sin(n*omega*time_wave))...
        +omega*n^3*cos(omega*time_wave).*(Gp_n(p*n+1)*cos(n*omega*time_wave)-Gpp_n(p*n+1)*sin(n*omega*time_wave));
end

%Reconstruct the Tangent vector
Rec_Tvect = Recon_Wave_dot;

%Reconstruct the Normal vector
Rec_Nvect = [(Recon_Wave_ddot(:,1).*(Recon_Wave_dot(:,2)).^2+Recon_Wave_ddot(:,1).*(Recon_Wave_dot(:,3)).^2-...
    Recon_Wave_ddot(:,2).*Recon_Wave_dot(:,1).*Recon_Wave_dot(:,2)-Recon_Wave_ddot(:,3).*Recon_Wave_dot(:,1).*Recon_Wave_dot(:,3))./...
    (((Recon_Wave_dot(:,1)).^2+(Recon_Wave_dot(:,2)).^2+(Recon_Wave_dot(:,3)).^2).^(3/2)),...
    (Recon_Wave_ddot(:,2).*(Recon_Wave_dot(:,1)).^2+Recon_Wave_ddot(:,2).*(Recon_Wave_dot(:,3)).^2-...
    Recon_Wave_ddot(:,1).*Recon_Wave_dot(:,1).*Recon_Wave_dot(:,2)-Recon_Wave_ddot(:,3).*Recon_Wave_dot(:,2).*Recon_Wave_dot(:,3))./...
    (((Recon_Wave_dot(:,1)).^2+(Recon_Wave_dot(:,2)).^2+(Recon_Wave_dot(:,3)).^2).^(3/2)),...
    (Recon_Wave_ddot(:,3).*(Recon_Wave_dot(:,1)).^2+Recon_Wave_ddot(:,3).*(Recon_Wave_dot(:,2)).^2-...
    Recon_Wave_ddot(:,1).*Recon_Wave_dot(:,1).*Recon_Wave_dot(:,3)-Recon_Wave_ddot(:,2).*Recon_Wave_dot(:,2).*Recon_Wave_dot(:,3))./...
    (((Recon_Wave_dot(:,1)).^2+(Recon_Wave_dot(:,2)).^2+(Recon_Wave_dot(:,3)).^2).^(3/2))];

%Reconstruct the Binormal vector
Rec_Bvect = [(Recon_Wave_ddot(:,3).*Recon_Wave_dot(:,2)-Recon_Wave_ddot(:,2).*Recon_Wave_dot(:,3))./...
    ((Recon_Wave_dot(:,1)).^2+(Recon_Wave_dot(:,2)).^2+(Recon_Wave_dot(:,3)).^2),...
    (Recon_Wave_ddot(:,1).*Recon_Wave_dot(:,3)-Recon_Wave_ddot(:,3).*Recon_Wave_dot(:,1))./...
    ((Recon_Wave_dot(:,1)).^2+(Recon_Wave_dot(:,2)).^2+(Recon_Wave_dot(:,3)).^2),...
    (Recon_Wave_ddot(:,2).*Recon_Wave_dot(:,1)-Recon_Wave_ddot(:,1).*Recon_Wave_dot(:,2))./...
    ((Recon_Wave_dot(:,1)).^2+(Recon_Wave_dot(:,2)).^2+(Recon_Wave_dot(:,3)).^2)];

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

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

SPPplus_print_v1(filename,M,omega,Recon_Wave,Gp_t,Gpp_t,G_speed,Gp_t_dot,...
    Gpp_t_dot,print_state(1),Rec_Tvect,Rec_Nvect,Rec_Bvect,'FOURIER',...
    'Data calculated via Fourier-transformation','Harmonics included up to ')

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

getfigs = [1,1]; %specifies to produce Standard SPP figure
                %(ex. SPP_LAOS_FOURIER_PLOT_sample_file1.jpg),
                % and Fourier Harmonics figure
                %(ex. SPP_LAOS_FOURIER_HARMONICS_sample_file1.jpg)

SPPplus_figures_v1(resp_wave,Recon_Wave,Gp_t,Gpp_t,G_speed,Gp_t_dot,...
    Gpp_t_dot,time_wave,f_domain,ft_amp_sign_fft_resp_s,M,W,print_state(2),...
    'FOURIER',filename,getfigs)

end