%%Attempt to solve a multiple scattering problem %Source term, num cylinders %Paramaters for crystal kwav=10; a=.47; %radius of cylinder (m) N=2; pho_scat=7.85; %density of steel (g/cm^3) nmax=4; pho_air=1.29*10^-3; %density of air (g/cm^3) nvals=[-nmax:nmax]; g_i=pho_scat/pho_air; %density ratio cl_scat=4.51e3; %long. sound vel. in steel (m/s) cl_air=343; %long. sound vel. in air (m/s) h_i=cl_scat/cl_air; %long.sound vel. ratio %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 %Set up for the cylinders----source is at rx=0,ry=0 rx=[2*ones(N,1)]; ry=[2*round([1:N]-((N+1)/2)),0]'; %Define Positions %rx=[2,3]'; ry=[0,0]'; ri=sqrt(rx(n).^2+ry(n).^2); %Get Magnitude of ri phi_ri=atan2(ry(n),rx(n)); %Get Angle of ri %phi_ri=atan(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_temp=atan(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; pause; clear('T_i_n','Gamma_i_n_mat','Gamma_i_n','Gijln_temp','Gijln','del_i_j','phi_ri_rj',... 'phi_ri_rj_temp'); %%Assuming I didn't screw up big time revx=[0:0.01:4]; revy=[-3:0.01:3]; 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); i_3rd_dim=zeros(1,1,num_columns); i_3rd_dim(1,1,:)=n; i3d=repmat(i_3rd_dim,[length width 1]); n_3rd_dim=zeros(1,1,num_columns); n_3rd_dim(1,1,:)=l_lr; n3d=repmat(n_3rd_dim,[length width 1]); Ain_3rd_dim=zeros(1,1,num_columns); Ain_3rd_dim(1,1,:)=Ain; Ain_3d=repmat(Ain_3rd_dim,[length width 1]); revX3d=repmat(revX,[1 1 num_columns]); revY3d=repmat(revY,[1 1 num_columns]); r_ri_3d=sqrt((revX3d-rx(i3d)).^2+(revY3d-ry(i3d)).^2); phi_ri_3d=atan2((revY3d-ry(i3d)),(revX3d-rx(i3d))); %phi_ri_3d=atan((revY3d-ry(i3d))./(revX3d-rx(i3d))); [val,uniqueis]=unique(n); inside_cyls=find(r_ri_3d(:,:,uniqueis) <= a); [rcylx,rcyly,icyl]=ind2sub(size(phi_ri_3d),inside_cyls); rcyls=sub2ind(size(revX),rcylx,rcyly); pause clear('icyl','i3d','revY3d','revX3d','revx','revy','revX','revY','rclyx','rcyly','icyl'); term_sum=sum(i*pi.*Ain_3d.*besselh(n3d,kwav*r_ri_3d).*exp(i.*n3d.*phi_ri_3d),3); pho_r=pho_source_r+term_sum; pho_r(rcyls)=-.5; imagesc([minrx,maxrx],[minry,maxry],real(pho_r));