%
% Place local Clusters
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% r..vector of the BS/MS/cluster
% d..distance in between the following objects
% a..main diagonal of the ellipsoid (delay)
% b..diagonal of the ellipsoid (azimuth)
% h..height of the unshifted ellipsoid
% Theta,Phi..orinentation of the ellipsoid in space

if N_C_local>0      % place the local cluster around the MS
    r_C_MS(1,:)=r_MS;
    r_C_BS(1,:)=r_MS;    % the local cluster is at the same position wether seen from the BS-side or the MS-side.
    h_C_MS(1)=d_BS_MS.*tan(phi_C_BS(1));    % height of local MS cluster
    h_C_BS(1)=h_C_MS(1);                    % height of local MS cluster
    a_C_MS(1)=d_tau(1);     % the local cluster is circular around the MS
    a_C_BS(1)=d_tau(1);
    b_C_MS(1)=d_tau(1);     % both, BS-side and MS-side cluster, are at exactly the same position
    b_C_BS(1)=d_tau(1);
    Phi_C_MS(1)=Phi_BS_MS-pi; % for later computations; to allow for the same procedure as for the twin clusters
    Phi_C_BS(1)=Phi_BS_MS;
    Theta_C_MS(1)=Theta_BS_MS-pi;
    Theta_C_BS(1)=Theta_BS_MS;
    tau_C_link(1)=0;        % there is no cluster-link delay for the local clusters
    k_local=1;
end % if N_C_local>0
if N_C_local==2     % place the local cluster around the BS
    r_C_MS(2,:)=r_BS;
    r_C_BS(2,:)=r_BS;    % the local cluster is at the same position wether seen from the BS-side or the MS-side.
    h_C_MS(2)=d_BS_MS.*tan(phi_C_BS(2));      % height of local MS cluster
    h_C_BS(2)=d_BS_MS.*tan(phi_C_BS(2));      % height of local MS cluster
    a_C_MS(2)=d_tau(2);
    a_C_BS(2)=d_tau(2);
    b_C_MS(2)=d_tau(2);
    b_C_BS(2)=d_tau(2);
    Phi_C_MS(2)=Phi_BS_MS; 
    Phi_C_BS(2)=Phi_BS_MS-pi; % for later computations; to allow for the same procedure as for the twin clusters
    Theta_C_MS(2)=Theta_BS_MS;
    Theta_C_BS(2)=Theta_BS_MS-pi;
    tau_C_link(2)=0;        % there is no cluster-link delay for the local clusters
    k_local=1:2;
end % if N_C_local==2
if N_C_local==0     % if no local clusters exist k_local has to be initialized
    k_local=[];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Place Twin clusters
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if N_C_twin>0
    k_twin=(N_C_local+1):(N_C_local+N_C_twin);    % k_twin keeps all the cluster numbers of twin clusters
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % the following lines should be replaced by the marginal pdfs of the
    % distributions
    Phi_C_MS(k_twin)=rand(N_C_twin,1)*1;  % use instead of marginal angular pdfs
    Phi_C_BS(k_twin)=rand(N_C_twin,1)*1;
    Theta_C_MS(k_twin)=rand(N_C_twin,1)*2*pi;
    Theta_C_BS(k_twin)=rand(N_C_twin,1)*2*pi;
    tau_C(k_twin)=tau_0+rand(N_C_twin,1)*5e-7;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    a_C_MS(k_twin)=d_tau(k_twin);
    a_C_BS(k_twin)=d_tau(k_twin);
    if 0    % if a distribution of cluster distances from BS/MS exists use it
        
    else    % if not the clusters are assumed to be circular in the delay-azimuth plane
        d_C_MS(k_twin)=d_tau(k_twin)./(2*tan(theta_C_MS(k_twin)));  % distance of cluster from MS
        d_C_BS(k_twin)=d_tau(k_twin)./(2*tan(theta_C_BS(k_twin)));  % distance of cluster from BS
        b_C_MS(k_twin)=a_C_MS(k_twin);
        b_C_BS(k_twin)=a_C_BS(k_twin);
    end
    [temp(:,1),temp(:,2),temp(:,3)]=sph2cart(Theta_C_MS(k_twin),Phi_C_MS(k_twin),d_C_MS(k_twin));
    r_C_MS(k_twin,:)=repmat(r_MS,N_C_twin,1)+temp;
    [temp(:,1),temp(:,2),temp(:,3)]=sph2cart(Theta_C_BS(k_twin),Phi_C_BS(k_twin),d_C_BS(k_twin));
    r_C_BS(k_twin,:)=repmat(r_BS,N_C_twin,1)+temp;
    clear temp
    h_C_MS(k_twin)=d_C_MS(k_twin).*tan(phi_C_MS(k_twin));   % height of MS-side cluster
    h_C_BS(k_twin)=d_C_BS(k_twin).*tan(phi_C_BS(k_twin));   % height of BS-side cluster
    
    tau_C_link(k_twin)=tau_C(k_twin)-d_C_MS(k_twin)/c0-d_C_BS(k_twin)/c0;
