SUBROUTINE trmatrix(mdo)
C   *****************************
      IMPLICIT REAL*8(A-H,O-Z), INTEGER (I-N)
      COMMON AMAT(250,6,27),NORLST(2500),
     1XPEL(250),YPEL(250),
     1EXPEL,N,LOGPAR(350)
      COMMON/INOUT/IIN,IOUT,ISOUT,ISO
      COMMON /INPUTT/ KODE(250),NAME(250,4),IADR(250),ELDAT(5000)
     +,MADR(250),KCOUNT,NA
	common /mat1/temp1(6,27),iflag1,nelstrt,nelend
      COMMON/MAT/TEMP(6,27),NORDER,MPRINT,IMAT,NMAT,IFITE,NELM,NOP,
     <NLIST(40)
	common /gplot/igplot,idev,ndispl
      COMMON/PRODCT/KODEPR,NEL,NOF
      COMMON/CONST/PI,TWOPI,CRDEG,CMAGEN,CLIGHT,EMASS,ERAD,ECHG
      COMMON/FOUT/OUTFL(300)
      COMMON/FUNC/BETAX,ALPHAX,ETAX,ETAPX,BETAY,ALPHAY,ETAY,ETAPY,
     <  STARTE,ENDE,DELTAE,DNU,MFPRNT,KOD,NLUM,NINT,NBUNCH
      COMMON/TWF/BETAOX,ALPHOX,ETAOX,ETAPOX,ANUX,
     >            BETAOY,ALPHOY,ETAOY,ETAPOY,ANUY,IE
      COMMON/LENGTH/TLENG,ALENG(250),ACLENG(2500)
      COMMON/LUM/ UO,TAUX,ALPHX,ALPHY,TAUY,
     <  ALPHE,TAUE,SIGE,SIGX,EPSX,
     <SIGY,EPSY,SIGXT,SIGYT,ENERGY,TAUREV
      COMMON/CBEAM/BSIG(6,6),BSIGF(6,6),MBPRT,nform
	common /bmax/bvarc,bharc,bvir,bhir,ir1,ir2,ir3,ir4
 	common /twiher/npnt,jfite(10),ipnt
	common /orbsave/xsave(569),ysave(569),partsave(569,4)
	data ipicb/2/,ipict/2/
	iflag1=0
	if(mdo.eq.0)then
	nint=1
C   SET INTEGRALS INITIALLY TO ZERO
C
      AI1X=0.0D0
      AI2X=0.0D0
      AI3X=0.0D0
      AI4X=0.0D0
      AI5X=0.0D0
      AI1Y=0.0D0
      AI2Y=0.0D0
      AI3Y=0.0D0
      AI4Y=0.0D0
      AI5Y=0.0D0
C
c  set up for partition number calculation with pretzel
	sumk2x=0.0
	sumk2x2=0.0
	sumk3x3=0.0
	hsumk3x3=0.0
	goto 840
	endif

c
	if(mdo.eq.1)then	!displaced orbit in quad
c
c
c    integrate quad times displaced orbit contribution to curlyD
c
c
	sumk2x=2*etaox*st*st*xsave(ie)*al+sumk2x
	sumk2x2=st*st*xsave(ie)**2*al+sumk2x2
	sumk3x3=st**3*xsave(ie)**3*al+sumk3x3
	hsumk3x3=dabs(hav*st**3*xsave(ie)**3*al)+sumk3x3
c	type *,' ie,etaox,st,al,xsave ',ie,etaox,st,al,xsave(ie)
	goto 840
	endif
c
c
C   IS THE NEXT ELEMENT A BEND
C
	if(mdo.ne.3)goto 840	!not a bend
C
      NEXT=IE+1
      NEXTEL=NORLST(NEXT)
      KODN=KODE(NEXTEL)
C
C   GET THE PARAMETERS
C
      IAT=IADR(NEXTEL)
      ANG=ELDAT(IAT+1)*CRDEG
      AL=ELDAT(IAT)
      RAD=AL/ANG
      FIELDN=ELDAT(IAT+2)
      PHI1=ELDAT(IAT+4)*CRDEG
      PHI2=ELDAT(IAT+8)*CRDEG
      RAD2=RAD*RAD
      RAD3=RAD2*RAD
      AK2=(1.0D0-FIELDN)/RAD2
      AK=DSQRT(DABS(AK2))
      ARG=AK*AL
      AL2=AL*AL
      ARG2=AK2*AL*AL
      ARG3=ARG2*ARG
      ARG4=ARG2*ARG2
      ARG5=ARG3*ARG2
      IF(AK2)1,2,3
