function [E1s,E1b,E1cf] = different_expint(z) sz = size(z); z = z(:); # First series expansion gm = 0.577215664901532860606512090082402431042159335; E1s = -gm - log(z); res = -z; ssum = res; k = 1; fflag = ones(size(res)); rrange = find(fflag); while(any(fflag)) rrange = rrange(find(fflag)); res_tmp = res(rrange); res_tmp = res_tmp*k.*(-z(rrange))./((k+1)^2); ssum(rrange) = ssum(rrange) + res_tmp; k++; res(rrange) = res_tmp; fflag = abs(res_tmp) > eps*abs(ssum(rrange)); endwhile E1s = E1s - ssum; E1s = reshape(E1s,sz); # Asymptotic series expansion z_b = z; E1b = exp(-z_b)./z_b; res = ones(size(z_b)); oldres = res; ssum = res; k = 0; fflag = ones(size(z_b)); rrange = find(fflag); while (k<1e3 && any(fflag)) rrange = rrange(find(fflag)); res_tmp = res(rrange); oldres_tmp = oldres(rrange); oldres_tmp = res_tmp; res_tmp = res_tmp*(k+1)./(-z_b(rrange)); ssum(rrange) = ssum(rrange) + res_tmp; k++; res(rrange) = res_tmp; oldres(rrange) = oldres_tmp; fflag = abs(oldres_tmp)>abs(res_tmp); endwhile E1b = E1b.*ssum; E1b = reshape(E1b,sz); # continued fraction, Lent'z modified algorithm f_new = 2^(-100)*ones(size(z)); C_new = f_new; D_new = zeros(size(z)); k = 0; fflag = ones(size(z)); rrange = find(fflag); Delta = ones(size(z)); C_old = ones(size(z)); D_old = ones(size(z)); f_old = ones(size(z)); while (k<1e3 && any(fflag)) rrange = rrange(find(fflag)); f_new_tmp = f_new(rrange); C_new_tmp = C_new(rrange); D_new_tmp = D_new(rrange); f_old = f_new_tmp; C_old = C_new_tmp; D_old = D_new_tmp; b = z(rrange)*(mod(k,2)==0) + (mod(k,2)==1); a = exp(-z(rrange))*(k==0) + ceil((k)/2)*(k>=1); D_new_tmp = b + a.*D_old; D_new_tmp = D_new_tmp.*(D_new_tmp~=0) + 2^(-100)*(D_new_tmp==0); C_new_tmp = b + a./C_old; C_new_tmp = C_new_tmp.*(C_new_tmp~=0) + 2^(-100)*(C_new_tmp==0); D_new_tmp = 1./D_new_tmp; Delta(rrange) = C_new_tmp.*D_new_tmp; f_new(rrange) = f_old.*Delta(rrange); C_new(rrange) = C_new_tmp; D_new(rrange) = D_new_tmp; fflag = abs(Delta(rrange)-1) > 2^(-50); k++; endwhile E1cf = reshape(f_new,sz); endfunction