Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix for PTC 6D closed orbit seach with TOTALPATH=true; update of PTC; #551

Merged
merged 1 commit into from Jan 31, 2018
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
26 changes: 19 additions & 7 deletions libs/ptc/src/Ci_tpsa.f90
Expand Up @@ -9942,12 +9942,14 @@ end subroutine c_normal_radiation
subroutine c_stochastic_kick(m,ait,ki,eps)
!#general: stochastic kick
!# This routine creates a random kick
!# It diagonalises the beam envelope kick by noting
!# that this object is dual to Lie operators
implicit none
type(c_damap) , intent(inout) :: m
real(dp), intent(out) :: ait(6,6),ki(6)
real(dp), intent(in) :: eps
integer i,j
real(dp) norm,f(6,6),s(6,6), at(6,6),ai(6,6)
real(dp) norm,normy,f(6,6),s(6,6), at(6,6),ai(6,6)
real(dp) b(6,6),a(6,6),d(6,6)
type(c_vector_field) vf
type(c_damap) id
Expand Down Expand Up @@ -9980,18 +9982,24 @@ subroutine c_stochastic_kick(m,ait,ki,eps)
enddo
enddo

norm=0

!!!! In rings without errors and middle plane symmetry
!!! the y-part of the envelope is exactly zero
!!! So I add a y-part to make sure that my normal form
!!! does not crap out
normy=0
do i=1,6
do j=1,6
yhere=i==3.or.i==4.or.j==3.or.j==4
if(yhere) norm=norm +abs(f(i,j))
if(yhere) normy=normy +abs(f(i,j))
enddo
enddo
! write(6,*) " norm y",norm
if(norm<eps) then
if((normy/10)/norm<eps) then
f(3,3)=0.234567_dp*twopi
f(4,4)=0.234567_dp*twopi
endif
!!!!!!!!!!!!!!!
f=matmul(f,s)
do i=1,6
do j=1,6
Expand Down Expand Up @@ -10022,16 +10030,20 @@ subroutine c_stochastic_kick(m,ait,ki,eps)
enddo

! construct z_i =ki(i) * r_i
! then x=

! then x= Ait z is the appropriate kick
! in the original space.
!

if(global_verbose) then
write(6,*)" Stochastic kick"
do i=1,6
do j=1,6
if(d(i,j)/=0.d0) then
if(abs(d(i,j))>eps*norm) then
write(6,*) i,j,d(i,j)
endif
enddo
enddo
endif

ai=-matmul(matmul(s,at),s)
ait=transpose(ai)
Expand Down
24 changes: 14 additions & 10 deletions libs/ptc/src/Sc_euclidean.f90
Expand Up @@ -547,23 +547,27 @@ SUBROUTINE ROT_XZR(A,X,b,EXACT,ctime)
real(dp),INTENT(INOUT):: X(6)
real(dp) XN(6),PZ,PT
real(dp),INTENT(IN):: A,b
real(dp) sina, cosa, tana
LOGICAL(lp),INTENT(IN):: EXACT,ctime

IF(EXACT) THEN
COSA = COS(A)
SINA = SIN(A)
TANA = TAN(A)
if(ctime) then
PZ=ROOT(1.0_dp+2.0_dp*x(5)/b+X(5)**2-X(2)**2-X(4)**2)
PT=1.0_dp-X(2)*TAN(A)/PZ
XN(1)=X(1)/COS(A)/PT
XN(2)=X(2)*COS(A)+SIN(A)*PZ
XN(3)=X(3)+X(4)*X(1)*TAN(A)/PZ/PT
XN(6)=X(6)+X(1)*TAN(A)/PZ/PT*(1.0_dp/b+x(5))
PT=1.0_dp-X(2)*TANA/PZ
XN(1)=X(1)/COSA/PT
XN(2)=X(2)*COSA+SINA*PZ
XN(3)=X(3)+X(4)*X(1)*TANA/PZ/PT
XN(6)=X(6)+X(1)*TANA/PZ/PT*(1.0_dp/b+x(5))
else
PZ=ROOT((1.0_dp+X(5))**2-X(2)**2-X(4)**2)
PT=1.0_dp-X(2)*TAN(A)/PZ
XN(1)=X(1)/COS(A)/PT
XN(2)=X(2)*COS(A)+SIN(A)*PZ
XN(3)=X(3)+X(4)*X(1)*TAN(A)/PZ/PT
XN(6)=X(6)+(1.0_dp+X(5))*X(1)*TAN(A)/PZ/PT
PT=1.0_dp-X(2)*TANA/PZ
XN(1)=X(1)/COSA/PT
XN(2)=X(2)*COSA+SINA*PZ
XN(3)=X(3)+X(4)*X(1)*TANA/PZ/PT
XN(6)=X(6)+(1.0_dp+X(5))*X(1)*TANA/PZ/PT
endif
X(1)=XN(1)
X(2)=XN(2)
Expand Down
100 changes: 50 additions & 50 deletions libs/ptc/src/Sh_def_kind.f90

