%%Attempt to solve a multiple scattering problem with a TRIANGULAR LATTICE %General Options movie=2; %If movie=1, a .avi file is created intensedouble=2; %multiply intensity on right side %Angular Momentum value %Paramaters for crystal nmax=4; nvals=[-nmax:nmax]; a=2; %radius of cylinder (m) pho_scat=7.67; pho_air=1.29*10^-3; %density of steel/air (g/cm^3 g_i=pho_scat/pho_air; %density ratio cl_scat=6010; cl_air=340; %long. sound steel/vel. in air (m/s) symmetric=1; h_i=cl_scat/cl_air; %long.sound vel. ratio %User enters index of refraction nind=input('enter index of refraction:'); if nind==-.7 f=.2; %Filling Fraction lat_a=sqrt(2*pi*a^2/(sqrt(3)*f)); %Lattice Constant xo=5*lat_a; %First x value lambda=1.54*lat_a; %Wavelength kwav=2*pi/lambda; %Wavevector elseif nind==-1 f=.906; lat_a=sqrt(2*pi*a^2/(sqrt(3)*f)); xo=3.5*lat_a; lambda=2.74*lat_a; kwav=2*pi/lambda; end %Set up for the cylinders----source is at rx=0,ry=0 -----lattice constant a %This is a triangular lattice so that every other column is offset by a/2 %I is number of columns J is the number of rows %Total number of cylinders=floor(I/2)*J+ceil((I-1)/2)*(J-1) I=2; Imax=I; J=3; N=ceil(I/2)*J+ceil((I-1)/2)*(J-1); %Values for the columns index=(1:N); col1=ones(J,1); col2=zeros(J-1,1); subcols=cat(1,col1,col2); rx=xo+(sqrt(3)/2)*lat_a*(2*ceil(index/(2*J-1))'-repmat(subcols,I/2,1)); %This will add an extra column to make it symmetric about the middle if symmetric==1 rx=cat(1,rx,(I+1)*sqrt(3)/2*lat_a+xo*col1); end %Values for the Rows; row1=(lat_a)*(-floor(J/2):floor(J/2)); row2=row1(1:J-1)+lat_a/2; subrows=cat(1,row1',row2'); ry=repmat(subrows,I/2,1); %Add extra column y values if symmetric about the middle if symmetric==1 ry=cat(1,ry,row1'); end %If the extra column has been added N will be different, update anyway [N,meaningless]=size(rx); %Boundaries for Field xright=ceil(max(rx)+10*lat_a); ytop=ceil(max(ry)+2*lat_a); ybot=floor(min(ry)-2*lat_a); %Spot to evaluate Amplitude evalx=max(rx)+2*lat_a; evaly=0; %The Matrix will have Columns------ N*(2*nmax+1) % Rows--------- N*(2*nmax+1) num_columns=N*(2*nmax+1); j=[1:num_columns]; k=[1:num_columns]; %Index each column by l and each row by n l_lr=repmat(nvals,1,N); l_ud=l_lr'; n_ud=l_ud; n=ceil(j/(2*nmax+1))'; [J,I]=meshgrid(n,n); [l_ind,n_ind]=meshgrid(l_lr,l_ud); l_minus_n=l_ind-n_ind; del_i_j=I~=J; %Set up the radius vectors and the angles ri=sqrt(rx(n).^2+ry(n).^2); %Get Magnitude of ri phi_ri=atan2(ry(n),rx(n)); %Get Angle of ri ri_rj=sqrt((rx(I)-rx(J)).^2+(ry(I)-ry(J)).^2); %Calculate Mag Differences phi_ri_rj_temp=atan2(del_i_j.*(ry(I)-ry(J)),(rx(I)-rx(J))); %Get angles phi_ri_rj=zeros(num_columns); %Create array for phi vals phi_ri_rj(find(del_i_j==1))=phi_ri_rj_temp(find(del_i_j==1)); %Keep i!=j %Loop through various values of k to get a transmission curve %Arguments of the bessel function num_frequencies=100; Amplitudes=zeros(1,num_frequencies); kas=zeros(1,num_frequencies); reval=sqrt(evalx^2+evaly pho_source_r=i*pi*besselh(0, for knum=1:num_frequencies kwav=1000*knum+1000; ka=kwav.*a.*ones(num_columns,1); kas(knum)=kwav*a; kah=ka/h_i; gh=g_i*h_i*ones(num_columns,1); %Assumes that all cylinders are the same %Here are the coefficients for the matrix multiplication to solve % i i N nmax i,j j i %gamma A - sum sum G A = T % n n j,j~=i l=-nmax l,n l nr Gamma_i_n=(besselh(l_ud,ka).*.5.*(besselj(l_ud-1,kah)-besselj(l_ud+1,kah))-... gh.*.5.*(besselh(l_ud-1,ka)-besselh(l_ud+1,ka)).*besselj(l_ud,kah))./... (gh.*.5.*(besselj(l_ud-1,ka)-besselj(l_ud+1,ka)).*besselj(l_ud,kah)-... besselj(l_ud,ka).*.5.*(besselj(l_ud-1,kah)-besselj(l_ud+1,kah))); Gamma_i_n_mat=diag(Gamma_i_n); T_i_n=besselh(-n_ud,kwav*ri).*exp(-i.*n_ud.*phi_ri); Gijln_temp=besselh(l_minus_n,kwav.*ri_rj).*exp(i.*l_minus_n.*phi_ri_rj); Gijln=zeros(num_columns); Gijln(find(del_i_j==1))=Gijln_temp(find(del_i_j==1)); %Keep i!=j %%Solve for the Coeeficients Alj Ain=(Gamma_i_n_mat-Gijln)^-1*T_i_n; term_sum=0; for in=1: num_columns if mod(in,2*nmax+1)==1 r_ri=sqrt((evalx-rx(n(in))).^2+(evaly-ry(n(in))).^2); phi_ri=zeros(size(r_ri)); phi_ri=atan2((evaly-ry(n(in))),(evalx-rx(n(in)))); end term_sum=term_sum+i*pi.*Ain(in).*besselh(l_lr(in),kwav*r_ri).*... exp(i.*l_lr(in).*phi_ri); end pho_eval=pho_ Amplitudes(knum)=abs(term_sum)^2; end plot(kas,Amplitudes);