else
    k_twin=[];

end % if N_C_twin>0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Place single interacting clusters
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if N_C_single>0

    k_single=(N_C_local+N_C_twin+1):N_C;    % k_single keeps all the cluster numbers of the single interacting clusters
    
    d_C_BS(k_single)=r_min-sigma_r.*log(1-rand(N_C_single,1)); % distance of cluster from BS
    Theta_C_BS(k_single)=Theta_BS_MS+randn(N_C_single,1).*sigma_Theta_C_BS;
    [temp(:,1),temp(:,2),temp(:,3)]=sph2cart(Theta_C_BS(k_single),Phi_C_BS(k_single),d_C_BS(k_single));
    r_C_BS(k_single,:)=repmat(r_BS,N_C_single,1)+temp;
    clear temp
    r_C_MS(k_single,:)=r_C_BS(k_single,:); % the single interacting cluster is at the same position at the MS side as seen from the BS side.
        
    h_C_BS(k_single)=d_C_BS(k_single).*tan(phi_C_BS(k_single));      % height of cluster
    h_C_MS(k_single)=h_C_BS(k_single);      % height of cluster
    b_C_BS(k_single)=d_C_BS(k_single).*tan(theta_C_BS(k_single));
    b_C_MS(k_single)=b_C_BS(k_single);
    a_C_BS(k_single)=d_tau(k_single);
    a_C_MS(k_single)=a_C_BS(k_single);
    
    Phi_C_BS(k_single)=asin((-r_BS(3)+r_C_BS(k_single,3))./d_C_BS(k_single)); % points from the BS towards the SI cluster
    Phi_C_MS(k_single)=Phi_C_BS(k_single);   % again this value is only usedfor rotating the cluster but not for any positioning
    Theta_C_MS(k_single)=Theta_C_BS(k_single);       % again this value is only usedfor rotating the cluster but not for any positioning
    tau_C_link(k_single)=0;     % there is no cluster link delay for single interacting clusters
else
    k_single=[];
end % if N_C_single>0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Place MPCs
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k=k_twin
    r=randn(N_MPC,3);
    r_MPC_MS_rel(k,:,:)=-r*diag([a_C_MS(k); b_C_MS(k); h_C_MS(k)])*rotation_matrix(0,Phi_C_MS(k),Theta_C_MS(k));
    r_MPC_MS(k,:,:)=r_MPC_MS_rel(k,:,:)+repmat(permute(r_C_MS(k,:),[1 3 2]),[1 N_MPC 1]);
    r_MPC_BS_rel(k,:,:)=r*diag([a_C_BS(k); b_C_BS(k); h_C_BS(k)])*rotation_matrix(0,Phi_C_BS(k),Theta_C_BS(k));
    r_MPC_BS(k,:,:)=r_MPC_BS_rel(k,:,:)+repmat(permute(r_C_BS(k,:),[1 3 2]),[1 N_MPC 1]);
