diff --git a/examples/glacial_splitting/glacial_splitting.m b/examples/glacial_splitting/glacial_splitting.m index 232db01..97c0f69 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 @@ -122,6 +122,12 @@ 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)); + AVp = zeros(1,length(all_data)); + AVs_max = zeros(1,length(all_data)); j = 0; fprintf('\nLayer data (bottom to top) \n'); fprintf('---------------------------\n'); @@ -137,18 +143,79 @@ 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 (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) + % 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 - % The effective splitting calculation + + % Plot stratigraphic column + + figure(); + 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; + 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 (km)') + hold off; + + 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,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,6,5); + scatter(AVp,dtts./1000);ylim([0,4]); + ylabel('Depth (km)'); + 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) + + + % The effective splitting calculation for vertical ray path fprintf('\nEffective splitting calculation \n'); fprintf('--------------------------------\n'); % Set up inclination and azimuth and effective splitting @@ -158,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) @@ -174,11 +241,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)); @@ -208,6 +272,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 @@ -217,7 +389,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)); @@ -227,8 +399,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 @@ -242,6 +419,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 @@ -301,40 +482,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 - 'Tav', -25, ... % Average temperature (deg. C) + % 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', -24, ... % 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(3) = struct('tex_file', fullfile(dat_dir, 'V3321c.txt'), ... - 'dtb', 2900, ... - 'Tav', -15, ... + 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(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(17) = struct('tex_file', fullfile(dat_dir, 'V3311g.txt'), ... % A4 + 'dtb', 3460, ... + 'Tav', -1.9, ... 'azi', 0); - data(8) = struct('tex_file', fullfile(dat_dir, 'V3311g.txt'), ... - 'dtb', 3450, ... + 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', 0, ... + 'azi', 0); end @@ -347,13 +575,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(3,:); % phi1 - eulers(2,:) = data(4,:); % phi1 - eulers(3,:) = data(5,:); % 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 @@ -363,13 +606,84 @@ 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 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 + + +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. + % (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\n', 10, -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', 'pa', 'dunit', 'kgm3', 'Aij'); end + 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 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