Large diffs are not rendered by default.

165 changes: 87 additions & 78 deletions libs/ptc/src/Si_def_element.f90
Expand Up @@ -135,73 +135,7 @@ MODULE S_DEF_ELEMENT


CONTAINS

SUBROUTINE decode_element(EL)
IMPLICIT NONE
TYPE(ELEMENT),INTENT(INOUT):: EL

SELECT CASE(EL%KIND)
CASE(KIND0)
print*,"this is KIND0"
case(KIND1)
print*,"this is KIND1"
case(KIND2)
print*,"this is KIND2"
case(KIND3)
print*,"this is KIND3"
case(KIND4)
print*,"this is KIND4"
case(KIND5)
print*,"this is KIND5"
case(KIND6)
print*,"this is KIND6"
case(KIND7)
print*,"this is KIND7"
case(KIND8)
print*,"this is KIND8"
case(KIND9)
print*,"this is KIND9"
case(KIND10)
print*,"this is KIND11"
CASE(KIND11)
print*,"this is KIND12"
CASE(KIND12)
print*,"this is KIND13"
CASE(KIND13)
print*,"this is KIND14"
CASE(KIND14)
print*,"this is KIND11"
CASE(KIND15)
print*,"this is KIND15"
CASE(KIND16)
print*,"this is KIND16"
CASE(KIND18)
print*,"this is KIND18"
CASE(KIND19)
print*,"this is KIND19"
CASE(KIND20)
print*,"this is KIND20"
CASE(KIND21)
print*,"this is KIND21"
CASE(KIND22)
print*,"this is KIND22"
case(KINDWIGGLER)
print*,"this is KINDWIGGLER"
case(KINDPA)
print*,"this is KINDPA"
case(kindsuperdrift)
print*,"this is KINDSUPERDRIFT"
case(KINDABELL)
print*,"this is KINDABELL"

case default

write(6,'(1x,i4,a21)') el%kind," not supported decode_element"
! call !write_e(0)
END SELECT


end SUBROUTINE decode_element

SUBROUTINE TRACKR(EL,X,K,MID)
IMPLICIT NONE
Expand Down Expand Up @@ -2846,8 +2780,8 @@ SUBROUTINE ZERO_EL(EL,I)
! EL=DEFAULT;
! ANBN
CALL ZERO_ANBN(EL,I)
ALLOCATE(EL%FINT);EL%FINT=0.5_dp;
ALLOCATE(EL%HGAP);EL%HGAP=0.0_dp;
ALLOCATE(EL%FINT(2));EL%FINT(1)=0.5_dp;EL%FINT(2)=0.5_dp;
ALLOCATE(EL%HGAP(2));EL%HGAP(1)=0.0_dp;EL%HGAP(2)=0.0_dp;
ALLOCATE(EL%H1);EL%H1=0.0_dp;
ALLOCATE(EL%H2);EL%H2=0.0_dp;
ALLOCATE(EL%VA);EL%VA=0.0_dp;
Expand Down Expand Up @@ -3116,8 +3050,8 @@ SUBROUTINE ZERO_ELP(EL,I)
! EL=DEFAULT;
! ANBN
CALL ZERO_ANBN(EL,I)
ALLOCATE(EL%FINT);CALL ALLOC(EL%FINT);EL%FINT=0.5_dp;
ALLOCATE(EL%HGAP);CALL ALLOC(EL%HGAP);EL%HGAP=0.0_dp;
ALLOCATE(EL%FINT(2));CALL ALLOC(EL%FINT);EL%FINT(1)=0.5_dp;EL%FINT(2)=0.5_dp;
ALLOCATE(EL%HGAP(2));CALL ALLOC(EL%HGAP);EL%HGAP(1)=0.0_dp;EL%HGAP(2)=0.0_dp;
ALLOCATE(EL%H1);CALL ALLOC(EL%H1);EL%H1=0.0_dp;
ALLOCATE(EL%H2);CALL ALLOC(EL%H2);EL%H2=0.0_dp;
ALLOCATE(EL%VA);CALL ALLOC(EL%VA);EL%VA=0.0_dp;
Expand Down Expand Up @@ -3165,8 +3099,10 @@ SUBROUTINE copy_el_elp(ELP,EL)
ELP%vorname=EL%vorname
ELP%KIND=EL%KIND
ELP%L=EL%L
ELP%FINT=EL%FINT
ELP%HGAP=EL%HGAP
ELP%FINT(1)=EL%FINT(1)
ELP%FINT(2)=EL%FINT(2)
ELP%HGAP(1)=EL%HGAP(1)
ELP%HGAP(2)=EL%HGAP(2)
ELP%H1=EL%H1
ELP%H2=EL%H2
ELP%VA=EL%VA
Expand Down Expand Up @@ -3548,8 +3484,10 @@ SUBROUTINE copy_elp_el(ELP,EL)
ELP%vorname=EL%vorname
ELP%KIND=EL%KIND
ELP%L=EL%L
ELP%FINT=EL%FINT
ELP%HGAP=EL%HGAP
ELP%FINT(1)=EL%FINT(1)
ELP%FINT(2)=EL%FINT(2)
ELP%HGAP(1)=EL%HGAP(1)
ELP%HGAP(2)=EL%HGAP(2)
ELP%H1=EL%H1
ELP%H2=EL%H2
ELP%VA=EL%VA
Expand Down Expand Up @@ -3930,8 +3868,10 @@ SUBROUTINE copy_el_el(ELP,EL)
ELP%KIND=EL%KIND
ELP%PLOT=EL%PLOT
ELP%L=EL%L
ELP%FINT=EL%FINT
ELP%HGAP=EL%HGAP
ELP%FINT(1)=EL%FINT(1)
ELP%FINT(2)=EL%FINT(2)
ELP%HGAP(1)=EL%HGAP(1)
ELP%HGAP(2)=EL%HGAP(2)
ELP%H1=EL%H1
ELP%H2=EL%H2
ELP%VA=EL%VA
Expand Down Expand Up @@ -4307,8 +4247,10 @@ SUBROUTINE reset31(ELP)
ELP%knob=.FALSE.

