% raytrace_lab.m
%
% raytrace_lab is a modification of Mark Orzech's "raytrace.m" program for 
% use in the OC4213/OC4610 labs. 
%
% Differences between this program and raytrace.m:
%     - this version is interactive and prompts the user for offshore wave
%       direction and frequency.
%     - the user clicks the mouse on the offshore location from which to
%       start the refraction and can pick as many locations as desired.
%     - The NCEX stations sites 1-12,18, and 19 may be overlayed as
%       "targets" for the rays (Paul Jessen, December 2006)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Original documentation:
%
% Trace wave rays over NCEX bathymetry using "Snell's Law"-type refraction.
% Input is specified by modifying "initial values" in program below.
% Output includes plot of wave rays over bathymetry and
%     xsall = [N x NR] matrix of x-coordinates of all rays plotted (meters),
%             where N is total number of steps (of size "ds") taken
%             and NR is total number of rays (equal to length(ys) below)
%     ysall = [N x NR] matrix of corresponding y-coordinates (meters)
%     thetall = [N x NR] matrix of corresponding ray direction values (degrees).
%    
% Files required to run program:
%
%  sd*msfxyns_v3b.mat = bathymetry data file, containing xi,yi vectors and 
%                       zi matrix of depths (positive downward from MWL)
%  sdcoast.mat = data file for plotting coastline contour (xcoast, ycoast)
%   
% Mark Orzech 8/23/04
if exist('zi','var') == 0,
  clear
  pack
% load xy-bathymetry, if necessary
  bathfile='sd_decimate_fac2.mat';
%bathfile='sd037.5msfxyns_v3b.mat';             % Bathymetry filename
  disp('Loading bathymetry data ...');
  load(bathfile);
  minz=5;
  disp('Done.');
  %clear dcdx     % So that new wave speed values will be calculated below
  lenxi=length(xi); lenyi=length(yi);
end
%
% Prompt user for offshore wave direction and frequency
%
prcnt=0;
diropt = 0;
while diropt < 4
  diropt=menu('Choose Offshore Wave Direction','240','270','300','Quit');
  if diropt == 4, break, end
%osdir=input('Enter Offshore Wave Direction(200-340): ');
  per=menu('Choose Offshore Wave Period','5 seconds','10 seconds','20 seconds');
  freq=1/20;
  if per == 1, freq=1/5;, end
  if per == 2, freq=1/10;, end
  osdir=240;
  if diropt == 2, osdir=270;, end;
  if diropt == 3, osdir=300;, end;
  prcnt=prcnt+1;
%freq=input('Enter Offshore Wave Frequency (hz): ');
  %
  % convert direction and frequency to program units
  %
  theta=270-osdir; % Initial wave direction, degrees (CCW from 270)
  sig=2*pi*freq;  % Wave frequency, radians 
  %
  % Calculate wave speed and its gradients at all points on bathymetry, if req'd
  %
  if exist('dcdy','var') == 0,
    disp('Generating wave speed values ...');
    c=sig./waveno(sig,zi);
    dcdx=diff(c,1,2)/(xi(2)-xi(1));
    dcdy=diff(c,1,1)/(yi(2)-yi(1));
    dcdx=[dcdx dcdx(:,end)];
    dcdy=[dcdy; dcdy(end,:)];
    disp('Done.');
  end
  %
  % create a contour plot of the region
  %
  % Plot bathy and coastline
  %
  figure; 
  contour(xi,yi,zi,[50:50:max(max(zi))]);
  axis('square')
  hold on;
  if prcnt == 1, load sdcoast.mat;, end
  plot(xcoast,ycoast,'k*','markersize',2);
  title('NCEX Study Area')
  % Set some initial values, get axes limits
  ds=1.25;  %2.5;  %5;           % Step size along wave rays, in meters
  xlm=get(gca,'xlim');
  ylm=get(gca,'ylim');
  %
  % overlay station positions
  %
%   if exist('sta_x','var') == 0, load stapos, end
%   ph=plot(sta_x,sta_y,'ys');
%   text(sta_x+5,sta_y,num2str(id));
%
% overlay region boundaries
%
if prcnt == 1,
  if exist('region_bounds.mat','file') == 2,
    load('region_bounds.mat');
  end
end
if exist('x','var') == 1,
  for j=1:2:12
    line(x(j:j+1),y(j:j+1),'linewidth',2,'color','k')
  end
  for j=1:length(xl)
    text(xl(j),yl(j),num2str(j))
  end