1     C=DCOSH(ARG)
      S=DSINH(ARG)
      GO TO 4
2     C=1.0D0
      S=0.0D0
      F1=1.0D0
      F2=.5D0
      F3=1.0D0/6.0D0
      F4=1.0D0/3.0D0
      F5=1.0D0/60.0D0
      GO TO 5
3     C=DCOS(ARG)
      S=DSIN(ARG)
4     F1=S/ARG
      F2=(1.0D0-C)/ARG2
      F3=(ARG-S)/ARG3
      F4=(ARG-S*C)/ARG3
      F5=(3.0D0*ARG-4.0D0*S+S*C)/ARG5
C
C   FIND THE VALUES THAT CONTRIBUTE TO THE INTEGRALS
C
    5 IF(KODN.EQ.13)GOTO 870    !vertical bend
      ETAP1=ETAPOX+ETAOX*DTAN(PHI1)/RAD
      AVETA=ETAOX*F1+ETAP1*AL*F2+AL2*F3/RAD
      ETA2=ETAOX*C+ETAP1*S*AL/ARG+(AL2*(1.0D0-C))/(RAD*ARG2)
      AVEFOR=FIELDN*AVETA/RAD3+(ETAOX*DTAN(PHI1)+ETA2*DTAN(PHI2))/
     <  (2.0D0*AL*RAD2)
      ALPHA1=ALPHOX-BETAOX*DTAN(PHI1)/RAD
      GAMMA=(1.0D0+ALPHOX*ALPHOX)/BETAOX
      T1=GAMMA*ETAOX*ETAOX
      T2=2.0D0*ALPHA1*ETAOX*ETAP1
      T3=BETAOX*ETAP1*ETAP1
      T4=-(GAMMA*ETAOX+ALPHA1*ETAP1)*F3*AL
      T5=(ALPHA1*ETAOX+BETAOX*ETAP1)*F2
      T45=2.0D0*AL*(T4+T5)/RAD
      T6=GAMMA*F5*AL2/2.0D0
      T7=ALPHOX*F2*F2*AL
      T8=BETAOX*F4/2.0D0
      T678=AL2*(T6-T7+T8)/RAD2
      HAV=T1+T2+T3+T45+T678
C
C   COMPUTE THE INTEGRALS
C
      AI1X=AI1X+AL*AVETA/RAD
      AI2X=AI2X+AL/RAD2
      AI3X=AI3X+AL/(DABS(RAD))**3
      AI4X=AI4X+(AL*AVETA/RAD3-2.0D0*AL*AVEFOR)
      AI5X=AI5X+AL*HAV/(DABS(RAD))**3
      GO TO 840
  870 ETAP1=ETAPOY+ETAOY*DTAN(PHI1)/RAD
      AVETA=ETAOY*F1+ETAP1*AL*F2+AL2*F3/RAD
      ETA2=ETAOY*C+ETAP1*S*AL/ARG+(AL2*(1.0D0-C))/(RAD*ARG2)
      AVEFOR=FIELDN*AVETA/RAD3+(ETAOY*DTAN(PHI1)+ETA2*DTAN(PHI2))/
     <  (2.0D0*AL*RAD2)
      ALPHA1=ALPHOY-BETAOY*DTAN(PHI1)/RAD
      GAMMA=(1.0D0+ALPHOY*ALPHOY)/BETAOY
      T1=GAMMA*ETAOY*ETAOY
      T2=2.0D0*ALPHA1*ETAOY*ETAP1
      T3=BETAOY*ETAP1*ETAP1
      T4=-(GAMMA*ETAOY+ALPHA1*ETAP1)*F3*AL
      T5=(ALPHA1*ETAOY+BETAOY*ETAP1)*F2
      T45=2.0D0*AL*(T4+T5)/RAD
      T6=GAMMA*F5*AL2/2.0D0
      T7=ALPHOY*F2*F2*AL
      T8=BETAOY*F4/2.0D0
      T678=AL2*(T6-T7+T8)/RAD2
      HAV=T1+T2+T3+T45+T678
C
C   COMPUTE THE INTEGRALS
C
      AI1Y=AI1Y+AL*AVETA/RAD
      AI2Y=AI2Y+AL/RAD2
      AI3Y=AI3Y+AL/(DABS(RAD))**3
      AI4Y=AI4Y+(AL*AVETA/RAD3-2.0D0*AL*AVEFOR)
      AI5Y=AI5Y+AL*HAV/(DABS(RAD))**3

C
C
 840  continue
	if(mdo.ne.4)return
      AI1=AI1X+AI1Y
      AI2=AI2X+AI2Y
      AI3=AI3X+AI3Y
      AI4=AI4X+AI4Y
       AI5=AI5X+AI5Y
