C
C RD.F - Gateway function for RD
C
C This subroutine is the main gateway to MATLAB.  When a MEX function
C  is executed MATLAB calls the MEXFUNCTION subroutine in the corresponding
C  MEX file.  

      SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)

C-----------------------------------------------------------------------
C     (pointer) Replace integer by integer*8 on the DEC Alpha
C     64-bit platform
C
      INTEGER PLHS(*), PRHS(*)

C-----------------------------------------------------------------------
C

      INTEGER NLHS, NRHS
C
C-----------------------------------------------------------------------
C     (pointer) Replace integer by integer*8 on the DEC Alpha
C     64-bit platform
C
      INTEGER MXCREATEDOUBLEMATRIX, MXGETPR

C-----------------------------------------------------------------------
C

      INTEGER MXGETM, MXGETN
C
C KEEP THE ABOVE SUBROUTINE, ARGUMENT, AND FUNCTION DECLARATIONS FOR USE
C IN ALL YOUR FORTRAN MEX FILES.
C---------------------------------------------------------------------
C
C-----------------------------------------------------------------------
C     (pointer) Replace integer by integer*8 on the DEC Alpha
C     64-bit platform
C
      INTEGER RDlhs, X, Y, Z
      REAL*8 RRDlhs, RX, RY, RZ

C
C CHECK FOR PROPER NUMBER OF ARGUMENTS
C
      IF (NRHS .NE. 3) THEN
        CALL MEXERRMSGTXT('RD requires three input arguments')
      ELSEIF (NLHS .GT. 1) THEN
        CALL MEXERRMSGTXT('RD requires one output argument')
      ENDIF
C
C CHECK THE DIMENSIONS OF Y.  IT CAN BE 4 X 1 OR 1 X 4.
C
C      M = MXGETM(PRHS(2))
C      N = MXGETN(PRHS(2))
C
C      IF ((MAX(M,N) .NE. 4) .OR. (MIN(M,N) .NE. 1)) THEN
C        CALL MEXERRMSGTXT('YPRIME requires that Y be a 4 x 1 vector')
C      ENDIF
C
C CREATE A MATRIX FOR RETURN ARGUMENT
C
      PLHS(1) = MXCREATEDOUBLEMATRIX(1,1,0)
C
C ASSIGN POINTERS TO THE VARIOUS PARAMETERS
C
      RDlhs = MXGETPR(PLHS(1))
C
      X = MXGETPR(PRHS(1))
      Y = MXGETPR(PRHS(2))
      Z = MXGETPR(PRHS(3))
C
C COPY RIGHT HAND ARGUMENTS TO LOCAL ARRAYS OR VARIABLES
      CALL MXCOPYPTRTOREAL8(X, RX, 1)
      CALL MXCOPYPTRTOREAL8(Y, RY, 1)
      CALL MXCOPYPTRTOREAL8(Z, RZ, 1)
C
C DO THE ACTUAL COMPUTATIONS IN A SUBROUTINE
C       CREATED ARRAYS.  
C
      CALL RD2(RRDlhs,RX,RY,RZ)
C
C COPY OUTPUT WHICH IS STORED IN LOCAL ARRAY TO MATRIX OUTPUT
      CALL MXCOPYREAL8TOPTR(RRDlhs, RDlhs, 1)
C
      RETURN
      END
     
C
C
C
	SUBROUTINE RD2(RD,X,Y,Z)
CC COMPUTES CARLSON'S ELLIPTIC INTEGRAL OF THE SECOND KIND
CC X AND Y MUST BE NON-NEGATIVE AND ONLY ONE MAY BE ZERO
CC Z MUST BE POSITIVE
CC
CC  TAKEN DIRECTLY FROM NUMERICAL RECIPIES IN FORTRAN, CHAPTER 6 PAGE 257.
CC  MODIFIED IN THAT ITS DOUBLE PRECISION
	DOUBLE PRECISION RD, X, Y, Z, ERRTOL, TINY, BIG
	DOUBLE PRECISION C1, C2, C3, C4, C5, C6
	DOUBLE PRECISION ALAMB, AVE, DELX, DELY, DELZ, EA, EB, EC, ED
	DOUBLE PRECISION EE, FAC, SQRTX, SQRTY, SQRTZ, SUM, XT, YT, ZT
	PARAMETER (ERRTOL = 0.0015)
	PARAMETER (TINY = 1.D-25)
	PARAMETER (BIG = 4.5D21)
	PARAMETER (C1 = 3./14.)
	PARAMETER (C2 = 1./6.)
	PARAMETER (C3 = 9./22.)	
	PARAMETER (C4 = 3./26.)
	PARAMETER (C5 = 9./88.)
	PARAMETER (C6 = 9./52.)
CC
	IF ( MIN(X,Y) .LT. 0. .OR. MIN(X+Y,Z) .LT. TINY 
     &		.OR. MAX(X,Y,Z) .GT. BIG ) THEN
     		WRITE(*,*) 'INVALID ARGUMENTS IN RD'
     		WRITE(*,*) 'X = ',X
     		WRITE(*,*) 'Y = ',Y
     		WRITE(*,*) 'Z = ',Z
     		RD = -1.0
     		RETURN
     	END IF
CC
	XT = X
	YT = Y
	ZT = Z
	SUM = 0.
	FAC = 1.
 1	CONTINUE
 		SQRTX = SQRT(XT)
 		SQRTY = SQRT(YT)
 		SQRTZ = SQRT(ZT)
 		ALAMB = SQRTX * (SQRTY + SQRTZ) + SQRTY * SQRTZ
 		SUM = SUM + FAC / (SQRTZ * (ZT + ALAMB))
 		FAC = FAC / 4.
 		XT = (XT + ALAMB) / 4.
 		YT = (YT + ALAMB) / 4.
 		ZT = (ZT + ALAMB) / 4.
 		AVE = (XT + YT + 3.*ZT) / 5.
 		DELX = (AVE - XT) / AVE
 		DELY = (AVE - YT) / AVE
 		DELZ = (AVE - ZT) / AVE
 	IF ( MAX( ABS(DELX), ABS(DELY), ABS(DELZ)) .GT. ERRTOL ) GOTO 1
 	EA = DELX * DELY
 	EB = DELZ * DELZ
 	EC = EA - EB
 	ED = EA - 6. * EB
 	EE = ED + EC + EC
 	RD = 3. * SUM + FAC * 
     &		(1. + ED * (-C1 + C5 * ED - C6 * DELZ * EE)
     &		+ DELZ * (C2 * EE + DELZ * (-C3 * EC + DELZ * C4 * EA)))
     & 	/ (AVE * SQRT(AVE))
      	RETURN
      	END