SSH_Josephson.m

Go to the documentation of this file.

%   DC Josephson current through an SSH chain including a topological band inversion -  Based on EQuUs v4.7
%    Copyright (C) 2016 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/.
%%
function SSH_Josephson( filenum )

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




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

outfilename = [fncname,'_',num2str( filenum )];
outputdir   = [];
setOutputDir()

EF = 0.06;
width = 1;
height = 30;
silent = 1;

phi1 = 0;
phi2vec = [];

Edb = 1113;

inputXML1 = fullfile( directory, 'Basic_Input_SSH_Josephson_1.xml');
inputXML2 = fullfile( directory, 'Basic_Input_SSH_Josephson_2.xml');


Ribbons = cell(2,1);
combineribbons = [];

currentvec_continuum = [];
currentvec_continuum_scattering = [];
current_surf = [];
currentvec_ABS = [];
currentvec_NBS = [];
DeltaPhi_vec = [];
currentvec_1range = [];


% combined ribbon
setvectors();
currentVSDeltaPhi();
EgeszAbra();



%% set vectors    
    function setvectors()
        phi_db = 40;
        deltaPhi = pi/(phi_db+1);     
        phi2vec = 0:deltaPhi:pi;      
        %phi2vec = pi/3;        
        
    end
    
%% currentVSDeltaPhi
    function currentVSDeltaPhi()
                
            
    DeltaPhi_vec = phi2vec - phi1;
    
    createRibbons();
    combineribbons = CombineRibbons('Ribbons', Ribbons );
    
    % continuum contribution
    SNSJosephson_handles = SNSJosephson('ribbon', combineribbons, 'gfininvfromHamiltonian', 1);
    if isempty(currentvec_continuum_scattering)
%         [currentvec_continuum_scattering, current_surf] = SNSJosephson_handles.CurrentCalc_continuum( DeltaPhi_vec, 'Edb', Edb );
%         save( [outputdir,'/',outfilename,'.mat'], 'DeltaPhi_vec', 'currentvec_continuum', 'currentvec_continuum_scattering', 'current_surf', 'currentvec_ABS', 'currentvec_NBS', 'currentvec_1range', 'width', 'height', 'EF', 'Edb');
%         EgeszAbra();
    end
    if isempty(currentvec_continuum)
        currentvec_continuum = SNSJosephson_handles.CurrentCalc_discrete( DeltaPhi_vec, 'range', 'CONT', 'Edb', Edb );
        save( [outputdir,'/',outfilename,'.mat'], 'DeltaPhi_vec', 'currentvec_continuum', 'currentvec_continuum_scattering', 'current_surf', 'currentvec_ABS', 'currentvec_NBS', 'currentvec_1range', 'width', 'height', 'EF', 'Edb');
        EgeszAbra();
    end
    
    % contribution from Andreev bound states
    if isempty(currentvec_ABS)
        currentvec_ABS  = SNSJosephson_handles.CurrentCalc_discrete( DeltaPhi_vec, 'range', 'ABS', 'Edb', Edb );
        save( [outputdir,'/',outfilename,'.mat'], 'DeltaPhi_vec', 'currentvec_continuum', 'currentvec_continuum_scattering', 'current_surf', 'currentvec_ABS', 'currentvec_NBS', 'currentvec_1range', 'width', 'height', 'EF', 'Edb');
        EgeszAbra();
    end
    
    % contribution from normal bound states
    if isempty(currentvec_NBS)
        currentvec_NBS  = SNSJosephson_handles.CurrentCalc_discrete( DeltaPhi_vec, 'range', 'NBS', 'Edb', Edb );
        save( [outputdir,'/',outfilename,'.mat'], 'DeltaPhi_vec', 'currentvec_continuum', 'currentvec_continuum_scattering', 'current_surf', 'currentvec_ABS', 'currentvec_NBS', 'currentvec_1range', 'width', 'height', 'EF', 'Edb');
        EgeszAbra();
    end
    
    if isempty(currentvec_1range)
        currentvec_1range = SNSJosephson_handles.CurrentCalc_discrete( DeltaPhi_vec, 'range', 'ALL', 'Edb', Edb );
        save( [outputdir,'/',outfilename,'.mat'], 'DeltaPhi_vec', 'currentvec_continuum', 'currentvec_continuum_scattering', 'current_surf', 'currentvec_ABS', 'currentvec_NBS', 'currentvec_1range', 'width', 'height', 'EF', 'Edb');
    end
    
    
 
    
    end




