FigDir = ['../../fig_revision1'];

FlagExp0 = 0;
FlagExp1 = 0;
FlagExp2 = 0;
FlagExp3 = 0;
FlagExp4 = 1;
curdir        = pwd;

Cfile   = mfilename('fullpath');
Cfile   = [Cfile '.m'];
copyfile(Cfile,FigDir);

PaperSize = [12, 5];
PaperSizeME = [12 8];
FONTS      = 20;
fcaxis          = [1.48,5];

velfile = {'/scratch/slim/zfang/Model/BG2D/BG450kmTrue.mat', ...
    '/scratch/slim/zfang/Model/BG2D/Ini_Vel.mat'};

strfile = {'BG_TrueModel', ...
    'BG_InitialModel'};

strtitle = {'True Model', ...
    'Initial Model'};


for i = 1:length(velfile)
    [A n d o] = ReadAllData(velfile{i});
    [z ,x]    = odn2grid(o,d,n);
    figure; set(gca,'fontsize',FONTS);imagesc(x,z,A);
    xlabel('Lateral [m]'); ylabel('Depth [m]');
    set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
    %; title(strtitle{i})
    caxis(fcaxis);fidx=gcf;
    PrintFigure(fidx,PaperSize,[FigDir '/' strfile{i}]);
end

if FlagExp0 > 0;
    Velfile = {'/scratch/slim/zfang/Model/SimpleLayer/SimpleLayerkm.mat', ...
        '/scratch/slim/zfang/Model/SimpleLayer/Ini_Vel.mat', ...
        '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/Layer_WrongSrc_15Hz/m_2.mat', ...
        '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/Layer_TrueSrc2_15Hz/m_2.mat' };
    
    str = {'Ture', 'Ini', 'Wrong', 'Correct'};
    
    for i = 1:length(Velfile);
        [A ns ds os] = ReadAllData(Velfile{i});
        [zz,xx]     = odn2grid(os,ds,ns);
        figure;imagesc(xx,zz,A);caxis([1,3]);set(gca,'fontsize',FONTS);
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
        fidx=gcf;
        PrintFigure(fidx,PaperSize,[FigDir '/Layer_V' str{i} ]);
        
    end
    
    f    = 15;
    t0t = 0.3;
    t0w = 0.32;
    fw    = 14;
    dt    = 0.004;
    nt    = 501;
    
    wt   = getrick(f,t0t,dt,nt);
    ww   = getrick(fw,t0w,dt,nt);
    [tt]  = odn2grid(0,dt,nt);
    figure; plot(tt,[wt(:)],'linewidth',2,'linestyle','-.');hold on;plot(tt(:),ww(:),'linewidth',2,'linestyle','--')
    set(gca,'fontsize',FONTS);legend('Correct Source','Wrong Source');
    xlabel('Time [s]'); set(gca,'fontsize',FONTS);xlim([0,2]);
    fidx=gcf;
    PrintFigure(fidx,PaperSize,[FigDir '/Layer_SrcCmp'  ]);
    
    
    
end



