clear all; close all; clc;
format long;

%constantes
me = 9.10955*10^(-31);
eV = 1.6021895*10^(-19);
hb = 1.0545896*10^(-34);
mp = 1.672651*10^(-27);
a = 1*10^(-10);

N = 100;

%Rendre N pair
%N = 2*abs(N/2);




function [k] = kha(E,m,hb)
  k = sqrt((2*m*E)/hb^2);
endfunction


%Definition de k1
k1 = kha(1*eV,me,hb);



function [g] = G(k,xj,xl)
  im = sqrt(-1);
  den = 2*im*k;
  g = abs(xj-xl);
  g = im*k*g;
  g = exp(g);
  g = g/den;
endfunction

function [u] = U(m,hb,V,a,x)
  k = (2*m)/hb^2;
  if(abs(x)>a/2)
    u=0;
  else
    u = k*V;
  endif
endfunction

function [x] = Xdis(N,d)
  x = linspace(-1*d,d,N);
endfunction

function [x] = Xmis(N,d,p)
  x1 = linspace(-3*p*d,-1*d,N/2);
  x2 = linspace(1*d,p*d,N/2);
  x = [x1 x2];
endfunction

function [phi] = PhiO(k,x)
  im = sqrt(-1);
  phi = exp(-im*k*x);
endfunction


function [] = Proba(E,m,hb,V,a,N,s,p)
k1 = kha(E,m,hb);


% Calcul de B
XB = Xdis(N,a/2);
for i = 1:N
  for j = 1: N
    GBB(i,j) = G(k1, XB(i), XB(j));
    if(i == j)
      UBB(i,j) = U(m,hb,V,a,XB(i));
    else
      UBB(i,j) = 0;
    endif
  endfor
endfor

IB = eye(N);


pasXB = a/N;

M = IB - pasXB*GBB*UBB;
M = inv(M);


for i = 1:N
  PhiOB(i) = PhiO(k1,XB(i));
endfor

PhiB = M*PhiOB';

for i = 1:N
  Y(i) = PhiB(i)*conj(PhiB(i));
endfor


% Calcul de A
XA = Xmis(N,a/2,p);

for i = 1:N
  for j = 1: N
    GAB(i,j) = G(k1, XA(i), XB(j));
  endfor
endfor

for i = 1:N
  PhiOA(i) = PhiO(k1,XA(i));
endfor

PhiA = PhiOA' + pasXB*GAB*UBB*PhiB;

for i = 1:N
  YA(i) = PhiA(i)*conj(PhiA(i));
endfor

%Plot
figure('name', s);
plot(XA(1:N/2),YA(1:N/2),"b",XB,Y,"r",XA(N/2+1:N),YA(N/2+1:N),"b");

%Proba
disp(strcat('Probabilite de passage par effet tunnel,', s, ' :'));
proba = Y(N);
  if(proba > 10^(-2))
    prc = strcat(mat2str(round(proba*100000)/1000) , '%');
  else
    prc = strcat(mat2str(proba*100) , '%');
  endif
disp(prc);


%ouvre un fichier ou le créé
%fid = fopen('test.txt','w');
%fclose(fid);
%fid = fopen('test.txt','at');
%écrit dans ce fichier
%fprintf(fid,'%s\t','x1');
%fprintf(fid,'%s\t','y1');
%fprintf(fid,'%s\t','x1');
%fprintf(fid,'%s\t','y1');
%fprintf(fid,'%s\t','x2');
%fprintf(fid,'%s\t','y2');
%fprintf(fid,'%s\t','x3');
%fprintf(fid,'%s\n','y3');


%for i = 1:N/2  
%  fprintf(fid,'%i',XB(i));
%  fprintf(fid,'\t%i\t',Y(i));
%  fprintf(fid,'%i',XB(i+N/2));
%  fprintf(fid,'\t%i\t',Y(i+N/2));
%  fprintf(fid,'%i',XA(i));
%  fprintf(fid,'\t%i\t',YA(i));
%  fprintf(fid,'%i',XA(i+N/2));
%  fprintf(fid,'\t%i\n',YA(i+N/2));
%endfor

%fclose(fid);



endfunction

Proba(1*eV,me,hb,2*eV,a,N,"E=1eV, V0 = 2eV, Electron",5);

Proba(1*eV,mp,hb,2*eV,a,N,"E=1eV, V0 = 2eV, Proton",0.8);