function [JCharge,JSpin,JTIS,JFull,JEmpty] = diffuseCurrent( EkT,x )
% it calculates diffuse band current in absence of any standing wave
% electrons. It is simple derivitive of SpinStatistics in respect of spin
% polarization
%
% [JCharge,JSpin,JTIS,JFull,JEmpty] = diffuseCurrent( EkT,x );
% x is spin polarization  EkT  is electron energy /kT
% 
%JCharge  is the charge curent (spin detection)
%JSpin  is the spin curent or current of spin-polarized electrons (spin
%diffusion)
% polarization of the electron gas
% JTIS is the current of spin-unpolarized electrons, which occupies
% half-filled state
% stFull - is the current of spin-unpolarized electrons, which occupies
% full-filled state
%stEmpty is the current of states, which are not filled
%
% use for gragh EkT=-5:0.05:5; EkT=EkT';
%M(:,1)=EkT(:);M(:,2)=JCharge(:);M(:,3)=JSpin(:);M(:,4)=JTIS(:);M(:,5)=JFull(:);
%
%        


JCharge=zeros(length(EkT),length(x));
JSpin=zeros(length(EkT),length(x));
Jsp=zeros(length(EkT),length(x));
JspR=zeros(length(EkT),length(x));
JTIS=zeros(length(EkT),length(x));
JFull=zeros(length(EkT),length(x));
JEmpty=zeros(length(EkT),length(x));

dx=0.0001; % step to take derivation

for jE=1:length(EkT)

    for jx=1:length(x)
        if x(jx)<~0
            [ TIA1, TIS1, stFull1, ~, stEmpty1, ~] = SpinStatistic( EkT(jE),x(jx)+dx);
            [ TIA2, TIS2, stFull2, ~, stEmpty2, ~] = SpinStatistic( EkT(jE),x(jx));
            
            JTIA(jE,jx)=(TIA1-TIA2)/dx;
            JTIS(jE,jx)=(TIS1-TIS2)/dx;
            JFull(jE,jx)=(stFull1-stFull2)/dx;
            JEmpty(jE,jx)=(stEmpty1-stEmpty2)/dx;
        elseif x(jx)>~1
            [ TIA1, TIS1, stFull1, ~, stEmpty1, ~] = SpinStatistic( EkT(jE),x(jx)+dx);
            [ TIA2, TIS2, stFull2, ~, stEmpty2, ~] = SpinStatistic( EkT(jE),x(jx)-dx);
            
            JTIA(jE,jx)=-(TIA1-TIA2)/dx;
            JTIS(jE,jx)=-(TIS1-TIS2)/dx;
            JFull(jE,jx)=-(stFull1-stFull2)/dx;
            JEmpty(jE,jx)=-(stEmpty1-stEmpty2)/dx;

        else
[ TIA1, TIS1, stFull1, ~, stEmpty1, ~] = SpinStatistic( EkT(jE),x(jx)+dx);
[ TIA2, TIS2, stFull2, ~, stEmpty2, ~] = SpinStatistic( EkT(jE),x(jx)-dx);

JTIA(jE,jx)=-(TIA1-TIA2)/2/dx;
JTIS(jE,jx)=-(TIS1-TIS2)/2/dx;
JFull(jE,jx)=-(stFull1-stFull2)/2/dx;
JEmpty(jE,jx)=-(stEmpty1-stEmpty2)/2/dx;

        end
        
JCharge(jE,jx)=JTIA(jE,jx)+JTIS(jE,jx)+JFull(jE,jx);
JSpin(jE,jx)=JTIA(jE,jx);
%Jsp(jE,jx)=JSpin(jE,jx)/JCharge(jE,jx);
%JspR(jE,jx)=Jsp(jE,jx)/x(jx);

    end
end







end