if FlagExp1 > 0;
    
    RstTrueSrcDir = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_TrueSrc_20Hz';
    RstEstSrcDir     = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz';
    RstEstSrcDir2  = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_ModelSrc2';
    
    ves_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_TrueSrc_20Hz/m_11.mat';
    vts_file  = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz/m_11.mat';
    vesP_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_ModelSrc2/m_11.mat';
    est_src  = ReadAllData('/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz/EstSrcWave.mat');
    est_srcP = ReadAllData('/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_ModelSrc2/EstSrcWave.mat');
    
    ves          = ReadAllData(ves_file);
    vts           = ReadAllData(vts_file);
    vesP       = ReadAllData(vesP_file);
    
    freq        = [2:0.5:29.5];
    f0             = 20;
    t0             = 0.5;
    flagErr   = 0;
    flagModel = 1;
    flagGradient = 0;
    
    if flagGradient > 0
        gt_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz/G_t.mat';
        gw_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz/G_w.mat';
        ge_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz/G_e.mat';
        [gt, n, d, o] = ReadAllData(gt_file);
        [gw, n, d, o] = ReadAllData(gw_file);
        [ge, n, d, o] = ReadAllData(ge_file);
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,gt);
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);caxis([-2,2]*1e4);
        %; title(strtitle{i})
        fidx=gcf;
        PrintFigure(fidx,PaperSize,[FigDir '/G_truesrc' ]);
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,gw);
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);caxis([-2,2]*1e4);
        %; title(strtitle{i})
        fidx=gcf;
        PrintFigure(fidx,PaperSize,[FigDir '/G_wrongsrc' ]);
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,ge);
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);caxis([-2,2]*1e4);
        %; title(strtitle{i})
        fidx=gcf;
        PrintFigure(fidx,PaperSize,[FigDir '/G_estsrc' ]);
        
        
    end
    
    
    if flagModel > 0
        
        src_true = fwi_wavelet(freq,t0,f0);
        
        bb  = norm(est_src(1,:)) / norm(est_srcP);
        est_srcP = est_srcP*bb;
        
        
        
        figure;set(gca,'fontsize',20);plot(freq(:), abs([src_true(:)]), '+','linewidth',1,'color','black'); hold on;
        plot(freq(:), abs([est_srcP(:)]), 'x','linewidth',1,'color','black');
        plot(freq(:), abs(vec(est_src(1,:))),'o', 'linewidth',1,'color','black');
        h = legend('True source','Estimated source with WRI-SE-MS',  'Estimated source with WRI-SE-WS','Location','South');
        set(h,'fontsize',10);
        xlabel('Frequency [Hz]','fontsize',20);
        ylabel('Amplitude','fontsize',20);
        str = [FigDir '/Freq_SrcCmp_Amp_all'];
        set(gca,'fontsize',20);xlim([0 32])
        PrintFigure(gcf,[7,5],str);
        
        figure;set(gca,'fontsize',20);plot(freq(:), abs([src_true(:)]), '--+','linewidth',2); hold on;
        plot(freq(:), abs(vec(est_src(1,:))),'o', 'linewidth',2);
        h = legend('True source','Estimated source by VP','Location','SouthEast');
        set(h,'fontsize',10);
        xlabel('Frequency [Hz]','fontsize',20);
        ylabel('Amplitude','fontsize',20);
        str = [FigDir '/Freq_SrcCmp_Amp'];
        set(gca,'fontsize',20);xlim([0 32])
        PrintFigure(gcf,[7,5],str);
        
        aa = angle(vec(est_src(1,:)));
        ii   =find(aa>2);
        aa(ii) = aa(ii) - 2*pi;
       
        freq1         = [2:0.1:31];
        src_true1 = fwi_wavelet(freq1,t0,f0);
        %   plot(freq1(:), angle([src_true1(:)]), '--','linewidth',2,'color','red');hold on
        cc = angle([src_true(:)]);
        ii   =find(cc>2);
        cc(ii) = cc(ii) - 2*pi;
        bb = angle(est_srcP);
        id = find(bb >2);
        bb(id) = bb(id)-2*pi;
        id = find(bb < -4)
        bb(id) = bb(id) + 2*pi;
        id = find(bb >3);
        bb(id) = bb(id)-2*pi;
