temfwd.m

matlab中常用电法各种装置的电磁正演计算

function out = temfwd(varargin);

% out = temfwd(config, layers, cparams);
% 
% temfwd.m is part of the CR1Dmod forward modeling package, and contains
% the code used to calculate the response of a central-loop configuration
% over a 1D layered half space.
%
% Input:
% config :    structure defining the configuation (generated by cr1dmod) 
% layers :    structure defining layer-parameters (generated by cr1dmod)
% cparams:    structure defining calculation options (generated by compute)
%
% Output:
% out:        vector of calculated impedances
%
% The file contains an option to calculate the response of a square loop
% transmitter loop as measured by a vertical magnetic dipole receiver.
% So far only dipole location in the center of the transmitter loop is 
% allowed. The Sq-dip option is not included in the GUI.
%
% Written by:
% Thomas Ingeman-Nielsen
% The Arctic Technology Center, BYG
% Technical University of Denmark
% Email: [e-mail]

cpu_t=cputime;

global FDkern

if ( nargin < 3 )
    disp('A call to temfwd must provide three input arguments!');
    return
end

% Set up structures of parameters...
ok = 0;
if ( isstruct(varargin{1})         )
    config = varargin{1};
    if ( isstruct(varargin{2})     )
        layers = varargin{2};
        if ( isstruct(varargin{3}) )
            cparams = varargin{3};                    
            ok = 1;
        end
    end
end
if ~ok
    disp('Non-structure inputs have not been implemented yet!');
    return
end

config.Txz = 0; % assume Tx loop is on the ground (only option implemented)
config.Rxz = 0; % assume Rx loop is on the ground (only option implemented)
% The function is not capable of calculating with heights!

if ( ~isfield(config,'TEMtype') )  || ...
    ( isempty(config.TEMtype)   )
    config.TEMtype = 'circ-dip';       % Geometry assumption [circ-dip, sq-dip, sq-sq]
end

if ( ~isfield(config,'TxS')     )  || ...
        ( isempty(config.TxS)   )
    config.TxS = 40;
    config.TxR = sqrt(config.TxS.^2./pi);   
    disp('TxS set to 40m');
    beep;
end

if ( ~isfield(config,'TxR')     )  || ...
        ( isempty(config.TxR)   )
    config.TxR = sqrt(config.TxS.^2./pi);       % Tx equivalent circular loop radius used for normalization
end

% Prepares frequency domain spline or individual calculations by either
% calculating a spline function and returning a handle to a function that
% evaluates the spline, or returning a handle to the function that calculates
% the frequency domain response.
TEM_handle = Hz_TEM_driver(config, layers, cparams);

disp('Finished the frequency domain computations...');

switch  cparams.domain
    case {'TD'}
        if ( cparams.FDspline      )  && ...
                ( cparams.showFDsp )
            h1 = figure;
            semilogx(sqrt(FDkern.b), FDkern.kern(2:end-1), '-xb',       ...
                'markersize', 4);
            h2 = legend('Spline',3);
        end
        
        T = 2.*cparams.times./(1./layers(1).rho.*4.*pi.*1e-7.*          ...
            config.TxR.^2); % Define normalized time
        
        hw = waitbar(0,'Please wait...');
        for k=1:max(size(T))
            switch cparams.FtoTtype
                case {'FST'}
                    out(k) = -FCST('SI4', T(k), -1, 0, 0,            ...
                        cparams.FST_err, TEM_handle, config,            ...
                        layers, cparams);
                case {'NST'}                
                    out(k) = -NJCST('SIN', TEM_handle, T(k),             ...
                        cparams.Seg_tol, cparams.NST_tol,               ...
                        cparams.Max_seg, config, layers, cparams);
                case {'NCT'}
                    out(k) = NJCST('COS', TEM_handle, T(k),              ...
                        cparams.Seg_tol, cparams.NCT_tol,               ...
                        cparams.Max_seg, config, layers, cparams);
            end % switch
            waitbar(k/(length(T)),hw);
        end % k
        close(hw);
        
        switch config.TEMtype
            case {'circ-dip'}
                out = out.*2./pi.*config.RxA./                          ...
                    (1./layers(1).rho.*config.TxR.^3);
            case {'sq-dip'}
                out = out.*8.*config.Rxn.*config.RxR.^2./               ...
                    (config.TxR.^2.*1./layers(1).rho);
            case {'sq-sq'}
                % not implemented yet
        end        
        
        if ( cparams.FDspline      )  && ...
                ( cparams.showFDsp )  && ...
                ( ishandle(h1)     )
            figure(h1);
            hold on;
            est = unique(FDkern.estimate);
            semilogx(sqrt(est), ppval(FDkern,est), 'xr', 'markersize', 4);
            delete(h2);
            legend({'Spline'; 'Evaluations'},3);
            xlimits = xlim;
            omega = 2.*xlimits.^2./                                     ...
                (1./layers(1).rho.*4e-7.*pi.*config.TxR.^2);
            disp(omega);
            
        end
            
    case {'FD'}
        b = cparams.freq.*2.*pi./                                       ...
            layers(1).rho.*4e-7.*pi.*config.TxR.^2./2;
        out = feval(TEM_handle, b, config, layers, cparams);