end
for k=[k_local k_single]
    r=randn(N_MPC,3);
    r_MPC_MS_rel(k,:,:)=r*diag([a_C_MS(k); b_C_MS(k); h_C_MS(k)])*rotation_matrix(0,Phi_C_MS(k),Theta_C_MS(k));
    r_MPC_MS(k,:,:)=r_MPC_MS_rel(k,:,:)+repmat(permute(r_C_MS(k,:),[1 3 2]),[1 N_MPC 1]);
    r_MPC_BS_rel(k,:,:)=r*diag([a_C_BS(k); b_C_BS(k); h_C_BS(k)])*rotation_matrix(0,Phi_C_BS(k),Theta_C_BS(k));
    r_MPC_BS(k,:,:)=r_MPC_BS_rel(k,:,:)+repmat(permute(r_C_BS(k,:),[1 3 2]),[1 N_MPC 1]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute attenuation of MPCs
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Attenuation of clusters
P_C=max(exp(-k_tau*(tau_C-tau_0)*1e6),exp(-k_tau*(tau_B-tau_0)*1e6));

% Power factors of MPCs within a cluster; Ricean distribution
A_K=repmat(K_MPC,N_C,N_MPC);
P_S=abs(randn(N_C,N_MPC)*sigma_K+A_K+j*randn(N_C,N_MPC)*sigma_K);
P_S=P_S/sqrt(mean(abs(P_S(:).^2)));
% Shadow fading of each cluster
SF=10.^(0.1*randn(N_C,1)*sigma_S);
SF=repmat(SF,1,N_MPC);

% Attenuation amplitude of each MPC
A_MPC=sqrt(repmat(P_C,1,N_MPC).*P_S.*SF);
A_MPC=A_MPC/sqrt(sum(abs(A_MPC(:).^2)));


% Attenuation of the LOS component
A_LOS=mu_K+10.^(0.05*rand(1)*sigma_K);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute the impulse responses
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k2=1:No_of_imp_resp      % compute 5 impulse responses; impulse responses are only valid if teh movement of teh MS is not too large!
    % MPC distance from BS/MS plus angles for steering vectors
    temp_x=r_MPC_MS(:,:,1)-repmat(permute(r_MS(:,1),[3 2 1]),[N_C,N_MPC,1]);
    temp_y=r_MPC_MS(:,:,2)-repmat(permute(r_MS(:,2),[3 2 1]),[N_C,N_MPC,1]);
    temp_z=r_MPC_MS(:,:,3)-repmat(permute(r_MS(:,3),[3 2 1]),[N_C,N_MPC,1]);
    [Theta_MPC_MS,Phi_MPC_MS,d_MPC_MS]=cart2sph(temp_x,temp_y,temp_z);
    temp_x=r_MPC_BS(:,:,1)-repmat(permute(r_BS(:,1),[3 2 1]),[N_C,N_MPC,1]);
    temp_y=r_MPC_BS(:,:,2)-repmat(permute(r_BS(:,2),[3 2 1]),[N_C,N_MPC,1]);
    temp_z=r_MPC_BS(:,:,3)-repmat(permute(r_BS(:,3),[3 2 1]),[N_C,N_MPC,1]);
    [Theta_MPC_BS,Phi_MPC_BS,d_MPC_BS]=cart2sph(temp_x,temp_y,temp_z);
	
    [Theta_BS_MS,Phi_BS_MS,d_BS_MS]=cart2sph(r_MS(1)-r_BS(1),r_MS(2)-r_BS(2),r_MS(3)-r_BS(3));
	
	% Steering vectors for all MPCs and LOS for BS and MS antennas
    % Since a planar array in XY is assumed the elevation angle of the
    % paths have no influence on the steering vectors
    S_MS=dir_resp(Theta_MPC_MS,MSAnr,d_Ant_MS);
    S_BS=dir_resp(Theta_MPC_BS,BSAnr,d_Ant_BS);
    S_MS_LOS=dir_resp(Theta_BS_MS,MSAnr,d_Ant_MS);
    S_BS_LOS=dir_resp(Theta_BS_MS+pi,BSAnr,d_Ant_BS); % pi once more from above

	% Attenuation of clusters
	P_C=max(exp(-k_tau*(tau_C-tau_0)*1e6),exp(-k_tau*(tau_B-tau_0)*1e6));
    
    % compute the delay of all MPCs
    
    tau_MPC=(d_MPC_MS+d_MPC_BS)/c0+repmat(tau_C_link',1,N_MPC);
    A_MIMO=repmat([S_MS_LOS; S_MS],[1 1 BSAnr]).*repmat([A_LOS; A_MPC(:)],[1 MSAnr BSAnr]).*repmat(permute([S_BS_LOS; S_BS],[1 3 2]),[1 MSAnr 1]);
    tau_MIMO=repmat([tau_0; tau_MPC(:)],[1 MSAnr BSAnr]);
    A_MIMO_MPC=repmat([S_MS],[1 1 BSAnr]).*repmat([A_MPC(:)],[1 MSAnr BSAnr]).*repmat(permute([S_BS],[1 3 2]),[1 MSAnr 1]);
    tau_MIMO_MPC=repmat([tau_MPC(:)],[1 MSAnr BSAnr]);
    A_MIMO_LOS=repmat([S_MS_LOS],[1 1 BSAnr]).*repmat([A_LOS],[1 MSAnr BSAnr]).*repmat(permute([S_BS_LOS],[1 3 2]),[1 MSAnr 1]);
    tau_MIMO_LOS=repmat([tau_0],[1 MSAnr BSAnr]);

    imp_resp=zeros(ceil(max(max(tau_MPC))/t_sample+1),MSAnr,BSAnr); % assume a maximum path delay of 5km
    imp_resp_MPC=zeros(ceil(max(max(tau_MPC))/t_sample+1),MSAnr,BSAnr); % assume a maximum path delay of 5km
    imp_resp_LOS=zeros(ceil(max(max(tau_MPC))/t_sample+1),MSAnr,BSAnr); % assume a maximum path delay of 5km
    % Create impluse response
    % In this implementation all paths are reduced to the t_sample delay.
    % If t_sample is equal to the sampling rate of the receiver the
    % resulting model is simple tapp-delay model. The channel model should
    % be at least four times oversampled compared to the receiver since no
    % sinc-function is implemented. 
    for k=1:length(tau_MIMO)
        imp_resp((ceil(tau_MIMO(k,1,1)./t_sample+1)),:,:)=imp_resp((ceil(tau_MIMO(k,1,1)./t_sample+1)),:,:)+...
            A_MIMO(k,:,:).*exp(j*2*pi*mod(tau_MIMO(k,:,:).*f_c,1));
    end
    for k=1:length(tau_MIMO_MPC)
        imp_resp_MPC((ceil(tau_MIMO_MPC(k,1,1)./t_sample+1)),:,:)=imp_resp_MPC((ceil(tau_MIMO_MPC(k,1,1)./t_sample+1)),:,:)+...
            A_MIMO_MPC(k,:,:).*exp(j*2*pi*mod(tau_MIMO_MPC(k,:,:).*f_c,1));
    end
    for k=1
        imp_resp_LOS((ceil(tau_MIMO_LOS(k,1,1)./t_sample+1)),:,:)=imp_resp_LOS((ceil(tau_MIMO_LOS(k,1,1)./t_sample+1)),:,:)+...
            A_MIMO_LOS(k,:,:).*exp(j*2*pi*mod(tau_MIMO_LOS(k,:,:).*f_c,1));
    end
    imp_resp=imp_resp/sqrt(sum(abs(imp_resp(:)).^2)/MSAnr/BSAnr);   % normalize impulse response power to one.
    
    % filter the impulse response with a cosine filter. Only used if the 
    %   last line is uncommented!
    imp_filter = rcosflt(reshape(imp_resp,[length(imp_resp),MSAnr*BSAnr]), cosf.Fd, cosf.Fs, 'filter/Fs', cosf.ys);
    imp_filter=reshape(imp_filter,[length(imp_filter),MSAnr,BSAnr]);
  %  keyboard
            impn(k2,1:size(imp_resp,1),:,:)=permute(imp_resp,[4 1 2 3]);
            imp_LOS(k2,1:size(imp_resp_LOS,1),:,:)=permute(imp_resp_LOS,[4 1 2 3]);
            imp_MPC(k2,1:size(imp_resp_MPC,1),:,:)=permute(imp_resp_MPC,[4 1 2 3]);

            
            
    imp_cosine(k2,1:size(imp_filter(1:4:end,:,:),1),:,:)=permute(imp_filter(1:4:end,:,:),[4 1 2 3]);
    
   
    
    r_MS=r_MS+v_MS*delta_t;     % update the MS position. The MS is the only moving part in macrocells!

end
imp.imp=impn;
imp.imp_MPC=imp_MPC;
imp.imp_LOS=imp_LOS;
imp.imp_cosine=imp_cosine;
%imp.imp=imp_cosine;    % uncomment this line for a filtered impulse
%response (also given in imp.imp_cosine)