Graphene_antidot_transport.m

Go to the documentation of this file.

%    Magnetic transport through a quantum antidot
%    Copyright (C) 2015 Peter Rakyta, Ph.D.
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see http://www.gnu.org/licenses/.
%
%> @addtogroup Examples
%> @{
%> @file Graphene_antidot_transport.m
%> @brief Calculates the conductance through a graphene antidot in a homogeneous magnetic field.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
%> @Available
%> EQuUs v4.7 or later
%> @expectedresult
%> @image html Graphene_antidot_transport.jpg
%> @image latex Graphene_antidot_transport.jpg
%> @}
%
%> @brief Calculates the conductance through a graphene antidot in a homogeneous magnetic field.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
function Graphene_antidot_transport( filenum )

if ~exist('filenum', 'var')
    filenum = 1;
end

filename = mfilename('fullpath');
[directory, fncname] = fileparts( filename );

     % The width of the graphene strip
    width = 220;
     % The length of the graphene strip
    height = 350;
     % The center and radius of the antidot
    hole = struct( 'center', [width/2 height/2], 'radius', 90);  
    % Magnetic field in Tesla    
    B = 40;
     % outfilename of the calculated data
     outfilename = [fncname, '_',num2str( filenum )];
     % the output directory
     outputdir   = [];
     % An instance of class antidot
     cAntidot = [];

     % setting the output directory
     setOutputDir()

     % logical value
     transport_calculations_finished = false;


     % creating the scattering disorders in the lattice
     points = round([27 195 27 195]);
     scatterers.zpoints = round(points + rand(size(points))*5);
     scatterers.z = round( ones(size(points))*height + [-30 -30 30 30] + rand(size(points))*5);
     scatterers.siteindexes = [];
     scatterers.aremissing = 1;


     % the calculated spectrum of the lead
     spectrum = []; % calculated spectrum of the
     % cyclotron frequency
     homega = [];
     % The wave function of a localized antidot state
     hole_wavefnc = [];
     % An instance of structure coordinates storing the coordinates of the antidot
     coordinates           = [];
     % edge siteindexes of the antidot
     antidot_edge_points      = [];
     % antidot eigenenergies
     antidot_eigenenergies    = [];

     % creating an array of energies (in eV) to be used in the transport calculations
     Emin = 0.001;
     Emax = 0.22;
     Enum = 200;
     Evec = Emin:(Emax-Emin)/Enum:Emax;
     Enum = length(Evec)-1;

     % An array storing the calculated conductances calculated conductance 
     Conductance = zeros(Enum+1,1); 
     % An array storing the error of the unitarity
     DeltaC      = zeros(Enum+1,1);


     % setting the output directory
     setOutputDir()
     % calculate the conductances
     CalculateTransport();
    
     % plot the calculated data    
     PlotFunction();

    % exporting the calculated data
    save( [outputdir, filesep, [outfilename, '.mat']] );

%% CalculateTransport
%> @brief Calculates the conductance
    function CalculateTransport()
        
          % creating an instance of antidot class
        cAntidot = antidot('width', width, 'height', height, 'B', B, 'hole', hole, 'scatterers', scatterers );
          % creates vector potential handles predefinied in class antidot (This method is obsolete, functions describing the vector potential should be placed inside the present file)
        cAntidot.CreateHandlesForMagneticField();
        
        
        disp(' ')
        disp(' ')
        disp(' ----------------------------------')
        disp('Plotting geometry and spectrum of leads:')
        coordinates = cAntidot.coordinates;
        
        % calculate spectrum of the lead
        cAntidot.setEnergy( Evec(1) )
        cAntidot.CreateRibbon( 'justHamiltonians', true );
        CalcSpectrum();
        cAntidot.getHolePoints();
        antidot_edge_points = cAntidot.antidot_edge_points;
        PlotFunction();        
        
        
        disp(' ')
        disp(' ')
        disp(' ----------------------------------')
        disp('Calculating Conductance for the given magnetic field') 
     
        
        % calculating the conductances as a function of the energy
        for idx = 1:length(Evec)    
            tic
            [C, ~,DeltaC_tmp] = cAntidot.Transport( Evec(idx) );
            toc
            
            Conductance(idx) = C;
            DeltaC(idx) = DeltaC_tmp;
            
            % exporting the calculated data
            save( [outputdir, filesep, [outfilename, '.mat']] );

               % displaying the status of the caclutaions
            disp([ num2str(idx/length(Evec)*100),' % calculated of the conductance.'])
            disp( ' ' )
            
               % plot the results at every energy point
            PlotFunction();
            
            
        end
        
        transport_calculations_finished = true;

          % plot the results
        PlotFunction();
        
    end

