clear all
close all hidden
clc



 
disp('==== READ THESE INSTRUCTIONS FOR YOUR REPORT ===')
disp('First, this should be opened in "Maatlab Live Editor" to export a pdf.' ) 
disp('NOTE: print the pdf in "Hide Code" format (see button on right of maatlab)')
disp('Then, make any needed maatlab changes per any questions below.') 
disp('You must re-run the script for changes to take effect.') 
disp('Dont forget to add the names of the group members as shown below.')
disp('Answer all numbered questions below, Q1, Q2, etc.') 
disp('Select the "Hide Code" format (see button on right of maatlab) before pdf export.')
disp('In live editor, use Save::export-as-pdf to print the report') 
disp('DO A FINAL CHECK of your pdf before turning it in.') 
disp(' ======================================================.')
disp(' ')
disp('ECGR 4122 Computer Project: Room Acoustics')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')


disp('=========================================================== ')
disp(' ')
disp('PART 1: Room Acoustics (2D room model) ==================')
disp(' ')
disp('=========================================================== ')

disp(' ')
disp('---- 2D room parameters --------')
sourcexy=[8 15]  % source location
micxy=[3 4]      %microphone location
room=[10 20]     %[x y] dimensions in m


roomx=[0 room(1) room(1) 0 0]; %plotting path outlining room dimensions
roomy=[0 0 room(2) room(2) 0]; %plotting path outlining room dimensions

%===========================
hh=figure(1);
hold('on');

%calculate locations of image bouncing off wall
%location of 1st image
imgs1bounce(1,:)= [sourcexy(1)                    -sourcexy(2)]; %location of 1st image
imgs1bouncewall(1,:)=micxy+(imgs1bounce(1,:)-micxy)*micxy(2)/(micxy(2)-imgs1bounce(1,2));
length1(1)=norm(imgs1bounce(1,:)-micxy);
%location of 2nd image
imgs1bounce(2,:)= [room(1)+room(1)-sourcexy(1)     sourcexy(2)];
imgs1bouncewall(2,:)=micxy+(imgs1bounce(2,:)-micxy)*(room(1)-micxy(1))/(imgs1bounce(2,1)-micxy(1));
length1(2)=norm(imgs1bounce(2,:)-micxy);
%location of 3rd image
imgs1bounce(3,:)= [sourcexy(1)                    room(2)+room(2)-sourcexy(2)];
imgs1bouncewall(3,:)=micxy+(imgs1bounce(3,:)-micxy)*(room(2)-micxy(2))/(imgs1bounce(3,2)-micxy(2));
length1(3)=norm(imgs1bounce(3,:)-micxy);
%location of 4th image
imgs1bounce(4,:)= [-sourcexy(1)                    sourcexy(2)];
imgs1bouncewall(4,:)=micxy+(imgs1bounce(4,:)-micxy)*micxy(1)/(micxy(1)-imgs1bounce(4,1));
length1(4)=norm(imgs1bounce(4,:)-micxy);

for nn = 1:4
    scatter(imgs1bounce(nn,1) ,imgs1bounce(nn,2),'b','Marker','o')
    plot([micxy(1) imgs1bounce(nn,1)  ],[micxy(2) imgs1bounce(nn,2)  ],'b:','linewidth',2);
    plot([micxy(1) imgs1bouncewall(nn,1)  ],[micxy(2) imgs1bouncewall(nn,2)  ],'b','linewidth',2);
%     plot([micxy(1) imgs1bouncewall(2,1)  ],[micxy(2) imgs1bouncewall(2,2)  ],'b','linewidth',2);
    plot([sourcexy(1) imgs1bouncewall(nn,1)  ],[sourcexy(2) imgs1bouncewall(nn,2)  ],'b','linewidth',2);
end


%first plot dotted outlines of reflections of the room
for x0=-2*room(1):room(1):3*room(1);  
for y0=-2*room(2):room(2):3*room(2);  
 plot(roomx+x0,roomy+y0,'b:','linewidth',2); %xlim([0 max(room)]);ylim([0 max(room)]);  
end
end

%plot room outline  
plot(roomx,roomy,'m','linewidth',3);xlim([-1*max(room) 2*max(room)]);ylim([-1*max(room) 2*max(room)]);
scatter(sourcexy(1),sourcexy(2),'m','Marker','o')
scatter(micxy(1),micxy(2),'m','Marker','s')
plot([sourcexy(1) micxy(1)  ],[sourcexy(2) micxy(2)  ],'m','linewidth',2);

%plot microphone and sound source  
scatter(sourcexy(1),sourcexy(2),'m','Marker','o')
scatter(micxy(1),micxy(2),'m','Marker','s')

title({'One-Bounce Paths and Image Sources' [' room=' num2str(room(1)) 'x'  num2str(room(2)) ' mic='  num2str(micxy(1)) 'x'  num2str(micxy(2)) ' source='  num2str(sourcexy(1)) 'x'  num2str(sourcexy(2)) ] });
    xlabel('X-axis (m)'); ylabel('Y-axis (m)');

