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(400) 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 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(400) 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, 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, & monposrev(2500) 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,2),part(1,3),part(1,4),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