function J = einstein(z,E)
% function einstein(z,E)
%   This function needs to call function einstein_f1
%
% Written by Junke Guo, 07/05/2005

if abs(z-round(z)) > 5e-3       
    % Step 2 in the closure
    J1 = z.*pi./sin(z.*pi) - einstein_f1(z,E);
    
    % Step 3, the bracket in Eq. (22)
    n = ceil(z) + 2; k = 1:n; an = 1./(k-z); S1 = sum(an);
    gamma = 0.57721566490153286060651;
    fz = 1 - gamma - log(n-z) - 0.5./(n-z) + 1./12./(n-z).^2 - 1./120./(n-z).^4 + S1;
    
    % Step 4 
    S2 = 0; 
    for k = 1:n, S2 = S2 + (-1).^k./(z-k)./(z-k-1).*einstein_f1(z-k,E); end
    
    % Step 5
    J2 = -z.*pi./sin(z.*pi).*fz - einstein_f1(z,E).*(log(E) + 1./(z-1)) - z.*S2;

    J = [J1; J2];
    return
else
    n = round(z);
    J1 = (-1).^n.*(n.*log(E) - E + 1);
    J2 = (-1).^n.*(n./2.*log(E).^2 - E.*log(E) + E - 1);
    for k = 0:n-2
        J1 = J1 + (-1).^k.*nchoosek(n,k).*((E.^(k-n+1)-1)./(n-k-1));
        J2 = J2 + (-1).^k.*nchoosek(n,k).*(E.^(1+k-n).*log(E)./(n-k-1) + (E.^(1+k-n)-1)./(n-k-1).^2);
    end
    J = [J1; J2];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F1 = einstein_f1(z,E)
% This function is called by einstein.m
%
% Step 1 in the closure

E1 = E./(1-E);
k = 1:ceil(z)+4;
an = (-1).^k./(k-z).*E1.^(k-z);
F1 = E./E1.^z - z.*sum(an);