#include "fintrf.h"
C
#if 0
C     
C     
#endif
	  SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
	  IMPLICIT NONE
	  INTEGER PLHS(*), PRHS(*)
	  INTEGER MXGETM, MXGETN
      INTEGER NLHS, NRHS
	  INTEGER M, N


C      Check for proper number of arguments. 
      IF (NRHS .NE. 3) THEN
         CALL MEXERRMSGTXT('Three inputs required.')
      ENDIF

C	  Check that input #1 is a scalar.
      M = MXGETM(PRHS(1))
      N = MXGETN(PRHS(1))
      IF(M .NE. 1 .OR. N .NE. 1) THEN
         CALL MEXERRMSGTXT('Input #1 is not a scalar.')
      ENDIF

C	  Check that input #2 is a scalar.
      M = MXGETM(PRHS(2))
      N = MXGETN(PRHS(2))
      IF(M .NE. 1 .OR. N .NE. 1) THEN
         CALL MEXERRMSGTXT('Input #2 is not a scalar.')
      ENDIF

C	  Check that input #3 is a column vector.
      N = MXGETN(PRHS(3))
      IF( N .NE. 1) THEN
         CALL MEXERRMSGTXT('Input #3 is not a column vector.')
      ENDIF

	  M = MXGETM(PRHS(3))
	
C	CALL INTERMEDIATE GATEWAY:
	  CALL GATEWAY(PLHS,PRHS,M)

	  END 
C ------ GATEWAY SUBROUTINE

	  SUBROUTINE GATEWAY(PLHS, PRHS, ROWS)
	  IMPLICIT NONE
      INTEGER PLHS(*), PRHS(*)
	  INTEGER MXCREATEDOUBLEMATRIX, MXGETPR, MXGETPI
	  INTEGER MXISCOMPLEX
      INTEGER ROWS, ISCOMPLEX
	  DOUBLE PRECISION SIG2, T1(ROWS), ALPHA, OMEGA
	  DOUBLE PRECISION RUPSI(ROWS), IUPSI(ROWS)
      COMPLEX*16 GAMMA
      
      PLHS(1) = MXCREATEDOUBLEMATRIX(ROWS,1,1)
      CALL MXCOPYPTRTOREAL8(MXGETPR(PRHS(1)), ALPHA ,1)  

C ------  Convert MATLAB variables to FORTRAN variables  
      
      ISCOMPLEX = MXISCOMPLEX(PRHS(1))       
      IF(ISCOMPLEX.EQ.1) THEN
         CALL MXCOPYPTRTOREAL8(MXGETPI(PRHS(1)),OMEGA,1)          
      ELSE 
         OMEGA = 0
      END IF        
      CALL MXCOPYPTRTOREAL8(MXGETPR(PRHS(2)), SIG2 ,1)
      CALL MXCOPYPTRTOREAL8(MXGETPR(PRHS(3)), T1 ,ROWS)
      
      GAMMA = DCMPLX(ALPHA,OMEGA)

C ------ Call computational subroutine     

      CALL COMPUP(GAMMA,SIG2, T1, ROWS, RUPSI, IUPSI )   

C ------  Convert FORTRAN variables to MATLAB variables  

      CALL MXCOPYREAL8TOPTR(RUPSI,MXGETPR(PLHS(1)),ROWS)
      CALL MXCOPYREAL8TOPTR(IUPSI,MXGETPI(PLHS(1)),ROWS)

      RETURN 
      END

