C     *************************
      SUBROUTINE ANAL
C     *************************
      IMPLICIT REAL*8(A-H,O-Z), INTEGER (I-N)
      COMMON/MAT/TEMP(6,27),NORDER,MPRINT,IMAT,NMAT,IFITE,NELM,NOP,
     <NLIST(40)
	include 'analc.inc'
      COMMON/INOUT/IIN,IOUT,ISOUT,ISO
      COMMON/LENGTH/TLENG,ALENG(250),ACLENG(2500)
      COMMON/CONST/PI,TWOPI,CRDEG,CMAGEN,CLIGHT,EMASS,ERAD,ECHG
      DIMENSION C(16),ETA(4)
      NSTABX=0
      NSTABY=0
      NSTAB=0
      NWRNCP=0
      DO 561 IM=1,4
      ETA(IM)=-TEMP(IM,6)
      DO 561 JM=1,4
      IC=(IM-1)*4 +JM
      C(IC)=TEMP(JM,IM)
  561 IF(IM.EQ.JM)C(IC)=C(IC)-1.0D0
      CALL DSIMQ (C,ETA,4,KS)
      IF(KS.NE.1) GO TO 562
      WRITE(IOUT,563)
      IF(ISO.NE.0)WRITE(ISOUT,563)
563   FORMAT(/,' SINGULAR SET OF EQUATIONS SENT TO DSIMQ',/)
      NSTAB=1
      RETURN
  562 COSX=.5D0*(TEMP(1,1)+TEMP(2,2))
      IF(DABS( COSX).GE.1.0D0)GOTO 1
      SIN=DSQRT(1- COSX* COSX)
      SIN=DSIGN(SIN,TEMP(1,2))
      RMU0X=DATAN2(SIN, COSX)
      IF(RMU0X.LT.0.0D0)RMU0X=TWOPI+RMU0X
      CBETAX=TEMP(1,2)/SIN
      CALPHX=.5D0*(TEMP(1,1)-TEMP(2,2))/SIN
      RNU0X=RMU0X/TWOPI
      DET1=DABS(TEMP(1,3)*TEMP(2,4)-TEMP(1,4)*TEMP(2,3))
      DET2=DABS(TEMP(3,1)*TEMP(4,2)-TEMP(3,2)*TEMP(4,1))
      EPS=DET1+DET2
      IF(EPS.GT.1.0D-09)NWRNCP=1
      CETAX=ETA(1)
      CETAPX=ETA(2)
      CETAY=ETA(3)
      CETAPY=ETA(4)
      IF (TLENG.EQ.0.0D0)GOTO 10
      COMPF=(TEMP(5,1)*CETAX+TEMP(5,2)*CETAPX+TEMP(5,6))/TLENG
      GOTO 11
   10 WRITE(IOUT,99999)
99999 FORMAT(/,' LENGTH IS ZERO: COMPACTION FACTOR CANNOT BE COMPUTED',
     >/)
      COMPF=0.0D0
   11 IF (NORDER.EQ.1) GO TO 2
      A110=2.0D0*TEMP(1,7)*ETA(1)+TEMP(1,8)*ETA(2)+TEMP(1,9)*ETA(3)
     >+TEMP(1,10)*ETA(4) +TEMP(1,12)
      A214=TEMP(2,8)*ETA(1)+2.0D0*TEMP(2,13)*ETA(2)+TEMP(2,14)*ETA(3)
     >+TEMP(2,15)*ETA(4)+TEMP(2,17)
      A317=TEMP(3,9)*ETA(1)+TEMP(3,14)*ETA(2)+2.0D0*TEMP(3,18)*ETA(3)
     >+TEMP(3,19)*ETA(4) +TEMP(3,21)
      A419=TEMP(4,10)*ETA(1)+TEMP(4,15)*ETA(2)+TEMP(4,19)*ETA(3)
     >+2.0D0*TEMP(4,22)*ETA(4)+TEMP(4,24)
      A114=TEMP(1,8)*ETA(1)+2.0D0*TEMP(1,13)*ETA(2)+TEMP(1,14)*ETA(3)
     >+TEMP(1,15)*ETA(4)+TEMP(1,17)
      A319=TEMP(3,10)*ETA(1)+TEMP(3,15)*ETA(2)+TEMP(3,19)*ETA(3)
     >+2.0D0*TEMP(3,22)*ETA(4)+TEMP(3,24)
      RMU1X=-.5D0*(A110+A214)/SIN
      CHROMX=RMU1X/TWOPI
      BETA1X=(A114-CBETAX* COSX*RMU1X)/SIN
      ALPH1X=(.5D0*(A110-A214)-RMU1X*CALPHX* COSX)/SIN
    2 COSY=.5D0*(TEMP(3,3)+TEMP(4,4))
      IF(DABS(COSY).GE.1.0D0)GOTO 3

      SIN=DSQRT(1-COSY*COSY)
      SIN=DSIGN(SIN,TEMP(3,4))
      RMU0Y=DATAN2(SIN,COSY)
      IF(RMU0Y.LT.0.0D0)RMU0Y=TWOPI+RMU0Y
      CBETAY=TEMP(3,4)/SIN
      CALPHY=.5D0*(TEMP(3,3)-TEMP(4,4))/SIN
      RNU0Y=RMU0Y/TWOPI
      IF (NORDER.EQ.1) GO TO 4
      RMU1Y=-.5D0*(A317+A419)/SIN
      CHROMY=RMU1Y/TWOPI
      BETA1Y=(A319-CBETAY*COSY*RMU1Y)/SIN
      ALPH1Y=(.5D0*(A317-A419)-RMU1Y*CALPHY*COSY)/SIN
    4 continue
c	call dynamic
c	call fitpar
	RETURN
    1 NSTABX=1
      IF(COSX.EQ.1) GOTO 5
      ALX1 = COSX + DSQRT(COSX*COSX-1)
      ALX2 = COSX - DSQRT(COSX*COSX-1)
      A12=TEMP(1,2)
      A11=TEMP(1,1)
      VXP1=A11-ALX1
      DENOM=DSQRT(A12*A12+VXP1*VXP1)
      VX1=-A12/DENOM
      VXP1=VXP1/DENOM
      VXP2=A11-ALX2
      DENOM=DSQRT(A12*A12+VXP2*VXP2)
      VX2=-A12/DENOM
      VXP2=VXP2/DENOM
      GOTO 2
    5 NSTABX = 2
      GOTO 2
    3 NSTABY=1
      IF(COSY.EQ.1) GOTO 6
      ALY1 = COSY + DSQRT(COSY*COSY-1)
      ALY2 = COSY - DSQRT(COSY*COSY-1)
      A34=TEMP(3,4)
      A33=TEMP(3,3)
      VYP1=A33-ALY1
      DENOM=DSQRT(A34*A34+VYP1*VYP1)
      VY1=-A34/DENOM
      VYP1=VYP1/DENOM
      VYP2=A33-ALY2
      DENOM=DSQRT(A34*A34+VYP2*VYP2)
      VY2=-A34/DENOM
      VYP2=VYP2/DENOM
      GOTO 4
    6 NSTABY = 2
      GOTO 4
      END