The loop below is a standard loop for lattice functions. It computes two types of things:
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)