function plot_RAPTOR_LIUQE(L,vk,model, ... % not iterated
                LYiim1, gkiim1, xkiim1, xdotkm1,... % previous iteartion at the end of the subrelaxatin process
                LYii,   gkii_after_relax, xkii_after_relax, xdotkii_after_relax,  ... % corresponding to gkii, xkii, 
                        gkii_before_relax, xkii_before_relax, xdotkii_before_relax ) % computed before the relaxation.
persistent hax hf

% List of quantities to be plot
% The subplot are filled by raw
plot_requested = { 'pprime', 'ttprime', 'jphioR' , ...
      'pprime_diff', 'ttprime_diff',  'Q0Q_diff', 'Q1Q_diff', 'Q2Q_diff', 'Q3Q_diff', 'Q4Q_diff','iqQ_diff'};
nplots = numel(plot_requested);

nraw = 2;
ncolums = ceil(nplots/nraw);

% Initilize figure and get axis handle
if isempty(hf)
    % initialize axes
    hf=figure();
    set(hf,'name','RAPTOR_LIUQE iterations');
    set(hf,'units','inch','position',[2,2,11,9]) 
    
    % Subplot are filled by columns
    for ii=1:nplots
        hax(ii) = subplot(nraw, ncolums,ii,'parent',hf);
    end
end