%        plot(freq(:), cc, '+','linewidth',2);hold on;
%        plot(freq(:), bb, 'x','linewidth',2,'color','black');
%        plot(freq(:), aa,'o', 'linewidth',2,'color','red');
        figure;set(gca,'fontsize',20);
        plot(freq(:), cc, '+','linewidth',1,'color','black');hold on;
        plot(freq(:), bb, 'x','linewidth',1,'color','black');
        plot(freq(:), aa,'o', 'linewidth',1,'color','black');
        h = legend('True source','Estimated source with WRI-SE-MS',  'Estimated source with WRI-SE-WS','Location','NorthEast');
        set(h,'fontsize',10);
        xlabel('Frequency [Hz]','fontsize',20);
        ylabel('Phase','fontsize',20);
        str = [FigDir '/Freq_SrcCmp_Phase_all'];
        set(gca,'fontsize',20);xlim([0 32]);ylim([-4 4])
        PrintFigure(gcf,[7,5],str);
        
        figure;set(gca,'fontsize',20);
        plot(freq(:), cc, '+','linewidth',2);hold on;
        plot(freq(:), aa,'o', 'linewidth',2,'color','red');
        h = legend('True source','Estimated source by VP');
        set(h,'fontsize',10);
        xlabel('Frequency [Hz]','fontsize',20);
        ylabel('Phase','fontsize',20);
        str = [FigDir '/Freq_SrcCmp_Phase'];
        set(gca,'fontsize',20);xlim([0 32]);ylim([-4 4])
        PrintFigure(gcf,[7,5],str);
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,ves);
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
        %; title(strtitle{i})
        caxis(fcaxis);fidx=gcf;
        PrintFigure(fidx,PaperSize,[FigDir '/vel_freq_estsrc' ]);
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,vts); caxis(fcaxis);fidx=gcf;
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
        PrintFigure(fidx,PaperSize,[FigDir '/vel_freq_truesrc' ]);
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,vesP); caxis(fcaxis);fidx=gcf;
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
        PrintFigure(fidx,PaperSize,[FigDir '/vel_freq_estsrcP' ]);
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,abs(vesP-vts)./vts);
        xlabel('Lateral [m]'); ylabel('Depth [m]');fidx=gcf;
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');colormap(flipud(gray))
        PrintFigure(fidx,PaperSize,[FigDir '/P_reldiff_freq_truesrc' ]);
        acaxis = caxis;
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,abs(ves-vts)./vts);
        xlabel('Lateral [m]'); ylabel('Depth [m]');fidx=gcf;
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');colormap(flipud(gray));caxis(acaxis)
        PrintFigure(fidx,PaperSize,[FigDir '/reldiff_freq_truesrc' ]);
        
        
        
    end
    
    if flagErr > 0
        vt = ReadAllData('/scratch/slim/zfang/Model/BG2D/BG450kmTrue.mat');
        vi = ReadAllData('/scratch/slim/zfang/Model/BG2D/Ini_Vel.mat');
        cd(RstTrueSrcDir);
        k = 1;
        for i = 1:11
            cd(['xm' num2str(i)]);
            for j = 1:21
                xx         = load(['x_' num2str(j) '.dat']);
                METrue(k) = norm(1./sqrt(xx)-vt(:));
                k         = k+1;
            end
            cd('../')
        end
        
        cd(RstEstSrcDir);
        k = 1;
        for i = 1:11
            cd(['xm' num2str(i)]);
            for j = 1:21
                xx         = load(['x_' num2str(j) '.dat']);
                MEEst2(k)  = norm(1./sqrt(xx)-vt(:));
                k         = k+1;
            end
            cd('../')
        end
        cd(curdir)
        
        cd(RstEstSrcDir2);
        k = 1;
        for i = 1:11
            cd(['xm' num2str(i)]);
            for j = 1:21
                file_t   = ['./x_' num2str(j) '.dat'];
                if exist(file_t,'file') > 0
                    xx         = load(file_t);
                end
                MEEst1(k)  = norm(1./sqrt(xx(1:prod(n)))-vt(:));
                k         = k+1;
            end
            cd('../')
        end
        cd(curdir)
        
        METrue = METrue / norm(vi-vt);
        MEEst2 = MEEst2 / norm(vi-vt);
        MEEst1 = MEEst1 / norm(vi-vt);
        
        figure;set(gca,'fontsize',20);p=plot([METrue(:) MEEst1(:) MEEst2(:)],'linewidth',3);
        set(p(1),'color','black');set(p(1),'LineStyle','--');
        set(p(2),'color','black');set(p(2),'LineStyle','-');
        set(p(3),'color','black');set(p(3),'LineStyle',':');
        xlabel('Iteration');ylabel('Relative Model Error');legend('Ture Source Wavelet', 'WRI-SE-MS','WRI-SE-WS');set(gca,'fontsize',FONTS);xlim([0,length(METrue)])
        fidx = gcf;PrintFigure(fidx, PaperSizeME, [FigDir '/ModelErrorComp'])
        
    end
    
end


