% 3D BCCB SPM employing calibration with dielectric scattering probe
% Version date: 16/06/2017
% Author: Denys Shumakov
% This code is a specific application of the Scattered-power Mapping that
% employs a BCCB solver. BCCB stands for block circulant with circulant
% blocks. The code allows for choosing alternative solvers as well.

%CODE REQUIREMENTS:

%1. The S-parameter raw data are assumed to be in the (x,y,freq,z) format.
%2. CO scan area should be of the same size as that of the OUT. 
%   Note that in this code, number of samples along x and y is odd.
%3. Spatial sampling along x and y has to be the same. 
%   Otherwise, trace function is not applicable for the power maps generation. 
%   Note that solvers work regardless of the sampling rate (system matrix is always square).
%4. CVX (optional for solver 2) has to be downloaded, installed and a 
%   license has to be obtained from cvxr.com.
%5. Path to SMT has to be added (optional for preconditioning in solver 3).
%6. Note that the code works with the data of single precision (optional for faster speed).

%P.S.: For the details, see [Shumakov and Nikolova "Fast quantitative microwave imaging
%      with scattered-power maps", IEEE Trans. MTT, 2017]
%      All references are to this publication, NOT the book, as the BCCB
%      solver is not explained in the book.

%% load the data

clear;clc
tic
% load reflection coefficient at port 1 - calibration object (CO) acquisition
load 'S11_CO';
%[S11_CO_signal,S11_CO_noise]=CO_denoising(S11_CO,0);
S1_CO=single(S11_CO);
clearvars S11_CO S11_CO_noise
% load reflection coefficient at port 2 - CO acquisition
load 'S22_CO';
%[S22_CO_signal,S22_CO_noise]=CO_denoising(S22_CO,0);
S2_CO=single(S22_CO);
clearvars S22_CO S22_CO_noise
% load transmission coefficient from port 2 to port 1 - CO acquisition
load 'S12_CO';
%[S12_CO_signal,S12_CO_noise]=CO_denoising(S12_CO,0);
S3_CO=single(S12_CO);
clearvars S12_CO S12_CO_noise
% load transmission coefficient from port 1 to port 2 - CO acquisition
load 'S21_CO';
%[S21_CO_signal,S21_CO_noise]=CO_denoising(S21_CO,0);
S4_CO=single(S21_CO);
clearvars S21_CO S21_CO_noise
% load reflection coefficient at port 1 - object under test (OUT) acquisition
load 'S11_OUT';
S1_OUT=single(S11_OUT);
clearvars S11_OUT
% load reflection coefficient at port 2 - OUT acquisition
load 'S22_OUT';
S2_OUT=single(S22_OUT);
clearvars S22_OUT
% load transmission coefficient from port 2 to port 1 - OUT acquisition
load 'S12_OUT';
S3_OUT=single(S12_OUT);
clearvars S12_OUT
% load transmission coefficient from port 1 to port 2 - OUT acquisition
load 'S21_OUT';
S4_OUT=single(S21_OUT);
clearvars S21_OUT
% load reflection coefficient at port 1 - reference object (RO) acquisition
load 'S11_RO';
%[S11_RO_signal,S11_RO_noise]=RO_denoising(S11_RO);
S1_RO=single(S11_RO);
clearvars S11_RO S11_RO_noise
% load reflection coefficient at port 2 - (RO) acquisition
load 'S22_RO';
%[S22_RO_signal,S22_RO_noise]=RO_denoising(S22_RO);
S2_RO=single(S22_RO);
clearvars S22_RO S22_RO_noise
% load transmission coefficient from port 2 to port 1 - RO acquisition
load 'S12_RO';
%[S12_RO_signal,S12_RO_noise]=RO_denoising(S12_RO);
S3_RO=single(S12_RO);
clearvars S12_RO S12_RO_noise
% load transmission coefficient from port 1 to port 2 - RO acquisition
load 'S21_RO';
%[S21_RO_signal,S21_RO_noise]=RO_denoising(S21_RO);
S4_RO=single(S21_RO);
clearvars S21_RO S21_RO_noise

[xnum1, ynum1, fnum, znum]=size(S1_CO); %xnum1, ynum1 and znum are number of samples along x,y and z; 
%fnum is number of frequency samples

%zero-padding S_COs to allow for spatial shifting
SCOtemp=zeros(2*xnum1-1, 2*ynum1-1, fnum, znum,'single');
[xnum, ynum, fnum, znum]=size(SCOtemp);
SCOtemp(ceil(0.25*xnum):ceil(0.75*xnum), ceil(0.25*ynum):ceil(0.75*ynum), :,:)=S1_CO; S1_CO=SCOtemp;
SCOtemp(ceil(0.25*xnum):ceil(0.75*xnum), ceil(0.25*ynum):ceil(0.75*ynum), :,:)=S2_CO; S2_CO=SCOtemp;
SCOtemp(ceil(0.25*xnum):ceil(0.75*xnum), ceil(0.25*ynum):ceil(0.75*ynum), :,:)=S3_CO; S3_CO=SCOtemp;
SCOtemp(ceil(0.25*xnum):ceil(0.75*xnum), ceil(0.25*ynum):ceil(0.75*ynum), :,:)=S4_CO; S4_CO=SCOtemp;
clearvars SCOtemp

%% initialization

f_st=1; %number of freq sample to start with
f_stop=fnum; %number of freq sample to finish with
f_range=f_st:f_stop; %set the frequency range to be processed
ROcontr=1; %known relative permittivity of RO (background medium)
COcontr=0.1; %known contrast of scattering probe in CO (i.e., the difference between 
% the relative permittivity of the scattering probe and that of the RO)
COsize=10; %in mm
scan_step=10; %in mm

% Choose type of forward model: Born or Rytov
approx=1; %Born is 1, Rytov is 0
solver=1; %choose the solver: BCCB 1, CVX 2, PCG 3, LSQLIN 4

scaling_factor=(COsize^2)/(scan_step^2); %calculating scaling factor for cubic scattering probe
CO_cx=(xnum+1)/2; %finding the center of the imaged region along x
CO_cy=(ynum+1)/2; %finding the center of the imaged region along y

%% generating OUT power maps. See Section II.C

range_x=ceil(0.25*xnum):ceil(0.75*xnum);
xnumm=length(range_x);
range_y=ceil(0.25*ynum):ceil(0.75*ynum);
ynumm=length(range_y);
cmap_out=zeros(xnum,ynum,znum,'single');
cmap_out_tmp=zeros(xnum,ynum,znum,'single');
cmap_co_tmp=zeros(xnum,ynum,znum,'single');
S1_COcal=zeros(length(range_x),length(range_y),fnum,znum,'single');
S2_COcal=zeros(length(range_x),length(range_y),fnum,znum,'single');
S3_COcal=zeros(length(range_x),length(range_y),fnum,znum,'single');
S4_COcal=zeros(length(range_x),length(range_y),fnum,znum,'single');
wx=(0:(2*pi)/xnumm:(2*pi*(xnumm-1))/xnumm)'; %calculating angular frequencies along x
wy=0:(2*pi)/ynumm:(2*pi*(ynumm-1))/ynumm; %calculating angular frequencies along y
fftw('planner','exhaustive'); %optimization of fft running speed
dS1=cell(CO_cx*CO_cy,length(f_range),znum); %initializing the cell with delta S1_COs
dS2=cell(CO_cx*CO_cy,length(f_range),znum); %initializing the cell with delta S2_COs
dS3=cell(CO_cx*CO_cy,length(f_range),znum); %initializing the cell with delta S3_COs
dS4=cell(CO_cx*CO_cy,length(f_range),znum); %initializing the cell with delta S4_COs
co_max=zeros(1,fnum);
counttf=1; %frequency counter
for i=f_range
    counttz=1; %range counter
    for t=1:znum
        
        if approx == 1 %calculating scattered responses within the LQL approximation
            deltaS1=S1_OUT(:, :, i)-S1_RO(:, :, i); 
            deltaS2=S2_OUT(:, :, i)-S2_RO(:, :, i); 
            deltaS3=S3_OUT(:, :, i)-S3_RO(:, :, i); 
            deltaS4=S4_OUT(:, :, i)-S4_RO(:, :, i); 
            S1_COcal(:,:,i,t)=S1_CO(range_x, range_y, i, t)-S1_RO(:, :, i);
            S2_COcal(:,:,i,t)=S2_CO(range_x, range_y, i, t)-S2_RO(:, :, i); 
            S3_COcal(:,:,i,t)=S3_CO(range_x, range_y, i, t)-S3_RO(:, :, i);
            S4_COcal(:,:,i,t)=S4_CO(range_x, range_y, i, t)-S4_RO(:, :, i);
        else %calculating scattered responses within the LQLR approximation
            deltaS1=S1_RO(:, :, i).*log(S1_OUT(:, :, i)./S1_RO(:, :, i));
            deltaS2=S2_RO(:, :, i).*log(S2_OUT(:, :, i)./S2_RO(:, :, i));
            deltaS3=S3_RO(:, :, i).*log(S3_OUT(:, :, i)./S3_RO(:, :, i));
            deltaS4=S4_RO(:, :, i).*log(S4_OUT(:, :, i)./S4_RO(:, :, i));
            S1_COcal(:,:,i,t)=S1_RO(:, :, i).*log(S1_CO(range_x, range_y, i, t)./S1_RO(:, :, i));
            S2_COcal(:,:,i,t)=S2_RO(:, :, i).*log(S2_CO(range_x, range_y, i, t)./S2_RO(:, :, i));
            S3_COcal(:,:,i,t)=S3_RO(:, :, i).*log(S3_CO(range_x, range_y, i, t)./S3_RO(:, :, i));
            S4_COcal(:,:,i,t)=S4_RO(:, :, i).*log(S4_CO(range_x, range_y, i, t)./S4_RO(:, :, i));
        end
        %applying 2D FFT to the CO data
        S1COF_origin=fft2(S1_COcal(:,:,i,t)); 
        S2COF_origin=fft2(S2_COcal(:,:,i,t)); 
        S3COF_origin=fft2(S3_COcal(:,:,i,t)); 
        S4COF_origin=fft2(S4_COcal(:,:,i,t)); 
        counttxy=1; %lateral (x and y) counter
        for ix=-(xnumm-1)/2:(xnumm-1)/2
            for iy=-(ynumm-1)/2:(ynumm-1)/2
                phaseshift_x=exp(-1i*wx*ix); %phase shift along x
                phaseshift_y=exp(-1i*wy*iy); %phase shift along y
                phaseshift=phaseshift_x*phaseshift_y; %phase shift along x and y. See Section II.E
                
                S1_CO_shifted_1=phaseshift.*S1COF_origin; %shifting the original delta S1_CO
                dS1{counttxy,counttf,counttz}=ifft2(S1_CO_shifted_1); %return to spatial domain
                
                S2_CO_shifted_1=phaseshift.*S2COF_origin; %shifting the original delta S2_CO
                dS2{counttxy,counttf,counttz}=ifft2(S2_CO_shifted_1); %return to spatial domain
                
                S3_CO_shifted_1=phaseshift.*S3COF_origin; %shifting the original delta S3_CO
                dS3{counttxy,counttf,counttz}=ifft2(S3_CO_shifted_1); %return to spatial domain
                
                S4_CO_shifted_1=phaseshift.*S4COF_origin; %shifting the original delta S4_CO
                dS4{counttxy,counttf,counttz}=ifft2(S4_CO_shifted_1); %return to spatial domain
                
                %calculating the OUT and CO power maps for each range location and each frequency sample                
                cmap_out_tmp(xnumm+ix, ynumm+iy, t)=cmap_out_tmp(xnumm+ix, ynumm+iy, t)+trace(deltaS1*[dS1{counttxy,counttf,counttz}]')+trace(deltaS2*[dS2{counttxy,counttf,counttz}]')+trace(deltaS3*[dS3{counttxy,counttf,counttz}]')+trace(deltaS4*[dS4{counttxy,counttf,counttz}]');
                cmap_co_tmp(xnumm+ix, ynumm+iy, t)=cmap_co_tmp(xnumm+ix, ynumm+iy, t)+trace([dS1{counttxy,counttf,counttz}]*[dS1{counttxy,counttf,counttz}]')+trace([dS2{counttxy,counttf,counttz}]*[dS2{counttxy,counttf,counttz}]')+trace([dS3{counttxy,counttf,counttz}]*[dS3{counttxy,counttf,counttz}]')+trace([dS4{counttxy,counttf,counttz}]*[dS4{counttxy,counttf,counttz}]');
                counttxy=counttxy+1; %updating the counter
            end
        end               
        counttz=counttz+1; %updating the counter
    end
    co_max(:,i)=max(max(max(abs(cmap_co_tmp)))); %finding the CO power map maximum over all range locations per each frequency sample    
    for t=1:znum
        cmap_out(:,:,t)=cmap_out(:,:,t)+cmap_out_tmp(:,:,t)./co_max(:,i); %frequency combination of the normalized OUT power maps. See Section II.D
    end
    cmap_out_tmp=zeros(xnum,ynum,znum,'single');
    counttf=counttf+1; %updating the counter
end

clearvars S1_COcal S2_COcal phaseshift cmap_out_tmp S1COF_origin S2COF_origin %freeing the memory

%% CO power maps generator. See Section II.C

CO_Cmap_temp=zeros(xnum, ynum, znum,'single');
cmap_co_tmp=zeros(xnum,ynum,znum,'single');
CO_Cmap_Store1=cell(znum,1);
ixx=-(xnumm-1)/2; %shifting value along x
iyy=-(ynumm-1)/2; %shifting value along y
countz=1; %counter along z

for tt=1:znum       
    
    disp('Creating Map')
    disp(['ixx=' num2str(ixx) ', iyy=' num2str(iyy)])
    counttf=1; %frequency counter
    for i=f_range
        
        deltaS1=[dS1{1,counttf,countz}]; %using the previously calculated CO power maps       
        deltaS2=[dS2{1,counttf,countz}];  %using the previously calculated CO power maps       
        deltaS3=[dS3{1,counttf,countz}]; %using the previously calculated CO power maps         
        deltaS4=[dS4{1,counttf,countz}]; %using the previously calculated CO power maps 
        
        counttz=1; %range counter
        for t=1:znum            
            counttxy=1; %lateral (x and y) counter
            for ix=-(xnumm-1)/2:(xnumm-1)/2
                for iy=-(ynumm-1)/2:(ynumm-1)/2
                    
                    %calculating the CO power maps for each range location and each frequency sample                
                    cmap_co_tmp(xnumm+ix, ynumm+iy, t)=cmap_co_tmp(xnumm+ix, ynumm+iy, t)+trace(deltaS1*[dS1{counttxy,counttf,counttz}]')+trace(deltaS2*[dS2{counttxy,counttf,counttz}]')+trace(deltaS3*[dS3{counttxy,counttf,counttz}]')+trace(deltaS4*[dS4{counttxy,counttf,counttz}]');
                    counttxy=counttxy+1; %updating the counter
                end
            end
            counttz=counttz+1; %updating the counter
        end
                
        CO_Cmap_temp=CO_Cmap_temp+cmap_co_tmp./co_max(:,i); %frequency combination of the normalized CO power maps. See Section II.D
        cmap_co_tmp=zeros(xnum,ynum,znum,'single');
        counttf=counttf+1; %updating the counter
    end    
    CO_Cmap_Store1{tt}=CO_Cmap_temp(ceil(0.25*xnum):ceil(0.75*xnum), ceil(0.25*ynum):ceil(0.75*ynum),:); %storing the CO power maps
    CO_Cmap_temp=zeros(xnum, ynum, znum,'single');
    disp('Map is Done')
    countz=countz+1; %updating the counter
end

clearvars CO_Cmap_temp cmap_co_tmp deltaS1 deltaS2 %freeing the memory

%% forming matrix A (only first column) and data vector b. See Section II.C

CO_Cmap1=zeros(xnumm,ynumm,znum,CO_cx*CO_cy*znum,'single');
counter1=0; %counter initialization
for t=1:znum
    for mm=ceil(0.25*xnum):ceil(0.75*xnum)
        for nn=ceil(0.25*ynum):ceil(0.75*ynum)
            counter1=counter1+1; %counter update                       
            CO_Cmap1(:,:,:,counter1)=[CO_Cmap_Store1{t}]; %arranging the CO power maps' values in accordance with a counter
        end
    end
end

A=zeros(CO_cx*CO_cy*znum,CO_cx*CO_cy*znum,'single');
for hh=1 %this counter is 1 since we need only the first column. See Section II.F
    ff=1;
    for t=1:znum        
        for m=1:CO_cx
            for n=1:CO_cy                
                A(ff, hh)=CO_Cmap1(m, n, t, hh); %forming the first column of system matrix
                ff=ff+1;
            end
        end
    end
end

b=zeros(CO_cx*CO_cy*znum,1,'single');
ff=1;
for t=1:znum
    for m=1:CO_cx
        for n=1:CO_cy
            b(ff,1)=cmap_out(m+(xnum-1)/4, n+(ynum-1)/4, t); %forming the data vector 
            ff=ff+1;
        end
    end
end

clearvars CO_Cmap1

%% solving linear system of equations Ax=b

if solver == 1 %BCCB solver. See Section II.F
    disp('BCCB solver used')       
    AA=reshape(A(:,1),[CO_cx CO_cy znum]);
    bb=reshape(b(:,1),[CO_cx CO_cy znum]);
    kernel=fftn(bb)./fftn(AA);
    proper_sol=scaling_factor.*COcontr.*flipud(rot90(ifftn(kernel),-3))+ROcontr; %flipud and rot90 are used to orient the quantitative images accordingly to the qualitative one   
    
elseif solver == 2 %cvx solver (allows for constraints). Taken from: cvxr.com
    disp('CVX solver used')      
    Fmtx_x=dftmtx(CO_cx);
    Fmtx_y=dftmtx(CO_cy);
    Fmtx_z=dftmtx(znum);
    Fmtx=kron(kron(Fmtx_z,Fmtx_y),Fmtx_x);
    Fmt_adj=Fmtx'./(CO_cx*CO_cy*znum);
    BCCB=Fmt_adj*diag(reshape(fftn(reshape(A(:,1),[CO_cx CO_cy znum])),[CO_cx*CO_cy*znum 1]))*Fmtx; %forming the full system matrix A
    
    %re-formulating the problem to incorporate physicality constraints
    e_ro=ones(length(b),1)*ROcontr;
    m_hat=double(BCCB)*e_ro+scaling_factor*COcontr*double(b);
    N=[real(double(BCCB)) -imag(double(BCCB));imag(double(BCCB)) real(double(BCCB))];
    K=[real(m_hat); imag(m_hat)];
    Nv=length(b);
    AAA=[-eye(Nv) zeros(Nv); zeros(Nv) eye(Nv)];
    bbb=[-ones(Nv,1); zeros(Nv,1)];
    
    run('\\srv-cem-fs\RedirectedFolders$\shumakds\Desktop\cvx-w64\cvx\cvx_startup.m')
    cvx_begin
    cvx_solver Mosek %SDPT3 and sedumi are other available solvers within CVX package
    variable lsqcoe_solution(2*Nv)
    minimize (norm(N*lsqcoe_solution-K))
    subject to
    AAA*lsqcoe_solution <= bbb %using physicality constraints. See Section III.C
    cvx_end
    solution=lsqcoe_solution(1:Nv)+1i*lsqcoe_solution(Nv+1:end);
    proper_sol=reshape(solution,xnum1,ynum1,znum);
    
elseif solver == 3  %preconditioned or not conjugate gradients solver  
    disp('Conjugate gradients used')
    Fmtx_x=dftmtx(CO_cx);
    Fmtx_y=dftmtx(CO_cy);
    Fmtx_z=dftmtx(znum);
    Fmtx=kron(kron(Fmtx_z,Fmtx_y),Fmtx_x);
    Fmt_adj=Fmtx'./(CO_cx*CO_cy*znum);
    BCCB=Fmt_adj*diag(reshape(fftn(reshape(A(:,1),[CO_cx CO_cy znum])),[CO_cx*CO_cy*znum 1]))*Fmtx; %forming the full system matrix A. See Section II.F
        
    oldpath=addpath('H:\MATLAB\smt','-end'); %enabling SMT
    L=smtcprec('superopt',double(BCCB)); %using circulant preconditioner
    [pcg_solution, flag, res, iter]=pcg(double(BCCB),double(b),[],500,L);    
    proper_sol=reshape(scaling_factor*COcontr*pcg_solution+ROcontr,xnum1,ynum1,znum);
    
else %linear least-squares solver (allows for constraints)
    disp('Linear least-squares used')
    Fmtx_x=dftmtx(CO_cx);
    Fmtx_y=dftmtx(CO_cy);
    Fmtx_z=dftmtx(znum);
    Fmtx=kron(kron(Fmtx_z,Fmtx_y),Fmtx_x);
    Fmt_adj=Fmtx'./(CO_cx*CO_cy*znum);
    BCCB=Fmt_adj*diag(reshape(fftn(reshape(A(:,1),[CO_cx CO_cy znum])),[CO_cx*CO_cy*znum 1]))*Fmtx; %forming the full system matrix A. See Section II.F
    
    %re-formulating the problem to incorporate physicality constraints
    e_ro=ones(length(b),1)*ROcontr;
    m_hat=double(BCCB)*e_ro+scaling_factor*COcontr*double(b);
    N=[real(double(BCCB)) -imag(double(BCCB));imag(double(BCCB)) real(double(BCCB))];
    K=[real(m_hat); imag(m_hat)];
    Nv=length(b);
    AAA=[-eye(Nv) zeros(Nv); zeros(Nv) eye(Nv)];
    bbb=[-ones(Nv,1); zeros(Nv,1)];
        
    tmp=double(b)./diag(double(BCCB)); %setting initial point according to Sheng's paper
    start=[real(tmp);imag(tmp)];
    options = optimset('MaxIter',5000);
    [lsqcoe_solution,resnorm,residual,exitflag,output] =lsqlin(N, K, AAA, bbb,[],[],[],[],start,options); %using physicality constraints. See Section III.C
    disp(exitflag)
    solution=lsqcoe_solution(1:Nv)+1i*lsqcoe_solution(Nv+1:end);
    proper_sol=reshape(solution,xnum1,ynum1,znum);
    
end

toc

%% Plotting the results

int=0; %interpolation yes 1, no 0. Note that the colormap values may change during certain interpolation settings
int_rate=0.2; %set interpolation rate

if znum == 1 %plotting 2D images (this section can be also run to get the 3D images slice by slice separately: that is why loop over range locations is still there)
    
    for t=1:znum
        figure('Units', 'centimeters', 'OuterPosition',[0 10 70 14]);
        
        subplot(1, 3, 1) %qualitative outcome
        drawmap=cmap_out(ceil(0.25*xnum):ceil(0.75*xnum), ceil(0.25*ynum):ceil(0.75*ynum), t);
        
        V=abs(drawmap/max(max(drawmap))); %normalized modulus
        imagesc(0:(xnum1-1)*scan_step/10, 0:(ynum1-1)*scan_step/10,V);
        colorbar
        colormap parula
        shading interp
        alpha(0.9);
        
        subplot(1, 3, 2) %real part of the relative dielectric permittivity distribution        
        drawmap=proper_sol(:, :, t);        
        if int == 1
            [X,Y]=meshgrid(1:CO_cx,1:CO_cy);
            [XI,YI]=meshgrid(1:int_rate:CO_cx,1:int_rate:CO_cy);
            V2=interp2(X,Y,drawmap,XI,YI,'spline');
            imagesc(0:(xnum1-1)*scan_step/10, 0:(ynum1-1)*scan_step/10,real(V2));
        else
            V=real(drawmap);
            imagesc(0:(xnum1-1)*scan_step/10, 0:(ynum1-1)*scan_step/10,V);
        end
        colorbar
        colormap parula
        shading interp
        alpha(0.9);
        
        subplot(1, 3, 3) %imaginary part of the relative dielectric permittivity distribution
        if int == 1
            [X,Y]=meshgrid(1:CO_cx,1:CO_cy);
            [XI,YI]=meshgrid(1:int_rate:CO_cx,1:int_rate:CO_cy);
            V2=interp2(X,Y,drawmap,XI,YI,'spline');
            imagesc(0:(xnum1-1)*scan_step/10, 0:(ynum1-1)*scan_step/10,imag(V2));
        else
            V=imag(drawmap);
            imagesc(0:(xnum1-1)*scan_step/10, 0:(ynum1-1)*scan_step/10,V);
        end
        colorbar
        colormap parula
        shading interp
        alpha(0.9);
        
        set(gcf, 'color','w');
    end
    
else % plotting 3D images
    
    if int == 1
        Z=ones(xnum1,ynum1);
        [X,Y]=meshgrid(1:xnum1,1:ynum1);
        [X1,Y1]=meshgrid(linspace(0,(xnum1-1)*scan_step/10,xnum1),linspace(0,(ynum1-1)*scan_step/10,ynum1));
        
        figure
        for k=1:znum
            [XI,YI]=meshgrid(1:int_rate:xnum1,1:int_rate:ynum1);
            V3(:,:,k)=interp2(X,Y,proper_sol(:, :, k), XI,YI);
            surface(X1,Y1,Z.*k,real(V3(:,:,k)),'FaceColor','texturemap','CDataMapping','scaled')
        end
        colormap parula
        view(3), box on, axis tight square
        shading flat
        set(gca,'YDir','reverse','ZLim',[1 znum])
        colorbar
        
        figure
        for k=1:znum
            [XI,YI]=meshgrid(1:int_rate:xnum1,1:int_rate:ynum1);
            V3(:,:,k)=interp2(X,Y,proper_sol(:, :, k), XI,YI);
            surface(X1,Y1,Z.*k,imag(V3(:,:,k)),'FaceColor','texturemap','CDataMapping','scaled')
        end
        colormap parula
        view(3), box on, axis tight square
        shading flat
        set(gca,'YDir','reverse','ZLim',[1 znum])
        colorbar
        
    else
        Z=ones(xnum1,ynum1);
        [X,Y]=meshgrid(1:xnum1,1:ynum1);
        figure
        for k=1:znum
            surface(X,Y,Z.*k,real(proper_sol(:,:,k)),'FaceColor','texturemap','CDataMapping','scaled')
        end
        colormap parula
        view(3), box on, axis tight square
        shading flat
        set(gca,'YDir','reverse','ZLim',[1 znum])
        colorbar
        
        figure
        for k=1:znum
            surface(X,Y,Z.*k,imag(proper_sol(:,:,k)),'FaceColor','texturemap','CDataMapping','scaled')
        end
        colormap parula
        view(3), box on, axis tight square
        shading flat
        set(gca,'YDir','reverse','ZLim',[1 znum])
        colorbar
    end
end