end
%
% put wave directions and frequencies on map
%
tstr1=['Wave Direction: ',num2str(osdir),' deg'];
tstr2=['Wave period: ',sprintf('%4.1f',1/freq),' sec'];
text(-1500,-2200,tstr1);
text(-1500,-2500,tstr2);
  %
  % now have main loop where user will click on position to start refraction
  %
  while 1
    disp('Click on offshore starting position (click outside map to quit)...')
    [xs,ys]=ginput(1); % Initial x/y coordinates of rays
    if xs < xlm(1) || xs > xlm(2), break, end
    if ys < ylm(1) || ys > ylm(2), break, end
    %
    % check to see if this depth is too shallow
    %
    [vl,aix]=min(abs(xi-xs));
    [vl,aiy]=min(abs(yi-ys));
    osd=zi(aiy,aix);
    if osd < minz
      errordlg('Offshore depth too shallow, try again');
    else
      % Prepare for calculations
      thet=theta*pi/180; % convert initial offshore direction to radians
      xsall=xs; ysall=ys; thetall=thet; 
      % Progress wave rays through the bathymetry
      disp('Calculating wave rays ...');
      while 1
        [minx,ix]=min(abs(xi*ones(size(xs))-ones(size(xi))*xs));
        [miny,iy]=min(abs(yi'*ones(size(ys))-ones(size(yi'))*ys));
        idat=lenyi*(ix-1)+iy;
        % Calculate new locn and angle values
        xs=xs + ds*cos(thet);
        ys=ys + ds*sin(thet);
        thet=thet + ds*(sin(thet).*dcdx(idat)-cos(thet).*dcdy(idat))./c(idat);
        % If any results are imaginary set those values to NaN
        idx=find(imag(xs)~=0 | imag(ys)~=0 | imag(thet)~=0);
        xs(idx)=NaN; ys(idx)=NaN; thet(idx)=NaN;
        % If any rays outside data region, set those values to NaN
        idx=find(xs<min(xi) | xs>max(xi) | ys<min(yi) | ys>max(yi));
        xs(idx)=NaN; ys(idx)=NaN; thet(idx)=NaN;
        % If any rays near zero depth or out of water (<=0), set those values to NaN
        idx=find(zi(idat)<=0);
        xs(idx)=NaN; ys(idx)=NaN; thet(idx)=NaN;
        % If all NaN in one of the variables, quit.
        if all(isnan(xs)) | all(isnan(ys)) | all(isnan(thet))
          break;
        end
        % Append new results to output matrix
        xsall=[xsall; xs];
        ysall=[ysall; ys];
        thetall=[thetall; thet];
      end % while
      % Convert theta values back to degrees
      thetall=thetall*180/pi;
      % Add path results to plot
      plot(xsall,ysall,'k*','markersize',1);
      % If desired, set axes to correct relative scaling
  %    axis equal; axis([min(xi) max(xi) min(yi) max(yi)]);
      %
      %  now see if the ray came close to any of the stations
      %
  %     if exist('sta_x','var') == 1,
  %       for j=1:length(sta_x)
  %         xdst=xsall-sta_x(j);
  %         ydst=ysall-sta_y(j);
  %         dst=sqrt(xdst.^2+ydst.^2);
  %         [vldst(j),dind(j)]=min(dst); 
  %       end
  %       %
  %       % find points within 15 meters of one of the stations
  %       %
  %       cp=find(vldst<=15);
  %       if isempty(cp) == 0,
  %         disp('Ray Intersected one or more instrument sites!');
  %         fprintf(1,'\nSite   Distance(m)    Wave Angle');
  %         for j=1:length(cp)
  %           fprintf('\n %d %f %f',id(cp(j)), vldst(cp(j)), thetall(cp(j)));
  %         end
  %         fprintf(1,'\n');
  %       end
  %     end
  %
    end
  end
  %
  % save the figure if desired
  %
  qstr=questdlg('Save figure?');
  if strcmp(qstr,'Yes') == 1,
    tp=round(1/freq);
    onme=[num2str(osdir),'deg_',num2str(tp),'sec'];
    eval(['print -djpeg99 ',onme]);
    fprintf('\n Figure saved as: %s.jpg\n',onme);
  end
  %
  % clear the figure for the next run
  %
  clear dcdx dcdy
  close(gcf)
end
disp('Done.');