if FlagExp2 > 0;
    
    RstTrueSrcDir = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/TrueSrc';
    RstEstSrcDir     = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEst';
    RstEstSrcDir2   = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_ModelSrc_TimeData2';
    
    ves_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEst_TimeData_50pml/m_11.mat';
    vts_file  = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/TrueSrc/m_11.mat';
    vesP_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_ModelSrc_TimeData2/m_11.mat';
    est_src  = ReadAllData('/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEst_TimeData_50pml/EstSrcWave.mat');
    est_srcP = ReadAllData('/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_ModelSrc_TimeData2/EstSrcWave.mat');
    
    ves          = ReadAllData(ves_file);
    vesP       = ReadAllData(vesP_file);
    %    vts           = ReadAllData(vts_file);
    
    freq        = [2:0.5:29.5];
    f0             = 20;
    t0             = 0.5;
    flagErr   = 0;
    
    src_true = -fwi_wavelet(freq,t0,f0);
    est_src   = est_src(1,1:56);
    
    
    aa                = norm(vec(est_src)) / norm(src_true);
    
    figure;set(gca,'fontsize',20);plot(freq(:), abs([aa*src_true(:)]), '--+','linewidth',2);hold on;
    plot(freq(:), abs(vec(est_src(1,:))),'o', 'linewidth',2);
    h = legend('True source', 'Estimated source by VP' ,'Location','SouthEast');
    set(h,'fontsize',10);
    xlabel('Frequency [Hz]','fontsize',20);
    ylabel('Amplitude','fontsize',20);
    str = [FigDir '/Time_SrcCmp_Amp'];
    set(gca,'fontsize',20);xlim([0 32])
    PrintFigure(gcf,[7,5],str);
    
    
    bb = norm(est_src)/norm(est_srcP);
    est_srcP = est_srcP*bb;
    figure;set(gca,'fontsize',20);plot(freq(:), abs(aa*[src_true(:)]), '+','linewidth',1,'color','black'); hold on;
    plot(freq(:), abs([est_srcP(:)]), 'x','linewidth',1,'color','black');
    plot(freq(:), abs(vec(est_src(1,:))),'o', 'linewidth',1,'color','black');
    h = legend('True source','Estimated source with WRI-SE-MS',  'Estimated source with WRI-SE-WS','Location','South');
    set(h,'fontsize',10);
    xlabel('Frequency [Hz]','fontsize',20);
    ylabel('Amplitude','fontsize',20);
    str = [FigDir '/Time_SrcCmp_Amp_all'];
    set(gca,'fontsize',20);xlim([0 32])
    PrintFigure(gcf,[7,5],str);
    
    a_est_src =  angle(vec(est_src(1,:)));
    a_est_src(1) = a_est_src(1) - 2*pi;
    a_true_src     = angle(src_true);
    a_true_src(49) = a_true_src(49) - 2*pi;
    figure;set(gca,'fontsize',20);
    freq1         = [2:0.1:31];
    src_true1 = fwi_wavelet(freq1,t0,f0);
    %   plot(freq1(:), angle([src_true1(:)]), '--','linewidth',2,'color','red');hold on
    plot(freq(:), [a_true_src(:)], '+','linewidth',2);hold on
    plot(freq(:),a_est_src,'o', 'linewidth',2,'color','red');
    h = legend('True source', 'Estimated source by VP');
    set(h,'fontsize',10);
    xlabel('Frequency [Hz]','fontsize',20);
    ylabel('Phase','fontsize',20);
    str = [FigDir '/Time_SrcCmp_Phase'];
    set(gca,'fontsize',20);xlim([0 32]);ylim([-4 4])
    PrintFigure(gcf,[7,5],str);
    
    aa = angle(vec(est_src(1,:)));
    ii   =find(aa>2);
    aa(ii) = aa(ii) - 2*pi;
    figure;set(gca,'fontsize',20);
    freq1         = [2:0.1:31];
    src_true1 = fwi_wavelet(freq1,t0,f0);
    %   plot(freq1(:), angle([src_true1(:)]), '--','linewidth',2,'color','red');hold on
    cc = angle([src_true(:)]);
    ii   =find(cc>2);
    cc(ii) = cc(ii) - 2*pi;
    bb = angle(est_srcP);
    id = find(bb >2);
    bb(id) = bb(id)-2*pi;
    id = find(bb < -4)
    bb(id) = bb(id) + 2*pi;
    id = find(bb >3);
    bb(id) = bb(id)-2*pi;
    plot(freq(:), cc, '+','linewidth',1,'color','black');hold on;
    plot(freq(:), bb, 'x','linewidth',1,'color','black');
    plot(freq(:), aa,'o', 'linewidth',1,'color','black');
    h = legend('True source','Estimated source with WRI-SE-MS',  'Estimated source with WRI-SE-WS','Location','NorthEast');
    set(h,'fontsize',10);
    xlabel('Frequency [Hz]','fontsize',20);
    ylabel('Phase','fontsize',20);
    str = [FigDir '/Time_SrcCmp_Phase_all'];
    set(gca,'fontsize',20);xlim([0 32]);ylim([-4 4])
    PrintFigure(gcf,[7,5],str);
    
    
    
    figure; set(gca,'fontsize',FONTS);imagesc(x,z,ves);
    xlabel('Lateral [m]'); ylabel('Depth [m]');
    set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
    %; title(strtitle{i})
    caxis(fcaxis);fidx=gcf;
    PrintFigure(fidx,PaperSize,[FigDir '/vel_time_estsrc' ]);
    
    figure; set(gca,'fontsize',FONTS);imagesc(x,z,vesP);
    xlabel('Lateral [m]'); ylabel('Depth [m]');
    set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
    %; title(strtitle{i})
    caxis(fcaxis);fidx=gcf;
    PrintFigure(fidx,PaperSize,[FigDir '/vel_time_estsrc_P' ]);
    
end

