*   Plots longitudinal phase space trajectories
*   of particles in an accelerator

*   Usage: PAW> exec long_dynamics

*   Author: Giulio Stancari (stancari@fnal.gov)
*   Date: 26 February 2000


MACRO MAIN

SET *
OPT *

MESSAGE
MESSAGE '===== PARTICLE LONGITUDINAL DYNAMICS ====='
MESSAGE

* guess what THIS is
IF ($VEXIST(PI).NE.0) THEN
  V/DEL PI
ENDIF
V/CREA PI(1) R 3.141592654

MESSAGE 'PARTICLE DATA'

IF ($VEXIST(MP).NE.0) THEN
  V/DEL MP
ENDIF
MP=938.272
READ MP 'Rest energy [MeV]'
V/CREA MP(1) R [MP]

IF ($VEXIST(C).NE.0) THEN
  V/DEL C
ENDIF
C=1.
READ C 'Charge in units of |e|'
V/CREA C(1) R [C]

IF ($VEXIST(ES).NE.0) THEN
  V/DEL ES
ENDIF
ES=8938.272
READ ES 'Total energy [MeV]'
V/CREA ES(1) R [ES]

MESSAGE
MESSAGE 'RF AND MACHINE DATA'

IF ($VEXIST(V).NE.0) THEN
  V/DEL V
ENDIF
V=5.
READ V 'Voltage [MV]'
V/CREA V(1) R [V]

IF ($VEXIST(EV).NE.0) THEN
  V/DEL EV
ENDIF
V/CREA EV(1) R
SIGMA EV=[C]*V

IF ($VEXIST(F).NE.0) THEN
  V/DEL F
ENDIF
F=0.053
READ F 'Frequency [GHz]'
V/CREA F(1) R [F]

IF ($VEXIST(H).NE.0) THEN
  V/DEL H
ENDIF
H=90.
READ H 'Harmonic number'
V/CREA H(1) R [H]

IF ($VEXIST(PHIS).NE.0) THEN
  V/DEL PHIS
ENDIF
INPUT:
PHIS=3.1416
READ PHIS 'Synchronous phase [0, 2*pi) [rad]'
IF ( ([PHIS].LT.0.).OR.($SIGMA(2.*PI).LE.[PHIS]) ) THEN
  MESSAGE ' *** Invalid value '//[PHIS]
  GOTO INPUT
ENDIF
V/CREA PHIS(1) R [PHIS]

IF ($VEXIST(ETA).NE.0) THEN
  V/DEL ETA
ENDIF
ETA=0.006
READ ETA 'Slip factor eta'
V/CREA ETA(1) R [ETA]

IF ($VEXIST(DENEW).NE.0) THEN
  V/DEL DENEW
ENDIF
V/CREA DENEW(1)

IF ($VEXIST(PHINEW).NE.0) THEN
  V/DEL PHINEW
ENDIF
V/CREA PHINEW(1)

* useful quantities
SIGMA GAMMA = ES/MP
SIGMA PS = SQRT( ES*ES - MP*MP )
SIGMA BETA = PS/ES
SIGMA GAMMAT = 1. / SQRT( ETA + 1./GAMMA/GAMMA )

* synchrotron tune
IF $SIGMA((ETA*COS(PHIS))).GT.0 THEN
  MESSAGE
  MESSAGE ' *** Warning: Motion is unstable!'
  MESSAGE
ENDIF
SIGMA NU=SQRT(ABS(ETA*H*EV*COS(PHIS)/2./PI/BETA/BETA/ES))
* bucket width in ns
SIGMA N=INT(ABS(PHIS)/PI)
SIGMA DPHI = MIN(ABS(PHIS-PI/2.),ABS(PHIS-3./2.*PI))
SIGMA DT = 2.*DPHI/2./PI/F
* this is the stationary bucket area (phi=0)
SIGMA A0 = 8.*BETA/PI/F*SQRT(EV*ES/2./PI/H/ABS(ETA))
EXEC INTEGR
* this is the stationary bucket area
SIGMA A = A0 * I($SIGMA(INT(NN(1)))) / I0($SIGMA(INT(NN0(1))))

MESSAGE
MESSAGE 'KINEMATIC AND DYNAMIC QUANTITIES'
MESSAGE 'Charge in units of |e| = '//$FORMAT(C(1),G8.6)
MESSAGE 'Rest energy [MeV] =      '//$FORMAT(MP(1),G8.6)
MESSAGE 'Relativistic beta =      '//$FORMAT(BETA(1),G8.6)
MESSAGE 'Relativistic gamma =     '//$FORMAT(GAMMA(1),G8.6)
MESSAGE 'Momentum [MeV] =         '//$FORMAT(PS(1),G8.6)
MESSAGE 'Energy [MeV] =           '//$FORMAT(ES(1),G8.6)

