%%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; %filling factor nvals=[-nmax:nmax]; a=.1; %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=.47; %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 pause; %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); %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; %Arguments of the bessel function ka=kwav.*a.*ones(num_columns,1); kah=ka/h_i; gh=g_i*h_i*ones(num_columns,1); %Assumes that all cylinders are the same 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 %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; clear('T_i_n','Gamma_i_n_mat','Gamma_i_n','Gijln_temp','Gijln','del_i_j','phi_ri_rj',... 'phi_ri_rj_temp','J','K','j','k'); %%Assuming I didn't screw up big time revx=[-1:0.02:xright]; revy=[ybot:0.02:ytop]; maxrx=max(revx); maxry=max(revy); minrx=min(revx); minry=min(revy); [revX,revY]=meshgrid(revx,revy); rev=(sqrt(revX.^2+revY.^2)); pho_source_r=i*pi*besselh(0,rev); %Create the 3 dimensional array to get this done. Summing Ani terms % i designates num of cyliner n designates order of hankel function [length,width]=size(revX); %Use for loop this time because computer can't handle the memory storage term_sum=zeros(size(revX)); inside_cyls=[]; for in=1: num_columns if mod(in,2*nmax+1)==1 r_ri=sqrt((revX-rx(n(in))).^2+(revY-ry(n(in))).^2); phi_ri=zeros(size(r_ri)); %phi_ri=atan2((revY-ry(n(in))),(revX-rx(n(in)))); %Don't care about values inside cyls inside_cyls=cat(1,inside_cyls,find(r_ri <=a)); outside_cyls=find(r_ri >a); phi_ri(outside_cyls)=atan2((revY(outside_cyls)-ry(n(in))),... (revX(outside_cyls)-rx(n(in)))); end term_sum(outside_cyls)=term_sum(outside_cyls)+... i*pi.*Ain(in).*besselh(l_lr(in),kwav*r_ri(outside_cyls)).*... exp(i.*l_lr(in).*phi_ri(outside_cyls)); end pho_r=pho_source_r+term_sum; pho_r(inside_cyls)=-.5; %For pretty picture purposes---make length x width x 3 true color image inside_cyls=repmat(inside_cyls,[1 1 3]); pho_r_scaled=(real(pho_r)-min(real(pho_r(:))))/(max(real(pho_r(:)))-min(real(pho_r(:)))); pho_r3_scaled=repmat(pho_r_scaled,[1 1 3]); pho_r3_scaled(inside_cyls(:,:,1))=1; pho_r3_scaled(inside_cyls(:,:,2))=0; pho_r3_scaled(inside_cyls(:,:,3))=0; imagesc([minrx,maxrx],[minry,maxry],real(pho_r3_scaled)); pause %Attempt at a movie? if movie==1 avistring=['/nfs/slac/g/ki/ki08/lsst/CPanalysis/Phononic_Crys/','MST_Tri_Lat_lat_a',num2str(lat_a),'_f_',num2str(f),... '_lambda_',num2str(lambda),'.avi']; aviobj=avifile(avistring); end tmax=80; omega=.65; set(gcf,'doublebuffer','on'); for t=0:tmax pho_movie=exp(-i*omega*t).*pho_r; pho_movie(inside_cyls)=-.5; if intensedouble==1 pho_movie(find(revX > (Imax+1)*sqrt(3)/2*lat_a+xo+a))=... 10*pho_movie(find(revX > (Imax+1)*sqrt(3)/2*lat_a+xo+a)); end imagesc([minrx,maxrx],[minry,maxry],real(pho_movie),[-15 15]);colormap(gray); if movie==1 frame=getframe(gca); for frames=1:10 aviobj=addframe(aviobj,frame); end end pause(0.1); end if movie==1 aviobj=close(aviobj); end