% Programme pour calculer l influence mutuelle
% entre atomes polarisables
%
% Licence de Physique, UL5A
% Universite Paul Sabatier
% Licence de Physique UL5, TD numeriques, 2003-2004
% Version 08/02/2004 LK

% ------------ Initialisation  ------------
clear all
clf
nom='??'
 
% ------- Champ electrique exterieur ------
E0y = ??
E0z = ??

% ------- polarisabilite des atomes -------
alp1 = 3.e-27;
alp2 = 3.e-27;

% ---------- Position des atomes ----------
y1 = 0;
z1 = 0;
y2 = ??
z2 = ??

%  ------ Calcul de la distance entre atomes ------
r12y = y2-y1;
r12z = z2-z1;
r12 = ??

% ------- Initialisation du champ electrique ------
   E1y = E0y;
   E1z = E0z;

   E2y = E0y;
   E2z = E0z;

%  -------------- Iterations  ---------------------
for iter = 1:10
    E2(iter) = sqrt (E2y^2 + E2z^2);  % module 
    E1(iter) = sqrt (E1y^2 + E1z^2);  % module

    % Moment dipolaire de l atome 1 

    p1y = ??        % selon y
    p1z = ??        % selon z
    p1(iter) = sqrt (p1y^2+p1z^2);   % module

    % Champ produit en 2

    p1r12 = p1y * r12y + p1z * r12z;
    E2yp = E0y +(3 * r12y * p1r12 - p1y * r12^2) / r12^5;
    E2zp = E0z +(3 * r12z * p1r12 - p1z * r12^2) / r12^5;
    E2p = sqrt (E2yp^2 + E2zp^2);   % module

    DelE2 (iter) = abs (E2p - E2(iter));
    E2y = E2yp;
    E2z = E2zp;

   % Moment dipolaire de l atome 2 

    p2y = ??          % selon y
    p2z = ??          % selon z
    p2(iter) = sqrt (p2y^2 + p2z^2);

    % Champ produit en 1
    
    p2r12=p2y*r12y+p2z*r12z;
    E1yp = E0y +(3* r12y * p2r12 - p2y * r12^2) / r12^5;
    E1zp = E0z +(3* r12z * p2r12 - p2z * r12^2) / r12^5;
    E1p = sqrt (E1yp^2 + E1zp^2);  % module
    
    DelE1 (iter) = abs (E1p - E1(iter));
    E1y = E1yp;
    E1z = E1zp;
end

%  ---------------- Partie graphique ---------------------
semilogy (DelE1)
title ('DelE1')
pause

plot (E1)
title ('E1')
pause

semilogy (DelE2)
title ('DelE2')
pause

plot (E2)
title('E2')
pause

plot (p1)
title('p1')
pause

plot (p2)
title ('p2')
