cr1dmod.m

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

function varargout = cr1dmod(varargin)

% CR1DMOD is the main function in a 1D forward modelling package which has
% the capability of modelling Complex Resistivity effects.
% A call to CR1DMOD will open the main gui interface, from which layered
% earth models can be created and measurement configurations defined.
%
% The program is written for and tested in Matlab 6.5 R13.
%
% The program is described in:
% Ingeman-Nielsen, T. and Baumgartner, F. (2005): CR1Dmod: A Matlab program
% to model 1D Complex Resistivity effects in electrical and electromagnetic
% surveys, submitted to Computers & Geosciences.
%
% The CR1DMOD package includes the files:
% cr1dmod.m:          This file
% cr1dmod.fig:        Gui layout of the main window
% topview.fig:        Gui layout of the topview window, which allows the
%                       user to place electrodes and wires arbitrarily on
%                       the surface of the half-space.
% topview.m:          Functions relating to the topview GUI
% compute.fig:        Gui layout of the compute window, from which
%                       calculation speciffic parameters are defined and the
%                       actual computation started.
% compute.m:          Functions relating to the compute GUI
% dcgsafwd.m:         Forward modelling code for the DC response of a
%                       general surface electrode array
% emgsafwd.m:         Forward modelling code for the EM response of a
%                       general surface electrode array
% temfwd.m:           Forward modelling code for the frequency and transient
%                       response of the central loop configuration
% fdemfwd.m:          Forward modelling code for the frequency response of a
%                       horizontal coplanar loop configuration
% NJCST.m:            Routine to perform numerical hankel and harmonic
%                       transforms
% zeros_J.mat:        Precalculated values of the zeros of the bessel
%                       functions used in the NJCST hankel transform
% FJCST.m:            Routine to perform hankel and harmonic transforms
%                       using the digital filter method
% FCST.m:             Similar to FJCST but calculates only harmonic
%                       transforms
% filters.mat:        Filter coefficients needed in FJCST and FCST
% Z_CR.m:             Function to calculate the resistivity dispersion based
%                       on the cole-cole model
% neg_loglog.m:       Function to create a loglog plot where positive data
%                       values plot in blue, and negative data values in red.
%
% Written by:
% Thomas Ingeman-Nielsen
% The Arctic Technology Center, BYG
% Technical University of Denmark
% Email: [e-mail]


% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @cr1dmod_OpeningFcn, ...
    'gui_OutputFcn',  @cr1dmod_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --------------------------------------------------------------------
% --- Executes just before cr1dmod is made visible.
function cr1dmod_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to cr1dmod (see VARARGIN)

% Choose default command line output for cr1dmod
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

if nargin > 3
    for k = 1:nargin-3
        if isstr(varargin{k}) && strcmpi(varargin{k},'debug')
            setappdata(0, 'debug', 1);  % 1: go through debug code
            disp('The code will display debug information!');
            break;
        end
    end
end     