%% sets the output directory
    function setOutputDir()
        resultsdir = 'results';
        resultsdir = fullfile( directory, resultsdir );
        mkdir(resultsdir );
        outputdir = resultsdir;        
    end


%% createRibbons
% Fills ribbon list
function createRibbons()
     Ribbons{1} = Ribbon('E', 0, 'EF', EF, 'width', width, 'height', height, 'silent', silent, 'phi', [phi1, phi1+DeltaPhi_vec(1)],  'filenameIn', inputXML1, 'filenameOut', fullfile( outputdir, [outfilename,'_1.xml']) );
     Ribbons{2} = Ribbon('E', 0, 'EF', EF, 'width', width, 'height', height, 'silent', silent, 'phi', [phi1, phi1+DeltaPhi_vec(1)],  'filenameIn', inputXML2, 'filenameOut', fullfile( outputdir, [outfilename,'_2.xml']) );
end






%% plotfunction
    function EgeszAbra()
        
        figure1 = figure();
        
        fontsize = 16;
          
        
                
        Position = [0.18 0.18 0.53 0.53];
        axes1 = axes('Parent',figure1, 'Position', Position,...
                'Visible', 'on',...
                'FontSize', fontsize,...                 'xlim', x_lim,...                    'ylim', y_lim,...                'XTick', XTick,...                'YTick', YTick,...
                'Box', 'on',...
                'FontName','Times New Roman');
        hold on; 
        
        % Create xlabel
        xlabel('\Delta\phi/\pi','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes1);

        % Create ylabel
        ylabel('I','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes1);
        
        legend_labels = [];
        
        if ~isempty( currentvec_continuum )
            indexes = logical( currentvec_continuum );
            plot(DeltaPhi_vec(indexes)/pi, currentvec_continuum(indexes), 'b-', 'Parent', axes1)
               legend_labels{end+1} = 'CONT';
        end

        if ~isempty( currentvec_continuum_scattering )
            indexes = logical( currentvec_continuum_scattering );
            plot(DeltaPhi_vec(indexes)/pi, currentvec_continuum_scattering(indexes), 'b--', 'Parent', axes1)
               legend_labels{end+1} = 'CONT S-matrix';
        end
        
        if ~isempty( currentvec_ABS )
            indexes = ~isnan( currentvec_ABS );
            plot(DeltaPhi_vec(indexes)/pi, currentvec_ABS(indexes), 'Linestyle', '-', 'Color', [0 0.83 0], 'Parent', axes1) 
               legend_labels{end+1} = 'Andreev bound states';
        end
        
        if ~isempty( currentvec_NBS )
            indexes = ~isnan( currentvec_NBS );
            plot(DeltaPhi_vec(indexes)/pi, currentvec_NBS(indexes), 'r-', 'Parent', axes1)
               legend_labels{end+1} = 'Normal bound states';
        end
        
        if (~isempty( currentvec_NBS )) && (~isempty( currentvec_continuum )) && (~isempty( currentvec_NBS ))
            total = currentvec_continuum+currentvec_ABS+currentvec_NBS;
            indexes = ~isnan( total );
            plot(DeltaPhi_vec(indexes)/pi, total, 'k-', 'Parent', axes1)
               legend_labels{end+1} = 'Total current';
        end          
        
        
        if ~isempty( currentvec_1range )
            indexes = ~isnan( currentvec_1range );
            plot(DeltaPhi_vec(indexes)/pi, currentvec_1range(indexes), 'c-', 'Parent', axes1)
               legend_labels{end+1} = 'ALL';
        end
        
               

        fig_legend = legend(axes1, legend_labels);
        set(fig_legend, 'FontSize', fontsize, 'FontName','Times New Roman', 'Box', 'off', 'Location', 'NorthEast')
    
    
        xlabel('\Delta\phi/\pi')
        ylabel('I')        
               

        print('-depsc2', fullfile(outputdir, [outfilename,'.eps']))
        close(figure1)
        
        
        
        
    end

end