% Written by Min-Hua Huang
% June 9, 2011
% This is a file to calculate the "h" function.
% input: h(s,c,n2); s is the base parameter; c is the power parameter;
% n2 is the number of summation terms for calculating h when s is a
% negative integer.
% Due to the limited computational capacity, when abs(c)>30,
% the result may be invalid.
function fv=h(s,c,n2)
format long;  n1=250;  % n1 is the number of h's expansion terms. 
%if abs(c)>30
%   fprintf('Due to the limited computational capacity, when abs(c)>30,\nthe result may be invalid!\n');
%end
%************** Case 1: s is a negative integer  *******************
if sign(s)<0 && mod(s,1)==0
   sum2=log(abs(c)); d=s;
     for j=0:1:n2
         s=j;
         start1=1/(s+1); sum1=start1;
         for i=1:1:n1
             start1=start1*(c)/(s+1+i);
             sum1=sum1+start1;
         end         
         sum2=sum2+sum1;
     end
   value=sum2*exp(-c);
   %********************     s=-1   ********************************** 
   if d==-1
      hvalue=value;
   else
   %******************** s=-2,-3,-4,...... ***************************
      start2=1/(d+1); sum3=start2;
      for i=1:1:(-d-2)
          start2=start2*(-c)/(d+1+i);
          sum3=sum3+start2;
      end
      hvalue=sum3+start2*(-c)*value;
   end 
else
%************** Case 2: s is a not negative integer ****************
   start1=1/(s+1); sum2=start1;
    for i=1:1:n1
        start1=start1*(-c)/(s+1+i);
        sum2=sum2+start1;
    end
   hvalue=sum2; 
end
fv=hvalue;