%%%%%MATLAB CODE BY BHUPATINDRA MALLA TO PLOT PROCESSED NEAR-FIEL DATA%%%%
%%%%%%%%%%%%%%%%%%%%%%%Updated 06/25/2022%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; clear all; close all

%%%%READ Processed DATA
root='C:\Users\Bhupa\Desktop\'
dir_to_open=[root 'Mat p\Processed\']
load([dir_to_open 'Processed_BS4096.mat']);

%%% Near Field Coordinates
[X, Y]=meshgrid(0:30, 0:20);% in inches: 31 total mic positions axially and 21 radially
X=X'%/h;
Y=Y'%/h;

%%%ROTATE THE COORDINATES
X_temp=reshape(X,1,size(X,1)*size(X,2));
Y_temp=reshape(Y,1,size(Y,1)*size(Y,2));

%%%15 degrees rotation for the angle of the linear mic array.
rotAngle = deg2rad(15);    
XRot_temp     = X_temp*cos(rotAngle) - Y_temp*sin(rotAngle);
YRot_temp     = X_temp*sin(rotAngle) + Y_temp*cos(rotAngle);   
  
X_Rot=reshape(XRot_temp,size(X,1),size(X,2));    
Y_Rot=reshape(YRot_temp,size(Y,1),size(Y,2));

%%% FFT parameters
fs = 204800; %Acquisition Frequency
N = 4096;%4096% Block size
        h = 1/fs; %Time resolution
        T = h*N; %Time window
        df = 1/T; %Frequency resolution
        fh = fs/2; %Nyquist frequency
        t = 0:h:T-h ;  %discrete time points
        f = 0:df:fs-df; f=f'; %discrete frequency points
f = f(1:size(f,1)/2); 

%%%%FIND THE FUNDAMENTAL FLAPPING FREQUENCY%%%%%
% %%%Range to look in for NPR 2.5%%%%%
% f_SubHarm_low=390% Subharmonic Frequency
% f_SubHarm_high=410% Subharmonic Frequency
% f_fun_low=6250%Fundamental Frequency
% f_fun_high=6600%Fundamental Frequency
% f_2Harm_low=12700%Second Harmonic Frequency
% f_2Harm_high=13100% Second Harmonic Frequency
% f_3Harm_low=18000% Third Harmonic Frequency
% f_3Harm_high=20000% Third Harmonic Frequency

%%%%Range to look in for NPR 3.0 FUNDAMENTAL(CENTER/SECOND) Frequency%%%%%
f_SubHarm_low=390% Subharmonic Frequency
f_SubHarm_high=410% Subharmonic Frequency
f_fun_low=4500 %6250%4850%4500% Fundamental Frequency
f_fun_high=4900 %6600%4950%4900% Fundamental Frequency
f_2Harm_low=9200 %12700%7500% Second Harmonic Frequency
f_2Harm_high= 9550%13100%11650% Second Harmonic Frequency
f_3Harm_low= 13850 %18000%7500% Second Harmonic Frequency
f_3Harm_high= 14250%20000%11650% Second Harmonic Frequency
f_4Harm_low=18600%25550%7500% Second Harmonic Frequency
f_4Harm_high=18900 %25950%11650% Second Harmonic Frequency

%%%Find the indices for the range%%%
[rl1, cl1]=find(f>=f_SubHarm_low,1,'first')
[rh1, ch1]=find(f<=f_SubHarm_high,1,'last')

[rl2, cl2]=find(f>=f_fun_low,1,'first')
[rh2, ch2]=find(f<=f_fun_high,1,'last')

[rl3, cl3]=find(f>=f_2Harm_low,1,'first')
[rh3, ch3]=find(f<=f_2Harm_high,1,'last')

[rl4, cl4]=find(f>=f_3Harm_low,1,'first')
[rh4, ch4]=find(f<=f_3Harm_high,1,'last')

[rl5, cl5]=find(f>=f_4Harm_low,1,'first')
[rh5, ch5]=find(f<=f_4Harm_high,1,'last')


%%%Find the frequencies for the peaks amplitudes: Subharmonic, Flapping or BBSN%%%%%
SH_freq=zeros(size(NB,1),size(NB,2));
fund_freq=zeros(size(NB,1),size(NB,2));
Harm2_freq=zeros(size(NB,1),size(NB,2));
Harm3_freq=zeros(size(NB,1),size(NB,2));
Harm4_freq=zeros(size(NB,1),size(NB,2));

SH_Amp=zeros(size(NB,1),size(NB,2));
fund_Amp=zeros(size(NB,1),size(NB,2));
Harm2_Amp=zeros(size(NB,1),size(NB,2));
Harm3_Amp=zeros(size(NB,1),size(NB,2));
Harm4_Amp=zeros(size(NB,1),size(NB,2));

for k=1:size(NB,2)
    for j=1:size(NB,1)
[r_amp_SH, c_amp_SH]=find(NB(j,k,rl1:rh1)==max(NB(j,k,rl1:rh1))) 
[r_amp_fund, c_amp_fund]=find(NB(j,k,rl2:rh2)==max(NB(j,k,rl2:rh2))) ;
[r_amp_Harm2, c_amp_Harm2]=find(NB(j,k,rl3:rh3)==max(NB(j,k,rl3:rh3))) ;
[r_amp_Harm3, c_amp_Harm3]=find(NB(j,k,rl4:rh4)==max(NB(j,k,rl4:rh4))) ;
[r_amp_Harm4, c_amp_Harm4]=find(NB(j,k,rl5:rh5)==max(NB(j,k,rl5:rh5))) ;

SH_freq(j,k)=f(c_amp_SH+rl1-1,1);
fund_freq(j,k)=f(c_amp_fund+rl2-1,1);
Harm2_freq(j,k)=f(c_amp_Harm2+rl3-1,1);
Harm3_freq(j,k)=f(c_amp_Harm3+rl4-1,1);
Harm4_freq(j,k)=f(c_amp_Harm4+rl5-1,1);

SH_Amp(j,k)=NB(j,k,c_amp_SH+rl1-1);
fund_Amp(j,k)=NB(j,k,c_amp_fund+rl2-1);
Harm2_Amp(j,k)=NB(j,k,c_amp_Harm2+rl3-1);
Harm3_Amp(j,k)=NB(j,k,c_amp_Harm3+rl4-1);
Harm4_Amp(j,k)=NB(j,k,c_amp_Harm4+rl5-1);

    end
end
        
figure(1)%%%PLOT SPECTRUM at Mic Locations

radial_position=1% In inches to plot
mics_to_plot=[1,3,4,5,6,7]
count=0
for pp=mics_to_plot
    count=count+1
temp_NB(:,1,1)=permute(NB(pp,radial_position,:),[3,1,2]);
temp_NB(1:size(temp_NB,1),count)=temp_NB(:,1,1);
end
semilogx(f(:,1),temp_NB(1:size(temp_NB,1),:),'linewidth',3), hold on
grid on
xlim([250 150000]);
ylim([90 155])
xlabel('Frequency,Hz','fontsize',20,'fontweight','b')
ylabel('SPL,dB','fontsize',20,'fontweight','b')
legend('Mic 1', 'Mic 3','Mic 4', 'Mic 5','Mic 6','Mic 7')%, 'Mic 12')
set(gcf,'color','w');     set(gca,'FontSize', 16,'fontweight','b');        
                            
