/******************************************************************
!*      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 <stdio.h>
#include <math.h>

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