Aide Modélisation-> Diffusion quantique par un potentiel unidimensionnel

Version imprimable

2 pièce(s) jointe(s)

Bonjour à tous.
Si je viens vous demander votre aide c'est que là je suis vraiment bloqué... Ca fait une bonne dizaine d'heures que je travail sur ce projet et ça avançait plutôt bien jusqu'à présent mais là je me retrouve dans une impasse.
Que je vous explique mon problème: Dans mon UE de méthodes numériques il m'est demandé de réaliser la modélisation de la fonction d'onde d'une particule à travers un potentiel (énoncé en annexe.).
La partie interessante de l'énoncé est la suivante:
Pour résoudre le problème numériquement, il est moins couteux en temps de calcul de prendre le problème à
l’envers :
– on prend en x > L une fonction d’onde eikx ;
– on écrit les conditions aux limites en L ;
– on résout l’équation différentielle entre L et 0 ;
– les conditions aux limites en 0 permettent de calculer la solution en x < 0 : Aeikx + Be−ikx ;
– enfin,A=1/T etB=R/T.
L'équation de Schrodinger dans un état stationnaire à une dimension est donné par l'équation différentielle de type y"=-a(x)y.
J'ai donc créé un code me permettent de résoudre cette équation avec a(x) variant suivant si on se trouve en dedans ou en dehors du potentiel.
Et....ca marche! (enfin ça m'a l'air). Petit graphe sous root en annexe.
Mais là où je bloque c'est pour déterminer le coefficient de transmission et le coefficient de réflexion. Suivant les conditions limite en x=0 on a psi(0)=A+B et psi(0)'=k(A-B). Donc on en tire A et B et donc R et T. Mais là je trouve des trucs complètement aberrants! Sachant que R+T=1 normalement. (mon code en annexe)
Est-ce que le ciel (ou vous) peut me venir en aide?! Merci:zoubi:.
PS: n'hésitez pas si vous avez des questions, j'ai conscience de ne pas avoir été forcément très clair sur tout les points^^

Cordialement, Jeremy.

Annexe:
Enoncé:
Pièce jointe 208779

Graphe:
Pièce jointe 208780

Code:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
#ifndef __EDO2_H__
#define __EDO2_H__
#include "Vecteur.h"
#include <cmath>
class EDO2 {
protected:
    int N; double E;
    Vecteur yf; // Valeurs de y(xf)  stockees ici par rk4()
    double *xtab;
    Vecteur *ytab, *yptab;
    void stocke(int i,double x, Vecteur y, Vecteur yp);
 
public:
    double  L=21e-9;
    double V=1.6e-19,m=1.9e-31,hb=1.04082e-33, k1=sqrt(2*m*E)/hb, k2=sqrt(2*m*(V-E))/hb;
 
    EDO2(int Ni, long double Ei);
    ~EDO2();
    Vecteur evalyp(double x, Vecteur y);
 
    // integration par Runge-Kutta d'ordre 4
    Vecteur rk4(double x0,Vecteur y0, double xf);
 
    void ecrirexyz(const char * nom_fichier);
    void ecrirexyzROOT(const char * nom_fichier);
};
 
#endif

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
#include "EDO2.h"
#include <fstream>
#include <iomanip>
#include <iostream>
using namespace std;
#include "TGraph.h"
#include "TFile.h"
#include "TString.h"
 
// On veux résoudre y"=-ay, on pose y1=y et y2=y' => y1'=y2 y2'=-ay1
 
EDO2::EDO2(int Ni, long double Ei) : N(Ni), E(Ei)
{
    yf.redim(2);
    // On alloue les tableaux servant a stocker x, y, y', z, z' pour chaque point.
    xtab=new double [N+1]; // [0..N] => N+1
    ytab=new Vecteur [N+1];
    yptab=new Vecteur [N+1];
    for (int i=0;i<N+1;i++) {
        ytab[i].redim(2);
        yptab[i].redim(2);
    }
}
EDO2::~EDO2()
{
    delete[] xtab; delete[] ytab; delete[] yptab;
}
 
