📄 compute.m
字号:
% ylabel('Real(H_z)');
% subplot(2,1,2);
% semilogx(cparams.freq, imag(v), '-xb', 'Markersize',5);
% ylabel('Imag(H_z)');
% xlabel('Frequency (Hz)');
% end
% --------------------------------------------------------------------
% function plot_hcp(result, config, cparams, fig)
%
% if min(size(result)) == 1
% figure(fig, 'visible', 'on');
% for k = 1:length(result)
% subplot(2,1,1);
% semilogx(cparams.freq, real(result(k).H./result(k).H_prim));
% hold on;
% subplot(2,1,2);
% semilogx(cparams.freq, imag(result(k).H./result(k).H_prim));
% hold on;
% end % k
% subplot(2,1,1);
% ylabel('Real(Z/Z_0)');
% hold off;
% subplot(2,1,2);
% ylabel('Imag(Z/Z_0)');
% xlabel('Freq. (Hz)');
% hold off;
% else
% disp('I don'' know how to plot these results!');
% end % if
% --------------------------------------------------------------------
% function plot_emgsa(result, config, cparams, fig)
%
% if min(size(result)) == 1
% figure(fig);
% set(fig, 'visible', 'on');
% for k = 1:length(result)
% if ~isfield(result(k), 'G_factor')
% result(k).G_factor = 0;
% if isfield(result(k), 'Aspac')
% result(k).G_factor = pi.*result(k).Aspac ...
% .*(result(k).Nspac) ...
% .*(result(k).Nspac+1) ...
% .*(result(k).Nspac+2);
% end
% if isempty(result(k).G_factor) || result(k).G_factor == 0
% G_factor = 1;
% end % if
% end
% subplot(2,1,1);
% % halldat = - 6.7e-011*cparams.freq.^3 + ...
% % 1.1e-006*cparams.freq.^2 - 0.01*cparams.freq - 1;
% semilogx(cparams.freq, abs(result(k).Z.*result(k).G_factor));
% hold on;
% subplot(2,1,2);
% semilogx(cparams.freq, ...
% -angle(result(k).Z).*1000);%+halldat.');
% hold on;
% end % k
% subplot(2,1,1);
% ylabel('\rho_a (\Omega\cdot m)');
% hold off;
% subplot(2,1,2);
% ylabel('neg. phase (mrad)');
% xlabel('Freq. (Hz)');
% hold off;
% else
% % if you want results in a 3D matrix for D-D data
% % Z = reshape([result.Z],length(result(1).Z), ...
% % size(result,1),size(result,2));
% disp('I don'' know how to plot these results!');
% end % if
% --------------------------------------------------------------------
% function plot_dcgsa(result, config, cparams, fig)
%
% if min(size(result)) == 1 && max(size(result)) > 1
% figure(fig, 'visible', 'on');
% yval = squeeze(reshape([result.Z], size(result,1), size(result,2)));
% G_factor = squeeze(reshape([result.G_factor], size(result,1), ...
% size(result,2)));
%
% for k = 1:length(yval)
% if G_factor(k) == 0
% G_factor(k) = 1;
% end % if
% switch config.type
% case {'Dipole-Dipole', '*Capacitance*'}
% if size(result,1) > 1
% xval = squeeze(reshape( ...
% [result.Aspac], ...
% size(result,1), 1));
% loglog(xval, yval.*G_factor);
% else
% xval = squeeze(reshape( ...
% [result.Nspac], ...
% 1, size(result,2)));
% semilogy(xval, yval.*G_factor);
% end % if
% case {'Wenner'}
% xval = squeeze(reshape( ...
% [result.Aspac], ...
% size(result,1), 1));
% semilogy(xval, yval.*G_factor);
% case {'Schlumberger'}
% if size(result,1) > 1
% xval = squeeze(reshape( ...
% [result.OA], ...
% size(result,1), 1));
% loglog(xval, yval.*G_factor);
% else
% xval = squeeze(reshape( ...
% [result.OM], ...
% 1, size(result,2)));
% semilogy(xval, yval.*G_factor);
% end % if
% end
% hold on;
% end % k
% hold off;
% ylabel('Apparent Resistivity (\Omegam)');
% switch config.type
% case {'Dipole-Dipole', '*Capacitance*'}
% if size(result,1) > 1
% xlabel('A-spacing (m)');
% else
% xlabel('N-spacing (m)');
% end % if
% case {'Wenner'}
% xlabel('A-spacing (m)');
% case {'Schlumberger'}
% if size(result,1) > 1
% xlabel('OA (m)');
% else
% xlabel('OM (m)');
% end % if
% end % switch
% end % if
% --------------------------------------------------------------------
function saveforward(model, varargin)
if getappdata(0, 'debug')
disp('compute.m:saveforwad:debug: Entered saveforward routine');
end
if nargin == 2
[pathname,filename,ext] = fileparts(varargin{1});
filterindex = [];
switch lower(ext)
case {'','.mat'}
filterindex = 1;
case {'.dat'}
filterindex = 2;
otherwise
filterindex = 3;
end
else
if getappdata(0, 'compiled')
% !!! for standalone
[filename, pathname] = uiputfile( ...
{'*.mat','MAT-files (*.mat)'; ...
'*.dat','ASCII file (*.dat)'; ...
'*.*', 'All Files (*.*)'}, ...
'Save result as...', 'model_result.mat');
if isequal(filename,0), return; end
filterindex = [];
[tmp1,filename,ext] = fileparts(filename);
switch lower(ext)
case {'','.mat'}
filterindex = 1;
case {'.dat'}
filterindex = 2;
otherwise
filterindex = 3;
end
else
% !!! for matlab mode
[filename, pathname, filterindex] = uiputfile( ...
{'*.mat','MAT-files (*.mat)'; ...
'*.dat','ASCII file (*.dat)'; ...
'*.*', 'All Files (*.*)'}, ...
'Save result as...', 'model_result.mat');
if isequal(filename,0), return; end
[tmp1,filename,ext] = fileparts(filename);
end
end
if filterindex == 1 || filterindex == 3
if isempty(ext), ext = '.mat'; end
save(fullfile(pathname, [filename ext] ),'model','-mat');
elseif filterindex==2
fid = fopen(fullfile(pathname, [filename ext]),'w');
fprintf(fid,'Configuration: %s\n',model.config.type);
switch model.config.type
case {'Dipole-Dipole', '*Capacitance*'}
fprintf(fid,'A-spacing: ');
fprintf(fid,'%f\t',model.config.Aspac);
fprintf(fid,'\nN-spacing: ');
fprintf(fid,'%f\t',model.config.Nspac);
case {'Schlumberger'}
fprintf(fid,'OA: %f\n',model.config.OA);
fprintf(fid,'OM: %f\n',model.config.OM);
case {'Wenner'}
fprintf(fid,'A-spacing: %f\n',model.config.a);
case {'HCP FDEM (HLEM)'}
fprintf(fid,'R-spacing: %f\n',model.config.Rspac);
case {'TEM Central Loop'}
fprintf(fid,'Tx-side: %f\n',model.config.TxS);
fprintf(fid,'Rx-radius: %f\n',model.config.RxR);
end
fprintf(fid,'\nModel:\n');
fprintf(fid,'Rho DC: ');
fprintf(fid,'%f\t',[model.layers.rho]);
fprintf(fid,'\nThickness: ');
fprintf(fid,'%f\t',[model.layers.thickness]);
fprintf(fid,'\ndepth_to_top: ');
fprintf(fid,'%f\t',[model.layers.depth_to_top]);
fprintf(fid,'\n\nCalculation: %s\n',model.calc_type);
fprintf(fid,'Transform type: %s\n',model.hank_type);
if strcmp(model.hank_type,'cconvol')
fprintf(fid,'cconvol tol: %s\n',model.cconvolerr);
else
fprintf(fid,'Hank tol: %s\n',model.Hank_tol);
end
fprintf(fid,'\nResult:\n');
switch model.calc_type
case{'DC'}
switch model.config.type
case {'Dipole-Dipole', '*Capacitance*'}
fprintf(fid,'Aspac\t\tNspac\t\tRhoA\n');
for k = 1:length(model.result)
fprintf(fid,'%f\t%f\t%f', ...
model.result(k).Aspac, ...
model.result(k).Nspac, ...
model.result(k).Z);
end
case {'Schlumberger'}
fprintf(fid,'OA\t\tOM\t\tRhoA\n');
for k = 1:length(model.result)
fprintf(fid,'%f\t%f\t%f', ...
model.result(k).OA, ...
model.result(k).OM, ...
model.result(k).Z);
end
case {'Wenner'}
for k = 1:length(model.result)
fprintf(fid,'%f\t%f', ...
model.result(k).Aspac, ...
model.result(k).Z);
end
case {'General Surface Array','*GSA capacitance*'}
for k = 1:length(model.result)
fprintf(fid,'%i\t%f', ...
k, model.result(k).Z);
end
end
case {'FD'}
switch model.config.type
case {'HCP FDEM (HLEM)'}
for k = 1:length(model.result)
fprintf(fid,'%f\t%f', ...
model.result(k).Rspac, ...
model.result(k).Z);
end
end
otherwise
disp('I don''t know how to save these data...!');
end
fclose(fid);
end
% --- Executes on button press in Reusefig_checkbox.function Reusefig_checkbox_Callback(hObject, eventdata, handles)% hObject handle to Reusefig_checkbox (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Hint: get(hObject,'Value') returns toggle state of Reusefig_checkbox
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -