%refrac.m , an m-file for oc4213 to plot refractions rays and calculate
%Ks and Kr
%!rm sxsy.dat;
%!rm box.dat;
%!rm input.dat;
%!rm boxes.dat;
%!rm RKSK.dat;
%!rm XY;
% 
% Create refract.output file if it doesn't exist
if ~exist('refract.output','file')
   fid=fopen('refract.output','w'); fclose(fid);
end

skr=[];sks=[];sxf=[];syf=[];saf=[];
!rinit.exe
clf
hold off
load pshore.dat
load mbdepth
load box.dat
contour(flipud(depth),[-2400:200:-200 -100 -50 -20])
hold on;
plot(shorex,shorey); 
title('Monterey Bay')
xbox=box(:,1);ybox=box(:,2);
dbx=xbox(2)-xbox(1); xbox([1 4 5])=xbox([1 4 5])+dbx/2.5;
xbox(2:3)=xbox(2:3)-dbx/2.5;
dby=ybox(3)-ybox(2); ybox([1 2 5])=ybox([1 2 5])+dby/2.5;
ybox(3:4)=ybox(3:4)-dby/2.5;
plot(xbox,ybox); hold on; fill(xbox,ybox,'y');
disp(['    '])
disp(['When you want to stop plotting refraction rays type 0.'])
disp(['Otherwise just hit RETURN once before each new ray plotted.'])
disp(' ');
for i = 1:100
 X = 0;        
 Y = 0;
 stop = input('Hit RETURN to plot, enter 1 to clear plot, or enter 0 to quit:  ');
 if ~isempty(stop)
  if stop == 0
   break
  end
 end
 if ~isempty(stop)
  if stop == 1
  
   % Replot Monterey Bay
   clf
   hold off
   load pshore.dat
   load mbdepth
   load box.dat
   contour(flipud(depth))
   hold on;
   plot(shorex,shorey); 
   title('Monterey Bay')
   xbox=box(:,1);ybox=box(:,2);
   dbx=xbox(2)-xbox(1); xbox([1 4 5])=xbox([1 4 5])+dbx/2.5;
   xbox(2:3)=xbox(2:3)-dbx/2.5;
   dby=ybox(3)-ybox(2); ybox([1 2 5])=ybox([1 2 5])+dby/2.5;
   ybox(3:4)=ybox(3:4)-dby/2.5;
   plot(xbox,ybox); hold on; fill(xbox,ybox,'y');
  
  end
  
 else
  
  % Add another ray
   [X,Y]=ginput(1);
   save XY X Y /ascii;
   !refract.exe
   load sxsy.dat
   load ideep.dat
   if ideep==-9999       %if not in deep water try again
     disp('Must start in deep water - try again')
   else
     sx=sxsy(:,1);sy=sxsy(:,2);
     plot(sx,sy)
     load boxes.dat
     if ~isempty(boxes)
      x13=boxes(1);xp1=boxes(2);xp2=boxes(3);y13=boxes(4);yp1=boxes(5);yp2=boxes(6);
%      [x13 xp1 xp2 y13 yp1 yp2]
      if x13>=xp1 & x13<=xp2 & y13>=yp1 & y13<=yp2,
       load RKSK.dat
       xf=RKSK(1);yf=RKSK(2);af=RKSK(3);kr=RKSK(4);ks=RKSK(5);
       sxf=[sxf xf];syf=[syf yf];saf=[saf af];
       skr=[skr kr];
       sks=[sks ks];
       numberhits=length(skr);
       disp(['HIT TARGET! Total number of hits: ' int2str(numberhits)]);
       delete boxes.dat
      end
     end % if ~isempty(boxes)
   end % if ideep
   
 end % if ~isempty(stop) 
 
end % for i

numberhits=length(skr);
if (numberhits>0)
  disp('Mean values of coefficients obtained:  ');
  Kr=mean(skr);
  Ks=mean(sks); 
  save KrKs numberhits  sxf syf saf skr sks Kr Ks /ascii
  !rfin.exe
  disp('Results saved to file refract.output.');
end
