%    penney_price.m
%
%    Construct an animation of the diffraction of an incident wave train
%    by a semi-infinite, planar breakwater, following Penney and Price 195?.
%
%    Angle of incidence 'alf' is arbitrary.

%    Define alf (the angle of incidence) first

     alf = input('input incident wave angle alf in radians:  ');

%    Dimensions of Cartesian grid in x,y.

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

%    Set up dimensions for numerical arrays.

     dx = 12*pi/(nptsx-1);

     x = -6*pi + (0:(nptsx-1))*dx;
     y=x;

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

%    We want angle theta to vary from 0 to 2pi rather 
%    than -pi to pi

     for index1=1:nptsx,
         for index2=1:nptsx,
             if ( theta(index1,index2)<0 ) ...
                theta(index1,index2)=theta(index1,index2)+2*pi;
             end
         end
     end
     
%    Define phases for time stepping - 12 points per period.

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

%    Define sig, sig' in theory.

     sig = 2 * sqrt(rad/pi) .* sin( (theta-alf)/2 );
     sigp = -2 * sqrt(rad/pi) .* sin( (theta+alf)/2 );

%    Define spatial phases

     phase1 = -rad.*cos(theta - alf);
     phase2 = -rad.*cos(theta + alf);

%    Arguments for complex error functions

     arg1 = sqrt( (i*pi)/2 ) * sig;
     arg2 = sqrt( (i*pi)/2 ) * sigp;

%    Complex error functions

     array1=mfun('erf',arg1);
     array2=mfun('erf',arg2);


%    Define surface displacements

     eta1 = exp( i * phase1) .* (1 + array1 )/sqrt(2);

     eta2 = exp( i * phase2) .* (1 + array2 )/sqrt(2);

     eta = eta1 + eta2;


%    Incident wave field:

%    Make a movie looping over the 12 values of phase

     for index=1:12,

     figure(1)

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

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

     M1(index)=getframe;

     end

%    Scattered wave field.

%    Make a movie looping over the 12 values of phase

     for index=1:12,

     figure(2)

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

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

     M2(index)=getframe;

     end 

%    Total wave field.

%    Make a movie looping over the 12 values of phase

     for index=1:12,

     figure(3)

     wave=real(eta*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_pp.avi')
     movie2avi(M2, 'M2_pp.avi')
     movie2avi(M3, 'M3_pp.avi')