% Numerische und grafische Auswertung der Funktion Imag(1/kappa).
% "kappa" = dielektrische Funktion eines homogenen Elektronengases 
% nach Lindhard = random phase approximation RPA
%
% Impulse in Einheiten Fermi-Impulse
% Frequenzen in Einheiten Fermi-Frequenzen
%
% Realteil und Imaginaerteil werden getrennt berechnet und
% graphisch dargestellt.
%
% H. Sormann   6- 9-2001
% Last update 12-11-2003
%

  clc

  clear all

%  ndelay=input('ndelay (Vorschlag: 10000) = ');
  ndelay=15000;
%
% Eingabe von rs in Bohr (ab) = Definition der Elektronendichte bzw. von (abkf)
  fprintf('RS=2.0 mit AL-EXP.!!\n');
  rs=input('rs in Einheiten Bohr (rs>=3.0) ');
  fprintf('RS IN BOHR = %5.2f\n', rs);

  abkf=(9*pi/4)^(1/3) / rs;
  const=4.0/pi/abkf;

% Im Fall rs=2 : Einlesen exp. Werte fuer Al (Batson et al 1976):
  if(rs==2.0)
    al_exp=load('Im1kappa.Batson_76');
  end  

% Definition des (omega)-Darstellungsbereiches:
%  omega:  0  bis ommax
%  ommax=input('omega-max Vorschlag: 6) = ');
  ommax=6.0;

% nomega Werte im omega-Raum:
  nomega=400;

  delom=ommax/nomega;

 omegaf= 0.0 : delom : ommax;
  for iom=1:nomega
    omega=omegaf(iom);
    if(omega==1) omegaf(iom)=1.0001;
    end
  end

  max_omegaf=zeros(1,72);
  max_wertef=zeros(1,72);
  kf=zeros(1,72);

  knull=0;
  nullstf=0;
  istop=1;
  fprintf('ENTER fuer Programmstart\n');
  pause
  for ik=1:72 

    k=0.025*ik;
    fprintf('       K/K-FERMI = %6.3f\n',k);

    if(k==0.5) k=0.4999;
    end
    if(k==1.0) k=0.9999;
    end
    kf(ik)=k;

    for iom=1:nomega+1
      omega=omegaf(iom);

% Realteil von kappa:
      vec1=omega-2*k-k^2;
      vec2=omega+2*k-k^2;
      vec3=omega+2*k+k^2;
      vec4=omega-2*k+k^2;

      if(vec1==0)
        vec1=0.0001;
      end
      if(vec2==0)
        vec2=0.0001;
      end
      if(vec3==0)
        vec3=0.0001;
      end
      if(vec4==0)
        vec4=0.0001;
      end
    
      lnarg1=abs(vec1 / vec2);
      lnarg2=abs(vec3 / vec4);
      ln1=log(lnarg1);
      ln2=log(lnarg2);
      vec5=1 - (omega-k^2)^2/(k^2)/4;
      vec6=1 - (omega+k^2)^2/(k^2)/4;

      kappa1(iom)=1 + const /k^2 *(0.5 + (vec5 * ln1 + vec6 * ln2) /k/4);

% Imaginaerteil von kappa:
      kappa2(iom)=0;
      if (k<=2)
        if(omega<= 2*k-k^2) 
          kappa2(iom)=omega/(abkf*k^3);
        elseif(omega<= 2*k+k^2) 
          kappa2(iom)=(1 - (omega-k^2)^2/4/k^2)/abkf/k^3;
        end
      elseif(k>2)
        if((omega >= k^2-2*k) & (omega<= 2*k+k^2))
          kappa2(iom)=(1 - (omega-k^2)^2/4/k^2)/abkf/k^3;
        end
      end
          
    end

    werte=kappa2./(kappa1.^2 + kappa2.^2);
    plot(omegaf',werte','r');
    legend('electron-hole scattering');
    wertmax=6;
    axis([0 ommax 0 wertmax]);
    xlabel('OMEGA');
    ylabel('IMAG.(1/KAPPA)');
    title('-IMAG.TEIL 1/LINDHARD-FKT.');
    hold on
    for iii=1:length(werte)
      if(werte(iii)==max(werte))
        max_omega=omegaf(iii);
      end
    end
    plot(max_omega,max(werte),'ro')
    legend('electron-hole-scattering','Maximum der el-hole-Kurve');
    max_omegaf(ik)=max_omega;
    max_wertef(ik)=max(werte);