C
C
      WRITE(IOUT,1100)AI1X,AI2X,AI3X,AI4X,AI5X
 1100 FORMAT(/,'   SYNCHROTRON INTEGRALS ARE  ',//,
     >3X,'INTEGRAL X1=',E18.11,13X,'INTEGRAL X2=',E18.11,/,
     <       3X,'INTEGRAL X3=',E18.11,13X,'INTEGRAL X4=',E18.11,/,
     <       3X,'INTEGRAL X5=',E18.11,/)
      WRITE(IOUT,1101)AI1Y,AI2Y,AI3Y,AI4Y,AI5Y
1101  FORMAT(/,3X,'INTEGRAL Y1=',E18.11,13X,'INTEGRAL Y2=',E18.11,/,
     <       3X,'INTEGRAL Y3=',E18.11,13X,'INTEGRAL Y4=',E18.11,/,
     <       3X,'INTEGRAL Y5=',E18.11,/)
5000  COMPAC = AI1/TLENG
      IF(DABS(AI2).LT.1.0E-30)RETURN
      AJX = 1.0D0-(AI4X/AI2)
c  sumk2x2 is the integral over effective G2 of displaced orbit in quads
c  sumk2x2 is the integral of quad k and effective g in quads
	CURLYDX=AI4X/AI2
	DELDEN=SUMK2X2/AI2
	DELNUM=SUMK2X/AI4X
	DXPLUS=(AI4X+SUMK2X)/(AI2+SUMK2X2)
	DXMINUS=(AI4X-SUMK2X)/(AI2+SUMK2X2)
	AJXPLUS=1.D0-DXPLUS
	AJXMINUS=1.D0-DXMINUS
      AJY = 1.0D 0-(AI4Y/AI2)
	DY=AI4Y/(AI2+SUMK2X2)
	AJYPRETZ=1.D0-DY
      AJE = 2.0D0+(AI4/AI2)
	DZPLUS=(AI4+SUMK2X)/(AI2+SUMK2X2)
	DZMINUS=(AI4-SUMK2X)/(AI2+SUMK2X2)
	AJEPLUS=2.D0+DZPLUS
	AJEMINUS=2.D0+DZMINUS
      TAUREV = NINT*TLENG/CLIGHT
      FREV = 1.0D0/TAUREV
      ENERGY = STARTE
C
C   OUTPUT THE RESULTS
C
      WRITE(IOUT,734)COMPAC,AJX,AJY,AJE,TAUREV,FREV
734   FORMAT(/,' MOMENTUM COMPACTION=',E15.7,3X,'JX=',F8.5,3X,
     < 'JY=',F8.5,3X,
     <  'JE=',F8.5,/,' TAUREV=',E15.7,'(SEC)',3X,'FREV=',E15.7,
     <  '(HZ)'/)
	WRITE(IOUT,1737)AJXPLUS,AJXMINUS,AJYPRETZ,AJEPLUS,AJEMINUS
1737	FORMAT(/,' JXPLUS=',F8.5,5X,'JXMINUX=',F8.5,
     &	  /,' JYPRETZ=',F8.5,3X,'JEPLUS='F8.5,3X,'JEMINUS=',F8.5)
	WRITE(IOUT,1734)DELDEN,DELNUM,CURLYDX,DXPLUS,DXMINUS
1734	FORMAT(' DELTA DEN=',E15.7,3X,'DELTA NUM=',E15.7,/,
     &	  ' CURLYDX=',E15.7,3X,' DXPLUS=',E15.7,3X,' DXMINUS=',E15.7)
