% Grafische Darstellung der Funktion

%       -IMAG(1/KAPPA(Impuls q, Frequenz omega)

% eines homogenen Elektronengases nach Lindhard =
% random phase approximation RPA
%
% Impulse in Einheiten Fermi-Impulse
% Frequenzen in Einheiten Fermi-Frequenzen
%
% H. Sormann   6- 9-2001
% Last update 16-11-2004
%

clear all

warning off MATLAB:divideByZero

%
% Eingabe von rs in Bohr (ab) = Definition der Elektronendichte bzw. von (abkf)
  rs=input('rs in Einheiten Bohr (Vorschlag: 4.0) ');
%  rs=4.0;

  abkf=(9*pi/4)^(1/3) / rs;
%  pause

  const=4.0/pi/abkf;

% Definition des (omega)-Darstellungsbereiches:
%  omega:  0  bis ommax
%  ommax=input('omega-max Vorschlag: 6) = ');
  ommax=6.0;
% nomega Werte im omega-Raum:
  nomega=input('nomega = (Vorschlag: 100) ');
%  nomega=50;
  delom=ommax/nomega;
  omegaf= 0.0 : delom : ommax;
  for iom=1:nomega
    omega=omegaf(iom);
    if(omega==1) omegaf(iom)=1.0001;
    end
  end

% Definition des (k)-Darstellungsbereiches
%  k:  kmin bis kmax
%  kmin=input('kmin = ');
%  kmax=input('kmax = ');
  kmin=.01;
  kmax=2.0;
% nk Werte im k-Raum:
%  nk=input('nk = (Vorschlag: 100) ');
  nk=nomega;
  delk=(kmax-kmin)/nk;
  kf=kmin : delk : kmax;  

  for ik=1 : length(kf)
    k=kf(ik);
    if(k==0.5) k=0.4999;
    end
    if(k==1.0) k=0.9999;
    end
    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;
    
      lnarg1=abs(vec1 / vec2);
      lnarg2=abs(vec3 / vec4);
%    [lnarg1 ; lnarg2]'
%    pause
      ln1=log(lnarg1);
      ln2=log(lnarg2);
%    [ln1 ; ln2]'
%    pause
      vec5=1 - (omega-k^2)^2/(k^2)/4;
      vec6=1 - (omega+k^2)^2/(k^2)/4;
%    [k ; vec5 ; vec6]'
%    pause

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

% Imaginaerteil von kappa:
      kappa2(ik,iom)=0;
      if (k<=2)
        if(omega<= 2*k-k^2) 
          kappa2(ik,iom)=omega/(abkf*k^3);
        elseif(omega<= 2*k+k^2) 
          kappa2(ik,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(ik,iom)=(1 - (omega-k^2)^2/4/k^2)/abkf/k^3;
        end
      end
          
      werte(ik,iom)=kappa2(ik,iom)/(kappa1(ik,iom)^2 + kappa2(ik,iom)^2);
    end
  end

% Approx. Bestimmung der Maxima:
  for ik=1:length(kf)
    for iom=1:nomega
      if(kappa1(ik,iom)*kappa1(ik,iom+1)<0)
        if(kappa1(ik,iom)<kappa1(ik,iom+1))
          omnull=iom;
        else
          omnull=iom+1;
        end
      end
    end
%    fprintf('ik  nullst  2k-k^2   %5i   %9.6f   %9.6f\n',ik,omegaf(omnull),2*kf(ik)+kf(ik)^2);
%    pause
    if(omegaf(omnull) > 2*kf(ik)+kf(ik)^2)
%      werte(ik,omnull)=0.2;
      werte(ik,omnull)=0;
    end
    omnullf(ik)=omegaf(omnull);
  end

  [K,OM]=meshgrid(kf,omegaf);
  mesh(K,OM,werte');
  view(-67.5,40.0);
  hold on
  
  axis([kmin kmax 0.0 ommax 0.0 1]);
  xlabel('q')
  ylabel('omega')
  zlabel('IM(1/kappa)')
  title('INVERSE DER DIELEKTRISCHEN FUNKTION IN DER RPA');
  
  for i=1:length(omnullf)
    if(omnullf(i)==omnullf(end))
      break
    end
  end
  pause

  
  plot(kf(1:i),omnullf(1:i),'r*');
  hold off

%  title('Klicken Sie im Bildfenster: Tools --> Rotate-3D')