% Bestimmung der Nullstellen des Realteils:

    for iom=1:length(kappa1)-1
      if(kappa1(iom)*kappa1(iom+1) <= 0.0)
        if(abs(kappa1(iom))<abs(kappa1(iom+1)))
          nullst=omegaf(iom);
          fkappa2=kappa2(iom);
        else
          nullst=omegaf(iom+1);
          fkappa2=kappa2(iom+1);
        end
        if(k<2.0)
        knull(ik)=k;
        nullstf(ik)=nullst;
        end
        if(fkappa2==0) 
          vertikal=[nullst 0 ; nullst wertmax];
          plot(vertikal(:,1),vertikal(:,2),'b');
          legend('electron-hole-scattering','Maximum der el-hole-Kurve','plasmon excitation');
        end
      end
    end
    if(abs(nullst-(2*k+k^2))<0.008)
      fprintf('************************** KRITISCHES K/K-FERMI = %5.3f\n', k);
      fprintf('*** WEITER MIT ENTER.\n');
      title('KRITISCHE WELLENZAHL');
      krit=k;
      pause  
    end
    if(abs(nullst-(2*k+k^2))<0.5)
      if(istop==1)
        fprintf('************************************  K/K-FERMI = %5.3f\n', k);
        fprintf('*** PROGRAMM-FORTSETZUNG:\n');
        fprintf('    Command window anklicken und ENTER.\n'); 
        pause
        istop=0;
      end  
    end
    kfeld(ik)=k;

    for idelay=1:ndelay
      delay=sin(0.1*idelay);
    end
    drawnow
    hold off
  end

  fprintf('************************************  K/K-FERMI = %5.3f\n', k);
  fprintf('*** WEITER ZU OMEGA-K-GRAFIK MIT ENTER.\n');
  pause

  knull(length(knull))=[];
  nullstf(length(nullstf))=[];

  plot(knull,nullstf,'b',kfeld,-2.*kfeld+kfeld.^2,'r--',kfeld,2.*kfeld+kfeld.^2,'r--');
  axis([-inf  inf  0   inf]);
  text(1.2,2.0,'Imag(KAPPA) UNGLEICH NULL');
  title('NULLSTELLEN-DISKUSSION  LINDHARD-FKT.');
  xlabel('WELLENZAHL IN EINH. K-FERMI');
  ylabel('FREQUENZ');
  hold on
  kc=[krit 0 ; krit 11];
  plot(kc(:,1),kc(:,2),'k-.');
  text(krit-0.22,9.0,'krit. Wellenzahl');
  text(0.05,3.0,'Plasmonen');
  fprintf('ENTER FUER DIE MAXIMUMSKURVE DER ELEK-HOLE-PROFILE.\n');
  pause
  plot(kf,max_omegaf,'r.')
  title('ROTE PUNKTE = FREQUENZEN DER MAXIMA DER ELEK-HOLE-PROFILE')
  hold off

  if(rs==2.0)
    fprintf('*** WEITER ZU: VERGLEICH MIT AL-EXPERIMENT\n');
    pause

    plot(al_exp(:,1),1.0+0.2*al_exp(:,2),'ro')
    legend('EXP: Batson et al, 1976)');
    hold on
    plot(knull.^2,nullstf,'b');
    axis([0 1.5 0 3.5]);
    grid on
    title('NULLSTELLEN-DISKUSSION  LINDHARD-FKT.');
    xlabel('(WELLENZAHL IN EINH. K-FERMI)^2');
    ylabel('FREQUENZ');
    kc=[krit 0 ; krit 11];
    plot(kc(:,1).^2,kc(:,2),'k-.');
    text(0.05,3.2,'THEORIE = RPA');
    text(0.05,2.7,'RS = 2.0 (ALUMINIUM)');
    text(0.05,1.7,'Plasmonen');
    fprintf('ENTER FUER DIE MAXIMUMSKURVE DER ELEK-HOLE-PROFILE.\n');
    plot(kf.^2,max_omegaf,'r.')
    title('ROTE PUNKTE = FREQUENZEN DER MAXIMA DER ELEK-HOLE-PROFILE')
    hold off
  end