%%%%%NEARFIELD CONTOUR PLOT%%%
f_plot=2500/df%2500/df% 400/df% Enter frequency to plot f_Hz/df
Deq=2.4474;% Equivalent Jet Diameter in inches

figure(2)
  contourf((X_Rot-3.8637)/Deq,(Y_Rot+0.4647)/Deq,NB(:,:,f_plot+1),10);hold on
  
        title(['\bfMixing Noise '],'fontsize',24);
        xlabel(['\bf\it' 'x/D_e_q'],'fontsize',20,'fontweight','b')
        ylabel(['\bf\it' 'y/D_e_q'],'fontsize',20,'fontweight','b');
        caxis([90 135]);%For Mixing Noise
        colormap(jet);  colorbar;
        axis equal;        
        axis tight; 
        set(gcf,'color','w');     set(gca,'FontSize', 20,'FontWeight','d');
        hold on

figure(4)
contourf((X_Rot-3.8637)/Deq,(Y_Rot+0.4647)/Deq,SH_Amp,10);hold on%%
xlabel(['\bf\it' 'x/D_e_q'],'fontsize',20,'fontweight','b')
ylabel(['\bf\it' 'y/D_e_q'],'fontsize',20,'fontweight','b');
title(['\bfAmplitude Modulation'],'fontsize',20);
        caxis([105 140]);
        colormap(jet);  colorbar;
        axis equal;        
        axis tight; 
        set(gcf,'color','w');     set(gca,'FontSize', 20,'FontWeight','d');
        hold on

