function zhf=hypergeom(a,b,c,z)
%This program is a direct conversion of the corresponding Fortran program in
%S. Zhang & J. Jin "Computation of Special Functions" (Wiley, 1996).%online: http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html%%Converted by f2matlab open source project:%online: https://sourceforge.net/projects/f2matlab/% written by Ben Barrowes (barrowes@alum.mit.edu)%

%     ============================================================
%     Purpose: This program computes hypergeometric function for
%     a complex argument, F(a,b,c,z), using subroutine
%     HYGFZ
%     Input :  a --- Parameter
%     b --- Parameter
%     c --- Parameter,  c <> 0,-1,-2,...
%     x --- Real part of complex argument z
%     y --- Imaginary part of complex argument z
%     ( z = x+iy )
%     Output:  ZHF --- F(a,b,c,z)
%     Examples:
%     a     b     c    z = x+ iy             F(a,b,c,z)
%     --------------------------------------------------------------
%     3.2   1.8   6.7   1.0+0.0 i    .54689992D+01+.00000000D+00 i

%     3.2  -1.8   6.7   1.0+0.0 i    .33750635D+00+.00000000D+00 i
%     -5.0   3.3   6.7   5.2+4.8 i    .11682745D+03+.60389104D+03 i
%     3.3  -6.0   3.7   5.2-4.8 i    .17620425D+05+.38293812D+05 i
%     -7.0   3.3  -3.7   5.2-4.8 i   -.11772779D+11-.14382286D+11 i
%     4.3  -8.0  -3.7   5.2+4.8 i    .13161188D+13-.10129870D+12 i
%     3.3   5.8   6.7   0.2+0.1 i    .17330557D+01+.63401030D+00 i
%     3.5  -2.4   6.7   0.2+0.5 i    .64762241D+00-.52110507D+00 i
%     3.3   4.3   6.7   0.8+0.3 i   -.14830086D+01+.83744258D+01 i
%     7.0   5.0   4.1   3.0-1.0 i   -.40376095D-02-.29566326D-02 i
%     5.0   7.0   4.1   3.0-1.0 i   -.40376095D-02-.29566326D-02 i
%     3.5   1.2   9.7   0.6+0.9 i    .10343044D+01+.54473814D+00 i
%     2.1   5.4   9.7   0.5+0.7 i    .68850442D+00+.12274187D+01 i
%     8.7   3.2   6.7   0.5+0.7 i   -.90046505D+00-.11198900D+01 i
%     8.7   2.7   6.7   0.6+0.9 i   -.46083890D+00-.54575701D+00 i
%     ============================================================



%     ======================================================
%     Purpose: Compute the hypergeometric function for a
%     complex argument, F(a,b,c,z)
%     Input :  a --- Parameter
%     b --- Parameter
%     c --- Parameter,  c <> 0,-1,-2,...
%     z --- Complex argument
%     Output:  ZHF --- F(a,b,c,z)
%     Routines called:
%     (1) MYGAMMA for computing gamma function
%     (2) PSI for computing psi function
%     ======================================================



