% clear
addpath 'D:\EM_work\MT3D\bin'
addpath 'D:\EM_work\MT3D\bin\MT1D'
close all

freq = logspace(4,7,21);
% freq = linspace(1e6,1e8,21);
% freq = logspace(2,3,11);
nfreq = length(freq);

% zNode = (-20:1:max(mesh.node2coord(:,3)))';
zNode = (-0:1:100)';
% zNode =zNode /30;
centerZ  = zNode(1:end-1)+0.5*diff(zNode) ; % 层中心Z坐标
% zNode = zNode/10;
id = find(abs(zNode)== min(abs(zNode)));

sig1D = 1e-12*ones(size(zNode,1)-1,1);
sig1D(centerZ>0,:)  = 1/1e3;
% sig1D(centerZ>10,:) = 1/1e2 ;

epsr  = ones(size(zNode,1)-1,1);
epsr(centerZ>0,:) = 3;
% epsr(centerZ>3,:) = 1;

mur = ones(size(zNode,1)-1,1);
% mur(centerZ>0,:) = 1;

%


figure(1)
subplot(1,2,1);
stairs( [sig1D(1,1);sig1D(:,1)],zNode,'k',"LineWidth",2)
title('Conductivity model');
ylabel('Z(km)') ;xlabel('sigma(S/m)');axis tight;box on;
set(gca,'YDir','reverse','Xscale','log');
subplot(1,2,2);
stairs( [epsr(1,1);epsr(:,1)],zNode,'k',"LineWidth",2)
title('Permittivity model');
ylabel('Z(km)') ;xlabel('epsilon');axis tight;box on;
set(gca,'YDir','reverse','Xscale','log');

%


data = table;
col = jet(nfreq);
f2 = figure(2);
l={};
for ifreq =1:nfreq
    [eField,hField]         =  mt1DAnalyticFieldsEpsMu(freq(ifreq),sig1D,epsr,mur,zNode,1,'E', 1);
    % [eField,hField]         =  mt1DAnalyticFields(freq(ifreq),sig1D,zNode,1,'E', 1);
    Ex1 = eField(id).';
    Ey1 = Ex1*0 ;
    Hx1 =  Ex1*0 ;
    Hy1 =  hField(id).' ;
    Hz1 = Ex1*0 ;
    Ex2 = Ey1 ;
    Ey2 = Ex1 ;
    Hx2 = Hy1 ;
    Hy2 = Ex1*0 ;
    Hz2 = Ex1*0 ;

    mtr = get_MTR(Ex1,Ey1,Hx1,Hy1,Hz1,Ex2,Ey2,Hx2,Hy2,Hz2,freq(ifreq));
    Frequency  =  repmat(freq(ifreq),1,1);
    RX   =(1:1)';

    datatemp = table(Frequency, RX,Ex1,Ey1,Hx1,Hy1,Hz1,Ex2,Ey2,Hx2,Hy2,Hz2,...
        mtr.Zxx, mtr.Zxy, mtr.Zyx, mtr.Zyy,mtr.rhoXX, mtr.rhoXY, mtr.rhoYX, mtr.rhoYY,...
        mtr.phaXX, mtr.phaXY, mtr.phaYX, mtr.phaYY,...
        'VariableNames',{'Frequency';'RX';'Ex1';'Ey1';'Hx1';'Hy1';'Hz1';'Ex2';...
        'Ey2';'Hx2';'Hy2';'Hz2';'Zxx';'Zxy';'Zyx';'Zyy';'rhoXX';'rhoXY';'rhoYX';'rhoYY'; ...
        'phaXX';'phaXY';'phaYX';'phaYY'});
    data = cat(1,data,datatemp);

    figure(f2);
    subplot(1,4,1);hold on;

    plot(real(eField),zNode,'-','LineWidth',1.5,'Color',col(ifreq,:));
    ylabel('Z(km)') ;xlabel('real(E)');
    % title('Real part of  electric field')
    set(gca,'YDir','reverse','XDir','reverse');
    axis tight;box on;

    subplot(1,4,2);hold on;
    plot(imag(eField),zNode,'-','LineWidth',1.5,'Color',col(ifreq,:));
    ylabel('Z(km)') ;xlabel('imag(E)');
    % title('Imaginary part of  electric field')
    set(gca,'YDir','reverse','XDir','reverse');
    axis tight;box on;

    subplot(1,4,3);hold on;
    plot(abs(eField),zNode,'-','LineWidth',1.5,'Color',col(ifreq,:));
    ylabel('Z(km)') ;xlabel('$|\mathbf{E}|$','Interpreter','latex');
    % title('Amplitude of electric field')
    set(gca,'YDir','reverse','XDir','reverse');
    axis tight;box on;
    % xlim([0,1])
    l{ifreq} = [num2str(freq(ifreq)/1000),'kHz'] ;
    legend(l,'FontSize',6,"Position", [0.722,0.392,0.15,0.53],Box="off");
    
    % MakeGif('my.Gif',ifreq);


    drawnow ;
    pause(0.3);
end

%

figure;
% yyaxis left;
p=plot(freq,data.rhoXY,'-s','LineWidth',2);
c =get(p,'Color');
xlabel('frequency/Hz','LineWidth',1,'FontWeight','bold');
ylabel('Apparent resistivity(ohm.m)','LineWidth',1,'FontWeight','bold');
% hold on;
% loglog(freq,rhoaXY_3D,'.-','LineWidth',2);
xlabel('frequency(Hz)','FontWeight','bold');
% ylim([500,2000]);
set(gca,'LineWidth',1,'XScale','log');
box on;axis tight;

%

figure;
% yyaxis left;
p=plot(freq,data.phaYX,'-s','LineWidth',2);
c =get(p,'Color');
xlabel('frequency/Hz','LineWidth',1,'FontWeight','bold');
ylabel('Phass(degree)','LineWidth',1,'FontWeight','bold');
% hold on;
% loglog(freq,rhoaXY_3D,'.-','LineWidth',2);
xlabel('frequency(Hz)','FontWeight','bold');
% ylim([500,2000]);
set(gca,'LineWidth',1,'XScale','log');
box on;axis tight;