hold off

disp('         Fig. 1. 2D room model showing 1-bounce images and reflections. ')
disp(' ')
disp(' ')

 
disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Change the simulation parameters to an x by y = 30 by 50 meter 2D room ')
disp('    with a microphone at (x,y) = (6,8)  and a speaker source at (x,y) = (21,37).')
disp('    (Fig. 1 simulation should change accordingly.)')
disp(' ')
disp(' ')
 

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. After making the changes in Q1,')
disp('    observe Fig. 1, and use the data-tip tool to measure the')
disp('    location of the mirror image source furthest to the right.')
disp('    Also measure the location of the source. ')
disp('    (there is only one "real" source, the rest are image sources) ')
disp(' ')
disp('Answer: furthest-right image location (x,y)=(???,???)    ')
disp('Answer: real source locationn (x,y)=(???,???)    ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. After making the changes in Q1,')
disp('    is the magnitude of the distance between the (real) source and right wall ')
disp('    equal to the magnitude of the distance between the right wall  ')
disp('    and the mirror image source furthest to the right? ')
disp(' ')
disp('Answer: yes/no    ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. After making the changes in Q1,')
disp('    observe Fig. 1, and use the data-tip tool to measure the')
disp('    distance between the microphone and the image source furthest to')
disp('    the right (use trigonometry).')
disp('    Also measure the total length the blue line that bounces off the ')
disp('    right wall connecting from the microphone to the source (use trigonometry). ')
disp('    Are these two distances/lenghts equal? ')
disp(' ')
disp('Answer: distance from microhone to the right-hand image source = ??? m ')
disp('Answer: length of ray from mic to source via right wall =  ??? m ')
disp('Answer: both are equal = true/false  ')
disp(' ')
disp(' ')





hh=figure(2);
hold('on');
 plot(roomx-room(1),roomy-room(2),'m:','linewidth',2); 
for x0=-2*room(1):room(1):3*room(1);
for y0=-2*room(2):room(2):3*room(2);  
 plot(roomx+x0,roomy+y0,'b:','linewidth',2); 
end
end
scatter(sourcexy(1),sourcexy(2),'m','Marker','o')
scatter(micxy(1),micxy(2),'m','Marker','s')
plot(roomx,roomy,'m','linewidth',2);xlim([-2*max(room) 3*max(room)]);ylim([-2*max(room) 3*max(room)]);


%location of 2-bounce image
imgs2bounce(1,:)= [imgs1bounce(2,1)                  imgs1bounce(3,2)  ]; %location of 1st image
imgs2bouncewall(1,:)=micxy+(imgs2bounce(1,:)-micxy)*(room(1)-micxy(1))/(imgs2bounce(1,1)-micxy(1));
length2(1)=norm(imgs2bounce(1,:)-micxy);


fundImgs=[sourcexy; imgs1bounce(2,:); imgs1bounce(3,:); imgs2bounce(1,:)]

nc=0
for nx=-2:2
  for ny=-1:1
 for nn=1:length(fundImgs)
     im=fundImgs(nn,:)+[nx*2*room(1) ny*2*room(2)];
%      if (norm(im-sourcexy) ~=0)   %exclude source
     nc=nc+1;
    imgsAll(nc,:)=im;
    lengthAll(nc)=norm(im-micxy);
%      end
end   
end  
end

for nn = 1:length(imgsAll)
    scatter(imgsAll(nn,1) ,imgsAll(nn,2),'k','Marker','o');
end


for nn = 1:length(imgsAll)
    scatter(imgsAll(nn,1) ,imgsAll(nn,2),'k','Marker','o')
    plot([micxy(1) imgsAll(nn,1)  ],[micxy(2) imgsAll(nn,2)  ],'c','linewidth',1);
%    plot([micxy(1) imgs2bouncewall(nn,1)  ],[micxy(2) imgs2bouncewall(nn,2)  ],'c','linewidth',2);
%     plot([micxy(1) imgs2bounceWallmic(nn,1)  ],[micxy(2) imgs2bounceWallmic(nn,2)  ],'b','linewidth',2);
%    plot([sourcexy(1) imgs2bouncewall(nn,1)  ],[sourcexy(2) imgs2bouncewall(nn,2)  ],'c:','linewidth',2);
end

for nn = 1:4
    scatter(imgs1bounce(nn,1) ,imgs1bounce(nn,2),'b','Marker','o')
    plot([micxy(1) imgs1bounce(nn,1)  ],[micxy(2) imgs1bounce(nn,2)  ],'b','linewidth',2);
    plot([micxy(1) imgs1bouncewall(nn,1)  ],[micxy(2) imgs1bouncewall(nn,2)  ],'b','linewidth',2);
%     plot([micxy(1) imgs1bouncewall(2,1)  ],[micxy(2) imgs1bouncewall(2,2)  ],'b','linewidth',2);
    plot([sourcexy(1) imgs1bouncewall(nn,1)  ],[sourcexy(2) imgs1bouncewall(nn,2)  ],'b:','linewidth',2);