void EDO2::stocke(int i,double x, Vecteur y, Vecteur yp)
{
    xtab[i]=x;
    ytab[i]=y;
    yptab[i]=yp;
}
void EDO2::ecrirexyz(const char * nom_fichier)
{
    ofstream fichier( nom_fichier );
    // On ecrit tout en 5 colonnes dans un fichier
    for (int i = 0 ; i < N+1 ; i++) {
        fichier << setw(20) << setprecision(12) << xtab[i];
        for ( int d=1; d<=1; d++){
            fichier << setw(20) << setprecision(12) << ytab[i](d)<< setw(20) << setprecision(12) << yptab[i](d);
        }
        fichier << endl;}
    fichier.close();
}
void EDO2::ecrirexyzROOT(const char * nom_fichier)
{
    TFile fichier( nom_fichier , "RECREATE" );
    TGraph* gr=new TGraph [2]; TGraph* grp=new TGraph [2]; for (int d=0; d<1; d++)
    {
        TString s("Fonction "); TString name("F");
        gr[d].Set(N+1); // mets le nombre de points dans les graphes s+=(d+1);name+=(d+1);
        gr[d].SetNameTitle(name,s);
        grp[d].Set(N+1); // mets le nombre de points dans les graphes des derivees
        s+=" derivee"; name+="p";
        grp[d].SetNameTitle(name,s);
        for (int i = 0 ; i < N+1 ; i++)
        {
            gr[d].SetPoint(i,xtab[i],ytab[i](d+1));
            grp[d].SetPoint(i,xtab[i],yptab[i](d+1));
        }
        gr[d].Write();
        grp[d].Write();
    }
    delete [] gr; delete [] grp;
    fichier.Close();
}
 
Vecteur EDO2::evalyp(double x, Vecteur y)
{
    long double K, Ki=-((V-E)*2*m)/(hb*hb), Ke=-((-E)*2*m)/(hb*hb); //a=2m/h
    if (x>=0 && x<=L) K=Ki;
    else K=Ke;
    Vecteur v(2);
    v(1)=y(2);
    v(2)=-K*y(1);
    return v;
}
 
Vecteur EDO2::rk4(double x0, Vecteur y0,double xf)
{double x;
    Vecteur y(2);
    Vecteur yp(2);
    Vecteur k1(2), k2(2), k3(2), k4(2);
    y=y0;
 
    double h = (xf - x0)/N;
    for (int i = 0 ; i<N ; i++) {
        x = x0 + i*h;
        yp= evalyp(x, y);
        stocke(i,x,y,yp);
        k1 = h * yp;
        k2 = h * evalyp(x+h/2, y+k1/2 );
        k3 = h * evalyp(x+h/2, y+k2/2);
        k4 = h * evalyp(x+h,   y+k3 );
        y=y+(k1+2*k2+2*k3+k4)/6;
            }
    yf=y;
    stocke(N,xf,y,evalyp(xf,y));
    return yf; }

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include "EDO2.h"
using namespace std;
#include "TGraph.h"
#include "TFile.h"
#include <cmath>
 
int main() {
    cout << endl << "EDO2" << endl;
    cout << endl << "N=1000" << endl;
    EDO2 T1(1000, 1.1e-19);
 
    Vecteur yres(2);
    double y00[2]={0,cos(T1.k1*100e-9)};
    Vecteur y0(2,y00);
 
    yres=T1.rk4(100e-9,y0,-100e-9);
 
    T1.ecrirexyz("Diffusion.dat");
    T1.ecrirexyzROOT("Diffusion.root");
 
 
    //k1=k3=sqrt((2*m*E)/(h*h))
    //k2=sqrt((2*m*(E-V))/(h*h))
 
    // 1/T + R/T = phi0
    // (1/T - R/T)k1 = phi0'
 
    double psi0, psi0p;
    psi0=T1.rk4(100e-9,y0,0)(1);
    psi0p=T1.rk4(100e-9,y0,0)(2);
    cout<<"psi0="<<psi0<<endl;
    cout<<"psi0p="<<psi0p<<endl;
 
    double A,B,R,T;
    A=0.5*(psi0+(psi0p/T1.k1));
    B=psi0-A;
    cout<<"A="<<A<<endl;
    cout<<"B="<<B<<endl;
 
    T=1/A;
    R=B*T;
    cout<<"T="<<T<<endl;
    cout<<"R="<<R<<endl;
    cout<<"R+T="<<R+T<<endl;
    cout<<"k(A-B)="<<T1.k1*((1/T)-(R/T))<<endl;
    return 0;}