MESSAGE
MESSAGE 'RF AND MACHINE QUANTITIES'
MESSAGE 'RF Voltage [MV] =         '//$FORMAT(V(1),G8.6)
MESSAGE 'RF Frequency [GHz] =      '//$FORMAT(F(1),G8.6)
MESSAGE 'RF Period [ns] =          '//$FORMAT(1./F(1),G8.6)
MESSAGE 'Harmonic number =         '//$FORMAT(H(1),G8.6)
MESSAGE 'Synchronous phase [rad] = '//$FORMAT(PHIS(1),G8.6)
MESSAGE 'Synchrotron tune =        '//$FORMAT(NU(1),G8.6)
MESSAGE '1 / (synchrotron tune) =  '//$FORMAT(1./NU(1),G8.6)
MESSAGE 'Synchrotron freq. [Hz] =  '//$FORMAT(NU(1)*F(1)/H(1)*1E9,G8.6)
MESSAGE 'Bucket width [ns] =       '//$FORMAT(DT(1),G8.6)
MESSAGE 'Bucket area [MeV ns] =    '//$FORMAT(A(1),G8.6)
MESSAGE 'Eta =                     '//$FORMAT(ETA(1),G8.6)
MESSAGE 'Gamma_t =                 '//$FORMAT(GAMMAT(1),G8.6)

OPT=1
XMIN=0; XMAX=20
DEMIN=-600.
DEMAX=600.

DRAWFRAME:

MESSAGE
MESSAGE 'PLOTTING OPTIONS'

CHOOSE:
READ OPT 'X-axis = phase (1) or time (2)?'
IF (([OPT].NE.1).AND.([OPT].NE.2)) THEN
  MESSAGE ' *** Invalid option '//[OPT]
  GOTO CHOOSE
ENDIF

CASE [OPT] IN
(1) VA='[f]   '; TX='(rad)'
(2) VA='t   ';TX='(ns)'
ENDCASE

READ XMIN 'X-axis lower limit '//[TX]
READ XMAX 'X-axis upper limit '//[TX]

READ DEMIN 'Y-axis lower limit [MeV]'
READ DEMAX 'Y-axis upper limit [MeV]'

OPT GRID
NULL [XMIN] [XMAX] [DEMIN] [DEMAX]
ATIT [VA]//[TX] '[D]E   (MeV)'

COL=1
PHINEW=3.14
DTNEW=4.5
DENEW=100.
NP=144

ANOTHER:

MESSAGE
CASE [OPT] IN
(1) READ PHINEW 'Delta phi [rad]'; V/INP PHINEW(1) [PHINEW]
(2) READ DTNEW 'Delta t [ns]'; SIGMA PHINEW=2.*PI*F*[DTNEW]
ENDCASE
SIGMA PHIOLD=PHINEW

READ DENEW 'Delta E [MeV]'
V/INP DENEW(1) [DENEW]
SIGMA DEOLD=DENEW

READ NP 'Number of turns'

COL=[COL]+1
IF ([COL].GT.7) THEN
  COL=2
ENDIF
SET PMCI [COL]

MRKR=20
SIZE=0.03

DO I=1,[NP]
  CASE [OPT] IN
  (1) GRA/HPL/SYM PHINEW DENEW 1 [MRKR] [SIZE]
  (2) SIGMA DT=PHINEW/2./PI/F; GRA/HPL/SYM DT DENEW 1 [MRKR] [SIZE]
  ENDCASE
  V/COP PHINEW PHIOLD
  V/COP DENEW DEOLD
  SIGMA DENEW = DEOLD + EV * ( SIN(PHIOLD) - SIN(PHIS) )
  SIGMA PHINEW = PHIOLD + DENEW * 2.*PI*H*ETA/BETA/BETA/ES
ENDDO

ANS='N'
READ ANS 'Change plot options? [y/n]'
IF (([ANS]='y').OR.([ANS]='Y')) GOTO DRAWFRAME

ANS='Y'
READ ANS 'Another particle? [y/n]'
IF (([ANS]='y').OR.([ANS]='Y')) GOTO ANOTHER

RETURN



MACRO INTEGR

IF ($VEXIST(X1).NE.0) THEN
  V/DEL X1
ENDIF
IF ($VEXIST(X2).NE.0) THEN
  V/DEL X2
ENDIF
IF ($VEXIST(NN0).NE.0) THEN
  V/DEL NN0
ENDIF
IF ($VEXIST(NN).NE.0) THEN
  V/DEL NN
ENDIF
V/CREA X1(1) R
V/CREA X2(1) R
V/CREA NN0(1) R 1000.
V/CREA NN(1) R 1000.

SIGMA X2=PI/2.
SIGMA X1=X2-DPHI

SIGMA PHI0=ARRAY($SIGMA(NN0(1)),0.#$SIGMA(X2(1)))
SIGMA PHI=ARRAY($SIGMA(NN(1)),$SIGMA(X1(1))#$SIGMA(X2(1)))

SIGMA STEP0 = (X2-0.)/(NN0-1.)
SIGMA STEP  = (X2-X1)/(NN-1.)

SIGMA Y0=SQRT(COS(PHI0))
SIGMA Y=SQRT(ABS( COS(PHI) + PHI*SIN(PHIS(1)) - PI/2.*SIN(PHIS) ))

SIGMA I0=QUAD(Y0,STEP0)
SIGMA I=QUAD(Y,STEP)

RETURN