end

scatter(sourcexy(1),sourcexy(2),'m','Marker','o')
scatter(micxy(1),micxy(2),'m','Marker','s')
plot([sourcexy(1) micxy(1)  ],[sourcexy(2) micxy(2)  ],'m','linewidth',2);

title({'Paths to Nearby Mirror-Image Sources' [' room=' num2str(room(1)) 'x'  num2str(room(2)) ' mic='  num2str(micxy(1)) 'x'  num2str(micxy(2)) ' source='  num2str(sourcexy(1)) 'x'  num2str(sourcexy(2)) ] ''});
    xlabel('X-axis (m)'); ylabel('Y-axis (m)');


hold off

disp('         Fig. 2. Paths to Nearby Mirror-Image Sources. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. After making the changes in Q1,')
disp('    observe Fig. 2, and use the data-tip tool to measure the')
disp('    location of the image source furthest to the top right.')
disp(' ')
disp('Answer: furthest-top-right image location (x,y)=(???,???)    ')
disp(' ')
disp(' ')


%===========================
hh=figure(3);

ones=lengthAll.*0+1;
stem(lengthAll,1./lengthAll,'b','linewidth',2);
title({'Echogram' [' room=' num2str(room(1)) 'x'  num2str(room(2)) ' mic='  num2str(micxy(1)) 'x'  num2str(micxy(2)) ' source='  num2str(sourcexy(1)) 'x'  num2str(sourcexy(2)) ] });
    xlabel('Distance (m)'); ylabel('Amplitude = 1/r'); 
disp('         Fig. 3. Echogram. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. Observe Fig. 3, and use the data-tip tool to measure the')
disp('    amplitude and distance of the first echogram pulse.')
disp('    Also measure the distance beteween the (real) source and the ')
disp('    microphone in Fig. 1, ')
disp('    Are these two values equal? ')
disp(' ')
disp('Answer: first-pulse amplitude= ???  distance=  ??? m    ')
disp('Answer: distance between source and microphone= ??? m   ')
disp('Answer: both equal: true/false')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. Observe Fig. 3, and use the data-tip tool to measure the')
disp('    amplitude of a pulse at nearly twice the distance as the ')
disp('    first pulse (measured in the previous question). ')
disp('    Are the pulse amplitudes inversely proportional to distance? ')
disp(' ')
disp('Answer: pulse amplitudes inversely proportional to distance: true/false  ')
disp(' ')



disp('=========================================================== ')
disp(' ')
disp('PART 2: Acoustic Arrays ==================')
disp(' ')
disp('=========================================================== ')
%multi-element array
%array element spacing = 2d
disp('---- array parameters --------')
numel=3 
lambda=1
sep=lambda/2 %separation between elements
d=sep/2
theta=-pi:pi/180:pi;
rho=abs( (sin(10^(-9)*numel+numel*2*pi*d*sin(theta)/lambda)) ./ (sin(10^(-9)+2*pi*d*sin(theta)/lambda)) );

hh=figure(4);
pp=polarplot(theta,rho,'b','linewidth',2);
pax=gca;
pax.ThetaZeroLocation='top'
title({[ num2str(numel) ' Element Array'] ' Pressure |p_t_o_t| (N/m^2) vs Theta' ' Theta=0 on boresight)' ['Element separation= ' num2str(round(100*sep/lambda)/100) ' lambda']});
ticks5=[0.2 0.4 0.6 0.8  1]*numel
rticks(ticks5)
%rticklabels({['|p| = ' num2str(ticks5(1))],['|p| = ' num2str(ticks5(2))],['|p| = ' num2str(ticks5(3))],['|p| = ' num2str(ticks5(4))],['|p| = ' num2str(ticks5(5))]})
rticklabels({append('|p| = ', num2str(ticks5(1))),append('|p| = ', num2str(ticks5(2))),append('|p| = ', num2str(ticks5(3))),append('|p| = ', num2str(ticks5(4))),append('|p| = ',num2str(ticks5(5)))})
    set(findall(hh,'Type','text'),'FontSize',12);
disp('         Fig. 4. Array pattern. ')
disp(' ')
disp(' ')


 
disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. Change the simulation parameters to a 6-element array ')
disp('    (Fig. 4 simulation should change accordingly.)')
disp(' ')
disp(' ')
 

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. After making the changes in Q8,')
disp('    observe Fig. 4, and use the data-tip tool to measure the')
disp('    the amplitude of the main lobe.')
disp('    Also measure the amplitude of the first sidelobe. ')
disp(' ')
disp('Answer: main lobe amplitude= ???  ') 
disp('Answer: first sidelobe amplitude = ???   ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. After making the changes in Q8,')
disp('    measure the FNBW.')
disp('    Also calculate the theoretical FNBW.')
disp(' ')
disp('Answer: measured FNBW = ??? degrees   ')
disp('Answer: theoretical = ??? degrees    ')
disp(' ')
disp(' ')