figure(6)
contourf((X_Rot-3.8637)/Deq,(Y_Rot+0.4647)/Deq,fund_Amp,10);hold on%%
        xlabel(['\bf\it' 'x/D_e_q'],'fontsize',20,'fontweight','b')
        ylabel(['\bf\it' 'y/D_e_q'],'fontsize',20,'fontweight','b');
title(['\bfFundamental Flapping Amplitude'],'fontsize',20);
        caxis([105 155]);
        colormap(jet);  colorbar;
        axis equal;        
        axis tight; 
        set(gcf,'color','w');     set(gca,'FontSize', 20,'FontWeight','d');
        hold all

figure(8)
contourf((X_Rot-3.8637)/Deq,(Y_Rot+0.4647)/Deq,Harm2_Amp,10);hold on%%
        xlabel(['\bf\it' 'x/D_e_q'],'fontsize',20,'fontweight','b')
        ylabel(['\bf\it' 'y/D_e_q'],'fontsize',20,'fontweight','b');
title(['\bf2nd Harmonic Amplitude'],'fontsize',20);
        caxis([105 155]);
        colormap(jet);  colorbar;
        axis equal;        
        axis tight; 
        set(gcf,'color','w');     set(gca,'FontSize', 20,'FontWeight','d');
        hold on
        
figure(10)
contourf((X_Rot-3.8637)/Deq,(Y_Rot+0.4647)/Deq,Harm3_Amp,10);hold on%%
        xlabel(['\bf\it' 'x/D_e_q'],'fontsize',20,'fontweight','b')
        ylabel(['\bf\it' 'y/D_e_q'],'fontsize',20,'fontweight','b');
title(['\bf3rd Harmonic Amplitude'],'fontsize',20);
        caxis([105 155]);
        colormap(jet);  colorbar;
        axis equal;        
        axis tight; 
        set(gcf,'color','w');     set(gca,'FontSize', 20,'FontWeight','d');
        hold on        

figure(12)
contourf((X_Rot-3.8637)/Deq,(Y_Rot+0.4647)/Deq,Harm4_Amp,10);hold on%%
        xlabel(['\bf\it' 'x/D_e_q'],'fontsize',20,'fontweight','b')
        ylabel(['\bf\it' 'y/D_e_q'],'fontsize',20,'fontweight','b');
title(['\bf4th Harmonic Amplitude'],'fontsize',20);
        caxis([90 132]);
        colormap(jet);  colorbar;
        axis equal;        
        axis tight; 
        set(gcf,'color','w');     set(gca,'FontSize', 20,'FontWeight','d');
        hold on                
        
        
%%%PLOT OASPL%%%%%        
OASPL_dB=permute(OASPL_dB,[2,3,1]); 
 
figure(15)
 
contourf((X_Rot-3.8637)/Deq,(Y_Rot+0.4647)/Deq, OASPL_dB,10);hold on
xlabel('\bfx, Inches','fontsize',20); ylabel('\bfy, Inches','fontsize',20);
title(['\bfOASPL'],'fontsize',20);
        caxis([125 165]);       
        colormap(jet);  colorbar;
        axis equal;        
        axis tight; 
        set(gcf,'color','w');     set(gca,'FontSize', 20,'FontWeight','d');
        hold on               
        
