PROGRAM trntr3
C  Comment lines found in trntr1 are not repeated here, and should be read
C  first.
C ------------------------ WARNING ---------------------------
C This program uses quadruple precision (REAL*16) which is only supported
C by a few compilers.  It was the program actually used to obtain the 
C data files.
C
      REAL*16 P(0:2307277),V(0:2307277),Q(0:2307277),FAC1(151,151),XN
      REAL*16 FCTR1,FCTR2,FCTR3,FCTR4,FCTRQ,FCTRV,FACTOR
      REAL*16 PSUM,VSUM,QSUM,FAC,XS
      INTEGER MMAX,MAXS,NM,N,IS,M,NQ,IQ,INDX1,JS,IP,KS
      COMMON/ ARRAYS/ P,V,Q,FAC1
C This COMMON block reduces disk space requirements with many compilers
C  and helps memory management on some computers.
      MAXS=300
      MMAX=MAXS+1
C Tabulate FAC1.
      DO 10 N=1,(MMAX+1)/2
        XN=qext(N)
        FAC1(N,1)=1q0
        DO 10 IS=2,(MMAX+1)/2
          XS=qext(IS)
10    FAC1(N,IS)=FAC1(N,IS-1)*(XN+XS)/(XS-1q0)
C We first set the initial conditions for translation.
      P(INDX1(1,0,0))=1q0
      V(INDX1(1,0,0))=1q0
      Q(INDX1(1,0,0))=0q0
      DO 500 M=1,MMAX
       NM=M/2+1
       IF(M.EQ.MMAX)NM=2
       DO 400 N=1,NM
        NQ=M-2*N+2
        XN=qext(N)
        FCTR1=(XN+.5q0)/(XN+1q0)
        FCTR2=XN*(XN-.5q0)/(XN+1q0)
        FCTR3=XN*(2q0*XN*XN-.5q0)/(XN+1q0)
        FCTR4=(8q0*XN*XN-2q0)/(3q0*XN+3q0)
        FCTRV=2q0*XN/((XN+1q0)*(2q0*XN+3q0))
        FCTRQ=3q0/(2q0*XN*(XN+1q0))
        DO 300 IQ=0,NQ
         IP=M-IQ-N+1
         JS=(IQ+1)/2
         KS=IQ/2
         PSUM=0q0
         VSUM=0q0
         QSUM=0q0
C If JS is less than 1, we skip the loop on IS and set the elements to 0
         IF (JS.LT.1) GO TO 250
         DO 200 IS=1,JS
          XS=qext(IS)
          FAC=FAC1(N,IS)
C
          FACTOR=2q0*(XN*XS+1q0)**2/(XN+XS)
          PSUM=PSUM+FAC*FCTR1*(4q0+XN*XS-FACTOR)
     &     *P(INDX1(IS,IQ-IS,IP-N+1))/(XS*(2q0*XS-1q0))
          IF(IP.GE.N) THEN
            QSUM=QSUM-FAC*FCTRQ*P(INDX1(IS,IQ-IS,IP-N))/XS
          ENDIF
          IF(IP.GT.N) THEN
            PSUM=PSUM+FAC*FCTR2*P(INDX1(IS,IQ-IS,IP-N-1))
            VSUM=VSUM+FAC*FCTRV*P(INDX1(IS,IQ-IS,IP-N-1))
          ENDIF
          IF(IS.LE.KS) THEN
            PSUM=PSUM-FAC*FCTR4*Q(INDX1(IS,IQ-IS-1,IP-N+1))
            IF(IP.GE.N) THEN
              QSUM=QSUM+FAC*XS*Q(INDX1(IS,IQ-IS-1,IP-N))/(XN+1q0)
             ENDIF
          ENDIF
          IF(IS.LT.JS) THEN
            PSUM=PSUM+FAC*FCTR3*V(INDX1(IS,IQ-IS-2,IP-N+1))/(2q0*XS+1q0)
          ENDIF
  200    CONTINUE
  250    P(INDX1(N,IP,IQ))=PSUM
         V(INDX1(N,IP,IQ))=PSUM+VSUM
  300   Q(INDX1(N,IP,IQ))=QSUM
  400  CONTINUE
  500 CONTINUE
C
C The full set of coefficients having now been calculated, I collect the
C  ones I wish to keep in the fronts of the arrays ready for writing.
C  I also enforce non-dimensionalizations and other consequences
C  of the definitions of the resistance functions.
C
C I map the RYA coefficients using: P(INDX1(1,m-q,q)) -> V( 1+m(m+1)/2+q ).
C I map the RYB ones using: 2Q(INDX1(1,m-q,q)) -> Q( 1+m(m+1)/2+q ).
C
      DO 600 M=0,MAXS
        NM=(M*(M+1))/2
        DO 600 IQ=0,M
          Q(NM+IQ)=2q0*Q(INDX1(1,M-IQ,IQ))
          V(NM+IQ)=P(INDX1(1,M-IQ,IQ))
  600   CONTINUE  
C
C I map the RYG ones using: (3/4) P(INDX1(2,m-q,q)) -> P( 1+m(m+1)/2+q ).
C  The term P(INDX1(2,0,1)) was not stored and is set separately.
C
      P(0)=0q0
      DO 620 M=1,MAXS
        NM = (M*(M+1))/2
        DO 610 IQ=0,(M-1)
  610   P(NM+IQ)=0.75q0*P( INDX1(2,M-IQ,IQ) )
  620 P(NM+M) = 0q0
C
C Now write the output files
C
      MMAX=((MAXS+1)*(MAXS+2))/2 -1
700   FORMAT(D22.16)
      OPEN (1,FILE='ryacof.dat',STATUS='NEW')
      WRITE(1,710)MAXS
710   FORMAT(' RYA',I5)
      DO 720 M=0,MMAX
720   WRITE (1,700) V(M)
      CLOSE (1)
C
      OPEN (2,FILE='rybcof.dat',STATUS='NEW')
      WRITE(2,740)MAXS
 740  FORMAT(' RYB',I5)
      DO 750 M=0,MMAX
750   WRITE (2,700) Q(M)
      CLOSE (2)
C
      OPEN (3,FILE='rygcof.dat',STATUS='NEW')
      WRITE(3,760)MAXS
 760  FORMAT(' RYG',I5)
      DO 770 M=0,MMAX
770   WRITE (3,700) P(M)
      CLOSE (3)
C
      STOP
      END
C
      INTEGER FUNCTION INDX1(N,IP,IQ)
C ---------------------------- Arguments --------------------------- 
      INTEGER N,IP,IQ
C ------------------------- Local Variables ------------------------
      INTEGER M,M2,LM
      M=IP+IQ+N-1
      M2=M/2
      LM=( M2*(M2+1) )/2
      LM=LM*(M2+1)+ ( LM*(M2+2) )/3
C LM is calculated as above to avoid overflow.
C If MMAX=200, then LM = 681750
      LM=LM+(M-2*M2)*(M2+1)*(M2+1)
      INDX1=LM+(N-1)*(M-N+3)+IQ
      RETURN
      END