From 4dbb57947f2f0970f9ff97b8e0301bdf450f772b Mon Sep 17 00:00:00 2001 From: Piotr Date: Wed, 31 Jan 2018 12:26:30 +0100 Subject: [PATCH] Fix for PTC 6D closed orbit seach with TOTALPATH=true; update of PTC; --- libs/ptc/src/Ci_tpsa.f90 | 26 +- libs/ptc/src/Sc_euclidean.f90 | 24 +- libs/ptc/src/Sh_def_kind.f90 | 100 +-- libs/ptc/src/Si_def_element.f90 | 165 ++--- libs/ptc/src/Sl_family.f90 | 44 +- libs/ptc/src/Sn_mad_like.f90 | 2 +- libs/ptc/src/Sp_keywords.f90 | 28 +- libs/ptc/src/Sra_fitting.f90 | 601 ++++++++++++++++-- libs/ptc/src/a_def_all_kind.inc | 24 +- libs/ptc/src/a_def_element_fibre_layout.inc | 7 +- src/madx_ptc_module.f90 | 181 +++++- src/madx_ptc_twiss.f90 | 3 +- tests/share/ALS/als.seqx | 4 +- .../test-ptc-twiss-6D-ALS.madx | 3 + .../test-ptc-twiss-6D-ALS.ref | 160 +---- 15 files changed, 933 insertions(+), 439 deletions(-) diff --git a/libs/ptc/src/Ci_tpsa.f90 b/libs/ptc/src/Ci_tpsa.f90 index 04b356a10..9f278c707 100644 --- a/libs/ptc/src/Ci_tpsa.f90 +++ b/libs/ptc/src/Ci_tpsa.f90 @@ -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 @@ -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(normeps*norm) then write(6,*) i,j,d(i,j) endif enddo enddo +endif ai=-matmul(matmul(s,at),s) ait=transpose(ai) diff --git a/libs/ptc/src/Sc_euclidean.f90 b/libs/ptc/src/Sc_euclidean.f90 index 992ca8c10..ec9f908f1 100644 --- a/libs/ptc/src/Sc_euclidean.f90 +++ b/libs/ptc/src/Sc_euclidean.f90 @@ -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) diff --git a/libs/ptc/src/Sh_def_kind.f90 b/libs/ptc/src/Sh_def_kind.f90 index b5644f52e..493171dcd 100644 --- a/libs/ptc/src/Sh_def_kind.f90 +++ b/libs/ptc/src/Sh_def_kind.f90 @@ -4911,7 +4911,7 @@ SUBROUTINE EDGER(EL,BN,H1,H2,FINT,HGAP,I,X,k) IMPLICIT NONE logical(lp) :: doneitt=.true. real(dp),INTENT(INOUT):: X(6) - real(dp),INTENT(IN)::FINT,HGAP,H1,H2 + real(dp),INTENT(IN)::FINT(2),HGAP(2),H1,H2 real(dp),INTENT(IN),dimension(:)::BN TYPE(MAGNET_CHART),INTENT(IN):: EL INTEGER, INTENT(IN) :: I @@ -4924,7 +4924,7 @@ SUBROUTINE EDGER(EL,BN,H1,H2,FINT,HGAP,I,X,k) call ROT_XZ(EL%EDGE(1),X,EL%BETA0,DONEITT,k%TIME) CALL FACE(EL,BN,H1,X,k) endif - CALL FRINGE_dipole(EL,BN,FINT,HGAP,I,X,k) + CALL FRINGE_dipole(EL,BN,FINT(I),HGAP(I),I,X,k) IF(I==2) then !doubling exit angle if second half CALL FACE(EL,BN,H2,X,k) x(1)=x(1)+EL%LC*SIN((EL%EDGE(2)-EL%EDGE(1))*0.5_dp) @@ -4936,7 +4936,7 @@ SUBROUTINE EDGER(EL,BN,H1,H2,FINT,HGAP,I,X,k) x(1)=x(1)+EL%DIR*EL%LC*SIN((EL%EDGE(2)-EL%EDGE(1))*0.5_dp) CALL FACE(EL,BN,H2,X,k) endif - CALL FRINGE_dipole(EL,BN,FINT,HGAP,I,X,k) + CALL FRINGE_dipole(EL,BN,FINT(I),HGAP(I),I,X,k) IF(I==1) then !doubling entrance angle if first half (1/2 magnet for designer) CALL FACE(EL,BN,H1,X,k) call ROT_XZ(EL%EDGE(1),X,EL%BETA0,DONEITT,k%TIME) @@ -4959,14 +4959,14 @@ SUBROUTINE EDGER(EL,BN,H1,H2,FINT,HGAP,I,X,k) IF(EL%BEND_FRINGE.and.(.NOT.((I==1.AND.EL%KILL_ENT_FRINGE).OR.(I==2.AND.EL%KILL_EXI_FRINGE)))) THEN - X(4)=X(4)-TAN(EL%EDGE(I)-EL%DIR*EL%CHARGE*2.0_dp*FINT*HGAP*(1.0_dp+SIN(EL%EDGE(I))**2)*BN(1)/COS(EL%EDGE(I))) & + X(4)=X(4)-TAN(EL%EDGE(I)-EL%DIR*EL%CHARGE*2.0_dp*FINT(I)*HGAP(I)*(1.0_dp+SIN(EL%EDGE(I))**2)*BN(1)/COS(EL%EDGE(I))) & & *EL%DIR*EL%CHARGE*BN(1)*X(3) ! SECTOR WEDGE (PROT) + FRINGE !!! coming to you with sadistic delight fsad=0.0_dp - if(fint*hgap/=0.0_dp) then - fsad=1.d0/(fint*hgap*2)/36.0_dp + if(fint(i)*hgap(i)/=0.0_dp) then + fsad=1.d0/(fint(i)*hgap(I)*2)/36.0_dp endif c3=bn(1)**2*fsad x(4)=x(4)-4*c3*x(3)**3 @@ -4976,11 +4976,11 @@ SUBROUTINE EDGER(EL,BN,H1,H2,FINT,HGAP,I,X,k) else IF(EL%BEND_FRINGE.and.(.NOT.((I==1.AND.EL%KILL_ENT_FRINGE).OR.(I==2.AND.EL%KILL_EXI_FRINGE)))) THEN - CALL FRINGE_dipole(EL,BN,FINT,HGAP,I,X,k) + CALL FRINGE_dipole(EL,BN,FINT(I),HGAP(I),I,X,k) fsad=0.0_dp - if(fint*hgap/=0.0_dp) then - fsad=1.d0/(fint*hgap*2)/36.0_dp + if(fint(I)*hgap(I)/=0.0_dp) then + fsad=1.d0/(fint(I)*hgap(I)*2)/36.0_dp endif c3=bn(1)**2*fsad x(4)=x(4)-4*c3*x(3)**3 @@ -5003,7 +5003,7 @@ END SUBROUTINE EDGER SUBROUTINE EDGEP(EL,BN,H1,H2,FINT,HGAP,I,X,k) IMPLICIT NONE TYPE(REAL_8),INTENT(INOUT):: X(6) - TYPE(REAL_8),INTENT(IN)::FINT,HGAP,H1,H2 + TYPE(REAL_8),INTENT(IN)::FINT(2),HGAP(2),H1,H2 TYPE(REAL_8),INTENT(IN),dimension(:)::BN TYPE(MAGNET_CHART),INTENT(IN):: EL INTEGER, INTENT(IN) :: I @@ -5016,7 +5016,7 @@ SUBROUTINE EDGEP(EL,BN,H1,H2,FINT,HGAP,I,X,k) ! call ROT_XZ(EL%EDGE(1),X,EL%BETA0,DONEITT,k%TIME) !DONE IN TRUE_PARALLEL ROUTINE!!!!! CALL FACE(EL,BN,H1,X,k) endif - CALL FRINGE_dipole(EL,BN,FINT,HGAP,I,X,k) + CALL FRINGE_dipole(EL,BN,FINT(I),HGAP(I),I,X,k) IF(I==2) then !doubling exit angle if second half CALL FACE(EL,BN,H2,X,k) ! x(1)=x(1)+EL%LC*SIN((EL%EDGE(2)-EL%EDGE(1))*half) @@ -5028,7 +5028,7 @@ SUBROUTINE EDGEP(EL,BN,H1,H2,FINT,HGAP,I,X,k) ! x(1)=x(1)+EL%DIR*EL%LC*SIN((EL%EDGE(2)-EL%EDGE(1))*half) CALL FACE(EL,BN,H2,X,k) endif - CALL FRINGE_dipole(EL,BN,FINT,HGAP,I,X,k) + CALL FRINGE_dipole(EL,BN,FINT(I),HGAP(I),I,X,k) IF(I==1) then !doubling entrance angle if first half (1/2 magnet for designer) CALL FACE(EL,BN,H1,X,k) ! call ROT_XZ(EL%EDGE(1),X,EL%BETA0,DONEITT,k%TIME) @@ -5050,14 +5050,14 @@ SUBROUTINE EDGEP(EL,BN,H1,H2,FINT,HGAP,I,X,k) IF(EL%BEND_FRINGE.and.(.NOT.((I==1.AND.EL%KILL_ENT_FRINGE).OR.(I==2.AND.EL%KILL_EXI_FRINGE)))) THEN - X(4)=X(4)-TAN(EL%EDGE(I)-EL%DIR*EL%CHARGE*2.0_dp*FINT*HGAP*(1.0_dp+SIN(EL%EDGE(I))**2)*BN(1)/COS(EL%EDGE(I))) & + X(4)=X(4)-TAN(EL%EDGE(I)-EL%DIR*EL%CHARGE*2.0_dp*FINT(I)*HGAP(I)*(1.0_dp+SIN(EL%EDGE(I))**2)*BN(1)/COS(EL%EDGE(I))) & & *EL%DIR*EL%CHARGE*BN(1)*X(3) ! SECTOR WEDGE (PROT) + FRINGE !!! coming to you with sadistic delight call alloc(fsad,c3) fsad=0.0_dp - if(fint*hgap/=0.0_dp) then - fsad=1.d0/(fint*hgap*2)/36.0_dp + if(fint(I)*hgap(I)/=0.0_dp) then + fsad=1.d0/(fint(I)*hgap(I)*2)/36.0_dp endif c3=bn(1)**2*fsad x(4)=x(4)-4*c3*x(3)**3 @@ -5067,11 +5067,11 @@ SUBROUTINE EDGEP(EL,BN,H1,H2,FINT,HGAP,I,X,k) else IF(EL%BEND_FRINGE.and.(.NOT.((I==1.AND.EL%KILL_ENT_FRINGE).OR.(I==2.AND.EL%KILL_EXI_FRINGE)))) THEN call alloc(fsad,c3) - CALL FRINGE_dipole(EL,BN,FINT,HGAP,I,X,k) + CALL FRINGE_dipole(EL,BN,FINT(I),HGAP(I),I,X,k) fsad=0.0_dp - if(fint*hgap/=0.0_dp) then - fsad=1.d0/(fint*hgap*2)/36.0_dp + if(fint(I)*hgap(I)/=0.0_dp) then + fsad=1.d0/(fint(I)*hgap(I)*2)/36.0_dp endif c3=bn(1)**2*fsad x(4)=x(4)-4*c3*x(3)**3 @@ -11565,7 +11565,7 @@ SUBROUTINE fringe_TEAPOTr(EL,X,k,J) IF(EL%P%EDGE(1)/=0.0_dp) THEN CALL ROT_XZ(EL%P%EDGE(1),X,EL%P%BETA0,DONEITT,k%TIME) CALL FACE(EL%P,EL%BN,EL%H1,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) IF(k%FRINGE.or.EL%P%permfringe/=0) then IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) @@ -11579,7 +11579,7 @@ SUBROUTINE fringe_TEAPOTr(EL,X,k,J) CALL WEDGE(-EL%P%EDGE(1),X,k,EL2=EL) ELSE CALL FACE(EL%P,EL%BN,EL%H1,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) ENDIF @@ -11597,13 +11597,13 @@ SUBROUTINE fringe_TEAPOTr(EL,X,k,J) x(2)=x(2)+EL%P%EDGE(2)*el%bn(2)*(x(1)**2-x(3)**2) x(4)=x(4)-EL%P%EDGE(2)*el%bn(2)*(2.0_dp*x(1)*x(3)) endif - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) CALL FACE(EL%P,EL%BN,EL%H2,X,k) CALL ROT_XZ(EL%P%EDGE(2),X,EL%P%BETA0,DONEITT,k%TIME) ELSE IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) CALL FACE(EL%P,EL%BN,EL%H2,X,k) ENDIF ENDIF ! J=2 @@ -11615,7 +11615,7 @@ SUBROUTINE fringe_TEAPOTr(EL,X,k,J) IF(EL%P%EDGE(2)/=0.0_dp) THEN CALL ROT_XZ(EL%P%EDGE(2),X,EL%P%BETA0,DONEITT,k%TIME) CALL FACE(EL%P,EL%BN,EL%H2,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) IF(k%FRINGE.or.EL%P%permfringe/=0) then IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) @@ -11628,7 +11628,7 @@ SUBROUTINE fringe_TEAPOTr(EL,X,k,J) CALL WEDGE(-EL%P%EDGE(2),X,k,EL2=EL) ELSE CALL FACE(EL%P,EL%BN,EL%H2,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3) CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) ENDIF @@ -11646,13 +11646,13 @@ SUBROUTINE fringe_TEAPOTr(EL,X,k,J) x(2)=x(2)-EL%P%EDGE(1)*el%bn(2)*(x(1)**2-x(3)**2) x(4)=x(4)+EL%P%EDGE(1)*el%bn(2)*(2.0_dp*x(1)*x(3)) endif - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) CALL FACE(EL%P,EL%BN,EL%H1,X,k) CALL ROT_XZ(EL%P%EDGE(1),X,EL%P%BETA0,DONEITT,k%TIME) ELSE IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) CALL FACE(EL%P,EL%BN,EL%H1,X,k) ENDIF @@ -11678,7 +11678,7 @@ SUBROUTINE fringe_TEAPOTP(EL,X,k,J) IF(EL%P%EDGE(1)/=0.0_dp) THEN CALL ROT_XZ(EL%P%EDGE(1),X,EL%P%BETA0,DONEITT,k%TIME) CALL FACE(EL%P,EL%BN,EL%H1,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) IF(k%FRINGE.or.EL%P%permfringe/=0) then IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) @@ -11692,7 +11692,7 @@ SUBROUTINE fringe_TEAPOTP(EL,X,k,J) CALL WEDGE(-EL%P%EDGE(1),X,k,EL2=EL) ELSE CALL FACE(EL%P,EL%BN,EL%H1,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) ENDIF @@ -11710,13 +11710,13 @@ SUBROUTINE fringe_TEAPOTP(EL,X,k,J) x(2)=x(2)+EL%P%EDGE(2)*el%bn(2)*(x(1)**2-x(3)**2) x(4)=x(4)-EL%P%EDGE(2)*el%bn(2)*(2.0_dp*x(1)*x(3)) endif - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) CALL FACE(EL%P,EL%BN,EL%H2,X,k) CALL ROT_XZ(EL%P%EDGE(2),X,EL%P%BETA0,DONEITT,k%TIME) ELSE IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) CALL FACE(EL%P,EL%BN,EL%H2,X,k) ENDIF ENDIF ! J=2 @@ -11728,7 +11728,7 @@ SUBROUTINE fringe_TEAPOTP(EL,X,k,J) IF(EL%P%EDGE(2)/=0.0_dp) THEN CALL ROT_XZ(EL%P%EDGE(2),X,EL%P%BETA0,DONEITT,k%TIME) CALL FACE(EL%P,EL%BN,EL%H2,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) IF(k%FRINGE.or.EL%P%permfringe/=0) then IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) @@ -11741,7 +11741,7 @@ SUBROUTINE fringe_TEAPOTP(EL,X,k,J) CALL WEDGE(-EL%P%EDGE(2),X,k,EL2=EL) ELSE CALL FACE(EL%P,EL%BN,EL%H2,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3) CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) ENDIF @@ -11759,13 +11759,13 @@ SUBROUTINE fringe_TEAPOTP(EL,X,k,J) x(2)=x(2)-EL%P%EDGE(1)*el%bn(2)*(x(1)**2-x(3)**2) x(4)=x(4)+EL%P%EDGE(1)*el%bn(2)*(2.0_dp*x(1)*x(3)) endif - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) CALL FACE(EL%P,EL%BN,EL%H1,X,k) CALL ROT_XZ(EL%P%EDGE(1),X,EL%P%BETA0,DONEITT,k%TIME) ELSE IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) CALL FACE(EL%P,EL%BN,EL%H1,X,k) ENDIF @@ -12956,13 +12956,13 @@ SUBROUTINE fringe_STREXr(EL,X,k,J) ANGH=EL%P%B0*EL%P%LD*0.5_dp-EL%P%EDGE(1) CALL ROT_XZ(EL%P%EDGE(1),X,EL%P%BETA0,DONEITT,k%TIME) CALL FACE(EL%P,EL%BN,EL%H1,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) CALL WEDGE(ANGH,X,k,EL1=EL) ELSE - CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT,EL%HGAP,1,X,k) + CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT(1),EL%HGAP(1),1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3) CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) ENDIF @@ -12975,13 +12975,13 @@ SUBROUTINE fringe_STREXr(EL,X,k,J) CALL WEDGE(ANGH,X,k,EL1=EL) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) CALL FACE(EL%P,EL%BN,EL%H2,X,k) CALL ROT_XZ(EL%P%EDGE(2),X,EL%P%BETA0,DONEITT,k%TIME) ELSE IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) - CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT,EL%HGAP,2,X,k) + CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT(2),EL%HGAP(2),2,X,k) ENDIF ENDIF ! J @@ -12995,12 +12995,12 @@ SUBROUTINE fringe_STREXr(EL,X,k,J) ANGH=EL%P%B0*EL%P%LD*0.5_dp-EL%P%EDGE(2) CALL ROT_XZ(EL%P%EDGE(2),X,EL%P%BETA0,DONEITT,k%TIME) CALL FACE(EL%P,EL%BN,EL%H2,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) CALL WEDGE(ANGH,X,k,EL1=EL) ELSE - CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT,EL%HGAP,2,X,k) + CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT(2),EL%HGAP(2),2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) ENDIF @@ -13013,13 +13013,13 @@ SUBROUTINE fringe_STREXr(EL,X,k,J) CALL WEDGE(ANGH,X,k,EL1=EL) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) CALL FACE(EL%P,EL%BN,EL%H1,X,k) CALL ROT_XZ(EL%P%EDGE(1),X,EL%P%BETA0,DONEITT,k%TIME) ELSE IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) - CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT,EL%HGAP,1,X,k) + CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT(1),EL%HGAP(1),1,X,k) ENDIF ENDIF ! J @@ -13047,13 +13047,13 @@ SUBROUTINE fringe_STREXP(EL,X,k,J) ANGH=EL%P%B0*EL%P%LD*0.5_dp-EL%P%EDGE(1) CALL ROT_XZ(EL%P%EDGE(1),X,EL%P%BETA0,DONEITT,k%TIME) CALL FACE(EL%P,EL%BN,EL%H1,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) CALL WEDGE(ANGH,X,k,EL1=EL) ELSE - CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT,EL%HGAP,1,X,k) + CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT(1),EL%HGAP(1),1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3) CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) ENDIF @@ -13066,13 +13066,13 @@ SUBROUTINE fringe_STREXP(EL,X,k,J) CALL WEDGE(ANGH,X,k,EL1=EL) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) CALL FACE(EL%P,EL%BN,EL%H2,X,k) CALL ROT_XZ(EL%P%EDGE(2),X,EL%P%BETA0,DONEITT,k%TIME) ELSE IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) - CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT,EL%HGAP,2,X,k) + CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT(2),EL%HGAP(2),2,X,k) ENDIF ENDIF ! J @@ -13086,12 +13086,12 @@ SUBROUTINE fringe_STREXP(EL,X,k,J) ANGH=EL%P%B0*EL%P%LD*0.5_dp-EL%P%EDGE(2) CALL ROT_XZ(EL%P%EDGE(2),X,EL%P%BETA0,DONEITT,k%TIME) CALL FACE(EL%P,EL%BN,EL%H2,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,2,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(2),EL%HGAP(2),2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) CALL WEDGE(ANGH,X,k,EL1=EL) ELSE - CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT,EL%HGAP,2,X,k) + CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT(2),EL%HGAP(2),2,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,2,X,k) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,2,X,k) ENDIF @@ -13104,13 +13104,13 @@ SUBROUTINE fringe_STREXP(EL,X,k,J) CALL WEDGE(ANGH,X,k,EL1=EL) IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) - CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT,EL%HGAP,1,X,k) + CALL FRINGE_dipole(EL%P,EL%BN,EL%FINT(1),EL%HGAP(1),1,X,k) CALL FACE(EL%P,EL%BN,EL%H1,X,k) CALL ROT_XZ(EL%P%EDGE(1),X,EL%P%BETA0,DONEITT,k%TIME) ELSE IF(el%p%permfringe==2.or.el%p%permfringe==3) CALL FRINGE2QUAD(EL%P,EL%bn(2),EL%an(2),EL%VA,EL%VS,1,X,k) IF(k%FRINGE.or.el%p%permfringe==1.or.el%p%permfringe==3)CALL MULTIPOLE_FRINGE(EL%P,EL%AN,EL%BN,1,X,k) - CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT,EL%HGAP,1,X,k) + CALL EDGE_TRUE_PARALLEL(EL%P,EL%BN,EL%H1,EL%H2,EL%FINT(1),EL%HGAP(1),1,X,k) ENDIF ENDIF ! J diff --git a/libs/ptc/src/Si_def_element.f90 b/libs/ptc/src/Si_def_element.f90 index 482c8825a..cef7549ba 100644 --- a/libs/ptc/src/Si_def_element.f90 +++ b/libs/ptc/src/Si_def_element.f90 @@ -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 @@ -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; @@ -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; @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/libs/ptc/src/Sl_family.f90 b/libs/ptc/src/Sl_family.f90 index 9869d1349..e3650f3f0 100644 --- a/libs/ptc/src/Sl_family.f90 +++ b/libs/ptc/src/Sl_family.f90 @@ -951,8 +951,8 @@ recursive SUBROUTINE MISALIGN_FIBRE(S2,S1,OMEGA,BASIS,ADD,preserve_girder) ! S2%MAG%R(I)=S1(3+I); S2%MAGP%R(I)=S1(3+I); ! ENDDO DO I=1,3 - D(I)=S1(I); D(I)=S1(I); - R(I)=S1(3+I); R(I)=S1(3+I); + D(I)=S1(I); !D(I)=S1(I); + R(I)=S1(3+I); !R(I)=S1(3+I); ENDDO S2%CHART%D_IN=0.0_dp;S2%CHART%D_OUT=0.0_dp; S2%CHART%ANG_IN=0.0_dp;S2%CHART%ANG_OUT=0.0_dp; @@ -986,7 +986,10 @@ recursive SUBROUTINE MISALIGN_FIBRE(S2,S1,OMEGA,BASIS,ADD,preserve_girder) CALL ROTATE_FRAME(F,OMEGAT,ANGLE,1,BASIS=BASIST) IF(PRESENT(BASIS)) THEN ! MUST ROTATE THAT FRAME AS WELL FOR CONSISTENCY IN DEFINITION WHAT A MISALIGNMENT IS IN PTC - CALL GEO_ROT(BASIST,ANGLE,1) +! CALL GEO_ROT(BASIST,ANGLE,1) + e1=basist + CALL GEO_ROT(BASIST,e2,ANGLE,e1) ! 2018 correction: agrees with successive rotations followed by rotations + basist=e2 ELSE BASIST=F%MID ! ALREADY ROTATED ENDIF @@ -1048,43 +1051,62 @@ recursive SUBROUTINE MISALIGN_FIBRE(S2,S1,OMEGA,BASIS,ADD,preserve_girder) END SUBROUTINE MISALIGN_FIBRE + SUBROUTINE MAD_MISALIGN_FIBRE(S2,S1) ! MISALIGNS FULL FIBRE; FILLS IN CHART AND MAGNET_CHART IMPLICIT NONE REAL(DP),INTENT(IN):: S1(6) TYPE(FIBRE),target,INTENT(INOUT):: S2 REAL(DP) ENT(3,3),ENT1(3,3),ENT2(3,3),T(3),MAD_ANGLE(3),T_GLOBAL(3),ANGLE(3),MIS(6) + real(dp) r1(3,3),r2(3,3),r3(3,3) ent=S2%CHART%F%ent T(1)=S1(1);T(2)=S1(2);T(3)=S1(3); MAD_ANGLE(1)=-S1(4) MAD_ANGLE(2)=-S1(5) MAD_ANGLE(3)=S1(6) - CALL CHANGE_BASIS(T,ENT,T_GLOBAL,GLOBAL_FRAME) - ANGLE=0.0_dp; ANGLE(3)=MAD_ANGLE(3) + + + ANGLE=0.0_dp; ANGLE(1)=MAD_ANGLE(1) ent1=ent ent2=ent CALL GEO_ROT(ENT1,ENT,ANGLE,ENT2) - ANGLE=0.0_dp; ANGLE(1)=MAD_ANGLE(1) + + ANGLE=0.0_dp; ANGLE(2)=MAD_ANGLE(2) ent1=ent ent2=ent CALL GEO_ROT(ENT1,ENT,ANGLE,ENT2) - ANGLE=0.0_dp; ANGLE(2)=MAD_ANGLE(2) + + ANGLE=0.0_dp; ANGLE(3)=MAD_ANGLE(3) ent1=ent ent2=ent CALL GEO_ROT(ENT1,ENT,ANGLE,ENT2) + + + + - CALL CHANGE_BASIS(T_GLOBAL,GLOBAL_FRAME,T,ENT) CALL COMPUTE_ENTRANCE_ANGLE(S2%CHART%F%ent,ENT,ANGLE) - MIS(1:3)=T + + + ent1=S2%CHART%F%ent + ent2=ent + CALL CHANGE_BASIS(T,ENT1,T_GLOBAL,ent2) + + + + MIS(1:3)=T_GLOBAL MIS(4:6)=ANGLE ENT=S2%CHART%F%ent T=S2%CHART%F%A - call MISALIGN_SIAMESE(S2,MIS,T,ENT) - ! call MISALIGN_FIBRE(S2,MIS,S2%CHART%F%A,S2%CHART%F%ent) + + call MISALIGN_FIBRE(S2,MIS,T,ENT) END SUBROUTINE MAD_MISALIGN_FIBRE + + + ! NEW ROUTINES TO CHANGE LAYOUT using only magnets!!!! SUBROUTINE TRANSLATE_girder(S2,D,ORDER,BASIS,PATCH,PREC) ! TRANSLATES A fibre diff --git a/libs/ptc/src/Sn_mad_like.f90 b/libs/ptc/src/Sn_mad_like.f90 index 219d986a6..0a5744a27 100644 --- a/libs/ptc/src/Sn_mad_like.f90 +++ b/libs/ptc/src/Sn_mad_like.f90 @@ -63,7 +63,7 @@ module Mad_like real(dp) T1,T2,B0 real(dp) volt,freq0,harmon,lag,DELTA_E,BSOL real(dp) tilt - real(dp) FINT,hgap,h1,h2,X_COL,Y_COL + real(dp) FINT,hgap,FINT2,hgap2,h1,h2,X_COL,Y_COL real(dp) thin_h_foc,thin_v_foc,thin_h_angle,thin_v_angle,hf,vf,ls ! highly illegal additions by frs CHARACTER(120) file CHARACTER(120) file_rev diff --git a/libs/ptc/src/Sp_keywords.f90 b/libs/ptc/src/Sp_keywords.f90 index 6964b0aeb..917b711be 100644 --- a/libs/ptc/src/Sp_keywords.f90 +++ b/libs/ptc/src/Sp_keywords.f90 @@ -2915,12 +2915,14 @@ subroutine el_el0(f,dir,mf) ele0%VOLT_FREQ_PHAS=0.0_dp ele0%B_SOL=0.0_dp - ele0%fint_hgap_h1_h2_va_vs(1)=f%fint - ele0%fint_hgap_h1_h2_va_vs(2)=f%hgap - ele0%fint_hgap_h1_h2_va_vs(3)=f%h1 - ele0%fint_hgap_h1_h2_va_vs(4)=f%h2 - ele0%fint_hgap_h1_h2_va_vs(5)=f%va - ele0%fint_hgap_h1_h2_va_vs(6)=f%vs + ele0%fint_hgap_h1_h2_va_vs(1)=f%fint(1) + ele0%fint_hgap_h1_h2_va_vs(2)=f%fint(2) + ele0%fint_hgap_h1_h2_va_vs(3)=f%hgap(1) + ele0%fint_hgap_h1_h2_va_vs(4)=f%hgap(2) + ele0%fint_hgap_h1_h2_va_vs(5)=f%h1 + ele0%fint_hgap_h1_h2_va_vs(6)=f%h2 + ele0%fint_hgap_h1_h2_va_vs(7)=f%va + ele0%fint_hgap_h1_h2_va_vs(8)=f%vs ele0%L=f%L @@ -2994,12 +2996,14 @@ subroutine el_el0(f,dir,mf) f%bn(1:f%p%nmul)=ele0%bn(1:f%p%nmul) endif -f%fint= ele0%fint_hgap_h1_h2_va_vs(1) -f%hgap= ele0%fint_hgap_h1_h2_va_vs(2) -f%h1 = ele0%fint_hgap_h1_h2_va_vs(3) -f%h2 = ele0%fint_hgap_h1_h2_va_vs(4) -f%va = ele0%fint_hgap_h1_h2_va_vs(5) -f%vs = ele0%fint_hgap_h1_h2_va_vs(6) +f%fint(1)= ele0%fint_hgap_h1_h2_va_vs(1) +f%fint(2)=ele0%fint_hgap_h1_h2_va_vs(2) +f%hgap(1)= ele0%fint_hgap_h1_h2_va_vs(3) +f%hgap(2)= ele0%fint_hgap_h1_h2_va_vs(4) +f%h1 = ele0%fint_hgap_h1_h2_va_vs(5) +f%h2 = ele0%fint_hgap_h1_h2_va_vs(6) +f%va = ele0%fint_hgap_h1_h2_va_vs(7) +f%vs = ele0%fint_hgap_h1_h2_va_vs(8) if(f%kind==kind4.or.f%kind==kind21) then ! cavities diff --git a/libs/ptc/src/Sra_fitting.f90 b/libs/ptc/src/Sra_fitting.f90 index d71d393fe..7c033e43b 100644 --- a/libs/ptc/src/Sra_fitting.f90 +++ b/libs/ptc/src/Sra_fitting.f90 @@ -3071,10 +3071,11 @@ subroutine lattice_fit_bump_min_rcs(R0,R1,EPSF,poly,np,sca) end subroutine lattice_fit_bump_min_rcs - SUBROUTINE FIND_ORBIT_LAYOUT_da(RING,FIX,STATE,TURNS,fibre1,node1) ! Finds orbit without TPSA in State or compatible state + SUBROUTINE FIND_ORBIT_LAYOUT_da(RING,FIX,STATE,TURNS,fibre1,node1,total) ! Finds orbit without TPSA in State or compatible state IMPLICIT NONE TYPE(layout),target,INTENT(INOUT):: RING real(dp) , intent(inOUT) :: FIX(6) + real(dp), optional :: total INTEGER , optional,intent(in) :: TURNS,node1,fibre1 TYPE(INTERNAL_STATE),optional, intent(in) :: STATE type(fibre), pointer :: object_fibre1 @@ -3086,21 +3087,22 @@ SUBROUTINE FIND_ORBIT_LAYOUT_da(RING,FIX,STATE,TURNS,fibre1,node1) ! Finds orbit do i=1,fibre1-1 object_fibre1=>object_fibre1%next enddo - call FIND_ORBIT_LAYOUT_da_object(FIX,STATE,TURNS,fibre1=object_fibre1) + call FIND_ORBIT_LAYOUT_da_object(FIX,STATE,TURNS,fibre1=object_fibre1,total=total) else object_node1=>ring%t%start do i=1,node1-1 object_node1=>object_node1%next enddo - call FIND_ORBIT_LAYOUT_da_object(FIX,STATE,TURNS,node1=object_node1) + call FIND_ORBIT_LAYOUT_da_object(FIX,STATE,TURNS,node1=object_node1,total=total) endif end SUBROUTINE FIND_ORBIT_LAYOUT_da - SUBROUTINE FIND_ORBIT_LAYOUT_da_object(FIX0,STATE,TURNS,fibre1,node1) ! Finds orbit without TPSA in State or compatible state + SUBROUTINE FIND_ORBIT_LAYOUT_da_object(FIX0,STATE,TURNS,fibre1,node1,total) ! Finds orbit without TPSA in State or compatible state IMPLICIT NONE real(dp) , intent(inOUT) :: FIX0(6) + real(dp), optional :: total INTEGER , optional,intent(in) :: TURNS type(fibre), optional, pointer :: fibre1 type(integration_node), optional, pointer :: node1 @@ -3249,13 +3251,34 @@ SUBROUTINE FIND_ORBIT_LAYOUT_da_object(FIX0,STATE,TURNS,fibre1,node1) ! Finds or c=>c%next i=i+1 enddo - - if(freq_redefine) then + + if(ring%harmonic_number==0) then + t=>ring%t%start + tot=0 + + do i=1,ring%t%n + tot= tot + t%ds_ac + t=>t%next + enddo + + if(state%time) tot=tot/c%beta0 + + if(freq_redefine) then + ring%harmonic_number=tot*freq/twopi + else + ring%harmonic_number=tot*freq/CLIGHT + endif + + else + + if(freq_redefine) then tot=RING%HARMONIC_NUMBER*twopi/FREQ - else + else tot=RING%HARMONIC_NUMBER*CLIGHT/FREQ - endif + endif endif + + endif @@ -3390,14 +3413,16 @@ SUBROUTINE FIND_ORBIT_LAYOUT_da_object(FIX0,STATE,TURNS,fibre1,node1) ! Finds or ! FIX(6)=FIX(6)+freq*turns0 c_%APERTURE_FLAG=APERTURE fix0=fix + if(present(total)) total=tot END SUBROUTINE FIND_ORBIT_LAYOUT_da_object - SUBROUTINE FIND_ORBIT_LAYOUT_noda(RING,FIX,STATE,eps,TURNS,fibre1,node1) ! Finds orbit without TPSA in State or compatible state + SUBROUTINE FIND_ORBIT_LAYOUT_noda(RING,FIX,STATE,eps,TURNS,fibre1,node1,total) ! Finds orbit without TPSA in State or compatible state IMPLICIT NONE TYPE(layout),target,INTENT(INOUT):: RING real(dp) , intent(inOUT) :: FIX(6) + real(dp), optional :: total INTEGER , optional,intent(in) :: TURNS,node1,fibre1 real(dp) eps,TOT,freq,t6 TYPE(INTERNAL_STATE),optional, intent(in) :: STATE @@ -3410,22 +3435,23 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda(RING,FIX,STATE,eps,TURNS,fibre1,node1) ! Finds do i=1,fibre1-1 object_fibre1=>object_fibre1%next enddo - call FIND_ORBIT_LAYOUT_noda_object(FIX,STATE,eps,TURNS,fibre1=object_fibre1) + call FIND_ORBIT_LAYOUT_noda_object(FIX,STATE,eps,TURNS,fibre1=object_fibre1,total=total) else object_node1=>ring%t%start do i=1,node1-1 object_node1=>object_node1%next enddo - call FIND_ORBIT_LAYOUT_noda_object(FIX,STATE,eps,TURNS,node1=object_node1) + call FIND_ORBIT_LAYOUT_noda_object(FIX,STATE,eps,TURNS,node1=object_node1,total=total) endif end SUBROUTINE FIND_ORBIT_LAYOUT_noda - SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Finds orbit without TPSA in State or compatible state + SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1,total) ! Finds orbit without TPSA in State or compatible state IMPLICIT NONE TYPE(layout),pointer :: RING real(dp) , intent(inOUT) :: FIX0(6) + real(dp), optional :: total INTEGER , optional,intent(in) :: TURNS type(fibre), optional, pointer :: fibre1 type(integration_node), optional, pointer :: node1 @@ -3434,8 +3460,8 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi TYPE(INTERNAL_STATE) stat real(dp) DIX(6),xdix,xdix0,tiny,beta1,t6,freqmin - real(dp) X(6),Y(6),MX(6,6),sxi(6,6),SX(6,6),dt,dl,fix(6) - integer NO1,ND2,I,IU,ITE,ier,j,ITEM,try + real(dp) X(6),Y(6),MX(6,6),sxi(6,6),SX(6,6),fix(6), Y6start + integer NO1,ND2,I,IU,ITE,ier,j,ITEM !,try TYPE (fibre), POINTER :: C TYPE (integration_node), POINTER :: t logical(lp) APERTURE,use_bmad_units_temp @@ -3446,7 +3472,7 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi fix=fix0 didPhaseJump = .false. - + tot=0 if(present(fibre1)) then ring=>fibre1%parent_layout @@ -3487,6 +3513,7 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi endif call convert_bmad_to_ptc(fix,beta1,STATE%TIME) endif + use_bmad_units=.false. TURNS0=1 @@ -3554,9 +3581,6 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi if((stat%totalpath==1).and.(.not.stat%nocavity)) then - - if(global_verbose) print*,"Totpath and cavity: Looking for Frequency" - C=>RING%START freq=0.0_dp i=1 @@ -3572,15 +3596,45 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi c=>c%next i=i+1 enddo - - if(freq_redefine) then - tot=RING%HARMONIC_NUMBER*twopi/FREQ - if(global_verbose) print*,"Totpath and cavity: f_redefine TOT=",TOT - else - tot=RING%HARMONIC_NUMBER*CLIGHT/FREQ - if(global_verbose) print*,"Totpath and cavity: TOT=",TOT - endif - endif + + + if(global_verbose) write(6,*) " Using frequency ", freq + + + if(ring%harmonic_number==0) then + t=>ring%t%start + tot=0 + do i=1,ring%t%n + tot= tot + t%ds_ac + t=>t%next + enddo + + if(state%time) tot=tot/c%beta0 + + if(global_verbose) write(6,*) " Integrated tot of the machine ", tot + + if(freq_redefine) then + ring%harmonic_number=tot*freq/twopi + else + ring%harmonic_number=tot*freq/CLIGHT + endif + + if(global_verbose) write(6,*) " Corresponding harmonic number ", ring%harmonic_number + + else + + if(freq_redefine) then + tot=RING%HARMONIC_NUMBER*twopi/FREQ + else + tot=RING%HARMONIC_NUMBER*CLIGHT/FREQ + endif + + if(global_verbose) write(6,*) " Harmonic number already defined", ring%harmonic_number + if(global_verbose) write(6,*) " Corresponding tot ", tot + + endif + + endif @@ -3595,58 +3649,58 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi X=FIX DO I=1,TURNS0 - ! CALL TRACK(RING,X,LOC,STAT) - ! trackflag=TRACK_flag(RING,X,LOC,STAT) - !! xs%x=x - - if(global_verbose) print*,"ITEM=",ITEM," x=",x - call TRACK_probe_X(x,stat,fibre1=fibre1,node1=node1) if(.not.check_stable) then messagelost(len_trim(messagelost)+1:255)=" -> Unstable tracking guessed orbit " c_%APERTURE_FLAG=APERTURE - if(try>0) goto 1111 + ! if(try>0) goto 1111 return endif - ENDDO - ! write(6,*) x + + if(global_verbose) then + write(6,*) "#############################################" + write(6,*) "ITERATION ", ITEM + write(6,*) "" + write(6,*) "FIX start :", FIX(:) + write(6,*) "FIX end :", X(:) + write(6,*) "DIFF :", FIX - X + endif + mx=0.0_dp DO J=1,ND2 Y=FIX Y(J)=FIX(J)+EPS + Y6start=Y(6) DO I=1,TURNS0 - ! CALL TRACK(RING,Y,LOC,STAT) - !! xs%x=y call TRACK_probe_X(Y,stat,fibre1=fibre1,node1=node1) if(.not.check_stable) then messagelost(len_trim(messagelost)+1:255)=" -> Unstable while tracking small rays around the guessed orbit " ! fixed_found=my_false c_%APERTURE_FLAG=APERTURE - if(try>0) goto 1111 + ! if(try>0) goto 1111 return endif ENDDO - + if(stat%totalpath==1) then + y(6)=y(6)-TURNS0*gettot((y(6)/TURNS0) - Y6start , freq) + endif + do i=1,ND2 - if(i==6) then MX(I,J)=Y(i)/2/eps+MX(I,J) - ! MX(I,J)=(Y(i)-tot-X(i))/eps - else - ! MX(I,J)=(Y(i)-X(i))/eps - MX(I,J)=Y(i)/2/eps+MX(I,J) - endif enddo + Y=FIX Y(J)=FIX(J)-EPS + Y6start=Y(6) DO I=1,TURNS0 ! CALL TRACK(RING,Y,LOC,STAT) !! xs%x=y @@ -3656,34 +3710,47 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi messagelost(len_trim(messagelost)+1:255)=" -> Unstable while tracking small rays around the guessed orbit " ! fixed_found=my_false c_%APERTURE_FLAG=APERTURE - if(try>0) goto 1111 + ! if(try>0) goto 1111 return endif ENDDO + if(stat%totalpath==1) then + y(6)=y(6)-TURNS0*gettot(y(6)/TURNS0 - Y6start, freq) + endif do i=1,ND2 - if(i==6) then MX(I,J)=-Y(i)/2/eps+MX(I,J) - ! MX(I,J)=(Y(i)-tot-X(i))/eps - else - ! MX(I,J)=(Y(i)-X(i))/eps - MX(I,J)=-Y(i)/2/eps+MX(I,J) - endif enddo ENDDO + if(global_verbose) then + write(6,*) "" + write(6,*) "MX" + do i=1,ND2 + write(6,*) " ", MX(I,:) + enddo + endif + SX=MX; DO I=1,nd2 ! 6 before SX(I,I)=MX(I,I)-1.0_dp ENDDO + if(global_verbose) then + write(6,*) "" + write(6,*) "SX" + do i=1,ND2 + write(6,*) " ", SX(I,:) + enddo + endif + DO I=1,ND2 DIX(I)=FIX(I)-X(I) - if(i==6) dix(6)=dix(6)+tot + if(i==6 .and. stat%totalpath==1) dix(6)=dix(6)+gettot(X(6) - FIX(6),freq) enddo CALL matinv(SX,SXI,ND2,6,ier) @@ -3700,6 +3767,8 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi enddo enddo dix=x + + DO I=1,ND2 FIX(I)=FIX(I)+DIX(I) ENDDO @@ -3708,6 +3777,14 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi do iu=1,ND2 xdix=abs(dix(iu))+xdix enddo + + if(global_verbose) then + write(6,*) "" + write(6,*) "DIX (corr): ", DIX + write(6,*) "PENALTY F : ", XDIX + endif + + ! write(6,*) " Convergence Factor = ",nd2,xdix,deps_tracking ! pause 123321 ! if(verbose) write(6,*) " Convergence Factor = ",xdix @@ -3731,7 +3808,8 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi check_stable=my_false ! ENDIF ITE=0 - if(try>0) goto 1111 + ! if(try>0) goto 1111 + return endif ! write(6,*) item,xdix,xdix0,tiny @@ -3743,8 +3821,8 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi if (ND2 == 6.and.check_longitudinal) then isStableFixPoint = is_ORBIT_STABLE(FIX,EPS,STAT,fibre1,node1) - - if (isStableFixPoint .eqv. .false. ) then + + if (isStableFixPoint .eqv. .false.) then if (didPhaseJump ) then messagelost= "Found unstable fixed point" @@ -3754,15 +3832,385 @@ SUBROUTINE FIND_ORBIT_LAYOUT_noda_object(FIX0,STATE,eps,TURNS,fibre1,node1) ! Fi else if(global_verbose) print*,"Orbit seemed to be unstable in longitudinal" + + fix = fix0 + fix(6)= fix(6)+ clight/freqmin/2 + + didPhaseJump = .true. ! it is protection against infinite lopp + + goto 1111 + endif + endif + endif + + + if(use_bmad_units_temp) then + + call convert_ptc_to_bmad(fix,beta1,STATE%TIME) + endif + use_bmad_units=use_bmad_units_temp + ! FIX(6)=FIX(6)+freq*turns0 + c_%APERTURE_FLAG=APERTURE + fix0=fix + if(present(total)) total=tot + END SUBROUTINE FIND_ORBIT_LAYOUT_noda_object + + + + + SUBROUTINE FIND_ORBIT_LAYOUT_noda_object_orig(FIX0,STATE,eps,TURNS,fibre1,node1,total) ! Finds orbit without TPSA in State or compatible state + IMPLICIT NONE + TYPE(layout),pointer :: RING + real(dp) , intent(inOUT) :: FIX0(6) + real(dp), optional :: total + INTEGER , optional,intent(in) :: TURNS + type(fibre), optional, pointer :: fibre1 + type(integration_node), optional, pointer :: node1 + real(dp) eps,freq,tot + TYPE(INTERNAL_STATE),optional, intent(in) :: STATE + TYPE(INTERNAL_STATE) stat + + real(dp) DIX(6),xdix,xdix0,tiny,beta1,t6,freqmin + real(dp) X(6),Y(6),MX(6,6),sxi(6,6),SX(6,6),dt,dl,fix(6) + integer NO1,ND2,I,IU,ITE,ier,j,ITEM !,try + TYPE (fibre), POINTER :: C + TYPE (integration_node), POINTER :: t + logical(lp) APERTURE,use_bmad_units_temp + logical isStableFixPoint, didPhaseJump + INTEGER TURNS0,trackflag - fix = fix0 - fix(6)= fix(6)+ clight/freqmin/2 - didPhaseJump = .true. ! it is protection against + fix=fix0 + + didPhaseJump = .false. + + tot=0 + if(present(fibre1)) then + ring=>fibre1%parent_layout + else + ring=>node1%parent_fibre%parent_layout + endif + + + if (.not.STATE%NOCAVITY.and.check_longitudinal) then + freqmin=1.d38 + t6=0 + c=>ring%start + do i=1,ring%n + if(c%mag%kind==kind4.or.c%mag%kind==kind21) then + if(abs(c%mag%freq)>0) then + if( abs(c%mag%freq)<=freqmin) freqmin=c%mag%freq + endif + endif + t6=t6+c%mag%p%ld + c=>c%next + enddo + if(state%time) t6=t6/c%beta0 + if(global_verbose) write(6,*) " Harmonic # = ", freqmin*t6/clight + if(global_verbose) write(6,*) " freqmin , dt", freqmin,clight/freqmin/2 + endif - goto 1111 + ! fixed_found=my_true + !! type(probe) xs + if(.not.associated(RING%t)) call MAKE_NODE_LAYOUT(ring) + !! xs%x=zero + !! xs%s%x=zero + use_bmad_units_temp=use_bmad_units + if(use_bmad_units_temp) then + if(present(fibre1)) then + beta1=fibre1%mag%p%beta0 + else + beta1=node1%parent_fibre%mag%p%beta0 + endif + call convert_bmad_to_ptc(fix,beta1,STATE%TIME) + endif + + use_bmad_units=.false. + + TURNS0=1 + trackflag=0 + IF(PRESENT(TURNS)) TURNS0=TURNS + + APERTURE=c_%APERTURE_FLAG + c_%APERTURE_FLAG=.false. + + !! call move_to(ring,c,loc) + !! loct=c%t1%pos + + + Nullify(C); + + + dix(:)=0.0_dp + tiny=1e-40_dp + xdix0=1e4_dp*DEPS_tracking + NO1=1 + + if(.not.present(STATE)) then + IF(default%NOCAVITY) THEN + ! ND1=2 + stat=default+ only_4d + ELSE + ! ND1=3 + STAT=default + C=>RING%START + do i=1,RING%n + if(C%magp%kind==kind4.OR.C%magp%kind==kind21) goto 101 + C=>C%NEXT + enddo + + messagelost= " No Cavity in the Line " + check_stable=.false. + return + + ENDIF + else ! (.not.present(STATE)) t + IF(STATE%NOCAVITY) THEN + ND2=4 + STAT=STATE+only_4d0 + if(state%radiation) then + check_stable=.false. + + messagelost= " Cavity needed when radiation present " + return + endif + ELSE + ND2=6 + STAT=STATE + C=>RING%START + do i=1,RING%n + if(C%magp%kind==kind4.OR.C%magp%kind==kind21) goto 101 + C=>C%NEXT + enddo + check_stable=.false. + messagelost= " State present; no cavity: FIND_ORBIT_LAYOUT will crash => exiting" + return + + ENDIF + endif +101 continue + + + if((stat%totalpath==1).and.(.not.stat%nocavity)) then + C=>RING%START + freq=0.0_dp + i=1 + xdix=0.0_dp + do while(i<=RING%n) + if(associated(c%magp%freq)) then + IF(FREQ==0.0_dp) THEN + freq=c%magp%freq + ELSEIF(c%magp%freq/=0.0_dp.AND.c%magp%freqc%next + i=i+1 + enddo + + + if(global_verbose) write(6,*) " Using frequency ", freq + + + if(ring%harmonic_number==0) then + t=>ring%t%start + tot=0 + do i=1,ring%t%n + tot= tot + t%ds_ac + t=>t%next + enddo + + if(state%time) tot=tot/c%beta0 + + if(global_verbose) write(6,*) " Integrated tot of the machine ", tot + + if(freq_redefine) then + ring%harmonic_number=tot*freq/twopi + else + ring%harmonic_number=tot*freq/CLIGHT + endif + + if(global_verbose) write(6,*) " Corresponding harmonic number ", ring%harmonic_number + + else + + if(freq_redefine) then + tot=RING%HARMONIC_NUMBER*twopi/FREQ + else + tot=RING%HARMONIC_NUMBER*CLIGHT/FREQ endif + + if(global_verbose) write(6,*) " Harmonic number already defined", ring%harmonic_number + if(global_verbose) write(6,*) " Corresponding tot ", tot + + endif + + endif + + + + + +1111 continue + + + ITEM=0 +3 continue + ITEM=ITEM+1 + X=FIX + + DO I=1,TURNS0 + call TRACK_probe_X(x,stat,fibre1=fibre1,node1=node1) + + if(.not.check_stable) then + messagelost(len_trim(messagelost)+1:255)=" -> Unstable tracking guessed orbit " + c_%APERTURE_FLAG=APERTURE + ! if(try>0) goto 1111 + return + endif + + + ENDDO + + mx=0.0_dp + DO J=1,ND2 + Y=FIX + Y(J)=FIX(J)+EPS + DO I=1,TURNS0 + call TRACK_probe_X(Y,stat,fibre1=fibre1,node1=node1) + + if(.not.check_stable) then + messagelost(len_trim(messagelost)+1:255)=" -> Unstable while tracking small rays around the guessed orbit " + ! fixed_found=my_false + c_%APERTURE_FLAG=APERTURE + ! if(try>0) goto 1111 + return + endif + + + ENDDO + + if(stat%totalpath==1) then + y(6)=y(6)-TURNS0*tot + endif + do i=1,ND2 + MX(I,J)=Y(i)/2/eps+MX(I,J) + enddo + Y=FIX + Y(J)=FIX(J)-EPS + DO I=1,TURNS0 + ! CALL TRACK(RING,Y,LOC,STAT) + !! xs%x=y + call TRACK_probe_X(Y,stat,fibre1=fibre1,node1=node1) + + if(.not.check_stable) then + messagelost(len_trim(messagelost)+1:255)=" -> Unstable while tracking small rays around the guessed orbit " + ! fixed_found=my_false + c_%APERTURE_FLAG=APERTURE + ! if(try>0) goto 1111 + return + endif + + ENDDO + + if(stat%totalpath==1) then + y(6)=y(6)-TURNS0*tot + endif + + do i=1,ND2 + MX(I,J)=-Y(i)/2/eps+MX(I,J) + enddo + + + ENDDO + + SX=MX; + DO I=1,nd2 ! 6 before + SX(I,I)=MX(I,I)-1.0_dp + ENDDO + + DO I=1,ND2 + DIX(I)=FIX(I)-X(I) + if(i==6) dix(6)=dix(6)+tot + enddo + + CALL matinv(SX,SXI,ND2,6,ier) + IF(IER==132) then + messagelost= " Inversion failed in FIND_ORBIT_LAYOUT_noda" + check_stable=.false. + return + endif + + x=0.0_dp + do i=1,nd2 + do j=1,nd2 + x(i)=sxi(i,j)*dix(j)+x(i) + enddo + enddo + dix=x + DO I=1,ND2 + FIX(I)=FIX(I)+DIX(I) + ENDDO + + xdix=0.0_dp + do iu=1,ND2 + xdix=abs(dix(iu))+xdix + enddo + ! write(6,*) " Convergence Factor = ",nd2,xdix,deps_tracking + ! pause 123321 + ! if(verbose) write(6,*) " Convergence Factor = ",xdix + if(xdix.gt.deps_tracking) then + ite=1 + else + if(xdix.ge.xdix0.or.xdix<=tiny) then + ite=0 + else + ite=1 + xdix0=xdix + endif + endif + + if(iteM>=MAX_FIND_ITER) then + ! C_%stable_da=.FALSE. + ! IF(iteM==MAX_FIND_ITER+100) THEN + ! write(6,*) " Unstable in find_orbit without TPSA" + messagelost= "Maximum number of iterations in find_orbit without TPSA" + xlost=fix + check_stable=my_false + ! ENDIF + ITE=0 + ! if(try>0) goto 1111 + return + endif + ! write(6,*) item,xdix,xdix0,tiny + + if(ite.eq.1.or.item null() real(dp), POINTER ::L => null() real(dp), DIMENSION(:), POINTER :: AN => null(),BN => null() !Multipole component - real(dp), POINTER:: FINT => null(),HGAP => null() !FRINGE FUDGE FOR MAD + real(dp), DIMENSION(:), POINTER:: FINT => null(),HGAP => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: H1 => null(),H2 => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: VA => null(),VS => null() !valishev-like multipole integer, pointer :: f => null() @@ -56,7 +56,7 @@ TYPE DKD2P TYPE(MAGNET_CHART), POINTER:: P => null() TYPE(REAL_8), POINTER ::L => null() TYPE(REAL_8), DIMENSION(:), POINTER :: AN => null(),BN => null() !Multipole component - TYPE(REAL_8), POINTER:: FINT => null(),HGAP => null() !FRINGE FUDGE FOR MAD + TYPE(REAL_8), DIMENSION(:), POINTER:: FINT => null(),HGAP => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: H1 => null(),H2 => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: VA => null(),VS => null() !valishev-like multipole integer, pointer :: f => null() @@ -170,7 +170,7 @@ TYPE SOL5 real(dp), POINTER ::L => null() real(dp), POINTER ::B_SOL => null() real(dp), DIMENSION(:), POINTER :: AN => null(),BN => null() !Multipole component - real(dp), POINTER:: FINT => null(),HGAP => null() !FRINGE FUDGE FOR MAD + real(dp), DIMENSION(:), POINTER:: FINT => null(),HGAP => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: H1 => null(),H2 => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: VA => null(),VS => null() ! sad f1,f2 real(dp), POINTER:: dx => null(),dy => null(),pitch_x => null(),pitch_y => null() ! SADISTIC @@ -181,7 +181,7 @@ TYPE SOL5P TYPE(REAL_8), DIMENSION(:), POINTER :: AN => null(), BN => null() !Multipole component TYPE(REAL_8), POINTER ::L => null() TYPE(REAL_8), POINTER ::B_SOL => null() - TYPE(REAL_8), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD + TYPE(REAL_8), DIMENSION(:), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: H1 => null(), H2 => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: VA => null(), VS => null() !valishev-like multipole real(dp), POINTER:: dx => null(),dy => null(),pitch_x => null(),pitch_y => null() ! SADISTIC @@ -195,7 +195,7 @@ TYPE KTK real(dp), DIMENSION(:,:), POINTER :: MATX => null(), MATY => null() !LINEAR MATRIX !frs real(dp), DIMENSION(:), POINTER :: lx(:) => null(), ly(:) => null() real(dp), DIMENSION(:), POINTER :: lx => null(), ly => null() - real(dp), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD + real(dp), DIMENSION(:), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: H1 => null(), H2 => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: VA => null(), VS => null() !valishev-like multipole END TYPE KTK @@ -207,7 +207,7 @@ TYPE KTKP TYPE(REAL_8), DIMENSION(:,:), POINTER :: MATX => null(), MATY => null() !LINEAR MATRIX !frs TYPE(REAL_8), DIMENSION(:), POINTER :: lx(:) => null(), ly(:) => null() TYPE(REAL_8), DIMENSION(:), POINTER :: lx => null(), ly => null() - TYPE(REAL_8), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD + TYPE(REAL_8), DIMENSION(:), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: H1 => null(), H2 => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: VA => null(), VS => null() !valishev-like multipole END TYPE KTKP @@ -225,7 +225,7 @@ TYPE TKTF real(dp), DIMENSION(:), POINTER :: Rlx => null() ! real(dp), DIMENSION(:), POINTER :: dx(:) => null() ! real(dp), DIMENSION(:), POINTER :: dy(:) => null() - real(dp), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD + real(dp), DIMENSION(:), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: H1 => null(), H2 => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: VA => null(), VS => null() !valishev-like multipole integer, pointer :: f => null() @@ -244,7 +244,7 @@ TYPE TKTFP TYPE(REAL_8), DIMENSION(:), POINTER :: Rlx => null() ! real(dp), DIMENSION(:), POINTER :: dx(:) => null() ! real(dp), DIMENSION(:), POINTER :: dy(:) => null() - TYPE(REAL_8), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD + TYPE(REAL_8), DIMENSION(:), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: H1 => null(), H2 => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: VA => null(), VS => null() !valishev-like multipole integer, pointer :: f => null() @@ -278,7 +278,7 @@ TYPE TEAPOT real(dp), DIMENSION(:), POINTER :: bf_x => null(),bf_y => null() ! B field polynomial ! INTEGER,POINTER :: SECTOR_NMUL => null() logical(lp), POINTER :: DRIFTKICK => null() ! Split flag - real(dp), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD + real(dp), DIMENSION(:), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: H1 => null(), H2 => null() !FRINGE FUDGE FOR MAD integer, pointer :: f => null() real(dp), POINTER:: VA => null(), VS => null() !valishev-like multipole @@ -295,7 +295,7 @@ TYPE TEAPOTP TYPE(REAL_8), DIMENSION(:), POINTER :: bf_x => null(),bf_y => null() ! B field polynomial ! INTEGER,POINTER :: SECTOR_NMUL => null() logical(lp), POINTER :: DRIFTKICK => null() - TYPE(REAL_8), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD + TYPE(REAL_8), DIMENSION(:), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: H1 => null(), H2 => null() !FRINGE FUDGE FOR MAD integer, pointer :: f => null() TYPE(REAL_8), POINTER:: VA => null(), VS => null() !valishev-like multipole @@ -360,7 +360,7 @@ TYPE STREX real(dp), POINTER ::L => null() real(dp), DIMENSION(:), POINTER :: AN => null(), BN => null() !Multipole component logical(lp), POINTER :: DRIFTKICK => null(), LIKEMAD => null() - real(dp), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD + real(dp), DIMENSION(:), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: H1 => null(), H2 => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: VA => null(), VS => null() !valishev-like multipole integer, pointer :: f => null() @@ -371,7 +371,7 @@ TYPE STREXP TYPE(REAL_8), POINTER ::L => null() TYPE(REAL_8), DIMENSION(:), POINTER :: AN => null(), BN => null() !Multipole component logical(lp), POINTER :: DRIFTKICK => null(), LIKEMAD => null() - TYPE(REAL_8), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD + TYPE(REAL_8), DIMENSION(:), POINTER:: FINT => null(), HGAP => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: H1 => null(), H2 => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: VA => null(), VS => null() !valishev-like multipole integer, pointer :: f => null() diff --git a/libs/ptc/src/a_def_element_fibre_layout.inc b/libs/ptc/src/a_def_element_fibre_layout.inc index 5e06a96c5..6af1f68df 100644 --- a/libs/ptc/src/a_def_element_fibre_layout.inc +++ b/libs/ptc/src/a_def_element_fibre_layout.inc @@ -182,7 +182,7 @@ TYPE ELEMENT real(dp), POINTER :: L => null() ! Length of integration often same as LD ! real(dp), DIMENSION(:), POINTER:: AN => null(),BN => null() !Multipole component - real(dp), POINTER:: FINT => null(),HGAP => null() !FRINGE FUDGE FOR MAD + real(dp), DIMENSION(:), POINTER:: FINT => null(),HGAP => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: H1 => null(),H2 => null() !FRINGE FUDGE FOR MAD real(dp), POINTER:: VA => null(),VS => null() ! use in quad fringe from sad ! @@ -248,7 +248,7 @@ TYPE ELEMENTP ! TYPE(REAL_8), POINTER :: L => null() ! LENGTH OF INTEGRATION OFTEN SAME AS LD, CAN BE ZERO TYPE(REAL_8), DIMENSION(:), POINTER :: AN => null(),BN => null() !MULTIPOLE COMPONENT - TYPE(REAL_8), POINTER:: FINT => null(),HGAP => null() !FRINGE FUDGE FOR MAD + TYPE(REAL_8), DIMENSION(:), POINTER:: FINT => null(),HGAP => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: H1 => null(),H2 => null() !FRINGE FUDGE FOR MAD TYPE(REAL_8), POINTER:: VA => null(),VS => null() !valishev-like multipole ! @@ -342,7 +342,8 @@ END TYPE FIBRE TYPE LAYOUT CHARACTER(120), POINTER :: NAME => null()! IDENTIFICATION - INTEGER, POINTER :: INDEX,HARMONIC_NUMBER => null()! IDENTIFICATION, CHARGE SIGN + INTEGER, POINTER :: INDEX => null()! IDENTIFICATION, CHARGE SIGN + REAL(DP), POINTER :: HARMONIC_NUMBER => null()! logical(lp),POINTER ::CLOSED => null() INTEGER, POINTER :: N => null() ! TOTAL ELEMENT IN THE CHAIN INTEGER,POINTER ::NTHIN => null() ! NUMBER IF THIN LENSES IN COLLECTION (FOR SPEED ESTIMATES) diff --git a/src/madx_ptc_module.f90 b/src/madx_ptc_module.f90 index fffccd2b3..41f2e317a 100644 --- a/src/madx_ptc_module.f90 +++ b/src/madx_ptc_module.f90 @@ -52,6 +52,7 @@ MODULE madx_ptc_module type(fibre), pointer :: p => null() end type fibreptr + integer, private, parameter:: maxelperclock = 10 ! maximum 10 ac dipols with given clock type clockdef real(dp) :: tune = -1 ! negative means inactive, in fact it is tune, left like this for backward compatibility, can be changed during LS2 @@ -784,7 +785,10 @@ subroutine ptc_input() endif endif case(3) ! PTC accepts mults watch out sector_nmul defaulted to 4 + + if (getdebug()>2) print*,"Translating SBEND" if(l.eq.zero) then + if (getdebug()>2) print*,"Length zero -> translating as MARKER" key%magnet="marker" goto 100 endif @@ -846,6 +850,18 @@ subroutine ptc_input() call augment_count('errors_dipole ') endif endif + + if (getdebug()>2) then + print*,"B0=", key%list%b0 + print*,"K=", key%list%k + print*,"KS=", key%list%ks + print*,"TILT=", key%tiltd + print*,"T1=", key%list%t1 + print*,"T2=", key%list%t2 + print*,"H1=", key%list%h1 + print*,"H2=", key%list%h2 + endif + case(5) key%magnet="quadrupole" !VK @@ -950,6 +966,7 @@ subroutine ptc_input() !================================================================ case(8) + if (getdebug()>2) print*,"Translating MULTIPOLE" key%magnet="multipole" !---- Multipole components. F_ERRORS = zero @@ -1032,6 +1049,14 @@ subroutine ptc_input() endif endif + if (getdebug()>2) then + print*,"thin_h_angle=", key%list%thin_h_angle + print*,"thin_v_angle=", key%list%thin_v_angle + print*,"K=", key%list%k + print*,"KS=", key%list%ks + print*,"TILT=", key%tiltd + endif + case(9) ! PTC accepts mults key%magnet="solenoid" @@ -1458,6 +1483,7 @@ subroutine ptc_input() endif + if(advance_node().ne.0) goto 10 @@ -1503,7 +1529,7 @@ subroutine ptc_input() write(6,*) "Before start: ",my_ring%start%chart%f%a write(6,*) "Before end: ",my_ring%end%chart%f%b endif - + call survey(my_ring) if (getdebug() > 0) then @@ -1882,6 +1908,7 @@ subroutine ptc_align() type(fibre), pointer :: f !--------------------------------------------------------------- + j=restart_sequ() j=0 f=>my_ring%start @@ -1891,60 +1918,150 @@ subroutine ptc_align() n_align = node_al_errors(al_errors) if (n_align.ne.0) then if (getdebug() > 3) then + write(6,*) " ----------------------------------------------- " write(6,*) f%mag%name," Translation Error " write(6,'(3f11.8)') al_errors(1:3) write(6,*) f%mag%name," Rotation Error " write(6,'(3f11.8)') al_errors(4:6) - - write(6,*) - write(6,*) "Ac:", f%chart%f%a - write(6,*) "Oc:", f%chart%f%o - write(6,*) "Bc:", f%chart%f%b - write(6,*) "Am:", f%mag%p%f%a - write(6,*) "Om:", f%mag%p%f%o - write(6,*) "Bc:", f%mag%p%f%b - write(6,*) - write(6,*) f%mag%p%f%ent(1,:) - write(6,*) f%mag%p%f%ent(2,:) - write(6,*) f%mag%p%f%ent(3,:) - + call print_elframes(f) endif - ! this routine is buggy + ! this routine is buggy -> fixed by Etienne on 2018.01.29 ! call mad_misalign_fibre(f,al_errors(1:6)) - ! this is PTC original, but it shifts in frame after roations + ! this is PTC original, but it first rotates and than shifts (in frame after rotations) ! call misalign_fibre(f,al_errors(1:6)) !our new routine call misalign_element(f,al_errors) + !a workaround to handle misaligment of thin dipoles that is not handled by PTC + call misalign_thindipole(f,al_errors) + + if (getdebug() > 3) then - write(6,*) - write(6,*) "Ac:", f%chart%f%a - write(6,*) "Oc:", f%chart%f%o - write(6,*) "Bc:", f%chart%f%b - - write(6,*) "Am:", f%mag%p%f%a - write(6,*) "Om:", f%mag%p%f%o - write(6,*) "Bc:", f%mag%p%f%b - - - write(6,*) - write(6,*) f%mag%p%f%ent(1,:) - write(6,*) f%mag%p%f%ent(2,:) - write(6,*) f%mag%p%f%ent(3,:) - + write(6,*) " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv " + call print_elframes(f) endif endif f=>f%next if(advance_node().ne.0) goto 10 + END subroutine ptc_align !_________________________________________________________________ + + subroutine print_elframes(f) + implicit none + TYPE(FIBRE),target,INTENT(INOUT):: f + + write(6,*) "Ac:", f%chart%f%a + write(6,*) "Oc:", f%chart%f%o + write(6,*) "Bc:", f%chart%f%b + write(6,*) + write(6,*) "Am:", f%mag%p%f%a + write(6,*) "Om:", f%mag%p%f%o + write(6,*) "Bm:", f%mag%p%f%b + write(6,*) + write(6,*) "entm(1,:) :", f%mag%p%f%ent(1,:) + write(6,*) "entm(2,:) :", f%mag%p%f%ent(2,:) + write(6,*) "entm(3,:) :", f%mag%p%f%ent(3,:) + write(6,*) + write(6,*) "midm(1,:) :", f%mag%p%f%mid(1,:) + write(6,*) "midm(2,:) :", f%mag%p%f%mid(2,:) + write(6,*) "midm(3,:) :", f%mag%p%f%mid(3,:) + write(6,*) + write(6,*) "exim(1,:) :", f%mag%p%f%exi(1,:) + write(6,*) "exim(2,:) :", f%mag%p%f%exi(2,:) + write(6,*) "exim(3,:) :", f%mag%p%f%exi(3,:) + write(6,*) + write(6,*) "ang_in :", f%CHART%ang_in + write(6,*) "ang_out :", f%CHART%ang_out + + end subroutine print_elframes + + !_____________________________________________________ + ! Routine to handle thin dipole misalignments + ! Etienne says that KIND3 is "highly illigal" and he does not support it + ! therefore we need to tweak ut ourselves + ! Patching from rotation errors is done using FIBRE%CHART%ANG_OUT and FIBRE%CHART%ANG_OUT + ! Here we calculate these angles + subroutine misalign_thindipole(f,al_errors) + use twiss0fi, only: align_max + TYPE(FIBRE),target,INTENT(INOUT):: f + REAL(DP),INTENT(IN) :: al_errors(align_max) + REAL(DP) :: Fent(3,3), F0ent(3,3), Fexi(3,3), F0exi(3,3), A + REAL(DP) :: F1(3,3), F2(3,3), F3(3,3) + + if (f%mag%kind /= kind3 ) return + + !in MADX vertical bend is tilted horizontal bend, so thin_v_angle is always zero + if ( abs(f%mag%K3%thin_h_angle) .lt. 1e-12 ) return; + + if (getdebug() > 2) then + write(6,*) "misalign_thindipole angle = ",f%mag%K3%thin_h_angle + endif + + + A = f%mag%K3%thin_h_angle + ! identity + F0ent = zero + F0ent(1,1) = 1 + F0ent(2,2) = 1 + F0ent(3,3) = 1 + + !!!!!!!!!!!!! + ! depends only on the bend angle + F0exi = zero + F0exi(1,1) = cos(A) + F0exi(1,3) = sin(A) + F0exi(2,2) = 1 + F0exi(3,1) = -sin(A) + F0exi(3,3) = cos(A) + + !!!!!!!!!!!!! + ! Lumped error matrix + f1 = zero + f2 = zero + f3 = zero + + f1(1,1) = 1 + f1(2,2) = cos(al_errors(4)) + f1(2,3) =-sin(al_errors(4)) + f1(3,2) = sin(al_errors(4)) + f1(3,3) = cos(al_errors(4)) + + f2(1,1) = cos(al_errors(5)) + f2(1,3) =-sin(al_errors(5)) + f2(2,2) = 1 + f2(3,1) = sin(al_errors(5)) + f2(3,3) = cos(al_errors(5)) + + ! change of sign + f3(1,1) = cos(-al_errors(6)) + f3(1,2) =-sin(-al_errors(6)) + f3(2,1) = sin(-al_errors(6)) + f3(2,2) = cos(-al_errors(6)) + f3(3,3) = 1 + + Fent = matmul(f2, f1) + Fent = matmul(f3,Fent) + + !____________________________________ + + !lumped bend rot with alignment rotation + Fexi = matmul(F0exi,Fent) + + ! Calculate athe angles needed for tracking (a PTC routine) + CALL COMPUTE_ENTRANCE_ANGLE(F0ENT,FENT,f%CHART%ANG_IN) + CALL COMPUTE_ENTRANCE_ANGLE(FEXI,F0EXI,f%CHART%ANG_OUT) + + end subroutine misalign_thindipole + + !_____________________________________________________ subroutine misalign_element(f,al_errors) use twiss0fi, only: align_max TYPE(FIBRE),target,INTENT(INOUT):: f @@ -3331,6 +3448,4 @@ subroutine acdipoleramping(t) end subroutine acdipoleramping - - END MODULE madx_ptc_module diff --git a/src/madx_ptc_twiss.f90 b/src/madx_ptc_twiss.f90 index 66920812a..34c7735f6 100644 --- a/src/madx_ptc_twiss.f90 +++ b/src/madx_ptc_twiss.f90 @@ -566,7 +566,6 @@ subroutine ptc_twiss(tab_name,summary_tab_name) logical(lp) :: doTMtrack ! true if rmatrix==true .and. isRing==true . do not track theTransferMap and save time ! .or. already tracked form closed solution search logical(lp) :: usertableActive = .false. ! flag to mark that there was something requested with ptc_select - integer :: countSkipped character(48) :: summary_table_name character(12) :: tmfile='transfer.map' @@ -1035,7 +1034,7 @@ subroutine ptc_twiss(tab_name,summary_tab_name) do i=1,MY_RING%n - if (getdebug() > 1) then + if (getdebug() > 2) then write(6,*) "" write(6,*) "##########################################" write(6,'(i4, 1x,a, f10.6)') i,current%mag%name, suml diff --git a/tests/share/ALS/als.seqx b/tests/share/ALS/als.seqx index dd47cc643..8d1ef1c81 100644 --- a/tests/share/ALS/als.seqx +++ b/tests/share/ALS/als.seqx @@ -52,7 +52,9 @@ BEND1 : RBEND,L= LBEND, ANGLE=ALPHA, k1=-0.778741; CAVM:MARKER; rfvolt = 0.2d0; -CAV:RFCAVITY,L=0.2000,VOLT:=rfvolt,FREQ=500., LAG=0.25; +rffreq = 500; +rflag = 0.25; +CAV:RFCAVITY,L=0.2000,VOLT:=rfvolt,FREQ:=rffreq, LAG:=rflag; sfline: line = (1*sf); sdline: line = (1*sd); diff --git a/tests/test-ptc-twiss-6D-ALS/test-ptc-twiss-6D-ALS.madx b/tests/test-ptc-twiss-6D-ALS/test-ptc-twiss-6D-ALS.madx index cc9081a75..eee60022b 100644 --- a/tests/test-ptc-twiss-6D-ALS/test-ptc-twiss-6D-ALS.madx +++ b/tests/test-ptc-twiss-6D-ALS/test-ptc-twiss-6D-ALS.madx @@ -1,4 +1,7 @@ +!example of 6D calculation on ALS (synchrotron light source) +option, -echo; call, file="../share/ALS/als.seqx"; +option, echo; beam, particle = electron, energy = sqrt(1.5*1.5+emass*emass); use, period=ALS; diff --git a/tests/test-ptc-twiss-6D-ALS/test-ptc-twiss-6D-ALS.ref b/tests/test-ptc-twiss-6D-ALS/test-ptc-twiss-6D-ALS.ref index b2a4b3642..f48c47a99 100644 --- a/tests/test-ptc-twiss-6D-ALS/test-ptc-twiss-6D-ALS.ref +++ b/tests/test-ptc-twiss-6D-ALS/test-ptc-twiss-6D-ALS.ref @@ -3,165 +3,11 @@ + MAD-X 5.03.07 (64 bit, Linux) + + Support: mad@cern.ch, http://cern.ch/mad + + Release date: 2017.10.20 + - + Execution date: 2018.01.26 14:45:42 + + + Execution date: 2018.01.31 12:23:26 + ++++++++++++++++++++++++++++++++++++++++++++ -call, file="../share/ALS/als.seqx"; - -L1 : drift, L= 2.832695; - -L2 : drift, L= 0.45698; - -L3 : drift, L= 0.08902; - -L4 : drift, L= 0.2155; - -L5 : drift, L= 0.219; - -L6 : drift, L= 0.107078; - -L7 : drift, L= 0.105716; - -L8 : drift, L= 0.135904; - -L9 : drift, L= 0.2156993; - -L10: drift, L= 0.089084; - -L11: drift, L= 0.235416; - -L12: drift, L= 0.1245; - -L13: drift, L= 0.511844; - -L14: drift, L= 0.1788541; - -L15: drift, L= 0.1788483; - -L16: drift, L= 0.511849; - -L17: drift, L= 0.1245; - -L18: drift, L= 0.235405; - -L19: drift, L= 0.089095; - -L20: drift, L= 0.2157007; - -L21: drift, L= 0.177716; - -L22: drift, L= 0.170981; - -L23: drift, L= 0.218997; - -L24: drift, L=0.215503; - -L25: drift, L=0.0890187; - -L26: drift, L=0.45698; - -L27: drift, L=2.832696; - -L27c: drift, L=2.832696-0.2; - -ds : drift, L=0.1015; - - - -QF1 : QUADRUPOLE,L=0.344, K1= 2.2474D0+6.447435260914397e-03; - -QF2 : QUADRUPOLE,L=0.344, K1= 2.2474; - -QD1 : QUADRUPOLE,L=0.187, K1= -2.3368D0-2.593018157427161e-02; - -QD2 : QUADRUPOLE,L=0.187, K1= -2.3368; - -QFA1: QUADRUPOLE,L=0.448, K1= 2.8856; - -QFA2: QUADRUPOLE,L=0.448, K1= 2.8856; - - - -!!! 1/2 mad-x value - -ksf=-41.3355516397069748d0; - -ksd=56.2564709584745489d0; - - - -sf:sextupole,l=2*0.1015d0, K2= ksf; - -sd: sextupole,l=2*0.1015d0, K2= ksd; - - - -VC5:marker; - -ALPHA=0.17453292519943295769236907684886d0; - - - -LBEND=0.86621d0; - - - -BEND : RBEND, L=LBEND, ANGLE=ALPHA, k1=-0.778741; - -BEND1 : RBEND,L= LBEND, ANGLE=ALPHA, k1=-0.778741; - - - -CAVM:MARKER; - -rfvolt = 0.2d0; - -CAV:RFCAVITY,L=0.2000,VOLT:=rfvolt,FREQ=500., LAG=0.25; - - - - sfline: line = (1*sf); - - sdline: line = (1*sd); - - - -SUP1: line =(L1+L2+L3+QF1+VC5+L4+L5+QD1+L6+L7+L8+VC5+BEND+VC5+L9+sfline+L10+ - - L11+QFA1+L12+sdline+L13+ - - L14+BEND+L15+L16+sdline+L17+ - - QFA2+L18+L19+sfline+L20+BEND+L21+ - - L22+QD2+L23+L24+QF2+L25+ - - L26+VC5+L27+cavm); - - - -SUPb: line=(L1+L2+L3+QF1+VC5+L4+L5+QD1+L6+L7+L8+VC5+BEND+VC5+L9+sfline+L10+ - - L11+QFA1+L12+sdline+L13+ - - L14+BEND+L15+L16+sdline+L17+ - - QFA2+L18+L19+sfline+L20+BEND1+L21+ - - L22+QD2+L23+L24+QF2+L25+ - - L26+VC5+L27c+cav); - - - - - -SSTART: marker; !to get compatible with lattice build in F90 code (pure PTC example) - - - -ALS: line = (SSTART + 11*sup1 + supb); - +!example of 6D calculation on ALS (synchrotron light source) +option, -echo;