if FlagExp3 > 0;
    
    RstTrueSrcDir = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/TrueSrc';
    RstEstSrcDir     = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEst';
    RstEstSrcDir2   = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_ModelSrc_TimeData2';
    
    ves_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEst_TimeDataNoise1e0New4/Noise1e1/m_11.mat';
    vts_file  = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/TrueSrc/m_11.mat';
    vesP_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_ModelSrc_Noise2/m_11.mat';
    est_src  = ReadAllData('/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEst_TimeDataNoise1e0New4/Noise1e1/EstSrcWave.mat');
    estP_src = ReadAllData('/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_ModelSrc_Noise2/EstSrcWave.mat');
    
    ves          = ReadAllData(ves_file);
    vesP       = ReadAllData(vesP_file);
    %    vts           = ReadAllData(vts_file);
    
    freq        = [2:0.5:29.5];
    f0             = 20;
    t0             = 0.5;
    flagErr   = 0;
    
    est_src   = est_src(1,1:56);
    
    src_true = -fwi_wavelet(freq,t0,f0);
    
    
    aa                = norm(vec(est_src)) / norm(src_true);
    bb                = norm(est_src)/norm(estP_src);
    est_srcP    = bb * estP_src;
    
    figure;set(gca,'fontsize',20);plot(freq(:), abs([aa*src_true(:)]), '--+','linewidth',2);hold on
    plot(freq(:), abs(vec(est_src(1,:))),'o', 'linewidth',2);
    h = legend('True source', 'Estimated source','Location','SouthEast');
    set(h,'fontsize',10);
    xlabel('Frequency [Hz]','fontsize',20);
    ylabel('Amplitude','fontsize',20);
    str = [FigDir '/Time_SrcCmp_Amp_Noise'];
    set(gca,'fontsize',20);xlim([0 32]);ylim([0 800])
    PrintFigure(gcf,[7,5],str);
    
    figure;set(gca,'fontsize',20);plot(freq(:), abs(aa*[src_true(:)]), '+','linewidth',1,'color','black'); hold on;
    plot(freq(:), abs([est_srcP(:)]), 'x','linewidth',1,'color','black');
    plot(freq(:), abs(vec(est_src(1,:))),'o', 'linewidth',1,'color','black');
    h = legend('True source','Estimated source with WRI-SE-MS',  'Estimated source with WRI-SE-WS','Location','South');
    set(h,'fontsize',10);
    xlabel('Frequency [Hz]','fontsize',20);
    ylabel('Amplitude','fontsize',20);
    str = [FigDir '/Time_SrcCmp_Amp_Noise_all'];
    set(gca,'fontsize',20);xlim([0 32])
    PrintFigure(gcf,[7,5],str);
    
    a_est_src =  angle(vec(est_src(1,:)));
    aa               = find(a_est_src > 2);
    a_est_src(aa) = a_est_src(aa) - 2*pi;
    a_true_src       = angle(src_true);
    aa                        = find(a_true_src > 2);
    a_true_src(aa) = a_true_src(aa) - 2*pi;
    
    figure;set(gca,'fontsize',20);
    freq1         = [2:0.5:29.5];
    src_true1 = fwi_wavelet(freq1,t0,f0);
    %   plot(freq1(:), angle([src_true1(:)]), '--','linewidth',2,'color','red');hold on
    plot(freq(:), [a_true_src(:)], '+','linewidth',2);;hold on
    plot(freq(:),a_est_src,'o', 'linewidth',2,'color','red')
    h = legend('True source', 'Estimated source' );
    set(h,'fontsize',10);
    xlabel('Frequency [Hz]','fontsize',20);
    ylabel('Phase','fontsize',20);
    str = [FigDir '/Time_SrcCmp_Phase_Noise'];
    set(gca,'fontsize',20);xlim([0 32]);ylim([-4 4])
    PrintFigure(gcf,[7,5],str);
    
    figure;set(gca,'fontsize',20);
    freq1         = [2:0.5:29.5];
    src_true1 = fwi_wavelet(freq1,t0,f0);
    %   plot(freq1(:), angle([src_true1(:)]), '--','linewidth',2,'color','red');hold on
    cc = angle([src_true(:)]);
    ii   =find(cc>2);
    cc(ii) = cc(ii) - 2*pi;
    bb = angle(est_srcP);
    id = find(bb >2);
    bb(id) = bb(id)-2*pi;
    id = find(bb < -4)
    bb(id) = bb(id) + 2*pi;
    id = find(bb >3);
    bb(id) = bb(id)-2*pi;
    aa = angle(vec(est_src(1,:)));
    ii   =find(aa>2);
    aa(ii) = aa(ii) - 2*pi;
    plot(freq(:), cc, '+','linewidth',1,'color','black');hold on;
    plot(freq(:), bb, 'x','linewidth',1,'color','black');
    plot(freq(:), aa,'o', 'linewidth',1,'color','black');
    h = legend('True source','Estimated source with WRI-SE-MS',  'Estimated source with WRI-SE-WS','Location','NorthEast');
    set(h,'fontsize',10);
    xlabel('Frequency [Hz]','fontsize',20);
    ylabel('Phase','fontsize',20);
    str = [FigDir '/Time_SrcCmp_Phase_Noise_all'];
    set(gca,'fontsize',20);xlim([0 32]);ylim([-4 4])
    PrintFigure(gcf,[7,5],str);
    
    figure; set(gca,'fontsize',FONTS);imagesc(x,z,ves);
    xlabel('Lateral [m]'); ylabel('Depth [m]');
    set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
    %; title(strtitle{i})
    caxis(fcaxis);fidx=gcf;
    PrintFigure(fidx,PaperSize,[FigDir '/vel_time_estsrc_noise' ]);
    
    figure; set(gca,'fontsize',FONTS);imagesc(x,z,vesP);
    xlabel('Lateral [m]'); ylabel('Depth [m]');
    set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
    %; title(strtitle{i})
    caxis(fcaxis);fidx=gcf;
    PrintFigure(fidx,PaperSize,[FigDir '/vel_time_estsrc_P_noise' ]);
    
