⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 compute.m

📁 matlab中常用电法各种装置的电磁正演计算
💻 M
📖 第 1 页 / 共 5 页
字号:
%     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 + -