There are all sorts of resonances around us, in the world, in our culture, and in our technology. A tidal resonance causes the 55 foot tides in the Bay of Fundy. Mechanical and acoustical resonances and their control are at the center of practically every musical instrument that ever existed. Even our voices and speech are based on controlling the resonances in our throat and mouth. Technology is also a heavy user of resonance. All clocks, radios, televisions, and gps navigating systems use electronic resonators at their very core. Doctors use magnetic resonance imaging or MRI to sense the resonances in atomic nuclei to map the insides of their patients. In spite of the great diversity of resonators, they all share many common properties. In this blog, we will delve into their various aspects. It is hoped that this will serve both the students and professionals who would like to understand more about resonators. I hope all will enjoy the animations.

For a list of all topics discussed, scroll down to the very bottom of the blog, or click here.

Origins of Newton's laws of motion

Non-mathematical introduction to relativity

Three types of waves: traveling waves, standing waves and rotating waves new

History of mechanical clocks with animations
Understanding a mechanical clock with animations
includes pendulum, balance wheel, and quartz clocks

Water waves, Fourier analysis


Saturday, October 21, 2017

Octave programs used in Emitter Games

Octave programs used in Emitter Games

Sample of some of the images in the previous posting



The Octave programs listed below were used to make the images and movies shown in the previous posting. Some of the lower postings are quite lengthy because they were used to make several types of images of the same wavefield. That is to say, a specific output command was active while the output commands were commented out for each specific image type.

It is recommended that the viewer read the previous posting before studying this posting. The methods and math are discussed in that posting.



Program for snap shots of emitter waves 
% snapShotSimpleEmitter.m
clear;
waveLen=2;
k=2*pi/waveLen;
A=1;% strength of central emitter
Nx=200;% number of points in x direction
Ny=100;% number of points in y direction
x=linspace(-10,10,Nx);
y=linspace(-5,5,Ny);
[xx,yy]=meshgrid(x,y);
rr=sqrt( xx.^2 + yy.^2 );
pphi=atan2(yy, xx);
multiply=sqrt(2*A./(pi*k*rr));
Phi=multiply.*(besselh(1,1,k*rr).*cos(pphi));

pcolor(xx(:,:),yy(:,:),real(Phi))
shading interp
%colorbar
axis("equal","tight")
cmax=0.1;
caxis([0,cmax]);
This program produces a snap shot (a 2D color mapping) of the wavefield around a dipole emitter as shown in Fig. 4a of Emitter games.

By changing the "Phi=multiply.*..." line into "Phi=multiply.*besselh(0,1,k*rr).*cos(pphi));" we will have the field around a monopole.

By changing the same line to
"Phi=multiply.*(besselh(0,1,k*rr) + I*besselh(1,1,k*rr).*cos(pphi));" we will have the field around a directional emitter.

By changing the next line, "pcolor(...", to "pcolor(xx(:,:),yy(:,:),abs(Phi));" we will show the amplitude of the wave field everywhere such as in Fig. 4b of Emitter games.




Program for making movies of wave fields 
% framesSimpleEmitter.m
clear;
n=0;% number of emitters in addition 
  %         to center one. Can be 0.
Nemitt=2*n+1;
waveLen=2;
k=2*pi/waveLen;
a=0.5*waveLen;% spacing between emitters
Ao=1;% strength of central emitter
An=1; % strength of other emitters
Nx=120;% number of points in x direction
Ny=40;% number of points in y direction
x=linspace(-5,15,Nx);
y=linspace(-10,10,Ny);
if (n==0) nM=0;yEmit=0;AA=Ao;
  else nM=-n:n; 
  AA=An*ones(1,1,Nemitt); 
  AA(1,1,n+1)=Ao; 
  endif
AA(1,1,:);
yEmitt=nM*a;
[xx,yy,yyEmitt]=meshgrid(x,y,yEmitt);
rr=sqrt( xx.^2 + (yy-yyEmitt).^2 );
pphi=atan2( (yy-yyEmitt), xx);
multiply1=(sqrt(2*AA./(pi*k*rr)))/sqrt(Nemitt);
Phi=sum(multiply1.*(besselh(0,1,k*rr)\
     +I*(besselh(1,1,k*rr)).*cos(pphi)),3);

iN=99;% normally set to 99, 
% total count as if all frames are to be printed
iStart=1;%iN; normally set this to 1
iend=iN; % iN normally set this to iN
maxRadians=12*pi;% set range of wt in 
%            radians, 12pi is good
dwt=maxRadians/iN;
wt=dwt*(iStart-1);

for i=iStart:iend%[12,33,38,56,85,86,91]   
  wt=i*dwt;
  PhiOut=real(Phi).*cos(wt)+imag(Phi).*sin(wt);
  pcolor(xx,yy,PhiOut);
  shading interp;
  cmax=0.1;
  caxis([-cmax cmax]);
  grid
  %colorbar("EastOutside");
  axis('equal','tight');
  pause(0.4);
  stringNm=sprintf('dirEmSingle/wav10%02d.png',i);
  %print(stringNm,"-dpng","-S335,335","-F:5");
endfor
This program produces frames for a movie of the waves around a directional emitter as shown in Fig. 5 of Emitter games.

This program takes a considerable time to produce all 99 frames, with the bulk of the time in the "print" statement at the end. For that reason is recommended that the program be fine tuned with the print statement commented out. Also, an appropriate subdirector should be created to receive the frames. I found that some of the frames were improperly created and needed to be recreated. For that I changed iStart and iend both to equal the defective frames index.

To make the movie from the frames I used the command line program ffmpeg for my linux mint operating system. The actual command was:
ffmpeg -i wav10%2d.png movieName.mp4 Excution of this was almost instantaneous. I then uploaded the movie to youtube and linked it to the html document.




Program for making 2D color maps of linear arrays of emitters 
%  linearArrayOfEmitter2Dmap.m
clear;
n=5;% number of emitters in addition to center one. Can be 0.
Nemitt=2*n+1
waveLen=2;
k=2*pi/waveLen;
a=0.25*waveLen;% spacing between emitters
Ao=1;% strength of central emitter
An=1; % strength of other emitters
Nx=200;% number of points in x direction
Ny=200;% number of points in y direction
x=linspace(-5,15,Nx);
y=linspace(-10,10,Ny);
if (n==0) nM=0;yEmit=0;AA=Ao;
else nM=-n:n; AA=An*ones(1,1,Nemitt); AA(1,1,n+1)=Ao; endif
AA(1,1,:);
yEmitt=nM*a;
[xx,yy,yyEmitt]=meshgrid(x,y,yEmitt);
krr=k*sqrt( xx.^2 + (yy-yyEmitt).^2 );
pphi=atan2( (yy-yyEmitt), xx);
multiply1=AA.*sqrt(2./pi)/sqrt(Nemitt);
Phi=0.5*sum(multiply1.*(besselh(0,1,krr)+I*(besselh(1,1,krr)).*cos(pphi)),3);
%Phi=sum(multiply1.*(besselh(0,1,krr)),3);
%Phi=sum(multiply1.*(besselh(1,1,krr)).*cos(pphi),3);

subplot (1,2,1);
pcolor(xx(:,:,1),yy(:,:,1),real(Phi))
shading interp
%colorbar
axis("equal","tight")
cmax=0.4;
caxis([0,cmax]);

subplot (1,2,2);
pcolor(xx(:,:,1),yy(:,:,1),abs(Phi))
shading interp
%colorbar
axis("equal","tight")
cmax=0.4;
caxis([0,cmax]);
This program creates the false color mappings of wave fields from linear arrays of emitters as shown in Figs. 7, 8 and 9.




Program for making Fraunhofer diffraction patterns 
%  FraunhoferDiffractionGraph.m
clear;
n=20;% number of emitters in addition
%          to center one. Can be 0.
Nemitt=2*n+1;
waveLen=4;
k=2*pi/waveLen;
a=5*waveLen% slit width
spacing=a/(Nemitt-1)% hide this statement if n=0
Ao=1;% strength of central emitter
An=1; % strength of other emitters
Nx=2;% number of points in x direction
Ny=400;% number of points in y direction
x=linspace(500,501,Nx);
ymax=400;
y=linspace(-ymax,ymax,Ny);
if (n==0) nM=0;
  yEmit=0;AA=Ao;
  else nM=-n:n; 
  AA=An*ones(1,1,Nemitt); 
  AA(1,1,n+1)=Ao; 
