Skip to content

Commit

Permalink
Fix for PTC 6D closed orbit seach with TOTALPATH=true; update of PTC;
Browse files Browse the repository at this point in the history
  • Loading branch information
@piotrskowronski
piotrskowronski committed Jan 31, 2018
1 parent 9187481 commit 4dbb579
Show file tree
Hide file tree
Showing 15 changed files with 933 additions and 439 deletions.
26 changes: 19 additions & 7 deletions libs/ptc/src/Ci_tpsa.f90
View file
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
View file
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
View file

Large diffs are not rendered by default.

165 changes: 87 additions & 78 deletions libs/ptc/src/Si_def_element.f90
View file
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

0 comments on commit 4dbb579

Please sign in to comment.