PROGRAM rax 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 as set in the DATA statements. C C The program writes its output to FORTRAN-77 i/o unit *, which C is standard output, typically a console. For output to a file or printer C the unit can be changed, or the output re-directed by the operating system. C DOUBLE PRECISION DLAMDA(9),RAX,RAX11,RAX12,RAX21,RAX22,W1 INTEGER J DATA DLAMDA/.5D0,1D0,2D0,3D0,4D0,5D0,10D0,20D0,100D0/ WRITE(6,10) 10 FORMAT(1X, & ' LAMBDA RAX11 RAX12 RAX21 RAX22 W1'/) DO 100 J=1,9 RAX11=RAX(1,1,DLAMDA(J)) RAX12=RAX(1,2,DLAMDA(J)) RAX21=RAX(2,1,DLAMDA(J)) RAX22=RAX(2,2,DLAMDA(J)) W1 = RAX11 +(1D0+DLAMDA(J))*RAX12 + DLAMDA(J)*RAX22 WRITE(*,50)DLAMDA(J),RAX11,RAX12,RAX21,RAX22,W1 50 FORMAT(1X,F7.3, 5(F11.5) ) 100 CONTINUE STOP END C DOUBLE PRECISION FUNCTION RAX(IALPHA,IBETA,DLAMDA) C ----------------------------- Arguments of function ----------------- INTEGER IALPHA, IBETA DOUBLE PRECISION DLAMDA C --------------------------------------------------------------------- C X C This subprogram calculates the functions A (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.22)-(3.23), although the form in which the C singularities are expanded follows Jeffrey & Corless (1988) C PhysicoChemical Hydrodynamics, Vol 10, 461-470. C The subprogram works generally the same as rxa.f where further comments C can be found. C C In the program the coefficients are read in only once and checked C to make certain the signature is correct (cf axitr1). C The coefficients are expected in rxa300.dat in the 95 format. C Thus MAXS = 200, RXACOF(0:20300) C MAXS = 250, RXACOF(0:31625) C MAXS = 300, RXACOF(0:45450) C C ------------------------ Local variables ---------------------------- DOUBLE PRECISION RXACOF(0:45450), G1, G2, G3, G1SUM, G2SUM, G3SUM DOUBLE PRECISION ZL, COEF, 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 IF(I.EQ.J) THEN RAX=1D0 - G1/4D0 - G3 MMIN=2 MMAX=MAXS ELSE RAX=G1/4D0+G2*DLOG(4D0)-2D0*G3 MMIN=1 MMAX=MAXS-1 ENDIF 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 ZL=1D0 COEF=0D0 C Some computers will underflow in 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) RAX=RAX+COEF-G1SUM-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M) 100 CONTINUE IF(I.EQ.J)RETURN RAX=-RAX*2D0/XLP1 IF(I.EQ.2)RAX=RAX*XL RETURN END