endif
AA(1,1,:);
yEmitt=nM*spacing;
[xx,yy,yyEmitt]=meshgrid(x,y,yEmitt);
krr=k*sqrt( xx.^2 + (yy-yyEmitt).^2 );
pphi=atan2( (yy-yyEmitt), xx);
multiply1=AA.*sqrt(2./pi)/sqrt(Nemitt);
Phi3=0.5*sum(multiply1.*(besselh(0,1,krr)...
   +I*(besselh(1,1,krr)).*cos(pphi)),3);
Phi1=sum(multiply1.*(besselh(0,1,krr)),3);
Phi2=sum(multiply1.*(besselh(1,1,krr)).*cos(pphi),3);
PhiFraun1=(0.0025)*sum(exp(-I*k*(yy(:,1,:).* ...
    yyEmitt(:,1,:))./(x(1,1))),3);
xFraun=k*a*y/(2*500*pi); % the pi is because
%      the octave sinc function has a pi in it
PhiFraun2=0.1*abs(sinc(xFraun));
xFraun2=k*a*y./(2*pi*sqrt(500^2+y.^2));
PhiFraun3=0.1*abs(sinc(xFraun2));
delP=0.006;
% subplot(1,1,1);
plot(y,abs(Phi1(:,1)),y,abs(Phi2(:,1))...
     +delP,y,abs(Phi3(:,1))+2*delP,...
     y,abs(PhiFraun3)+3*delP,y,PhiFraun2+4*delP,...
     y,abs(PhiFraun1)+5*delP);
grid
This program creates the graph of diffraction patterns shown in Fig. 12 of "Emitter games". This represents a small modification from the program just above used to make Figs. 7, 8 and 9.




Program for making reflections via emitters 
%   mapping2DwithMirror.m

clear;
n=30;% number of emitters in addition 
%      to center one. 1 or larger
Nemitt=2*n+1;
waveLen=3;
xmax=50;
xmin=-5;
ymax=20;
ymin=-40;
k=2*pi/waveLen;
a=20;% slit width of all emitters
Nx=200;% number of points in x direction
Ny=260;% number of points in y direction

x=linspace(xmin,xmax,Nx);
y=linspace(ymin,ymax,Ny);
dyEm=a/(Nemitt-1);
nEm0=1:Nemitt;
nEm(1,1,:)=nEm0;
yEm= -(a/2)+(nEm-1).*dyEm;
xEm=ones(1,1,Nemitt)+5;
emDen=(Nemitt-1)/a;
[xx1,yy1,nnEm]=meshgrid(x,y,nEm0);
rEm=k*sqrt((xx1-xEm).^2+(yy1-yEm).^2);

%Next do reflector
Nmir=60;
xStartMirror=16;
slope=-1.;
widthMir=1.5*a;
totMirWidth=widthMir/cos(atan(slope));
mirEmDensity=(Nmir-1)/totMirWidth;
inverseMirDensity=1/mirEmDensity;
dxMir=widthMir/(Nmir-1);
nMir0=1:Nmir;  % index for mirror counter emitters
nMir(1,1,:)=nMir0;
xMir=xStartMirror+dxMir*(nMir-1);
xMiddle=xStartMirror+(widthMir/2);
yMir=(xMir-xMiddle)*slope;
[xx2,yy2,nnMir]=meshgrid(x,y,nMir0);
rMir=k*sqrt((xx2-xMir).^2+(yy2-yMir).^2);
xMir1=shiftdim(xMir);
yMir1=shiftdim(yMir);
xEm1=(shiftdim(xEm))';
yEm1=(shiftdim(yEm))';