C ----- Computational subroutine
      SUBROUTINE COMPUP(GAMMA, SIG2, T1, ROWS, 
     +   RUPSI, IUPSI)        
	  INTEGER I,J, ROWS
      DOUBLE PRECISION RUPSI(ROWS), IUPSI(ROWS)
	  DOUBLE PRECISION SIG2, T1(ROWS), DEV, REALZ1, REALZ2 
	  COMPLEX *16 GAMMA, WOFZ1, WOFZ2
      COMPLEX *16 Z1, Z2, UPSI
	  
      
      DEV = DSQRT(SIG2)
      Z2 = DEV*GAMMA/2      
      REALZ2 = DREAL(Z2)

      IF(REALZ2 .GE.0) THEN    
         CALL WOFZ(DCMPLX(0,1)*Z2, WOFZ2)
      ELSE
         CALL WOFZ(DCMPLX(0,-1)*Z2, WOFZ2)
      END IF  

      DO 20, I=1,ROWS                
             Z1 = T1(I)/DEV - DEV*GAMMA/2   
             REALZ1 = DREAL(Z1)
             IF(REALZ1.GE.0) THEN    
                CALL WOFZ(DCMPLX(0,1)*Z1, WOFZ1)
             ELSE
                CALL WOFZ(DCMPLX(0,-1)*Z1, WOFZ1)
             END IF
C		real(z1) >= 0 and real(z2)>=0
             IF(REALZ1. GE. 0. AND .REALZ2. GE .0  ) THEN
                UPSI = 2*EXP(SIG2*(GAMMA**2)/4 - 
     +              GAMMA*(T1(I))) - EXP(-(T1(I))**2/SIG2
     +         + LOG(WOFZ1)) - EXP( - GAMMA*T1(I) 
     +         + LOG(WOFZ2))
             END IF
C		real(z1) < 0 and real(z2)>=0
		   IF(REALZ1. LT. 0. AND .REALZ2. GE .0  ) THEN
                UPSI =  EXP(-(T1(I))**2/SIG2
     +         + LOG(WOFZ1)) - EXP(- GAMMA*T1(I) 
     +         + LOG(WOFZ2))
             END IF
C		real(z1) >= 0 and real(z2)<0
		   IF(REALZ1. GE. 0. AND .REALZ2. LT .0  ) THEN
                UPSI = -  EXP(-(T1(I))**2/SIG2
     +         + LOG(WOFZ1)) + EXP(- GAMMA*T1(I) 
     +         + LOG(WOFZ2))
             END IF
C		real(z1) < 0 and real(z2)<0             
		   IF(REALZ1. LT. 0. AND .REALZ2. LT .0  ) THEN
                UPSI = - 2*EXP(SIG2*(GAMMA**2)/4 - 
     +              GAMMA*(T1(I))) + EXP(-(T1(I))**2/SIG2
     +         + LOG(WOFZ1)) + EXP( - GAMMA*T1(I) 
     +         + LOG(WOFZ2))
             END IF
           RUPSI(I) = DREAL(UPSI) 
           IUPSI(I) = DIMAG(UPSI)
C            RUPSI(I,J) = T1(1)
C            IUPSI(I,J) = T2(1)
10      CONTINUE
20    CONTINUE
      RETURN  
      END
      
C-----  Complex error function due to Hui et al.

      SUBROUTINE WOFZ(Z, W)       
      DOUBLE PRECISION A0, A1, A2, A3, A4, A5, A6
      DOUBLE PRECISION B0, B1, B2, B3, B4, B5, B6
      COMPLEX*16  Z, T, W
      DOUBLE PRECISION X,Y
      A0 = 122.607931777104326
      A1 = 214.382388694706425 
      A2 = 181.928533092181549 
      A3 = 93.155580458138441 
      A4 = 30.180142196210589 
      A5 = 5.912626209773153 
      A6 = 0.564189583562615 
      B0 = 122.607931773875350
      B1 = 352.730625110963558
      B2 = 457.334478783897737
      B3 = 348.703917719495792
      B4 = 170.354001821091472
      B5 = 53.992906912940207
      B6 = 10.479857114260399
      X = DREAL(Z)
      Y = DIMAG(Z)  
      T = DCMPLX(Y, -X)
      W  = ((((((A6*T + A5)*T + A4)*T + A3)*T + A2)*T + A1)*T+A0)
     +    /(((((((T + B6)*T + B5)*T + B4)*T + B3)*T + B2)*T + B1)*T+B0)
      RETURN
      END
