%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Program Burgers
%
% Translated into MATLAB by Ronald Kumon
% 
% First Translated:  21 Oct 1995
% Last Revised    :  21 Oct 1995
% Uses routines   :  F.m, PconvP.m, abstv.m, absrel.m, absbl.m
%
% This program solves the Burgers equation using a frequency domain
% technique.  Arbitrary absorption and dispersion may be included
% in the algorithm.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Declaration of global variables
global a n nhar tv d omega b

% Input parameters
clc
fprintf( '***********************************************\n');
fprintf( '*    Solution of Spectral Burgers Equation:   *\n');
fprintf( '*             Propagation Curves              *\n');
fprintf( '***********************************************\n');
step=input('Enter sigma step size (.01): ');
nout=input('Enter number of steps per output: ');
nhar=input('Enter number of harmonics in calculation (50): ');
xmax=input('Enter the max sigma range for output: ');
tv  =input('Enter the thermoviscous loss coefficient: ');
clc

% Define constants and initialize variables
x=0;				% Initial position
omega=1;			% Fundamental frequency
n=1:nhar;			% Harmonic amplitude index vector
p=zeros(1,nhar);		% Initialize harmonic amp. vector
R=zeros(1,nhar);		% Runge-Kutta intermediate vectors
S=zeros(1,nhar);		%	
a=abstv(n);			% Absorption coefficient vector

% Initial condition at source
p(1)=-j;			% Pure sine wave source condition 

% Second order Runge-Kutta integration routine
Nsteps=fix(xmax/step);		% Number of steps in spatial integration
x=0:step:xmax;			% Position vector
ndisplay=3; 			% Display ndisplay harmonics
pmatrix=p; 			% Cumulative harmonic amplitude matrix
pfubinimatrix=zeros(1,ndisplay); % Cumulative Fubini coefficient matrix
pfubinimatrix(1)=1;
pdiffmatrix  =zeros(1,ndisplay); % Cumumlative percent difference matrix
  for k=1:Nsteps; 		% Begin integration routine
  R=step*F(p);			% Runge-Kutta intermediate vectors
  S=step*F(p+R); 		%
  p=p+(R+S)/2;			% Update harmonic amplitude vector
  
  if ( rem(k,nout)==0 )		% Print results every nout steps
    clc 			% Clear screen 
    fprintf('x = %6.4f\n\n',x(k+1))
    for index=1:ndisplay 			
      y             =index*x(k+1); 	% Compute Fubini coefficients
      pfubini(index)=2/y*bessel(index,y);
      pdiff(index)  =(abs(p(index))-pfubini(index))/pfubini(index)*100;
      indexstr=num2str(index);
      fprintf(['|P(',indexstr,')|= %10.8f   B(',indexstr,')= %10.8f', ...
	       '   Percent diff= %10.8f\n'],abs(p(index)),pfubini(index), ...
               pdiff(index)) 
     end
     pmatrix=[pmatrix;p]; 		% Update cumulative matrices
     pfubinimatrix =[pfubinimatrix;pfubini];
     pdiffmatrix   =[pdiffmatrix;pdiff];
     ptotal =sum(abs(p).^2); 		% Sum of square mod of ampl. comp.
     pfubinitotal=sum(pfubini.^2);      % Sum of squares of Fubini coeff.
     fprintf('\nSum of |P(n)|^2= %10.8f   Sum of B(n)^2= %10.8f\n\n', ...
	 ptotal,pfubinitotal)
     pause
     fprintf('Working...\n')
  end
end 					% End integration routine

% Save key results
last=fix(Nsteps/nout)+1;
amp=abs(pmatrix(last,1:ndisplay));

% Plot graph
figure(1)
clf
hold on
xincr=0:(nout*step):xmax;
plot(xincr,abs(pmatrix(:,1)),'y') 	% Plot harmonic amplitudes (numerical)
plot(xincr,abs(pmatrix(:,2)),'r')
plot(xincr,abs(pmatrix(:,3)),'g')
plot(xincr,pfubinimatrix(:,1),'y--') 	% Plot harmonic amplitudes (Fubini)
plot(xincr,pfubinimatrix(:,2),'r--')
plot(xincr,pfubinimatrix(:,3),'g--')
axis([0 1 0 1])
title(['Problem #2     Parameters:   N = ',num2str(nhar), ...
       '  A_1 = ',num2str(tv)])
xlabel('sigma')
ylabel('|P_n|')
grid on
hold off

figure(2)
clf
hold on
plot(xincr,pdiffmatrix(:,1),'y')
plot(xincr,pdiffmatrix(:,2),'r--')
plot(xincr,pdiffmatrix(:,3),'g:')
title(['Problem #1     Parameters:   N = ',num2str(nhar), ...
       '  A_1 = ',num2str(tv)])
xlabel('sigma')
ylabel('Percent difference')
grid on
hold off

% Print graph in encapsulated Postscript format, if desired
%prt=input('Print graph? (1=yes,0=no): ');
%if prt==1
%  filename=input('Input filename in single quotes: ');
%end
%if prt==1 
%  eval(['print ',filename,' -deps']);
%end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of Program Burgers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%