rEmMir=k*sqrt((xMir1-xEm1).^2 + (yMir1-yEm1).^2);
Phi12 = inverseMirDensity*shiftdim(sum(besselh...
    (0,1,rEmMir).'),-1);
Phi13 = -((2/pi))*sum((Phi12.*besselh(0,1,rMir)),3);

Phi1=sum(besselh(0,1,rEm),3);
Phi3=Phi1+Phi13;
PPhi1=real(Phi3);
PPhi2=imag(Phi3);

pcolor(xx2(:,:,1),yy2(:,:,1),PPhi2)
shading interp
%colorbar
axis("equal","tight")
cmax=5;
caxis([-cmax cmax]);
%print("pics/mirrorStraight.png","-dpng",
%    "-S480,233","-F:5")
 This Octave program created the images of reflections via emitters shown in Figs. 14 and 15 of "Emitter games".



Program for modeling a circular resonant cavity 
% resCavityMethod2_findingResonance_16Sept17.m
%  expected Jo resonances at 
%      wavelen= 2.6a, 1.14a, 0.726a

clear;
xSource=0;
ySource=0;
radiusCavity=1;
radiusSample=1.4*radiusCavity;
nCavity=15;
nSample=nCavity;
dphiCavity=2*pi/nCavity;
dphiSample=2*pi/nSample;
iCavity=1:nCavity;
iSample=1:nSample;
phiCavity=(iCavity-1)*dphiCavity;
phiSample=(iSample-1)*dphiSample;
xCavity=radiusCavity*cos(phiCavity);
yCavity=radiusCavity*sin(phiCavity);
xSample=radiusSample*cos(phiSample);
ySample=radiusSample*sin(phiSample);
rCavitySample=sqrt((xCavity-xSample.')
      .^2+(yCavity-ySample.').^2);
rSourceSample=sqrt((xSource-xSample.')
       .^2+(ySource-ySample.').^2);

% SETUP MATRIX
nXMatrix=100;
nYMatrix=85;
xmax=1.5*radiusCavity;
xmin=-xmax;
ymax=xmax;
ymin=-ymax;
xMatrix=linspace(xmin,xmax,nXMatrix);
yMatrix=linspace(ymin,ymax,nYMatrix);
xSource1=shiftdim(xSource,-1);
ySource1=shiftdim(ySource,-1);
rSourceMatrix=sqrt((xSource1.-xMatrix')
        .^2 + (ySource1.-yMatrix).^2);
xCavity1=shiftdim(xCavity,-1);
   % makes a 1x1xN array where cavity variable is last
yCavity1=shiftdim(yCavity,-1);
rCavityMatrix=sqrt((xCavity1-xMatrix').^2
        + (yCavity1-yMatrix).^2);

minWavelen=radiusCavity*2.611;
     %1.13;%2.5; resonance at 1.1375
maxWavelen=radiusCavity*2.616;
      %1.145;%2.7  resonance at 2.6136
nWavelen=40;
delWavelen=(maxWavelen-minWavelen)/(nWavelen-1);
iWavelen=1:nWavelen;
wavelen1=(iWavelen-1)*delWavelen + minWavelen;
fieldSliceX=zeros(1,nWavelen);
response=zeros(1,nWavelen);
iMiddleY=floor(nYMatrix/2);
iMiddleX=floor(nXMatrix/2);
[xx,yy]=meshgrid(xMatrix,yMatrix);

for i=1:nWavelen
  wavelen=(i-1)*delWavelen + minWavelen;
  k1=2*pi/wavelen;
  MCavitySample1=(besselh(0,1,k1*rCavitySample)).';
  S1 = -besselh(0,1,k1*rSourceSample);
  Phi1=linsolve(MCavitySample1,S1,struct ());
  Phi2=shiftdim(Phi1,-2);
  fieldCavity1=sum(Phi2.*besselh(0,1,k1*rCavityMatrix),3);
  fieldSource1=besselh(0,1,k1*rSourceMatrix); 
  totalField1=fieldCavity1+fieldSource1;
  pause(0.1);
  response(1,i)=totalField1(iMiddleY,iMiddleX);
  pcolor(xx,yy,abs(totalField1'));
  shading interp
  caxis([0 200]);%1000,  100
 % colorbar
 % print("pics/secondResonance2D.png","-dpng",
 %   "-S380,250","-F:5");
endfor 


plot(wavelen1,abs(response));
xlabel("wavelength");
axis([minWavelen maxWavelen 0 800]);% 70
print("pics/firstResonance.png","-dpng",
    "-S320,233","-F:5");
 This Octave program created the 2D color mapping and resonance curves for the emulated circular resonant cavities in Figs. 17 and 18 of Emitter games.



Program for finding resonance for circular cavity 
% resCavityMethod2_findingResonance_16Sept17.m
%  expected Jo resonances at wavelen= 2.6a, 1.14a, 0.726a

clear;
xSource=0;
ySource=0;
radiusCavity=1;
radiusSample=1.4*radiusCavity;
nCavity=15;
nSample=nCavity;
dphiCavity=2*pi/nCavity;
dphiSample=2*pi/nSample;
iCavity=1:nCavity;
iSample=1:nSample;
phiCavity=(iCavity-1)*dphiCavity;
phiSample=(iSample-1)*dphiSample;% add a little offset from cavity points
xCavity=radiusCavity*cos(phiCavity);
yCavity=radiusCavity*sin(phiCavity);
xSample=radiusSample*cos(phiSample);
ySample=radiusSample*sin(phiSample);
rCavitySample=sqrt((xCavity-xSample.').^2+(yCavity-ySample.').^2);
rSourceSample=sqrt((xSource-xSample.').^2+(ySource-ySample.').^2);

% SETUP MATRIX
nXMatrix=100;
nYMatrix=85;
xmax=1.5*radiusCavity;
xmin=-xmax;
ymax=xmax;
ymin=-ymax;
xMatrix=linspace(xmin,xmax,nXMatrix);
yMatrix=linspace(ymin,ymax,nYMatrix);
xSource1=shiftdim(xSource,-1);
ySource1=shiftdim(ySource,-1);

rSourceMatrix=sqrt((xSource1.-xMatrix').^2 + (ySource1.-yMatrix).^2);
xCavity1=shiftdim(xCavity,-1);% makes a 1x1xN array where cavity variable is last
yCavity1=shiftdim(yCavity,-1);
rCavityMatrix=sqrt((xCavity1-xMatrix').^2 + (yCavity1-yMatrix).^2);

minWavelen=radiusCavity*2.611;%1.13;%2.5; resonance at 1.1375
maxWavelen=radiusCavity*2.616;%1.145;%2.7  resonance at 2.6136
nWavelen=40;
delWavelen=(maxWavelen-minWavelen)/(nWavelen-1);
iWavelen=1:nWavelen;
wavelen1=(iWavelen-1)*delWavelen + minWavelen;
fieldSliceX=zeros(1,nWavelen);
response=zeros(1,nWavelen);
iMiddleY=floor(nYMatrix/2);
iMiddleX=floor(nXMatrix/2);
[xx,yy]=meshgrid(xMatrix,yMatrix);

for i=1:nWavelen
  wavelen=(i-1)*delWavelen + minWavelen;
  k1=2*pi/wavelen;
  MCavitySample1=(besselh(0,1,k1*rCavitySample)).';
  S1 = -besselh(0,1,k1*rSourceSample);
  Phi1=linsolve(MCavitySample1,S1,struct ());
  Phi2=shiftdim(Phi1,-2);
  fieldCavity1=sum(Phi2.*besselh(0,1,k1*rCavityMatrix),3);
  fieldSource1=besselh(0,1,k1*rSourceMatrix); 
  totalField1=fieldCavity1+fieldSource1;
  pause(0.1);
  response(1,i)=totalField1(iMiddleY,iMiddleX);
endfor 

plot(wavelen1,abs(response));
xlabel("wavelength");
axis([minWavelen maxWavelen 0 800]);% 70
grid
%print("pics/firstResonanceCircular.png","-dpng",
%    "-S320,233","-F:5");



Octave program for making a movie of a rotating wave mode in a resonant cavity 
% % resonantCavityRotating_21Sept17.m
clear;
radiusCavity=1;
xSource11=0.5*radiusCavity;
ySource11=0;
xSource22=0.35*radiusCavity;
ySource22=0.35*radiusCavity;
radiusSample=1.4*radiusCavity;
nCavity=22;
nSample=nCavity;
dphiCavity=2*pi/nCavity;
dphiSample=2*pi/nSample;
iCavity=1:nCavity;
iSample=1:nSample;
phiCavity=(iCavity-1)*dphiCavity;
phiSample=(iSample-1)*dphiSample;
xCavity=radiusCavity*cos(phiCavity);
yCavity=radiusCavity*sin(phiCavity);
xSample=radiusSample*cos(phiSample);
ySample=radiusSample*sin(phiSample);
rCavitySample=sqrt((xCavity-xSample.')
      .^2+(yCavity-ySample.').^2);
rSource1Sample=sqrt((xSource11-xSample.')
       .^2+(ySource11-ySample.').^2);
rSource2Sample=sqrt((xSource22-xSample.')
       .^2+(ySource22-ySample.').^2);

% SETUP MATRIX
nXMatrix=100;
nYMatrix=85;
xmax=1.2*radiusCavity;
xmin=-xmax;
ymax=xmax;
ymin=-ymax;
xMatrix=linspace(xmin,xmax,nXMatrix);
yMatrix=linspace(ymin,ymax,nYMatrix);
xSource1=shiftdim(xSource11,-1);
ySource1=shiftdim(ySource11,-1);
rSource1Matrix=sqrt((xSource11.-xMatrix')
        .^2 + (ySource11.-yMatrix).^2);
rSource2Matrix=sqrt((xSource22.-xMatrix')
        .^2 + (ySource22.-yMatrix).^2);
xCavity1=shiftdim(xCavity,-1);
   % makes a 1x1xN array where 
   %cavity variable is last
yCavity1=shiftdim(yCavity,-1);
rCavityMatrix=sqrt((xCavity1-xMatrix').^2
        + (yCavity1-yMatrix).^2);

[xx,yy]=meshgrid(xMatrix,yMatrix);
wavelen = 1.223*radiusCavity;
k1=2*pi/wavelen;
MCavitySample1=(besselh(0,1,k1*rCavitySample)).';
S1 = -besselh(0,1,k1*rSource1Sample);
S2 = -besselh(0,1,k1*rSource2Sample);
Phi1a=linsolve(MCavitySample1,S1,struct ());
Phi1b=linsolve(MCavitySample1,S2,struct ());
Phi2a=shiftdim(Phi1a,-2);
Phi2b=shiftdim(Phi1b,-2);
fieldCavity1a=sum(Phi2a.*besselh(0,1,
    k1*rCavityMatrix),3);
fieldCavity1b=sum(Phi2b.*besselh(0,1,
    k1*rCavityMatrix),3);
fieldSource1a=besselh(0,1,k1*rSource1Matrix);
fieldSource1b=besselh(0,1,k1*rSource2Matrix);
totalField1a=fieldCavity1a+fieldSource1a;
totalField1b=fieldCavity1b+fieldSource1b;
totalField=totalField1a + I*totalField1b;
 
iN=99;% normally set to 99, total count as 
%   if all frames are to be printed
iStart=1;%iN; normally set this to 1
iend=5; % iN normally set this to iN
maxRadians=12*pi;% set range of wt in radians,
%                    12pi is good
dwt=maxRadians/iN;
wt=dwt*(iStart-1);

for i=iStart:iend%   
  wt=i*dwt;
  PhiOut=real(totalField*exp(-I*wt));
  %PhiOut=real(Phi).*cos(wt)+imag(Phi).*sin(wt);
  pcolor(xx,yy,PhiOut.');
  shading interp;
  cmax=270;
  caxis([-cmax cmax]);
  grid
  %colorbar("EastOutside");
  axis('equal','tight');
  pause(0.4);
  stringNm=sprintf(
        'cavityRotatingMovie/rot%02d.png',i);
%  print(stringNm,"-dpng","-S330,261","-F:5"); 
endfor
This Octave program created the 2D color mapping of a circular cavity excited in the (2,1) rotating mode. The output appeared in "Emitter games".



Program for showing the spectra response of an L-shaped resonator. It can also be used for plotting the layout of the cavity, emitter points, exciter, and sample points. 
% LresonatorUniformlySpacedFindRes_7Oct17.m

% make wall emitter arrayfun
clear;
xVertex=[0,0,-10,-10,10,10];
yVertex=[0,20,20,-10,-10,0];
n=columns(xVertex);
sideLengths=zeros(1,n);
XlengthEachSide(1,1:n-1)=xVertex(1,2:n)-xVertex(1,1:n-1);
XlengthEachSide(1,n)=xVertex(1,1)-xVertex(1,n);
YlengthEachSide(1,1:n-1)=yVertex(1,2:n)-yVertex(1,1:n-1);
YlengthEachSide(1,n)=yVertex(1,1)-yVertex(1,n);
sideLengths=sqrt((XlengthEachSide).^2+(YlengthEachSide).^2);
totalLength=sum(sideLengths,2); 
approxPoints=30;
ptsPerLength=approxPoints/totalLength;
pointsPerEachSide=round(sideLengths*ptsPerLength); 
totalPoints=sum(pointsPerEachSide,2); 
i=1:totalPoints;
delXeachSide=XlengthEachSide./pointsPerEachSide;
delYeachSide=YlengthEachSide./pointsPerEachSide;
xCavity=zeros(1,totalPoints);
yCavity=zeros(1,totalPoints);

lineNo=zeros(1,n);
indexi=zeros(1,totalPoints);
indexj=zeros(1,totalPoints);
icounts=1;
for i4=1:n % indexj counts 0 to emitters in each side
           %  indexi has a bunch of 1's then a bunch of 2's, etc.
   for j4=1:pointsPerEachSide(1,i4)
      indexj(icounts)=j4;
      indexi(icounts)=i4;
      icounts += 1;
   endfor
endfor
xCavity(1,i)=xVertex(1,indexi)+(indexj).*delXeachSide(1,indexi);
yCavity(1,i)=yVertex(1,indexi)+(indexj).*delYeachSide(1,indexi);

gap=totalLength/totalPoints;
mult=1.4;
xSample(1,i)=(xCavity<0).*(1.4.*xCavity(1,i)+2)+ ...
    (xCavity>=0).*(1..*xCavity(1,i)+2) ;
ySample(1,i)=(yCavity<0).*(1.3*yCavity(1,i)+1)+ ...
     (yCavity>=0).*(1.01.*yCavity(1,i)+2) ;
   
xSource= -5;
ySource=5;

rCavitySample=sqrt((xCavity-xSample.').^2+(yCavity-ySample.').^2);
rSourceSample=sqrt((xSource-xSample.').^2+(ySource-ySample.').^2);

xSource1=shiftdim(xSource,-1);
ySource1=shiftdim(ySource,-1);   

xCavity1=shiftdim(xCavity,-1);
      % makes a 1x1xN array where cavity variable is last
yCavity1=shiftdim(yCavity,-1);

xResponse=0;
yResponse=-5;
rResponseCavity=sqrt((xResponse-xCavity).^2 + (yResponse-yCavity).^2);
rResponseSource=sqrt((xResponse-xSource).^2 + (yResponse-ySource).^2);

minWaveLen=12.; % resonances at 20.065, 17.862,  15.643
maxWaveLen=25;
NwaveLen=1000;
wavelen=linspace(minWaveLen,maxWaveLen,NwaveLen);% n=2 res at 17, 
waveLresponse=zeros(1,NwaveLen);
response=zeros(1,NwaveLen);
response1=zeros(1,NwaveLen);
response2=zeros(1,NwaveLen);

for (iwave=1:NwaveLen)

k1=2*pi./wavelen(iwave);
krCavitySample=k1.*rCavitySample;
MCavitySample1=besselh(0,1,krCavitySample);
S1 = -besselh(0,1,k1.*rSourceSample);
Phi1=linsolve(MCavitySample1,S1,struct ());
response(1,iwave)=sum((abs(Phi1)).');
response1(1,i)=median(abs(Phi1));
resp1=sum((besselh(0,1,k1*rResponseCavity)).*Phi1.');
response2(1,i)= resp1 + (besselh(0,1,k1*rResponseSource));

endfor

plot(xVertex,yVertex,"-o",xCavity,yCavity,"x",xSample,ySample,".");
axis([-15,20 -15,35],"equal","tight")

semilogy(wavelen,abs(response));
xlabel("wavelength");
axis("tight");
grid

%print("pics/LresonatorSpectra11.png","-dpng",
 %  "-S380,250","-F:5");



Program for showing the fields in a resonance of an L-shaped cavity 
% LresonatorUniformlySpaced_3Oct17.m

% make wall emitter arrayfun
clear;
xVertex=[0,0,-10,-10,10,10];
yVertex=[0,20,20,-10,-10,0];
n=columns(xVertex);
sideLengths=zeros(1,n);
XlengthEachSide(1,1:n-1)=xVertex(1,2:n)-xVertex(1,1:n-1);
XlengthEachSide(1,n)=xVertex(1,1)-xVertex(1,n);
YlengthEachSide(1,1:n-1)=yVertex(1,2:n)-yVertex(1,1:n-1);
YlengthEachSide(1,n)=yVertex(1,1)-yVertex(1,n);
sideLengths=sqrt((XlengthEachSide).^2+(YlengthEachSide).^2);
totalLength=sum(sideLengths,2); 
approxPoints=30;
ptsPerLength=approxPoints/totalLength;
pointsPerEachSide=round(sideLengths*ptsPerLength); 
totalPoints=sum(pointsPerEachSide,2); 
i=1:totalPoints;
delXeachSide=XlengthEachSide./pointsPerEachSide;
delYeachSide=YlengthEachSide./pointsPerEachSide;
xCavity=zeros(1,totalPoints);
yCavity=zeros(1,totalPoints);

lineNo=zeros(1,n);
indexi=zeros(1,totalPoints);
indexj=zeros(1,totalPoints);
icounts=1;
for i4=1:n % indexj counts 0 to emitters in each side
           %  indexi has a bunch of 1's then a bunch of 2's, etc.
   for j4=1:pointsPerEachSide(1,i4)
      indexj(icounts)=j4;
      indexi(icounts)=i4;
      icounts += 1;
   endfor
endfor
xCavity(1,i)=xVertex(1,indexi)+(indexj).*delXeachSide(1,indexi);
yCavity(1,i)=yVertex(1,indexi)+(indexj).*delYeachSide(1,indexi);

gap=totalLength/totalPoints;
mult=1.4;
xSample(1,i)=(xCavity<0).*(1.4.*xCavity(1,i)+2)+ ...
    (xCavity>=0).*(1..*xCavity(1,i)+2) ;
ySample(1,i)=(yCavity<0).*(1.3*yCavity(1,i)+1)+ ...
     (yCavity>=0).*(1.01.*yCavity(1,i)+2) ;
   
xSource= -5;
ySource=5;

rCavitySample=sqrt((xCavity-xSample.').^2+(yCavity-ySample.').^2);
rSourceSample=sqrt((xSource-xSample.').^2+(ySource-ySample.').^2);

xSource1=shiftdim(xSource,-1);
ySource1=shiftdim(ySource,-1);   

xCavity1=shiftdim(xCavity,-1);
      % makes a 1x1xN array where cavity variable is last
yCavity1=shiftdim(yCavity,-1);

xResponse=0;
yResponse=-5;
rResponseCavity=sqrt((xResponse-xCavity).^2 + (yResponse-yCavity).^2);
rResponseSource=sqrt((xResponse-xSource).^2 + (yResponse-ySource).^2);

wavelen=17.862;   %20.065, 17.862, 15.643

k1=2*pi/wavelen;
MCavitySample1=(besselh(0,1,k1*rCavitySample));
S1 = -besselh(0,1,k1*rSourceSample);
Phi1=linsolve(MCavitySample1,S1,struct ());
Phi2=shiftdim(Phi1,-2);
response(1,i)=sum((abs(Phi1)).');
response1(1,i)=median(abs(Phi1));
resp1=sum((besselh(0,1,k1*rResponseCavity)).*Phi1.');
response2(1,i)= resp1 + (besselh(0,1,k1*rResponseSource));

%  Below shows 2D color mapping of the above solution

nXMatrix=50;
nYMatrix=60;
hSide=20;
xmax=15;%15.325;
xmin=-15;%-15.118;
ymax=hSide;%25.492;
ymin=-15;%-15.124;
xMatrix=linspace(xmin,xmax,nXMatrix);
yMatrix=linspace(ymin,ymax,nYMatrix);

rSourceMatrixx=sqrt((xSource1.-xMatrix').^2 + (ySource1.-yMatrix).^2);
rSourceMatrix=(rSourceMatrixx<1).*1 + (rSourceMatrixx>=1).*rSourceMatrixx;
xCavity1=shiftdim(xCavity,-1);% makes a 1x1xN array where cavity variable is last
yCavity1=shiftdim(yCavity,-1);
rCavityMatrixx=sqrt((xCavity1-xMatrix').^2 + (yCavity1-yMatrix).^2);
minDist=.00001;
rCavityMatrix=(rCavityMatrixx=minDist).*rCavityMatrixx;
[xx,yy]=meshgrid(xMatrix,yMatrix);


fieldCavity1=sum(Phi2.*besselh(0,1,k1*rCavityMatrix),3);
fieldSource1=besselh(0,1,k1*rSourceMatrix); 
totalField1=fieldCavity1+fieldSource1;


pcolor(xx,yy,real(totalField1'));
shading interp
%caxis([-700 1100]);%0,  100
colorbar

%print("pics/firstLresonance.png","-dpng",
%   "-S380,250","-F:5");


Posted by P. Ceperley at 12:43 PM