Commit b21f4151 authored by Thomas Planche's avatar Thomas Planche
Browse files

committing zgoubi source files

parent a3de9fd2
C ZGOUBI, a program for computing the trajectories of charged particles
C in electric and magnetic fields
C Copyright (C) 1988-2007 Franois Mot
C
C This program is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation; either version 2 of the License, or
C (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this program; if not, write to the Free Software
C Foundation, Inc., 51 Franklin Street, Fifth Floor,
C Boston, MA 02110-1301 USA
C
C Franois Mot <fmeot@bnl.gov>
C Brookhaven National Laboratory
C C-AD, Bldg 911
C Upton, NY, 11973
C -------
SUBROUTINE ACCEN(JJ,
> RATIN)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C--------------------------------------------------------------------------
C Computes the rms ellipse with emittance AY matched to a particle set,
C and gives particle rate within that ellipse
C--------------------------------------------------------------------------
INCLUDE "C.CDF.H" ! COMMON/CDF/ IES,LF,LST,NDAT,NRES,NPLT,NFAI,NMAP,NSPN,NLOG
INCLUDE "C.CONST2.H" ! COMMON/CONST2/ ZERO, UN
INCLUDE 'MXLD.H'
INCLUDE "C.DON.H" ! COMMON/DON/ A(MXL,MXD),IQ(MXL),IP(MXL),NB,NOEL
INCLUDE "MAXCOO.H"
INCLUDE "MAXTRA.H"
LOGICAL AMQLU(5),PABSLU
INCLUDE "C.FAISC.H" ! COMMON/FAISC/ F(MXJ,MXT),AMQ(5,MXT),DP0(MXT),IMAX,IEX(MXT),
C $ IREP(MXT),AMQLU,PABSLU
INCLUDE "C.UNITS.H" ! COMMON/UNITS/ UNIT(MXJ)
SAVE EPSPI
Compute rms ellipse
CALL LPSFIT(JJ,
> EMIT,ALP,BET,XM,XPM)
Compute number of particles alive and number inside ellipse
CALL CNTINL(JJ,EPSPI,ALP,BET,XM,XPM,
> NLIV,NINL)
C GAM = (1.D0+ALP*ALP) / BET
CC------- JJ = 1, 2 or 3
C J1=2*JJ
C J2=J1+1
C IF(J1.EQ.6) THEN
C J1=1
C J2=6
C ENDIF
C NINW=0
C NOFF = 0
C DO 1 I=1,IMAX
C IF(IEX(I) .LT. -1) THEN
C NOFF = NOFF+1
C GOTO 1
C ENDIF
C Y2 = F(J1,I)*UNIT(J1-1)
C T2 = F(J2,I)*UNIT(J2-1)
C YT = Y2*T2
C Y2 = Y2*Y2
C T2 = T2*T2
C IF( GAM*Y2+2.D0*ALP*YT+BET*T2 .LE. EPSPI) NINW=NINW+1
C 1 CONTINUE
C----- RATIN is of course at maximum 1.
C RATIN = DBLE(NINW)/DBLE(IMAX-NOFF)
C RATIN = DBLE(NINL)/DBLE(NLIV)
RATIN = DBLE(NINL)/DBLE(IMAX)
RETURN
ENTRY ACCENW(EPSPIN)
EPSPI = EPSPIN
RETURN
END
C ZGOUBI, a program for computing the trajectories of charged particles
C in electric and magnetic fields
C Copyright (C) 1988-2007 Franois Mot
C
C This program is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation; either version 2 of the License, or
C (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this program; if not, write to the Free Software
C Foundation, Inc., 51 Franklin Street, Fifth Floor,
C Boston, MA 02110-1301 USA
C
C Franois Mot <fmeot@bnl.gov>
C Brookhaven National Laboratory
C C-AD, Bldg 911
C Upton, NY, 11973
C -------
SUBROUTINE ACCEP(JJ,
> EMIT,AMA,BMA,XMA,XPMA,NLIV,MXINL)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C----------------------------------------------------------------------------
C Looks for ellipse with emittance EPSPI and other parameters taken in
C the vicinity of the rms ellipse, that has the maximum particle content.
C----------------------------------------------------------------------------
INCLUDE "C.CDF.H" ! COMMON/CDF/ IES,LF,LST,NDAT,NRES,NPLT,NFAI,NMAP,NSPN,NLOG
INCLUDE "C.CONST2.H" ! COMMON/CONST2/ ZERO, UN
INCLUDE 'MXLD.H'
INCLUDE "C.DON.H" ! COMMON/DON/ A(MXL,MXD),IQ(MXL),IP(MXL),NB,NOEL
INCLUDE "MAXCOO.H"
INCLUDE "MAXTRA.H"
LOGICAL AMQLU(5),PABSLU
INCLUDE "C.FAISC.H" ! COMMON/FAISC/ F(MXJ,MXT),AMQ(5,MXT),DP0(MXT),IMAX,IEX(MXT),
C $ IREP(MXT),AMQLU,PABSLU
INCLUDE "C.UNITS.H" ! COMMON/UNITS/ UNIT(MXJ)
PARAMETER (RANGE=0.2D0,ISTP=3)
SAVE EPSPI
Compute rms ellipse
CALL LPSFIT(JJ,
> EMIT,ALP,BET,XM,XPM)
Compute number of particles alive and number inside ellipse
DX = XM*RANGE/ISTP
IF(DX.EQ.0.D0) DX = 1.D-10
DXP = XPM*RANGE/ISTP
IF(DXP.EQ.0.D0) DXP = 1.D-10
DA = ALP*(2.D0*RANGE)/ISTP
IF(DA.EQ.0.D0) DA = 1.D-10
DB = BET*(2.D0*RANGE)/ISTP
IF(DB.EQ.0.D0) DB = 1.D-10
MXINL = -1
DO 10 IXX = -ISTP,ISTP
XX = XM + IXX*DX
DO 10 IXXP = -ISTP,ISTP
XXP = XPM + IXXP*DXP
DO 10 IAA = -ISTP,ISTP
AA = ALP + IAA*DA
DO 10 IBB = -ISTP,ISTP
BB = BET + IBB*DB
CALL CNTINL(JJ,EPSPI,AA,BB,XX,XXP,
> NLIV,NINL)
IF(NINL .GT. MXINL) THEN
MXINL = NINL
AMA = AA
BMA = BB
XMA = XX
XPMA = XXP
ENDIF
10 CONTINUE
C RATIN = DBLE(MXINL)/DBLE(IMAX)
EMIT = EPSPI
C write(*,*) ' sbr accep nliv, ninl :', nliv, ninl
RETURN
ENTRY ACCEPW(EPSPIN)
EPSPI = EPSPIN
RETURN
END
C ZGOUBI, a program for computing the trajectories of charged particles
C in electric and magnetic fields
C Copyright (C) 1988-2007 Franois Mot
C
C This program is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation; either version 2 of the License, or
C (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this program; if not, write to the Free Software
C Foundation, Inc., 51 Franklin Street, Fifth Floor,
C Boston, MA 02110-1301 USA
C
C Franois Mot <fmeot@bnl.gov>
C Brookhaven National Laboratory
C C-AD, Bldg 911
C Upton, NY, 11973
C -------
SUBROUTINE ADPOL(ID,
> IK,B,DB,DDB,D3BX,D3BY,D4BX,D4BY,BT)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C--------------------------------------------------------------------
C Add up pole components of multipole lenses. Can be E or B field
C--------------------------------------------------------------------
DIMENSION B(5,3),DB(3,3),DDB(3,3,3)
DIMENSION D3BX(3,3,3), D3BY(3,3,3)
DIMENSION D4BX(3,3,3,3) ,D4BY(3,3,3,3)
DIMENSION BT(5,15)
GOTO(1,2,3) IK
1 CONTINUE
IF(ID.GE.0) THEN
BT(1,1)= B(1,1)
BT(1,2)= B(1,2)
BT(1,3)= B(1,3)
IF(ID.GE.1) THEN
BT(2,1)= DB(1,1)
BT(2,2)= DB(2,1)
BT(2,3)= DB(3,1)
BT(2,4)= DB(2,2)
BT(2,5)= DB(3,2)
IF(ID.GE.2) THEN
BT(3,1)= DDB(1,1,1)
BT(3,2)= DDB(2,1,1)
BT(3,3)= DDB(3,1,1)
BT(3,4)= DDB(2,2,1)
BT(3,5)= DDB(3,2,1)
BT(3,7)= DDB(2,2,2)
BT(3,8)= DDB(3,2,2)
IF(ID.GE.3) THEN
BT(4,1)= D3BX(1,1,1)
BT(4,2)= D3BX(2,1,1)
BT(4,3)= D3BX(3,1,1)
BT(4,4)= D3BX(2,2,1)
BT(4,5)= D3BX(3,2,1)
BT(4,7)= D3BX(2,2,2)
BT(4,8)= D3BX(3,2,2)
BT(4,11)= D3BY(2,2,2)
BT(4,12)= D3BY(3,2,2)
IF(ID.GE.4) THEN
BT(5,1)= D4BX(3,1,1,1)
BT(5,2)= D4BX(3,2,1,1)
BT(5,3)= D4BX(3,2,2,1)
BT(5,4)= D4BX(3,3,3,1)
BT(5,5)= D4BX(3,2,2,2)
BT(5,6)= D4BX(3,3,3,2)
BT(5,7)= D4BY(3,2,2,2)
BT(5,8)= D4BY(3,3,3,2)
ENDIF
C ID GE 4
ENDIF
C ID GE 3
ENDIF
C ID GE 2
ENDIF
C ID GE 1
ENDIF
C ID GE 0
C FM june 2015. Already in chamc CALL RAZDRV(KFL)
IK=2
GOTO 99
2 CONTINUE
IF(ID.GE.0) THEN
BT(1,1)=BT(1,1) + B(1,1)
BT(1,2)=BT(1,2) + B(1,2)
BT(1,3)=BT(1,3) + B(1,3)
IF(ID.GE.1) THEN
BT(2,1)=BT(2,1) + DB(1,1)
BT(2,2)=BT(2,2) + DB(2,1)
BT(2,3)=BT(2,3) + DB(3,1)
BT(2,4)=BT(2,4) + DB(2,2)
BT(2,5)=BT(2,5) + DB(3,2)
IF(ID.GE.2) THEN
BT(3,1)=BT(3,1) + DDB(1,1,1)
BT(3,2)=BT(3,2) + DDB(2,1,1)
BT(3,3)=BT(3,3) + DDB(3,1,1)
BT(3,4)=BT(3,4) + DDB(2,2,1)
BT(3,5)=BT(3,5) + DDB(3,2,1)
BT(3,7)=BT(3,7) + DDB(2,2,2)
BT(3,8)=BT(3,8) + DDB(3,2,2)
IF(ID.GE.3) THEN
BT(4,1)=BT(4,1) + D3BX(1,1,1)
BT(4,2)=BT(4,2) + D3BX(2,1,1)
BT(4,3)=BT(4,3) + D3BX(3,1,1)
BT(4,4)=BT(4,4) + D3BX(2,2,1)
BT(4,5)=BT(4,5) + D3BX(3,2,1)
BT(4,7)=BT(4,7) + D3BX(2,2,2)
BT(4,8)=BT(4,8) + D3BX(3,2,2)
BT(4,11)=BT(4,11) + D3BY(2,2,2)
BT(4,12)=BT(4,12) + D3BY(3,2,2)
IF(ID.GE.4) THEN
BT(5,1)=BT(5,1) + D4BX(3,1,1,1)
BT(5,2)=BT(5,2) + D4BX(3,2,1,1)
BT(5,3)=BT(5,3) + D4BX(3,2,2,1)
BT(5,4)=BT(5,4) + D4BX(3,3,3,1)
BT(5,5)=BT(5,5) + D4BX(3,2,2,2)
BT(5,6)=BT(5,6) + D4BX(3,3,3,2)
BT(5,7)=BT(5,7) + D4BY(3,2,2,2)
BT(5,8)=BT(5,8) + D4BY(3,3,3,2)
ENDIF
C ID GE 4
ENDIF
C ID GE 3
ENDIF
C ID GE 2
ENDIF
C ID GE 1
ENDIF
C ID GE 0
C FM june 2015. Already in chamc CALL RAZDRV(KFL)
GOTO 99
3 CONTINUE
IF(ID.GE.0) THEN
B(1,1)=BT(1,1)
B(1,2)=BT(1,2)
B(1,3)=BT(1,3)
IF(ID.GE.1) THEN
DB(1,1)=BT(2,1)
DB(2,1)=BT(2,2)
DB(3,1)=BT(2,3)
DB(2,2)=BT(2,4)
DB(3,2)=BT(2,5)
IF(ID.GE.2) THEN
DDB(1,1,1)=BT(3,1)
DDB(2,1,1)=BT(3,2)
DDB(3,1,1)=BT(3,3)
DDB(2,2,1)=BT(3,4)
DDB(3,2,1)=BT(3,5)
DDB(2,2,2)=BT(3,7)
DDB(3,2,2)=BT(3,8)
c write(89,*) b(1,3),DB(3,1),DB(3,2),DDB(3,2,2),DDB(3,1,1),
c > ' sbr adpol '
c write(89,*) DB(3,1), DB(3,2), DB(1,1), DB(2,1),' sbr adpol '
IF(ID.GE.3) THEN
D3BX(1,1,1)=BT(4,1)
D3BX(2,1,1)=BT(4,2)
D3BX(3,1,1)=BT(4,3)
D3BX(2,2,1)=BT(4,4)
D3BX(3,2,1)=BT(4,5)
D3BX(2,2,2)=BT(4,7)
D3BX(3,2,2)=BT(4,8)
D3BY(2,2,2)=BT(4,11)
D3BY(3,2,2)=BT(4,12)
IF(ID.GE.4) THEN
D4BX(3,1,1,1)=BT(5,1)
D4BX(3,2,1,1)=BT(5,2)
D4BX(3,2,2,1)=BT(5,3)
D4BX(3,3,3,1)=BT(5,4)
D4BX(3,2,2,2)=BT(5,5)
D4BX(3,3,3,2)=BT(5,6)
D4BY(3,2,2,2)=BT(5,7)
D4BY(3,3,3,2)=BT(5,8)
ENDIF
C ID GE 4
ENDIF
C ID GE 3
ENDIF
C ID GE 2
ENDIF
C ID GE 1
ENDIF
C ID GE 0
99 CONTINUE
RETURN
END
C ZGOUBI, a program for computing the trajectories of charged particles
C in electric and magnetic fields
C Copyright (C) 1988-2007 François Méot
C
C This program is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation; either version 2 of the License, or
C (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this program; if not, write to the Free Software
C Foundation, Inc., 51 Franklin Street, Fifth Floor,
C Boston, MA 02110-1301 USA
C
C François Meot <fmeot@bnl.gov>
C Brookhaven National Laboratory
C C-AD, Bldg 911
C Upton, NY, 11973
C -------
SUBROUTINE AGSBLW(MOD2,NOEL,DEV,NBLW,NB,WN,WA,
> BM)
C ----------------------------------------------
C NBLW : # of blwg ; WN(i=1,nblw) : # of turns ;
C WA(i=1,nblw) : amperes.
C ----------------------------------------------
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION WN(*),WA(*),BM(NB)
INCLUDE 'MXFS.H'
INCLUDE 'MXLD.H'
PARAMETER (LBLSIZ=20)
CHARACTER(LBLSIZ) LABEL
INCLUDE "C.LABEL.H" ! COMMON/LABEL/ LABEL(MXL,2)
PARAMETER (KSIZ=10)
CHARACTER(KSIZ) FAM
CHARACTER(LBLSIZ) LBF
INCLUDE "C.SCALT.H" ! COMMON/SCALT/ FAM(MXF),LBF(MXF,MLF)
CHARACTER(LBLSIZ) MMNM
PARAMETER(TURNMM=16.D0)
INTEGER DEBSTR, FINSTR
MMNM = LABEL(NOEL,1)
MMNM = MMNM(DEBSTR(MMNM):FINSTR(MMNM))
IF (MOD2 .EQ.0) THEN
C User defined values.
ELSEIF(MOD2 .EQ.1) THEN
C Actual AGS values
c/========================================================================
c/ Cold Snake a16-19_blw bump. Backlegs on A16, A17, -A18, -A19 in series.
c/========================================================================
C WN(2) = 0.D0
IF (MMNM .EQ. 'MM_A16AD') THEN
NBLW = 1
SIGN = 1.D0
NWL = 10
WN(1) = DBLE(NWL) * SIGN
ELSEIF(MMNM .EQ. 'MM_A17CF') THEN
NBLW = 1
SIGN = 1.D0
NWL = 10
WN(1) = DBLE(NWL) * SIGN
ELSEIF(MMNM .EQ. 'MM_A18CF') THEN
NBLW = 1
SIGN = -1.D0
NWL = 10
WN(1) = DBLE(NWL) * SIGN
ELSEIF(MMNM .EQ. 'MM_A19BD') THEN
NBLW = 1
SIGN = -1.D0
NWL = 12
WN(1) = DBLE(NWL) * SIGN
C/================================================================
C/ COLD SNAKE A20 BLW. BACKLEG WIRED ONLY on A20, 12 turns.
c/ Note: There are three power supplies wired on A20, but the controls
c/ system (currently) combines them into one total current which
c/ we use here.
c/================================================================
ELSEIF(MMNM .EQ. 'MM_A20BD') THEN
NBLW = 1
SIGN = 1.D0
NWL = 12
WN(1) = DBLE(NWL) * SIGN
c/====================================================================
c/ Cold Snake b2-5_blw bump. Backlegs wired on B2, B3, -B4, -B5 in series
c/====================================================================
c B02 magnet has 2 windings
ELSEIF(MMNM .EQ. ' MM_B02BF') THEN
NBLW = 2
SIGN = 1.D0
NWL = 12
WN(1) = DBLE(NWL) * SIGN
SIGN = 1.D0
NWL = 6
WN(2) = DBLE(NWL) * SIGN
ELSEIF(MMNM .EQ. 'MM_B03CD') THEN
NBLW = 1
SIGN = 1.D0
NWL = 10
ELSEIF(MMNM .EQ. 'MM_B04CD ') THEN
NBLW = 1
SIGN = -1.D0
NWL = 10
ELSEIF(MMNM .EQ. 'MM_B05AF') THEN
NBLW = 1
SIGN = -1.D0
NWL = 10
c/====================================================================
c/ L20 Position Bump
c/ BLWs on K19, K20, L13, L14, A7, A8, B1, B2
c/ with 6, 6, -5, -5, -5, -5, 6, 6 turns respectively.
c/====================================================================
ELSEIF(MMNM .EQ. 'MM_K19BD ') THEN
NBLW = 1
SIGN = 1.D0
NWL = 6
WN(1) = DBLE(NWL) * SIGN
ELSEIF(MMNM .EQ. 'MM_K20BD') THEN
NBLW = 1
SIGN = 1.D0
NWL = 6
WN(1) = DBLE(NWL) * SIGN
ELSEIF(MMNM .EQ. 'MM_L13CF ') THEN
NBLW = 1
SIGN = -1.D0
NWL = 5
WN(1) = DBLE(NWL) * SIGN
ELSEIF(MMNM .EQ. 'MM_L14CF') THEN
NBLW = 1
SIGN = -1.D0
NWL = 5
WN(1) = DBLE(NWL) * SIGN
ELSEIF(MMNM .EQ. 'MM_A07CD ') THEN
NBLW = 1
SIGN = -1.D0
NWL = 5
WN(1) = DBLE(NWL) * SIGN
ELSEIF(MMNM .EQ. 'MM_A08CD') THEN
NBLW = 1
SIGN = -1.D0
NWL = 5
WN(1) = DBLE(NWL) * SIGN
ELSEIF(MMNM .EQ. 'MM_B01BF') THEN
NBLW = 1
SIGN = 1.D0
NWL = 6
WN(1) = DBLE(NWL) * SIGN
C BLIL20=BLI
c/===============================================================
c/ L20 Angle Bump (1 lambda)
c/ BLW windings around L6, L7, A14, A15
c/ with 5, 5, -5, -5 turns, respectively.
c/===============================================================
ELSEIF(MMNM .EQ. 'MM_L06AF') THEN
NBLW = 1
SIGN = 1.