gc=[];gcab=[];gca=[];gcb=[];g1=[];g2=[];g3=[];ga=[];gb=[];gam=[];gbm=[];pa=[];pb=[];gabc=[];gab=[];gba=[];pca=[];pac=[];gm=[];t0=[];g0=[];gcbk=[];
zw=0.0;
x=real(z);
y=imag(z);
eps=1.0d-15;
l0=c == fix(c)&c < 0.0d0;
l1=abs(1.0d0-x) < eps&y == 0.0d0&c-a-b <= 0.0d0;
l2=abs(z+1.0d0) < eps&abs(c-a+b-1.0d0) < eps;
l3=a == fix(a)&a < 0.0d0;
l4=b == fix(b)&b < 0.0d0;
l5=c-a == fix(c-a)&c-a <= 0.0d0;
l6=c-b == fix(c-b)&c-b <= 0.0d0;
aa=a;
bb=b;
a0=abs(z);
if (a0 > 0.95d0) eps=1.0d-8; end;
pi=3.141592653589793d0;
el=.5772156649015329d0;
if (l0|l1);
'the hypergeometric series is divergent',
return;
end;
if (a0 == 0.0d0|a == 0.0d0|b == 0.0d0);
zhf=complex(1.0d0,0.0d0);
elseif (z == 1.0d0&c-a-b > 0.0d0);
[gc]=mygamma(c,gc);
[gcab]=mygamma(c-a-b,gcab);
[gca]=mygamma(c-a,gca);
[gcb]=mygamma(c-b,gcb);
zhf=gc.*gcab./(gca.*gcb);
elseif (l2);
g0=sqrt(pi).*2.0d0.^(-a);
[g1]=mygamma(c,g1);
[g2]=mygamma(1.0d0+a./2.0d0-b,g2);
[g3]=mygamma(0.5d0+0.5d0.*a,g3);
zhf=g0.*g1./(g2.*g3);
elseif (l3|l4);
if (l3) nm=fix(abs(a)); end;
if (l4) nm=fix(abs(b)); end;
zhf=complex(1.0d0,0.0d0);
zr=complex(1.0d0,0.0d0);
for  k=1:nm;
zr=zr.*(a+k-1.0d0).*(b+k-1.0d0)./(k.*(c+k-1.0d0)).*z;
zhf=zhf+zr;
end;
elseif (l5|l6);
if (l5) nm=fix(abs(c-a)); end;
if (l6) nm=fix(abs(c-b)); end;
zhf=complex(1.0d0,0.0d0);
zr=complex(1.0d0,0.0d0);
for  k=1:nm;
zr=zr.*(c-a+k-1.0d0).*(c-b+k-1.0d0)./(k.*(c+k-1.0d0)).*z;
zhf=zhf+zr;
end;
zhf=(1.0d0-z).^(c-a-b).*zhf;
elseif (a0 <= 1.0d0);
if (x < 0.0d0);
z1=z./(z-1.0d0);
if (c > a&b < a&b > 0.0);
a=bb;
b=aa;
end;
zc0=1.0d0./((1.0d0-z).^a);
zhf=complex(1.0d0,0.0d0);
zr0=complex(1.0d0,0.0d0);
for  k=1:500;
zr0=zr0.*(a+k-1.0d0).*(c-b+k-1.0d0)./(k.*(c+k-1.0d0)).*z1;
zhf=zhf+zr0;
if (abs(zhf-zw) < abs(zhf).*eps) break; end;
zw=zhf;
end;
     zhf=zc0.*zhf;
elseif (a0 >= 0.90d0);
gm=0.0d0;
mcab=fix(c-a-b+eps.*(abs(1.0d0).*sign(c-a-b)));
if (abs(c-a-b-mcab) < eps);
m=fix(c-a-b);
[ga]=mygamma(a,ga);
[gb]=mygamma(b,gb);
[gc]=mygamma(c,gc);
[gam]=mygamma(a+m,gam);
[gbm]=mygamma(b+m,gbm);
[pa]=mypsi(a,pa);
[pb]=mypsi(b,pb);
if (m ~= 0) gm=1.0d0; end;
for  j=1:abs(m)-1;
gm=gm.*j;
end;
rm=1.0d0;
for  j=1:abs(m);
rm=rm.*j;
end;
zf0=complex(1.0d0,0.0d0);
zr0=complex(1.0d0,0.0d0);
zr1=complex(1.0d0,0.0d0);
sp0=0.d0;
sp=0.0d0;
if (m >= 0);
zc0=gm.*gc./(gam.*gbm);
zc1=-gc.*(z-1.0d0).^m./(ga.*gb.*rm);
for  k=1:m-1;
zr0=zr0.*(a+k-1.d0).*(b+k-1.d0)./(k.*(k-m)).*(1.d0-z);
zf0=zf0+zr0;
end;
for  k=1:m;
sp0=sp0+1.0d0./(a+k-1.0d0)+1.0./(b+k-1.0d0)-1.d0./k;
end;
zf1=pa+pb+sp0+2.0d0.*el+log(1.0d0-z);
for  k=1:500;
sp=sp+(1.0d0-a)./(k.*(a+k-1.0d0))+(1.0d0-b)./ (k.*(b+k-1.0d0));
sm=0.0d0;
for  j=1:m;
sm=sm+(1.0d0-a)./((j+k).*(a+j+k-1.0d0)) +1.0d0./(b+j+k-1.0d0);
end;
zp=pa+pb+2.0d0.*el+sp+sm+log(1.0d0-z);
zr1=zr1.*(a+m+k-1.0d0).*(b+m+k-1.0d0)./(k.*(m+k)) .*(1.0d0-z);
zf1=zf1+zr1.*zp;
if (abs(zf1-zw) < abs(zf1).*eps) break; end;
zw=zf1;
end;
       zhf=zf0.*zc0+zf1.*zc1;
