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);

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

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

function [a] = alph(E,m,hb,V)
  a = sqrt((2*m*(V-E))/hb^2);
endfunction

%rho(3*eV,me,hb)

%kha(3*eV,me,hb,2*eV)

%alph(2*eV,me,hb,4*eV)

function [r] = Ref(E,m,hb,V,a)
  if(E == V)
    p = rho(E,m,hb);
    r = p*p*a*a/4;
    r = r/(1+r);
  else
    p = rho(E,m,hb);
    k = kha(E,m,hb,V);
    r = (p^2-k^2)/(2*p*k);
    r = r*r;
    r = (r*(sin(k*a))^2)/(1+r*(sin(k*a))^2);
  endif
endfunction

function [t] = Tran(E,m,hb,V,a)
  if(E == V)
    p = rho(E,m,hb);
    t = p*p*a*a/4;
    t = 1/(1+t);
  else
    p = rho(E,m,hb);
    k = kha(E,m,hb,V);
    t = (p^2-k^2)/(2*p*k);
    t = t*t;
    t = 1/(1+t*(sin(k*a))^2);
  endif
endfunction

function [r] = Refc(E,m,hb,V,a)
  im = sqrt(-1);
  if(E == 0)
    r = 1;
  elseif(E == V)
    p = rho(E,m,hb);
    r = p*p*a*a/4;
    r = r/(1+r);
  else
    p = rho(E,m,hb);
    k = alph(E,m,hb,V);
    r = (p^2-(im*k)^2)/(2*p*(im*k));
    r = r*r;
    r = (r*(im*sinh(k*a))^2)/(1+r*(im*sinh(k*a))^2);
  endif
endfunction

function [t] = Tranc(E,m,hb,V,a)
  im = sqrt(-1);
  if(E == 0)
    t = 0;
  elseif(E == V)
    p = rho(E,m,hb);
    t = p*p*a*a/4;
    t = 1/(1+t);
  else
    p = rho(E,m,hb);
    k = alph(E,m,hb,V);
    t = (p^2-(im*k)^2)/(2*p*(im*k));
    t = t*t;
    t = 1/(1+t*(im*sinh(k*a))^2);
  endif
endfunction


N = 200;

V0 = 2*eV;
E = linspace(0*V0,10*V0,N);
Rc = zeros(N);
Tc = zeros(N);

for i = 1:N
  if(E(i) < V0)
    Rc(i) = Refc(E(i),me,hb,V0,a);
    %Tc(i) = Tranc(E(i),me,hb,V0,a);
  else
    Rc(i) = Ref(E(i),me,hb,V0,a);
    %Tc(i) = Tran(E(i),me,hb,V0,a);
  endif
  Tc(i) = 1 - Rc(i);
endfor

Ex = E/eV;

l1 = length(Ex);
l2 = length(Rc);


%coefficients R et T
%figure;
%plot(Ex,Rc,Ex,Tc)

% Fonctions de calcul des coefficients
function [aa] = AA(E,m,hb,V,a)
  p = rho(E,m,hb);
  k = kha(E,m,hb,V);
  aa = i*exp(-i*p*a);
  aa = aa * (k*k - p*p)*sin(k*a);
  den = 2*p*k*cos(k*a) - i*(p*p + k*k)*sin(k*a);
  aa = aa/den;
endfunction

function [ba] = BA(E,m,hb,V,a,ind)
  p = rho(E,m,hb);
  k = kha(E,m,hb,V);
  ba = ind*p*(p+ind*k)*exp(-ind*i*a*(k+ind*p)/2);
  den = 2*p*k*cos(k*a)-i*(p*p+k*k)*sin(k*a);
  ba = ba/den;
endfunction

function [ca] = CA(E,m,hb,V,a)
  p = rho(E,m,hb);
  k = kha(E,m,hb,V);
  ca = 2*p*k*exp(-i*p*a);
  ca = ca * (2*p*k*cos(k*a) + i*(k*k+ p*p)*sin(k*a));
  den = 4*p*p*k*k*cos(k*a)*cos(k*a)+(p*p+k*k)^2 * sin(k*a)*sin(k*a);
  ca = ca/den;
endfunction

function [ca] = Proba(E,m,hb,V,a,n)
%définition des constantes
p3 = rho(E,m,hb);
k3 = kha(E,m,hb,V);
%calcul des coefficients
AAmp = AA(E,m,hb,V,a);
BAp = BA(E,m,hb,V,a,1);
BAm = BA(E,m,hb,V,a,-1);
CAp = CA(E,m,hb,V,a);


%Définition de l'intervalle
X1x = linspace(-5*a,-a/2,n);
X2x = linspace(-a/2,a/2,n);
X3x = linspace(a/2,5*a,n);
Xfx = [X1x X2x X3x];

for i = 1:n
  im = sqrt(-1);
  Y1y(i) = AAmp*exp(-1*im*p3*X1x(i)) + exp(im*p3*X1x(i));
  Y1y(i) = Y1y(i)*conj(Y1y(i));
  Yfy(i) = Y1y(i);
endfor


for i = 1:n
  im = sqrt(-1);
  Y2y(i) = BAm*exp(-1*im*k3*X2x(i)) + BAp*exp(im*k3*X2x(i));
  Y2y(i) = Y2y(i)*conj(Y2y(i));
  Yfy(i+n) = Y2y(i);
endfor


for i = 1:n
  Y3y(i) = CAp * conj(CAp);
  Yfy(i+2*n) = Y3y(i);
endfor

%norme de Phi au carré
figure;
plot(X1x,Y1y,"b",X2x,Y2y,"r",X3x,Y3y,"b");

endfunction

%cas E=3eV
Proba(3*eV,me,hb,2*eV,a,N);

%cas E=1eV
Proba(1*eV,me,hb,2*eV,a,N);


%cas E=10eV
Proba(20*eV,me,hb,2*eV,a,N);