CC---PROGRAM TO PROVIDE EXACT SOLUTIONS
      SUBROUTINE ABQMAIN
C
      IMPLICIT REAL*8(A-H,O-Z)
C*****MAIN PROGRAM
      DIMENSION P2(27),P3(27),DSIG1(27),DSIG2(27)
C
      DATA P2/30.,32.5,35.,37.5,40.,37.5,35.,32.5,
     1 30.,27.5,25.,22.5,20.,17.5,15.,12.5,10.,
     2 7.5,5.,2.5,0.,-2.5,-5.,-7.5,-10.,-12.5,-15./
      DATA P3/0.,0.,0.,0.,0.,2.5,5.,7.5,10.,12.5,
     1 15.,17.5,20.,22.5,25.,27.5,30.,32.5,35.,37.5,40.,42.5,45.,
     2 47.5,50.,52.5,55./
C*****P2 CONTAINS PRESCRIBED SIGMAX VALUE
C*****P3 CONTAINS PRESCRIBED SIGMAY VALUE
C*****INITIALIZE VARIABLES
      SIG0=3.E4
      SIG1=0.
      SIG2=0.
      ALF1=0.
      ALF2=0.
      E1=0.
      E2=0.
      EP1=0.
      EP2=0.
C*****COMPUTE STRESS INCREMENTS
      DSIG1(1)=30000.
      DSIG2(1)=0.
      DO 200 J=1,26
         DSIG1(J+1)=(P2(J+1)-P2(J))*1000.
         DSIG2(J+1)=(P3(J+1)-P3(J))*1000.
  200 CONTINUE
C*****DSIG1 CONTAINS PRESCRIBED X-STRESS INCREMENTS
C*****DSIG2 CONTAINS PRESCRIBED Y-STRESS INCREMENTS
      SUB=1./1000.
      DO 300 J=1,27
         DS1=DSIG1(J)*SUB
         DS2=DSIG2(J)*SUB
C*****DIVIDE EACH STEP INTO 1000 SUBINCREMENTS
         DO 400 I=1,1000
            BETA=0.
            SIG1=SIG1+DS1
            SIG2=SIG2+DS2
C*****TEST IF YIELDING
            TEST=(SIG1-ALF1)**2+(SIG2-ALF2)**2-(SIG1-ALF1)*(SIG2-ALF2)
            TEST=SQRT(TEST)
C*****BETA IS ZERO IF ELASTIC,  ONE IF PLASTIC
            IF(TEST.GT.1.0001*SIG0)BETA=1.
            SIG1=SIG1-DS1
            SIG2=SIG2-DS2
C*****CALL KINEMATIC HARDENING SUBROUTINE TO UPDATE VARIABLES
            CALL KINPLA(SIG1,SIG2,ALF1,ALF2,DS1,DS2,E1,E2,EP1,EP2,
     Q                  SIG0,BETA)
  400    CONTINUE
C*****WRITE OUT RESULTS AT END OF EACH MAJOR STRESS INCREMENT
         WRITE(6,9999)SIG1,SIG2,E1,E2,EP1,EP2,SIG0,BETA
 9999    FORMAT(1P8E14.4)
  300 CONTINUE
      STOP
      END
C
      SUBROUTINE KINPLA(SIG1,SIG2,ALF1,ALF2,DSIG1,DSIG2,E1,E2,EP1,EP2,
     1 SIG0,BETA)
      IMPLICIT REAL*8(A-H,O-Z)
C*****ON ENTRY TO THE SUBROUTINE THE VARIABLES SIG1,SIG2,ALF1,ALF2,
C*****E1,E2,EP1,EP2 ARE THE VALUES OF STRESS, YIELD SURFACE CENTRE,
C*****TOTAL STRAINS AND PLASTIC STRAINS AT THE BEGINNING OF THE
C*****INCREMENT.  ON EXIT FROM THE SUBROUTINE THESE VARIABLES CONTAIN
C*****THE VALUES AT THE END OF THE INCREMENT.  DSIG1 AND DSIG2 ARE THE
C*****PRESCRIBED STRESS INCREMENTS
      REAL*8 N1,N2
C*****AC IS THE STRESS VS. PLASTIC STRAIN SLOPE, E IS YOUNGS MODULUS,
C*****ANU IS POISSONS RATIO
      AC=1.589473684E6
      ANU=0.3
      E=30.E6
      AD=E/(1.-ANU*ANU)
      BD=ANU*AD
      N1=(SIG1-.5*SIG2-ALF1+.5*ALF2)/SIG0
      N2=(SIG2-0.5*SIG1-ALF2+.5*ALF1)/SIG0
      ANDN=AD*N1*N1+2.*BD*N1*N2+AD*N2*N2
      DENOM=ANDN+AC
      P=AD*N1+BD*N2
      Q=BD*N1+AD*N2
      A=AD-(P*P/DENOM)*BETA
      B=BD-(P*Q/DENOM)*BETA
      C=AD-(Q*Q/DENOM)*BETA
      DET=B*B-A*C
      R=-C/DET
      S=B/DET
      T=-A/DET
      DE1=R*DSIG1+S*DSIG2
      DE2=S*DSIG1+T*DSIG2
      DEP=(P*DE1+Q*DE2)/DENOM
      IF(BETA.EQ.0.)DEP=0.
      DEP1=DEP*N1
      DEP2=DEP*N2
      FAC1=AC*DEP/SIG0
      DALF1=(SIG1-ALF1)*FAC1
      DALF2=(SIG2-ALF2)*FAC1
      ANDEP=ALF1*DEP1+ALF2*DEP2
      IF(ANDEP.LT.0..AND.SIG0.EQ.3.E+4)WRITE(6,999)SIG1,SIG2
  999 FORMAT(10H CHANGE AT,1P2E15.4)
C*****TEST FOR ORNL YIELD SURFACE EXPANSION
      IF(ANDEP.LT.0.)SIG0=3.4E4
      SIG1=SIG1+DSIG1
      SIG2=SIG2+DSIG2
      ALF1=ALF1+DALF1
      ALF2=ALF2+DALF2
      E1=E1+DE1
      E2=E2+DE2
      EP1=EP1+DEP1
      EP2=EP2+DEP2
      RETURN
      END