elseif (m < 0);
m=-m;
zc0=gm.*gc./(ga.*gb.*(1.0d0-z).^m);
zc1=-(-1).^m.*gc./(gam.*gbm.*rm);
for  k=1:m-1;
zr0=zr0.*(a-m+k-1.0d0).*(b-m+k-1.0d0)./(k.*(k-m)) .*(1.0d0-z);
zf0=zf0+zr0;
end;
for  k=1:m;
sp0=sp0+1.0d0./k;
end;
zf1=pa+pb-sp0+2.0d0.*el+log(1.0d0-z);
for  k=1:500;
sp=sp+(1.0d0-a)./(k.*(a+k-1.0d0))+(1.0d0-b)./(k.* (b+k-1.0d0));
sm=0.0d0;
for  j=1:m;
sm=sm+1.0d0./(j+k);
end;
zp=pa+pb+2.0d0.*el+sp-sm+log(1.0d0-z);
zr1=zr1.*(a+k-1.d0).*(b+k-1.d0)./(k.*(m+k)).*(1.d0-z);
zf1=zf1+zr1.*zp;
if (abs(zf1-zw) < abs(zf1).*eps) break; end;
zw=zf1;
end;
       zhf=zf0.*zc0+zf1.*zc1;
end;
else;
[ga]=mygamma(a,ga);
[gb]=mygamma(b,gb);
[gc]=mygamma(c,gc);
[gca]=mygamma(c-a,gca);
[gcb]=mygamma(c-b,gcb);
[gcab]=mygamma(c-a-b,gcab);
[gabc]=mygamma(a+b-c,gabc);
zc0=gc.*gcab./(gca.*gcb);
zc1=gc.*gabc./(ga.*gb).*(1.0d0-z).^(c-a-b);
zhf=complex(0.0d0,0.0d0);
zr0=zc0;
zr1=zc1;
for  k=1:500;
zr0=zr0.*(a+k-1.d0).*(b+k-1.d0)./(k.*(a+b-c+k)).*(1.d0-z);
zr1=zr1.*(c-a+k-1.0d0).*(c-b+k-1.0d0)./(k.*(c-a-b+k)) .*(1.0d0-z);
zhf=zhf+zr0+zr1;
if (abs(zhf-zw) < abs(zhf).*eps) break; end;
zw=zhf;
end;
      zhf=zhf+zc0+zc1;
end;
else;
z00=complex(1.0d0,0.0d0);
if (c-a < a&c-b < b);
z00=(1.0d0-z).^(c-a-b);
a=c-a;
b=c-b;
end;
zhf=complex(1.0d0,0.d0);
zr=complex(1.0d0,0.0d0);
for  k=1:1500;
zr=zr.*(a+k-1.0d0).*(b+k-1.0d0)./(k.*(c+k-1.0d0)).*z;
zhf=zhf+zr;
if (abs(zhf-zw) <= abs(zhf).*eps) break; end;
zw=zhf;
end;
    zhf=z00.*zhf;
