PROGRAM TRAP
C     Trapezoidal rule solution to one-dimensional shock problem
C     using an analytical form of the doubly-asymptotic approximation. 
C     This program uses DAA2; to get a DAA1 solution
C     simply set daa2=0.0
C
C     Gamma is a constant used in some DAA formulations to offset
C     the computed pressure field in the fluid to account for the
C     shock arrival. This is set to 0 here.
C
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 Ms,Mf,Ks, M(2,2), C(2,2), K(2,2), F(2), U(2), UDOT(2),
     1     U2DOT(2), A(2,2), AINV(2,2), RHS(2), VEC1(2), VEC2(2),
     2     UOLD(2),UDOTOLD(2)
      DATA M/4*0.0/, C/4*0.0/, K/4*0.0/, F/2*0.0/,
     1     U/2*0.0/, UDOT/2*0.0/, U2DOT/2*0.0/
C
      daa2 = 1.0
      Pi = 1.0
      Ms = 5.0
      Ks = 197.39209
      Mf = 0.5611
      Af = 1.0
      G  = 1.0    
      rhocee = 1.0  
      Gamma = 0.0
      Ds  = rhocee*Af*G*G*Af/Ms
      Df1 = rhocee*Af*Af/Mf
      Df2 = rhocee*Df1*Af/Mf
      H1  = Ds - (Ds + (1.0-daa2)*Df1)*Gamma
      H2  = daa2*Df2*Gamma
      dt = 0.025
      dtsq = dt**2
      NINC = 400
C
      M(1,1) = Ms
      M(2,2) = Af
      C(1,2) = G*Af
      C(2,1) = -daa2*rhocee*Df1*G
      C(2,2) = Df1 + Ds
      K(1,1) = Ks
      K(2,1) = rhocee*Af*G*Ks/Ms
      K(2,2) = daa2*Df2
      F(1)   = -G*Af*(1.0-Gamma)*Pi
      U2DOT(1) = F(1)/M(1,1)
C
      DO I=1,2
        DO J=1,2
          A(I,J) = M(I,J)*4.0/dtsq + C(I,J)*2.0/dt + K(I,J)
        END DO
      END DO
      DET = A(1,1)*A(2,2) - A(1,2)*A(2,1)
      AINV(1,1) =  A(2,2)/DET
      AINV(1,2) = -A(1,2)/DET
      AINV(2,1) = -A(2,1)/DET
      AINV(2,2) =  A(1,1)/DET
      WRITE(6,*) 'Ks = ',Ks,'   Ms = ',Ms
      WRITE(6,*) 'daa2 = ',daa2
      WRITE(6,*) 'DET = ',DET              
      WRITE(6,*) 'A matrix ', A(1,1),A(1,2)
      WRITE(6,*) '         ', A(2,1),A(2,2)
      DISCRIM = SQRT((AINV(1,1)-AINV(2,2))**2 + AINV(1,2)*AINV(2,1))
      EIG1 = (AINV(1,1) + AINV(2,2) + DISCRIM)/2.0
      EIG2 = (AINV(1,1) + AINV(2,2) - DISCRIM)/2.0
      WRITE(6,*) 'AINV EIGS = ', EIG1,EIG2
C
      T = 0.0
      Ptot = (1.0-Gamma)*Pi + UDOT(2)
      WRITE(6,100)
  100 FORMAT('   T          X             dX/dt       ',
     1 '  d2X/dt2          Ptot            qM            dqM/dt',
     2 '       d2qM/dt2'/
     3       ' ---------------------------------------',
     4 '-------------------------------------------------------',
     5 '---------------')  
      WRITE(6,200) T, U(1),UDOT(1),U2DOT(1), Ptot,
     1                U(2),UDOT(2),U2DOT(2)
  200 FORMAT(1X,F6.3,', ',6(G13.5,', '),G13.5)
C
      DO KINC=1,NINC
         T = T + dt
         F(2) = -H1*Pi + H2*Pi*T
         DO I=1,2
           VEC1(I) = U(I)*2.0/dt + UDOT(I)
           VEC2(I) = U(I)*4.0/dtsq + UDOT(I)*4.0/dt + U2DOT(I)
           RHS(I)  = F(I)
         END DO
         DO I=1,2
           DO J=1,2
             RHS(I) = RHS(I) + C(I,J)*VEC1(J) + M(I,J)*VEC2(J)
           END DO
         END DO
         DO I=1,2
           UOLD(I) = U(I)
           UDOTOLD(I) = UDOT(I)
           U(I) = 0.0
           DO J=1,2
             U(I) = U(I) + AINV(I,J)*RHS(J)
           END DO
         END DO
         DO I=1,2
           UDOT(I)  = (U(I)-UOLD(I))*2.0/dt - UDOTOLD(I)
           U2DOT(I) = (UDOT(I)-UDOTOLD(I))*2.0/dt - U2DOT(I)
         END DO
         Ptot = (1.0-Gamma)*Pi + UDOT(2)
        WRITE(6,200) T, U(1),UDOT(1),U2DOT(1), Ptot,
     1                  U(2),UDOT(2),U2DOT(2)
      END DO
      STOP
      END