!=================================================================
subroutine Track_a_ray_with_spin
use madx_ptc_module
use pointer_lattice
implicit none
type(layout), pointer :: ALS
REAL(dp) x0(6),X(6),sig0(6),BETA_XX,BETA_XY,theta0,N0(3),DN0_DX(3),DN0_DX2(3),dx
type(probe_8) xs
type(probe) xs0,XST,xs0m,XSTm, xs02,XST2,xs0m2,XSTm2
type(real_8) ray8(6)
TYPE(INTERNAL_STATE) STATE
type(fibre), pointer :: p
type(damap) id,one_turn,A
integer i,mf,n,nh
type(normalform) nf
type(damapspin) id_s,one_turn_S,a_S,A_f,A_spin,A_l,A_nl,NORMALSPIN
type(normal_spin) nf_S
TYPE(spinor_8) N_AXIS
TYPE(spinor) N_AXIS_avep,N_AXIS_avem
als=>M_U%start
STATE=DEFAULT0+only_4d0
STATE=STATE+SPIN0
WRITE(6,*) "SPIN IS ON "
als=>m_U%start
sig0=0.00001d0
call MESS_UP_ALIGNMENT(als,SIG0,5.d0)
x0=0.d0
state=default0+nocavity0+DELTA0
CALL FIND_ORBIT(ALS,x0,1,STATE,c_1d_5)
write(6,*) "CLOSED ORBIT "
WRITE(6,'(6(1x,E15.8))') x0
CALL INIT(STATE,1,0)
call alloc(ray8)
call alloc(id,one_turn,A)
call alloc(nf)
id=1
ray8=x0+id
call track(als,ray8,1,state)
one_turn=ray8
nf=one_turn
write(6,*) " TUNES "
WRITE(6,'(2(1x,E15.8))') nf%tune(1:2)
RAY8=x0+NF%A_T
write(6,*) " Regular Twiss in file LATTICE.TXT "
CALL KANALNUMMER(MF,"LATTICE.TXT")
P=>ALS%START
A=RAY8
BETA_XX=(A%V(1).SUB.'1')**2+(A%V(1).SUB.'01')**2
BETA_XY=(A%V(1).SUB.'001')**2+(A%V(1).SUB.'0001')**2
WRITE(MF,*) P%MAG%NAME,BETA_XX,BETA_XY
DO I=1,ALS%N
call track(als,ray8,I,I+1,state)
A=RAY8
BETA_XX=(A%V(1).SUB.'1')**2+(A%V(1).SUB.'01')**2
BETA_XY=(A%V(1).SUB.'001')**2+(A%V(1).SUB.'0001')**2
P=>P%NEXT
WRITE(MF,'(a8,1x,2(1x,E15.8))') P%MAG%NAME(1:8),BETA_XX,BETA_XY
ENDDO
CLOSE(MF)
call kill(ray8)
call kill(id)
call kill(nf)
state=default0+DELTA0+SPIN0
XS0=x0
CALL FIND_ORBIT_probe_spin(ALS,XS0,STATE,c_1d_5,FIBRE1=1,theta0=theta0)
write(6,*) "CLOSED ORBIT "
WRITE(6,'(6(1x,E15.8))') XS0%X
write(6,*) " N0 "
WRITE(6,*) XS0%S(1)%X
WRITE(6,*) " THETA0"
WRITE(6,'(1x,E15.8)') THETA0
N0=XS0%S(1)%X
WRITE(6,*)" Setting the second spin axis to e_y : XS0%S(2)=2 "
XS0%S(2)=2
CALL TRACK_PROBE(ALS,XS0,STATE,FIBRE1=1)
WRITE(6,*)" tracking n0: should be returning to iself "
WRITE(6,'(3(1x,E15.8))') XS0%S(1)%X
WRITE(6,*)" tracking y_y: should not be invariant "
WRITE(6,'(3(1x,E15.8))') XS0%S(2)%X
WRITE(6,*)" "
WRITE(6,*)" Doing a Spin Normal Form for n-vector "
CALL INIT(STATE,2,0)
call alloc(ray8)
call alloc(id_S,one_turn_S,A_S,A_f,A_spin,A_l,A_nl,NORMALSPIN)
call alloc(nf_S)
ID_S=1
XS=XS0+ID_S
WRITE(6,*)" Computing a one-turn spin taylor map "
CALL TRACK_PROBE(ALS,XS,STATE,FIBRE1=1)
ID_S=XS
WRITE(6,*)" Doing a Spin-Orbital Normal Form: NF_S=ID_S "
NF_S=ID_S
A_S=NF_S%A_T
call factor(A_S,A_f,A_spin,A_l,A_nl)
WRITE(16,*) " FIXED POINT "
CALL PRINT(A_f,16)
WRITE(16,*) " SPIN PART "
CALL PRINT(A_spin,16)
WRITE(16,*) " ORBITAL LINEAR "
CALL PRINT(A_L,16)
WRITE(16,*) " ORBITAL NONLINEAR "
CALL PRINT(A_NL,16)
NORMALSPIN=A_S**(-1)*ID_S*A_S
WRITE(16,*) " NORMAL FORM "
CALL PRINT(NORMALSPIN,16)
N_AXIS=2
WRITE(16,*) " NORMALIZED N-AXIS "
CALL PRINT(N_AXIS,16)
WRITE(16,*) " N-AXIS "
N_AXIS= A_spin*N_AXIS
CALL PRINT(N_AXIS,16)
N0(1)=N_AXIS%X(1).SUB.'0'
N0(2)=N_AXIS%X(2).SUB.'0'
N0(3)=N_AXIS%X(3).SUB.'0'
WRITE(6,*) " n0 computed from N-AXIS "
WRITE(6,'(3(1x,E15.8))') N0
DN0_DX(1)=N_AXIS%X(1).SUB.'001'
DN0_DX(2)=N_AXIS%X(2).SUB.'001'
DN0_DX(3)=N_AXIS%X(3).SUB.'001'
DN0_DX2(1)=2*N_AXIS%X(1).SUB.'002'
DN0_DX2(2)=2*N_AXIS%X(2).SUB.'002'
DN0_DX2(3)=2*N_AXIS%X(3).SUB.'002'
write(6,*) " Spin-orbit Twiss in file LATTICE_SPIN.TXT "
CALL KANALNUMMER(MF,"LATTICE_SPIN.TXT")
P=>ALS%START
XS=XS0+NF_S%A_T ! Probe= closed orbit + A (A as in M=AoRoA**(-1) )
DO I=1,ALS%N
! send probe across each element EXACTLY
CALL TRACK_PROBE(ALS,XS,STATE,FIBRE1=I,FIBRE2=I+1)
!!!!!! USER DEFINED
A_S=XS ! turn the exact probe into a Berz-like obscenity
! Compute an obscenity here. Here it is the Ripken beta function
BETA_XX=(A_S%M%V(1).SUB.'1')**2+(A_S%M%V(1).SUB.'01')**2
BETA_XY=(A_S%M%V(1).SUB.'001')**2+(A_S%M%V(1).SUB.'0001')**2
!!!!!!! END OF USER DEFINED
!!!!!! EXTRA USER DEFINED
call factor(A_S,A_f,A_spin,A_l,A_nl) ! A_S= A_f o A_spin o A_l o A_nl
n_axis=2; ! creates e_y=(0,1,0) FPP arbitrary choice for the normal spin n-axis
n_axis=A_spin*n_axis ! The spin part of A_total acts on e_y. That gives you the n-field
!!!!!!! END OF EXTRA USER DEFINED
P=>P%NEXT
DN0_DX(1)=N_AXIS%X(1).SUB.'001'
DN0_DX(2)=N_AXIS%X(2).SUB.'001'
DN0_DX(3)=N_AXIS%X(3).SUB.'001'
DN0_DX2(1)=2*N_AXIS%X(1).SUB.'002'
DN0_DX2(2)=2*N_AXIS%X(2).SUB.'002'
DN0_DX2(3)=2*N_AXIS%X(3).SUB.'002'
WRITE(MF,'(a8,1x,8(1x,E15.8))') P%MAG%NAME(1:8),BETA_XX,BETA_XY,DN0_DX,DN0_DX2
ENDDO
CLOSE(MF)
WRITE(6,*) " !!! STROBOSCOPIC AVERAGE A LA HOFFSTADER !!!! "
!!! STROBOSCOPIC AVERAGE AS HOFFSTADER !!!!
n=10000
dx=1.d-4
nh=n*2 ! To avoid printing intermediate results
XS0=x0
XS0%X(3)=XS0%X(3)+dx
XS0%S(1)=1
XS0%S(2)=2
XS0%S(3)=3
XST=0
XS0m=x0
XS0m%X(3)=XS0m%X(3)-dx
XS0m%S(1)=1
XS0m%S(2)=2
XS0m%S(3)=3
XSTm=0
9 format('(a20,3(1x,E15.8))')
2001 call stroboscopic_average(als,xs0,xst,1,STATE,n,nh,N_AXIS_avep)
2002 call stroboscopic_average(als,xs0m,xstm,1,STATE,n,nh,N_AXIS_avem)
write(6,*)" $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ "
WRITE(6,9) " Numerical dn/dy =",(N_AXIS_avep%x(:)-N_AXIS_avem%x(:))/dx/2.d0
WRITE(6,9) " Analytical dn/dy =", DN0_DX
WRITE(6,9) " Numerical d2n/dy2 =",(N_AXIS_avep%x(:)+N_AXIS_avem%x(:)-2*n0(:))/dx**2
WRITE(6,9) "Analytical d2n/dy2 =", DN0_DX2
write(6,*) " more -> 1 for yes"
read(5,*)i
if(i==1) goto 2001
end subroutine Track_a_ray_with_spin