05/05/2016, 02h42
koala01
Salut,

En toute honnêteté, la diffusion quantique et moi, ca fait deux... Je ne vais donc pas m'étendre sur l'algorithme, car je suis totalement incapable de le corriger ;)

Mais un rapide survol de ton code m'a fait littéralement bondir sur ma chaise pour une raison toute simple : tu code exactement comme si tu codais en C (hormis le fait que tu utilises new et non malloc). Et ca, ca pourrait te mériter des baffes!

Un exemple frappant :
Code:

1 2 3 xtab=new double [N+1]; // [0..N] => N+1 ytab=new Vecteur [N+1]; yptab=new Vecteur [N+1];
(et encore, je te fais grace de la fonction redim() que je n'ai pas vérifiée)
Sais tu que tu donne trois occasions à ton code de lever une exception en faisant cela :question:
Sais tu que la présence des pointeurs xtab, ytab et yptab tend à transformer ta classe en bombe à retardement si tu essaye de la copier ou d'utiliser l'opérateur = :question:

Si au moins tu utilisais la classe std::vector (car je présumes que la taille de ces trois tableaux est définie à l'utilisation), tu pourrais surement gagner énormément en "fiabilité" de la part de ta classe...

Ta fonction ecrirexyz :
Evite de te taper sur la tête avec un marteau pour dire que ca "fait tellement de bien quand ca s'arrête" et transmet le nom du fichier sous la forme d'une référence constante vers un objet de type std::string... (d'ailleurs, plus simplement, utilise std::string à chaque fois que tu devras manipuler des chaines de caractères : tu ne t'en sentiras que mieux ;) )
l'ouverture d'un fichier peut échouer. ta variable fichier serait alors égale à false, si cela devait arriver... Tu serais bien inspiré de tester l'ouverture du fichier sous la forme de
Code:

1 2 3 4 5 6 7 std::ofstream fichier(nom); if(fichier){ // OK, il est ouvert, on peut faire ce qu'on veut ici }else{ // oopps, pourquoi il n'est pas ouvert??? }
De plus, le destructeur de std::ofstream va automatiquement fermer le fichier lorsqu'il sera appelé... y a donc absolument aucune raison de le fermer explicitement (la ligne fichier.close(); ne sert absolument à rien ;) )

la fonction ecrirexyzROOT :
pourquoi utiliser l'allocation dynamique de la mémoire pour gr et grp alors que leur taille est clairement définie :question: Pourquoi ne pas avoir recours à la classe (proposée par la bibliothèque standard depuis la sortie de C++11) qui fonctionne comme un tableau de taille connue à la compilation : std::array :question:
Là, tu donnes deux occasions à ton code de lancer une exception, alors que tu pourrais carrément éviter tout risque en n'ayant simplement pas recours à l'allocation dynamique de la mémoire ;)
N, si je ne m'abuse, représente la taille de tes tableaux, ce qui implique qu'un accès correct ne devrait aller que de 0 à N inclus... Pourquoi utilises tu N+1 à la ligne gr[d].Set(N+1); // mets le nombre de points dans les graphes s+=(d+1);name+=(d+1); (meme chose pour la ligne grp[d].Set(N+1); // mets le nombre de points dans les graphes des derivees)

Il se peut que je me trompe (je l'ai dit, je n'ai fait que parcourir ton code en diagonale), mais je ne serais pas étonné si le fait d'utiliser N+1 au lieu de N ne venait à occasionner un accès "hors bornes", avec pour conséquences le fait que tu essayerais d'aller chercher un élément de plus que ce que tes tableaux ne contiennent. Du coup, le dernier élément aurait très vraisemblablement une valeur totalement aberrante correspondant aux "résidus d'une utilisation antérieure de la mémoire", et c'est, j'en suis presque sur, ce qui tes résultats aberrants ;)
Tu utilises std::ofstream dans ecrirexyz et un TFile dans ecrirexyzROOT, c'est normal :question: D'ailleurs, c'est quoi cette classe TFile :question:

la fonction evalyp:
Heu, tu pourrais essayer de l'indenter un tout petit peu correctement ? Car je me demande si tu sais ce qui se passe réellement dans cette fonction ou pas ;)

Et, de manière générale:
Je ne te dis pas tout le temps qu'il m'a fallu pour me rendre compte qu'il n'y avait pas deux, mais trois variables déclarées dans la ligne de code long double K, Ki=-((V-E)*2*m)/(hb*hb), Ke=-((-E)*2*m)/(hb*hb); //a=2m/h :P
les quatre qualités essentielles d'un code sont :
1. d'être facilie à lire
2. d'être facile à écrire
3. d'être facile à comprendre
4. d'être facile à faire évoluer
Cela ne fera aucune différence au niveau de l'exécutable généré si tu prend la peine d'écrire une seule (et unique!) instruction par ligne. Et cela ne changera rien non plus si tu décide d'utiliser des noms clairement explicites.

Peut être que des noms comme xtab, ytab yptab, gr ou grp font partie du langage spécifique à ton domaine, et, si tu as la certitude que n'importe quel collègue arrivera à les comprendre sans devoir venir te demander à quoi cela correspond, tant mieux pour toi, mais pour ma part, ces noms ne me disent absolument rien. Et comme tu n'es absolument pas limité par le nombre de caractères de tes identifiants, si tu pouvais trouver des nom qui indiquent clairement à quoi servent ces variables, tout le monde t'en serait sans doute reconnaissant ;)

Salut,
Oui, je suis désolé pour ce genre de désagrément, ça ne fait que 7-8 mois que je code. xtab, ytab, yptab étant définit comme protected dans ma classe il n'y a pas de problème non? C'est marrant que tu me dise que les fonctions ecrirexyz et ecrirexyzROOT soit écrits avec les pieds car ce n'est pas moi qui les ai écrites mais mes profs hahaha! Voyant leurs utilité et un peu près leurs fonctionnement j'ai fais un bête copier-coller d'est corrigé de TP.

la fonction evalyp:
plus propre ;) :
Code:

1 2 3 4 5 6 7 8 9 10 Vecteur EDO2::evalyp(double x, Vecteur y) {long double K=-(2*m)/(hb*hb), Kiterieur=K*(V-E), Kexterieur=K*-E; if (x>=0 && x<=L) K=Kiterieur; else K=Kexterieur; Vecteur v(2); v(1)=y(2); v(2)=-K*y(1); return v;}
Pour résoudre une équation différentielle d'ordre N il faut créer un système à N équations différentielle du premier ordre.
D'ou y"=-a(x)y => changement de variable on pose; => v=y' => v'=y"
On le définit comme un vecteur v(1)=y' et v(2)=y"
On définit notre fonction elle aussi comme un vecteur=> y(1)=y et y(2)=y'
On applique ça à notre système v(1)=y(2) et v(2)=-a(x)y(1)
Or a(x) dépends de si on se trouve dans le potentiel ou non, d'où mon if(x>...

Citation:
Envoyé par momovelo
Salut,
la fonction evalyp:
plus propre ;) :
Code:

1 2 3 4 5 6 7 8 9 10 Vecteur EDO2::evalyp(double x, Vecteur y) {long double K=-(2*m)/(hb*hb), Kiterieur=K*(V-E), Kexterieur=K*-E; if (x>=0 && x<=L) K=Kiterieur; else K=Kexterieur; Vecteur v(2); v(1)=y(2); v(2)=-K*y(1); return v;}
la fonction evalyp plus propre :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Vecteur EDO2::evalyp(double x, Vecteur y) { long double K = -(2 * m) / (hb * hb); long double Kiterieur = K * (V - E); long double Kexterieur = K * -E; if (x>=0 && x<=L) { K=Kiterieur; } else { K=Kexterieur; } Vecteur v(2); v(1) = y(2); v(2) = -K * y(1); return v; }
(Je n'ai pas renommé les variables, ne sachant pas à quoi elles correspondent...)

1) Moi, j'aime bien quand mon code respire.
2) Un if sur une ligne, c'est mal, un if/else/for/while... sans accolades, c'est error prone.
3) Les déclarations multiples séparées par une virgule c'est à éviter, si on veut comprendre le code.

05/05/2016, 10h16
Bousk

Attention tu as changé les types de Kiterieur et Kexterieur de long double en long
05/05/2016, 10h26
dragonjoker59

Bien vu, j'ai corrigé (Ca ne fait que confirmer ce que je pense, au final... J'avais pas vu le "double")
05/05/2016, 12h07
momovelo

Personnellement j'aime bien quand c'est clair et aéré pour des parties de code dont je ne suis pas sûr. Au contraire, pour les parties dont je suis sûr j'aime bien avoir quelque chose de très compacte qui prends le moins de lignes possible.

Enfin bref, ça ne m'aide pas beaucoup tout ça :mrgreen:
05/05/2016, 19h08
koala01

C'est une très mauvaise habitude...

La raison est simple : un code est beaucoup plus souvent lu qu'il n'est écrit, et il y a de très fortes chances pour que tu doive relire ton code TRES longtemps après avoir eu largement le temps d'oublier ce qu'il est sensé faire.

Aère ton code, fais en sorte qu'il soit le plus facilement compréhensible par TOUT LE MONDE, que l'on ne passe pas un "temps infernal" à se poser des questions idiotes (du genre de : où telle variable est-elle définie ? combien de variable il me déclare sur telle ligne ? ce if, il englobe quelle instruction ? et ce else, alors ?)

Déjà, cela simplifiera la vie de ceux qui voudront bien t'aider, mais, en plus, cela simplifiera ta propre vie, car tu peux me croire : lorsque tu reviendras sur le code après parfois trois ou six mois (ou plus), tu seras exactement dans la même situation que quelqu'un qui n'a pas la moindre idée de ce que peut faire ton code ;)
06/05/2016, 12h02
Prof

Bonjour.

Il y a quelque chose qui me chifonne dans ce programme.

Les fonctions d'ondes sont des fonctions complexes.
D'ailleurs, $Formule mathématique$ est visiblement un complexe.

Par contre, dans votre code, vous traitez la fonction inconnue y comme une fonction réelle.

Le i des complexes a disparu de vos calculs, ce qui vous amène à écrire

Citation:

psi(0)'=k(A-B)

alors que cela devrait être : psi(0)'=ik(A-B).

Et plus loin, dans le test, vous introduisez une onde monochromatique plane sous la forme d'une fonction réelle du type y = cos(kx).

Le problème, c'est qu'une telle fonction ne vérifie pas y'= ky.
06/05/2016, 15h53
momovelo

C'est exactement ce que je me suis dit cette nuit. J'ai donc essayé de traité d'un côté la partie réelle et de l'autre la partie imaginaire. Mais ça ne fonctionne toujours pas. Le plus simple serait de directement travailler avec les complexes. Je vais voir s'il existe des librairies gsl qui me le permettent.
06/05/2016, 16h16
Bousk

Je sais pas si ça répond au sujet, mais il existe std::complex

Malheureusement ce n'est pas si simple...
Si j'initialise mon y0 comme ceci (y0 étant la première valeur utilisé pour ma méthode rk4):
Code:

1 2 3 4 5 complex<double>monComplex(cos(T0.k1*x0), sin(T0.k1*x0 )); complex<double>monComplexP(-T0.k1*sin(T0.k1*x0), T0.k1*cos(T0.k1*x0 )); complex<double> y00[2]={monComplex,monComplexP}; Vecteur y0(2,y00);
je me retrouve avec ce genre d'erreur:

Citation:

main.cpp:27:13: error: no matching constructor for initialization of 'Vecteur'
Vecteur y0(2,y00);
^ ~~~~~
./Vecteur.h:21:3: note: candidate constructor not viable: no known conversion
from 'complex<double> [2]' to 'const double *' for 2nd argument
Vecteur(int dim_i, const double * Vtab_i); // Prend un tableau en val. init.
^
./Vecteur.h:22:3: note: candidate constructor not viable: no known conversion
from 'complex<double> [2]' to 'const gsl_vector *' for 2nd argument
Vecteur(int dim_i, const gsl_vector * Vtab_i); // Prend un gsl_vector ...
^
./Vecteur.h:19:3: note: candidate constructor not viable: allows at most single
argument 'dim_i', but 2 arguments were provided
Vecteur(int dim_i = 1); // Initialise les composantes à 0
^
./Vecteur.h:23:3: note: candidate constructor not viable: requires single
argument 'v', but 2 arguments were provided
Vecteur(const Vecteur & v); // Constructeur par copie
^

06/05/2016, 17h35
Bousk

Je ne vois là que des erreurs de compilation dû à l'introduction d'un nouveau type qu'il faut prendre en compte. :calim2:

Au fait, j'ai réussi ^^
Comme vous vous en doutiez (moi aussi) le problème était d'ordre mathématique. En fait les coefficient A et B sont complexe!
Si on reprend nos conditions au limites:
psi(0)=Aexp(ik*0)+Bexp(-ik*0) => psi(0)=A+B
psi'(0)= ikAexp(ik*0)-ikBexp(-ik*0) =>psi'(0)=ik(A-B)

Maintenant on pose A=a+ia' et B= b+ib'

psi(0)= a+b +i(a'+b') => d'où Re.psi(0)=a+b et Im.psi(0)a'+b'
On déroule on retrouve A et B et enfin T et R!
Voici mon code si ça vous intéresse:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
Vecteur EDO2::rk4(double x0,double xf, double E)
{
    double k=(sqrt(2*m*E))/hb;
    double y00Re[2]={cos(k*x0),-k*sin(k*x0)};
    Vecteur y0Re(2,y00Re); //condition initiale Re
 
    double y00Im[2]={sin(k*x0),k*cos(k*x0)};
    Vecteur y0Im(2,y00Im); //condition initiale Im
 
 
    double x;
    Vecteur yRe(2),yIm(2);
    Vecteur ypRe(2),ypIm(2);
    Vecteur k1(2), k2(2), k3(2), k4(2);
    yRe=y0Re;
    yIm=y0Im;
    double ynorm;
    double h = (xf - x0)/N;
    for (int i = 0 ; i<N ; i++) {
        x = x0 + i*h;
        ynorm=yRe(1)*yRe(1)+yIm(1)*yIm(1);
        ypRe= evalyp(x,E, yRe);
        ypIm= evalyp(x,E, yIm);
        stocke(i,x*(1/nm),ynorm,yRe,ypRe,yIm,ypIm);
 
       //Re
        k1 = h * ypRe;
        k2 = h * evalyp(x+h/2,E ,yRe+k1/2 );
        k3 = h * evalyp(x+h/2,E, yRe+k2/2);
        k4 = h * evalyp(x+h,E ,  yRe+k3 );
        yRe=yRe+(k1+2*k2+2*k3+k4)/6;
 
        //Im
        k1 = h * ypIm;
        k2 = h * evalyp(x+h/2,E, yIm+k1/2 );
        k3 = h * evalyp(x+h/2,E, yIm+k2/2);
        k4 = h * evalyp(x+h, E,  yIm+k3 );
        yIm=yIm+(k1+2*k2+2*k3+k4)/6;
 
    yf(1)=yRe(1);
    yf(2)=yRe(2);
    yf(3)=yIm(1);
    yf(4)=yIm(2);
 
    }
    stocke(N,xf*(1/nm),ynorm,yRe,evalyp(xf,E,yRe),yIm,evalyp(xf,E,yIm));
    return yf; }
 
void EDO2::CalculT(double x0, double E0,double Ef)//Calcul de T pour différentes energies
{double h = (Ef - E0)/N;
    double E,k, a,b,ia,ib,t,it,r,ir,T,TAnaly,K;
    Vecteur psi0;
    for(int i=0; i<=N; i++)
    {E = E0 + i*h;
        psi0=rk4(x0,0,E);
 
        k=(sqrt(2*m*E))/hb;//issu de la physique quantique
        ib=0.5*(psi0(3)+(psi0(2)/k));
        ia=0.5*(psi0(3)-(psi0(2)/k));
        a=0.5*(psi0(1)+(psi0(4)/k));
        b=0.5*(psi0(1)-(psi0(4)/k));
 
        t=a/(a*a+ia*ia);
        it=-ia/(a*a+ia*ia);
        T=pow(t,2)+pow(it,2);
        K=sqrt(2*m*(V-E))/hb;
        TAnaly=(4*K*K*k*k)/(pow((K*K+k*k),2)*pow(sinh(K*L),2)+(4*K*K*k*k));
        stockeT(i,(E*(1/eV)),T,TAnaly);
        }
    }