% Compute the range independent modal (rim) acoustic field using the 
% Legendre-Galerkin method: Edit the generic input file rimLG.dat;
% Read the imput file, run rimLG and plot complex horizontal wave
% numbers, modes , acoustic field and transmission loss versus range.
% The transmisson loss is written on the text file rimLG.tl for post 
% processing.
%
% Calls read_rimLG, rimLG, pkh_plot, mode_plot, tl_plot and fld_plot.
%
% See rimLG_readme.txt for user instructions. See Legendre_Galerkin.pdf 
% for documentation on the computational method.
% (R. B. Evans, N. Stonington, CT, Nov. 17, 2016)
%
clear
reply_edit='Y';
display('Starts by opening rimLG.dat in the editor');
display('Press any key to continue');
pause
% Edit the data file
while ((reply_edit == 'Y')||(reply_edit == 'y'))
edit 'rimLG.dat' % Brings editor back if changes have not been saved.
% Make changes and save rimLG.dat for the desired caclulation.
display(' ');
display('Make changes to rimLG.dat, if any.');
display('Save rimLG.dat for the desired calculation.');
display(' ');
display(['*** Click in the command window and ' ...
                          'press any key to continue ***']);
pause
% Confirm that any required changes have been saved
reply_save = input('Have any required changes been saved? y/n :','s');
   if ((reply_save == 'Y')||(reply_save == 'y'))
      reply_edit='N';
   end
end
% The editor may be closed at this point.
display(' ');
display('You may close the editor now.');
% Give option to quit
display(' ');
reply_run = input('Do you want to run rimLG? y/n :','s');
if ((reply_run == 'Y')||(reply_run == 'y')) % Proceed with calculation
%
% Define text data file
data_file='rimLG.dat';
% Read the data file
[case_name,FREQ,ZS,M,MS,NL,NDEP,ZMIN,ZINC,RMIN,RMAX,RINC,...
 HW,HB,depw,svpw,dbpwlw,rhow,NW,depb,svpb,dbpwlb,rhob,NB]= ...
                                         read_rimLG(data_file);
% Compute complex pressure on the depth-range grid
display('rimLG calculation started: please wait.');
[zr,cdfr,rkm,kh,cpr,NRNG]= ...
    rimLG(FREQ,ZS,NL,MS,NDEP,ZMIN,ZINC,RMIN,RMAX,RINC, ...
    HW,HB,depw,svpw,dbpwlw,rhow,NW,depb,svpb,dbpwlb,rhob,NB);
display('rimLG calculation completed.');
% Plot the wave numbers if requested
reply_kh = input('Do you want to plot the wave numbers? y/n :','s');
if((reply_kh == 'Y')||(reply_kh == 'y'))
% Plot the horizontal wave numbers in the complex plane.
% pkh are the wave numbers used in the sum to form the complex pressure
   pkh=kh(1:MS);
   pkh_plot(pkh,FREQ,case_name)
end;
% Plot a mode functions if requested
reply_mode = input('Do you want to plot a mode? y/n :','s');
%
while((reply_mode == 'Y')||(reply_mode == 'y'))
  mode_num = input('What mode number do you want to plot?:');
  display('Mode plot started.');
  % Plot mode
  mode_plot(mode_num,kh,cdfr,FREQ,zr,case_name)
  reply_mode = input('Do you want to plot another mode? y/n :','s');
end
% Define text transmission loss file
tl_file='rimLG.tl';
tl_fid = fopen(tl_file,'wt');
fprintf(tl_fid,'%s\n',case_name);
today11=date;
today=[today11(4:6) today11(1:2) '_' today11(8:11)];
NRNG0=0;
head=[FREQ ZS NRNG0 NDEP];
% Write header (spaces make spaces) and receiver depths
fprintf(tl_fid,'%10s %13.7f %13.7f %5d %5d\n',today,head); 
fprintf(tl_fid,'%13.7f',zr);
fprintf(tl_fid,'\n');
% Compute transmission loss
display('Transmission loss calculation started.');
tl=zeros(NDEP,NRNG);
   for n=1:NRNG
      for j=1:NDEP
      tl(j,n)=-20.0*log10(abs(cpr(j,n)*(4*pi)));
      end
   rkm_tl=[rkm(n); tl(1:NDEP,n)];
   fprintf(tl_fid,'%13.7f',rkm_tl);
      if n~=NRNG
      fprintf(tl_fid,'\n'); 
      end
   end
fclose(tl_fid);
display('Transmission loss calculation completed.');
%
reply_fld = input('Do you want to plot the acoustic field? y/n :','s');
if((reply_fld == 'Y')||(reply_fld == 'y'))
   % Plot tranmission loss field
   dbmin = input('What min loss (dB) do you want to plot?:');
   dbmax = input('What max loss (dB) do you want to plot?:');
   fld_plot(tl,rkm,zr,dbmin,dbmax,case_name)
 end
% Plot transmisson loss if requested
reply_tl = input('Do you want to plot loss vrs. range? y/n :','s');
%
while((reply_tl == 'Y')||(reply_tl == 'y'))
  rec_dep = input('What receiver depth (m) do you want to plot?:');
  % Plot transmission loss at zr nearest to rec_dep
  dep_dif=abs(rec_dep-zr);
  min_dif=min(dep_dif);
  j=find(dep_dif == min_dif);
  NR=j(1);
  display('Transmission loss plot started.');
  tl_plot(NR,NDEP,NRNG,tl,rkm,ZS,FREQ,zr,case_name)
  reply_tl = input('Do you want to plot another receiver? y/n :','s');
end
else
display('Calculation avioded.');
end % run rimLG