for ii=1:nplots
    ax = hax(ii);
    cla(ax);
    switch plot_requested{ii}
        case 'psicontour'
            FxN = (LYii.Fx-LYii.FA)/(LYii.FB-LYii.FA);
            contourf(ax,L.rx,L.zx,FxN,21); hold(ax,'on');
            contour(ax,L.rx,L.zx,FxN,[1 1],'w','linewidth',2);
            plot(ax,LYii.rA,LYii.zA,'ow');
            axis(ax,'equal');
            
        case 'pprime'
            hold(ax,'on')
            
            [pprime_LIUQE_input, ~ ,rhopsiN_LIUQE_input] = calc_pp_ttp( xkiim1,gkiim1,vk,model);

            pprime_LIUQE_output = LYii.PpQ; 
            pprime_RAPTOR_iim1 = eval_pprime(xkiim1,[],vk,model,false);
            pprime_RAPTOR_ii_before_relax = eval_pprime(xkii_before_relax,[],vk,model,false);
            pprime_RAPTOR_ii_afert_relax = eval_pprime(xkii_after_relax,[],vk,model,false);
             
            psiN_LIUQE_input = rhopsiN_LIUQE_input.^2;
            psiN_LIUQE_output = L.pQ.^2;
            psiN_RAPTOR_iim1 =  eval_psiN(xkiim1,[],[],model,false);
            psiN_RAPTOR_ii_before_relax = eval_psiN(xkii_before_relax,[],[],model,false);
            psiN_RAPTOR_ii_afert_relax = eval_psiN(xkii_after_relax,[],[],model,false);
             
            
            plot(ax, psiN_LIUQE_input, pprime_LIUQE_input, 'DisplayName', 'LIUQE input');
            plot(ax, psiN_LIUQE_output, pprime_LIUQE_output, 'DisplayName','LIUQE output');
            plot(ax, psiN_RAPTOR_iim1, pprime_RAPTOR_iim1, 'DisplayName', 'RAPTOR iim1');
            plot(ax, psiN_RAPTOR_ii_before_relax, pprime_RAPTOR_ii_before_relax, 'DisplayName', 'RAPTOR ii before relax');
            plot(ax, psiN_RAPTOR_ii_afert_relax, pprime_RAPTOR_ii_afert_relax, 'DisplayName', 'RAPTOR ii after relax');

    
            xlim(ax, [-0.01, 1])
            ylabel(ax,'pprime')
            xlabel(ax,'psiN')
            legend(ax,'show', 'location', 'south')
            
        case 'pprime_zoom_axes'
            hold(ax,'on')
            
            [pprime_LIUQE_input, ~ ,rhopsiN_LIUQE_input] = calc_pp_ttp( xkiim1,gkiim1,vk,model);
            
            pprime_LIUQE_output = LYii.PpQ;
            pprime_RAPTOR_iim1 = eval_pprime(xkiim1,[],vk,model,false);
            pprime_RAPTOR_ii_before_relax = eval_pprime(xkii_before_relax,[],vk,model,false);
            pprime_RAPTOR_ii_afert_relax = eval_pprime(xkii_after_relax,[],vk,model,false);
            
            psiN_LIUQE_input = rhopsiN_LIUQE_input.^2;
            psiN_LIUQE_output = L.pQ.^2;
            psiN_RAPTOR_iim1 =  eval_psiN(xkiim1,[],[],model,false);
            psiN_RAPTOR_ii_before_relax = eval_psiN(xkii_before_relax,[],[],model,false);
            psiN_RAPTOR_ii_afert_relax = eval_psiN(xkii_after_relax,[],[],model,false);
            
            
            plot(ax, psiN_LIUQE_input, pprime_LIUQE_input, 'DisplayName', 'LIUQE input');
            plot(ax, psiN_LIUQE_output, pprime_LIUQE_output, 'DisplayName','LIUQE output');
            plot(ax, psiN_RAPTOR_iim1, pprime_RAPTOR_iim1, 'DisplayName', 'RAPTOR iim1');
            plot(ax, psiN_RAPTOR_ii_before_relax, pprime_RAPTOR_ii_before_relax, 'DisplayName', 'RAPTOR ii before relax');
            plot(ax, psiN_RAPTOR_ii_afert_relax, pprime_RAPTOR_ii_afert_relax, 'DisplayName', 'RAPTOR ii after relax');
            
            xlim(ax,[-0.009 0.1])
            ylabel(ax,'pprime')
            xlabel(ax,'psiN')
            legend(ax,'show', 'location', 'south')
            
            
        case 'ttprime'
            hold(ax,'on')

            [~, ttprime_LIUQE_input,rhopsiN_LIUQE_input] = calc_pp_ttp( xkiim1,gkiim1,vk,model);
            
            ttprime_LIUQE_output = LYii.TTpQ;
            ttprime_RAPTOR_iim1 = eval_ttprime(xkiim1,gkiim1,vk,model,false);
            ttprime_RAPTOR_ii_before_relax = eval_ttprime(xkii_before_relax,gkii_before_relax,vk,model,false);
            ttprime_RAPTOR_ii_afert_relax = eval_ttprime(xkii_after_relax,gkii_after_relax,vk,model,false);
            
            psiN_LIUQE_input = rhopsiN_LIUQE_input.^2;
            psiN_LIUQE_output = L.pQ.^2;
            psiN_RAPTOR_iim1 =  eval_psiN(xkiim1,[],[],model,false);
            psiN_RAPTOR_ii_before_relax = eval_psiN(xkii_before_relax,[],[],model,false);
            psiN_RAPTOR_ii_afert_relax = eval_psiN(xkii_after_relax,[],[],model,false);
            
            
            plot(ax, psiN_LIUQE_input, ttprime_LIUQE_input, 'DisplayName', 'LIUQE input');
            plot(ax, psiN_LIUQE_output, ttprime_LIUQE_output, 'DisplayName','LIUQE output');
            plot(ax, psiN_RAPTOR_iim1, ttprime_RAPTOR_iim1, 'DisplayName', 'RAPTOR iim1');
            plot(ax, psiN_RAPTOR_ii_before_relax, ttprime_RAPTOR_ii_before_relax, 'DisplayName', 'RAPTOR ii before relax');
            plot(ax, psiN_RAPTOR_ii_afert_relax, ttprime_RAPTOR_ii_afert_relax, 'DisplayName', 'RAPTOR ii after relax');
            
            xlim(ax,[-0.01 1.])            
            ylabel(ax,'ttprime')
            xlabel(ax,'psiN')
            legend(ax,'show', 'location', 'south')
            
                 
        case 'ttprime_zoom_axes'
            hold(ax,'on')
            
            [~, ttprime_LIUQE_input,rhopsiN_LIUQE_input] = calc_pp_ttp( xkiim1,gkiim1,vk,model);
            
            ttprime_LIUQE_output = LYii.TTpQ;
            ttprime_RAPTOR_iim1 = eval_ttprime(xkiim1,gkiim1,vk,model,false);
            ttprime_RAPTOR_ii_before_relax = eval_ttprime(xkii_before_relax,gkii_before_relax,vk,model,false);
            ttprime_RAPTOR_ii_afert_relax = eval_ttprime(xkii_after_relax,gkii_after_relax,vk,model,false);
            
            psiN_LIUQE_input = rhopsiN_LIUQE_input.^2;
            psiN_LIUQE_output = L.pQ.^2;
            psiN_RAPTOR_iim1 =  eval_psiN(xkiim1,[],[],model,false);
            psiN_RAPTOR_ii_before_relax = eval_psiN(xkii_before_relax,[],[],model,false);
            psiN_RAPTOR_ii_afert_relax = eval_psiN(xkii_after_relax,[],[],model,false);
            
            
            plot(ax, psiN_LIUQE_input, ttprime_LIUQE_input, 'DisplayName', 'LIUQE input');
            plot(ax, psiN_LIUQE_output, ttprime_LIUQE_output, 'DisplayName','LIUQE output');
            plot(ax, psiN_RAPTOR_iim1, ttprime_RAPTOR_iim1, 'DisplayName', 'RAPTOR iim1');
            plot(ax, psiN_RAPTOR_ii_before_relax, ttprime_RAPTOR_ii_before_relax, 'DisplayName', 'RAPTOR ii before relax');
            plot(ax, psiN_RAPTOR_ii_afert_relax, ttprime_RAPTOR_ii_afert_relax, 'DisplayName', 'RAPTOR ii after relax');
            
            
            xlim(ax,[-0.009 0.1])            
            ylabel(ax,'ttprime')
            xlabel(ax,'psiN')
            legend(ax,'show', 'location', 'south')
        case 'jphioR'
            hold(ax,'on')
            %TODO
            % Add LIUQE plot 

            jphioR_RAPTOR_iim1 = eval_jphioR(xkiim1,gkiim1,vk,model,false);
            jphioR_RAPTOR_ii_before_relax = eval_jphioR(xkii_before_relax,gkii_before_relax,vk,model,false);
            jphioR_RAPTOR_ii_afert_relax = eval_jphioR(xkii_after_relax,gkii_after_relax,vk,model,false);
            

            psiN_RAPTOR_iim1 =  eval_psiN(xkiim1,[],[],model,false);
            psiN_RAPTOR_ii_before_relax = eval_psiN(xkii_before_relax,[],[],model,false);
            psiN_RAPTOR_ii_afert_relax = eval_psiN(xkii_after_relax,[],[],model,false);
            
            

            plot(ax, psiN_RAPTOR_iim1, jphioR_RAPTOR_iim1, 'DisplayName', 'RAPTOR iim1');
            plot(ax, psiN_RAPTOR_ii_before_relax, jphioR_RAPTOR_ii_before_relax, 'DisplayName', 'RAPTOR ii before relax');
            plot(ax, psiN_RAPTOR_ii_afert_relax, jphioR_RAPTOR_ii_afert_relax, 'DisplayName', 'RAPTOR ii after relax');
            
            xlim(ax,[-0.01 1])
            ylabel(ax,'jphioR')
            xlabel(ax,'psiN')
            legend(ax,'show', 'location', 'south')
            
        case 'g23or'
            hold(ax,'on')
            %TODO
            % Add LIUQE plot 

            g23or_RAPTOR_iim1 = eval_g23or(xkiim1,gkiim1,vk,model,false);
            g23or_RAPTOR_ii_before_relax = eval_g23or(xkii_before_relax,gkii_before_relax,vk,model,false);
            g23or_RAPTOR_ii_afert_relax = eval_g23or(xkii_after_relax,gkii_after_relax,vk,model,false);
            

            psiN_RAPTOR_iim1 =  eval_psiN(xkiim1,[],[],model,false);
            psiN_RAPTOR_ii_before_relax = eval_psiN(xkii_before_relax,[],[],model,false);
            psiN_RAPTOR_ii_afert_relax = eval_psiN(xkii_after_relax,[],[],model,false);
            
            

            plot(ax, psiN_RAPTOR_iim1, g23or_RAPTOR_iim1, 'DisplayName', 'RAPTOR iim1');
            plot(ax, psiN_RAPTOR_ii_before_relax, g23or_RAPTOR_ii_before_relax, 'DisplayName', 'RAPTOR ii before relax');
            plot(ax, psiN_RAPTOR_ii_afert_relax, g23or_RAPTOR_ii_afert_relax, 'DisplayName', 'RAPTOR ii after relax');
            
            xlim(ax,[-0.01 1])            
            ylabel(ax,'g23or')
            xlabel(ax,'psiN')
            legend(ax,'show', 'location', 'south')
            
            
        case  'pprime_diff'
            hold(ax,'on')
            
            pprime_RAPTOR_iim1 = eval_pprime(xkiim1,[],vk,model,false);
            pprime_RAPTOR_ii_afert_relax = eval_pprime(xkii_after_relax,[],vk,model,false);

            plot(ax, model.geom.rho, pprime_RAPTOR_iim1-pprime_RAPTOR_ii_afert_relax , 'DisplayName', 'RAPTOR iim1');
            
            xlim(ax, [-0.01, 1])
            ylabel(ax,'pprime_{im1} - pprime_{i}')
            xlabel(ax,'rhotor')
            legend(ax,'show', 'location', 'south')
            
            
            
        
        case  'ttprime_diff'
            hold(ax,'on')
            
            
            ttprime_RAPTOR_iim1 = eval_ttprime(xkiim1,gkiim1,vk,model,false);
            ttprime_RAPTOR_ii_afert_relax = eval_ttprime(xkii_after_relax,gkii_after_relax,vk,model,false);
            
            plot(ax, model.geom.rho, ttprime_RAPTOR_iim1-ttprime_RAPTOR_ii_afert_relax, 'DisplayName', 'RAPTOR iim1');
            
            xlim(ax,[-0.01 1.])
            ylabel(ax,'ttprime_{im1}-ttprime_{i}')
            xlabel(ax,'rhotor')
            legend(ax,'show', 'location', 'south')
                        
        case  'Q0Q_diff'
            try
            hold(ax,'on')
            
            plot(ax, L.pQ , LYiim1.Q0Q- LYii.Q0Q, 'DisplayName', 'RAPTOR iim1');
            
            xlim(ax,[-0.01 1.])
            ylabel(ax,'Q0Qiim1 - Q0Qii')
            xlabel(ax,'rhotor')
            legend(ax,'show', 'location', 'south')
            catch
            end
        case  'Q1Q_diff'
            try
            hold(ax,'on')
            
            plot(ax, L.pQ , LYiim1.Q1Q- LYii.Q1Q, 'DisplayName', 'RAPTOR iim1');
            
            xlim(ax,[-0.01 1.])
            ylabel(ax,'Q1Qiim1 - Q1Qii')
            xlabel(ax,'rhotor')
            legend(ax,'show', 'location', 'south')
            catch
            end
        case  'Q2Q_diff'
            try
            hold(ax,'on')
            
            plot(ax, L.pQ , LYiim1.Q2Q- LYii.Q2Q, 'DisplayName', 'RAPTOR iim1');
            
            xlim(ax,[-0.01 1.])
            ylabel(ax,'Q2Qiim1 - Q2Qii')
            xlabel(ax,'rhotor')
            legend(ax,'show', 'location', 'south')
            catch
            end
        case  'Q3Q_diff'
            try
            hold(ax,'on')
            
            plot(ax, L.pQ , LYiim1.Q3Q- LYii.Q3Q, 'DisplayName', 'RAPTOR iim1');
            
            xlim(ax,[-0.01 1.])
            ylabel(ax,'Q3Qiim1 - Q3Qii')
            xlabel(ax,'rhotor')
            legend(ax,'show', 'location', 'south')
            catch
            end
        case  'Q4Q_diff'
            try
            hold(ax,'on')
            
            plot(ax, L.pQ , LYiim1.Q4Q- LYii.Q4Q, 'DisplayName', 'RAPTOR iim1');
            
            xlim(ax,[-0.01 1.])
            ylabel(ax,'Q4Qiim1 - Q4Qii')
            xlabel(ax,'rhotor')
            legend(ax,'show', 'location', 'south')
            catch
            end
        case  'iqQ_diff'
            try
            hold(ax,'on')
            
            plot(ax, L.pQ , LYiim1.iqQ- LYii.iqQ, 'DisplayName', 'RAPTOR iim1');
            
            xlim(ax,[-0.01 1.])
            ylabel(ax,'iqQiim1 - iqQii')
            xlabel(ax,'rhotor')
            legend(ax,'show', 'location', 'south')
            catch
            end
    
        otherwise
            warning('\n Plot requested not recognised\n \n ')
            
    end
    
    
    
    
end

drawnow