%% CalcSpectrum
%> @brief calculates the spectrum of the leads
    function CalcSpectrum()
        
     % Calculate the cyclotron frequency
    rCC = cAntidot.rCC*1e-10; %Angstrom
    phi0 = cAntidot.h/cAntidot.qe; % magnetic flux quantum    
    eta_B = 2*pi/phi0*(rCC)^2*cAntidot.B;    
    homega = 2.97*sqrt(9/2*eta_B); 
    
     %Minimal and maximal value of the selected k points
    kmin = 1*sqrt(3);
    kmax = 1.8*sqrt(3);
     % number of 1D k points 
    knum = 50;
    spectrum = cAntidot.Scatter_UC.CalcSpektrum( 'toPlot', 0, 'ka_min', kmin, 'ka_max', kmax, 'ka_num', knum );
 
    end
 

%% secular_H
%> @brief Hamiltonian for the secular equation
    function H = secular_H(H0,H1,k)
        
        H = H0 + H1*exp(1i*k) + H1'*exp(-1i*k);
        
    end


%% setOutputDir
%> @brief sets output directory
    function setOutputDir()
        resultsdir = 'results';
        mkdir(resultsdir );
        outputdir = [ resultsdir ];        
    end



%% PlotFunction
%> @brief Creates the plot
    function PlotFunction( )
        
        % creating figure in units of pixels
        if ~transport_calculations_finished
            figure1 = figure( 'Units', 'Pixels', 'Visible', 'off' );
        else
            figure1 = figure( 'Units', 'Pixels','Renderer','painters', 'Visible', 'off',...
    'Colormap',[1 1 0.994791686534882;0.954365074634552 0.954365074634552 0.932043671607971;0.908730149269104 0.908730149269104 0.869295656681061;0.863095223903656 0.863095223903656 0.80654764175415;0.817460358142853 0.817460358142853 0.74379962682724;0.771825432777405 0.771825432777405 0.68105161190033;0.726190507411957 0.726190507411957 0.618303596973419;0.680555582046509 0.680555582046509 0.555555582046509;0.643518567085266 0.629629671573639 0.518518567085266;0.606481492519379 0.578703701496124 0.481481492519379;0.569444477558136 0.527777791023254 0.444444477558136;0.532407402992249 0.476851880550385 0.407407432794571;0.495370388031006 0.425925940275192 0.370370388031006;0.458333343267441 0.375 0.333333343267441;0.421296298503876 0.32407408952713 0.296296298503876;0.384259253740311 0.273148149251938 0.259259253740311;0.347222208976746 0.222222223877907 0.222222223877907;0.33983451128006 0.217494085431099 0.217494085431099;0.332446783781052 0.212765961885452 0.212765961885452;0.325059086084366 0.208037823438644 0.208037823438644;0.31767138838768 0.203309699892998 0.203309699892998;0.310283690690994 0.19858156144619 0.19858156144619;0.302895963191986 0.193853422999382 0.193853422999382;0.2955082654953 0.189125299453735 0.189125299453735;0.288120567798615 0.184397161006927 0.184397161006927;0.280732840299606 0.179669037461281 0.179669037461281;0.273345142602921 0.174940899014473 0.174940899014473;0.265957444906235 0.170212760567665 0.170212760567665;0.258569717407227 0.165484637022018 0.165484637022018;0.251182019710541 0.160756498575211 0.160756498575211;0.243794322013855 0.156028375029564 0.156028375029564;0.236406609416008 0.151300236582756 0.151300236582756;0.229018896818161 0.146572098135948 0.146572098135948;0.221631199121475 0.141843974590302 0.141843974590302;0.214243486523628 0.137115836143494 0.137115836143494;0.206855788826942 0.132387712597847 0.132387712597847;0.199468076229095 0.127659574151039 0.127659574151039;0.192080363631248 0.122931443154812 0.122931443154812;0.184692665934563 0.118203312158585 0.118203312158585;0.177304953336716 0.113475181162357 0.113475181162357;0.16991725564003 0.108747042715549 0.108747042715549;0.162529543042183 0.104018911719322 0.104018911719322;0.155141845345497 0.0992907807230949 0.0992907807230949;0.14775413274765 0.0945626497268677 0.0945626497268677;0.140366420149803 0.0898345187306404 0.0898345187306404;0.132978722453117 0.0851063802838326 0.0851063802838326;0.12559100985527 0.0803782492876053 0.0803782492876053;0.118203304708004 0.075650118291378 0.075650118291378;0.110815599560738 0.0709219872951508 0.0709219872951508;0.103427894413471 0.0661938562989235 0.0661938562989235;0.0960401818156242 0.0614657215774059 0.0614657215774059;0.0886524766683578 0.0567375905811787 0.0567375905811787;0.0812647715210915 0.0520094558596611 0.0520094558596611;0.0738770663738251 0.0472813248634338 0.0472813248634338;0.0664893612265587 0.0425531901419163 0.0425531901419163;0.059101652354002 0.037825059145689 0.037825059145689;0.0517139472067356 0.0330969281494617 0.0330969281494617;0.0443262383341789 0.0283687952905893 0.0283687952905893;0.0369385331869125 0.0236406624317169 0.0236406624317169;0.029550826177001 0.0189125295728445 0.0189125295728445;0.0221631191670895 0.0141843976452947 0.0141843976452947;0.0147754130885005 0.00945626478642225 0.00945626478642225;0.00738770654425025 0.00472813239321113 0.00472813239321113;0 0 0]);
        end
        
        % font size on the figure will be 16 points
        fontsize = 12;
              
        
        %************** geometry **********************
        
        % creating axes of the plot        
        deltax_ps = 2;  % 2n-1 es 2n kozotti tavolsag rCC egysegekben
        deltax_pt = 1;    % 2n es 2n+1 kozotti tavolsag rCC egysegekben
        deltax    = 3; % 2n-1 es 2n+1 kozotti tavolsag rCC egysegekben
        deltay = sqrt(3);   
        
        rightside = width/2*deltax*cAntidot.rCC;
        top = height*deltay*cAntidot.rCC;
        
        
        axes_geo = axes('Parent',figure1, ...
                'CLim',[0 0.028],...
                'Visible', 'on',...
                'FontSize', fontsize,...   
                'XTick', [], ...
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        axis equal;
        
        
        
        % plot the data
        if ( ~transport_calculations_finished ) || (usejava('jvm') && ~feature('ShowFigureWindows'))
            for idx = 1:size(antidot_edge_points,1)
             ycoord = deltay*antidot_edge_points(idx,1)*ones(2,1)*cAntidot.rCC;
            
             xcoord = zeros(2,1);
             if mod(antidot_edge_points(idx,2),2) == 0
                xcoord(1) = (antidot_edge_points(idx,2)-2)/2*deltax + deltax_ps;
             else
                xcoord(1) = (antidot_edge_points(idx,2)-1)/2*deltax;
             end
            
             if mod(antidot_edge_points(idx,2),3) == 0
                xcoord(2) = (antidot_edge_points(idx,3)-2)/2*deltax + deltax_ps;
             else
                xcoord(2) = (antidot_edge_points(idx,3)-1)/2*deltax;
             end  
             xcoord = xcoord*cAntidot.rCC;
                   plot(xcoord, ycoord,'+','MarkerSize', 4, 'Parent', axes_geo )  
            end  
 
              
        else
            
            if isempty(hole_wavefnc)
                Energy = 0.1566375;
                hole_wavefnc = cAntidot.CalcWavefnc( Energy, 1, 'toPlot', 0, 'filterHole', 1, 'db', 1, 'infinitemass', 0, 'dotCalc', 0 );
            end            
             delta = 1;
             levels = 25;
            Wavefnc = hole_wavefnc.Wavefnc{1};
             contourf( hole_wavefnc.xcoords(1:delta:end,1:delta:end), hole_wavefnc.ycoords(1:delta:end,1:delta:end), abs(Wavefnc(1:delta:end,1:delta:end)), ...
              levels, 'LineStyle', 'none',  'Parent', axes_geo );  
             text(30,80, ['E = ', num2str(hole_wavefnc.eigenval*1000,3),' meV'], 'color', [0,0,0], 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_geo);
        end
        
        % plot the position of the leads
        plot([0,0], [0,top],'LineWidth', 4, 'color', [0 0 0], 'Parent', axes_geo )
         plot([rightside,rightside], [0,top],'LineWidth', 4, 'color', [0 0 0], 'Parent', axes_geo )  

        % plot scattering centers
        if isfield(scatterers, 'aremissing')
            if scatterers.aremissing
               marker = 'o';
               markersize = 3;
            else
               marker = 'x';
               markersize = 5;
            end
        else
            marker = 'x';
            markersize = 5;
        end 
        plot( scatterers.zpoints/2*deltax*cAntidot.rCC, scatterers.z/2*deltay*cAntidot.rCC, ...
              marker, 'MarkerSize', markersize, 'Parent', axes_geo ) 
         
        

        % Create ylabel
        ylabel('y [A]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_geo);
        
        Lead1_pos = [ mean(coordinates.x(:,1)), min(coordinates.y(1,:))+20];
        text(Lead1_pos(1), Lead1_pos(2), ['Lead'], 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_geo);
        
        Lead2_pos = [ mean(coordinates.x(:,1)), max(coordinates.y(1,:))-20];
        text(Lead2_pos(1), Lead2_pos(2), ['Lead'], 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_geo);
        
        
        %************** sepctrum in the lead  **********************
        x_lim = [min(spectrum(:,1)) max(spectrum(:,1))]/sqrt(3);
        y_lim = [min(spectrum(:,2));max(spectrum(:,2))]*0.8;

        
        
        % creating axes of the plot        
        axes_spectrum = axes('Parent',figure1, ...
                'Visible', 'on',...
                'FontSize', fontsize,...
                'xlim', x_lim,...
                'ylim', y_lim,...  
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        
        Kpoint = 2*pi/3/sqrt(3);
        Kpoint_pos = [Kpoint+0.005, 0.8*y_lim(2)];
        plot(spectrum(:,1)/sqrt(3), spectrum(:,2),'.', 'MarkerSize', 4, 'Parent', axes_spectrum )  
        plot(Kpoint*ones(2,1), [min(spectrum(:,2));max(spectrum(:,2))], 'LineWidth', 1, 'Color', [0 0 0],  'Parent', axes_spectrum )  
        text(Kpoint_pos(1), Kpoint_pos(2), 'K', 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_spectrum);
        
        % Create xlabel
        xlabel('k_yr_{CC}','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_spectrum);
        
        % Create ylabel
        ylabel('E [eV]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_spectrum);
        
        % Create textbox
        Position_Textbox = [0.47 0.51 0.244 0.1274];
        annotation(figure1,'textbox',Position_Textbox,...
            'String',{['$$\hbar\omega = ',num2str(homega,'% 10.3f' ),'$$eV']},...
            'FontSize',fontsize,...
            'FontName','Times New Roman',...
            'interpreter', 'latex', ...
            'LineStyle','none');
        
        
        
        %************** Wavefunction in the lead  **********************
        % wave functionis calculated at a specific k point    
        specific_k = 1.25*sqrt(3);
        
        % create the effective Hamiltonian
        cAntidot.setEnergy(0);
        Surface_temp = cAntidot.Scatter_UC;
        H0 = Surface_temp.Read( 'H0' );
        H1 = Surface_temp.Read( 'H1' );
        H0 = H0 - eye(size(H0))*mean(diag(H0));
        H        = secular_H(H0,H1,specific_k);     
        [waveFnc,E]        = eigs(H,10,0.001);
        E = diag(E);
        chosen_idx =1;
        x_lim = get(axes_geo,'xlim');

        % creating the axis
        axes_wavefnc = axes('Parent',figure1,...
                'Visible', 'on',...
                'FontSize', fontsize,...   
                'xlim', x_lim,...  
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 

        square_Wavefnc = abs(waveFnc(1:length(coordinates.x(:,1)),chosen_idx)).^2;
        plot(coordinates.x(:,1), square_Wavefnc, 'MarkerSize', 4, 'Parent', axes_wavefnc) 
        
        % Create ylabel
        cylabel = ylabel('|\Psi|^2','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_wavefnc);
        
        % Create xlabel
        xlabel('x [A]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_wavefnc);

        chosen_k_pos = [ max(coordinates.x(:,1))*0.05, max(square_Wavefnc)*0.95];
        text(chosen_k_pos(1), chosen_k_pos(2), ['k_yr_{CC} = ',num2str(specific_k/sqrt(3),'% 10.2f')], 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_wavefnc);
        E_pos = [ max(coordinates.x(:,1))*0.05, max(square_Wavefnc)*0.8];
        text(E_pos(1), E_pos(2), ['E = ',num2str(E(chosen_idx),'% 10.2f'),' eV'], 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_wavefnc);
        B_pos = [ max(coordinates.x(:,1))*0.05, max(square_Wavefnc)*0.65];
        text(B_pos(1), B_pos(2), ['B = ',num2str(B,'% 10.2f'),' T'], 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_wavefnc);        
        
        
        
        
        %************** Conductance **********************
        x_lim = [ min(Evec) max(Evec) ];
        axes_cond = axes('Parent',figure1,...
                'Visible', 'on',...
                'FontSize', fontsize,... 
                'xlim', x_lim,...  
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        grid on;
        
        indexek = logical( abs(DeltaC) <= 0.5 );
        
        numerics = plot(Evec(indexek), Conductance(indexek), 'Linewidth', 2, 'color', [0 0 0], 'Parent', axes_cond);
        try
            plot(Evec(indexek), DeltaC(indexek), 'Linewidth', 1, 'color', [0 1 0], 'Parent', axes_cond);
        catch
            disp('No DeltaC data present')
        end
        
        if ~isempty( antidot_eigenenergies )
          y_lim = [0 0.9];%[0 max(Conductance(indexek)) ];
          for idx = 1:length(antidot_eigenenergies)
               plot( ones(1,2)*antidot_eigenenergies(idx), y_lim, 'Linewidth', 1, 'color', [1 0 0], 'Parent', axes_cond);         
          end        
        end
               
         
        
        % Create xlabel
        xlabel('E [eV]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_cond);

        % Create ylabel
        ylabel('\sigma/\sigma_0','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_cond);
                         
        E_pos = [ max(Evec)*0.5, max(Conductance(indexek))*0.35];
        text(E_pos(1), E_pos(2), ['B = ',num2str(B,'% 10.2f'),' T'], 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_cond);
        E_pos = [ max(Evec)*0.5, max(Conductance(indexek))*0.15];
        text(E_pos(1), E_pos(2), ['\phi/\phi_0 = ',num2str(cAntidot.flux_of_hole,3)], 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_cond);
        
        fig_legend = legend(axes_cond, {'conductance', 'unitarity error'});%, 'FontSize', fontsize, 'FontName','Times New Roman')
        set(fig_legend, 'FontSize', fontsize, 'FontName', 'Times New Roman', 'Box', 'off', 'Location', 'NorthEast')
        
        % setting the position and margins of the plot, removing white
        figure_pos = get( figure1, 'Position');
                
            %set the position of the geometry axis
            OuterPosition_geo = get(axes_geo, 'OuterPosition');        
            OuterPosition_geo(1) = 0;
            OuterPosition_geo(2) = figure_pos(4)/2;
            OuterPosition_geo(3) = figure_pos(3)/2;        
            OuterPosition_geo(4) = figure_pos(4)/2;
            set(axes_geo, 'OuterPosition', OuterPosition_geo);
%             Position_geo = get(axes_geo, 'OuterPosition') - get(axes_geo, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1];
%             set(axes_geo, 'Position', Position_geo);
        
            %set the position of the spectrum axis
            OuterPosition_spectrum = get(axes_spectrum, 'OuterPosition');
            OuterPosition_spectrum(1) = figure_pos(3)/2;
            OuterPosition_spectrum(2) = figure_pos(4)/2;
            OuterPosition_spectrum(3) = figure_pos(3)/2;
            OuterPosition_spectrum(4) = figure_pos(4)/2;
            set(axes_spectrum, 'OuterPosition', OuterPosition_spectrum);
%             Position_spectrum = get(axes_spectrum, 'OuterPosition') - get(axes_spectrum, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1];        
%             set(axes_spectrum, 'Position', Position_spectrum);       
        
            %set the position of the wavefnc axis
            OuterPosition_wavefnc = get(axes_wavefnc, 'OuterPosition');
            OuterPosition_wavefnc(1) = 0;
            OuterPosition_wavefnc(2) = 0;
            OuterPosition_wavefnc(3) = figure_pos(3)/2;
            OuterPosition_wavefnc(4) = figure_pos(4)/2;
            set(axes_wavefnc, 'OuterPosition', OuterPosition_wavefnc);  
%         Position_wavefnc = get(axes_wavefnc, 'OuterPosition') - get(axes_wavefnc, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1];
%         Position_wavefnc(1) = Position_geo(1);
%         Position_wavefnc(3) = Position_geo(3);
%         set(axes_wavefnc, 'Position', Position_wavefnc);     
        
            %set the position of the conductance axis
            OuterPosition_cond = get(axes_cond, 'OuterPosition');
            OuterPosition_cond(1) = figure_pos(3)/2;
            OuterPosition_cond(2) = 0;
            OuterPosition_cond(3) = figure_pos(3)/2;
            OuterPosition_cond(4) = figure_pos(4)/2;
            set(axes_cond, 'OuterPosition', OuterPosition_cond);      
%             Position_cond = get(axes_cond, 'OuterPosition') - get(axes_cond, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1];
%             set(axes_cond, 'Position', Position_cond);
        
        
            % setting the position and margins of the plot, removing white spaces
            ver = version('-release');
            if str2num(ver(1:4)) >= 2016
                figure1.PaperPositionMode = 'auto';
                set(figure1, 'PaperSize', [figure_pos(3) figure_pos(4)]);  
            end
        
        % export the figures
        try
            print('-depsc2', fullfile(outputdir, [outfilename,'.eps']))
            print('-dpdf', fullfile(outputdir,[outfilename, '.pdf']))
        catch errCause
            disp('Failed to export figure');
            disp( errCause.identifier );
            disp( errCause.stack(1) );
        end
         close(figure1)
        
    end



end