/****************************************************************** !* Purpose: This program computes the psi function psi(z) * !* for a complex argument using subroutine CPSI * !* Input : x --- Real part of z * !* y --- Imaginary part of z * !* Output: PSR --- Real part of psi(z) * !* PSI --- Imaginary part of psi(z) * !* Examples: * !* x y Re[psi(z)] Im[psi(z)] * !* ------------------------------------------- * !* 3.0 2.0 1.16459152 .67080728 * !* 3.0 -2.0 1.16459152 -.67080728 * !* -3.0 2.0 1.39536075 2.62465344 * !* -3.0 -2.0 1.39536075 -2.62465344 * !* * !* -------------------------------------------------------------- * !* REFERENCE: "Fortran Routines for Computation of Special * !* Functions jin.ece.uiuc.edu/routines/routines.html" * !* * !* C++ Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !*****************************************************************/ #include #include void CPSI(double X, double Y, double *PSR, double *PSI) { /* ============================================= ! Purpose: Compute the psi function for a ! complex argument ! Input : x --- Real part of z ! y --- Imaginary part of z ! Output: PSR --- Real part of psi(z) ! PSI --- Imaginary part of psi(z) ! ============================================= */ double CT2,PI,RI,RR,TH,TM,TN,X0,X1,Y1,Z0,Z2; double A[9]; int K, N; A[1]=-0.83333333333333333e-01; A[2]=0.83333333333333333e-02; A[3]=-0.39682539682539683e-02; A[4]=0.41666666666666667e-02; A[5]=-0.83333333333333333e-01; A[6]=0.21092796092796093e-01; A[7]=-0.83333333333333333e-01; A[8]=0.4432598039215686; PI=3.141592653589793; if (Y == 0.0 && X == int(X) && X <= 0.0) { *PSR=1.0e+100; *PSI=0.0; } else { if (X < 0.0) { X1=X; Y1=Y; X=-X; Y=-Y; } X0=X; if (X < 8.0) { N=8-floor(X); X0=X+N; } if (X0 == 0.0 && Y != 0.0) TH=0.5*PI; if (X0 != 0.0) TH=atan(Y/X0); Z2=X0*X0+Y*Y; Z0=sqrt(Z2); *PSR=log(Z0)-0.5*X0/Z2; *PSI=TH+0.5*Y/Z2; for (K=1; K<9; K++) { *PSR=(*PSR)+A[K]*pow(Z2,(-K))*cos(2.0*K*TH); *PSI=(*PSI)-A[K]*pow(Z2,(-K))*sin(2.0*K*TH); } if (X < 8.0) { RR=0.0; RI=0.0; for (K=1; K<=N; K++) { RR=RR+(X0-K)/(pow((X0-K),2)+Y*Y); RI=RI+Y/(pow((X0-K),2)+Y*Y); } *PSR=(*PSR)-RR; *PSI=(*PSI)+RI; } if (X1 < 0.0) { TN=tan(PI*X); TM=tanh(PI*Y); CT2=TN*TN+TM*TM; *PSR=(*PSR)+X/(X*X+Y*Y)+PI*(TN-TN*TM*TM)/CT2; *PSI=(*PSI)-Y/(X*X+Y*Y)-PI*TM*(1.0+TN*TN)/CT2; X=X1; Y=Y1; } } } void main() { double X,Y,PSR,PSI; printf("\n x y Re[Psi(z)] Im[Psi(z)] \n"); printf(" ---------------------------------------------------\n"); X=3.0; Y=2.0; CPSI(X,Y,&PSR,&PSI); printf(" %f %f %e %e\n", X, Y, PSR,PSI); X=3.0; Y=-2.0; CPSI(X,Y,&PSR,&PSI); printf(" %f %f %e %e\n", X, Y, PSR,PSI); X=-3.0; Y=2.0; CPSI(X,Y,&PSR,&PSI); printf(" %f %f %e %e\n", X, Y, PSR,PSI); X=-3.0; Y=-2.0; CPSI(X,Y,&PSR,&PSI); printf(" %f %f %e %e\n", X, Y, PSR,PSI); } // end of file mcpsi.cpp