end;
elseif (a0 > 1.0d0);
mab=fix(a-b+eps.*(abs(1.0d0).*sign(a-b)));
if (abs(a-b-mab) < eps&a0 <= 1.1d0) b=b+eps; end;
if (abs(a-b-mab) > eps);
[ga]=mygamma(a,ga);
[gb]=mygamma(b,gb);
[gc]=mygamma(c,gc);
[gab]=mygamma(a-b,gab);
[gba]=mygamma(b-a,gba);
[gca]=mygamma(c-a,gca);
[gcb]=mygamma(c-b,gcb);
zc0=gc.*gba./(gca.*gb.*(-z).^a);
zc1=gc.*gab./(gcb.*ga.*(-z).^b);
zr0=zc0;
zr1=zc1;
zhf=complex(0.0d0,0.0d0);
for  k=1:500;
zr0=zr0.*(a+k-1.0d0).*(a-c+k)./((a-b+k).*k.*z);
zr1=zr1.*(b+k-1.0d0).*(b-c+k)./((b-a+k).*k.*z);
zhf=zhf+zr0+zr1;
if (abs((zhf-zw)./zhf) <= eps) break; end;
zw=zhf;
end;
    zhf=zhf+zc0+zc1;
else;
if (a-b < 0.0d0);
a=bb;
b=aa;
end;
ca=c-a;
cb=c-b;
nca=fix(ca+eps.*(abs(1.0d0).*sign(ca)));
ncb=fix(cb+eps.*(abs(1.0d0).*sign(cb)));
if (abs(ca-nca) < eps|abs(cb-ncb) < eps) c=c+eps; end;
[ga]=mygamma(a,ga);
[gc]=mygamma(c,gc);
[gcb]=mygamma(c-b,gcb);
[pa]=psi(a,pa);
[pca]=psi(c-a,pca);
[pac]=psi(a-c,pac);
mab=fix(a-b+eps);
zc0=gc./(ga.*(-z).^b);
[gm]=mygamma(a-b,gm);
zf0=gm./gcb.*zc0;
zr=zc0;
for  k=1:mab-1;
zr=zr.*(b+k-1.0d0)./(k.*z);
t0=a-b-k;
[g0]=mygamma(t0,g0);
[gcbk]=mygamma(c-b-k,gcbk);
zf0=zf0+zr.*g0./gcbk;
end;
if (mab == 0) zf0=complex(0.0d0,0.0d0); end;
zc1=gc./(ga.*gcb.*(-z).^a);
sp=-2.0d0.*el-pa-pca;
for  j=1:mab;
sp=sp+1.0d0./j;
end;
zp0=sp+log(-z);
sq=1.0d0;
for  j=1:mab;
sq=sq.*(b+j-1.0d0).*(b-c+j)./j;
end;
zf1=(sq.*zp0).*zc1;
zr=zc1;
rk1=1.0d0;
sj1=0.0d0;
for  k=1:10000;
zr=zr./z;
rk1=rk1.*(b+k-1.0d0).*(b-c+k)./(k.*k);
rk2=rk1;
for  j=k+1:k+mab;
rk2=rk2.*(b+j-1.0d0).*(b-c+j)./j;
end;
sj1=sj1+(a-1.0d0)./(k.*(a+k-1.0d0))+(a-c-1.0d0)./ (k.*(a-c+k-1.0d0));
sj2=sj1;
for  j=k+1:k+mab;
sj2=sj2+1.0d0./j;
end;
zp=-2.0d0.*el-pa-pac+sj2-1.0d0./(k+a-c) -pi./tan(pi.*(k+a-c))+log(-z);
zf1=zf1+rk2.*zr.*zp;
ws=abs(zf1);
if (abs((ws-w0)./ws) < eps) break; end;
w0=ws;
end;
    zhf=zf0+zf1;
end;
end;
  a=aa;
b=bb;
if (k > 150) ; end;
  %format(1x,'warning% you should check the accuracy');