if strcmp(get(hObject,'Visible'),'off')
    testVar = 1;
    if exist('testVar', 'var')
        % Executes when the code is running in the matlab environment
        setappdata(0, 'compiled', 0);  % 1: bypass non-compilable code
        handles.mfilePath = fileparts(mfilename('fullpath'));
        I = strmatch(handles.mfilePath,path);
        if isempty(I)
            addpath(handles.mfilePath,'-begin');
            handles.mfilePathRemove = 1;
        else
            handles.mfilePathRemove = 0;
        end
    else
        % Executes when the code is compiled        
        setappdata(0, 'compiled', 1);  % 1: bypass non-compilable code
        handles.mfilePathRemove = 0;
    end
    
    axes(handles.axes1);
    
    list = {'Dipole-Dipole',...
            ...%'*Capacitance*',...
            'Wenner',...
            'Schlumberger',...
            'General Surface Array',...
            ...%'*GSA capacitance*',...
            'TEM Central Loop',...
            'HCP FDEM (HLEM)'};
    %         'Pole-Dipole',...
    %         'Pole-Pole',...
    %         'Borehole',...
    %         'Free',...
    
    set(handles.Config_popup,'string',list);
    
    handles.config.type = 'Dipole-Dipole';
    handles.config.Aspac = 100;
    handles.config.Nspac = 1;
    handles.config.Rspac = 100;
    handles.config.OA = 100;
    handles.config.OM = 1;
    handles.config.Cwire = [];
    handles.config.Pwire = [];
    handles.config.TxS = 40;
    handles.config.RxA = 100;
    
    % Update handles structure
    guidata(hObject, handles);
    
    handles.layers(1).depth_to_top = 0;
    handles.layers(1).thickness = inf;
    handles.layers(1).rho = 100;
    handles.layers(1).m = 0;
    handles.layers(1).tau = 0;
    handles.layers(1).c = 0;
    handles.layers(1).mu = '0';   % NB this is the magnetic susceptibility!
    handles.layers(1).eps_r = 1;
    
    handles.model = [];
    
    set(handles.axes1,'ylim',[-50 20]);
    hold on;
    handles.surface = plot(get(handles.axes1,'xlim'),[0 0],'-k','linewidth',2);
    
    handles = update_config(handles);
    handles = update_model(handles);
    update_layer_param(handles);
    
    % Update handles structure
    guidata(hObject, handles);
    
end

% UIWAIT makes cr1dmod wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --------------------------------------------------------------------
% --- Outputs from this function are returned to the command line.
function varargout = cr1dmod_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% *****************************************************************
% * Callback functions related to model plot                      *
% *****************************************************************

% --------------------------------------------------------------------
% --- Executes on button press in axes1.
function axes1_Callback(hObject, eventdata, handles)

loc = get(hObject, 'currentpoint');
depth = round(loc(2,2)*10)/10;
if depth<0
    % find all interfaces above the one to be made
    index = find([handles.layers.depth_to_top]<-depth);    
    
    layer_num = index(end); 
    handles.layers(layer_num+1:end+1) = handles.layers(layer_num:end);
    if layer_num == 1
        handles.layers(layer_num).thickness = -depth;
    else
        handles.layers(layer_num).thickness =                           ...
            -depth-handles.layers(layer_num).depth_to_top;
    end
    handles.layers(layer_num+1).thickness =                             ...
        handles.layers(layer_num+1).thickness -                         ...
        handles.layers(layer_num).thickness;     
    handles.layers(layer_num+1).depth_to_top = -depth;
    handles.layers(layer_num).interface_handle = [];
    % ^ otherwise it points to the interface of layer_num+1!
    handles.layers(layer_num).prop_lab_handle = [];
    handles.layers(layer_num).label_handle = [];
    
    handles = update_model(handles);
    update_layer_param(handles);
    
    % Update handles structure
    guidata(hObject, handles);
    
    set(gcbf,'WindowButtonUpFcn',@interface_ButtonUp_Callback);
    set(gcbf,'WindowButtonMotionFcn',{@interface_Move_Callback, layer_num+1}); 
end


% --------------------------------------------------------------------
% --- Executes on button press on an interface.
function interface_ButtonDown_Callback(hObject, eventdata)

handles = guidata(hObject);
interface = find([handles.layers.interface_handle]==hObject)+1;
click_type = get(gcbf,'SelectionType');
switch click_type
    case{'normal'}   % left click
        set(gcbf,'WindowButtonUpFcn',@interface_ButtonUp_Callback);
        set(gcbf,'WindowButtonMotionFcn',                               ...
            {@interface_Move_Callback, interface});        
    case{'extend'}   % Shift - left
    case{'alt'}      % Ctrl - left    
    case{'open'}     % Double click
end


% --------------------------------------------------------------------
% --- Executes on button release on an interface.
function interface_ButtonUp_Callback(hObject, eventdata)

set(gcbf,'WindowButtonMotionFcn','');
handles = guidata(hObject);
ylimmits = ylim;

