From 615e0a35162342d91315b11bd79ca879c292f1d1 Mon Sep 17 00:00:00 2001 From: Alan Baird Date: Tue, 15 Jul 2014 14:42:36 +0100 Subject: [PATCH 01/12] added full ebsd data from Obbard and Baker (2007) --- .../glacial_splitting/glacial_splitting.m | 108 ++++++++++++++---- 1 file changed, 84 insertions(+), 24 deletions(-) diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index 232db01..88572fc 100644 --- a/examples/glacial_splitting/glacial_splitting.m +++ b/examples/glacial_splitting/glacial_splitting.m @@ -147,6 +147,19 @@ report_layer(i, all_data(i).tex_file, Cs(:,:,j), rhos(j), ... all_data(i).dtb, dtt, all_data(i).Tav, plot_pole) end + + % Plot stratigraphic column + + figure(); + bar([fliplr(thickness);nan(size(thickness))],'stack');xlim([0.5,1.5]); + set(gca,'YDir','Reverse'); + hold on; + samples=[97 248 399 650 1201 1500 1849 2082 2749 2874 3300 3311 3321 3329 3399 3416]./1000.0; + %scatter(ones(size(samples)),samples,'MarkerEdgeColor','k','MarkerFaceColor','w') + scatter(linspace(0.7,1.3,max(size(samples))),samples,'MarkerEdgeColor','k','MarkerFaceColor','w') + hold off; + + % The effective splitting calculation fprintf('\nEffective splitting calculation \n'); @@ -301,40 +314,87 @@ % Deal with file seperators in a cross platform way. dat_dir = fullfile('~', 'Dropbox', 'Ice_EBSD_data'); - % FIXME: AFB - A comment here as to where these layers - % come from is critical - data(1) = struct('tex_file', fullfile(dat_dir, 'V97-248R.txt'), ... - 'dtb', 450, ... % depth to base of layer m + % EBSD data from Obbard and Baker (2007). Depths to deeper layers from Lipenkov and Barkov (1998) + data(1) = struct('tex_file', fullfile(dat_dir, 'V97.txt'), ... % R-A + 'dtb', 200, ... % depth to base of layer m 'Tav', -25, ... % Average temperature (deg. C) 'azi', 0); % Horizonal rotation for texture - data(2) = struct('tex_file', fullfile(dat_dir, 'V3311g.txt'), ... - 'dtb', 2850, ... - 'Tav', -20, ... + data(2) = struct('tex_file', fullfile(dat_dir, 'V248.txt'), ... % R-B + 'dtb', 300, ... + 'Tav', -23, ... + 'azi', 0); + data(3) = struct('tex_file', fullfile(dat_dir, 'V399-3.txt'), ... % R-C + 'dtb', 454, ... + 'Tav', -22, ... + 'azi', 0); + data(4) = struct('tex_file', fullfile(dat_dir, 'V650.txt'), ... % A0-A + 'dtb', 900, ... + 'Tav', -19, ... + 'azi', 0); + data(5) = struct('tex_file', fullfile(dat_dir, 'V1201.txt'), ... % A0-B + 'dtb', 1300, ... + 'Tav', -16, ... + 'azi', 0); + data(6) = struct('tex_file', fullfile(dat_dir, 'V1500.txt'), ... % A0-C + 'dtb', 1700, ... + 'Tav', -14, ... + 'azi', 0); + data(7) = struct('tex_file', fullfile(dat_dir, 'V1849.txt'), ... % A0-D + 'dtb', 1900, ... + 'Tav', -12, ... + 'azi', 0); + data(8) = struct('tex_file', fullfile(dat_dir, 'V2082-M.txt'), ... % A0-E + 'dtb', 2700, ... + 'Tav', -7, ... 'azi', 0); - data(3) = struct('tex_file', fullfile(dat_dir, 'V3321c.txt'), ... - 'dtb', 2900, ... - 'Tav', -15, ... + data(9) = struct('tex_file', fullfile(dat_dir, 'V2749-132.txt'), ... % A1 + 'dtb', 2835, ... + 'Tav', -6, ... 'azi', 0); - data(4) = struct('tex_file', fullfile(dat_dir, 'V3311g.txt'), ... - 'dtb', 3150, ... - 'Tav', -10, ... + data(10) = struct('tex_file', fullfile(dat_dir, 'V2874-P9.txt'), ... % B1 + 'dtb', 2905, ... + 'Tav', -5, ... 'azi', 0); - data(5) = struct('tex_file', fullfile(dat_dir, 'V3321c.txt'), ... - 'dtb', 3225, ... - 'Tav', -9, ... + data(11) = struct('tex_file', fullfile(dat_dir, 'V2082-M.txt'), ... % A2 ------ girdle + 'dtb', 3140, ... + 'Tav', -4, ... 'azi', 0); - data(6) = struct('tex_file', fullfile(dat_dir, 'V3311g.txt'), ... - 'dtb', 3300, ... - 'Tav', -7, ... + data(12) = struct('tex_file', fullfile(dat_dir, 'V3321c.txt'), ... % B2 ------- cluster + 'dtb', 3220, ... + 'Tav', -3.5, ... 'azi', 0); - data(7) = struct('tex_file', fullfile(dat_dir, 'V3321c.txt'), ... + data(13) = struct('tex_file', fullfile(dat_dir, 'V3300_M.txt'), ... % A3 + 'dtb', 3310, ... + 'Tav', -3, ... + 'azi', 0); + data(14) = struct('tex_file', fullfile(dat_dir, 'V3311g.txt'), ... % CA + 'dtb', 3330, ... + 'Tav', -2.8, ... + 'azi', 0); + data(15) = struct('tex_file', fullfile(dat_dir, 'V3321c.txt'), ... % CB 'dtb', 3350, ... - 'Tav', -5, ... + 'Tav', -2.6, ... + 'azi', 0); + data(16) = struct('tex_file', fullfile(dat_dir, 'V3311g.txt'), ... % CA + 'dtb', 3370, ... + 'Tav', -2.5, ... 'azi', 0); - data(8) = struct('tex_file', fullfile(dat_dir, 'V3311g.txt'), ... - 'dtb', 3450, ... + data(17) = struct('tex_file', fullfile(dat_dir, 'V3311g.txt'), ... % A4 + 'dtb', 3460, ... + 'Tav', -1.9, ... + 'azi', 0); + data(18) = struct('tex_file', fullfile(dat_dir, 'V3311g.txt'), ... % DA+B + 'dtb', 3540, ... + 'Tav', -1.4, ... + 'azi', 0); + data(19) = struct('tex_file', fullfile(dat_dir, 'V3321c.txt'), ... % E + 'dtb', 3605, ... 'Tav', -1, ... - 'azi', 0); + 'azi', 0); + data(20) = struct('tex_file', fullfile(dat_dir, 'V3321c.txt'), ... % ? + 'dtb', 3750, ... + 'Tav', -1, ... + 'azi', 0); end From 48429716c3589fdc32767cc0c4de5143919cc114 Mon Sep 17 00:00:00 2001 From: Alan Baird Date: Wed, 16 Jul 2014 12:12:12 +0100 Subject: [PATCH 02/12] added graphs of temperature and velocities vs depth --- .../glacial_splitting/glacial_splitting.m | 50 ++++++++++++++++--- 1 file changed, 43 insertions(+), 7 deletions(-) diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index 88572fc..680900a 100644 --- a/examples/glacial_splitting/glacial_splitting.m +++ b/examples/glacial_splitting/glacial_splitting.m @@ -122,6 +122,10 @@ thickness = zeros(1,length(all_data)); Cs = zeros(6,6,length(all_data)); rhos = zeros(1,length(all_data)); + ddts = zeros(1,length(all_data)); + ddbs = zeros(1,length(all_data)); + Vpiso = zeros(1,length(all_data)); + Vsiso = zeros(1,length(all_data)); j = 0; fprintf('\nLayer data (bottom to top) \n'); fprintf('---------------------------\n'); @@ -137,28 +141,60 @@ dtt = all_data(i-1).dtb; end + dtts(j) = dtt; + dtbs(j) = all_data(i).dtb; + temps(j) = all_data(i).Tav; + + + + % Single crystal elasticity of this layer [C, rho] = ice_cij(all_data(i).Tav); % Calculate poly xtal elasticity of this layer [Cs(:,:,j), rhos(j)] = Cs_from_EBSD_file(C, rho, ... all_data(i).tex_file); + + % Calculate average isotropic velocities + CR = MS_axes(Cs(:,:,j), 'nowarn') ; + [Ciso] = MS_decomp(CR) ; + Vpiso(j) = 1e-3.*sqrt((Ciso(3,3)*1e9)./rhos(j)); + Vsiso(j) = 1e-3.*sqrt((Ciso(4,4)*1e9)./rhos(j)); report_layer(i, all_data(i).tex_file, Cs(:,:,j), rhos(j), ... all_data(i).dtb, dtt, all_data(i).Tav, plot_pole) end + Vpiso + Vsiso + % Plot stratigraphic column figure(); - bar([fliplr(thickness);nan(size(thickness))],'stack');xlim([0.5,1.5]); + subplot(1,3,1); + bar([fliplr(thickness);nan(size(thickness))],'stack');xlim([0.5,1.5]);ylim([0,4]); set(gca,'YDir','Reverse'); hold on; samples=[97 248 399 650 1201 1500 1849 2082 2749 2874 3300 3311 3321 3329 3399 3416]./1000.0; %scatter(ones(size(samples)),samples,'MarkerEdgeColor','k','MarkerFaceColor','w') scatter(linspace(0.7,1.3,max(size(samples))),samples,'MarkerEdgeColor','k','MarkerFaceColor','w') + ylabel('depth') hold off; + subplot(1,3,2); + stairs(temps,dtts./1000);xlim([-25,0]);ylim([0,4]); + ylabel('depth'); + xlabel('Temperature'); + set(gca,'YDir','Reverse'); + + subplot(1,3,3); + scatter(Vpiso,dtts./1000);ylim([0,4]); + hold on; + stairs(Vsiso,dtts./1000);ylim([0,4]); + ylabel('depth'); + xlabel('Velocity'); + set(gca,'YDir','Reverse'); + hold off; % The effective splitting calculation @@ -317,7 +353,7 @@ % EBSD data from Obbard and Baker (2007). Depths to deeper layers from Lipenkov and Barkov (1998) data(1) = struct('tex_file', fullfile(dat_dir, 'V97.txt'), ... % R-A 'dtb', 200, ... % depth to base of layer m - 'Tav', -25, ... % Average temperature (deg. C) + 'Tav', -24, ... % Average temperature (deg. C) 'azi', 0); % Horizonal rotation for texture data(2) = struct('tex_file', fullfile(dat_dir, 'V248.txt'), ... % R-B 'dtb', 300, ... @@ -393,7 +429,7 @@ 'azi', 0); data(20) = struct('tex_file', fullfile(dat_dir, 'V3321c.txt'), ... % ? 'dtb', 3750, ... - 'Tav', -1, ... + 'Tav', 0, ... 'azi', 0); end @@ -411,9 +447,9 @@ data = fscanf(fid, '%f', [12 inf]); nxtl = length(data(1,:)); eulers = zeros(3,nxtl); - eulers(1,:) = data(3,:); % phi1 - eulers(2,:) = data(4,:); % phi1 - eulers(3,:) = data(5,:); % phi1 + eulers(1,:) = data(5,:); % phi1 + eulers(2,:) = data(6,:); % phi1 + eulers(3,:) = data(7,:); % phi1 fclose(fid); end @@ -423,7 +459,7 @@ function report_layer(layernum, filename, Cvrh, rho, dtb, dtt, tav, ... fprintf('\nLayer: %i\n', layernum); fprintf(['data file: %s, \ntop: %f m, base: %f m, \n'... - 'thickness: %f m, teperature: %f C\n'], ... + 'thickness: %f m, temperature: %f C\n'], ... filename, dtt, dtb, dtb-dtt, tav); if plot_pole From 165d0ee417f6bc17bb59add3f97169d0987e067a Mon Sep 17 00:00:00 2001 From: Andrew Walker Date: Thu, 24 Jul 2014 10:59:07 +0100 Subject: [PATCH 03/12] Fix option in MS_plot Should be 'quiet' not 'silent' --- examples/glacial_splitting/glacial_splitting.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index 680900a..b8cfeaa 100644 --- a/examples/glacial_splitting/glacial_splitting.m +++ b/examples/glacial_splitting/glacial_splitting.m @@ -465,7 +465,7 @@ function report_layer(layernum, filename, Cvrh, rho, dtb, dtt, tav, ... if plot_pole MS_plot(Cvrh, rho, 'wtitle', filename, 'fontsize', 11, ... 'avscontours', 0:0.2:10.16, 'pcontours', 3.60:0.01:3.80, ... - 'polsize', 0.18, 0.16, 2.0, 1.0, 'limitsonpol', 'silent'); + 'polsize', 0.18, 0.16, 2.0, 1.0, 'limitsonpol', 'quiet'); end end From 95f7445feb6a74917590716fb8ad868fd8090d70 Mon Sep 17 00:00:00 2001 From: Andrew Walker Date: Thu, 24 Jul 2014 11:25:25 +0100 Subject: [PATCH 04/12] Fix calculation of isorropic velocity Rather than calculating the isotropic velcoty from the isotropic component of the poly-xtal anisotropic tensor (the norm of which changes with texture) use simple VRH mean of the single xtal tensor. This gives values that are just a function of the temperature and do not vary with texture. Also split the velocity graph into two graphs. --- .../glacial_splitting/glacial_splitting.m | 43 +++++++++---------- 1 file changed, 21 insertions(+), 22 deletions(-) diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index b8cfeaa..6446cfd 100644 --- a/examples/glacial_splitting/glacial_splitting.m +++ b/examples/glacial_splitting/glacial_splitting.m @@ -145,9 +145,6 @@ dtbs(j) = all_data(i).dtb; temps(j) = all_data(i).Tav; - - - % Single crystal elasticity of this layer [C, rho] = ice_cij(all_data(i).Tav); @@ -155,47 +152,49 @@ [Cs(:,:,j), rhos(j)] = Cs_from_EBSD_file(C, rho, ... all_data(i).tex_file); - % Calculate average isotropic velocities - CR = MS_axes(Cs(:,:,j), 'nowarn') ; - [Ciso] = MS_decomp(CR) ; - Vpiso(j) = 1e-3.*sqrt((Ciso(3,3)*1e9)./rhos(j)); - Vsiso(j) = 1e-3.*sqrt((Ciso(4,4)*1e9)./rhos(j)); + % Calculate average isotropic (VRH) velocity for this layer + [ K_vrh, G_vrh ] = MS_polyaverage( C ); + Vpiso(j) = 1e-3.*sqrt(((K_vrh*1e9)+((4.0/3.0)*G_vrh)*1e9)./rhos(j)); + Vsiso(j) = 1e-3.*sqrt((G_vrh*1e9)./rhos(j)); + % We may need the upper and lower bounds (V and R, respectivly) report_layer(i, all_data(i).tex_file, Cs(:,:,j), rhos(j), ... all_data(i).dtb, dtt, all_data(i).Tav, plot_pole) end - - Vpiso - Vsiso + % Plot stratigraphic column figure(); - subplot(1,3,1); + subplot(1,4,1); bar([fliplr(thickness);nan(size(thickness))],'stack');xlim([0.5,1.5]);ylim([0,4]); set(gca,'YDir','Reverse'); hold on; samples=[97 248 399 650 1201 1500 1849 2082 2749 2874 3300 3311 3321 3329 3399 3416]./1000.0; %scatter(ones(size(samples)),samples,'MarkerEdgeColor','k','MarkerFaceColor','w') scatter(linspace(0.7,1.3,max(size(samples))),samples,'MarkerEdgeColor','k','MarkerFaceColor','w') - ylabel('depth') + ylabel('Depth (km)') hold off; - subplot(1,3,2); + subplot(1,4,2); stairs(temps,dtts./1000);xlim([-25,0]);ylim([0,4]); - ylabel('depth'); - xlabel('Temperature'); + ylabel('Depth (km)'); + xlabel('Temperature (C)'); set(gca,'YDir','Reverse'); - subplot(1,3,3); + subplot(1,4,3); scatter(Vpiso,dtts./1000);ylim([0,4]); - hold on; - stairs(Vsiso,dtts./1000);ylim([0,4]); - ylabel('depth'); - xlabel('Velocity'); + ylabel('Depth (km)'); + xlabel('Isotropic p-wave velocity (km/s)'); set(gca,'YDir','Reverse'); - hold off; + + subplot(1,4,4); + scatter(Vsiso,dtts./1000);ylim([0,4]); + ylabel('Depth (km)'); + xlabel('Isotropic s-wave velocity (km/s)'); + set(gca,'YDir','Reverse'); + % The effective splitting calculation fprintf('\nEffective splitting calculation \n'); From 6d20dcadf716625eca5457abbc30d92f1f0df003 Mon Sep 17 00:00:00 2001 From: Andrew Walker Date: Thu, 24 Jul 2014 13:17:27 +0100 Subject: [PATCH 05/12] Handle both the new and old file formats Choose the format based on the first word in the first line of the EBSD files. Data is the same in both, but colomn order changes. --- .../glacial_splitting/glacial_splitting.m | 29 ++++++++++++++----- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index 6446cfd..edc80ac 100644 --- a/examples/glacial_splitting/glacial_splitting.m +++ b/examples/glacial_splitting/glacial_splitting.m @@ -442,13 +442,28 @@ fid = fopen(filename); % Read - the default assert((fid~=-1), 'Could not open file %s', filename); - fgetl(fid); % Header line - ignore - data = fscanf(fid, '%f', [12 inf]); - nxtl = length(data(1,:)); - eulers = zeros(3,nxtl); - eulers(1,:) = data(5,:); % phi1 - eulers(2,:) = data(6,:); % phi1 - eulers(3,:) = data(7,:); % phi1 + head = fgetl(fid); % Header line - use to distiguish the two + % kinds of file (ID as the first word or Index) + if all(head(1:6) == ' ID') + % First data format - from Geoff + data = fscanf(fid, '%f', [12 inf]); + nxtl = length(data(1,:)); + eulers = zeros(3,nxtl); + eulers(1,:) = data(3,:); % phi1 + eulers(2,:) = data(4,:); % phi2 + eulers(3,:) = data(5,:); % phi3 + elseif all(head(1:5) == 'Index') + % New data format - from Rachel + data = fscanf(fid, '%f', [12 inf]); + nxtl = length(data(1,:)); + eulers = zeros(3,nxtl); + eulers(1,:) = data(5,:); % euler1 + eulers(2,:) = data(6,:); % euler2 + eulers(3,:) = data(7,:); % euler3 + else + % Don't know what to do + error('Could not identify the EBSD data format from the 1st line') + end fclose(fid); end From 2d89bc3f22777d8ca02955be07e52a70a058fc4d Mon Sep 17 00:00:00 2001 From: Andrew Walker Date: Thu, 24 Jul 2014 14:43:31 +0100 Subject: [PATCH 06/12] Deal with non-vertical rays These are still straight but have a inclination that is not 90 degrees. Plot as a function of source polrization (otherwise we end up with chaos as we sweep around the azimuth) and inclnation for different azimuths. Some cleanup is needed. --- .../glacial_splitting/glacial_splitting.m | 132 ++++++++++++++++-- 1 file changed, 123 insertions(+), 9 deletions(-) diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index edc80ac..cf9e73c 100644 --- a/examples/glacial_splitting/glacial_splitting.m +++ b/examples/glacial_splitting/glacial_splitting.m @@ -196,7 +196,7 @@ set(gca,'YDir','Reverse'); - % The effective splitting calculation + % The effective splitting calculation for vertical ray path fprintf('\nEffective splitting calculation \n'); fprintf('--------------------------------\n'); % Set up inclination and azimuth and effective splitting @@ -222,11 +222,8 @@ [fast_eff(f,s), tlag_eff(f,s)] = do_effective_splitting(Cs, ... rhos, thickness, inc, azi, freq(f), spol(s), ... - plot_waves, eff_split_mode, quiet); - - % FIXME: do we need to correct fast_eff here? We are working in - % ray frame at the momenet. - % fast(j) = MS_unwind_pm_90((azi+pol')) ; % geog. ref frame + plot_waves, eff_split_mode, quiet); + if ~quiet fprintf('Effective fast direction: %f (deg)\n', fast_eff(f,s)); fprintf('Effective delay time: %f (s)\n', tlag_eff(f,s)); @@ -256,6 +253,114 @@ legend(freq_names) xlabel('Initial polarisation (deg)') ylabel('Fast directions (deg)') + + + % The effective splitting calculation for inclined rays + fprintf('\nEffective splitting calculation - inclined rays\n'); + fprintf('------------------------------------------------\n'); + % Set up inclination and azimuth and effective splitting + % params. These two could be functions of spol, for eg. + inc = [60 70 80] ; + azi = [0 45 90] ; + % Loop over source polarizations + % and frequencies + spol = 0:1:180; % Deg + freq = 30.0; % Hz - no point looking lower for SWS + inc_c = ['r', 'g', 'b']; + inc_names = {'60 degrees', '70 degrees', '80 degrees'}; + fast_eff_2 = zeros(length(spol), length(inc), length(azi)); + tlag_eff_2 = zeros(length(spol), length(inc), length(azi)); + for i = 1:length(inc) + for a = 1:length(azi) + for s = 1:length(spol) + + if ~quiet + fprintf('\n'); + fprintf('For source polarization: %f (deg)\n', spol(s)); + fprintf('azimuth: %f (deg)\n', azi(a)); + fprintf('inclination: %f (deg)\n', inc(i)); + end + + % Recalculate thickness etc + + [fast_eff_2(s,i,a), tlag_eff_2(s,i,a)] = do_effective_splitting(Cs, ... + rhos, thickness, inc(i), azi(a), freq, spol(s), ... + plot_waves, eff_split_mode, quiet); + + + if ~quiet + fprintf('Effective fast direction geographical frame: %f (deg)\n', fast_eff_2(s,i,a)); + fprintf('Effective delay time: %f (s)\n', tlag_eff_2(s,i,a)); + end + end + end + end + + % Make a figure to plot the results + scrsz = get(0,'ScreenSize'); + figure('Position',[1 scrsz(4) scrsz(3)/3 scrsz(4)*0.9]) ; + % Plot lag times + subplot(2,3,1) + for i = 1:length(inc) + plot(spol, tlag_eff_2(:,i,1), inc_c(i)) + hold on + end + xlabel('Source polarization (deg)') + ylabel('Lag time (s)') + thetitle = sprintf('Splitting as a function of source polarization and inclination, azi = %f (deg)', azi(1)); + title(thetitle); + + % Plot fast directions + subplot(2,3,4) + for i = 1:length(inc) + plot(spol ,fast_eff_2(:,i,1), inc_c(i)) + hold on + end + legend(inc_names) + xlabel('Source polarization (deg)') + ylabel('Fast directions (deg)') + + + subplot(2,3,2) + for i = 1:length(inc) + plot(spol, tlag_eff_2(:,i,2), inc_c(i)) + hold on + end + xlabel('Source polarization (deg)') + ylabel('Lag time (s)') + thetitle = sprintf('Splitting as a function of source polarization and inclination, azi = %f (deg)', azi(2)); + title(thetitle); + + % Plot fast directions + subplot(2,3,5) + for i = 1:length(inc) + plot(spol ,fast_eff_2(:,i,2), inc_c(i)) + hold on + end + legend(inc_names) + xlabel('Source polarization (deg)') + ylabel('Fast directions (deg)') + + subplot(2,3,3) + for i = 1:length(inc) + plot(spol, tlag_eff_2(:,i,3), inc_c(i)) + hold on + end + xlabel('Source polarization (deg)') + ylabel('Lag time (s)') + thetitle = sprintf('Splitting as a function of source polarization and inclination, azi = %f (deg)', azi(3)); + title(thetitle); + + % Plot fast directions + subplot(2,3,6) + for i = 1:length(inc) + plot(spol ,fast_eff_2(:,i,3), inc_c(i)) + hold on + end + legend(inc_names) + xlabel('Source polarization (deg)') + ylabel('Fast directions (deg)') + end @@ -265,7 +370,7 @@ % Header line for table... if ~quiet - fprintf('time lag (s) fast direction (deg)\n'); + fprintf('time lag (s) fast direction in ray frame (deg)\n'); end % Loop over layers and calculate splitting parameters fast = zeros(1,length(rhos)); @@ -275,8 +380,13 @@ % Calculate the splitting parameters for this layer [ pol, ~, vs1, vs2, ~, ~, ~ ] = MS_phasevels( Cs(:,:,j), ... rhos(j), inc, azi ); - fast(j) = pol; % FIXME: do we need to convert to geog ref? - tlag(j) = thickness(j)/vs2 - thickness(j)/vs1 ; + fast(j) = pol; + % Calculate tlag, ajusting for inclination (vertical ray has + % inc=90, so this is just length = thickness/1. Gets thicker + % as ray becomes shallower + tlag(j) = (thickness(j)/cosd(90-inc))/vs2 - ... + (thickness(j)/cosd(90-inc))/vs1 ; + if ~quiet fprintf(' %7f %7f \n', tlag(j), fast(j)); end @@ -290,6 +400,10 @@ [fast_eff,tlag_eff]=MS_effective_splitting_N(freq,spol,fast,... tlag,'mode',eff_split_mode); end + + % Put fast orentation into geographic reference frame. For + % Vertical ray azi = 0, so no change. + fast_eff = MS_unwind_pm_90(azi+fast_eff); end From d8be9328b2261143112472ab743ea28c77aac2d3 Mon Sep 17 00:00:00 2001 From: Andrew Walker Date: Mon, 28 Jul 2014 13:04:36 +0100 Subject: [PATCH 07/12] Allow MS_save to write to an already open file If 'fname' is a real (rather than a string) write to that assuming it is an fid (and avoid opening a new file). In this case the file is not closed prior to exiting MS_save. This can be used when either several Cijs must be written to the same file or to create parts of larger reports. Note that 'fname' = 1 will write data to the terminal. --- msat/MS_save.m | 42 ++++++++++++++++++++++++++++++++++++++---- 1 file changed, 38 insertions(+), 4 deletions(-) diff --git a/msat/MS_save.m b/msat/MS_save.m index 566aa1b..26a3973 100644 --- a/msat/MS_save.m +++ b/msat/MS_save.m @@ -2,7 +2,7 @@ % % // Part of MSAT - The Matlab Seismic Anisotropy Toolkit // % -% MS_save( fname, C, rho, ... ) +% MS_save( fname, C, rho, ... ) or MS_save( fid, C, rho, ... ) % % Usage: % MS_save(fname, C, rho) @@ -10,6 +10,16 @@ % 'fname'. C must be provided in GPa and rho must be provided in % kg/m3. See file format below. % +% +% MS_save(fid, C, rho) +% Write the elastic constants, C, and density, rho, to already +% open file described by integer valued double 'fid'. Fid is +% usually returned by a previous call to fopen. The file is not +% closed on exit from this function and no checks are made that +% the file is opened and can be used for output. C must be +% provided in GPa and rho must be provided in kg/m3. +% See file format below. +% % MS_save(fname, C, rho, ..., 'format',fmt) % Specify format for file. See below for descriptions. Available % options are: @@ -205,7 +215,7 @@ function MS_write_simple( fname, C, rho) - fid = fopen(fname, 'wt'); % Should we use t here? + [fid, fid_state] = MS_open_write_fh(fname); for i = 1:6 for j = i:6 @@ -214,14 +224,38 @@ function MS_write_simple( fname, C, rho) end fprintf(fid, '%1i %1i %f\n', 7, 7, rho); - fclose(fid); + MS_close_write_fh(fid, fid_state); end function MS_write_ematrix( fname, C, rho) error('MS:SAVE:NotImplemented', 'Writing to ematrix format is not implemented') end - + +function [fid, fid_state] = MS_open_write_fh(fname) + % Open file handle and remember if it needs to be closed + + if ischar(fname) + % This is a string - open it and return fid + fid = fopen(fname, 'wt'); % Should we use t here? + fid_state = 1; % We will need to close this fif + elseif isreal(fname) + % Not a string, but a number, so assume this is an fid + % we should use, but not close later + fid = fname; + fid_state = 0; % We will not need to close this one + else + % This isn't going to work + error('Cannot open file for write\n'); + end +end + +function MS_close_write_fh(fid, fid_state) + % If *we* opened the file handle, close it. + if fid_state + fclose(fid); + end +end From 952a2f55ab933b6ae0e91a543f8c16599e35093c Mon Sep 17 00:00:00 2001 From: Andrew Walker Date: Mon, 28 Jul 2014 13:08:48 +0100 Subject: [PATCH 08/12] Teach the glacial_splitting example to write ATRAK Create an input file for MODSIMPLE - a model builder for the ATRAK anisotropic ray tracer. --- .../glacial_splitting/glacial_splitting.m | 38 +++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index cf9e73c..8c448d5 100644 --- a/examples/glacial_splitting/glacial_splitting.m +++ b/examples/glacial_splitting/glacial_splitting.m @@ -195,6 +195,10 @@ xlabel('Isotropic s-wave velocity (km/s)'); set(gca,'YDir','Reverse'); + + % Create model for ATRAK + create_modsimple_rsp('modsimple.rsp', dtts, dtbs, Cs, rhos) + % The effective splitting calculation for vertical ray path fprintf('\nEffective splitting calculation \n'); @@ -597,3 +601,37 @@ function report_layer(layernum, filename, Cvrh, rho, dtb, dtt, tav, ... end end + +function create_modsimple_rsp(filename, dtts, dtbs, Cs, rhos) + + model_name = 'icesheet.mod'; + assert(length(dtts)==length(dtbs), 'Tops and bottoms must agree') + assert(length(dtts)==length(rhos), 'Tops and density must agree') + assert(length(dtts)==length(Cs), 'Tops and Cijs must agree') + + num_layers = length(dtts); + + fid = fopen(filename, 'w'); + fprintf(fid, '%s\n', model_name); + fprintf(fid, '%d\n', num_layers); + + % Hard code X, Y and Z nodes here. We could make this + % an argument, but I don't know what the rules are. + fprintf(fid, '%d %f %f\n', 10, 0.0, 10000.0); + fprintf(fid, '%d %f %f\n', 10, -5000.0, -5000.0); + fprintf(fid, '%d %f %f\n', 10, 0.0, -dtbs(num_layers)); % Base of hole + + % Now add each layer in turn (remembering we need to go backwards) + j = 0; + for i = num_layers:-1:1 + j = j + 1; + fprintf(fid, '%d %s\n', j, 'flat'); + fprintf(fid, '%f %f\n', -dtts(i), -dtts(i)); + fprintf(fid, '%s\n', 'ani'); + fprintf(fid, '%d\n', 0); + fprintf(fid, '%d\n', 21); + MS_save(fid, Cs(:,:,i), rhos(i), 'eunit', 'mbar', 'dunit', 'gcc', 'Aij'); + end + fclose(fid); +end + From ed061265fa2b63758fa36aef16d14ecbdce5984f Mon Sep 17 00:00:00 2001 From: Andrew Walker Date: Fri, 18 Sep 2015 12:03:59 +0100 Subject: [PATCH 09/12] Fix units for modsimple Aparanlty we had the wrong units for modsimple. According to Mike this is correct. --- examples/glacial_splitting/glacial_splitting.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index 8c448d5..2c3ead9 100644 --- a/examples/glacial_splitting/glacial_splitting.m +++ b/examples/glacial_splitting/glacial_splitting.m @@ -630,7 +630,7 @@ function create_modsimple_rsp(filename, dtts, dtbs, Cs, rhos) fprintf(fid, '%s\n', 'ani'); fprintf(fid, '%d\n', 0); fprintf(fid, '%d\n', 21); - MS_save(fid, Cs(:,:,i), rhos(i), 'eunit', 'mbar', 'dunit', 'gcc', 'Aij'); + MS_save(fid, Cs(:,:,i), rhos(i), 'eunit', 'pa', 'dunit', 'kgm3', 'Aij'); end fclose(fid); end From 05982c934d58448bf1ebb0a3e6f93129c636492c Mon Sep 17 00:00:00 2001 From: Andrew Walker Date: Mon, 28 Sep 2015 12:11:08 +0100 Subject: [PATCH 10/12] Anisotropies for vertical section Add these to the summary plot. We should probably refactor MS_plot and make this an MS_summary or something. --- .../glacial_splitting/glacial_splitting.m | 65 +++++++++++++++++-- 1 file changed, 58 insertions(+), 7 deletions(-) diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index 2c3ead9..d0e2896 100644 --- a/examples/glacial_splitting/glacial_splitting.m +++ b/examples/glacial_splitting/glacial_splitting.m @@ -21,7 +21,7 @@ % Note that setting this to 'GaussianWavelet' results in a much % longer execution time. -% (C) Alan Baird and Andrew Walker, 2014 +% (C) Alan Baird and Andrew Walker, 2014, 2015 % % Redistribution and use in source and binary forms, % with or without modification, are permitted provided @@ -126,6 +126,8 @@ ddbs = zeros(1,length(all_data)); Vpiso = zeros(1,length(all_data)); Vsiso = zeros(1,length(all_data)); + AVp = zeros(1,length(all_data)); + AVs_max = zeros(1,length(all_data)); j = 0; fprintf('\nLayer data (bottom to top) \n'); fprintf('---------------------------\n'); @@ -158,6 +160,9 @@ Vsiso(j) = 1e-3.*sqrt((G_vrh*1e9)./rhos(j)); % We may need the upper and lower bounds (V and R, respectivly) + % Anisotropy data for plot + [AVp(j), AVs_max(j)] = anisotropic_data(Cs(:,:,j), rhos(j)); + report_layer(i, all_data(i).tex_file, Cs(:,:,j), rhos(j), ... all_data(i).dtb, dtt, all_data(i).Tav, plot_pole) end @@ -166,7 +171,7 @@ % Plot stratigraphic column figure(); - subplot(1,4,1); + subplot(1,6,1); bar([fliplr(thickness);nan(size(thickness))],'stack');xlim([0.5,1.5]);ylim([0,4]); set(gca,'YDir','Reverse'); hold on; @@ -176,25 +181,35 @@ ylabel('Depth (km)') hold off; - subplot(1,4,2); + subplot(1,6,2); stairs(temps,dtts./1000);xlim([-25,0]);ylim([0,4]); ylabel('Depth (km)'); xlabel('Temperature (C)'); set(gca,'YDir','Reverse'); - subplot(1,4,3); + subplot(1,6,3); scatter(Vpiso,dtts./1000);ylim([0,4]); ylabel('Depth (km)'); xlabel('Isotropic p-wave velocity (km/s)'); set(gca,'YDir','Reverse'); + subplot(1,6,4); + scatter(Vpiso,dtts./1000);ylim([0,4]); + ylabel('Depth (km)'); + xlabel('Isotropic p-wave velocity (km/s)'); + set(gca,'YDir','Reverse'); - subplot(1,4,4); - scatter(Vsiso,dtts./1000);ylim([0,4]); + subplot(1,6,5); + scatter(AVp,dtts./1000);ylim([0,4]); ylabel('Depth (km)'); - xlabel('Isotropic s-wave velocity (km/s)'); + xlabel('P-wave anisotropy (%)'); set(gca,'YDir','Reverse'); + subplot(1,6,6); + scatter(AVs_max,dtts./1000);ylim([0,4]); + ylabel('Depth (km)'); + xlabel('Maximum s-wave anisotropy (%)'); + set(gca,'YDir','Reverse'); % Create model for ATRAK create_modsimple_rsp('modsimple.rsp', dtts, dtbs, Cs, rhos) @@ -635,3 +650,39 @@ function create_modsimple_rsp(filename, dtts, dtbs, Cs, rhos) fclose(fid); end +function [AVp, AVS_max] = anisotropic_data(C, rh) + % This is basically copied from MS_plot. We should + % pull this out into MS_something. + + grid_spacing = 6; + + % Set up inc-az grids... + [INC,AZ] = meshgrid( 90:-grid_spacing:0 , 0:grid_spacing:360 ) ; + + % Invoke MS_phasevels to get wave velocities etc. + [~,~,vs1,vs2,vp, S1P] = MS_phasevels(C,rh,... + reshape(INC, ((360/grid_spacing+1)*(90/grid_spacing+1)),1), ... + reshape(AZ, ((360/grid_spacing+1)*(90/grid_spacing+1)),1)); + + % reverse so sph2cart() works properly + % AZ = -AZ; + + % Reshape results back to grids + VS1 = reshape(vs1,61,16); + VS2 = reshape(vs2,61,16); + VP = reshape(vp,61,16); + % VS1_x = reshape(S1P(:,1),61,16); + % VS1_y = reshape(S1P(:,2),61,16); + % VS1_z = reshape(S1P(:,3),61,16); + + % aVS data + dVS = (VS1-VS2) ; + VSmean = (VS1+VS2)./2.0 ; + AVS = 100.0*(dVS./VSmean) ; + + % avp data + AVS_max = max(max(AVS)); + AVp = (max(VP)- min(VP))/(max(VP)+min(VP)).*200.0; + +end + From 995c205f199617a19af8605ace10a12ebc12fe2c Mon Sep 17 00:00:00 2001 From: Andrew Walker Date: Tue, 27 Oct 2015 15:53:47 +0000 Subject: [PATCH 11/12] Change frequencies in glacial spliting example These should correspond to SKS, local events and icequakes. --- examples/glacial_splitting/glacial_splitting.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index d0e2896..5689978 100644 --- a/examples/glacial_splitting/glacial_splitting.m +++ b/examples/glacial_splitting/glacial_splitting.m @@ -225,9 +225,9 @@ % Loop over source polarizations % and frequencies spol = min_azi:del_azi:max_azi; % Deg - freq = [0.3, 3.0, 30.0]; % Hz + freq = [3, 30, 3000]; % Hz freq_c = ['r', 'g', 'b']; - freq_names = {'0.3 Hz', '3.0 Hz', '30.0 Hz'}; + freq_names = {'3 Hz', '30 Hz', '3000 Hz'}; fast_eff = zeros(length(freq),length(spol)); tlag_eff = zeros(length(freq),length(spol)); for f = 1:length(freq) From 63fc8c9b3edf109551623573cd2f428db2abca5c Mon Sep 17 00:00:00 2001 From: Andrew Walker Date: Fri, 9 Dec 2016 14:18:29 +0000 Subject: [PATCH 12/12] Fix modsimple generator for glacial splitting We do not need a depth for the top of the model - so just print the depth, not zero. --- examples/glacial_splitting/glacial_splitting.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index 5689978..97c0f69 100644 --- a/examples/glacial_splitting/glacial_splitting.m +++ b/examples/glacial_splitting/glacial_splitting.m @@ -632,9 +632,10 @@ function create_modsimple_rsp(filename, dtts, dtbs, Cs, rhos) % Hard code X, Y and Z nodes here. We could make this % an argument, but I don't know what the rules are. + % (10 is the discritisation of the model) fprintf(fid, '%d %f %f\n', 10, 0.0, 10000.0); fprintf(fid, '%d %f %f\n', 10, -5000.0, -5000.0); - fprintf(fid, '%d %f %f\n', 10, 0.0, -dtbs(num_layers)); % Base of hole + fprintf(fid, '%d %f\n', 10, -dtbs(num_layers)); % Base of hole % Now add each layer in turn (remembering we need to go backwards) j = 0;