%%Attempt to solve a multiple scattering problem 

%Source term, num cylinders                 %Paramaters for crystal
kwav=10;                                    a=.2;                   %radius of cylinder (m)
N=5;                                        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

%Set up for the cylinders----source is at rx=0,ry=0
rows=[1,2];                      cols=[1,1];
rx=1.5*ones(N,1)*rows;           ry=round([1:N]'-((N+1)/2))*cols;          %Define Positions
[irx,jrx]=size(rx);              [iry,jry]=size(ry);
rx=reshape(rx,irx*jrx,1);        ry=reshape(ry,iry*jry,1);
N=irx*jrx;                                                                 %redefine cylinders
clear('irx','jrx','iry','jry','i','j');

%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=[-2:0.01:6];          revy=[-4:0.01:4];
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
    phi_ri=atan2((revY-ry(n(in))),(revX-rx(n(in))));
    r_ri=sqrt((revX-rx(n(in))).^2+(revY-ry(n(in))).^2);
    if mod(in,2*nmax+1)==1
        inside_cyls=cat(1,inside_cyls,find(r_ri <=a));
    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_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));

%Attempt at a movie?
tmax=100;       omega=.1;

for t=0:tmax
    pho_movie=exp(-i*omega*t).*pho_r;
    pho_movie(inside_cyls)=-.5;
    imagesc([minrx,maxrx],[minry,maxry],real(pho_r3_scale));
end