if length(handles.layers) > 1
    set(handles.axes1, 'ylim', ([-1 2/5] * (9/7 *                       ...
        handles.layers(end).depth_to_top)));
end

if strcmp(handles.config.type,'TEM Central Loop')
    yscale = diff(get(handles.axes1,'ylim'));
    t = sin([0:pi/100:pi]).*yscale./70;
    set(handles.config.plot_handle(1),'ydata',t.*4);
    set(handles.config.plot_handle(2),'ydata',t.*1);
elseif strcmp(handles.config.type,'HCP FDEM (HLEM)')
    yscale = diff(get(handles.axes1,'ylim'));
    t = sin([0:pi/100:2.*pi]').*yscale./70;
    set(handles.config.plot_handle(1:2), 'ydata', t, 'linewidth',       ...
        2, 'color', 'b');
    set(handles.config.plot_handle(3:4), 'ydata',                       ...
        [-yscale +yscale].*1.5./70,'linewidth',2,'color','b');
end

handles = update_model(handles);



% --------------------------------------------------------------------
% --- Executes on mouse movement with button pressed on an interface.
function interface_Move_Callback(hObject, eventdata, interface)

handles = guidata(hObject);

CurrentPoint = mean(get(gca,'CurrentPoint'));

layer_num = interface-1;

ylimmits = ylim;
if -CurrentPoint(2)<=handles.layers(layer_num).depth_to_top
    CurrentPoint(2) = -handles.layers(layer_num).depth_to_top;
elseif -CurrentPoint(2) >= handles.layers(layer_num+1).depth_to_top +   ...
        handles.layers(layer_num+1).thickness
    CurrentPoint(2) = -handles.layers(layer_num+1).depth_to_top -       ...
        handles.layers(layer_num+1).thickness;
elseif CurrentPoint(2)<ylimmits(1)
    temppoint = CurrentPoint;
    if isfield(handles,'lastpoint')
        if handles.lastpoint(2)>CurrentPoint(2)
            set(handles.axes1,'ylim',(ylimmits-[1 -2/5]));
            CurrentPoint(2) = ylimmits(1)-1;
        else 
            CurrentPoint(2) = ylimmits(1);
        end    
    else
        set(handles.axes1,'ylim',(ylimmits-[1 -2/5]));
        CurrentPoint(2) = ylimmits(1)-1;
    end
    handles.lastpoint = temppoint;
end

handles.layers(layer_num).thickness = round((-CurrentPoint(2) -         ...
    handles.layers(layer_num).depth_to_top) * 10) / 10;
handles.layers(layer_num+1).thickness =                                 ...
    round((handles.layers(layer_num+1).thickness -                      ...
    (-CurrentPoint(2)-handles.layers(layer_num+1).depth_to_top))*10)/10;
handles.layers(layer_num+1).depth_to_top = round(-CurrentPoint(2)*10)/10;

if strcmp(handles.config.type,'TEM Central Loop')
    yscale = diff(get(handles.axes1,'ylim'));
    t = sin([0:pi/100:pi]).*yscale./70;
    set(handles.config.plot_handle(1),'ydata',t.*4);
    set(handles.config.plot_handle(2),'ydata',t.*1);
elseif strcmp(handles.config.type,'HCP FDEM (HLEM)')
    yscale = diff(get(handles.axes1,'ylim'));
    t = sin([0:pi/100:2.*pi]').*yscale./70;
    set(handles.config.plot_handle(1), 'ydata', t, 'linewidth', 2,      ...
        'color', 'b');
    set(handles.config.plot_handle(2), 'ydata', t, 'linewidth', 2,      ...
        'color', 'b');
    set(handles.config.plot_handle(3), 'ydata',                         ...
        [-yscale +yscale].*1.5./70, 'linewidth', 2, 'color', 'b');
    set(handles.config.plot_handle(4), 'ydata',                         ...
        [-yscale +yscale].*1.5./70, 'linewidth', 2, 'color', 'b');
end

handles = update_model(handles);

update_layer_param(handles);