%    maccamy-fuchs.m
%
%    MacCamy-Fuchs Theory.
%
%    This program is a conversion of an earlier mathematica
%    program written by Robert A. Dalrymple and heavily 
%    modified by Kirby.
%
%    November 13, 2003.


%    Enter dimensionaless cylinder radius.

     rho0 = input('input cylinder radius kR:  ');

%    number of terms in series for incident and scattered 
%    wave trains.

     Ntermsinc=40;
     Ntermscat=40;

%    Dimensions of Cartesian grid in x,y.

     nptsx = input('input number of points in x,y:  ');

%    Set up dimensions for numerical arrays.

     dx = 8*rho0/(nptsx-1);

     x=-4*rho0 + (0:(nptsx-1))*dx;
     y=x;

     [X,Y]=meshgrid(x,y);
     rad=sqrt( X.*X + Y.*Y);
     theta=atan2(Y,X);

%    Define phases for time stepping - 12 points per period.

     phase = exp( - i*(0:1:11)*(2*pi)/12); 

%    Coefficients of scattered wave field. Note Length of coefs 
%    is Ntermscat+1

     coefs=zeros(Ntermscat+1, 1);

     coefs(1)=besselj(1,rho0)/besselh(1,2,rho0);  %changed to H(2) here.

     for index=1:Ntermscat,

       coefs(index+1) = ( besselj(index-1,rho0) - ...
          besselj(index+1,rho0) ) / ...
         ( besselh(index-1,2,rho0) - besselh(index+1,2,rho0) );  
     end

%    We are going to look at three movies:  The incident wave, 
%    the scattered wave, and the total wave field

%    Define the wave potential functions.

     phii=zeros(size(X)); phis=zeros(size(X));

%    Incident wave field.

     phii = besselj(0,rad);

     for index=1:Ntermsinc, phii = phii + 2* (i)^index * ...
         besselj(index,rad) .* cos(index*theta);
     end

%    scattered wave field.

     phis = - coefs(1)*besselh(0,rad);

     for index=1:Ntermscat,  phis = phis  ...
          -2*(i^index)*coefs(index+1)*besselh(index,rad) ...
          .*cos(index*theta);
     end

%    Incident wave field:

     phi = phii;

%    Mask out the cylinder.

     for indexx=1:nptsx,
       for indexy=1:nptsx,
           if (rad(indexx,indexy) <= rho0 ) 
               phi(indexx,indexy)=NaN; 
           end
       end
     end

%    Now, make a movie looping over the 12 values of phase

     for index=1:12,

     figure(1)

     wave=real(phi*phase(index));

     pcolor(x,y,wave),shading flat, axis('equal')

     M1(index)=getframe;

     end

%    Scattered wave field.

     phi = phis;

%    Mask out the cylinder.

     for indexx=1:nptsx,
       for indexy=1:nptsx,
           if (rad(indexx,indexy) <= rho0 ) 
              phi(indexx,indexy)=NaN; 
           end
       end
     end

%    Now, make a movie looping over the 12 values of phase

     for index=1:12,

     figure(2)

     wave=real(phi*phase(index));

     pcolor(x,y,wave),shading flat, axis('equal')

     M2(index)=getframe;

     end 

%    Total wave field.

     phi = phii  + phis;

%    Mask out the cylinder.

     for indexx=1:nptsx,
       for indexy=1:nptsx,
           if (rad(indexx,indexy) <= rho0 ) 
              phi(indexx,indexy)=NaN; 
           end
       end
     end

%    Now, make a movie looping over the 12 values of phase

     for index=1:12,

     figure(3)

     wave=real(phi*phase(index));

     pcolor(x,y,wave),shading flat, axis('equal')

     M3(index)=getframe;

     end

%    Play back the movie.

     figure(1)

     movie(M1,10)

     figure(2)

     movie(M2,10)

     figure(3)

     movie(M3,10)

%   Save the movie.

     movie2avi(M1, 'M1.avi')
     movie2avi(M2, 'M2.avi')
     movie2avi(M3, 'M3.avi')