PROGRAM rxa DEM C This program demonstrates the use of the coefficients calculated C by the programs axitrx. The program will produce a table of values C for different values of LAMBDA and S as set in the DATA statements. C DOUBLE PRECISION DLAMDA(5),S(4),RXA,RXA11,RXA12,RXA21,RXA22 INTEGER I,J DATA DLAMDA/1D0,.5D0,.25D0,.125D0,2D0/ DATA S/2.01D0,2.1D0,2.5D0,3D0/ WRITE(*,20) 20 FORMAT(1X,'TABLE OF VALUES FOR RESISTANCE FUNCTIONS XA'/ & 1X,'CHECK PROGRAM LISTING FOR NUMBER OF TERMS USED') DO 100 J=1,5 WRITE(*,30)DLAMDA(J) 30 FORMAT(/1X,'LAMBDA =',F8.3) WRITE(*,40) 40 FORMAT(/7X,'S',10X, & 'RXA(1,1)',8X,'RXA(1,2)',8X,'RXA(2,1)',8X,'RXA(2,2)') DO 100 I=1,4 RXA11=RXA(1,1,S(I),DLAMDA(J)) RXA12=RXA(1,2,S(I),DLAMDA(J)) RXA21=RXA(2,1,S(I),DLAMDA(J)) RXA22=RXA(2,2,S(I),DLAMDA(J)) WRITE(*,50)S(I),RXA11,RXA12,RXA21,RXA22 50 FORMAT(F10.3,4F16.5) 100 CONTINUE STOP END C DOUBLE PRECISION FUNCTION RXA(IALPHA,IBETA,S,DLAMDA) C ---------------------- Arguments of function ------------------------ INTEGER IALPHA, IBETA DOUBLE PRECISION S, DLAMDA C --------------------------------------------------------------------- C A C This routine calculates the functions X (s,lambda) C ALPHA BETA C defined in Jeffrey & Onishi (1984) J. Fluid Mech. Vol 139, 261-290. C Coefficients calculated by one of the programs axitrx are summed using C essentially equations (3.20)-(3.21), although the form in which the C singularities are expanded follows Jeffrey & Corless (1988) C PhysicoChemical Hydrodynamics, Vol 10, 461-470. C The subprogram uses Jeffrey & Onishi equation (1.9a) to convert C any call to one in which lambda.LT.1. It goes forwards or backwards C through the coefficients ( ISTEP) depending upon whether I=1 or 2. C This technique is not explained in the paper. C Starting at equn (1.9) we see that X can be obtained from X C 22 11 C by inverting lambda. Now looking at equn (3.15) we substitute 1/lambda C and redefine q as k-q. Factoring the lambda**k, we see that lambda**q C now multiplies P(1,q,k-q), which means running 'backwards' through the C coefficients relative to the published case. C C The program reads the coefficients from rxaxxx.dat only once and checks C to make certain the signature is correct (cf axitr2). C The xxx refers to the maximum number of terms. The array RXACOF C must be dimensioned to be greater that the total number of terms. C RXACOF(0:((MAXS+1)*(MAXS+2))/2 -1) C Thus MAXS = 200, RXACOF(0:20300) C MAXS = 250, RXACOF(0:31625) C MAXS = 300, RXACOF(0:45450) C ------------------------ Local variables ---------------------------- DOUBLE PRECISION RXACOF(0:45450), G1, G2, G3, G1SUM, G2SUM, G3SUM DOUBLE PRECISION T, ZT, ZL, COEF, EPS, XL, XLP1, XL13 INTEGER I, J, ISTEP, IQ, MMIN, MMAX, M, NM, MAXS LOGICAL LOADED CHARACTER*4 ID DATA LOADED/.FALSE./ C ------------------------ Start of executable code ------------------- IF(.NOT.LOADED) THEN OPEN(1,FILE='rxa300.dat',STATUS='OLD') READ(1,10) ID,MAXS 10 FORMAT(A4,I5) IF(ID.NE.' RXA') THEN WRITE(*,20) 20 FORMAT(1X,'WRONG COEFFICIENTS FOR RXA') STOP ENDIF READ(1,30) RXACOF 30 FORMAT(D22.16) CLOSE(1) LOADED=.TRUE. ENDIF IF(DLAMDA.GT.1D0) THEN I=3-IALPHA J=3-IBETA XL=1D0/DLAMDA ELSE I=IALPHA J=IBETA XL=DLAMDA ENDIF XLP1=XL+1D0 XL13=XLP1**3 IF(I.EQ.1) THEN C These are the g_i for X^A_11 and X^A_12 (i=1..3 in paper) G1=2D0*XL**2/XL13 G2=(XL+7D0*XL**2+XL**3)/(5D0*XL13) G3=(1D0+18D0*XL-29D0*XL**2+18D0*XL**3+XL**4)/(42D0*XL13) ISTEP=1 ELSE G1=2D0*XL/XL13 G2=(1D0+7D0*XL+XL**2)/(5D0*XL13) G3=(1D0+18D0*XL-29D0*XL**2+18D0*XL**3+XL**4)/(42D0*XL*XL13) ISTEP=-1 ENDIF T=2D0/S EPS=0.25D0*S*S-1D0 IF(I.EQ.J) THEN ZT=T*T C Notice that we supply the 1 because we start the summation at 2. RXA=1D0+G1/EPS-(G2+G3*EPS)*DLOG( 1D0-4D0/(S*S) )-G3 MMIN=2 MMAX=MAXS ELSE ZT=T RXA=G1*S/(2D0*EPS)+(G2+G3*EPS)*DLOG((S+2D0)/(S-2D0))-G3*S MMIN=1 MMAX=MAXS-1 ENDIF T=T*T G1SUM=G1 G2SUM=2D0*G2 G3SUM=4D0*G3 DO 100 M=MMIN,MMAX,2 IF (ISTEP .GT. 0) THEN NM = (M*(M+1))/2 ELSE NM = (M*(M+1))/2 + M ENDIF C This loop adds up the powers of lambda to obtain the term f_m in the C paper. Notice that the 2**(-m) is not needed, because the coefficients C in the file already contain that factor. ZL=1D0 COEF=0D0 C Some computers will underflow on the next statements. DO 50 IQ=0,M COEF=COEF+RXACOF(NM)*ZL NM=NM+ISTEP ZL=ZL*XL 50 CONTINUE COEF=COEF*XLP1**(-M) RXA=RXA + (COEF-G1SUM-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M) )*ZT ZT=ZT*T 100 CONTINUE IF(I.EQ.J)RETURN RXA=-RXA*2D0/XLP1 IF(I.EQ.2)RXA=RXA*XL RETURN END