Optimization for Mermin, Lindhard and vos_g functions

src/lib/ELF/Lindhard.m

@@ -1,4 +1,4 @@
-function eps=Lindhard(q,omega,gamma,omega0)
+function [eps,epsimag]=Lindhard(q,omega,gamma,omega0)
 
 %     sq = numel(q);
 %     sw = numel(omega);
@@ -14,18 +14,23 @@
     chi = sqrt(1.0 / (pi*v_f));
 
     z1_1 = omega./(q*v_f);
+    z1_1(isnan(z1_1)) = 0;
 
-    z1 = bsxfun(@plus,bsxfun(@plus,z1_1,z),1j*gamma./(q*v_f));
-    d1 = vos_g(z1);
-    z2 = bsxfun(@plus,bsxfun(@minus,z1_1,z),1j*gamma./(q*v_f));
-    d2 = vos_g(z2);
+    gq = gamma./(q*v_f); gq(isnan(gq)) = 0;
+    [reD1, imD1] = vos_g(z1_1 + z, gq);
+    [reD2, imD2] = vos_g(z1_1 - z, gq);
 
-    zzz = z.^3;
-    red1_d2 = real(d1) - real(d2);
-    imd1_d2 = imag(d1) - imag(d2);
+    red1_d2 = reD1 - reD2;
+    imd1_d2 = imD1 - imD2;
 
-    epsreal = bsxfun(@plus,1,bsxfun(@times,red1_d2 ,chi^2./(4 * zzz)));
-    epsimag = bsxfun(@times,imd1_d2 ,chi^2./(4 * zzz));
-    eps = complex(epsreal,epsimag);
+    chizzz = chi^2./(z.*z.*z * 4);
+    epsreal = 1 + red1_d2 .* chizzz;
+    epsimag = imd1_d2 .* chizzz;
+    if nargout == 1
+        eps = complex(epsreal,epsimag);
+    else
+        eps = epsreal;
+        % epsimag = epsimag
+    end
 
 end

src/lib/ELF/Mermin.m

@@ -6,11 +6,12 @@
 
     g_over_w = gamma ./ omega;
     z1 = complex(1,g_over_w); % omega should be unequal 0
-    z2 = bsxfun(@minus,Lindhard(q, omega, gamma, om_at_q),1);
-    z3 = bsxfun(@minus,Lindhard(q, zeros(size(omega)), 0.0, om_at_q),1);
 
-    top = bsxfun(@times,z1,z2);
-    bottom = bsxfun(@plus,1,bsxfun(@times,complex(0,g_over_w),z2)./z3);
-    eps = bsxfun(@plus,1,top./bottom);
+    z2 = Lindhard(q, omega, gamma, om_at_q) - 1;
+    z3 = Lindhard(q, zeros(size(omega)), 0.0, om_at_q) - 1;
+
+    top = z1 .* z2;
+    bottom = complex(0,g_over_w).*z2./z3 + 1;
+    eps = 1 + top./bottom;
 
 end

src/lib/ELF/vos_g.m

@@ -1,24 +1,52 @@
-function eps=vos_g(z)
+function [out, outim]=vos_g(z, img_z)
 
-    zplus1 = z + ones(size(z));
-    zminus1 = z - ones(size(z));
-
-    dummy1 = log( abs(zplus1) ./ abs(zminus1) );
-    dummy2 = angle(zplus1) - angle(zminus1);
-
-    real_sq = (real(z)).^2;
-    imag_sq = (imag(z)).^2;
-
-    reim1 = bsxfun(@minus,1,real_sq - imag_sq);
-
-    outreal_1 = real(z) + 0.5*bsxfun(@times,reim1,dummy1);
-    outreal = outreal_1 + bsxfun(@times,bsxfun(@times,real(z),imag(z)),dummy2);
+    if nargin < 2
+        reZ = real(z);
+        imgZ = imag(z);
+    else
+        reZ = z;
+        imgZ = img_z;
+    end
+
+    zplus1 = z + 1;
+    zminus1 = z - 1;
 
-    outimag_1 = imag(z) + 0.5*bsxfun(@times,reim1,dummy2);
-    outimag = outimag_1 - bsxfun(@times,bsxfun(@times,real(z),imag(z)),dummy1);
+    if any(any(imgZ ~= 0))
+        imgZ2 = imgZ.*imgZ;
+        % dummy1 = log( abs(zplus1) ./ abs(zminus1) );
+        dummy1 = log( sqrt((zplus1.*zplus1+imgZ2) ./ (zminus1.*zminus1+imgZ2)) );
+    %     dummy2 = angle(zplus1) - angle(zminus1);
+        dummy2 = atan2(imgZ, zplus1) - atan2(imgZ, zminus1);
+
+        reim1 = 1 - (reZ.^2 - imgZ2);
+
+        outreal_1 = reZ + 0.5*reim1.*dummy1;
+        outreal = outreal_1 + reZ.*imgZ.*dummy2;
+
+        outimag_1 = imgZ + 0.5*reim1.*dummy2;
+        outimag = outimag_1 - reZ.*imgZ.*dummy1;
+    else
+%         imgZ2 = imgZ.*imgZ;
+        dummy1 = log( abs(zplus1) ./ abs(zminus1) );
+    %     dummy2 = angle(zplus1) - angle(zminus1);
+        dummy2 = atan2(0, zplus1) - atan2(0, zminus1);
+
+        reim1 = 1 - reZ.^2;
+
+        outreal_1 = reZ + 0.5*reim1.*dummy1;
+        outreal = outreal_1;%; + reZ.*imgZ.*dummy2;
+
+        outimag_1 = 0.5*reim1.*dummy2;% + imgZ
+        outimag = outimag_1;% - reZ.*imgZ.*dummy1;        
+    end
 
 %     outreal = real(z) + 0.5* (1.0 - (real(z)*real(z) - imag(z)*imag(z))) *dummy1 + real(z)*imag(z)*dummy2;
-%     outimag = imag(z) + 0.5*(1.0 - (real(z)*real(z) - imag(z)*imag(z)))*dummy2 - real(z)*imag(z)*dummy1;
+%     outimag = imag(z) + 0.5* (1.0 - (real(z)*real(z) - imag(z)*imag(z))) *dummy2 - real(z)*imag(z)*dummy1;
 
-    eps=complex(outreal,outimag);
+    if nargout == 1
+        out = complex(outreal,outimag);
+    else
+        out = outreal;
+        outim = outimag;
+    end
 end