end

if FlagExp4 > 0;
    
    RstTrueSrcDir = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/TrueSrc';
    RstEstSrcDir     = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEst';
    RstEstSrcDir2   = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_ModelSrc_TimeData2';
    
    ves_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_LinearIni44/m_11.mat';
    vts_file  = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/TrueSrc/m_11.mat';
    vesP_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_FWI_TrueSrc_20Hz_LinearModel36/m_11.mat';
    vini_file   = '/scratch/slim/zfang/Model/BG2D/BGIniLinear17.mat';
    vt_file   = '/scratch/slim/zfang/Model/BG2D/BG450kmTrue.mat';
    est_src  = ReadAllData('/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_LinearIni44/EstSrcWave.mat');
    estP_src = ReadAllData('/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_FWI_TrueSrc_20Hz_LinearModel36/EstSrcWave.mat');
    
    xps         = [101, 251, 401];
    
    ves          = ReadAllData(ves_file);
    vesP       = ReadAllData(vesP_file);
    vini         = ReadAllData(vini_file);
    vt         = ReadAllData(vt_file);
    %    vts           = ReadAllData(vts_file);
    
    
    
    
    
    freq        = [2:0.5:29.5];
    f0             = 20;
    t0             = 0.5;
    flagErr   = 0;
    flagsrc   = 1;
    flagmodel = 0;
    flagdiff   = 0;
    flagDataf  = 1;
    flagData   = 1;
    if flagsrc > 0
        f_slice    = [11:56];
        est_src   = est_src(1,f_slice);
        
        src_true = fwi_wavelet(freq,t0,f0);
        src_true = src_true(f_slice);
        estP_src = estP_src(f_slice);
        
        
        aa                = norm(vec(est_src)) / norm(src_true);
        bb                = norm(est_src)/norm(estP_src);
        est_srcP    = bb * estP_src;
        
        figure;set(gca,'fontsize',20);plot(freq(f_slice), abs([aa*src_true(:)]), '--+','linewidth',2);hold on
        plot(freq(f_slice), abs(vec(est_src(1,:))),'o', 'linewidth',2);
        h = legend('True source', 'Estimated source','Location','SouthEast');
        set(h,'fontsize',10);
        xlabel('Frequency [Hz]','fontsize',20);
        ylabel('Amplitude','fontsize',20);
        str = [FigDir '/Time_SrcCmp_WRI_Amp'];
        set(gca,'fontsize',20);xlim([6 31]);ylim([0 800])
        PrintFigure(gcf,[7,5],str);
        
        figure;set(gca,'fontsize',20);plot(freq(f_slice), abs(aa*[src_true(:)]), '+','linewidth',1,'color','black'); hold on;
        plot(freq(f_slice), abs([est_srcP(:)]), 'x','linewidth',1,'color','black');
        plot(freq(f_slice), abs(vec(est_src(1,:))),'o', 'linewidth',1,'color','black');
        h = legend('True source','Estimated source of FWI',  'Estimated source of WRI','Location','SouthEast');
        set(h,'fontsize',10);
        xlabel('Frequency [Hz]','fontsize',20);
        ylabel('Amplitude','fontsize',20);
        xlim([6,31])
        ylim([40 160])
        str = [FigDir '/Time_SrcCmp_FWIWRI_AMP'];
        set(gca,'fontsize',20);xlim([6 31])
        PrintFigure(gcf,[7,5],str);
        
        a_est_src =  angle(vec(est_src(1,:)));
        aa               = find(a_est_src > 2);
        a_est_src(aa) = a_est_src(aa) - 2*pi;
        a_true_src       = angle(src_true);
        aa                        = find(a_true_src > 2);
        a_true_src(aa) = a_true_src(aa) - 2*pi;
        
        figure;set(gca,'fontsize',20);
        freq1         = [2:0.5:29.5];
        src_true1 = fwi_wavelet(freq1,t0,f0);
        %   plot(freq1(:), angle([src_true1(:)]), '--','linewidth',2,'color','red');hold on
        plot(freq(f_slice), [a_true_src(:)], '+','linewidth',2);;hold on
        plot(freq(f_slice),a_est_src,'o', 'linewidth',2,'color','red')
        h = legend('True source', 'Estimated source' );
        set(h,'fontsize',10);
        xlabel('Frequency [Hz]','fontsize',20);
        ylabel('Phase','fontsize',20);
        xlim([6,31])
        str = [FigDir '/Time_SrcCmp_WRI_Est'];
        set(gca,'fontsize',20);xlim([6 31]);ylim([-4 4])
        PrintFigure(gcf,[7,5],str);
        
        figure;set(gca,'fontsize',20);
        freq1         = [2:0.5:29.5];
        src_true1 = fwi_wavelet(freq1,t0,f0);
        src_true1 = src_true1(f_slice);
        %   plot(freq1(:), angle([src_true1(:)]), '--','linewidth',2,'color','red');hold on
        cc = angle([src_true(:)]);
        ii   =find(cc>2);
        cc(ii) = cc(ii) - 2*pi;
        bb = angle(est_srcP);
        id = find(bb >2);
        bb(id) = bb(id)-2*pi;
        id = find(bb < -4)
        bb(id) = bb(id) + 2*pi;
        id = find(bb >3);
        bb(id) = bb(id)-2*pi;
        aa = angle(vec(est_src(1,:)));
        ii   =find(aa>2);
        aa(ii) = aa(ii) - 2*pi;
        plot(freq(f_slice), cc, '+','linewidth',1,'color','black');hold on;
        plot(freq(f_slice), bb, 'x','linewidth',1,'color','black');
        plot(freq(f_slice), aa,'o', 'linewidth',1,'color','black');
        h = legend('True source','Estimated source with FWI',  'Estimated source with WRI','Location','NorthEast');
        set(h,'fontsize',10);
        xlabel('Frequency [Hz]','fontsize',20);
        ylabel('Phase','fontsize',20);
        xlim([6,31])
        str = [FigDir '/Time_SrcCmp_FWIWRI_Phase'];
        set(gca,'fontsize',20);xlim([6 31]);ylim([-4 4])
        PrintFigure(gcf,[7,5],str);
    end
    
    if flagmodel > 0
        for i = 1:length(xps)
            j = xps(i);
            figure; plot(z,[vt(:,j)],'-k','linewidth', 2); hold on
            plot(z,[ves(:,j)],'--k','linewidth', 2); 
            plot(z,[vesP(:,j)],':k','linewidth', 2); 
            xlabel('Velocity [km/s]','fontsize',14)
            ylabel('Depth [m]','fontsize',14);set(gca,'fontsize',14);
            h=legend('True model', 'WRI result', 'FWI result');
            set(h,'fontsize',14);
            xlim([0,2000])
            view(90,90)
            fidx=gcf;
            
            PrintFigure(fidx,[5,8],[FigDir '/cross_' num2str((j-1)*d(1)) ]);
            
            
            
        end
        
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,ves);
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
        %; title(strtitle{i})
        caxis(fcaxis);fidx=gcf;
        PrintFigure(fidx,PaperSize,[FigDir '/vel_WRI' ]);
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,vesP);
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
        %; title(strtitle{i})
        caxis(fcaxis);fidx=gcf;
        PrintFigure(fidx,PaperSize,[FigDir '/vel_FWI' ]);
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,vini);
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
        %; title(strtitle{i})
        caxis(fcaxis);fidx=gcf;
        PrintFigure(fidx,PaperSize,[FigDir '/vel_badini' ]);
        
    end
    
    if flagdiff > 0
        diff1_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_LinearIni38/Diff1.mat';
        diff2_file = '/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_FWI_TrueSrc_20Hz_LinearModel36/Diff1.mat';
        vd1 = ReadAllData(diff1_file);
        vd2 = ReadAllData(diff2_file);
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,vd1);
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
        %; title(strtitle{i})
        caxis([-.05,.05]);fidx=gcf;
        PrintFigure(fidx,PaperSize,[FigDir '/diff_WRI' ]);
        
        figure; set(gca,'fontsize',FONTS);imagesc(x,z,vd2);
        xlabel('Lateral [m]'); ylabel('Depth [m]');
        set(gca,'fontsize',FONTS);hc=colorbar; set(hc,'fontsize',FONTS);title(hc,'km/s');
        %; title(strtitle{i})
        caxis([-.05,.05]*2);fidx=gcf;
        PrintFigure(fidx,PaperSize,[FigDir '/diff_FWI' ]);
    
    end 
    
     if flagData > 0
      file1 = ['/scratch/slim/zfang/Data/BG2D/UQ2017/BGDataForPaper2017badIni.mat'];
      file2 = ['/scratch/slim/zfang/Data/BG2D/DataForSrcEstPaper/WRIsrcFreq_Cent20.mat'];

      D1 = ReadAllData(file1);
      D2 = ReadAllData(file2);

      D1 = squeeze(D1);
      D2 = squeeze(D2);
      D1 = vec(D1(:,45,1));
      D2 = vec(D2(:,45,11));

      figure;plot(x, real([D2]), '--k','linewidth',2); hold on
      plot(x, real([D1]),':k', 'linewidth',2);
      h = legend('Observed data','Predicted data')
      set(h,'fontsize',14);
      xlabel('Receiver [m]','fontsize',14);
      ylabel('Real(\bf{d})','fontsize',14);
      xlim([0 4500])
      set(gca,'fontsize',14);
      PrintFigure(gcf,[6,5],[FigDir '/DataCmpIniReal']);

      figure;plot(x, imag([D2]), '--k','linewidth',2);hold on
      plot(x, imag([D1]),':k', 'linewidth',2);
      h = legend('Observed data','Predicted data')
      set(h,'fontsize',14);
      xlabel('Receiver [m]','fontsize',14);
      ylabel('Imag(\bf{d})','fontsize',14);
      xlim([0 4500])
      set(gca,'fontsize',14);
      PrintFigure(gcf,[6,5],[FigDir '/DataCmpIniImag']);


    end

    if flagDataf > 0
      file1 = ['/scratch/slim/zfang/Data/BG2D/DataForSrcEstPaper/WRIsrcFreq_Cent20.mat'];
      file2 = ['/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_EstSrc_20Hz_LinearIni44/FinalData.mat'];
      file3 = ['/scratch/slim/zfang/Result/WRI_SRCEST/PaperForJournal/WithSrcEstFreq_FWI_TrueSrc_20Hz_LinearModel36/FinalData.mat'];

      Dt = ReadAllData(file1);
      Dw = ReadAllData(file2);
      Df = ReadAllData(file3);

      Dt = squeeze(Dt);
      Dw = squeeze(Dw);
      Df = squeeze(Df);
      Dt = vec(Dt(:,45,11));
      Dw = vec(Dw(:,45,1));
      Df = vec(Df(:,45,1));

      figure;plot(x, real([Dt]), '--k', 'linewidth',2);hold on
      plot(x, real([Dw]),':k', 'linewidth',2);
      h = legend('Observed data','Predicted data')
      set(h,'fontsize',14);
      xlabel('Receiver [m]','fontsize',14);
      ylabel('Real(\bf{d})','fontsize',14);
      xlim([0 4500])
      set(gca,'fontsize',14);
      PrintFigure(gcf,[6,5],[FigDir '/DataCmpWRIReal']);

      figure;plot(x, imag([Dt]), '--k', 'linewidth',2);hold on
      plot(x, imag([Dw]),':k', 'linewidth',2);
      h = legend('Observed data','Predicted data')
      set(h,'fontsize',14);
      xlabel('Receiver [m]','fontsize',14);
      ylabel('Imag(\bf{d})','fontsize',14);
      xlim([0 4500])
      set(gca,'fontsize',14);
      PrintFigure(gcf,[6,5],[FigDir '/DataCmpWRIImag']);

      figure;plot(x, real([Dt]), '--k', 'linewidth',2);hold on
      plot(x, real([Df]),':k', 'linewidth',2);
      h = legend('Observed data','Predicted data')
      set(h,'fontsize',14);
      xlabel('Receiver [m]','fontsize',14);
      ylabel('Real(\bf{d})','fontsize',14);
      xlim([0 4500])
      set(gca,'fontsize',14);
      PrintFigure(gcf,[6,5],[FigDir '/DataCmpFWIReal']);

      figure;plot(x, imag([Dt]), '--k', 'linewidth',2);hold on
      plot(x, imag([Df]),':k', 'linewidth',2);
      h = legend('Observed data','Predicted data')
      set(h,'fontsize',14);
      xlabel('Receiver [m]','fontsize',14);
      ylabel('Imag(\bf{d})','fontsize',14);
      xlim([0 4500])
      set(gca,'fontsize',14);
      PrintFigure(gcf,[6,5],[FigDir '/DataCmpFWIImag']);


    end
    
end



