%%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);