%%  Device Setup
clear all; clc;
%define constants

%t has been given as a parameter
t = -1;
device = 1;
%input('Valid inputs are 1 or 2 : ');

fprintf(1,'Enter the configuration of device you want to simulate Armchair/Zigzag\n')
fprintf(1,'1 for Armchair\n');
fprintf(1,'2 for Zigzag\n');
type = input('Valid inputs are 1 or 2 : ');
if (isempty(type)||(~((type == 1)||(type == 2))))
    type = 1;
end
%Points in GNR/CNT width
fprintf(1,'Enter the number of Unit Cells in the width\n')
NW = input('Enter a number between 3 and 10 : ');
if isempty(NW)
    NW = 5;
end
%Points in GNR/CNT length
NL = 1;
%Points in GNR/CNT unit cell
NU = 4;
%Define any scatterer as a delta potential barrier
%No scatterer
%Define phase breaking potential
D = 0*1e-12;
%Convergence criteria
errmax = 1e-5;
%Total H will be of LxWxU size
NT = NL*NW*NU;
fprintf(1,'Simulating...\n');

%%  NEGF Setup
%Define Hamiltonians and alpha, beta matrices
fprintf(1,'Hamiltonian Setup...please wait\n');
alphau = t*diag(ones(1,NU-1),1)+t*diag(ones(1,NU-1),-1);
betaw = zeros(NU);
betal = zeros(NU);
%Armchair Structure for GNR
if ((type == 1))
    betaw(2,1) = t; betaw(3,4) = t;
    betal(4,1) = t;
elseif ((type == 2))
    %Zigzag Structure for GNR
    betaw(4,1) = t;
    betal(2,1) = t; betal(3,4) = t;
end

alpha = kron(diag(ones(1,NW)),alphau)+kron(diag(ones(1,NW-1),1),betaw)+kron(diag(ones(1,NW-1),-1),betaw');

beta = kron(diag(ones(1,NW)),betal);
betaTrans = beta';
%Define Energy grid for calculation of Transmission
%E = linspace(0.0,0.0,1);
eta = input('Enter the Energy eta  for contact selfenergy 0.0 < eta < 1E-1 : ');
E = input('Enter the Energy for contact selfenergy -1.0 < E < 1.0 : ');

zplus=1i*eta;
%Define the matrices for NEGF 
% This algorithm solves the equation
% betaTrans = transpose(conjg(beta));
% Imat is an Identity matrix
% zplus = 1i*eta;
% g = inv( (E+1i*eta)*Imat - alpha - betaTrans*g*beta);
% The above equation can alos be written as 
% A*X*B*X + C*X + I = 0
% where the unknown X = g
% A = betaTrans
% B = beta
% C = - ((E+1i*eta)*Imat - alpha)
% I = Imat
% alpha, beta, and Imat are square matrices of size N=NW*NU


T = zeros(1,length(E));
Imat = eye(NW*NU);
g1 = inv((E(1)+zplus)*eye(NW*NU)-alpha);
Fnorm = 100;
ic=0;
g1old = g1;

while(Fnorm>errmax)
    ic=ic+1;
    Ax = (E(1)+zplus)*eye(NW*NU)-alpha-betaTrans*g1*beta;
    g1 = Ax\Imat;
    g1 = (0.75*g1+0.25*g1old);
    Xnorm =  norm(g1-g1old,'fro'); % Difference of new solution and old solution 
    g1old = g1;
    Ax = (E(1)+zplus)*eye(NW*NU)-alpha-betaTrans*g1*beta;
    Qx = Ax*g1 - eye(NW*NU); % Residue
    Fnorm = norm(Qx,'fro');
    fprintf('iterno: %i FnormResidue is %E  FnormDiffX is %E\n ',ic,Fnorm,Xnorm)
end
