%%%%%   Jakes Model  %%%%%
function [model,time]=JakesModel(p,Fd,T)  

%%% fd         :  Max Doppler shift 
%%% T          :  Sampling interval 
%%% num        :  Sampling Point 
%%% P          :  Factor of Correlation Time/Sampling Interval(if want to observe the time selective fading, observe time must longer than Correlation Time, so, total number of sampling points must larger than P) 
%%% M          :  Order of Jakes Model  (initially set M = 16) 
%%% N          :  Number of the arrive phases that are uniformly distributed   (N = 4*M + 2)  
%%% a          :  Max arrive phase of the multipaths  (initially set a = 0)
%%% Begin      :  Simulation begining time 
%%% beta_n     :  Arrive phase of the nth. multipath 
%%% fn         :  Doppler shift of the nth. multipath 
%%% h_I        :  inphase component of the model 
%%% h_Q        :  quadrature component of the model 
%%% model      :  Jakes Model outcome 
%%% time       :  Simulation time 

%%%%%   Define Parameters  %%%%% 
fd = Fd; P = p; 
M = 8; 
N = 4*M + 2; 
a = 0; 
Begin = 1000; 
Tc = 1/fd;              % Correlation Time 
T = Tc/P;               % ampling Interval Time << Correlation Time 
num = 0; 

%%%%%  Model Design  %%%%% 
for t = Begin+T:T:(Begin+100*Tc)     % insure that the observe time is longer than a Correlation Time 
    num = num + 1; 
    h_I = 0;
    h_Q = 0; 
  for n = 1:M 
    beta_n = pi*n/M;  
    fn = fd*cos(2*pi*n/N); 
    h_I = h_I + sqrt(2)*cos(beta_n)*cos(2*pi*fn*t);
    h_Q = h_Q + sqrt(2)*sin(beta_n)*cos(2*pi*fn*t);
%     h_I = h_I + sqrt(2/(M+1))*cos(beta_n)*cos(2*pi*fn*t);  % Energy normolized 
%     h_Q = h_Q + sqrt(2/(M))*sin(beta_n)*cos(2*pi*fn*t);    % Energy normolized 
  end
  h_I = h_I + sqrt(1/(M+1))*cos(a)*cos(2*pi*fd*t);
  h_Q = h_Q + sqrt(1/(M))*sin(a)*cos(2*pi*fd*t); 
  Model(num) = h_I + h_Q*i; 
end 
  
Time = T*(1:length(Model)); 
 
%%%%%  Return  %%%%% 
model = Model; 
time = Time;