CALL resetpoly_R31(ELP%L) ! SHARED BY EVERYONE
CALL resetpoly_R31(ELP%FINT) ! SHARED BY EVERYONE
CALL resetpoly_R31(ELP%HGAP) ! SHARED BY EVERYONE
CALL resetpoly_R31(ELP%FINT(1)) ! SHARED BY EVERYONE
CALL resetpoly_R31(ELP%FINT(2)) ! SHARED BY EVERYONE
CALL resetpoly_R31(ELP%HGAP(1)) ! SHARED BY EVERYONE
CALL resetpoly_R31(ELP%HGAP(2)) ! SHARED BY EVERYONE
CALL resetpoly_R31(ELP%H1) ! SHARED BY EVERYONE
CALL resetpoly_R31(ELP%H2) ! SHARED BY EVERYONE
CALL resetpoly_R31(ELP%VA) ! SHARED BY EVERYONE
Expand Down Expand Up @@ -4560,4 +4502,71 @@ subroutine remove_aperture_elp(el)
end subroutine remove_aperture_elp


SUBROUTINE decode_element(EL)
IMPLICIT NONE
TYPE(ELEMENT),INTENT(INOUT):: EL

SELECT CASE(EL%KIND)
CASE(KIND0)
print*,"this is KIND0"
case(KIND1)
print*,"this is KIND1"
case(KIND2)
print*,"this is KIND2"
case(KIND3)
print*,"this is KIND3"
case(KIND4)
print*,"this is KIND4"
case(KIND5)
print*,"this is KIND5"
case(KIND6)
print*,"this is KIND6"
case(KIND7)
print*,"this is KIND7"
case(KIND8)
print*,"this is KIND8"
case(KIND9)
print*,"this is KIND9"
case(KIND10)
print*,"this is KIND11"
CASE(KIND11)
print*,"this is KIND12"
CASE(KIND12)
print*,"this is KIND13"
CASE(KIND13)
print*,"this is KIND14"
CASE(KIND14)
print*,"this is KIND11"
CASE(KIND15)
print*,"this is KIND15"
CASE(KIND16)
print*,"this is KIND16"
CASE(KIND18)
print*,"this is KIND18"
CASE(KIND19)
print*,"this is KIND19"
CASE(KIND20)
print*,"this is KIND20"
CASE(KIND21)
print*,"this is KIND21"
CASE(KIND22)
print*,"this is KIND22"
case(KINDWIGGLER)
print*,"this is KINDWIGGLER"
case(KINDPA)
print*,"this is KINDPA"
case(kindsuperdrift)
print*,"this is KINDSUPERDRIFT"
case(KINDABELL)
print*,"this is KINDABELL"

case default

write(6,'(1x,i4,a21)') el%kind," not supported decode_element"
! call !write_e(0)
END SELECT


end SUBROUTINE decode_element

END MODULE S_DEF_ELEMENT