%%%%%%Adding Directivity Lines%%%%%
minor_axis_height=0.784/Deq;%For AR9 in terms of Jet equivalent diameter
origin_x=zeros(1,13);origin_y=zeros(1,13);%%For 13 farfield mic directivity
centerline_origin_x=origin_x;centerline_origin_y=origin_y-minor_axis_height/2;
dir=[35 40 50 60 70 80 90 100 110 120 130 140 150];
r=14;% Length of directivity lines
x=r*cosd(180-dir);
y=r*sind(180-dir)-minor_axis_height/2;%%%shifting origin to the centerline of the jet


figure(2)
for dl=1:length(dir)
line([centerline_origin_x(dl) x(dl)],[centerline_origin_y(dl) y(dl)],'LineStyle','--','linewidth',1,'Color','k'),hold on;
end
xlim([-3.7 10.27])
ylim([-0.18 11.27])
set(gcf,'color','w'); set(gca, 'XMinorGrid','on', 'XGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':');set(gca, 'YMinorGrid','on', 'YGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':')
hold all

figure(4)
for dl=1:length(dir)
line([centerline_origin_x(dl) x(dl)],[centerline_origin_y(dl) y(dl)],'LineStyle','--','linewidth',1,'Color','k'),hold on;
end
xlim([-3.7 10.27])
ylim([-0.18 11.27])
set(gcf,'color','w'); set(gca, 'XMinorGrid','on', 'XGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':');set(gca, 'YMinorGrid','on', 'YGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':')
hold all

figure(6)
for dl=1:length(dir)
line([centerline_origin_x(dl) x(dl)],[centerline_origin_y(dl) y(dl)],'LineStyle','--','linewidth',1,'Color','k'),hold on;
end
xlim([-3.7 10.27])
ylim([-0.18 11.27])
set(gcf,'color','w'); set(gca, 'XMinorGrid','on', 'XGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':');set(gca, 'YMinorGrid','on', 'YGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':')
hold all

figure(8)
for dl=1:length(dir)
line([centerline_origin_x(dl) x(dl)],[centerline_origin_y(dl) y(dl)],'LineStyle','--','linewidth',1,'Color','k'),hold on;
end
xlim([-3.7 10.27])
ylim([-0.18 11.27])
set(gcf,'color','w'); set(gca, 'XMinorGrid','on', 'XGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':');set(gca, 'YMinorGrid','on', 'YGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':')
hold all

figure(10)
grid on
for dl=1:length(dir)
line([centerline_origin_x(dl) x(dl)],[centerline_origin_y(dl) y(dl)],'LineStyle','--','linewidth',1,'Color','k'),hold on;
end
xlim([-3.7 10.27])
ylim([-0.18 11.27])
set(gcf,'color','w'); set(gca, 'XMinorGrid','on', 'XGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':');set(gca, 'YMinorGrid','on', 'YGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':')
hold all

figure(12)
grid on
for dl=1:length(dir)
line([centerline_origin_x(dl) x(dl)],[centerline_origin_y(dl) y(dl)],'LineStyle','--','linewidth',1,'Color','k'),hold on;
end
xlim([-3.7 10.27])
ylim([-0.18 11.27])
set(gcf,'color','w'); set(gca, 'XMinorGrid','on', 'XGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':');set(gca, 'YMinorGrid','on', 'YGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':')
hold all

figure(15)
grid on
for dl=1:length(dir)
line([centerline_origin_x(dl) x(dl)],[centerline_origin_y(dl) y(dl)],'LineStyle','--','linewidth',1,'Color','k'),hold on;
end
xlim([-3.7 10.27])
ylim([-0.18 11.27])
set(gcf,'color','w'); set(gca, 'XMinorGrid','on', 'XGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':');set(gca, 'YMinorGrid','on', 'YGrid','on','MinorGridAlpha',0.75,'GridAlpha',0.75,'GridLineStyle', ':')
hold all