RFDFWI.m

% RFDFWI - A FWI code for FWI of GPR data By Mrinal
%[This code is a MATLAB script for Full-Waveform Inversion (FWI) of
%   Ground Penetrating Radar (GPR) data. FWI is a geophysical imaging
%   technique used to reconstruct subsurface properties
%   (such as electromagnetic permittivity and conductivity)
%      by iteratively comparing modeled and observed data. ]
%RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRef
% For references, see Lavoué et al.(2014); 
%        Layek & Sengupta (2019,2021,&2023)
% REF: Köhn, D., De Nil, D. and Rabbel, W. (2017) Tutorial: 
%      Introduction to frequency domain modelling and FWI of georadar data 
%      with GERMAINE, DOI:10.13140/RG.2.2.29354.03523
%                      ____________________________
%Layek, M. K.,Sengupta, P., Multi parameter imaging by finite difference 
% frequencydomain full waveform inversion of GPR data: A guide for  
% sedimentary architecture modeling, 14 June 2023, PREPRINT (Version 1) 
% available @ Research Square[https://doi.org/10.21203/rs.3.rs 3056830/v1]
% % --------------------------------------------------------
%   Copyright @ Mrinal Kanti Layek
%   Senior Research Fellow (Geophysics),
%   Department of Geology and Geophysics,
%   Indian Institute of Technology Kharagpur,
%   Kharagpur - 721302, INDIA
%   contact: layek.mk@gmail.com
%   https://www.researchgate.net/profile/Mrinal_Layek
%-----------------------------------------------------------------
%Check pdf file - Notes.PDF for details
% -----------------------------------------------------------------------
% % Codes required:
%  grad_obj_hessMKLnew.m, inp_GPRmodel1.m; inp_GPRmodel.m; readeps1.m;
%  readssig1.m; extend_eps1JAN.m;forward_shotNEW; PML_9p_CFS_stag_para.m;
%  imp_A_TE_9p_cfs_PML_stag1_para or imp_A_TE_9p_cfs_PML_stag2_para;
%  RHS_TE1.m;
%  extract_model_ext, copy_grad_frame_TE; norm_matrix_TE.m ;
%  Tikhonov_grad_TE.m; Tikhonov_cost_TE;steep_descent.m;  wolfe_TENEW.m
%
%Folder : \outputmyINV : for all RFDFWI output figs and .mat files
%|||||||||||||||||||||||||||||||||||||||||||||||||||||||
% Tested ok @ LENEVO Think Station w/ configuartions: Windows 10 Pro,
% AMD Ryzen Threadripper PRO 3955WX 16-Cores 3.90 GHz, RAM 32 GB, 64bit,

close all
clear
% clc;
tic;

% define input model parameters
model1 = inp_GPRmodel1;
model = inp_GPRmodel;


% Read and plot input models

figure
set(gcf,'position',[0 0, model1.screenx, model1.screeny]);
set(gcf,'PaperPositionMode','Auto');
subplot (2,1,1)
eps_model1 = readeps1();% read eps model
eps_model=eps_model1';
%
subplot (2,1,2)
sig_model1 = readsig1(); % read sig model
sig_model=sig_model1';
% -=-=-=---------------
eps_model=imgaussfilt (eps_model,6);%generate initial model from true model
sig_model=imgaussfilt (sig_model,6);% generate initial model

exportgraphics(gcf,'outputmyINV\initial_model.tiff','Resolution',600)

% Apply PML grids to the mesh
if (model1.npml>0)
    [eps_model,model1] = extend_eps1JAN(model1,eps_model);
    [sig_model,model] = extend_eps1JAN(model,sig_model);
end

figure;
subplot (2,1,1)
imagesc(eps_model./model1.eps0);

subplot (2,1,2)
imagesc(model.x,model.y,sig_model);
title ("check this")


exportgraphics(gcf,sprintf('outputmyINV/figs/#0initial_model.tiff'),'Resolution',600)
% saveas(gcf,sprintf('outputmyINV/figs/#0Initial_claylensModelsSurface.tiff'));

% model1.ACQMY=2;
% (((((((((((((((()))))))))))))))
% read source and receiver positions from file
acq = read_acq_jan(model1);
%    shift acquisition geometry if model.npml > 0
if (model1.npml>0)
    acq = shift_acq(model1,acq);

end

fprintf('   model1.nx= %d   model1.ny= %d nsrc=%d  ntr=%d  \n', model1.nx, model1.ny,acq.nsrc, acq.ntr);
%----------------- plot models -----------------------------------
fh = figure();
fh.WindowState = 'maximized';
subplot (1,2,1)
imagesc(eps_model./model1.eps0);
axis ij;
axis equal;
axis tight;
xlabel('Distance [m]');
ylabel('Depth [m]');
load 'seismic.map'
colormap (flipud(seismic));
cbar_handle=colorbar('EastOutside');
hold on
scatter(acq.nxshot,acq.nyshot,'b+','LineWidth',1.5);
hold on
scatter(acq.nxrec,acq.nyrec,'wo','LineWidth',1.5);

subplot (1,2,2)
imagesc(sig_model*1000);
axis ij;
axis equal;
axis tight;
xlabel('Distance [m]');
ylabel('Depth [m]');
load 'seismic.map'
colormap (flipud(seismic));
cbar_handle=colorbar('EastOutside');
hold on
scatter(acq.nxshot,acq.nyshot,'b+','LineWidth',1.5);
hold on
scatter(acq.nxrec,acq.nyrec,'wo','LineWidth',1.5);

clear fh

fprintf('   model1.nx= %d   model1.ny= %d nsrc=%d  ntr=%d  \n', model1.nx, model1.ny,acq.nsrc, acq.ntr);
%************************** INITIALISATION **************************

nshots=acq.nsrc;% total no. of shots
ntr=acq.ntr;% no. of traces
itermax=1500; %number_of_FDFDFWI_iterations

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Observed^^^^^^^^^^^^^^^^^^^^^^^^^^
% load observed/field data (inverse crime method)----------
precobs1= load('obs\obs_data.mat');
precobs2=full(precobs1.precrnew);
precobs=precobs2(:,:,:);

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

%nnnnnnnnnneeeeeeeeeeewwwwwwwwwwwwwwwwwwww dddddddeeeeeeeffffffffffffff
% %simultaneous method (Meles et al. 2010);
% [see Lavoue et al. (2014)/ Layek& Sengupta (2023) for details]
f=[50e6 60e6 70e6 80e6 90e6 100e6 125e6 150e6 175e6 200e6];
w=2*pi*f;
omega2=w.^2;
nw=length(w);
fprintf('model1.nx= %d; model1.ny= %d; nw= %d, <itermax=%d> \n', model1.nx,model1.ny,nw, itermax);
% --------------------------------------------------------------------------
%                  (A)  GRADIENT CALCULATION
% --------------------------------------------------------------------------

parpool('local',16)

for iter=1:itermax
    % -----CALCULATE wavefield, adjoint wavefiled, gradients
    % and hessian all together - best way possible------
     % -------------------------------o-o-o----------------------------------------------------
    [te,te_adj, L2sumgrad, grad_sigma, grad_epsilon, hess_sigma, hess_epsilon ] = grad_obj_hessMKLnew( precobs,eps_model, sig_model,acq,nshots,ntr,nw,w, model1, iter,0);
    % -------------------------------o-o-o----------------------------------------------------

    tepp=te(:,:,nshots,nw);
    %

    % % Plot te and te_adj

    te_adjpp=te_adj(:,:,nshots,nw);

    fh1 = figure();
    fh1.WindowState = 'maximized';
    % figure;
    subplot (1,2,1)
    figer1 = plot_wavefld( model1,tepp);%plot_wvfld( model1,tepp,1) for te
    title(['te iter= '  num2str(iter) 'th iteration']);
    subplot (1,2,2)
    figer2 = plot_wavefld( model1,te_adjpp);%plot_wvfld( model1,tepp,1) for te_adj
    title(['te_adj iter= '  num2str(iter) 'th iteration']);
    % saveas(gcf,sprintf('outputmyINV/figs/01Wavefields_it=_%d_iter=_%d.tiff',model1.nf,iter));
    exportgraphics(gcf,sprintf('outputmyINV/figs/01Wavefields_it=_%d_iter=_%d.tiff',model1.nf,iter),'Resolution',600)
    % clear figure1 figure2
    clear fh1 tepp tepa figer1 figer2
    % oooooooooooooooooooooooooo (3) Assemble Gradient oooooooooooooooooooooo


    fprintf('L2sumgrad= %e  .\n', L2sumgrad);

    %  oooooooooooooooooooooooooo (5)Calculate L2 norm of gradients  oooooooooooooooooooooooooo
    scale_grad_sig = norm_matrix_TE (grad_sigma, model1); %used in Wolfels -not used in my case
    scale_grad_eps = norm_matrix_TE (grad_epsilon, model1);

    disp (' -------------- Gradient plotting ----------------------');
    fh1 = figure();
    fh1.WindowState = 'maximized';
    subplot(1,2,1)
    imagesc(model1.x,model1.y, grad_sigma);
    % colormap (flipud(gray));
    load 'seismic.map'
    colormap (flipud(seismic));
    % clim([-model1.caxis1 model1.caxis1]);
    cbar_handle=colorbar('EastOutside');
    set(gca,'FontSize',12);
    grid on;
    axis ij;
    axis equal;
    axis tight;
    xlabel('Distance [m]');
    ylabel('Depth [m]');
    title(['grad sigma iter= ' num2str(iter) 'th iteration']);

    subplot(1,2,2)
    imagesc(model1.x,model1.y, grad_epsilon);
    % colormap (flipud(gray));
    load 'seismic.map'
    colormap (flipud(seismic));
    % clim([-model1.caxis1 model1.caxis1]);
    cbar_handle=colorbar('EastOutside');
    set(gca,'FontSize',12);
    grid on;
    axis ij;
    axis equal;
    axis tight;
    xlabel('Distance [m]');
    ylabel('Depth [m]');
    title(['grad epsilon iter= ' num2str(iter) 'th iteration']);
    exportgraphics(gcf,sprintf('outputmyINV/figs/02grad_iteration_%d_iter=_%d.tiff',nw,iter), 'Resolution',600 )
    % saveas(gcf,sprintf('outputmyINV/figs/02grad_iteration_%d_iter=_%d.tiff',nw,iter));
    clear fh1
    disp(' ----------------------- SAVE mat FILES ----------------------------Save');
    save('outputmyINV\models\eps_model_old_.mat','eps_model');
    save('outputmyINV\models\sig_model_old_.mat','sig_model');
    save('outputmyINV\wavefield\te_old.mat','te');
    save('outputmyINV\wavefield\te_adj_old.mat','te_adj');
    %     save('outputmyINV\wavefield\res_old.mat','res');
    save('outputmyINV\wavefield\L2sumgrad_old.mat','L2sumgrad');
    save('outputmyINV\grad\grad_sigma_old.mat','grad_sigma');
    save('outputmyINV\grad\grad_epsilon_old.mat','grad_epsilon');
    save('outputmyINV\wavefield\scale_grad_sig_old.mat','scale_grad_sig');
    save('outputmyINV\wavefield\scale_grad_eps_old.mat','scale_grad_sig');

    % --------------------------------------------------------------------------
    %                  (B)  Tikhonov Regularisation
    % --------------------------------------------------------------------------

    % oooooooooooooooooooooooooo (1) Cost function Regularisation ooooooooooooo
    L2_tikhonov = Tikhonov_cost_TE(L2sumgrad,sig_model,eps_model, model1);
    L2t(1)=L2_tikhonov;
    % oooooooooooooooooooooooooo (2) Gradient Regularisation oooooooooooooooooo
    [ grad_sigma_tikh, grad_epsilon_tikh] = Tikhonov_grad_TE(model1,sig_model,eps_model,grad_sigma,grad_epsilon);

    % ----------------------- SAVE mat FILES ----------------------------Save

    save('outputmyINV\grad_tikh\grad_sigma_tikh_old.mat','grad_sigma_tikh');
    save('outputmyINV\grad_tikh\grad_epsilon_tikh_old.mat','grad_epsilon_tikh');


    % ---------------------- GRAD_TIKH plotting ----------------------
    figure
    set(gcf,'position',[0 0, model1.screenx, model1.screeny]);
    set(gcf,'PaperPositionMode','Auto');
    subplot(1,2,1)
    imagesc(model1.x,model1.y, grad_sigma_tikh);
    % colormap (flipud(gray));
    load 'seismic.map'
    colormap (flipud(seismic));
    cbar_handle=colorbar('EastOutside');
    set(gca,'FontSize',20);
    grid on;
    axis ij;
    axis equal;
    axis tight;
    xlabel('Distance [m]');
    ylabel('Depth [m]');
    title(['grad sigma TIKH iter= ' num2str(iter) 'th iteration']);

    subplot(1,2,2)
    imagesc(model1.x,model1.y, grad_sigma_tikh);
    % colormap (flipud(gray));
    load 'seismic.map'
    colormap (flipud(seismic));
    cbar_handle=colorbar('EastOutside');
    set(gca,'FontSize',20);
    grid on;
    axis ij;
    axis equal;
    axis tight;
    xlabel('Distance [m]');
    ylabel('Depth [m]');
    title(['grad epsilon TIKH iter= '  num2str(iter) 'th iteration']);

    %     toc;
    % --------------------------------------------------------------------------
    %                  (C)  Hessian Calculation -
    %          not required calculated in the first step
    % --------------------------------------------------------------------------
    disp('|||||||||||||||| (C)  Hessian Calculation  ||||||||||||||||||||||||');
    %
    %     [te_hess,hess_sigma,hess_epsilon] = pseudo_hess_shin_TEmNEW(model1,sig_model,eps_model,acq,nshots,w,nw,iter);



    save('outputmyINV\Hgrad\hess_sigma_old.mat','hess_sigma');
    save('outputmyINV\Hgrad\hess_epsilon_old.mat','hess_epsilon');

    % oooooooooooooooooooooooooo (1) HEssian plotting ooooooooooooooooooooo
    figure
    set(gcf,'position',[0 0, model1.screenx, model1.screeny]);
    set(gcf,'PaperPositionMode','Auto');
    subplot(1,2,1)
    imagesc(model1.x,model1.y, hess_sigma);
    colormap (flipud(gray));
    load 'seismic.map'
    colormap (flipud(seismic));
    % cbar_handle=colorbar('EastOutside');
    set(gca,'FontSize',12);
    grid on;
    axis ij;
    axis equal;
    axis tight;
    xlabel('Distance [m]');
    ylabel('Depth [m]');
    title([' hess sigma iter= ' num2str(iter) 'th iteration']);

    subplot(1,2,2)
    imagesc(model1.x,model1.y, hess_epsilon);
    colormap (flipud(gray));
    load 'seismic.map'
    colormap (flipud(seismic));
    % cbar_handle=colorbar('EastOutside');
    set(gca,'FontSize',12);
    grid on;
    axis ij;
    axis equal;
    axis tight;
    xlabel('Distance [m]');
    ylabel('Depth [m]');
    title([' hess epsilon iter= '  num2str(iter) 'th iteration']);
    exportgraphics(gcf,sprintf('outputmyINV/figs/03HESS_iteration_%d_%d.tiff',nw,iter), 'Resolution',600 )
    % saveas(gcf,sprintf('outputmyINV/figs/03HESS_iteration_%d_%d.tiff',nw,iter));
    %
    %

    clear teppHess

    disp('oooooooo (Time required upto Hessian oooooooooooooooooo \n');
    toc;
    disp('oooooooooooooooooooooo oooooooooooooo oooooooooooooooooo \n');

    SCALE_GRAD_sig_PRECOND =norm_matrix_TE(grad_sigma,model1);% not used
    SCALE_GRAD_eps_PRECOND = norm_matrix_TE(grad_epsilon,model1); % not used

    % *** calculate descent directon gradm from gradients
    
    gradm_sigma=steep_descent(model1,grad_sigma);
    gradm_epsilon=steep_descent(model1,grad_epsilon);

    figure;
    subplot(2,1,1)
    figuea=plot_wavefld( model1,gradm_sigma);
    title(['gradm sigma ' num2str(iter) 'th itearation'])
    subplot(2,1,2)
    figueb=plot_wavefld( model1,gradm_epsilon);
    title([' gradm epsilon  ' num2str(iter) 'th itearation'])

    %     saveas(gcf,sprintf('outputmyINV/figs/Gradm_iteration_%d_iter=%d.tiff',nw,iter));

    % ----------------------- SAVE mat FILES ----------------------------Save

    save('outputmyINV\gradm\gradm_sigma_old.mat','gradm_sigma');
    save('outputmyINV\gradm\gradm_epsilon_old.mat','gradm_epsilon');

    % --------------------------------------------------------------------------
    %                  (D)  Gauss quasi-Newton (GQN): l-BFGS method
    % --------------------------------------------------------------------------
    fprintf('LBFGS iter= %d .\n',iter);
    [Hgrad_sigma, Hgrad_epsilon] = LBFGS_TEmNEW(model1,sig_model,eps_model,hess_sigma,hess_epsilon,grad_sigma,grad_epsilon,iter);
    %     [Hgrad_sigma, Hgrad_epsilon] = LBFGS_TEmNEWmkl(model1,sig_model,eps_model,grad_sigma_tikh,grad_epsilon_tikh,iter);
    % [Hgrad_sigma(:,:,iter), Hgrad_epsilon(:,:,iter)] = LBFGS_TEm(model1,sig_model(:,:,iter),eps_model(:,:,iter),hess_sigma,hess_epsilon,grad_sigma,grad_epsilon,iter);
    clear sigmar epsilonr
    % oooooooooooooooooooooooooo (1) HEssian plotting ooooooooooooooooooooo
    figure
    set(gcf,'position',[0 0, model1.screenx, model1.screeny]);
    set(gcf,'PaperPositionMode','Auto');
    subplot(1,2,1)
    imagesc(model1.x,model1.y,Hgrad_sigma);
    % colormap (flipud(gray));
    load 'seismic.map'
    colormap (flipud(seismic));
    cbar_handle=colorbar('EastOutside');
    set(gca,'FontSize',12);
    %     grid on;
    axis ij;
    axis equal;
    axis tight;
    xlabel('Distance [m]');
    ylabel('Depth [m]');
    %     title('  Hgrad sigma');
    title(['Hgrad sigma iter= '  num2str(iter) 'th iteration']);

    subplot(1,2,2)
    imagesc(model1.x,model1.y, Hgrad_epsilon);
    % colormap (flipud(gray));
    load 'seismic.map'
    colormap (flipud(seismic));
    cbar_handle=colorbar('EastOutside');
    set(gca,'FontSize',12);
    %     grid on;
    axis ij;
    axis equal;
    axis tight;
    xlabel('Distance [m]');
    ylabel('Depth [m]');
    %     title(' Hgrad epsilon ');
    title(['Hgrad epsilon iter= '  num2str(iter) 'th iteration']);
    saveas(gcf,sprintf('outputmyINV/figs/04Hgrad_iteration_%d.tiff',iter));

    disp ('oooooooooooooooooooooooooo (2) check descent direction ooooooooooooooooooooo');

    % check if search direction is a descent direction, otherwise reset L-BFGS history
    nlbfgs_vec1=model1.nlbfgs_vec.*model1.nlbfgs;
    p_descent = check_descent_TE( model1,Hgrad_sigma,Hgrad_epsilon,grad_sigma_tikh,grad_epsilon_tikh);
    if (p_descent<0.0)
        disp('Loss of descent, reset l-bfgs history \n');

        for i=1:nlbfgs_vec1
            y_LBFGS(i)=0.0;
            s_LBFGS(i)=0.0;

        end


        for i=1:model1.nlbfgs_vec
            q_LBFGS(i)=0.0;
            r_LBFGS(i)=0.0;
        end

        for i=1:model1.nlbfgs

            alpha_LBFGS(i)=0.0;
            beta_LBFGS(i)=0.0;
            rho_LBFGS(i)=0.0;

        end

    end
    disp ('oooooooooooooooooooooooooo END: check descent direction ooooooooooooooooooooo\n');

    % --------------------------------------------------------------------------
    %                  (E)  Line-search Method: Wolfel condition
    % --------------------------------------------------------------------------

    disp('=============(E)  Line-search Method: Wolfel condition =======\n');
    if (model1.wolfe==1)
        [ ALPHA_OLD,sig_model,eps_model]= wolfe_TENEW(model1,acq,sig_model,eps_model,nshots, ntr,nw, precobs,grad_sigma,grad_epsilon,Hgrad_sigma,Hgrad_epsilon,w,iter,L2_tikhonov);
    else
        ALPHA_OLD =1.0;
    end
    fprintf('ALPHA_OLD = %e .\n',ALPHA_OLD);
    % ----------------------
    % %         [ ALPHA_OLD,sig_model,eps_model] = wolfe_TEmkl(model1,acq,sig_model,eps_model,nshots, ntr,nw, precobs,grad_sigma,grad_epsilon,Hgrad_sigma,Hgrad_epsilon,w,iter,L2_tikhonov);
    %     ALPHA_OLD = wolfe_TEjanM(model1,acq,sig_model(:,:,iter),eps_model(:,:,iter),nshots, ntr,nw, precobs,grad_sigma,grad_epsilon,Hgrad_sigma,Hgrad_epsilon,w,iter,L2_tikhonov,scale_grad_eps,SCALE_GRAD_eps_PRECOND,scale_grad_sig,SCALE_GRAD_sig_PRECOND);
    L2_hist(iter)=L2t(1);
    % --------------------------------------------------------------------------
    disp('=============apply optimum material parameter update -----');
    % % normalize material parameters
    sigmar1 = scale_grad_TE(sig_model ,model1.beta_sig./model1.sig0, model1);
    epsilonr1 = scale_grad_TE(eps_model,model1.beta_eps./model1.eps0, model1);

    %     sigmar1(:,:,iter)=scale_grad_TE(sig_model(:,:,iter),1.0./(model1.beta_sig.*model1.sig0), model1);
    %     epsilonr1(:,:,iter)=scale_grad_TE(eps_model(:,:,iter),1.0./(model1.beta_eps.*model1.eps0), model1);
    %
    %     [eps_scale1, eps_scale2]=alpha_meles(model1,sig_model, eps_model,Hgrad_sigma,Hgrad_epsilon);

    %     fprintf('ALPHA_OLD = %e .\n',ALPHA_OLD);
    % UPDATE sigmar and epsilonr
    eps_scale1 = ALPHA_OLD;
    eps_scale2=-ALPHA_OLD;
    %         fprintf('ALPHA_OLD = %e .\n',ALPHA_OLD);
    %         eps_scale1 = 1.0;

    sigmar = calc_mat_change_wolfe_multi_para_TE (sigmar1, model1,Hgrad_sigma,eps_scale1,1);
    epsilonr = calc_mat_change_wolfe_multi_para_TE (epsilonr1, model1,Hgrad_epsilon,eps_scale2,2);

    fprintf('ALPHA_OLD = %e .\n',ALPHA_OLD);
    clear sigma epsilon sigmar1 epsilonr1

    % convert sigmar to sigma & epsilonr

    sigma=scale_grad_TE(sigmar,(model1.sig0./model1.beta_sig), model1);
    epsilon=scale_grad_TE(epsilonr,(model1.eps0./model1.beta_eps), model1);


    %     sigma(:,:,iter)=scale_grad_TE(sigmar(:,:,iter),(model1.beta_sig.*model1.sig0), model1);
    %     epsilon(:,:,iter)=scale_grad_TE(epsilonr(:,:,iter),(model1.beta_eps.*model1.eps0), model1);

    close all;

    % oooooooooooooooooooooooooo (1)Models plotting ooooooooooooooooooooooooooo

    fh = figure();
    fh.WindowState = 'maximized';
    subplot (1,2,1)
    % waitforbuttonpress
    imagesc(model1.x,model1.y,epsilon./model1.eps0);
    caxis([model1.cx1 model1.cxe1]);
    cbar_handle=colorbar('EastOutside');
    set(get(cbar_handle,'ylabel'),'string','Relative permitivity','fontsize',model1.FSize,'FontWeight',model1.Fweight);
    load 'seismic.map'
    colormap (flipud(seismic));
    % colormap (flipud(gray));
    axis ij;
    axis equal;
    axis tight;
    set(gca,'FontSize',model1.FSize);
    xlabel('Distance [m]');
    ylabel('Depth [m]');
    title(['Output permitivity model ' num2str(iter) 'th iteration']);

    % waitforbuttonpress
    subplot (1,2,2)

    imagesc(model.x,model.y,sigma.*1000);
    caxis([model1.cx1 model1.cxs1]);
    cbar_handle=colorbar('EastOutside');
    set(get(cbar_handle,'ylabel'),'string','Conductivity(mS/m)','fontsize',model1.FSize,'FontWeight',model1.Fweight);
    % colormap (flipud(gray));
    load 'seismic.map'
    colormap (flipud(seismic));
    axis ij;
    axis equal;
    axis tight;
    set(gca,'FontSize',model1.FSize);
    xlabel('Distance [m]');
    ylabel('Depth [m]');
    title(['Output conductivity model ' num2str(iter) 'th iteration'])
    %save figure with iteration no.
    exportgraphics(gcf, sprintf('outputmyINV/figs/000Output_model_at_iteration= %d.tiff',iter),'Resolution',600)
    % saveas(gcf,sprintf('outputmyINV/figs/000Output_claylens_multioffset_at_iteration= %d.tiff',iter));
    clear fh
    %------------------------------------------------------------------------
    % calculating difference of the actual L2 and before two iterations, dividing with L2_hist[iter-2] provide changing in percent
    if (iter>1)
        model1.lbfgs_pointer=model1.lbfgs_pointer+1;

    end
    %     save('outputmyINV\wavefield\L2sumgrad_old.mat','L2sumgrad');
    if (iter>2)
        diff=abs((L2_hist(iter-2)-L2_hist(iter))./L2_hist(iter-2));
        ratio=L2_hist(iter)./L2_hist(1);
        ratio1(iter)=ratio;
    else
        diff=1e20;
        ratio=1.0e-1;
    end
    fprintf('diff = %e ratio = %e.\n',diff,ratio);

    save('outputmyINV\figs\FinalL2.mat','L2_hist');

    fprintf('diff = %e ratio = %e.\n',diff,ratio);
    if (iter>2)
        %--------------------
        %         LM2= load('E:\GPRFDFWI\outputmyINV\figs\L2_hist_2Cstag2_withwolfe_m02multioffsetbb.mat');
        %         m02cl=LM2.L2_hist;
        %         [xli1 a]=size(m02cl');
        %         if (iter<=xli1)
        %             xli=iter;
        %         else
        %             xli=xli1;
        %         end
        %--------------------
        figure;
        subplot(2,1,1)
        semilogy(ratio1);
        set(gca,'FontSize',20);
        xlabel('Iteration number');
        ylabel('ratio');
        title('ratio');
        grid on;

        %--------------------
        subplot(2,1,2)
        semilogy(L2_hist);
        %                 hold on;

        %         semilogy (1:xli,(m02cl(1:xli)),'r-');
        set(gca,'FontSize',12);
        xlabel('Iteration number');
        ylabel('L2');
        title('L2_hist');
        grid on;
        legend({'running'},'Location','northeast','FontSize',10,'FontName','Arial');
        title('new stag2 newcode april2019')
        %--------------------
        saveas(gcf,sprintf('outputmyINV/figs/#Output_2CStag1_log ratio_L2_.tiff'));
    end
    fprintf('diff = %e ratio = %e.\n',diff,ratio);
    PRO=0.0;
    % are Convergence criteria satisfied??
    % eps_scale =ALPHA_OLD;
    if((ratio<=5.0e-5)||(iter==itermax)||(diff<=PRO))% conditions ====>>>>>

        if (model1.npml>0)
            [epsilon,model1]  = extract_model_ext( model1, epsilon);
            [sigma,model]  = extract_model_ext( model, sigma);

        end
        fh = figure();
        fh.WindowState = 'maximized';
        subplot (1,2,1)
        % waitforbuttonpress
        imagesc(model1.x,model1.y,epsilon./model1.eps0);
        caxis([model1.cx1 model1.cxe1]);
        cbar_handle=colorbar('EastOutside');
        set(get(cbar_handle,'ylabel'),'string','Relative permitivity','fontsize',model1.FSize,'FontWeight',model1.Fweight);
        load 'seismic.map'
        colormap (flipud(seismic));
        % colormap (flipud(gray));
        axis ij;
        axis equal;
        axis tight;
        set(gca,'FontSize',12);
        xlabel('Distance [m]');
        ylabel('Depth [m]');
        title(['Output  Final  permitivity model ' num2str(iter) 'th iteration']);

        % waitforbuttonpress
        subplot (1,2,2)

        imagesc(model.x,model.y,sigma.*1000);
        caxis([model1.cx1 model1.cxs1]);
        %         caxis([0.0 model1.caxis1]);
        cbar_handle=colorbar('EastOutside');
        set(get(cbar_handle,'ylabel'),'string','Conductivity(mS/m)','fontsize',model.FSize,'FontWeight',model.Fweight);
        % colormap (flipud(gray));
        load 'seismic.map'
        colormap (flipud(seismic));
        axis ij;
        axis equal;
        axis tight;
        set(gca,'FontSize',12);
        xlabel('Distance [m]');
        ylabel('Depth [m]');
        title(['Output Final conductivity model ' num2str(iter) 'th iteration'])
        %save figure with iteration no.
        exportgraphics(gcf,sprintf('outputmyINV/figs/#Output_FINAL_Converged_Models_at_iteration= %d.tiff',iter),'Resolution',600)
        % saveas(gcf,sprintf('outputmyINV/figs/#Output_FINAL_Converged_Models_at_iteration= %d.tiff',iter));
        clear fh

        break;
    end

    fprintf('ALPHA_OLD = %e .\n',ALPHA_OLD);

    fprintf('diff = %e ratio = %e.\n',diff,ratio);
    %------------------------------------------------------------------------
    clear eps_model sig_model
    sig_model=sigma;
    eps_model=epsilon;
    save(['outputmyINV\models\output_eps_model_iter=',num2str(iter),'_.mat'],'eps_model');
    save(['outputmyINV\models\output_sig_model_iter=', num2str(iter),'_.mat'],'sig_model');

    clear sig_modelw sigmar res wiener te te_adj grad_sigma grad_epsilon

    disp('//////////////time rquired=');
    t=toc;
    disp(datestr(datenum(0,0,0,0,0,t),'HH:MM:SS'))
    disp('========================================');
    fprintf('diff = %e ratio = %e.\n',diff,ratio);
    fprintf('ALPHA_OLD = %e .\n',ALPHA_OLD);


end % END of iter loop

% matlabpool close
poolobj = gcp('nocreate');
delete(poolobj);

load chirp
sound(y,Fs)

disp('//////////////time rquired=');
t=toc;
disp(datestr(datenum(0,0,0,0,0,t),'HH:MM:SS'))
disp('========================================');