function [F] = gen_seasonal_trigonometric_F(s,m,omega)
% GEN_SEASONAL_TRIGONOMETRIC_F - Gen(erate) Seasonal Trigonometric F matrix
%
% Input: 
%   s     : total number of seasons = seasonal indicators.
%   m    : the number of trigonometric components to keep/consider.
%   omega : relaxation coefficient
% 
% Our assumed model is a locally constant LINEAR model with trigonometric predictors given by 
% 
% z_{n+j} = \beta_0 + \beta_1 j + \sum_{i=1}^{m} \beta_{1i} \sin( f_i t) + \beta_{2i} \cos( f_i t ) + \varepsilon_{n+j} 
% 
% See Page 164 of the book by Abraham on the section "Loally Constant Seasonal Models Using Trigonmetric Functions" 
% 
% Note: rather than symbolically calculate the given sums we will numerically evaluate them.
%   
% Written by:
% -- 
% John L. Weatherwax                2009-06-01
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

if( omega>1 )
  fprint('In this implementation we require omega<1\n');
  F=[]; 
  return; 
end

f = 2*pi/s;

% this code creates the matrix F: 
% 
F = zeros( 2 + 2*m ); 
F(1,1) =       1/( 1-omega    ); 
F(1,2) = -omega/( (1-omega)^2 ); 
F(2,2) = omega*(1+omega)/( (1-omega)^3 );
F(2,1) = F(1,2); 

for ii=1:m
  f_ii = f*ii; 
  g_1  = 1 - 2 * omega * cos( f_ii ) + omega^2;
  F(1,1+2*ii) = - omega * sin( f_ii ) / g_1; 
  F(1,2+2*ii) = ( 1 - omega * cos( f_ii ) ) / g_1; 
  F(2,1+2*ii) = omega * ( 1 - omega^2 ) * sin( f_ii ) / ( g_1^2 ); 
  F(2,2+2*ii) = - ( omega * ( 1 + omega^2 ) * cos( f_ii ) - 2 * omega^2 ) / ( g_1^2 );
  
  F(1+2*ii,1) = F(1,1+2*ii); 
  F(2+2*ii,1) = F(1,2+2*ii); 
  F(1+2*ii,2) = F(2,1+2*ii); 
  F(2+2*ii,2) = F(2,2+2*ii); 
end

for ii=1:m
  f_ii = f * ii; 
  for kk=1:m
    f_kk = f * kk; 
    g_2  = 1 - 2 * omega * cos( f_ii - f_kk ) + omega^2; 
    g_3  = 1 - 2 * omega * cos( f_ii + f_kk ) + omega^2; 
    
    F(1+2*ii,1+2*kk) =  0.5 * ( ( 1 - omega * cos( f_ii - f_kk ) )/g_2 - ( 1 - omega * cos( f_ii + f_kk ) )/g_3 );  
    F(2+2*ii,2+2*kk) =  0.5 * ( ( 1 - omega * cos( f_ii - f_kk ) )/g_2 + ( 1 - omega * cos( f_ii + f_kk ) )/g_3 );  
    F(1+2*ii,2+2*kk) =  -0.5 * ( omega * sin( f_ii - f_kk ) /g_2 + omega * sin( f_ii + f_kk )/g_3 );
    
    F(1+2*kk,1+2*ii) = F(1+2*ii,1+2*kk);
    F(2+2*kk,2+2*ii) = F(2+2*ii,2+2*kk);
    F(2+2*kk,1+2*ii) = F(1+2*ii,2+2*kk);
  end
end