C     ***********************
      SUBROUTINE SEED(IEND)
C     ********************************
      IMPLICIT REAL*8(A-H,O-Z), INTEGER (I-N)
	include 'input.inc'
	include 'miss.inc'
       DIMENSION DATA(1)
       INTEGER CLOCK1
c dlr
	data nran/1111111/
c dlr
       NCHAR=0
       NDIM=1
       NDATA=1
       NPRINT=1
       CALL INPUT(ICHAR,NCHAR,DATA,NDIM,IEND,NDATA,NPRINT)
       ISEED = DATA(1)
       IF(ISEED.EQ.0)ISEED=ran(nran)
       ISEED=(2*(ISEED/2))+1
       WRITE(IOUT,10000)ISEED
10000 FORMAT(//,'   THE SEED SELECTED FOR THIS RUN IS :',I16,//)
      RETURN
      END
C     ***********************
      SUBROUTINE SETCOR(IEND)
C     ***********************
      IMPLICIT REAL*8 (A-H,O-Z)
	include 'input.inc'
      COMMON/CORR/CORVAL(600,4),ICRID(600),ICRPOS(600),ICRSET(600),
     >ICROPT(600),NCORR,NCURCR,ICRFLG,ICRCHK,ALMNEL,NPARC
	include 'elements.inc'
      DIMENSION ICHAR(4),OPLIST(6)
      DATA NINE/'9'/
    1 NOP = 0
      NCHAR=4
      INPRT=1
      NDIM=0
      CALL INPUT(ICHAR,NCHAR,OPLIST,NDIM,IEND,NOP,INPRT)
      IF((ICHAR(1).EQ.NINE).AND.(ICHAR(2).EQ.NINE))GO TO 99
      NOP = 6
      NCHAR=0
      INPRT=1
      NDIM=6
      CALL INPUT(ICHAR,NCHAR,OPLIST,NDIM,IEND,NOP,INPRT)
      CALL ELID(ICHAR,NELID)
      ICPOS=OPLIST(1)
      DO 2 ICORR=1,NCORR
      IF(ICPOS.LT.ICRPOS(ICORR))GOTO 3
      IF((NELID.EQ.ICRID(ICORR)).AND.(ICPOS.EQ.ICRPOS(ICORR)))GOTO 4
    2 CONTINUE
    3 WRITE(IOUT,10001)
10001 FORMAT(/,'  NO MATCH WAS FOUND FOR CORRECTOR ID AND POSITION',
     >' IN THE CORRECTOR LIST',/,' DEFINED IN THE CORRECTOR DEFINITION',
     >' OPERATION . RUN IS STOPPED',/)
      CALL HALT
      STOP
    4 ICRSET(ICORR)=1
      ICROPT(ICORR)=OPLIST(2)
      CORVAL(ICORR,1)=OPLIST(3)
      CORVAL(ICORR,2)=OPLIST(4)
      CORVAL(ICORR,3)=OPLIST(5)
      CORVAL(ICORR,4)=OPLIST(6)
      GOTO 1
   99 ICRFLG=1
      WRITE(IOUT,88888)(ICRPOS(IW),ICRID(IW),ICRSET(IW),ICROPT(IW),
     >(CORVAL(IW,JW),JW=1,4),IW=1,NCORR)
88888 FORMAT(' ',4I6,4E12.4)
      RETURN
      END
C     ***********************
      SUBROUTINE SETERR(IEND)
C     ***********************
      IMPLICIT REAL*8(A-H,O-Z), INTEGER (I-N)
	include 'input.inc'
      COMMON/ERR/ERRVAL(7,500),NERELE(500),NERPAR(7,500),NERR,NEROPT,
     > NERRE,MERSEL(500),NERNGE(500),MERNGE(2,10,500),MERFLG
      DIMENSION DATA(20),ICHAR(4)
C INITIALIZE TO 0 ALL ARRAYS DEFINED IN THIS ROUTINE
      DO 5 INER=1,500
      MERSEL(INER)=0
      NERNGE(INER)=0
      DO 5 JNER=1,10
      DO 5 KNER=1,2
    5 MERNGE(KNER,JNER,INER)=0
      NDIM=20
      NPRINT=1
      NCHAR=0
      NDATA=2
      CALL INPUT(ICHAR,NCHAR,DATA,NDIM,IEND,NDATA,NPRINT)
      NEROPT=DATA(1)
      NERRE=DATA(2)
      DO 1 IER=1,NERRE
      NCHAR=4
      NDATA=0
      NDIM=0
      CALL INPUT(ICHAR,NCHAR,DATA,NDIM,IEND,NDATA,NPRINT)
      CALL ELID(ICHAR,NELID)
      DO 3 JM=1,NERR
      IF(NELID.EQ.NERELE(JM))GOTO 4
    3 CONTINUE
      WRITE(IOUT,99994)
99994 FORMAT(/,'  NAME IS NOT FOUND IN THE LIST OF THE',
     >' ELEMENTS WITH ERROR AS DEFINED IN ERRDAT',/)
      STOP
    4 MERSEL(IER)=JM
      NDATA = 1
      NCHAR=0
      NDIM=20
      CALL INPUT(ICHAR,NCHAR,DATA,NDIM,IEND,NDATA,NPRINT)
      NERNGE(IER)=DATA(1)
      IF((DATA(1).EQ.0.0D0).OR.(DATA(1).EQ.-1.0D0))GOTO 1
      NDATA=NERNGE(IER)*2
      CALL INPUT(ICHAR,NCHAR,DATA,NDIM,IEND,NDATA,NPRINT)
      NFRNGE=NERNGE(IER)
      DO 2 IRNGE=1,NFRNGE
      MERNGE(1,IRNGE,IER)=DATA((2*IRNGE)-1)
    2 MERNGE(2,IRNGE,IER)=DATA(2*IRNGE)
    1 CONTINUE
      MERFLG=1
      IOPT=NEROPT+1
C     WRITE(IOUT,99000)((NERELE(IN),MERSEL(IN),NERNGE(IN),
C    >MERNGE(1,1,IN),MERNGE(2,1,IN)),IN=1,NERR)
99000 FORMAT(/,' IN ESET',5I5)
      GOTO(10,11,12,13),IOPT
   10 WRITE(IOUT,99990)
99990 FORMAT(/,' NO RANDOM GENERATOR IS USED IN SETTING UP THE'
     >,' ERRORS',/)
      RETURN
   11 WRITE(IOUT,99991)
99991 FORMAT(/,' THE ERROR VALUES ARE MULTIPLIED RANDOMLY BY',
     >/,' +1 AND -1',/)
      RETURN
   12 WRITE(IOUT,99992)
99992 FORMAT(/,' THE ERROR VALUES ARE MULTIPLIED BY THE VALUES',
     >/,' OF A UNIFORM RANDOM DISTRIBUTION WHOSE SIGMA IS 1',/)
      RETURN
   13 WRITE(IOUT,99993)
99993 FORMAT(/,' THE ERROR VALUES ARE MULTIPLIED BY THE VALUES',
     >/,' OF A GAUSSIAN DISTRIBUTION WHOSE SIGMA IS 1 AND THAT IS',/,
     >' TRUNCATED AT 2 SIGMAS',/)
      RETURN
      END
C     ***********************
      SUBROUTINE SETMIS(IEND)
C     ********************************
      IMPLICIT REAL*8(A-H,O-Z), INTEGER (I-N)
	include 'input.inc'
	include 'miss.inc'
      DIMENSION DATA(20),ICHAR(4)
      NDIM=20
      NPRINT=1
      NCHAR=0
      NDATA=2
      CALL INPUT(ICHAR,NCHAR,DATA,NDIM,IEND,NDATA,NPRINT)
      NOPT=DATA(1)
      NMISE=DATA(2)
      DO 1 IMS=1,NMISE
      NCHAR=4
      NDATA=0
      NDIM=0
      CALL INPUT(ICHAR,NCHAR,DATA,NDIM,IEND,NDATA,NPRINT)
      CALL ELID(ICHAR,NELID)
      DO 3 JM=1,NMIS
      IF(NELID.EQ.MISELE(JM))GOTO 4
    3 CONTINUE
      WRITE(IOUT,99994)
99994 FORMAT(/,'  NAME IS NOT FOUND IN THE LIST OF THE',
     >' MISALIGNED ELEMENTS DEFINED IN MISDAT',/)
      STOP
    4 MISSEL(IMS)=JM
      NDATA = 1
      NCHAR=0
      NDIM=20
      CALL INPUT(ICHAR,NCHAR,DATA,NDIM,IEND,NDATA,NPRINT)
      NMRNGE(IMS)=DATA(1)
      IF((DATA(1).EQ.0.0D0).OR.(DATA(1).EQ.-1.0D0))GOTO 1
      NDATA=NMRNGE(IMS)*2
      CALL INPUT(ICHAR,NCHAR,DATA,NDIM,IEND,NDATA,NPRINT)
      NFRNGE=NMRNGE(IMS)
      DO 2 IRNGE=1,NFRNGE
      MSRNGE(1,IRNGE,IMS)=DATA((2*IRNGE)-1)
    2 MSRNGE(2,IRNGE,IMS)=DATA(2*IRNGE)
    1 CONTINUE
      MISFLG=1
      IOPT=NOPT+1
      GOTO(10,11,12,13),IOPT
   10 WRITE(IOUT,99990)
99990 FORMAT(/,' NO RANDOM GENERATOR IS USED IN THE MISALIGNEMENT',/)
      RETURN
   11 WRITE(IOUT,99991)
99991 FORMAT(/,' THE MISALIGNEMENT VALUES ARE MULTIPLIED RANDOMLY BY',
     >/,' +1 AND -1',/)
      RETURN
   12 WRITE(IOUT,99992)
99992 FORMAT(/,' THE MISALIGNMENT VALUES ARE MULTIPLIED BY THE VALUES',
     >/,' OF A UNIFORM RANDOM DISTRIBUTION WHOSE SIGMA IS 1',/)
      RETURN
   13 WRITE(IOUT,99993)
99993 FORMAT(/,' THE MISALIGNMENT VALUES ARE MULTIPLIED BY THE VALUES',
     >/,' OF A GAUSSIAN DISTRIBUTION WHOSE SIGMA IS 1 AND THAT IS',/,
     >' TRUNCATED AT 2 SIGMAS',/)
      RETURN
      END
C     ***********************
      SUBROUTINE SETSYN(IEND)
C     ********************************
      IMPLICIT REAL*8(A-H,O-Z), INTEGER (I-N)
	include 'input.inc'
      COMMON/SYNC/ENOM,SYNDEL,ISYNFL,ISYNGO
      DIMENSION DATA(2),ICHAR(4)
      NDIM=2
      NDATA=2
      NCHAR=0
      NPRINT=1
      CALL INPUT(ICHAR,NCHAR,DATA,NDIM,IEND,NDATA,NPRINT)
      ENOM=DATA(1)
      ISYNFL=DATA(2)
      RETURN
      END
C     *************************
c      SUBROUTINE SETOUT
cC     ******************************
c      IMPLICIT REAL*8(A-H,O-Z), INTEGER (I-N)
c	include 'input.inc'
c      COMMON/FOUT/OUTFL(300)
c      COMMON/CBEAM/BSIG(6,6),BSIGF(6,6),MBPRT,nform
c      COMMON/MAT/TEMP(6,27),NORDER,MPRINT,IMAT,NMAT,IFITE,NELM,NOP,
c     <NLIST(40)
c      include 'analc.inc'
c      COMMON/TWF/BETAOX,ALPHOX,ETAOX,ETAPOX,ANUX,
c     >            BETAOY,ALPHOY,ETAOY,ETAPOY,ANUY,IE
c      COMMON/FITD/AVEBX,AVEAX,AVEBY,AVEAY,AVENUX,AVENUY,
c     >AVER11,AVER12,AVER21,AVER22,AVER33,AVER34,AVER43,AVER44
c      COMMON/LAYOUT/XHI,YHI,ZHI,THETAI,PHIHI,PSIHI,CONVH,
c     >SRH,XRH,YRH,ZRH,THETAR,PHIRH,PSIRH,NLAY
c      DIMENSION RELAY2(10),RELAY1(20),RELAY3(162)
c      DIMENSION RELAY4(6),RELAY5(8),RELAY6(8)
c      EQUIVALENCE(RELAY4(1),AVEBX)
c      EQUIVALENCE(RELAY5(1),AVER11)
c      EQUIVALENCE(RELAY6(1),SRH)
c      EQUIVALENCE(COMPF,RELAY1(1))
c      EQUIVALENCE(RELAY2(1),BETAOX)
c      EQUIVALENCE(TEMP(1,1),RELAY3(1))
c      IF((IE.EQ.IFITE).AND.(IE.NE.0))GOTO 4
c      DO 1 I=1,10
c      OUTFL(I+20)=RELAY2(I)
c    1 CONTINUE
c      DO 2 I=1,20
c      OUTFL(I)=RELAY1(I)
c    2 CONTINUE
c      DO 3 I = 1,162
c      OUTFL(I+100)=RELAY3(I)
c    3 CONTINUE
c      IBS=1
c      DO 6 IB=1,6
c      DO 6 JB=IB,6
c      OUTFL(40+IBS)=BSIGF(IB,JB)
c    6 IBS=IBS+1
c      DO 8 IR=1,6
c      OUTFL(92+IR)=RELAY4(IR)
c    8 CONTINUE
c      DO 9 IR=1,8
c      OUTFL(270+IR)=RELAY5(IR)
c    9 CONTINUE
c      DO 10 IR=1,8
c      OUTFL(280+IR)=RELAY6(IR)
c   10 CONTINUE
c      RETURN
c    4 DO 5 I=1,10
c    5 OUTFL(I+30)=RELAY2(I)
c      IBS=1
c      DO 7 IB=1,6
c      DO 7 JB=IB,6
c      OUTFL(70+IBS)=BSIGF(IB,JB)
c    7 IBS=IBS+1
c      RETURN
c      END
C     *************************
      SUBROUTINE SETUP
C     *************************
      IMPLICIT REAL*8(A-H,O-Z)
	include 'input.inc'
	character*59 name
      ISOUT=6
      IIN=31
	call strget(' Input file ',name)
	open(unit=iin,file=name,readonly,type='OLD')
c
	call strip_dirname(name,istart,iend,leftbrak,ightbrak)
c
	filename=name(istart:iend)
      IOUT=34
	type 1,name(istart:iend)//'_out'
1	format(1x,' Open file ',a,' to write general output.')
	open(unit=iout,file=name(istart:iend)//'_out',type='new')
      ISO=0
      RETURN
      END
C     ***********************
      SUBROUTINE SOLQUA
C     ***********************
      IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
	include 'input.inc'
      COMMON/CONST/PI,TWOPI,CRDEG,CMAGEN,CLIGHT,EMASS,ERAD,ECHG
      COMMON/PRODCT/KODEPR,NEL,NOF                                         
	include 'tempstor.inc'
	include 'elements.inc'
      IAD=IADR(N)
      AL =ELDAT(IAD)
      AKQ=ELDAT(IAD+1)
      AKS=ELDAT(IAD+2)
      DO 100 I=1,6
      DO 100 J=1,NOF
      AMAT(N,I,J)=0.0D0
      IF(I.EQ.J) AMAT(N,I,J)=1.0D0
  100 CONTINUE
      IF(AKQ.NE.0.0D0.OR.AKS.NE.0.0D0)GO TO 2
      AMAT(N,1,2)=AL
      AMAT(N,3,4)=AL
      IF(NORDER.EQ.1)RETURN
      AMAT(N,5,13)=AL/2.0D0
      AMAT(N,5,22)=AL/2.0D0
      RETURN
    2 IKQ=1
      IF(AKQ.LT.0.0D0)IKQ=-1
      AKQ4=AKQ*AKQ
      AKS2=AKS*AKS
      AKS3=AKS2*AKS
      AKS4=AKS2*AKS2
      Q2=DSQRT(AKQ4+4.0D0*AKS4)
      Q2I=1.0D0/Q2
      SK1=DSQRT(2.0D0*AKS2+Q2)
      SK3=DSQRT(DABS(Q2-2.0D0*AKS2))
      SS=DSIN(SK1*AL)
      SC=DCOS(SK1*AL)
      BS=DSINH(SK3*AL)
      BC=DCOSH(SK3*AL)
      SKP=0.5D0*(SK1+SK3)
      SKM=0.5D0*(SK1-SK3)
      SKP2=SKP*SKP
      SKP3=SKP2*SKP
      SKM2=SKM*SKM
      SKM3=SKM2*SKM
      AMAT(N,1,1)=Q2I*(SKP2*SC+SKM2*BC)
      AMAT(N,1,2)=Q2I*(SKP*SS-SKM*BS)
      AMAT(N,1,3)=Q2I*AKS*(SKM*SS+SKP*BS)
      AMAT(N,1,4)=Q2I*AKS*(BC-SC)
      AMAT(N,2,1)=Q2I*(-SKP3*SS-SKM3*BS)
      AMAT(N,2,2)=Q2I*(SKP2*SC+SKM2*BC)
      AMAT(N,2,3)=Q2I*AKS3*(SC-BC)
      AMAT(N,2,4)=Q2I*AKS*(SKP*SS-SKM*BS)
      AMAT(N,3,1)=Q2I*AKS*(SKM*BS-SKP*SS)
      AMAT(N,3,2)=Q2I*AKS*(SC-BC)
      AMAT(N,3,3)=Q2I*(SKM2*SC+SKP2*BC)
      AMAT(N,3,4)=Q2I*(SKM*SS+SKP*BS)
      AMAT(N,4,1)=Q2I*AKS3*(BC-SC)
      AMAT(N,4,2)=Q2I*AKS*(-SKM*SS-SKP*BS)
      AMAT(N,4,3)=Q2I*(SKP3*BS-SKM3*SS)
      AMAT(N,4,4)=Q2I*(SKM2*SC+SKP2*BC)
      IF(NORDER.EQ.1)GOTO 300
      DAKS=-AKS
      DAKQ=-AKQ
      DAKQ4=2.0D0*AKQ*DAKQ
      DAKS2=2.0D0*AKS*DAKS
      DAKS3=3.0D0*AKS2*DAKS
      DAKS4=4.0D0*AKS3*DAKS

      DQ2=0.5D0*Q2I*(DAKQ4+4.0D0*DAKS4)
      Q2I2=Q2I*Q2I
      DQ2I=-Q2I2*DQ2
      DSK1=(2.0D0*DAKS2+DQ2)/(2.0D0*SK1)
      DSK3=0.0D0
      IF(SK3.NE.0.0D0)DSK3=(DQ2-2.0D0*DAKS2)/(2.0D0*SK3)
      DSS=SC*DSK1*AL
      DSC=-SS*DSK1*AL
      DBS=BC*DSK3*AL
      DBC=BS*DSK3*AL
      DSKP=0.5D0*(DSK1+DSK3)
      DSKM=0.5D0*(DSK1-DSK3)
      DSKP2=2.0D0*SKP*DSKP
      DSKP3=3.0D0*SKP2*DSKP
      DSKM2=2.0D0*SKM*DSKM
      DSKM3=3.0D0*SKM2*DSKM
      AMAT(N,1,12)=DQ2I*(SKP2*SC+SKM2*BC)+
     > Q2I*(DSKP2*SC+SKP2*DSC+DSKM2*BC+SKM2*DBC)
      AMAT(N,1,17)=DQ2I*(SKP*SS-SKM*BS)+
     > Q2I*(DSKP*SS+SKP*DSS-DSKM*BS-SKM*DBS)
      AMAT(N,1,21)=(DQ2I*AKS+Q2I*DAKS)*(SKM*SS+SKP*BS)
     > +Q2I*AKS*(DSKM*SS+SKM*DSS+DSKP*BS+SKP*DBS)
      AMAT(N,1,24)=(DQ2I*AKS+Q2I*DAKS)*(BC-SC)
     > +Q2I*AKS*(DBC-DSC)
      AMAT(N,2,12)=DQ2I*(-SKP3*SS-SKM3*BS)
     > +Q2I*(-DSKP3*SS-SKP3*DSS-DSKM3*BS-SKM3*DBS)
      AMAT(N,2,17)=AMAT(N,1,12)
      AMAT(N,2,21)=(DQ2I*AKS3+Q2I*DAKS3)*(SC-BC)
     > +Q2I*AKS3*(DSC-DBC)
      AMAT(N,2,24)=(DQ2I*AKS+Q2I*DAKS)*(SKP*SS-SKM*BS)
     > +Q2I*AKS*(DSKP*SS+SKP*DSS-DSKM*BS-SKM*DBS)
      AMAT(N,3,12)=-AMAT(N,2,24)
      AMAT(N,3,17)=-AMAT(N,1,24)
      AMAT(N,3,21)=DQ2I*(SKM2*SC+SKP2*BC)
     > +Q2I*(DSKM2*SC+SKM2*DSC+DSKP2*BC+SKP2*DBC)
      AMAT(N,3,24)=DQ2I*(SKM*SS+SKP*BS)
     > +Q2I*(DSKM*SS+SKM*DSS+DSKP*BS+SKP*DBS)
      AMAT(N,4,12)=-AMAT(N,2,21)
      AMAT(N,4,17)=-AMAT(N,1,21)
      AMAT(N,4,21)=DQ2I*(SKP3*BS-SKM3*SS)
     > +Q2I*(DSKP3*BS+SKP3*DBS-DSKM3*SS-SKM3*DSS)
      AMAT(N,4,24)=AMAT(N,3,21)
      AISSSC=0.0D0
      IF(SK1.NE.0.0D0)AISSSC=SS*SS/(2.0D0*SK1)
      AIBSBC=0.0D0
      IF(SK3.NE.0.0D0)AIBSBC=BS*BS/(2.0D0*SK3)
      AISSBC=Q2I*(SK3*SS*BS-SK1*SC*BC+SK1)
      AISCBS=Q2I*(SK3*SC*BC+SK1*SS*BS-SK3)
      AISS2=0.0D0
      AISC2=AL
      AIBS2=0.0D0
      AIBC2=AL
      IF(SK1.NE.0.0D0)AISS2=0.5D0*(AL-SS*SC/SK1)
      IF(SK1.NE.0.0D0)AISC2=0.5D0*(AL+SS*SC/SK1)
      IF(SK3.NE.0.0D0)AIBS2=0.5D0*(BS*BC/SK3 - AL)
      IF(SK3.NE.0.0D0)AIBC2=0.5D0*(AL+BS*BC/SK3)
      AISSBS=Q2I*(SK3*SS*BC-SK1*SC*BS)
      AISCBC=Q2I*(SK1*SS*BC+SK3*SC*BS)
      AKS5=AKS4*AKS
      AKS6=AKS5*AKS
      SKP4=SKP3*SKP
      SKP5=SKP4*SKP
      SKP6=SKP5*SKP
      SKM4=SKM3*SKM
      SKM5=SKM4*SKM
      SKM6=SKM5*SKM
      Q4I=Q2I*Q2I
      AMAT(N,5,7)=Q4I*0.5D0*(SKP6*AISS2+SKM6*AIBS2-2.0D0*SKP3*SKM3
     >   *AISSBS+AKS6*(AISC2+AIBC2-2.0D0*AISCBC))
      AMAT(N,5,8)=Q4I*(-SKP5*AISSSC-SKP3*SKM2*AISSBC-SKM3*SKP2*AISCBS
     > -SKM5*AIBSBC-AKS4*(-SKM*AISSSC-SKP*AISCBS+SKM*AISSBC+SKP*AIBSBC))
      AMAT(N,5,9)=Q4I*AKS3*((SKP3-SKM3)*(AISSBC-AISSSC)
     >  +(SKP3+SKM3)*(AIBSBC-AISCBS))
      AMAT(N,5,10)=Q4I*(AKS*(-SKP4*AISS2-(SKP3*SKM+SKM3*SKP)*AISSBS
     >  +SKM4*AIBS2)+AKS3*(-SKM2*AISC2-(SKP2-SKM2)*AISCBC+SKP2*AIBC2))
      AMAT(N,5,13)=Q4I*0.5D0*(SKP4*AISC2+2.0D0*SKP2*SKM2*AISCBC
     >  +SKM4*AIBC2+AKS2*(SKM2*AISS2+2.0D0*SKP*SKM*AISSBS+SKP2*AIBS2))
      AMAT(N,5,14)=Q4I*(AKS3*(SKP2*AISC2-(SKP2-SKM2)*AISCBC
     >  -SKM2*AIBC2)-AKS*(-SKM4*AISS2+(SKM*SKP3-SKP*SKM3)*AISSBS+
     >  SKP4*AIBS2))
      AMAT(N,5,15)=Q4I*AKS*((SKP3-SKM3)*AISSSC-(SKP2*SKM+SKP*SKM2)
     >  *AISCBS+(SKM2*SKP-SKM*SKP2)*AISSBC-(SKM3+SKP3)*AIBSBC)
      AMAT(N,5,18)=Q4I*0.5D0*(AKS6*(AISC2-2.0D0*AISCBC+AIBC2)
     >  +SKM6*AISS2-2.0D0*SKM3*SKP3*AISSBS+SKP6*AIBS2)
      AMAT(N,5,19)=Q4I*(AKS4*(SKP*AISSSC-SKM*AISCBS-SKP*AISSBC
     >  +SKM*AIBSBC)-SKM5*AISSSC-SKM3*SKP2*AISSBC+SKP3*SKM2*AISCBS
     >  +SKP5*AIBSBC)
      AMAT(N,5,22)=Q4I*0.5D0*(AKS2*(SKP2*AISS2-2.0D0*SKP*SKM*AISSBS
     > +SKM2*AIBS2)+SKM4*AISC2+2.0D0*SKM2*SKP2*AISCBC+SKP4*AIBC2)
  300 IF(IKQ.GE.0)RETURN
C
C     SET UP THE 90 DEGREE KICK MATRIX IN TEMP
C
      DO 10 IX = 1,6
      DO 10 IY = 1,NOF
      TEMP(IX,IY)=0.0D0
10    CONTINUE
      TEMP(1,3)=1
      TEMP(2,4)=1.0D0
      TEMP(3,1)=-1.0D0
      TEMP(4,2)=-1.0D0
      TEMP(5,5)=1.0D0
      TEMP(6,6)=1.0D0
       KODEPR=15
      CALL PROMAT
C
C  SET -90 DEGREE KICK MATRIX IN AMAT
C
      DO 200 I=1,6
      DO 200 J=1,NOF
      AMAT(N,I,J)=0.0D0
 200  CONTINUE
      AMAT(N,1,3)=-1.0D0
      AMAT(N,2,4)=-1.0D0
      AMAT(N,3,1)=1.0D0
      AMAT(N,4,2)=1.0D0
      AMAT(N,5,5)=1.0D0
      AMAT(N,6,6)=1.0D0
      KODEPR=8
      CALL PROMAT
C
C     PUT TEMP INTO AMAT(N,6,27)
C
      DO 20 IX=1,6
      DO 20 IY=1,NOF
  20  AMAT(N,IX,IY) = TEMP(IX,IY)
      RETURN
      END
C     **********************
c      SUBROUTINE SQFCT(F)
cC     *************************
c      IMPLICIT REAL*8(A-H,O-Z), INTEGER (I-N)
c	include 'input.inc'
c	include 'fitl.inc'
c      COMMON/FOUT/OUTFL(300)
c      F=0.0D0
c      IF(ISTART.NE.0)GO TO 1
c      DO 20 I=1,NCOND
c      VAL0=OUTFL(NVAL(I))
c   20 RVAL(I)=VAL0-(VAL0-VALF(I))/NDIV
c	type 22
c22	format(/)
c	type 21,(outfl(nval(i)),i=1,ncond)
c21	format(5e12.4)
c	type 21,(valf(i),i=1,ncond)
c      ISTART=1
c    1 DO 10 I=1,NCOND
c   10 F=F+((OUTFL(NVAL(I))-RVAL(I))*WGHT(I))**2
c      F=F/WV
c      F=DSQRT(F)
c	type *,' FIG ',f
c      RETURN
c      END
C     **************************
      SUBROUTINE SYNPRE(IAD,NEL)
C     **************************
      IMPLICIT REAL*8(A-H,O-Z), INTEGER (I-N)
	include 'input.inc'
      COMMON/SYNC/ENOM,SYNDEL,ISYNFL,ISYNGO
	include 'track.inc'
	include 'elements.inc'
      COMMON/CONST/PI,TWOPI,CRDEG,CMAGEN,CLIGHT,EMASS,ERAD,ECHG
      ANG=CRDEG*ELDAT(IAD+1)
      SYNDEL=1.408D-05*(ENOM**3)*ANG*ANG/ALENG(NEL)
      ISYNGO=1
      RETURN
      END
C     ***********************
      SUBROUTINE TRACKT
C     **********************
        IMPLICIT REAL*8 (A-H,O-Z)
cdlr
	common /opname/noptype(4)
	include 'elements.inc'
	character*3 ans,anss
cdlr
	include 'analc.inc'
	common /xampyamp/xampmax,yampmax
	include 'input.inc'
      COMMON/DETL/DENER(15),NH,NV,NVH,NHVP(105),MDPRT,NDENER,
     1NUXS(45),NUX(45),NUYS(45),NUY(45),NCO,NHNVHV,MULPRT,NSIG
	include 'fitl.inc'
cdlr
      COMMON/FUNC/BETAX,ALPHAX,ETAX,ETAPX,BETAY,ALPHAY,ETAY,ETAPY,
     <  STARTE,ENDE,DELTAE,DNU,MFPRNT,KOD,NLUM,NINT,NBUNCH
cdlr
	common /middle/interm,intprt,iemitt,idec,emit(2),nesave(1000)
	include 'track.inc'
       COMMON/V/VV(27),X1,X2,X3,X4,X5,X6
      COMMON/PLT/
     1XMIN,XMAX,YMIN,YMAX,XPMIN,XPMAX,YPMIN,YPMAX,
     >DELMIN,DELMAX,DNUMIN,DNUMAX,DBMIN,DBMAX,
     2MXXPR,MYYPR,MXY,MALL,NPLOT,NCCUM,NGRAPH,NCOL,NLINE
      COMMON/CONST/PI,TWOPI,CRDEG,CMAGEN,CLIGHT,EMASS,ERAD,ECHG
      COMMON/TRI/WCO(15,6),GEN(5,4),PGEN(75,6),DIST,
     <A1(15),B1(15),A2(15),B2(15),XCOR(3,15),XPCOR(3,15),
     1AL1(15),AL2(15),MESH,NIT,NENER,NITX,NITS,NITM1,NANAL,NOPRT
	include 'geom.inc'
	include 'miss.inc'
	include 'tempstor.inc'
      COMMON/ERR/ERRVAL(7,500),NERELE(500),NERPAR(7,500),NERR,NEROPT,
     > NERRE,MERSEL(500),NERNGE(500),MERNGE(2,10,500),MERFLG
      COMMON/ERSAV/SAV(7),SAVMAT(6,27),IERSET,IER,IDE
      COMMON/MONIT/VALMON(12,4,3),MNAME(600,4),MONPOS(600),NMON
      COMMON/MONFIT/VALFA(12),WGHTA(12),ERRA(12),
     >AMULTA(12,6),ADDA(12,6),DELA(12),NPARA(12,6),NELFA(12,6),
     >NPVARA(12),INDA(12,6),VALR(12),
     >NMONA(12),NVALA(12),NVARA,NCONDA,IALFLG,MONFLG,MONLST,NOPTER,
     >IAFRST,ISDBEG,IMONSD,IMSBEG
      COMMON/CORR/CORVAL(600,4),ICRID(600),ICRPOS(600),ICRSET(600),
     >ICROPT(600),NCORR,NCURCR,ICRFLG,ICRCHK,ALMNEL,NPARC
      COMMON/ORBIT/SIZEX,SIZEY,RMSX,RMSY,RMSIX,RMSIY,
     >RTEMPX,RTEMPY,RMSPX(5),RMSPY(5),RPX,RPY,
     >RMAXX,RMAXY,RMINX,RMINY,MAXX,MAXY,MINX,MINY,PLENG,
     >IRNG,IRANGE(5),NPRORB,IORB,IREF,IPAGE,IPOINT
      COMMON/SYNC/ENOM,SYNDEL,ISYNFL,ISYNGO
	include 'partsave.inc'
      LOGICAL MXXPR(101,51),MYYPR(101,51),MXY(101,51)
	dimension lab(4)
      DIMENSION INDEL(6,18),INDOR2(13,6,8),INDORD(13,6,18)
      DIMENSION INDOR1(13,6,10)
	common/ detmax/detxmax,detymax,scmax,detxymax,kspotmax
      EQUIVALENCE (INDOR1(1,1,1),INDORD(1,1,1))
      EQUIVALENCE(INDOR2(1,1,1),INDORD(1,1,11))
	data lab/'K','8','4','5'/
      DATA INDEL/2,1,2,1,3,1,
     <12,12,8,8,13,1,4,4,4,4,7,1,8,7,6,5,3,1,10,10,8,8,7,1,
     <6*0,6*0,4,2,4,2,1,1, 4,2,
     <4,2,1,1,2,2,2,2,1,1,6*0,3,7,2,6,1,1,
     <6*0,  8,8,12,12,13,1,  6*0, 8,8,8,8,11,1,6*0,6*0/
      DATA INDOR1/1,2,11*0, 2,12*0, 3,4,11*0, 4,12*0, 5,13,22,10*0,
     +  6,12*0, 1,2,6,7,8,12,13,17,18,19,22,27,0,
     <  1,2,6,7,8,12,13,17,18,19,22,27,0,
     <  3,4,9,10,14,15,21,24,5*0, 3,4,9,10,14,15,21,24,5*0,
     <  1,2,5,6,7,8,12,13,17,18,19,22,27, 6,12*0,
     <  1,2,12,17,9*0, 1,2,12,17,9*0, 3,4,21,24,9*0, 3,4,21,24,9*0,
     <  5,7,8,13,18,19,22,6*0, 6,12*0, 1,2,7,8,13,18,19,22,5*0,
     <  2,7,8,13,18,19,22,6*0, 3,4,9,10,14,15,7*0, 4,9,10,14,15,8*0,
     <5,13,22,10*0,6,12*0, 1,2,7,8,12,13,17,18,19,22,3*0,
     <1,2,7,8,12,13,17,18,19,22,3*0,3,4,9,10,14,15,21,24,5*0,
     <3,4,9,10,14,15,21,24,5*0,5,7,8,13,18,19,22,6*0,6,12*0,
     <78*0,78*0,1,2,3,4,9*0,2,4,11*0,1,2,3,4,9*0,
     <2,4,11*0, 5,12*0,6,12*0,
     <  1,2,3,4,9*0, 2,4,11*0,
     <  1,2,3,4,9*0, 2,4,11*0, 5,12*0, 6,12*0, 1,2,11*0, 1,2,11*0,
     <  3,4,11*0, 3,4,11*0, 5,12*0, 6,12*0/
       DATA INDOR2/78*0, 1,7,18,10*0,
     <  1,2,7,8,12,18,19,6*0, 3,9,11*0, 3,4,9,10,14,21,7*0,
     <  5,12*0, 6,12*0,78*0,
     <1,2,9,10,12,14,15,17,5*0,  1,2,9,10,12,14,15,17,5*0,
     <3,4,6,7,8,13,18,19,21,22,24,27,0,
     <3,4,6,7,8,13,18,19,21,22,24,27,0,
     <3,4,5,6,7,8,13,18,19,21,22,24,27, 6,12*0,
     <78*0,1,2,3,4,12,17,21,24,5*0,1,2,3,4,12,17,21,24,5*0,
     >1,2,3,4,12,17,21,24,5*0,1,2,3,4,12,17,21,24,5*0,
     >5,7,8,9,10,13,14,15,18,19,22,2*0,6,12*0,
     >78*0,78*0/
C
C   PRINT INITIAL POSITIONS
C
      NORDER=2
c     	  call bastart
        IF(NPRINT.NE.-2)CALL TRAKPR(-1,1)
cdlr
	if(idlr.eq.0)then
	  type 3
3	format(' Begin tracking at element number ?',$)
	accept *,itrbeg
c	  type 1
1	  format(' Do you want to plot  trajectory ? ',$)
c	  accept 2,ans
	  ans='y'
2	  format(a)
	    idlr=1
	  endif
cdlr
      NXREQ = NPLOT +NCTURN
      LSTREQ = NTURN/NPLOT*NPLOT
C
C   TURN LOOP NCTURN -> NTURN
C
c
	xampmax=0.  !max x amplitude for particle 1
	yampmax=0.  !max y amplitude for particle 1
c
   10 NCTURN=NCTURN+1
      ISYNGO=0
      IF(IALFLG.EQ.0)GOTO 601
      IF(IALFLG.EQ.1)GOTO 602
      ITRBEG=IAFRST
      ITREND=MONLST
      IXS=ISDBEG
      IMONSD=IMSBEG
      GOTO 603
  601 continue
c	ITRBEG=1
      ITREND=NELM
      IXS=ISEED
      IMONSD=ISEED
      GOTO 603
  602 ITRBEG=1
      ITREND=IAFRST-1
      IXS=ISEED
      IMONSD=ISEED
  603 ILIST=1
C     WRITE(IOUT,33331)IALFLG,ISEED,IXS,ISDBEG,ITRBEG,ITREND
33331 FORMAT(' IN TRACKT IALFLG,ISEED,IXS,ISDBEG,ITRBEG,ITREND =',/,
     >I3,3I11,2I5)
      DO 110 IEdum=1,nelm
	ie=iedum+ITRBEG-1
	if(ie.gt.nelm)ie=ie-nelm
       IEP=IE
      NEL=NORLST(IE)
      MNEL=NEL
      IAD=IADR(NEL)
      KD=KODE(NEL)
      IF(((KD.EQ.1).OR.(KD.EQ.13).OR.(KD.EQ.14)).AND.(ISYNFL.NE.0))
     >CALL SYNPRE(IAD,NEL)
      IF(ICRFLG.EQ.1)CALL CORCHK(IE)
C     WRITE(IOUT,77771)ICRFLG,ICRCHK
77771 FORMAT('  IN TRACKT ICRFLG,ICRCHK ARE',2I6)
      IF(MISFLG.EQ.1)CALL MISCHK
cdlr
	if(misflg.eq.1)call missave(ie)
	if(misflg.eq.-1)call misget(ie)
cdlr
C     IF((IREF.EQ.1).AND.(MCHFLG.NE.0))WRITE(IOUT,77772)MCHFLG,IE
77772 FORMAT('  IN TRACKT MCHFLG IS ',2I6)
      IF((MERFLG.EQ.1).OR.(ICRCHK.EQ.2))CALL ESET
C     IF(MCHFLG.NE.0)WRITE(IOUT,77772)MCHFLG,IE
      IF(IALFLG.NE.0)CALL MONCHK(IE)
      N=MADR(NEL)
      NT=KODE(NEL)
      IF(NT.EQ.5)CALL MULTIT(IAD)
      IF(NT.EQ.12)CALL ARBIT(IAD,NEL)
      IF(NT.EQ.11) GO TO 111
C
C   PROCEs THE PARTICLES
C
      DO 65 I=1,NPART
      IF(.NOT.LOGPAR(I)) GO TO 65
      TEST=(PART(I,1)/XPEL(NEL))**2+(PART(I,3)/YPEL(NEL))**2
      TEST=DSQRT(TEST)
      IF(TEST.LT.EXPEL) GO TO 20
      WRITE(IOUT,10024)I,IE,(NAME(NEL,IZ),IZ=1,4),NCTURN
      IF(ISO.NE.0)WRITE(ISOUT,10024)I,IE,(NAME(NEL,IZ),IZ=1,4)
     <,NCTURN
10024 FORMAT(//,'  PARTICLE #',I6,' IS LOST BEFORE ELEMENT',I6,
     +  '(',4A1,')',
     +  ' DURING TURN:',I6,/,' ITS POSITION IS :',//)
      NCPART=NCPART-1
      LOGPAR(I)=.FALSE.
      WRITE(IOUT,10300)I,(PART(I,J),J=1,6)
      IF(ISO.NE.0)WRITE(ISOUT,10300)I,(PART(I,J),J=1,6)
10300 FORMAT(I4,6(E14.5))
      GO TO 65
C
C   SET UP PARTICLES AND THE ARRAY TO MULTIPLY
C
   20 CONTINUE
      IF(ICRCHK.EQ.1)CALL CSET
C     IF(MCHFLG.NE.0)WRITE(IOUT,77772)MCHFLG,IE
      IF(MCHFLG.EQ.1)CALL MSET
      X1=PART(I,1)
      X2=PART(I,2)
      X3=PART(I,3)
      X4=PART(I,4)
      X5=PART(I,5)
      X6=PART(I,6)
      IF(NT.EQ.12)GOTO 1200
      IF(NT.EQ.5)GOTO 500
      CALL COMPL
C
C   DO FAST MATRIX MULTIPLY
C
      IF(KODE(NEL).EQ.10) GO TO 40    
      IIND=NT+1
C
C   PROCESS CODES 0-9
C
      IF(NT.NE.8)GO TO 33
      IKI=DABS(ELDAT(IAD+8))+.01
      IF(MOD(NCTURN,IKI).NE.0)GO TO 60
  33     DO 30 IM=1,6
      PART(I,IM)=0.0D0
      JMF=INDEL(IM,IIND)
      IF(IIND.EQ.7)JMF=27
      DO 30 JM=1,JMF
      IF(IIND.EQ.7)then
	IJM=JM
	else
      IJM=INDORD(JM,IM,IIND)
	endif
   30 PART(I,IM)=PART(I,IM)+AMAT(N,IM,IJM)*VV(IJM)
      IF(NT.NE.8)GO TO 60
      DO 32  IK=1,6
   32 PART(I,IK)=PART(I,IK)+ELDAT(IAD+IK-1)
C
C TREAT THE CASE OF A SYNCHROTRON OSCILLATION KICK
C
c      MSYN = ELDAT(IAD+9)
c      NSYN = ABS(FLOAT(MSYN))
c      IF(MSYN.LT.0)PART(I,6)=DEL(I)*DCOS(TWOPI*NCTURN/NSYN)
	syn=eldat(iad+9)
	xsyn=abs(syn)
	if(syn.lt.0.)part(i,6)=x6+x5*xsyn  !syn=V*2*pi*f/c/E
      GO TO 60
C
C    TREAT CODE 5 OF MULTIPOLAR ELEMENT
C
  500 CALL MULTTR(X1,X2,X3,X4,X5,X6)
      PART(I,1)=X1
      PART(I,2)=X2
      PART(I,3)=X3
      PART(I,4)=X4
      PART(I,5)=X5
      PART(I,6)=X6
      GOTO 60
C
C  TREAT CODE 12 ARBITRARY ELEMENT
C
 1200 CALL TRAFCT(X1,X2,X3,X4,X5,X6,Y1,Y2,Y3,Y4,Y5,Y6,NCTURN)
      PART(I,1)=Y1
      PART(I,2)=Y2
      PART(I,3)=Y3
      PART(I,4)=Y4
      PART(I,5)=Y5
      PART(I,6)=Y6
      GOTO 60
C
C   DO CODE 10 : GENERAL MATRIX
C
   40 DO 50 IM=1,6
      PART(I,IM)=0.0D0
      DO 50 JM=1,27
      AMIJ=AMAT(N,IM,JM)
      IF (AMIJ.EQ.0D0)GOTO 50
      PART(I,IM)=PART(I,IM)+AMIJ*VV(JM)
   50 CONTINUE
   60 CONTINUE
      IF(MCHFLG.EQ.1)CALL MRESET
      IF(ICRCHK.EQ.1)CALL CRESET
      IF(ISYNGO.EQ.1)PART(I,6)=PART(I,6)-SYNDEL
   65 CONTINUE
      ISYNGO=0
C      IF(PART(1,6).NE.0.0D0)WRITE(IOUT,99999)NPART,IE,IEP,NEL,
C     >((PART(INP,JNP),JNP=1,6),INP=1,NPART)
      IF(NANAL.EQ.0)GO TO 63
      DO 64 IEN=1,NENER
      XCOR(NCTURN,IEN)=PART(IEN,1)
      XPCOR(NCTURN,IEN)=PART(IEN,2)
   64 CONTINUE
   63 IF (NJOB .EQ. 0) GO TO 62
      NCP = 1
      DO 61 IC = 1, NCASE
      XG (NCTURN,IC) = PART(NCP, 1) - XCO
      XPG(NCTURN,IC) = PART(NCP, 2) - XPCO
      IF (NJOB .EQ. 2) NCP = NCP + 1
      YG (NCTURN,IC) = PART(NCP, 3) - YCO
      YPG(NCTURN,IC) = PART(NCP, 4) - YPCO
      NCP=NCP+1
   61 CONTINUE
   62 CONTINUE
      IF(NCPART.EQ.0) RETURN
      IF(IREF.EQ.1)CALL PLPORB(IE,NEL,NELM)
      IF(MONFLG.GE.1)CALL DETLPR(IE,ILIST)
C     WRITE(IOUT,77777)MDPRT
      IF((MDPRT.EQ.-2).AND.(IFITD.NE.1))GOTO 76
      IF(MDPRT.EQ.0)GOTO75
      IF((MDPRT.LE.-1).AND.(IE.NE.NELM))GOTO76
      IF(IE.EQ.NELM)GOTO75
      CALL PRTTST(IE,ILIST,IPRT)
      IF(IPRT.NE.1)GOTO76
   75 CALL DETLPR(IE,ILIST)
   76 IF(NPRINT)100,70,90
C
C   PRINT AFTER EVERY ELEMENT
C
   70 CALL TRAKPR(0,IEP)
      GO TO 100
C
C   PRINT AFTER N TURNS AT M LOCATIONS
C
   90 IF(MOD(NCTURN,NPRINT).NE.0) GO TO 100
      CALL PRTTST(IE,ILIST,IPRT)
      IF(IPRT.NE.1)GOTO 100
      CALL TRAKPR(0,IEP)
C
C   PLOTTING AFTER EVERY ELEMENT?
C
  100	continue
	do i=1,4
	  partsave(ie,i)=part(1,i)
	  end do
c	if(ans.ne.'y'.and.ans.ne.'Y')goto 1011
	if(kod.eq.200.and.(nitx+1.ne.nit))goto 1011	!kod=200 means MOVEment
	if(kod.eq.100.or.kod.eq.110)goto 1011 !simple and least fitting
        if(kod.eq.200.or.kod.eq.300) !kod=300 TRACking
     &	  call xyplot(part(1,1),part(1,3),acleng(ie))
        xsave(ie)=part(1,1)
	ysave(ie)=part(1,3)
1011  IF(NPLOT.NE.0) GO TO 111
      NZERO = 0
      NCCUM = 1
      CALL PLOTPR(IEP,NZERO)
  111 IF(IERSET.EQ.1)CALL ERESET
      IF(MONFLG.EQ.2)GOTO112
	if(kod.ne.200)goto 110        !kod=200 means MOVEment
	if(interm.eq.0)goto 110	      !No intermediate transfer calculations
	if(interm.eq.4)goto 1111      !All intermediate transfer calcs
	if((interm.eq.1.or.interm.eq.3).and.kode(nel).eq.1)goto 1111  !Dipoles
        if((interm.eq.2.or.interm.eq.3).and.kode(nel).eq.2)goto 1111  !Quads
 	if(interm.ge.5.and.(nel.eq.norlst(ie+1)))goto 1111 	!symmetry point
	if(interm.ne.6)goto 114		!Horizontal separators	
	  if(name(norlst(ie),1).ne.lab(1)
     &	      .and.name(norlst(ie),1).ne.lab(1))goto 114
	  if(name(norlst(ie),3).eq.lab(2))goto 1111
	  if(name(norlst(ie),2).eq.lab(3)
     &		.and.name(norlst(ie),3).eq.lab(4))goto 1111
  114	if(ie.eq.itrend)goto 1111	!End of line
	if(interm.gt.0)goto 110
	do 113 jj=1,iabs(interm)
	  if(ie.eq.nesave(jj))goto 110
  113	continue
	goto 110
 1111	call transave(ie)
  110 CONTINUE
c save max amplitude
	if(cbetax.ne.0.and.cbetay.ne.0.)then
	    xamp=part(1,1)**2*(1.+calphx**2)/cbetax
     &	      +2.*part(1,2)*part(1,1)*calphx+part(1,2)**2*cbetax
	    yamp=part(1,3)**2*(1.+calphy**2)/cbetay
     &	      +2.*part(1,3)*part(1,4)*calphy+part(1,4)**2*cbetay
	xampmax=dmax1(xamp,xampmax) !max x amplitude for particle 1
	yampmax=dmax1(yamp,yampmax) !max y amplitude for particle 2
	endif
c
	if(kod.eq.200)then     !kod=200 means MOVEment
	write(iout,1101)detxmax,detymax,scmax
1101	format(1x,
     &	  'Max detx =',e13.4,' Max dety =',e13.4,' Max coup =',e13.4)
	write(iout,1102)kspotmax,detxymax
1102	format(1x,' Max det is at ',i3,' Max det =',e13.4)
	endif
c     	  call bastart
  112 IF(IALFLG.EQ.1)ISDBEG=IXS
      IF(IALFLG.EQ.1)CALL DETLPR(ITREND,ILIST)
      IF(IALFLG.EQ.1)IMSBEG=IMONSD
C
C   CHECK THE PLOT REQUESTED FOR THIS TURN
C
      IF(NCTURN.LT.NTURN) GO TO 300
      IF(NPRINT.NE.-2)CALL TRAKPR(0,NELM)
  300 IF(NPLOT.LE.0) GO TO 200
      IF(NCTURN.NE.NXREQ) GO TO 200
      NZERO = 1
      IF(NCTURN.EQ.NPLOT) NZERO = 0
      IF(NCCUM.EQ.1) NZERO = 0
      IF(NCTURN.EQ.LSTREQ) NCCUM = 1
      CALL PLOTPR(NELM,NZERO)
      NXREQ = NXREQ+NPLOT
  200 IF(NCTURN.LT.NTURN) GO TO 10
cdlr
	misflg=-1
cdlr
      RETURN
      END
C     ***********************
c      SUBROUTINE TRAFCT(XI,XPI,YI,YPI,ALI,DELI,
c     > XO,XPO,YO,YPO,ALO,DELO,NTURN)
C     **********************
c        IMPLICIT REAL*8 (A-H,O-Z)
C
C
C THIS ROUTINE IS TO BE PROGRAMMED BY THE USER
C THE INPUT COORDINATES HAVE A TRAILING I AND THE OUTPUT
C COORDINATES HAVE A TRAILING O.
C IF THE TRANSFORM THAT IS PROGRAMMED HERE IS NOT SYMPLECTIC
C OR DOES NOT SATISFY LIOUVILLE'S THEOREM SOME OTHER PARTS OF
C THE PROGRAM MAY NOT GIVE MEANINGFUL RESULTS !!!!!
C THE ONLY PARTS OF THE PROGRAMME AFFECTED BY THIS ELEMENT ARE
C THE MOVE(MENT) ANALYSIS OPERATION, THE  TRAC(KING) OPERATION AND
C THE  MODI(FICATION) OF PARAMETERS OPERATION.
C NTURN IN THE CURRENT TURN BEING PROCESSED AND CAN BE USED
C TO DESIGN AN ELEMENT WHOSE BEHAVIOUR VARIES WITH THE TURN.
C
C
c      COMMON/ARB/PARA(20)
c	xo=xi
c	xpo=xpi
c      ALO=ALI
c      DELO=DELI
c      YO=YI
c      YPO=YPI
c      DELO=DELI
c      ALO=ALI
c	if(para(4).lt.3.)then
c	xl=para(1)
c	rho=para(2)
c	e=para(3)
c	rhoenergy=rho*(deli+1)
c	thetaenergy=xl/rhoenergy
c	theta0=xl/rho
c	xpo=xpi-thetaenergy+theta0
c	xo=xi+(1.-cos(theta0))*rho-(1.-cos(thetaenergy))*rhoenergy
c	endif
c	if(para(4).gt.3.)then
c	xl=para(1)
c	xk0=para(2)
c	e=para(3)
c	xk2=xk0*1./(1.+deli)
c	if(xk2.gt.0.)then
c	  xk=sqrt(xk2)
c	  arg=xk*xl
c	  r11=cos(arg)
c	  r12=1./xk*sin(arg)
c	  r21=-xk2*r12
c	  r22=r11
c	  r33=cosh(arg)
c	  r34=1./xk*sinh(arg)
c	  r43=xk2*r34
c	  r44=r33
c	else
c	  xk=sqrt(-xk2)
c	  arg=xk*xl
c	  r11=cosh(arg)
c	  r12=1./xk*sinh(arg)
c	  r21=xk*xk*r12
c	  r22=r11
c	  r33=cos(arg)
c	  r34=1./xk*sin(arg)
c	  r43=-xk*xk*r34
c	  r44=r33
c	endif
c	xo=xi*r11+xpi*r12
c	xpo=xi*r21+xpi*r22
c	yo=yi*r33+ypi*r34
c	ypo=yi*+ypi*r44
c	endif
c     RETURN
c     END
C     ************************
      SUBROUTINE TRAKPR(ICODE,IE)
C     ************************
      IMPLICIT REAL*8 (A-H,O-Z)
	include 'input.inc'
	include 'elements.inc'
	include 'track.inc'
	character*1 geo
	data ngeomfile/0/
C
C   WHAT KIND OF PRINTING?
C
        NEL = NORLST(IE)
      IF(ICODE)10,30,30
C
C   INITIAL POSITIONS PRINTING
C
10    WRITE(IOUT,10301)
10301 FORMAT(//,' INITIAL POSITIONS OF PARTICLES ',/)
c
	if(ngeomfile.ge.1)close(unit=17)
	      ngeomfile=ngeomfile+1
	      encode(1,2114,geo)ngeomfile
	type *,'ngeomfile,geo ',ngeomfile,geo
2114	      format(i1)
	call strip_dirname(filename,istart,iend,leftbrak,ightbrak)
	    type 2112,filename(istart:iend)//'_geom_'//geo
2112	    format(1x,' Open file ',a,
     &	      ' to write turn by turn tracking data.')
	    open(unit=17,file=filename(istart:iend)//'_geom_'//geo,
     &	      type='new')
	write(17,2113)
2113	format(6x,'x',12x,'xp',12x,'y',12x,'yp')
c
      DO 20 I=1,NPART
      IF(.NOT.LOGPAR(I)) GO TO 20
      WRITE(IOUT,10300)I,(PART(I,K),K=1,6)
10300 FORMAT(I4,6(E14.5))
c
	write(17,10021)i,(part(i,k),k=1,4)
c
20    CONTINUE
      RETURN
C
C   OTHER PRINTING
C
30    WRITE(IOUT,10302)IE,(NAME(NEL,IZ),IZ=1,4),NCTURN
10302 FORMAT(/,' PARTICLE POSITIONS AFTER ELEMENT',2X,I4,'(',4A1,')',
     >  2X,'DURING TURN',I6,/)
      DO 40 I=1,NPART
      IF(.NOT.LOGPAR(I)) GO TO 40
      WRITE(IOUT,10300) I,(PART(I,K),K=1,6)
c
	write(17,10021)I,(part(I,K),K=1,4)
10021	format(1x,i,4e13.3)
c
40    CONTINUE
      RETURN
      END
C     ***********************
      SUBROUTINE TWISS
C     ***********************
      IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
	include 'input.inc'
	include 'elements.inc'
      COMMON/CONST/PI,TWOPI,CRDEG,CMAGEN,CLIGHT,EMASS,ERAD,ECHG
      COMMON/PRODCT/KODEPR,NEL,NOF
      IAD=IADR(N)
      N=MADR(N)
C
C     CLEAR AMAT(N,6,27)
C
      DO 10 I=1,6
      DO 10 J=1,NOF
      AMAT(N,I,J)=0.0D0
      IF(I.EQ.J) AMAT(N,I,J)=1.0D0
10    CONTINUE
      ANG=ELDAT(IAD)*CRDEG
      CH=DCOS(ANG)
      SH=DSIN(ANG)
      ANG=ELDAT(IAD+3)*CRDEG
      CV=DCOS(ANG)
      SV=DSIN(ANG)
      AMAT(N,1,1)=CH+ELDAT(IAD+2)*SH
      AMAT(N,1,2)=SH*ELDAT(IAD+1)
      AMAT(N,2,1)=-(1.0+ELDAT(IAD+2)*ELDAT(IAD+2))*SH/ELDAT(IAD+1)
      AMAT(N,2,2)=CH-ELDAT(IAD+2)*SH
      AMAT(N,3,3)=CV+ELDAT(IAD+5)*SV
      AMAT(N,3,4)=SV*ELDAT(IAD+4)
      AMAT(N,4,3)=-(1.0D0+ELDAT(IAD+5)*ELDAT(IAD+5))*SV/ELDAT(IAD+4)
      AMAT(N,4,4)=CV-ELDAT(IAD+5)*SV
      RETURN
      END
C  **********************
c      FUNCTION URAND(IX)
cC  **********************
c      IMPLICIT REAL*8(A-H,O-Z)
c      DATA M2/0/,ITWO/2/
c      IF(M2.NE.0)GOTO 20
c      M=1
c   10 M2=M
c      M=ITWO*M2
c      IF(M.GT.M2)GOTO 10
c      HALFM=M2
c      IA=8*IDINT(HALFM*DATAN(1.0D0)/8.0D0) + 5
c      IC=2*IDINT(HALFM*(0.5D0-DSQRT(3.0D0)/6.0D0))+1
c      S=0.5D0/HALFM
c   20 IX=IX*IA+IC
c      IF(IX/2.GT.M2)IX=(IX-M2)-M2
c      IF(IX.LT.0)IX=(IX+M2)+M2
c      URAND=DFLOAT(IX)*S
c      RETURN
c      END
C  **********************
c      SUBROUTINE VARY(JV,NVPAR,VARVAL)
cC  **********************
c      IMPLICIT REAL*8 (A-H,O-Z)
c	include 'input.inc'
c	include 'elements.inc'
c      NV = IABS(JV)
c	if(jv.lt.0.and.NAME(nv,1).eq.'L'.and.NVPAR.eq.2)then  !skew quad
c	  call skewcal(nv,varval)
c	  endif
c      NVEL = IADR(NV)
c      ELDAT(NVEL+NVPAR-1) = VARVAL
cC
cC     RECOMPUTE MATRIX
cC
c      CALL MATGEN(NV)
c      RETURN
c      END
C     ***********************
      SUBROUTINE VBEND
C     ***********************
      IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
	include 'elements.inc'
	include 'tempstor.inc'
      COMMON/PRODCT/KODEPR,NEL,NOF
      COMMON/BND/IAD
C
C     GET ADDRESS OF THE PARAMETERS AND THE MATRIX
C
      IAD = IADR(N)
      N = MADR(N)
C
C     CLEAR AMAT(N,6,27)
C
      DO 100 I=1,6
      DO 100 J=1,NOF
      AMAT(N,I,J)=0.0D0
 100  CONTINUE
C
C     SET UP THE 90 DEGREE KICK MATRIX IN TEMP
C
      DO 10 IX = 1,6
      DO 10 IY = 1,NOF
      TEMP(IX,IY)=0.0D0
10    CONTINUE
      TEMP(1,3)=1
      TEMP(2,4)=1.0D0
      TEMP(3,1)=-1.0D0
      TEMP(4,2)=-1.0D0
      TEMP(5,5)=1.0D0
      TEMP(6,6)=1.0D0
      CALL BEND
      KODEPR=8
C
C  SET -90 DEGREE KICK MATRIX IN AMAT
C
      DO 200 I=1,6
      DO 200 J=1,NOF
      AMAT(N,I,J)=0.0D0
 200  CONTINUE
      AMAT(N,1,3)=-1.0D0
      AMAT(N,2,4)=-1.0D0
      AMAT(N,3,1)=1.0D0
      AMAT(N,4,2)=1.0D0
      AMAT(N,5,5)=1.0D0
      AMAT(N,6,6)=1.0D0
      CALL PROMAT
C
C     PUT TEMP INTO AMAT(N,6,27)
C
      DO 20 IX=1,6
      DO 20 IY=1,NOF
  20  AMAT(N,IX,IY) = TEMP(IX,IY)
      RETURN
      END
C     ***********************
C     SPECIAL FUNCTIONS
C     ***********************
      FUNCTION DSINH(X)
      IMPLICIT REAL*8(A-H,O-Z)
      DSINH=(DEXP(X)-DEXP(-X))/2.0D0
      RETURN
      END
      FUNCTION DCOSH(X)
      IMPLICIT REAL*8(A-H,O-Z)
      DCOSH=(DEXP(X)+DEXP(-X))/2.0D0
      RETURN
      END
      FUNCTION DTAN(X)
      IMPLICIT REAL*8(A-H,O-Z)
      DTAN=DSIN(X)/DCOS(X)
      RETURN
      END