733   ALPH = 7.039346E-06*ENERGY**3*FREV*NINT
      UO = 1.4078692E-02*ENERGY**4*AI2*NINT
      ALPHX = ALPH*(AI2-AI4X)
	ALPHXPLUS=ALPH*AI2*(1.0D0-DXPLUS)
	ALPHXMINUS=ALPH*AI2*(1.0D0-DXMINUS)
      TAUX = 1.0D0/ALPHX
	TAUXPLUS=1.0D0/ALPHXPLUS
	TAUXMINUS=1.0D0/ALPHXMINUS
      ALPHY = ALPH*(AI2-AI4Y)
	ALPHYPRETZ=ALPH*(AI2+SUMK2X2)*(1.D0-DY)
      TAUY = 1.0D0/ALPHY
	TAUYPRETZ=1.D0/ALPHYPRETZ
      ALPHE = ALPH*(2.0D0*AI2+AI4)
	ALPHEPLUS=ALPH*(AI2+SUMK2X2)*(2.D0+DZPLUS)
	ALPHEMINUS=ALPH*(AI2+SUMK2X2)*(2.D0+DZMINUS)
      TAUE = 1.0D0/ALPHE
	TAUEPLUS=1.D0/ALPHEPLUS
	TAUEMINUS=1.D0/ALPHEMINUS
      SIG = 1.211335E-03*ENERGY
      SIGE = SIG*DSQRT(AI3/(2.0D0*AI2+AI4))
	SIGEPLUS=SIG*DSQRT((AI3+SUMK3X3)/((AI2+SUMK2X2)*(2.D0+DZPLUS)))
	SIGEMINUS=SIG*DSQRT((AI3+SUMK3X3)/((AI2+SUMK2X2)*(2.D0+DZMINUS)))
      EPSX = DABS(AI5X/(AI2-AI4X))*SIG**2
	EPSXPLUS=DABS((AI5X+HSUMK3X3)/((AI2+SUMK2X2)*(1.D0-DXPLUS)))*SIG**2
	EPSXMINUS=DABS((AI5X+HSUMK3X3)/((AI2+SUMK2X2)*(1.D0-DXMINUS)))*SIG**2
      SIGX = DSQRT(EPSX*BETAX)
      EPSY = DABS(AI5Y/(AI2-AI4Y))*SIG**2
      SIGY = DSQRT(EPSY*BETAY)
      SIGXT=DSQRT(SIGX**2+(ETAX*SIGE)**2)
      SIGYT=DSQRT(SIGY**2+(ETAY*SIGE)**2)
      WRITE(IOUT,735)ENERGY,UO,ALPHX,ALPHY,ALPHE,
     <  TAUX,TAUY,TAUE,SIGE,EPSX,SIGX,
     <SIGXT,EPSY,SIGY,SIGYT
	WRITE(IOUT,1735)ALPHXPLUS,ALPHXMINUS,TAUXPLUS,TAUXMINUS
1735	FORMAT(' ALPHXPLUS=',E15.7,3X,' ALPHXMINUS=',E15.7,/,
     &	  ' TAUXPLUS=',E15.7,3X,'TAUXMINUS=',E15.7)
	WRITE(IOUT,1736)ALPHYPRETZ,TAUYPRETZ,ALPHEPLUS,ALPHEMINUS,
     &	  TAUEPLUS,TAUEMINUS,SIGEPLUS,SIGEMINUS,EPSXPLUS,EPSXMINUS
1736	FORMAT(/,10X,' AYPRETZ = ',E13.5,'(/SEC)',5X,'TAUYPRETZ=',E13.5,
     &	  '(SEC)'
     &	  ,/,' AEPLUS = ',E13.5,'(/SEC)',5X,' AEMINUS = ',E13.5,'(/SEC)'
     &	  ,/,' TAUEPLUS=',E13.5,'(SEC)',5X,'TAUEMINUS=',E13.5,'(SEC)',
     &	  /,' SIGEPLUS=',E13.5,5X,'SIGEMINUS=',E13.5,
     &	  /,' EPSXPLUS=',E13.5,'(M-RAD)',5X,
     &	 'EPSXMINUS=',E13.5,'(M-RAD)')
      IF(ISO.NE.0)WRITE(ISOUT,735)ENERGY,UO,ALPHX
     <,ALPHY,ALPHE,
     <  TAUX,TAUY,TAUE,SIGE,EPSX,SIGX,
     <SIGXT,EPSY,SIGY,SIGYT
735   FORMAT(//,20X,' ENERGY = ',F8.4,'(GEV)',//,' ENERGY LOSS =',
     < E13.5,'(MEV/REV)'
     < ,/,' AX=',E13.5,'(/SEC)',5X,'AY=',E13.5,'(/SEC)',5X,'AE=',
     <  E13.5,'(/SEC)',/,' TAUX=',E13.5,'(SEC)',4X,'TAUY=',E13.5,'(SEC)'
     <  4X,'TAUE=',E13.5,'(SEC)',/,' SIGE=',E13.5,/,
     <  ' EPSX=',E13.5,'(M-RAD)',2X,'SIGX=',E13.5,'(M.)',5X,
     <'SIGXT=',E13.5,'(M.)',/,
     <  ' EPSY=',E13.5,'(M-RAD)',2X,'SIGY=',E13.5,'(M.)',5X,
     <'SIGYT=',E13.5,'(M.)',/)
c      IF(NLUM.EQ.2) CALL LUMIN
      ENERGY = ENERGY+DELTAE
      IF(ENERGY.LE.ENDE) GO TO 733
      RETURN
      END