SIO227A Homework 5: 6.2 (Eric Gallimore)
Ch 6 #2
% Values from IASP model, in [deg, s, s/deg] iasp = [29,361.41,8.89; 30,370.27,8.85; 31,379.10,8.81;59,601.38,6.95;60,608.29,6.88;61,615.13,6.8;89,776.67,4.69;90,781.35,4.66;91,786,4.64]; % surface velocity v_surf = 5.8; % km/s u1 = 1/v_surf; % we are interested in rays traveling from the surface to the surface, so r1 = 6378; % km, at source r2 = 6378; % km, at receiver % convert to radians, s/radian, make variables easier to read T = iasp(:,2); % time, seconds Delta = pi/180 .* iasp(:,1); % distance, radians (of sphere) p = 180/pi .* iasp(:,3); % slowness, s/rad % get dp/dDelta diffpD = diff(p) ./ diff(Delta); % Our values are actually off by half a delta step. Apply a moving average % to get them to align. Drop the first element because it's an artifact of % the startup of this method. dPdD = filter([0.5 0.5], 1, [diffpD(1); diffpD; diffpD(end)]); dPdD = dPdD(2:end); % find takeoff angle using p = u0*sin(theta) % note that the surface velocity is the same for both the source and % receiver, so we get the same angle. theta = atan(p./u1); % Plot it to make sure it looks OK % plot(Delta, theta); % Use 6.23 E_coeff = (p ./ ((4*pi*u1^2*r1^2*r2^2) .* (cos(theta).^2) .* sin(Delta))) .* abs(dPdD); % Take square root to get amplitudes amp_coeff = sqrt(E_coeff); % Look at the results for 30, 60, 90 deg. idx = ismember(iasp(:,1), [30, 60, 90]); amp = amp_coeff(idx); % Compare the values at 60 and 90 deg to 30 deg (make 30 unity) amp_scaled = (1/amp(1)) .* amp; % and print the results fprintf('Distance (deg)\tRelative Amplitude\n'); fprintf('30\t%.4f\n', amp_scaled(1)); fprintf('60\t%.4f\n', amp_scaled(2)); fprintf('90\t%.4f\n', amp_scaled(3));
Distance (deg) Relative Amplitude 30 1.0000 60 0.7132 90 0.2136