%%Attempt to solve a multiple scattering problem with a TRIANGULAR LATTICE

%General Options
movie=2;		           %If movie=1, a .avi file is created
intensedouble=1;		   %multiply intensity on right side 
yoffset=0;		           %Move offcenter=1
posneg=2;			       %1=peaks/troughs  2=weak/strong
lattice=2;                 %1=triangle       2=square
if lattice==1
    lattice_str='Tri';
else
    lattice_str='Square';
end

%Angular Momentum value                    %Paramaters for crystal
nmax=2;                                     
nvals=[-nmax:nmax];                         
a=.1;               			    %radius of cylinder (m)

 %densities of steel/air (g/cm^3
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)
h_i=cl_scat/cl_air;     %long.sound vel. ratio

symmetric=0;

%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=6*lat_a-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)

if lattice==1
   I=8;  			Imax=I;                     J=25;                  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+...
      (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-sqrt(3)/2*lat_a)*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
   %Move the whole array up one notch to set it off from zero
   if yoffset==1
      ry=ry+lat_a/2;
   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)+12*lat_a);
   ytop=ceil(max(ry)+2*lat_a);
   ybot=floor(min(ry)-2*lat_a);

else
%Square Lattice 
  I=1;                 J=1;               Imax=I;
  N=I*J;
  %Values for the columns
  index=(1:N);    col1=ones(J,1);
  rx=xo+lat_a*(ceil(index/(J)));
  rx=rx';
  %Values for the Rows;
  row1=(lat_a)*(-floor(J/2):floor(J/2));
  ry=repmat(row1,I,1);
  ry=ry';

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


%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,kwav.*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)=0;

%For pretty picture purposes---make length x width x 3 true color image
if posneg==1
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))=1;
else
inside_cyls=repmat(inside_cyls,[1 1 3]);
pho_r_scaled=(abs(pho_r)).^2/(max(abs(pho_r(:))).^2);
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))=1;
end

%Determine operating system
osroot=matlabroot;
if osroot(1)=='C'
   filename=['C:\Documents and Settings\Lance\Desktop\School\Qual\',lattice_string,...
             '_Lat_Space_N_',num2str(nmax),...
             '_r_',num2str(a),'_a_'];
   filename=[filename,num2str(lat_a),'_f_',num2str(f),'_xo_',num2str(xo),'_k_'];
   filename=[filename,num2str(kwav),'.tif'];

   imname=['C:\Documents and Settings\Lance\Desktop\School\Qual\im_MST_',...
           lattice_string,'_Lattice_N_',...
           num2str(nmax),'_r_',num2str(a)];
   imname=[imname,'_a_',num2str(lat_a),'_f_',num2str(f),'_xo_',num2str(xo)];
   imname=[imname,'_k_',num2str(kwav),'.tif'];
else
   filename=['figures/MST_',lattice_string','_Lattice_N_',num2str(nmax),'_r_',...
             num2str(a),'_a_'];
   filename=[filename,num2str(lat_a),'_f_',num2str(f),'_xo_',num2str(xo),'_k_'];
   filename=[filename,num2str(kwav),'.tif'];
   imname=['figures/im_MST_',lattice_string,'_Lattice_N_',num2str(nmax),'_r_',...
           num2str(a)];
   imname=[imname,'_a_',num2str(lat_a),'_f_',num2str(f),'_xo_',num2str(xo)];
   imname=[imname,'_k_',num2str(kwav),'.tif'];
end
newfig=figure;


%Generate an intensity distribution with mins and max's
%imagesc([minrx,maxrx],[minry,maxry],abs(pho_r3_scaled))^2;
imagesc([minrx,maxrx],[minry,maxry],pho_r3_scaled);
text(0,1.05,['\lambda=',num2str(lambda),'m'],'units','normalized');
text(0.25,1.03,['r_{o}=',num2str(a),'m'],'units','normalized');
text(0.5,1.03,['a_{lat}=',num2str(lat_a),'m'],'units','normalized');
text(0.75,1.05,['N=',num2str(nmax)],'units','normalized');


imwrite(abs(pho_r3_scaled),imname,'tif');
saveas(newfig,filename,'tif');


pause;
%Attempt at a movie?
%Attempt at a movie?i

if movie==1
   if osroot(1)=='C'
      avistring=['/nfs/slac/g/ki/ki08/lsst/CPanalysis/Phononic_Crys/','MST_',...
              lattice_string,'_Lat_lat_a',num2str(lat_a),'_f_',num2str(f),...
              '_lambda_',num2str(lambda),'.avi'];
      aviobj=avifile(avistring); 
   else
      avistring=['C:\Documents and Settings\Lance\Desktop\School\Qual\','MST_',...
              lattice_string,'_Lat_lat_a',num2str(lat_a),'_f_',num2str(f),...
              '_lambda_',num2str(lambda),'.avi'];
      aviobj=avifile(avistring);
  end
tmax=200;       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:20
           aviobj=addframe(aviobj,frame);
       end
    end
    pause(0.1);

  aviobj=close(aviobj);

end

end