end    

disp(['Time spent: ' num2str(cputime-cpu_t) ' s']);

% -------------------------------------------------------------------------
function funchan = Hz_TEM_driver(config, layers, cparams);

global FDkern estimate

if ( ~cparams.FDspline )     % In the case no spline interpolation is wanted
                             %     return handle to FD calculation routine
    estimate = [];                             
    switch cparams.domain
        case {'TD'}
            switch cparams.FtoTtype
                case {'NST','FST'}
                    switch config.TEMtype
                        case {'circ-dip'}
                            funchan = @imHz_CircDip;
                        case {'sq-dip'}
                            funchan = @imHz_SqDip;
                        case {'sq-sq'}
                            % not implemented yet
                    end                    
                case {'NCT','FCT'}
                    switch config.TEMtype
                        case {'circ-dip'}
                            funchan = @reHz_CircDip;
                        case {'sq-dip'}
                            funchan = @reHz_SqDip;
                        case {'sq-sq'}
                            % not implemented yet
                    end                  
            end % switch
        case {'FD'}
            switch config.TEMtype
                case {'circ-dip'}
                    funchan = @Hz_CircDip;
                case {'sq-dip'}
                    funchan = @Hz_SqDip;
                case {'sq-sq'}
                    % not implemented yet
            end
    end % switch
        
else                         % Else set up the spline function         
    if ( ~isfield(cparams,'FDsp_NDEC')   )  ||                          ...
            ( isempty(cparams.FDsp_NDEC) )
        cparams.FDsp_NDEC = 15;
    end
    
    if ( ~isfield(cparams,'FDsp_Bmin')   )  ||                          ...
            ( isempty(cparams.FDsp_Bmin) )
        cparams.FDsp_Bmin = 1e-3;
    end
    
    if ( ~isfield(cparams,'FDsp_Bmax')   )  ||                          ...
            ( isempty(cparams.FDsp_Bmax) )
        cparams.FDsp_Bmax = 1e3;
    end
    
    b = [10.^(2.*([log10(cparams.FDsp_Bmin):                            ...
                1./cparams.FDsp_NDEC:                                   ...
                log10(cparams.FDsp_Bmax)]))];
    
    switch config.TEMtype
        case {'circ-dip'}
            kern = Hz_CircDip(b, config, layers, cparams);
        case {'sq-dip'}
            kern = Hz_SqDip(b, config, layers, cparams);
        case {'sq-sq'}
            % not implemented yet
    end                    
    
    switch cparams.domain
        case {'TD'}
            switch cparams.FtoTtype
                case {'NST','FST'}
                    % Used for Sine transform
                    b = [0 b];
                    kern = [0 0 imag(kern) 0]; % Make the endslopes of the spline horizontal
                case {'NCT','FCT'}
                    % Used for Cosine transform
                    b = [0 b cparams.FDsp_Bmax.^2.*1e6];
                    kern = [0 1 real(kern) 0 0]; 
            end % switch
            
        case {'FD'}
            kern = [0 kern 0]; % Make the endslopes of the spline horizontal
    end % switch
    
    FDkern = spline(b,kern);
    FDkern.b = b;
    FDkern.kern = kern;
    FDkern.estimate = [];
    
    funchan = @FDspline_eval;
end

% ---------------------------------------------------------------------
function kern = FDspline_eval(b, config, layers, cparams);

global FDkern
kern = ppval(b, FDkern);
FDkern.estimate = [FDkern.estimate b];

% ---------------------------------------------------------------------
function kern = imHz_CircDip(b, config, layers, cparams);

global estimate

kern = imag(Hz_CircDip(b, config, layers, cparams));
estimate = [estimate, [b;kern]];

% ---------------------------------------------------------------------
function kern = imHz_SqDip(b, config, layers, cparams);

global estimate

kern = imag(Hz_SqDip(b, config, layers, cparams));
estimate = [estimate, [b;kern]];


% ---------------------------------------------------------------------
function kern = reHz_CircDip(b, config, layers, cparams);

global estimate

kern = real(Hz_CircDip(b, config, layers, cparams));
estimate = [estimate, [b;kern]];

% ---------------------------------------------------------------------
function kern = reHz_SqDip(b, config, layers, cparams);

global estimate

kern = real(Hz_SqDip(b, config, layers, cparams));
estimate = [estimate, [b;kern]];

% *************************************************************************
% * Code for circular transmitter loop and dipole reciever                *
% *************************************************************************
function [Hz] = Hz_CircDip(b, config, layers, cparams)

% [Hz] = Hz_CircDip(b, config, layers, cparams);
% Calculates the Hz field for the central loop configuration
% at a speciffic frequency.