%This script runs an ode solver on the eikonal equation like described in %the problem set. The velocity model is set by functions profile.m and c.m. %Initial angle is given by theta in degrees (0 to 90). %The solution is integrated until the ray hits the surface, the depth of 50, %or distance 100 or the time 40, whichever comes first. This can be changed by editing the files. %This program is not foolproof, mainly because the ode solver can run into %problems with integration. The ray can also get caught in a waveguide, but escape %again because of numerical instabillities. Comments are scarce and random, %most should be self explanatory if some thought is given to it. Plotrays.m %plots some of the solution %****** F R A G I L E ---- H A N D L E W I T H C A R E ****** %Written by Keli Karason, karason@mit.edu, 8 Nov. 1999 [zz,cc,theta]=getinput; %Initial and final times. to=0; tfinal=40; %Initial position in space xo=0; zo=0; theta=theta(:); %Initial slowness vector from initial angles pxo=sin(theta/180*pi)/veldep(zo); pzo=cos(theta/180*pi)/veldep(zo); %Getting lengths nth=length(theta); %number of rays to trace nr=100; %initializing variables t=zeros(100,nth); x=zeros(100,nth); z=zeros(100,nth); px=zeros(100,nth); pz=zeros(100,nth); %Running the numerical integration and interpolating to get even number of time steps options=odeset('Events','on'); for j=1:nth [tt,yy]=ode45('yprime',[to tfinal],[xo zo pxo(j) pzo(j)]',options); t(:,j)=linspace(0,tt(end),nr)'; x(:,j)= interp1(tt,yy(:,1),t(:,j),'linear'); z(:,j)= interp1(tt,yy(:,2),t(:,j),'linear'); px(:,j)=interp1(tt,yy(:,3),t(:,j),'linear'); pz(:,j)=interp1(tt,yy(:,4),t(:,j),'linear'); disp([num2str(j) ' out of ' num2str(nth) ' finished']); end %Getting values for wavefronts twf=(1:1:20); xwf=zeros(nth,length(twf)); zwf=zeros(nth,length(twf)); for j=1:nth, xwf(j,:)=interp1(t(:,j),x(:,j),twf); zwf(j,:)=interp1(t(:,j),z(:,j),twf); end plotrays