PROGRAM RYH DEM C This program demonstrates the use of the coefficients calculated C by the programs trnrox or trnstx. The program will produce a table of C values for different values of lambda and s as set in the DATA statements. INTEGER MAXS,MAXL PARAMETER (MAXL=6,MAXS=5) DOUBLE PRECISION DLAMDA(MAXL),S(MAXS),RYH,VAL(MAXL) INTEGER I,J DATA DLAMDA/1D0,2D0,3D0,4D0,5D0,0.5D0/ DATA S/2.01D0,2.02D0,2.05D0,2.1D0,2.5D0/ WRITE(6,10) 10 FORMAT(1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION YH11'/ & ' s l=1 l=2 l=3 l=4 l=5') DO 50 I=1,MAXS DO 40 J=1,MAXL 40 VAL(J)=RYH(1,1,S(I),DLAMDA(J)) 50 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL) WRITE(6,110) 110 FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION YH12'/ & ' s l=1 l=2 l=3 l=4 l=5') DO 150 I=1,MAXS DO 140 J=1,MAXL 140 VAL(J)=RYH(1,2,S(I),DLAMDA(J)) 150 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL) WRITE(6,210) 210 FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION YH21'/ & ' s l=1 l=2 l=3 l=4 l=5') DO 250 I=1,MAXS DO 240 J=1,MAXL 240 VAL(J)=RYH(2,1,S(I),DLAMDA(J)) 250 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL) WRITE(6,310) 310 FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION YH22'/ & ' s l=1 l=2 l=3 l=4 l=5') DO 350 I=1,MAXS DO 340 J=1,MAXL 340 VAL(J)=RYH(2,2,S(I),DLAMDA(J)) 350 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL) 1100 FORMAT(F8.3,8F10.4) STOP END C DOUBLE PRECISION FUNCTION RYH(IALPHA,IBETA,S,DLAMDA) C ----------------------- Arguments of function ----------------------- DOUBLE PRECISION S, DLAMDA INTEGER IALPHA, IBETA C --------------------------------------------------------------------- C This calculates the function RYH defined in Jeffrey, 'The calculation C of the low-Reynolds-number resistance functions for two unequal spheres', C Phys. Fluids (1991). C The coefficients generated by one of the programs trnrox or trnstx are C read in from the file ryh300.dat. The coefficients are summed using C equations (36a,36b). The subprogram works by converting any call to one C in which lambda.LT.1. It goes forwards or backwards through the coefficients C (the variable ISTEP) depending upon whether I=1 or 2. C This version was written during February 1991. C ------------------------- Local variables --------------------------- DOUBLE PRECISION RYHCOF(0:45450) DOUBLE PRECISION XL,XLP1,G2,G3,G2SUM,G3SUM DOUBLE PRECISION T,ZT,ZL,COEF,EPS 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(4,FILE='ryh300.dat',STATUS='OLD') READ(4,10) ID, MAXS 10 FORMAT( A4,I5) IF( ID.NE.' RYH') THEN WRITE(*,20) 20 FORMAT( 1X,'WRONG COEFFICIENTS FOR RHY') STOP ENDIF READ(4,2000)RYHCOF 2000 FORMAT(D22.16) CLOSE(4) 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 T=2D0/S EPS=1D0-4D0/(S*S) IF(I.EQ.J) THEN IF(I.EQ.1) THEN G2=(2D0*XL-XL**2)/(10D0*XLP1**2) G3=(16D0-61D0*XL+180D0*XL**2+2D0*XL**3)/(500D0*XLP1**2) ISTEP=1 ELSE G2=(2D0*XL-1D0)/(10D0*XLP1**2) G3=(16D0*XL**3-61D0*XL**2+180D0*XL+2D0)/(500D0*XL*XLP1**2) ISTEP=-1 ENDIF ZT=T*T RYH=-(G2+G3*EPS*S*S/4D0)*DLOG(EPS)-G3 MMIN=2 MMAX=MAXS ELSE IF(I.EQ.1) THEN G2=(XL**2+7D0*XL**3)/(20D0*XLP1**2) G3=(43D0*XL+147D0*XL**2-185D0*XL**3+221D0*XL**4)/ & (1000D0*XLP1**2) ISTEP=1 ELSE G2=(7D0*XL+XL**2)/(20D0 * XL**2 * XLP1**2) G3=(221D0-185D0*XL+147D0*XL**2+43D0*XL**3)/ & (1000D0*XL**2*XLP1**2) ISTEP=-1 ENDIF ZT=T RYH=(G2+G3*EPS*S*S/4D0)*DLOG((S+2D0)/(S-2D0)) RYH=RYH-G3*S MMIN=1 MMAX=MAXS-1 ENDIF T=T*T 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=XLP1**(-M) COEF=0D0 DO 50 IQ=0,M COEF=COEF+RYHCOF(NM)*ZL NM=NM+ISTEP ZL=ZL*XL 50 CONTINUE RYH=RYH+ZT*(COEF-((G2SUM-G3SUM/DBLE(M+2))/DBLE(M))) ZT=ZT*T 100 CONTINUE IF(I.EQ.J) RETURN RYH=RYH*8D0/XLP1**3 IF(I.EQ.2)RYH=RYH*XL**3 RETURN END