source: PSPA/madxPSPA/libs/ptc/src/Sma_multiparticle.f90 @ 430

Last change on this file since 430 was 430, checked in by touze, 11 years ago

import madx-5.01.00
File size: 64.5 KB
Line 
1
2!The Polymorphic Tracking Code
3!Copyright (C) Etienne Forest and CERN
4
5module ptc_multiparticle
6  !  use S_TRACKING  !,FRINGE_=>FRINGE__MULTI !,FACE=>FACE_MULTI
7  use beam_beam_ptc
8  implicit none
9  public
10  CHARACTER*27 CASE_NAME(-2:3)
11  PRIVATE fuzzy_eq,fuzzy_neq
12  PRIVATE TRACK_FIBRE_FRONTR,TRACK_FIBRE_FRONTP
13  PRIVATE TRACK_FIBRE_BACKR,TRACK_FIBRE_BACKP
14  PRIVATE TRACKR_NODE_SINGLE,TRACKP_NODE_SINGLE,TRACKV_NODE_SINGLE
15  private DRIFTr_BACK_TO_POSITION,DRIFTp_BACK_TO_POSITION  !,DRIFT_BACK_TO_POSITION
16  private MAKE_NODE_LAYOUT_2 !,DRIFT_TO_TIME
17  PRIVATE MODULATE_R,MODULATE_P
18  PRIVATE TRACK_MODULATION_R,TRACK_MODULATION_P
19
20  !  LOGICAL :: OLD_MOD=.TRUE.
21
22  logical(lp),private, parameter :: dobb=.true.
23  logical(lp),private, parameter :: aperture_all_case0=.false.
24  type(probe) :: xsm,xsm0
25  !real(dp) :: unit_time =1.0e-3_dp
26  REAL(dp) :: x_orbit_sync(6)= 0.0_dp,dt_orbit_sync=0.0_dp
27 
28  INTERFACE TRACK_NODE_SINGLE
29     MODULE PROCEDURE TRACKR_NODE_SINGLE     !@1  t,x,state,charge
30     MODULE PROCEDURE TRACKP_NODE_SINGLE     !@1  t,y,state,charge
31     MODULE PROCEDURE TRACKV_NODE_SINGLE     !@1  t,v,state,charge
32  END INTERFACE
33
34
35
36
37  INTERFACE DRIFT_BACK_TO_POSITION
38     MODULE PROCEDURE DRIFTr_BACK_TO_POSITION
39     MODULE PROCEDURE DRIFTp_BACK_TO_POSITION
40  END INTERFACE
41
42  INTERFACE TRACK_FIBRE_FRONT
43     MODULE PROCEDURE TRACK_FIBRE_FRONTR
44     MODULE PROCEDURE TRACK_FIBRE_FRONTP
45  END INTERFACE
46
47  INTERFACE TRACK_FIBRE_BACK
48     MODULE PROCEDURE TRACK_FIBRE_BACKR
49     MODULE PROCEDURE TRACK_FIBRE_BACKP
50  END INTERFACE
51
52  INTERFACE COPY
53     MODULE PROCEDURE COPY_BEAM
54  END INTERFACE
55
56  INTERFACE OPERATOR (.feq.)
57     MODULE PROCEDURE fuzzy_eq
58  END INTERFACE
59  INTERFACE OPERATOR (.fne.)
60     MODULE PROCEDURE fuzzy_neq
61  END INTERFACE
62
63  INTERFACE MODULATE
64     MODULE PROCEDURE MODULATE_R   ! PROPAGATE FAKE MODULATED (Q,P) ; 7TH AND 8TH VARIABLES
65     MODULE PROCEDURE MODULATE_P   !
66  END INTERFACE
67
68
69  INTERFACE TRACK_MODULATION
70     MODULE PROCEDURE TRACK_MODULATION_R   ! PROPAGATE FAKE MODULATED (Q,P) ; 7TH AND 8TH VARIABLES
71     MODULE PROCEDURE TRACK_MODULATION_P   !
72  END INTERFACE
73
74  type three_d_info
75     !   character(nlp),pointer :: name
76     real(dp)  a(3),b(3)   ! Centre of entrance and exit faces
77     real(dp)  ent(3,3),exi(3,3)  ! entrace and exit frames for drawing magnet faces
78     real(dp)  wx,wy ! width of box for plotting purposes
79     real(dp)  o(3),mid(3,3)   ! frames at the point of tracking
80     real(dp)  reference_ray(6)  !
81     real(dp) x(6)   ! ray tracked with reference_ray using a  type(beam)
82     real(dp) r0(3),r(3)  ! ray position global returned
83     real(dp) scale     !  magnification using reference_ray
84     logical(lp) u(2)   ! unstable flag for both ray and reference_ray
85  END type three_d_info
86
87!  real :: ttime0,ttime1,dt1=0.0,dt2=0.0;
88
89CONTAINS
90
91  SUBROUTINE MODULATE_R(C,XS,K)
92    IMPLICIT NONE
93    type(INTEGRATION_NODE), pointer :: C
94    type(probe), INTENT(INOUT) :: xs
95    TYPE(INTERNAL_STATE) K
96    TYPE(ELEMENT),POINTER :: EL
97    TYPE(ELEMENTP),POINTER :: ELp
98    REAL(DP) v,dv
99
100
101
102    EL=>C%PARENT_FIBRE%MAG
103    ELP=>C%PARENT_FIBRE%MAGP
104
105    IF(K%MODULATION) THEN
106
107       DV=(XS%AC%X(1)*COS(EL%theta_ac)-XS%AC%X(2)*SIN(EL%theta_ac))
108       V=EL%DC_ac+EL%A_ac*DV
109       DV=el%D_ac*DV
110     else
111       V=0.0_dp
112       DV=0.0_dp
113     endif
114
115    CALL transfer_ANBN(EL,ELP,VR=V,DVR=DV)
116
117
118  END   SUBROUTINE MODULATE_R
119 
120   SUBROUTINE do_ramping_R(C,t,K)
121    IMPLICIT NONE
122    type(INTEGRATION_NODE), pointer :: C
123    TYPE(INTERNAL_STATE) K
124    TYPE(ELEMENT),POINTER :: EL
125    TYPE(ELEMENTP),POINTER :: ELp
126    REAL(DP) v,dv,t
127
128
129
130    EL=>C%PARENT_FIBRE%MAG
131    ELP=>C%PARENT_FIBRE%MAGP
132    if(.not.associated(EL%ramp)) return
133
134       V=EL%DC_ac
135       DV=0.0_dp
136       call set_ramp(C,t)
137
138    CALL transfer_ANBN(EL,ELP,VR=V,DVR=DV)
139
140
141  END   SUBROUTINE do_ramping_R
142 
143   SUBROUTINE DO_Ramping_p(C,t,K)
144    IMPLICIT NONE
145    type(INTEGRATION_NODE), pointer :: C
146    TYPE(INTERNAL_STATE) K
147    TYPE(ELEMENT),POINTER :: EL
148    TYPE(ELEMENTP),POINTER :: ELP
149    TYPE(REAL_8) V,DV
150    real(dp) t
151
152    EL=>C%PARENT_FIBRE%MAG
153    ELP=>C%PARENT_FIBRE%MAGP
154   
155    if(.not.associated(EL%ramp)) return
156
157    CALL ALLOC(V)
158    CALL ALLOC(DV)
159
160
161       V=elp%DC_ac
162       DV=0.0_dp
163         call set_ramp(C,t)
164
165    CALL transfer_ANBN(EL,ELP,VP=V,DVP=DV)
166
167    CALL KILL(V)
168    CALL KILL(DV)
169
170
171  END   SUBROUTINE DO_Ramping_p
172
173   SUBROUTINE set_all_ramp(R)
174    IMPLICIT NONE
175    TYPE(layout), target  :: r
176    TYPE(fibre), POINTER  :: p
177    integer i
178    REAL(DP) v,dv
179    v=0.0_dp
180    dv=0.0_dp     
181     p=>r%start
182     do i=1,r%n
183     
184      if(associated(p%mag%ramp)) then
185       call set_ramp(p%t1,x_orbit_sync(6))
186       CALL transfer_ANBN(p%mag,p%magp,VR=V,DVR=DV)
187
188      endif
189     
190      p=>p%next
191     enddo     
192     
193    end SUBROUTINE set_all_ramp
194     
195  SUBROUTINE set_ramp(t,t0)
196    IMPLICIT NONE
197    TYPE(INTEGRATION_NODE), POINTER  :: T
198    integer i,it
199    real(dp) r,ti,rat,dtot
200    type(ramping), pointer :: a
201    real(dp) an,bn,t0
202   
203
204 !   if(t%pos_in_fibre==1) return
205   
206    a=>t%parent_fibre%mag%ramp
207   
208
209    dtot=(a%table(a%n)%time-a%table(1)%time)  !/(a%n-1)
210   
211 !   ti=XSM0%ac%t/clight/a%unit_time    ! time in milliseconds
212    ti=t0/clight    !/a%unit_time    ! time in milliseconds
213   
214!    if(ti>a%t_max.or.ti<a%table(1)%time) then
215!    if(ti>a%table(a%n)%time.or.ti<a%table(1)%time) then
216!     return
217!    endif
218     
219
220     
221    if(ti>=a%t_max.or.ti<a%table(1)%time) then
222!    if(ti>a%table(a%n)%time.or.ti<a%table(1)%time) then
223      if(ti>=a%t_max) then
224           a%table(0)%bn=0.0_dp
225           a%table(0)%an=0.0_dp
226          do i=1,size(a%table(0)%bn)
227           a%table(0)%bn(i)= a%table(a%n)%bn(i)*a%r
228           a%table(0)%an(i)= a%table(a%n)%an(i)*a%r
229          enddo
230            a%table(0)%b_t= a%table(a%n)%b_t
231           a=>t%parent_fibre%magp%ramp
232           a%table(0)%bn=0.0_dp
233           a%table(0)%an=0.0_dp
234          do i=1,size(a%table(0)%bn)
235           a%table(0)%bn(i)= a%table(a%n)%bn(i)*a%r
236           a%table(0)%an(i)= a%table(a%n)%an(i)*a%r
237          enddo
238          a%table(0)%b_t= a%table(a%n)%b_t
239      else
240           a%table(0)%bn=0.0_dp
241           a%table(0)%an=0.0_dp
242          do i=1,size(a%table(0)%bn)
243           a%table(0)%bn(i)= a%table(1)%bn(i)*a%r
244           a%table(0)%an(i)= a%table(1)%an(i)*a%r
245          enddo
246         a%table(0)%b_t= a%table(1)%b_t
247          a=>t%parent_fibre%magp%ramp
248           a%table(0)%bn=0.0_dp
249           a%table(0)%an=0.0_dp
250          do i=1,size(a%table(0)%bn)
251           a%table(0)%bn(i)= a%table(1)%bn(i)*a%r
252           a%table(0)%an(i)= a%table(1)%an(i)*a%r
253          enddo
254          a%table(0)%b_t= a%table(1)%b_t
255     
256      endif
257
258    else
259   
260         ti=ti-a%table(1)%time
261         ti=mod(ti,dtot)+a%table(1)%time
262          dtot=dtot/(a%n-1)
263          ti=(ti-a%table(1)%time)/dtot+1
264           
265          it=int(ti)
266!          it=idint(ti)
267         
268           rat=(ti-it)   
269           
270
271           a%table(0)%bn=0.0_dp
272           a%table(0)%an=0.0_dp
273          do i=1,size(a%table(0)%bn)
274           a%table(0)%bn(i)=((a%table(it+1)%bn(i)-a%table(it)%bn(i))*rat + a%table(it)%bn(i))*a%r
275           a%table(0)%an(i)= ((a%table(it+1)%an(i)-a%table(it)%an(i))*rat + a%table(it)%an(i))*a%r
276          enddo
277           a%table(0)%b_t=((a%table(it+1)%b_t-a%table(it)%b_t)*rat + a%table(it)%b_t)
278
279          a=>t%parent_fibre%magp%ramp
280           a%table(0)%bn=0.0_dp
281           a%table(0)%an=0.0_dp
282          do i=1,size(a%table(0)%bn)
283           a%table(0)%bn(i)=((a%table(it+1)%bn(i)-a%table(it)%bn(i))*rat + a%table(it)%bn(i))*a%r
284           a%table(0)%an(i)= ((a%table(it+1)%an(i)-a%table(it)%an(i))*rat + a%table(it)%an(i))*a%r
285          enddo
286           a%table(0)%b_t=((a%table(it+1)%b_t-a%table(it)%b_t)*rat + a%table(it)%b_t)
287         
288    endif
289  end SUBROUTINE set_ramp
290 
291  SUBROUTINE MODULATE_P(C,XS,K)
292    IMPLICIT NONE
293    type(INTEGRATION_NODE), pointer :: C
294    type(probe_8), INTENT(INOUT) :: xs
295    TYPE(INTERNAL_STATE) K
296    TYPE(ELEMENT),POINTER :: EL
297    TYPE(ELEMENTP),POINTER :: ELP
298    TYPE(REAL_8) V,DV
299
300    EL=>C%PARENT_FIBRE%MAG
301    ELP=>C%PARENT_FIBRE%MAGP
302
303    CALL ALLOC(V)
304    CALL ALLOC(DV)
305
306    IF(K%MODULATION) THEN
307       DV=(XS%AC%X(1)*COS(ELP%theta_ac)-XS%AC%X(2)*SIN(ELP%theta_ac))
308       V=ELP%DC_ac+ELP%A_ac*DV
309       DV=elp%D_ac*DV
310
311       else  ! ramp
312          V=0.0_dp
313          DV=0.0_dp
314       endif
315 
316    CALL transfer_ANBN(EL,ELP,VP=V,DVP=DV)
317
318    CALL KILL(V)
319    CALL KILL(DV)
320
321
322  END   SUBROUTINE MODULATE_P
323
324 
325  SUBROUTINE TRACK_MODULATION_R(C,XS,K)
326    IMPLICIT NONE
327    type(INTEGRATION_NODE), pointer :: C
328    type(probe), INTENT(INOUT) :: xs
329    TYPE(INTERNAL_STATE) K
330    real(dp) xt
331    real(dp),pointer :: beta0
332
333    if(k%time) then
334       beta0=>C%PARENT_FIBRE%beta0
335       xs%ac%t=c%DS_AC/beta0+xs%ac%t
336       xt = cos(XS%AC%om * c%DS_AC/beta0) *XS%AC%X(1) + sin(XS%AC%om * c%DS_AC/beta0) *XS%AC%X(2)
337       XS%AC%X(2) = -sin(XS%AC%om * c%DS_AC/beta0) *XS%AC%X(1) + cos(XS%AC%om * c%DS_AC/beta0) *XS%AC%X(2)
338       XS%AC%X(1) = xt
339    else
340       xt = cos(XS%AC%om * c%DS_AC) *XS%AC%X(1) + sin(XS%AC%om * c%DS_AC) *XS%AC%X(2)
341       XS%AC%X(2) = -sin(XS%AC%om * c%DS_AC) *XS%AC%X(1) + cos(XS%AC%om * c%DS_AC) *XS%AC%X(2)
342       XS%AC%X(1) = xt
343       xs%ac%t=c%DS_AC+xs%ac%t
344    endif
345
346  END   SUBROUTINE TRACK_MODULATION_R
347
348  SUBROUTINE TRACK_MODULATION_P(C,XS,K)
349    IMPLICIT NONE
350    type(INTEGRATION_NODE), pointer :: C
351    type(probe_8), INTENT(INOUT) :: xs
352    TYPE(INTERNAL_STATE) K
353    TYPE(REAL_8) xt
354    real(dp),pointer :: beta0
355
356    CALL ALLOC(XT)
357
358    if(k%time) then
359       beta0=>C%PARENT_FIBRE%beta0
360       xs%ac%t=c%DS_AC/beta0+xs%ac%t
361       xt = cos(XS%AC%om * c%DS_AC/beta0) *XS%AC%X(1) + sin(XS%AC%om * c%DS_AC/beta0) *XS%AC%X(2)
362       XS%AC%X(2) = -sin(XS%AC%om * c%DS_AC/beta0) *XS%AC%X(1) + cos(XS%AC%om * c%DS_AC/beta0) *XS%AC%X(2)
363       XS%AC%X(1) = xt
364    else
365       xt = cos(XS%AC%om * c%DS_AC) *XS%AC%X(1) + sin(XS%AC%om * c%DS_AC) *XS%AC%X(2)
366       XS%AC%X(2) = -sin(XS%AC%om * c%DS_AC) *XS%AC%X(1) + cos(XS%AC%om * c%DS_AC) *XS%AC%X(2)
367       XS%AC%X(1) = xt
368       xs%ac%t=c%DS_AC+xs%ac%t
369    endif
370
371    CALL KILL(XT)
372
373  END   SUBROUTINE TRACK_MODULATION_P
374
375
376
377  FUNCTION fuzzy_eq( S1, S2 )
378    implicit none
379    logical(lp) fuzzy_eq
380    real(dp), INTENT (IN) :: S1,S2
381    fuzzy_eq=.false.
382
383    if(abs(s1-s2)<=c_%eps_pos) fuzzy_eq=.true.
384
385  end FUNCTION fuzzy_eq
386
387  FUNCTION fuzzy_neq( S1, S2 )
388    implicit none
389    logical(lp) fuzzy_neq
390    real(dp), INTENT (IN) :: S1,S2
391    fuzzy_neq=.false.
392
393    if(abs(s1-s2)>c_%eps_pos) fuzzy_neq=.true.
394
395  end FUNCTION fuzzy_neq
396
397  SUBROUTINE move_to_s( L,s,current,i,ds ) ! Moves position s
398    implicit none
399    TYPE (INTEGRATION_NODE), POINTER :: Current
400    TYPE (NODE_LAYOUT) L
401    real(dp) s,sp,ds
402    integer i,k
403    logical(lp) DOIT !,track_it
404
405    !  track_it=.false.
406    sp=mod(s,L%END%S(3))
407
408    if(sp==0.0_dp.and.s/=0.0_dp) then
409       current=>l%end
410       i=l%n+1
411       ds=0.0_dp
412       !  track_it=.true.
413       return
414    endif
415
416    if(sp==0.0_dp) then
417       current=>l%start
418       i=1
419       ds=0.0_dp
420       return
421    endif
422
423
424    nullify(current);
425    Current => L%LAST
426
427    k=L%LASTPOS
428    I=K
429    ds=0.0_dp
430
431    IF(SP>CURRENT%S(3) ) then
432
433       do i=k,l%n-1
434          if(current%next%s(3)>=sp) exit
435          current=>current%next
436       enddo
437
438       if(current%next%s(3)/=sp) ds=sp-current%s(3)
439       if(current%next%s(3)==sp) THEN
440          ds=0.0_dp
441          CURRENT=>current%next
442          I=I+1
443       ENDIF
444
445    elseif(SP<CURRENT%S(3)) then
446
447       do i=k-1,1,-1
448          current=>current%previous
449          if(current%s(3)<=sp) exit
450       enddo
451
452       ds=sp-current%s(3)
453
454    endif
455
456    L%LASTPOS=I; L%LAST => Current;
457
458    if(ds>0.0_dp) then
459       if(CURRENT%S(4)-ds.feq.0.0_dp) then
460          ds=0.0_dp
461          current=>Current%next
462          i=i+1
463          L%LAST => Current;
464       ELSEIF(ds.feq.0.0_dp) THEN
465          DS=0.0_dp
466       ENDIF
467    endif
468
469
470    DOIT=.TRUE.
471    !DOIT=.FALSE.
472
473    if(iabs(CURRENT%cas)==0.OR.iabs(CURRENT%cas)==1) then
474       do while(DS==0.0_dp.AND.DOIT)    ! PUTS AT BEGINNING IF DS=ZERO
475          CURRENT=>CURRENT%PREVIOUS
476          IF(ASSOCIATED(CURRENT)) THEN
477
478             IF((SP.FNE.CURRENT%S(3)).or.CURRENT%cas==-2) THEN
479                CURRENT=>CURRENT%NEXT
480                DOIT=.FALSE.
481             ELSE
482                I=I-1
483             ENDIF
484
485          ELSE
486             CURRENT=>CURRENT%NEXT
487             DOIT=.FALSE.
488          ENDIF
489
490       enddo
491    elseif(iabs(CURRENT%cas)==2) then
492       do while(DS==0.0_dp.AND.DOIT)    ! PUTS AT BEGINNING IF DS=ZERO
493          CURRENT=>CURRENT%next
494          IF(ASSOCIATED(CURRENT)) THEN
495
496             IF((SP.FNE.CURRENT%S(3)).or.CURRENT%cas==1) THEN
497                CURRENT=>CURRENT%previous
498                DOIT=.FALSE.
499             ELSE
500                I=I+1
501             ENDIF
502
503          ELSE
504             CURRENT=>CURRENT%previous
505             DOIT=.FALSE.
506          ENDIF
507
508       enddo
509    endif
510
511
512
513    L%LASTPOS=I; L%LAST => Current;
514
515    IF(I/=L%LAST%POS) THEN
516       WRITE(6,*) " ERROR IN move_to_s ",I,L%LAST%POS
517       STOP 999
518    ENDIF
519
520  END SUBROUTINE move_to_s
521
522  ! tracking one steps in the body
523
524
525
526  ! MULTIPARTICLE AT THE FIBRE LEVEL
527
528  ! front patch/misaglinments/tilt
529  SUBROUTINE TRACK_FIBRE_FRONTR(C,X,K)
530    implicit none
531    logical(lp) :: doneitt=.true.
532    TYPE(FIBRE),TARGET,INTENT(INOUT):: C
533    !    TYPE(BEAM),TARGET,INTENT(INOUT):: B
534    real(dp), INTENT(INOUT) :: X(6)
535    TYPE(INTERNAL_STATE)  K
536    !    TYPE(INTERNAL_STATE), INTENT(IN) :: K
537    logical(lp) ou,patch
538    INTEGER(2) PATCHT,PATCHG,PATCHE
539    TYPE (fibre), POINTER :: CN
540    real(dp), POINTER :: P0,B0
541
542    !FRONTAL PATCH
543    !    IF(ASSOCIATED(C%PATCH)) THEN
544    PATCHT=C%PATCH%TIME ;PATCHE=C%PATCH%ENERGY ;PATCHG=C%PATCH%PATCH;
545    !    ELSE
546    !       PATCHT=0 ; PATCHE=0 ;PATCHG=0;
547    !    ENDIF
548
549    ! PUSHING BEAM
550    !
551
552
553
554    IF(PATCHE/=0.AND.PATCHE/=2) THEN
555       NULLIFY(P0);NULLIFY(B0);
556       CN=>C%PREVIOUS
557       IF(ASSOCIATED(CN)) THEN ! ASSOCIATED
558          !          IF(.NOT.CN%PATCH%ENERGY) THEN     ! No need to patch IF PATCHED BEFORE
559          IF(CN%PATCH%ENERGY==0) THEN     ! No need to patch IF PATCHED BEFORE
560             P0=>CN%MAG%P%P0C
561             B0=>CN%MAG%P%BETA0
562
563             X(2)=X(2)*P0/C%MAG%P%P0C
564             X(4)=X(4)*P0/C%MAG%P%P0C
565             IF(k%TIME.or.recirculator_cheat)THEN
566                X(5)=root(1.0_dp+2.0_dp*X(5)/B0+X(5)**2)  !X(5) = 1+DP/P0C_OLD
567                X(5)=X(5)*P0/C%MAG%P%P0C-1.0_dp !X(5) = DP/P0C_NEW
568                X(5)=(2.0_dp*X(5)+X(5)**2)/(root(1.0_dp/C%MAG%P%BETA0**2+2.0_dp*X(5)+X(5)**2)+1.0_dp/C%MAG%P%BETA0)
569             ELSE
570                X(5)=(1.0_dp+X(5))*P0/C%MAG%P%P0C-1.0_dp
571             ENDIF
572          ENDIF ! No need to patch
573       ENDIF ! ASSOCIATED
574
575    ENDIF
576
577    ! The chart frame of reference is located here implicitely
578    IF(PATCHG==1.or.PATCHG==3) THEN
579       patch=ALWAYS_EXACT_PATCHING.or.C%MAG%P%EXACT
580       CALL PATCH_FIB(C,X,k,PATCH,MY_TRUE)
581    ENDIF
582
583    IF(PATCHT/=0.AND.PATCHT/=2.AND.(K%TOTALPATH==0)) THEN
584      if(K%time) then
585       X(6)=X(6)-C%PATCH%a_T/c%beta0
586      else
587       X(6)=X(6)-C%PATCH%a_T
588      endif
589    ENDIF
590
591    CALL DTILTD(C%DIR,C%MAG%P%TILTD,1,X)
592    ! The magnet frame of reference is located here implicitely before misalignments
593
594    !      CALL TRACK(C,X,EXACTMIS=K%EXACTMIS)
595    IF(C%MAG%MIS) THEN
596       ou = ALWAYS_EXACTMIS  !K%EXACTMIS.or.
597       CALL MIS_FIB(C,X,k,OU,DONEITT)
598    ENDIF
599
600  END SUBROUTINE TRACK_FIBRE_FRONTR
601
602  SUBROUTINE TRACK_FIBRE_FRONTP(C,X,K)
603    implicit none
604    logical(lp) :: doneitt=.true.
605    TYPE(FIBRE),TARGET,INTENT(INOUT):: C
606    !    TYPE(BEAM),TARGET,INTENT(INOUT):: B
607    TYPE(REAL_8), INTENT(INOUT) :: X(6)
608    TYPE(INTERNAL_STATE)  K
609    !    TYPE(INTERNAL_STATE), INTENT(IN) :: K
610    logical(lp) ou,patch
611    INTEGER(2) PATCHT,PATCHG,PATCHE
612    TYPE (fibre), POINTER :: CN
613    real(dp), POINTER :: P0,B0
614
615
616    !FRONTAL PATCH
617    !    IF(ASSOCIATED(C%PATCH)) THEN
618    PATCHT=C%PATCH%TIME ;PATCHE=C%PATCH%ENERGY ;PATCHG=C%PATCH%PATCH;
619    !    ELSE
620    !       PATCHT=0 ; PATCHE=0 ;PATCHG=0;
621    !    ENDIF
622
623    ! PUSHING BEAM
624    !
625
626
627
628    IF(PATCHE/=0.AND.PATCHE/=2) THEN
629       NULLIFY(P0);NULLIFY(B0);
630       CN=>C%PREVIOUS
631       IF(ASSOCIATED(CN)) THEN ! ASSOCIATED
632          !          IF(.NOT.CN%PATCH%ENERGY) THEN     ! No need to patch IF PATCHED BEFORE
633          IF(CN%PATCH%ENERGY==0) THEN     ! No need to patch IF PATCHED BEFORE
634             P0=>CN%MAGP%P%P0C
635             B0=>CN%MAGP%P%BETA0
636
637             X(2)=X(2)*P0/C%MAGP%P%P0C
638             X(4)=X(4)*P0/C%MAGP%P%P0C
639             IF(k%TIME.or.recirculator_cheat)THEN
640                X(5)=SQRT(1.0_dp+2.0_dp*X(5)/B0+X(5)**2)  !X(5) = 1+DP/P0C_OLD
641                X(5)=X(5)*P0/C%MAGP%P%P0C-1.0_dp !X(5) = DP/P0C_NEW
642                X(5)=(2.0_dp*X(5)+X(5)**2)/(SQRT(1.0_dp/C%MAGP%P%BETA0**2+2.0_dp*X(5)+X(5)**2)+1.0_dp/C%MAGP%P%BETA0)
643             ELSE
644                X(5)=(1.0_dp+X(5))*P0/C%MAGP%P%P0C-1.0_dp
645             ENDIF
646          ENDIF ! No need to patch
647       ENDIF ! ASSOCIATED
648
649    ENDIF
650
651    ! The chart frame of reference is located here implicitely
652    IF(PATCHG==1.or.PATCHG==3) THEN
653       patch=ALWAYS_EXACT_PATCHING.or.C%MAGP%P%EXACT
654       CALL PATCH_FIB(C,X,k,PATCH,MY_TRUE)
655    ENDIF
656
657    IF(PATCHT/=0.AND.PATCHT/=2.AND.(K%TOTALPATH==0)) THEN
658      if(K%time) then
659       X(6)=X(6)-C%PATCH%a_T/c%beta0
660      else
661       X(6)=X(6)-C%PATCH%a_T
662      endif
663    ENDIF
664
665    CALL DTILTD(C%DIR,C%MAGP%P%TILTD,1,X)
666    ! The magnet frame of reference is located here implicitely before misalignments
667
668    !      CALL TRACK(C,X,EXACTMIS=K%EXACTMIS)
669    IF(C%MAGP%MIS) THEN
670       ou = ALWAYS_EXACTMIS   !K%EXACTMIS.or.
671       CALL MIS_FIB(C,X,k,OU,DONEITT)
672    ENDIF
673
674
675  END SUBROUTINE TRACK_FIBRE_FRONTP
676
677
678  ! back patch/misaglinments/tilt
679
680  SUBROUTINE TRACK_FIBRE_BACKR(C,X,K)
681    implicit none
682    logical(lp) :: doneitf=.false.
683    TYPE(FIBRE),TARGET,INTENT(INOUT):: C
684    !    TYPE(BEAM),TARGET,INTENT(INOUT):: B
685    real(dp), INTENT(INOUT) :: X(6)
686    TYPE(INTERNAL_STATE)  K
687    !    TYPE(INTERNAL_STATE), INTENT(IN) :: K
688    logical(lp) ou,patch
689    INTEGER(2) PATCHT,PATCHG,PATCHE
690    TYPE (fibre), POINTER :: CN
691    real(dp), POINTER :: P0,B0
692
693
694    IF(ASSOCIATED(C%PATCH)) THEN
695       PATCHT=C%PATCH%TIME ;PATCHE=C%PATCH%ENERGY ;PATCHG=C%PATCH%PATCH;
696    ELSE
697       PATCHT=0 ; PATCHE=0 ;PATCHG=0;
698    ENDIF
699
700
701
702    IF(C%MAG%MIS) THEN
703       ou = ALWAYS_EXACTMIS  !K%EXACTMIS.or.
704       CALL MIS_FIB(C,X,k,OU,DONEITF)
705    ENDIF
706    ! The magnet frame of reference is located here implicitely before misalignments
707    CALL DTILTD(C%DIR,C%MAG%P%TILTD,2,X)
708
709    IF(PATCHT/=0.AND.PATCHT/=1.AND.(K%TOTALPATH==0)) THEN
710      if(K%time) then
711       X(6)=X(6)-C%PATCH%b_T/c%beta0
712      else
713       X(6)=X(6)-C%PATCH%b_T
714      endif
715    ENDIF
716
717    IF(PATCHG==2.or.PATCHG==3) THEN
718       patch=ALWAYS_EXACT_PATCHING.or.C%MAG%P%EXACT
719       CALL PATCH_FIB(C,X,k,PATCH,MY_FALSE)
720    ENDIF
721
722    ! The CHART frame of reference is located here implicitely
723
724    IF(PATCHE/=0.AND.PATCHE/=1) THEN
725       NULLIFY(P0);NULLIFY(B0);
726       CN=>C%NEXT
727       IF(.NOT.ASSOCIATED(CN)) CN=>C
728       !       P0=>CN%MAG%P%P0C
729       !       B0=>CN%MAG%P%BETA0
730       P0=>CN%MAG%P%P0C
731       B0=>CN%BETA0
732       X(2)=X(2)*C%MAG%P%P0C/P0
733       X(4)=X(4)*C%MAG%P%P0C/P0
734       IF(k%TIME.or.recirculator_cheat)THEN
735          X(5)=root(1.0_dp+2.0_dp*X(5)/C%MAG%P%BETA0+X(5)**2)  !X(5) = 1+DP/P0C_OLD
736          X(5)=X(5)*C%MAG%P%P0C/P0-1.0_dp !X(5) = DP/P0C_NEW
737          X(5)=(2.0_dp*X(5)+X(5)**2)/(root(1.0_dp/B0**2+2.0_dp*X(5)+X(5)**2)+1.0_dp/B0)
738       ELSE
739          X(5)=(1.0_dp+X(5))*C%MAG%P%P0C/P0-1.0_dp
740       ENDIF
741    ENDIF
742
743
744
745
746  END SUBROUTINE TRACK_FIBRE_BACKR
747
748  SUBROUTINE TRACK_FIBRE_BACKP(C,X,K)
749    implicit none
750    logical(lp) :: doneitf=.false.
751    TYPE(FIBRE),TARGET,INTENT(INOUT):: C
752    !    TYPE(BEAM),TARGET,INTENT(INOUT):: B
753    type(real_8), INTENT(INOUT) :: X(6)
754    TYPE(INTERNAL_STATE)  K
755    !    TYPE(INTERNAL_STATE), INTENT(IN) :: K
756    logical(lp) ou,patch
757    INTEGER(2) PATCHT,PATCHG,PATCHE
758    TYPE (fibre), POINTER :: CN
759    real(dp), POINTER :: P0,B0
760
761
762    IF(ASSOCIATED(C%PATCH)) THEN
763       PATCHT=C%PATCH%TIME ;PATCHE=C%PATCH%ENERGY ;PATCHG=C%PATCH%PATCH;
764    ELSE
765       PATCHT=0 ; PATCHE=0 ;PATCHG=0;
766    ENDIF
767
768
769
770    IF(C%MAGP%MIS) THEN
771       ou = ALWAYS_EXACTMIS   !K%EXACTMIS.or.
772       CALL MIS_FIB(C,X,k,OU,DONEITF)
773    ENDIF
774    ! The magnet frame of reference is located here implicitely before misalignments
775    CALL DTILTD(C%DIR,C%MAGP%P%TILTD,2,X)
776
777    IF(PATCHT/=0.AND.PATCHT/=1.AND.(K%TOTALPATH==0)) THEN
778      if(K%time) then
779       X(6)=X(6)-C%PATCH%b_T/c%beta0
780      else
781       X(6)=X(6)-C%PATCH%b_T
782      endif
783    ENDIF
784
785    IF(PATCHG==2.or.PATCHG==3) THEN
786       patch=ALWAYS_EXACT_PATCHING.or.C%MAGP%P%EXACT
787       CALL PATCH_FIB(C,X,k,PATCH,MY_FALSE)
788    ENDIF
789
790    ! The CHART frame of reference is located here implicitely
791
792    IF(PATCHE/=0.AND.PATCHE/=1) THEN
793       NULLIFY(P0);NULLIFY(B0);
794       CN=>C%NEXT
795       IF(.NOT.ASSOCIATED(CN)) CN=>C
796       !       P0=>CN%MAGP%P%P0C
797       !       B0=>CN%MAGP%P%BETA0
798       P0=>CN%MAGP%P%P0C
799       B0=>CN%BETA0
800       X(2)=X(2)*C%MAGP%P%P0C/P0
801       X(4)=X(4)*C%MAGP%P%P0C/P0
802       IF(k%TIME.or.recirculator_cheat)THEN
803          X(5)=sqrt(1.0_dp+2.0_dp*X(5)/C%BETA0+X(5)**2)  !X(5) = 1+DP/P0C_OLD
804          X(5)=X(5)*C%MAGP%P%P0C/P0-1.0_dp !X(5) = DP/P0C_NEW
805          X(5)=(2.0_dp*X(5)+X(5)**2)/(sqrt(1.0_dp/B0**2+2.0_dp*X(5)+X(5)**2)+1.0_dp/B0)
806       ELSE
807          X(5)=(1.0_dp+X(5))*C%MAGP%P%P0C/P0-1.0_dp
808       ENDIF
809    ENDIF
810
811  END SUBROUTINE TRACK_FIBRE_BACKP
812
813
814  ! thin lens tracking
815
816  SUBROUTINE TRACKV_NODE_SINGLE(T,V,K) !!
817    implicit none
818    TYPE(INTEGRATION_NODE),POINTER :: T
819    TYPE(INTERNAL_STATE)  K
820    REAL(DP) SC,reference_ray(6),x(6)
821    type(three_d_info),intent(INOUT) ::  v
822    TYPE(INTEGRATION_NODE),POINTER:: mag_in,mag_out
823
824    IF(.NOT.CHECK_STABLE) return
825    !       CALL RESET_APERTURE_FLAG
826    !    endif
827
828    if(abs(x(1))+abs(x(3))>absolute_aperture) then
829       messageLOST="exceed absolute_aperture in TRACKV_NODE_SINGLE"
830       lost_node=>t
831       lost_fibre=>t%parent_fibre
832       xlost=x
833       CHECK_STABLE=.false.
834    endif
835
836    x=V%X
837    reference_ray=V%reference_ray
838
839    CALL TRACK_NODE_SINGLE(T,V%X,K)
840
841    IF(.NOT.CHECK_STABLE)      V%U(1)=.TRUE.
842
843    CALL TRACK_NODE_SINGLE(T,V%reference_ray,K)
844
845    IF(.NOT.CHECK_STABLE)   V%U(2)=.TRUE.
846    IF(V%U(1).OR.V%U(2)) then
847       write(6,*) " Unstable in TRACKV_NODE_SINGLE ", V%U
848       RETURN
849    endif
850    IF(.NOT.ASSOCIATED(T%B)) THEN
851       WRITE(6,*) " NO FRAMES IN INTEGRATION NODES "
852       STOP 101
853    ENDIF
854
855    SC=1.0_dp
856    IF(v%SCALE/=0.0_dp) SC=v%SCALE
857    !      t=>B%POS(1)%NODE%previous
858
859    V%r0=t%A+(reference_ray(1)-SC*reference_ray(1))*t%ENT(1,1:3)+ SC*X(1)*t%ENT(1,1:3)
860    V%r0=v%r0+(reference_ray(3)-SC*reference_ray(3))*t%ENT(2,1:3)+ SC*X(3)*t%ENT(2,1:3)
861
862    V%r=t%B+(V%reference_ray(1)-SC*V%reference_ray(1))*t%EXI(1,1:3)+ SC*V%X(1)*t%EXI(1,1:3)
863    V%r=v%r+(V%reference_ray(3)-SC*V%reference_ray(3))*t%EXI(2,1:3)+ SC*V%X(3)*t%EXI(2,1:3)
864    mag_in=>t%parent_fibre%t1%next%next
865    mag_out=>t%parent_fibre%t2%previous%previous
866    v%a=mag_in%a
867    v%ent=mag_in%ent
868    v%b=mag_in%b
869    v%exi=mag_in%exi
870    v%o=t%B
871    v%mid=t%exi
872
873
874    IF(MAG_IN%PREVIOUS%CAS/=CASE1) STOP 201
875    IF(MAG_OUT%NEXT%CAS/=CASE2) STOP 202
876
877
878  END SUBROUTINE TRACKV_NODE_SINGLE
879
880  SUBROUTINE TRACKR_NODE_SINGLE(T,X,K) !!
881    ! This routines tracks a single thin lens
882    ! it is supposed to reproduce plain PTC
883    implicit none
884    TYPE(INTEGRATION_NODE), TARGET, INTENT(INOUT):: T
885    REAL(DP),INTENT(INOUT):: X(6)
886    TYPE(INTERNAL_STATE)  K
887    !    TYPE(INTERNAL_STATE), INTENT(IN) :: K
888    type(element),pointer :: el
889
890    IF(.NOT.CHECK_STABLE) return
891    !       CALL RESET_APERTURE_FLAG
892    !    endif
893
894    if(abs(x(1))+abs(x(3))>absolute_aperture) then   !.or.(.not.CHECK_MADX_APERTURE)) then
895       messageLOST="exceed absolute_aperture in TRACKR_NODE_SINGLE"
896       lost_node=>t
897       lost_fibre=>t%parent_fibre
898       xlost=x
899       CHECK_STABLE=.false.
900    endif
901
902
903    !  T%PARENT_FIBRE%MAG=K
904    T%PARENT_FIBRE%MAG%P%DIR=>T%PARENT_FIBRE%DIR
905    T%PARENT_FIBRE%MAG%P%beta0=>T%PARENT_FIBRE%beta0
906    T%PARENT_FIBRE%MAG%P%GAMMA0I=>T%PARENT_FIBRE%GAMMA0I
907    T%PARENT_FIBRE%MAG%P%GAMBET=>T%PARENT_FIBRE%GAMBET
908    T%PARENT_FIBRE%MAG%P%MASS=>T%PARENT_FIBRE%MASS
909    T%PARENT_FIBRE%MAG%P%CHARGE=>T%PARENT_FIBRE%CHARGE
910
911
912    !call cpu_time(ttime1)
913
914    !dt1=ttime1-ttime0+dt1
915
916
917    SELECT CASE(T%CAS)
918    CASE(CASEP1)
919       CALL TRACK_FIBRE_FRONT(T%PARENT_FIBRE,X,K)
920       if(associated(T%PARENT_FIBRE%MAG%p%aperture)) call CHECK_APERTURE(T%PARENT_FIBRE%MAG%p%aperture,X)
921    CASE(CASEP2)
922       CALL TRACK_FIBRE_BACK(T%PARENT_FIBRE,X,K)
923
924    CASE(CASE1,CASE2)
925       el=>T%PARENT_FIBRE%MAG
926       if(s_aperture_CHECK.and.associated(el%p%A).AND.CHECK_MADX_APERTURE.and.t%cas==case2) &
927            call check_S_APERTURE_out(el%p,t%POS_IN_FIBRE-2,x)
928
929       SELECT CASE(EL%KIND)
930       CASE(KIND0:KIND1,KIND3,KIND8:KIND9,KIND11:KIND15,KIND18:KIND19,kind22)
931       case(KIND2)
932          CALL TRACK_FRINGE(EL=EL%K2,X=X,k=k,J=T%CAS)
933       case(KIND4)
934          IF(T%CAS==CASE1) THEN
935             CALL ADJUST_TIME_CAV4(EL%C4,X,k,1)
936             CALL FRINGECAV(EL%C4,X,k=k,J=1)
937          ELSE
938             CALL FRINGECAV(EL%C4,X,k=k,J=2)
939             CALL ADJUST_TIME_CAV4(EL%C4,X,k,2)
940          ENDIF
941       case(KIND5)
942          CALL TRACK_FRINGE(EL5=EL%S5,X=X,k=k,J=T%CAS)
943       case(KIND6)
944          CALL TRACK_FRINGE(EL6=EL%T6,X=X,k=k,J=T%CAS)
945       case(KIND7)
946          CALL TRACK_FRINGE(EL7=EL%T7,X=X,k=k,J=T%CAS)
947       case(KIND10)
948          CALL FRINGE_teapot(EL%TP10,X,k=k,j=T%CAS)
949       case(KIND16,KIND20)
950          CALL fringe_STREX(EL%K16,X,k,T%CAS)
951       case(KIND17)
952          STOP 317
953       case(KIND21)
954          CALL FRINGE_CAV_TRAV(EL%CAV21,X=X,k=k,J=T%CAS)
955          CALL ADJUST_TIME_CAV_TRAV_OUT(EL%CAV21,X,k,T%CAS)   ! ONLY DOES SOMETHING IF J==2
956       case(KINDWIGGLER)
957          CALL ADJUST_WI(EL%WI,X,T%CAS)   ! ONLY DOES SOMETHING IF J==2
958       case(KINDPA)
959          CALL ADJUST_PANCAKE(EL%PA,X,k,T%CAS)
960       CASE DEFAULT
961          WRITE(6,*) "NOT IMPLEMENTED ",EL%KIND
962          stop 666
963       END SELECT
964
965    CASE(CASE0)
966       el=>T%PARENT_FIBRE%MAG
967       if(s_aperture_CHECK.and.associated(el%p%A).AND.CHECK_MADX_APERTURE)  &
968            call check_S_APERTURE(el%p,t%POS_IN_FIBRE-2,x)
969       if(associated(t%bb).and.dobb.and.do_beam_beam) then
970
971          if(t%bb%patch) call PATCH_BB(t%bb,X,k,EL%p%BETA0,ALWAYS_EXACT_PATCHING.or.EL%P%EXACT,my_true)
972          call BBKICK(t%bb,X)
973          if(t%bb%patch)call PATCH_BB(t%bb,X,k,EL%p%BETA0,ALWAYS_EXACT_PATCHING.or.EL%P%EXACT,my_false)
974
975       endif
976
977       SELECT CASE(EL%KIND)
978       CASE(KIND0)
979       case(KIND1)
980          CALL TRACK_SLICE(EL%D0,X,K)
981       case(KIND2)
982          CALL TRACK_SLICE(EL%K2,X,K,t%POS_IN_FIBRE-2)
983       case(KIND3)
984          CALL TRACK(EL%K3,X,K)
985       case(KIND4)
986          CALL TRACK_SLICE(EL%C4,X,K,t%POS_IN_FIBRE-2)
987       case(KIND5)
988          CALL TRACK_SLICE(EL%S5,X,K)
989       case(KIND6)
990          CALL TRACK_SLICE(EL%T6,X,K)
991       case(KIND7)
992          CALL TRACK_SLICE(EL%T7,X,K,t%POS_IN_FIBRE-2)
993       case(KIND8)
994          CALL TRACK(EL%S8,X,K)
995       case(KIND9)
996          CALL TRACK(EL%S9,X,K)
997       case(KIND10)
998          CALL TRACK_SLICE(EL%TP10,X,K,t%POS_IN_FIBRE-2)
999       case(KIND11:KIND14)
1000          CALL MONTI(EL%MON14,X,k,t%POS_IN_FIBRE-2)
1001          !          CALL TRACK_SLICE(EL%MON14,X,K)
1002       case(KIND15)
1003          call SEPTTRACK(EL%SEP15,X,k,t%POS_IN_FIBRE-2)
1004          !          CALL TRACK_SLICE(EL%SEP15,X,K)
1005       case(KIND16,KIND20)
1006          CALL TRACK_SLICE(EL%K16,X,K,t%POS_IN_FIBRE-2)
1007       case(KIND17)
1008          STOP 317
1009       case(KIND18)
1010          call RCOLLIMATORI(EL%RCOL18,X,k,t%POS_IN_FIBRE-2)
1011       case(KIND19)
1012          CALL ECOLLIMATORI(EL%ECOL19,X,k,t%POS_IN_FIBRE-2)
1013          !          CALL TRACK_SLICE(EL%ECOL19,X,K)
1014       case(KIND21)
1015          CALL TRACK_SLICE(EL%CAV21,X,k,t%POS_IN_FIBRE-2)
1016       case(KIND22)
1017          CALL TRACK_SLICE(EL%he22,X,k,t%POS_IN_FIBRE-2)
1018       case(KINDWIGGLER)
1019          CALL TRACK_SLICE(EL%WI,X,k,t%POS_IN_FIBRE-2)
1020       case(KINDPA)
1021          CALL TRACK_SLICE(EL%PA,X,k,T%POS_IN_FIBRE-2)
1022
1023       CASE DEFAULT
1024          WRITE(6,*) "NOT IMPLEMENTED ",EL%KIND
1025          stop 999
1026       END SELECT
1027       if(associated(T%PARENT_FIBRE%MAG%p%aperture).and.aperture_all_case0) &
1028            call CHECK_APERTURE(T%PARENT_FIBRE%MAG%p%aperture,X)
1029
1030
1031    case(CASET)
1032       if(associated(t%bb).and.dobb.and.do_beam_beam) then
1033
1034          if(t%bb%patch) call PATCH_BB(t%bb,X,k,EL%p%BETA0,ALWAYS_EXACT_PATCHING.or.EL%P%EXACT,my_true)
1035          call BBKICK(t%bb,X)
1036          if(t%bb%patch)call PATCH_BB(t%bb,X,k,EL%p%BETA0,ALWAYS_EXACT_PATCHING.or.EL%P%EXACT,my_false)
1037       endif
1038       IF(ASSOCIATED(T%T)) CALL TRACK(T%T,X)
1039    case(CASETF1,CASETF2)
1040
1041       IF(ASSOCIATED(T%T)) CALL TRACK(T%T,X)
1042
1043
1044    END SELECT
1045    ! CASE(CASE100)  ! FAKE BEAM BEAM CAKE AT SOME S
1046
1047
1048    !    T%PARENT_FIBRE%MAG=DEFAULT
1049    if(wherelost==2.and.(.not.check_stable)) then
1050       t%lost=t%lost+1
1051    endif
1052  END SUBROUTINE TRACKR_NODE_SINGLE
1053
1054
1055  SUBROUTINE TRACKP_NODE_SINGLE(T,X,K) !!
1056    ! This routines tracks a single thin lens
1057    ! it is supposed to reproduce plain PTC
1058    implicit none
1059    TYPE(INTEGRATION_NODE), TARGET, INTENT(INOUT):: T
1060    TYPE(REAL_8),INTENT(INOUT):: X(6)
1061    TYPE(INTERNAL_STATE)  K
1062    !    TYPE(INTERNAL_STATE), INTENT(IN) :: K
1063    type(elementp),pointer :: el
1064    logical(lp) BN2,L
1065    logical(lp) CHECK_KNOB
1066    integer(2), pointer,dimension(:)::AN,BN
1067
1068    IF(.NOT.CHECK_STABLE) return
1069    !       CALL RESET_APERTURE_FLAG
1070    !    endif
1071
1072    if(abs(x(1))+abs(x(3))>absolute_aperture) then
1073       messageLOST="exceed absolute_aperture in TRACKP_NODE_SINGLE"
1074       lost_node=>t
1075       lost_fibre=>t%parent_fibre
1076       xlost=x
1077       CHECK_STABLE=.false.
1078    endif
1079
1080    !   T%PARENT_FIBRE%MAGP=K
1081    IF(K%PARA_IN ) KNOB=.TRUE.
1082
1083    T%PARENT_FIBRE%MAGP%P%DIR=>T%PARENT_FIBRE%DIR
1084    T%PARENT_FIBRE%MAGP%P%beta0=>T%PARENT_FIBRE%beta0
1085    T%PARENT_FIBRE%MAGP%P%GAMMA0I=>T%PARENT_FIBRE%GAMMA0I
1086    T%PARENT_FIBRE%MAGP%P%GAMBET=>T%PARENT_FIBRE%GAMBET
1087    T%PARENT_FIBRE%MAGP%P%MASS=>T%PARENT_FIBRE%MASS
1088    T%PARENT_FIBRE%MAGP%P%CHARGE=>T%PARENT_FIBRE%CHARGE
1089
1090
1091
1092    SELECT CASE(T%CAS)
1093    CASE(CASEP1)
1094       CALL TRACK_FIBRE_FRONT(T%PARENT_FIBRE,X,K)
1095       if(associated(T%PARENT_FIBRE%MAGP%p%aperture)) call CHECK_APERTURE(T%PARENT_FIBRE%MAGP%p%aperture,X)
1096    CASE(CASEP2)
1097       !    if(abs(x(1))+abs(x(3))>absolute_aperture.or.(.not.CHECK_MADX_APERTURE)) then ! new 2010
1098       !       CHECK_STABLE=.false.
1099       !    endif
1100       CALL TRACK_FIBRE_BACK(T%PARENT_FIBRE,X,K)
1101
1102    CASE(CASE1,CASE2)
1103       el=>T%PARENT_FIBRE%MAGP
1104       if(s_aperture_CHECK.and.associated(el%p%A).AND.CHECK_MADX_APERTURE.and.t%cas==case2) &
1105            call check_S_APERTURE_out(el%p,t%POS_IN_FIBRE-2,x)
1106
1107
1108       SELECT CASE(EL%KIND)
1109       CASE(KIND0:KIND1,KIND3,KIND8:KIND9,KIND11:KIND15,KIND18:KIND19,kind22)
1110       case(KIND2)
1111          CALL TRACK_FRINGE(EL=EL%K2,X=X,k=k,J=T%CAS)
1112       case(KIND4)
1113          IF(T%CAS==CASE1) THEN
1114             CALL ADJUST_TIME_CAV4(EL%C4,X,k,1)
1115             CALL FRINGECAV(EL%C4,X,k=k,J=1)
1116          ELSE
1117             CALL FRINGECAV(EL%C4,X,k=k,J=2)
1118             CALL ADJUST_TIME_CAV4(EL%C4,X,k,2)
1119          ENDIF
1120       case(KIND5)
1121          CALL TRACK_FRINGE(EL5=EL%S5,X=X,k=k,J=T%CAS)
1122       case(KIND6)
1123          CALL TRACK_FRINGE(EL6=EL%T6,X=X,k=k,J=T%CAS)
1124       case(KIND7)
1125          CALL TRACK_FRINGE(EL7=EL%T7,X=X,k=k,J=T%CAS)
1126       case(KIND10)
1127          CALL FRINGE_teapot(EL%TP10,X,k,T%CAS)
1128       case(KIND16,KIND20)
1129          CALL fringe_STREX(EL%K16,X,k,T%CAS)
1130       case(KIND17)
1131          STOP 317
1132       case(KIND21)
1133          CALL FRINGE_CAV_TRAV(EL%CAV21,X=X,k=k,J=T%CAS)
1134          CALL ADJUST_TIME_CAV_TRAV_OUT(EL%CAV21,X,k,T%CAS)   ! ONLY DOES SOMETHING IF J==2
1135       case(KINDWIGGLER)
1136          CALL ADJUST_WI(EL%WI,X,T%CAS)   ! ONLY DOES SOMETHING IF J==2
1137       case(KINDPA)
1138          CALL ADJUST_PANCAKE(EL%PA,X,k,T%CAS)   ! ONLY DOES SOMETHING IF J==2
1139       CASE DEFAULT
1140          WRITE(6,*) "NOT IMPLEMENTED ",EL%KIND
1141          stop 666
1142       END SELECT
1143
1144    CASE(CASE0)
1145       el=>T%PARENT_FIBRE%MAGP
1146       if(s_aperture_CHECK.and.associated(el%p%A).AND.CHECK_MADX_APERTURE) &
1147            call check_S_APERTURE(el%p,t%POS_IN_FIBRE-2,x)
1148       if(associated(t%bb).and.dobb.and.do_beam_beam) then
1149
1150          if(t%bb%patch) call PATCH_BB(t%bb,X,k,EL%p%BETA0,ALWAYS_EXACT_PATCHING.or.EL%P%EXACT,my_true)
1151          call BBKICK(t%bb,X)
1152          if(t%bb%patch)call PATCH_BB(t%bb,X,k,EL%p%BETA0,ALWAYS_EXACT_PATCHING.or.EL%P%EXACT,my_false)
1153
1154       endif
1155       SELECT CASE(EL%KIND)
1156       CASE(KIND0)
1157       case(KIND1)
1158          CALL TRACK_SLICE(EL%D0,X,K)
1159       case(KIND2)
1160          CALL TRACK_SLICE(EL%K2,X,K,t%POS_IN_FIBRE-2)
1161       case(KIND3)
1162          CALL TRACK(EL%K3,X,K)
1163       case(KIND4)
1164          CALL TRACK_SLICE(EL%C4,X,K,t%POS_IN_FIBRE-2)
1165       case(KIND5)
1166          CALL TRACK_SLICE(EL%S5,X,K)
1167       case(KIND6)
1168          CALL TRACK_SLICE(EL%T6,X,K)
1169       case(KIND7)
1170          IF((EL%T7%BN(2)%KIND==3.OR.EL%T7%L%KIND==3).AND.KNOB) THEN
1171             CALL GETMAT7(EL%T7)                                      ! RECOMPUTES ONLY IF KNOB (SPEED)
1172          ENDIF
1173          CALL TRACK_SLICE(EL%T7,X,K,t%POS_IN_FIBRE-2)
1174          IF(KNOB) THEN
1175             BN2=.FALSE.
1176             L=.FALSE.
1177             IF(EL%T7%BN(2)%KIND==3) THEN
1178                BN2=.TRUE.
1179             ENDIF
1180             IF(EL%T7%L%KIND==3) THEN
1181                L=.TRUE.
1182             ENDIF
1183             IF(BN2.OR.L) THEN
1184                EL%T7%BN(2)%KIND=1
1185                EL%T7%L%KIND=1
1186                CALL KILL(EL%T7)                               ! RECOMPUTES ONLY IF KNOB (SPEED)
1187                CALL ALLOC(EL%T7)                               ! KNOB IS REMOVED THE SLOW WAY(SPEED)
1188                CALL GETMAT7(EL%T7)
1189                IF(BN2) EL%T7%BN(2)%KIND=3
1190                IF(L)  EL%T7%L%KIND=3
1191             ENDIF
1192          ENDIF
1193       case(KIND8)
1194          CALL TRACK(EL%S8,X,K)
1195       case(KIND9)
1196          CALL TRACK(EL%S9,X,K)
1197       case(KIND10)
1198          CALL MAKEPOTKNOB(EL%TP10,CHECK_KNOB,AN,BN)
1199          CALL TRACK_SLICE(EL%TP10,X,K,t%POS_IN_FIBRE-2)
1200          CALL UNMAKEPOTKNOB(EL%TP10,CHECK_KNOB,AN,BN)
1201       case(KIND11:KIND14)
1202          CALL MONTI(EL%MON14,X,k,t%POS_IN_FIBRE-2)
1203          !          CALL TRACK_SLICE(EL%MON14,X,K)
1204       case(KIND15)
1205          call SEPTTRACK(EL%SEP15,X,k,t%POS_IN_FIBRE-2)
1206          !          CALL TRACK_SLICE(EL%SEP15,X,K)
1207       case(KIND16,KIND20)
1208          CALL TRACK_SLICE(EL%K16,X,K,t%POS_IN_FIBRE-2)
1209       case(KIND17)
1210          STOP 317
1211       case(KIND18)
1212          call RCOLLIMATORI(EL%RCOL18,X,k,t%POS_IN_FIBRE-2)
1213          !          CALL TRACK_SLICE(EL%RCOL18,X,K)
1214       case(KIND19)
1215          CALL ECOLLIMATORI(EL%ECOL19,X,k,t%POS_IN_FIBRE-2)
1216          !          CALL TRACK_SLICE(EL%ECOL19,X,K)
1217       case(KIND21)
1218          CALL TRACK_SLICE(EL%CAV21,X,k,t%POS_IN_FIBRE-2)
1219       case(KINDWIGGLER)
1220          CALL TRACK_SLICE(EL%WI,X,k,t%POS_IN_FIBRE-2)
1221       case(KIND22)
1222          CALL TRACK_SLICE(EL%he22,X,k,t%POS_IN_FIBRE-2)
1223       case(KINDPA)
1224          CALL TRACK_SLICE(EL%PA,X,k,T%POS_IN_FIBRE-2)
1225       CASE DEFAULT
1226          WRITE(6,*) "NOT IMPLEMENTED ",EL%KIND
1227          stop 999
1228       END SELECT
1229       if(associated(T%PARENT_FIBRE%MAG%p%aperture).and.aperture_all_case0) &
1230            call CHECK_APERTURE(T%PARENT_FIBRE%MAG%p%aperture,X)
1231
1232    case(CASET)
1233       if(associated(t%bb).and.dobb.and.do_beam_beam) then
1234
1235          if(t%bb%patch) call PATCH_BB(t%bb,X,k,EL%p%BETA0,ALWAYS_EXACT_PATCHING.or.EL%P%EXACT,my_true)
1236          call BBKICK(t%bb,X)
1237          if(t%bb%patch)call PATCH_BB(t%bb,X,k,EL%p%BETA0,ALWAYS_EXACT_PATCHING.or.EL%P%EXACT,my_false)
1238       endif
1239       IF(ASSOCIATED(T%T)) CALL TRACK(T%T,X)
1240    case(CASETF1,CASETF2)
1241
1242       IF(ASSOCIATED(T%T)) CALL TRACK(T%T,X)
1243
1244
1245
1246    END SELECT
1247    ! CASE(CASE100)  ! FAKE BEAM BEAM CAKE AT SOME S
1248    !    T%PARENT_FIBRE%MAGP=DEFAULT
1249    ! KNOB IS RETURNED TO THE PTC DEFAULT
1250    ! NEW STUFF WITH KIND=3
1251    KNOB=.FALSE.
1252    ! END NEW STUFF WITH KIND=3
1253
1254  END SUBROUTINE TRACKP_NODE_SINGLE
1255
1256
1257
1258
1259
1260
1261
1262
1263  !  STUFF ABOUT LIST AND STRUCTURES
1264
1265  SUBROUTINE MAKE_NODE_LAYOUT( R) !
1266    ! Creates the thin layout and puts in R%T by calling MAKE_NODE_LAYOUT_2
1267    ! At this point large patches would be problematic; i.e. |d(3)|>>>0 .
1268    implicit none
1269    TYPE (LAYOUT), TARGET :: R
1270
1271    call MAKE_NODE_LAYOUT_2( R,R%T )
1272
1273  end SUBROUTINE MAKE_NODE_LAYOUT
1274
1275  SUBROUTINE sum_ds_ac( R,ds_ac_tot) !
1276    ! computes the total length for AC-modulation
1277    implicit none
1278    TYPE (LAYOUT), TARGET :: R
1279    TYPE (NODE_LAYOUT), pointer :: L
1280    INTEGER I
1281    TYPE(INTEGRATION_NODE), POINTER :: T
1282    REAL(DP), intent(inout) :: ds_ac_tot
1283
1284    t=>r%t%start
1285    ds_ac_tot=0.0_dp
1286    do i=1,r%t%n
1287
1288       ds_ac_tot=ds_ac_tot+t%ds_ac
1289
1290       t=>t%next
1291    enddo
1292  end SUBROUTINE sum_ds_ac
1293
1294  SUBROUTINE MAKE_NODE_LAYOUT_2( R,L ) !
1295    ! Creates a thin layout.
1296    ! At this point large patches would be problematic; i.e. |d(3)|>>>0 .
1297
1298    implicit none
1299    TYPE (LAYOUT), TARGET :: R
1300    TYPE (NODE_LAYOUT), pointer :: L
1301    TYPE(FIBRE), POINTER :: P
1302    INTEGER I,J,k,TEAPOT_LIKE
1303    REAL(DP) S,DLD,DL,LI,SL
1304    LOGICAL(LP) CIRCULAR
1305    TYPE(INTEGRATION_NODE), POINTER :: T1,T2,TM
1306
1307    CASE_NAME(CASEP1)="THE ENTRANCE PATCH"
1308    CASE_NAME(CASEP2)="THE EXIT PATCH  "
1309    CASE_NAME(CASE1)= "THE ENTRANCE FRINGE"
1310    CASE_NAME(CASE2)="THE EXIT FRINGE"
1311    CASE_NAME(CASE0)="A STEP OF INTEGRATION BODY"
1312    !    CASE_NAME(CASE3)="BEAM BEAM PANCAKE"
1313
1314    if(associated(L)) then
1315       CALL kill_NODE_LAYOUT(L)
1316       NULLIFY(L);
1317    endif
1318
1319    allocate(L)
1320    CALL Set_Up_NODE_LAYOUT( L )
1321    S=0.0_dp
1322    SL=0.0_dp  !  INTEGRATION LENGTH
1323    P=>R%START
1324    k=1
1325    DO I=1,R%N
1326
1327       TEAPOT_LIKE=0
1328       IF(P%MAG%P%B0/=0.0_dp) TEAPOT_LIKE=1
1329       IF(P%MAG%KIND==KIND16.OR.P%MAG%KIND==KIND16)TEAPOT_LIKE=0
1330       IF(P%MAG%KIND==KIND0.AND.P%MAG%P%NST/=1) THEN
1331          WRITE(6,*) "MARKER SHOULD HAVE NST=1 OTHERWISE PROBLEMS "
1332          WRITE(6,*) "WILL OCCUR WITH THE WORM AND THE NODE_LAYOUT SURVEY "
1333          STOP 500
1334       ENDIF
1335       IF(P%DIR==1) THEN
1336          LI=0.0_dp;
1337       ELSE
1338          LI=P%MAG%L;
1339       ENDIF
1340       DL=P%DIR*P%MAG%L/P%MAG%P%NST
1341       DLD=P%MAG%P%LD/P%MAG%P%NST
1342       CALL APPEND_EMPTY_THIN( L )
1343       L%END%TEAPOT_LIKE=TEAPOT_LIKE
1344       L%END%S(1)=S;L%END%S(2)=LI;L%END%S(3)=SL;L%END%S(4)=0.0_dp;    ! s(1) total ld
1345       L%END%S(5)=0.0_dp;L%END%DS_AC=0.0_dp;
1346       T1=>L%END                            ! s(2) local integration distance
1347       ! s(3) total integration distance
1348       L%END%CAS=CASEP1                       ! s(4) end of step =  DL
1349       L%END%pos_in_fibre=1
1350       L%END%pos=k;k=k+1;
1351       L%END%PARENT_NODE_LAYOUT=>L
1352       L%END%PARENT_FIBRE=>P
1353
1354       CALL APPEND_EMPTY_THIN( L )
1355       L%END%TEAPOT_LIKE=TEAPOT_LIKE
1356       L%END%S(1)=S;L%END%S(2)=LI;L%END%S(3)=SL;L%END%S(4)=0.0_dp;L%END%S(5)=0.0_dp;L%END%DS_AC=0.0_dp;
1357       L%END%CAS=CASE1
1358       L%END%pos_in_fibre=2
1359       L%END%pos=k;k=k+1;
1360       L%END%PARENT_NODE_LAYOUT=>L
1361       L%END%PARENT_FIBRE=>P
1362
1363       DO J=1,P%MAG%P%NST
1364          CALL APPEND_EMPTY_THIN( L )
1365          L%END%TEAPOT_LIKE=TEAPOT_LIKE
1366          L%END%S(1)=S;L%END%S(2)=LI;L%END%S(3)=SL;L%END%S(4)=DL;L%END%S(5)=DLD;L%END%ds_ac=DLD;
1367          L%END%CAS=CASE0
1368          L%END%pos_in_fibre=J+2
1369          L%END%pos=k;k=k+1;
1370          L%END%PARENT_NODE_LAYOUT=>L
1371          L%END%PARENT_FIBRE=>P
1372          S=S+DLD
1373          LI=LI+DL
1374          SL=SL+P%DIR*DL
1375          IF(MOD(P%MAG%P%NST,2)==0) THEN
1376             IF(J==P%MAG%P%NST/2) TM=>L%END    !+1
1377          ELSE
1378             IF(J==1) TM=>T1
1379          ENDIF
1380       ENDDO
1381
1382       CALL APPEND_EMPTY_THIN( L )
1383       L%END%TEAPOT_LIKE=TEAPOT_LIKE
1384       L%END%S(1)=S;L%END%S(2)=LI;L%END%S(3)=SL;L%END%S(4)=0.0_dp;L%END%S(5)=0.0_dp;L%END%DS_AC=0.0_dp;
1385       L%END%CAS=CASE2
1386       L%END%pos_in_fibre=P%MAG%P%NST+3
1387       L%END%pos=k;k=k+1;
1388       L%END%PARENT_FIBRE=>P
1389
1390       CALL APPEND_EMPTY_THIN( L )
1391       L%END%TEAPOT_LIKE=TEAPOT_LIKE
1392       L%END%S(1)=S;L%END%S(2)=LI;L%END%S(3)=SL;L%END%S(4)=0.0_dp;L%END%S(5)=0.0_dp;L%END%DS_AC=0.0_dp;
1393       L%END%CAS=CASEP2
1394       L%END%pos_in_fibre=P%MAG%P%NST+4
1395       L%END%pos=k;k=k+1;
1396       L%END%PARENT_NODE_LAYOUT=>L
1397       L%END%PARENT_FIBRE=>P
1398       T2=>L%END
1399
1400       P%T1=>T1
1401       P%T2=>T2
1402       P%TM=>TM
1403
1404       P=>P%NEXT
1405    ENDDO
1406    L%N=k-1
1407
1408    L%PARENT_LAYOUT=>R
1409
1410    IF(R%CLOSED) THEN
1411       l%closed=.true.
1412       CIRCULAR=.TRUE.
1413       CALL RING_L_THIN(L,CIRCULAR)
1414    ENDIF
1415
1416   if(lielib_print(12)==1)  call stat_NODE_LAYOUT(l)
1417
1418  END SUBROUTINE MAKE_NODE_LAYOUT_2
1419
1420
1421  SUBROUTINE stat_NODE_LAYOUT( L )
1422    implicit none
1423    TYPE (NODE_LAYOUT), pointer :: L
1424
1425    WRITE(6,*)  " PARENT LAYOUT NAME :", L%PARENT_LAYOUT%NAME(1:len_trim(L%PARENT_LAYOUT%NAME))
1426    WRITE(6,*) " NUMBER OF ORIGINAL LAYOUT ELEMENTS :", L%PARENT_LAYOUT%N
1427    WRITE(6,*) " NUMBER OF THIN OBJECTS :", L%N
1428    WRITE(6,*) " TOTAL IDEAL LENGTH OF STRUCTURE :", L%END%S(1)
1429    WRITE(6,*) " TOTAL INTEGRATION LENGTH OF STRUCTURE (mad8 style survey) :", L%END%S(3)
1430
1431  end SUBROUTINE stat_NODE_LAYOUT
1432
1433
1434  SUBROUTINE DRIFT_TO_TIME(T,YL,DT,X,k)
1435    ! Drifts to a given time using either a regular or TEAPOT drift. (Cylindrical coordinate drift)
1436    ! The amount of the drift YL is computed to achieve a time DT.
1437
1438    IMPLICIT NONE
1439    real(dp),INTENT(INOUT):: X(6)
1440    real(dp),INTENT(INOUT):: YL
1441    real(dp),INTENT(IN):: DT
1442    TYPE(INTEGRATION_NODE), pointer :: T
1443    TYPE(magnet_chart), pointer :: p
1444    real(dp) PZ
1445    real(dp)  b
1446    TYPE(INTERNAL_STATE) k !,OPTIONAL :: K
1447
1448    P=>T%PARENT_FIBRE%MAG%P
1449
1450    IF(k%TIME) then
1451       B=P%BETA0
1452    ELSE
1453       B=1.0_dp
1454    ENDIF
1455
1456
1457    !   if(P%B0/=zero.AND.T%TEAPOT_LIKE==1) then
1458    !
1459
1460    !       PZ=ROOT(one+two*x(5)/b+X(5)**2-X(2)**2-X(4)**2)
1461    !       R=one/P%B0
1462    !       A=(X(1)+R)*(one/b+x(5))/PZ
1463    !       A=ATAN(DT/(A+DT*X(2)/PZ) )
1464    !       YL=A/P%B0
1465    !       PT=one-X(2)*TAN(A)/PZ
1466    !       XN(1)=(X(1)+R)/COS(A)/PT-R
1467    !       XN(2)=X(2)*COS(A)+SIN(A)*PZ
1468    !       XN(3)=X(3)+X(4)*(X(1)+R)*TAN(A)/PZ/PT
1469    !       XN(6)=X(6)+(X(1)+R)*TAN(A)/PZ/PT*(one/b+x(5))
1470    !       X(1)=XN(1)
1471    !       X(2)=XN(2)
1472    !       X(3)=XN(3)
1473    !       X(6)=XN(6)
1474    !    else
1475
1476
1477    PZ=ROOT(1.0_dp+2.0_dp*X(5)/b+x(5)**2-X(2)**2-X(4)**2)
1478
1479    YL=DT/(1.0_dp/b+X(5))*PZ
1480
1481    X(1)=X(1)+YL*X(2)/PZ
1482    X(3)=X(3)+YL*X(4)/PZ
1483    X(6)=X(6)+YL*(1.0_dp/b+X(5))/PZ
1484
1485
1486
1487    !    endif
1488
1489
1490
1491  END SUBROUTINE DRIFT_TO_TIME
1492
1493  SUBROUTINE DRIFTr_BACK_TO_POSITION(T,YL,X,k)
1494    ! This is a regular drift
1495    ! It is used in time tracking to project back to the beginning of the thin lens
1496    ! and it is used in S tracking to drift in the middle of a step.
1497
1498    IMPLICIT NONE
1499    real(dp),INTENT(INOUT):: X(6)
1500    real(dp),INTENT(IN):: YL
1501    TYPE(INTEGRATION_NODE), pointer :: T
1502    !    TYPE(magnet_chart), pointer :: p
1503    !    TYPE(magnet_chart), pointer :: p
1504    real(dp) PZ
1505    real(dp)  b
1506    TYPE(INTERNAL_STATE) k !,OPTIONAL :: K
1507
1508    !    P=>T%PARENT_FIBRE%MAG%P
1509
1510    IF(k%TIME) then
1511       B=T%PARENT_FIBRE%BETA0
1512    ELSE
1513       B=1.0_dp
1514    ENDIF
1515
1516
1517    !    if(P%B0/=zero.AND.T%TEAPOT_LIKE==1) then
1518
1519
1520
1521    !      PZ=ROOT(one+two*x(5)/b+X(5)**2-X(2)**2-X(4)**2)
1522    !      R=one/P%B0
1523    !      A=-YL*P%B0
1524
1525    !       PT=one-X(2)*TAN(A)/PZ
1526
1527    !       XN(1)=(X(1)+R*(two*sin(a/two)**2+X(2)*sin(A)/PZ))/COS(A)/PT
1528    !       XN(2)=X(2)*COS(A)+SIN(A)*PZ
1529    !       XN(3)=X(3)+X(4)*(X(1)+R)*TAN(A)/PZ/PT
1530    !       XN(6)=X(6)+(X(1)+R)*TAN(A)/PZ/PT*(one/b+x(5))
1531    !          WRITE(6,*) "XN(6)-X(6)-DT , DT",DT,XN(6)-X(6)-DT
1532    !      X(1)=XN(1)
1533    !      X(2)=XN(2)
1534    !      X(3)=XN(3)
1535    !      X(6)=XN(6)
1536    !   else
1537    !       CALL DRIFT(YL,DL,P%beta0,1,P%EXACT,P%TIME,X)
1538    PZ=ROOT(1.0_dp+2.0_dp*X(5)/b+x(5)**2-X(2)**2-X(4)**2)
1539
1540
1541    X(1)=X(1)-YL*X(2)/PZ
1542    X(3)=X(3)-YL*X(4)/PZ
1543    X(6)=X(6)-YL*(1.0_dp/b+X(5))/PZ
1544
1545    !  endif
1546
1547
1548
1549  END SUBROUTINE DRIFTr_BACK_TO_POSITION
1550
1551
1552  SUBROUTINE DRIFTp_BACK_TO_POSITION(T,YL,X,k)
1553    ! This is a regular drift
1554    ! It is used in time tracking to project back to the beginning of the thin lens
1555    ! and it is used in S tracking to drift in the middle of an step.
1556
1557    IMPLICIT NONE
1558    type(real_8),INTENT(INOUT):: X(6)
1559    real(dp),INTENT(IN):: YL
1560    TYPE(INTEGRATION_NODE), pointer :: T
1561    TYPE(magnet_chart), pointer :: p
1562    type(real_8) XN(6),PZ,PT
1563    real(dp)  b
1564    TYPE(INTERNAL_STATE) k !,OPTIONAL :: K
1565
1566    P=>T%PARENT_FIBRE%MAG%P
1567
1568    IF(k%TIME) then
1569       B=P%BETA0
1570    ELSE
1571       B=1.0_dp
1572    ENDIF
1573
1574    call alloc(xn,6); call alloc(pz,pt);
1575
1576    !    if(P%B0/=zero.AND.T%TEAPOT_LIKE==1) then
1577
1578
1579
1580    !       PZ=sqrt(one+two*x(5)/b+X(5)**2-X(2)**2-X(4)**2)
1581    !       R=one/P%B0
1582    !       A=-YL*P%B0
1583
1584    !      PT=one-X(2)*TAN(A)/PZ
1585
1586    !      XN(1)=(X(1)+R*(two*sin(a/two)**2+X(2)*sin(A)/PZ))/COS(A)/PT
1587    !      XN(2)=X(2)*COS(A)+SIN(A)*PZ
1588    !      XN(3)=X(3)+X(4)*(X(1)+R)*TAN(A)/PZ/PT
1589    !      XN(6)=X(6)+(X(1)+R)*TAN(A)/PZ/PT*(one/b+x(5))
1590    !          WRITE(6,*) "XN(6)-X(6)-DT , DT",DT,XN(6)-X(6)-DT
1591    !      X(1)=XN(1)
1592    !      X(2)=XN(2)
1593    !      X(3)=XN(3)
1594    !      X(6)=XN(6)
1595    !   else
1596    !       CALL DRIFT(YL,DL,P%beta0,1,P%EXACT,P%TIME,X)
1597    PZ=sqrt(1.0_dp+2.0_dp*X(5)/b+x(5)**2-X(2)**2-X(4)**2)
1598
1599
1600    X(1)=X(1)-YL*X(2)/PZ
1601    X(3)=X(3)-YL*X(4)/PZ
1602    X(6)=X(6)-YL*(1.0_dp/b+X(5))/PZ
1603
1604    !   endif
1605
1606    call kill(xn,6); call kill(pz,pt);
1607
1608
1609
1610  END SUBROUTINE DRIFTp_BACK_TO_POSITION
1611
1612
1613  !  Survey still worm like
1614
1615  SUBROUTINE FILL_SURVEY_DATA_IN_NODE_LAYOUT(R)
1616    ! THIS SUBROUTINE ALLOCATES NODE FRAMES IF NEEDED
1617    ! IT SURVEYS THE NODES USING THE OLD REAL WORMS
1618    ! SHOULD BE CALLED AFTER MISALIGNMENTS OR MOVING PART OF LATTICE
1619
1620    IMPLICIT NONE
1621    type(layout),target:: r
1622    type(fibre), pointer ::c
1623    type(INTEGRATION_NODE), pointer ::t
1624    type(worm) vers
1625    integer k,ic,j
1626    real(dp) x(6),ent(3,3),a(3)
1627    LOGICAL(LP) APER
1628    aper=APERTURE_FLAG
1629    APERTURE_FLAG=.FALSE.
1630
1631    if(.not.associated(r%t)) call MAKE_NODE_LAYOUT(r)
1632
1633    CALL  allocate_node_frame( R)
1634
1635    call survey(r)
1636
1637    CALL ALLOC(vers,r)
1638    C=>r%START
1639    CALL XFRAME(vers%E,C%chart%f%ent,C%chart%f%A,-7)  ! initializes the survey part of worm
1640    vers%E%L(-1)=0.d0 !Starts beam line at z=0   fake distance along ld for cheap work
1641
1642    do k=1,r%n
1643       x=0.0_dp
1644       CALL TRACK(r,x,k,k+1,default,vers)
1645       if(.not.check_stable) then
1646          Write(6,*) " fake instability at ",c%mag%name, " ",k
1647          check_stable=.true.
1648          CALL RESET_APERTURE_FLAG
1649       endif
1650
1651       t=>c%t1
1652       j=-6
1653       call gMID(vers,x,j)
1654       call G_FRAME(vers%e,ENT,A,j)
1655       t%ent=ent
1656       t%a=a
1657
1658       t=>t%next
1659       if(t%cas/=case1) then
1660          write(6,*)" error in fill_survey_data_in_NODE_LAYOUT",j,t%cas
1661          stop 665
1662       endif
1663       j=vers%POS(2)
1664       call gMID(vers,x,j)
1665       call G_FRAME(vers%e,ENT,A,j)
1666       t%ent=ent
1667       t%a=a
1668       t%previous%exi=ent
1669       t%previous%b=a
1670       t=>t%next
1671       ic=0
1672       DO J=vers%POS(2)+1,vers%POS(3)-1     ! pos(2)+1 to pos(3)-1 inside the magnet
1673
1674          ic=ic+1
1675
1676          call gMID(vers,x,j)
1677          call G_FRAME(vers%e,ENT,A,j)
1678
1679          if(j/=vers%POS(2)+1) then
1680             t%previous%exi=ent
1681             t%previous%b=a
1682             if(t%previous%cas/=case0) then
1683                write(6,*)" error in fill_survey_data_in_NODE_LAYOUT",j,t%previous%cas
1684                stop 666
1685             endif
1686          else
1687             t%previous%exi=ent
1688             t%previous%b=a
1689             if(t%previous%cas/=case1) then
1690                write(6,*)" error in fill_survey_data_in_NODE_LAYOUT",j,t%previous%cas
1691                stop 664
1692             endif
1693          endif
1694
1695          if(j/=vers%POS(3)-1) then
1696             t%ent=ent
1697             t%a=a
1698             if(t%cas/=case0) then
1699                write(6,*)" error in fill_survey_data_in_NODE_LAYOUT",j,t%cas
1700                stop 666
1701             endif
1702          else
1703             t%ent=ent
1704             t%a=a
1705             if(t%cas/=case2) then
1706                write(6,*)" error in fill_survey_data_in_NODE_LAYOUT",j,t%cas
1707                write(6,*)t%POS,T%PARENT_FIBRE%MAG%NAME
1708                write(6,*)T%PARENT_FIBRE%T1%POS,T%PARENT_FIBRE%T2%POS
1709                stop 668
1710             endif
1711          endif
1712
1713
1714          !  omega(1)= a(1)+scale*(xr(1)*ent(1,1)+xr(3)*ent(2,1))
1715          !  omega(2)= a(2)+scale*(xr(1)*ent(1,2)+xr(3)*ent(2,2))
1716          !  omega(3)= a(3)+scale*(xr(1)*ent(1,3)+xr(3)*ent(2,3))
1717          !  r1(1)=omega0(3)
1718          !  r1(2)=omega0(1)
1719          !  r2(1)=omega(3)
1720          !  r2(2)=omega(1)
1721          !  if(abs(r1(1))>1.d6) r1=r2
1722          !  call gMoveTo2D(r1(1),r1(2))
1723          !  call  gDrawLineTo2D(r2(1),r2(2))
1724          !  omega0=omega
1725          t=>t%next
1726       enddo
1727       j=vers%POS(3)
1728       call gMID(vers,x,j)
1729       call G_FRAME(vers%e,ENT,A,j)
1730       t%previous%exi=ent
1731       t%previous%b=a
1732       t%ent=ent
1733       t%a=a
1734       if(t%previous%cas/=case2) then
1735          write(6,*)" error in fill_survey_data_in_NODE_LAYOUT",j,t%cas
1736          stop 669
1737       endif
1738       !      t=>t%next
1739
1740       j=vers%nst
1741       call gMID(vers,x,j)
1742       call G_FRAME(vers%e,ENT,A,j)
1743
1744       t%exi=ent
1745       t%b=a
1746
1747       if(t%cas/=casep2) then
1748          write(6,*)" error in fill_survey_data_in_NODE_LAYOUT",j,t%cas
1749          stop 670
1750       endif
1751
1752
1753       if(ic/=c%mag%p%nst+1) then
1754          write(6,*)" error in fill_survey_data_in_NODE_LAYOUT"
1755          write(6,*)k, ic,c%mag%name,c%mag%p%nst
1756          stop 888
1757       endif
1758       c=>c%next
1759    enddo
1760
1761    CALL kill(vers)
1762
1763    APERTURE_FLAG=aper
1764
1765  end  subroutine fill_survey_data_in_NODE_LAYOUT
1766
1767
1768  ! BEAM STUFF
1769
1770  subroutine create_beam(B,N,CUT,SIG,A,C,T)
1771    USE gauss_dis
1772    implicit none
1773    INTEGER N,I,J,K
1774    REAL(DP) CUT,SIG(6),X
1775    TYPE(BEAM) B
1776    REAL(DP), OPTIONAL :: A(6,6),C(6)
1777    REAL(DP)  AS(6,6),CL(6),XT(6)
1778    TYPE (INTEGRATION_NODE),optional,target::  T
1779
1780    IF(.NOT.ASSOCIATED(B%N)) THEN
1781       CALL ALLOCATE_BEAM(B,N)
1782    ELSEIF(B%N/=N) THEN
1783       CALL KILL_BEAM(B)
1784       CALL ALLOCATE_BEAM(B,N)
1785    ENDIF
1786
1787    CL=0.0_dp; AS=0.0_dp;
1788    DO I=1,6
1789       AS(I,I)=1.0_dp
1790    ENDDO
1791
1792    IF(PRESENT(A)) AS=A
1793    IF(PRESENT(C)) CL=C
1794
1795    DO I=1,N
1796       DO J=1,6
1797          CALL GRNF(X,cut)
1798          XT(J)=X*SIG(J)
1799       ENDDO
1800       B%X(I,1:6)=CL(:)
1801       DO J=1,6
1802          DO K=1,6
1803             B%X(I,J)=AS(J,K)*XT(K)+B%X(I,J)
1804          ENDDO
1805
1806       ENDDO
1807    ENDDO
1808
1809
1810    if(present(t)) then
1811       DO I=1,N
1812          ! if(associated(B%POS(I)%NODE))then
1813          B%POS(I)%NODE=>T
1814          ! endif
1815       ENDDO
1816
1817    endif
1818
1819  end    subroutine create_beam
1820
1821  subroutine create_PANCAKE(B,N,CUT,SIG,T,A)
1822    USE gauss_dis
1823    implicit none
1824    INTEGER N,I,J
1825    REAL(DP) CUT,SIG(6),X,Y(LNV),beta(2)
1826    TYPE(BEAM) B
1827    TYPE (INTEGRATION_NODE),optional,target::  T
1828    TYPE (DAMAP),OPTIONAL :: A
1829    TYPE (tree) monkey
1830
1831    IF(.NOT.ASSOCIATED(B%N)) THEN
1832       CALL ALLOCATE_BEAM(B,N)
1833    ELSEIF(B%N/=N) THEN
1834       CALL KILL_BEAM(B)
1835       CALL ALLOCATE_BEAM(B,N)
1836    ENDIF
1837    write(6,*) n," particles created"
1838    Y=0.0_dp
1839    IF(.not.PRESENT(A)) THEN
1840
1841       DO I=1,N
1842          DO J=1,6
1843             CALL GRNF(X,cut)
1844             B%X(I,J)=X*SIG(J)
1845          ENDDO
1846          B%X(I,7)=0.0_dp
1847       enddo
1848    ELSE
1849       call alloc(monkey)
1850       beta(1)=(a%v(1).sub.'1')**2+(a%v(1).sub.'01')**2
1851       beta(2)=(a%v(3).sub.'001')**2+(a%v(3).sub.'0001')**2
1852       write(6,*) " Betas in create_PANCAKE ",beta
1853       monkey=A
1854       DO I=1,N
1855          DO J=1,C_%ND
1856             CALL GRNF(X,cut)
1857             Y(2*j-1)=X*sqrt(SIG(j)/2.0_dp)
1858             CALL GRNF(X,cut)
1859             Y(2*j)=X*sqrt(SIG(j)/2.0_dp)
1860          ENDDO
1861          y=monkey*Y
1862          B%X(I,1:C_%ND2)=y(1:c_%nd2)
1863
1864          DO J=C_%ND2+1,6
1865             CALL GRNF(X,cut)
1866             B%X(I,J)=X*SIG(J)
1867          ENDDO
1868          B%X(I,7)=0.0_dp
1869       enddo
1870       CALL KILL(MONKEY)
1871    ENDIF
1872
1873
1874
1875
1876    if(present(t)) then
1877       DO I=1,N
1878          if(associated(B%POS(I)%NODE))then
1879             B%POS(I)%NODE=>T
1880          endif
1881       ENDDO
1882
1883    endif
1884  end    subroutine create_PANCAKE
1885
1886  subroutine copy_beam(B1,B2)
1887    implicit none
1888    INTEGER I
1889    TYPE(BEAM), INTENT(INOUT) :: B1,B2
1890
1891    IF(.NOT.ASSOCIATED(B2%N)) THEN
1892       CALL ALLOCATE_BEAM(B2,B1%N)
1893    ELSEIF(B1%N/=B2%N) THEN
1894       CALL KILL_BEAM(B2)
1895       CALL ALLOCATE_BEAM(B2,B1%N)
1896    ENDIF
1897
1898    B2%X=B1%X
1899    B2%U=B1%U
1900    B2%N=B1%N
1901    !    B2%CHARGE=B1%CHARGE
1902    B2%LOST=B1%LOST
1903    DO I=0,B1%N
1904       if(associated(B1%POS(I)%NODE))then
1905          B2%POS(I)%NODE=>B1%POS(I)%NODE
1906       endif
1907    ENDDO
1908
1909  END subroutine copy_beam
1910
1911  subroutine READ_beam_raw(B,MF)
1912    implicit none
1913    INTEGER k,mf
1914    TYPE(BEAM), INTENT(IN):: B
1915
1916    DO K=1,b%n
1917       IF(.not.B%U(K)) THEN
1918          if(associated(b%pos(k)%NODE)) then
1919             WRITE(MF,100) B%X(K,1:6),b%pos(k)%NODE%s(3)+B%X(K,7)
1920          else
1921             WRITE(MF,100) B%X(K,1:6),B%X(K,7)
1922          endif
1923       ENDIF
1924    ENDDO
1925100 FORMAT(7(1x,e13.6))
1926  END subroutine READ_beam_raw
1927
1928  subroutine PRINT_beam_raw(B,MF)
1929    implicit none
1930    INTEGER k,mf
1931    TYPE(BEAM), INTENT(IN):: B
1932
1933    DO K=1,b%n
1934       IF(.not.B%U(K)) THEN
1935          if(associated(b%pos(k)%NODE)) then
1936             WRITE(MF,100) B%X(K,1:6),b%pos(k)%NODE%s(3)+B%X(K,7)
1937          else
1938             WRITE(MF,100) B%X(K,1:6),B%X(K,7)
1939          endif
1940       ENDIF
1941    ENDDO
1942100 FORMAT(7(1x,e13.6))
1943  END subroutine PRINT_beam_raw
1944
1945  subroutine stat_beam_raw(B,n,MF,xm)
1946    implicit none
1947    INTEGER i,j,k,mf,NOTlost,N
1948    TYPE(BEAM), INTENT(IN):: B
1949    real(dp), optional :: xm(6)
1950    real(dp), allocatable :: av(:,:)
1951    real(dp) em(2),beta(2),xma(6)
1952    allocate(av(n,n))
1953    av=0.0_dp
1954    notlost=0
1955    xma(:)=-1.0_dp
1956    DO K=1,b%n
1957       IF(.not.B%U(K)) THEN
1958          do i=1,6
1959             if(abs(b%x(k,i))>xma(i)) xma(i)=abs(b%x(k,i))
1960          enddo
1961
1962          do i=1,n
1963             do j=i,n
1964                av(i,j)= b%x(k,i)*b%x(k,j)+av(i,j)
1965             enddo
1966          enddo
1967          notlost=notlost+1
1968       ENDIF
1969    ENDDO
1970    IF(NOTLOST==0) THEN
1971       if(mf/=6) then
1972          WRITE(mf,*) " ALL PARTICLES ARE LOST "
1973          WRITE(mf,*) " NO STATISTICS "
1974       else
1975          WRITE(6,*) " ALL PARTICLES ARE LOST "
1976          WRITE(6,*) " NO STATISTICS "
1977       endif
1978       deallocate(av)
1979       RETURN
1980    ENDIF
1981    if(notlost/=b%n-b%lost) then
1982       Write(6,*) " Error keeping track of lost particles "
1983       stop 999
1984    endif
1985
1986    WRITE(MF,*) " NUMBER LEFT ",B%N-B%LOST
1987    if(mf/=6)WRITE(6,*) " NUMBER LEFT ",B%N-B%LOST
1988    WRITE(MF,*) " LOST ",B%LOST
1989    if(mf/=6)WRITE(6,*) " LOST ",B%LOST
1990    av=av/notlost
1991    em(1)=2.0_dp*sqrt(av(1,1)*av(2,2)-av(1,2)**2)
1992    em(2)=2.0_dp*sqrt(av(3,3)*av(4,4)-av(3,4)**2)
1993    beta(1)=2.0_dp*av(1,1)/em(1)
1994    beta(2)=2.0_dp*av(3,3)/em(2)
1995
1996    write(mf,*) " average arrays "
1997    write(mf,*) "betas ",beta
1998    write(mf,*) "emittances ",em
1999    if(mf/=6) then
2000       write(6,*) " average arrays "
2001       write(6,*) "betas ",beta
2002       write(6,*) "emittances ",em
2003    endif
2004    write(6,*) " limits "
2005    write(6,*) xma(1:2)
2006    write(6,*) xma(3:4)
2007    write(6,*) xma(5:6)
2008    if(present(xm)) xm=xma
2009
2010100 FORMAT(7(1x,e13.6))
2011
2012    deallocate(av)
2013  END subroutine stat_beam_raw
2014
2015
2016  subroutine PRINT_beam(B,MF,I)
2017    implicit none
2018    INTEGER K,MF,I1,I2
2019    INTEGER,OPTIONAL:: I
2020    TYPE(BEAM), INTENT(IN):: B
2021    TYPE(INTEGRATION_NODE),POINTER::T
2022    TYPE(FIBRE),POINTER::F
2023
2024    I1=1
2025    I2=B%N
2026
2027    IF(PRESENT(I)) THEN
2028       I1=I
2029       I2=I
2030    ENDIF
2031    !    IF(B%TIME_INSTEAD_OF_S) THEN
2032    !       WRITE(MF,*) "____________________________ TIME TRACKED BEAM __________________________________"
2033    !    ELSE
2034    WRITE(MF,*) "_________________ POSITION TRACKED BEAM (AS IN PTC PROPER)_______________________"
2035    !    ENDIF
2036
2037    DO K=I1,I2
2038       IF(B%U(K)) THEN
2039          WRITE(MF,*) " PARTICLE # ",K, " IS LOST "
2040       ELSE
2041          T=>B%POS(K)%NODE
2042          F=>T%PARENT_FIBRE
2043          WRITE(MF,*) "_________________________________________________________________________"
2044          WRITE(MF,*) " PARTICLE # ",K, " IS LOCATED AT SLICE # ",T%POS," IN FIBRE  ",F%MAG%NAME
2045          WRITE(MF,*) " IN THE FIBRE POSITION  ",T%pos_in_fibre
2046          WRITE(MF,*) " IN ",CASE_NAME(T%CAS)
2047          IF(T%CAS==CASE0)WRITE(MF,*) " AT THE STEP NUMBER ",T%pos_in_fibre-2
2048
2049          WRITE(MF,*) "........................................................................."
2050          !          IF(B%TIME_INSTEAD_OF_S) THEN
2051          !             WRITE(MF,*) " TIME AND POSITION AFTER THIN SLICE = ",B%X(K,6:7)
2052          !          ELSE
2053          WRITE(MF,*) " TIME AND POSITION  = ",B%X(K,6:7)
2054          !          ENDIF
2055          WRITE(MF,*) " X,Y = ",B%X(K,1),B%X(K,3)
2056          WRITE(MF,*) " PX,PY = ",B%X(K,2),B%X(K,4)
2057          WRITE(MF,*) " ENERGY VARIABLE = ",B%X(K,5)
2058       ENDIF
2059       WRITE(MF,*) "_________________________________________________________________________"
2060    ENDDO
2061
2062  END subroutine PRINT_beam
2063
2064
2065
2066  SUBROUTINE NULLIFY_BEAM(B)
2067    IMPLICIT NONE
2068    TYPE(BEAM) , INTENT (INOUT) :: B
2069    NULLIFY(B%N,B%LOST)
2070    !    NULLIFY(B%Y)
2071    NULLIFY(B%X)
2072    NULLIFY(B%U)
2073    NULLIFY(B%POS)
2074    !    NULLIFY(B%CHARGE)
2075    !    NULLIFY(B%TIME_INSTEAD_OF_S)
2076    !    NULLIFY(B%SIGMA)
2077    !    NULLIFY(B%DX,B%ORBIT)
2078    !    NULLIFY(B%BBPAR,B%BEAM_BEAM,B%BBORBIT)
2079  END SUBROUTINE NULLIFY_BEAM
2080
2081  SUBROUTINE NULLIFY_BEAMS(B)
2082    IMPLICIT NONE
2083    TYPE(BEAM) , INTENT (INOUT) :: B(:)
2084    INTEGER I
2085    DO I=1,SIZE(B)
2086       CALL NULLIFY_BEAM(B(i))
2087    ENDDO
2088
2089  END SUBROUTINE NULLIFY_BEAMS
2090
2091  subroutine alloc_three_d_info(v)
2092    IMPLICIT NONE
2093    TYPE(three_d_info) , INTENT (INOUT) :: V
2094    v%a=0.0_dp
2095    v%b=0.0_dp
2096    v%o=0.0_dp
2097    v%ent=global_frame
2098    v%exi=global_frame
2099    v%mid=global_frame
2100    v%reference_ray=0.0_dp
2101    v%r0=0.0_dp
2102    v%r=0.0_dp
2103    v%x=0.0_dp
2104    v%scale=1.0_dp
2105    v%u=my_false
2106    v%wx=0.1_dp
2107    v%wy=0.1_dp
2108
2109  end subroutine alloc_three_d_info
2110
2111  SUBROUTINE ALLOCATE_BEAM(B,N)
2112    IMPLICIT NONE
2113    TYPE(BEAM) , INTENT (INOUT) :: B
2114    INTEGER , INTENT (IN) :: N
2115    INTEGER I
2116
2117    ALLOCATE(B%N,B%LOST)
2118
2119    B%N=N
2120    B%LOST=0
2121    !    NULLIFY(B%Y)
2122    !    IF(PRESENT(POLYMORPH)) THEN
2123    !       IF(POLYMORPH) then
2124    !          ALLOCATE(B%Y(6))
2125    !          CALL ALLOC(B%Y)
2126    !       endif
2127    !    ENDIF
2128    ALLOCATE(B%X(N,7))
2129    ALLOCATE(B%U(0:N))
2130    ALLOCATE(B%POS(0:N))
2131    !    ALLOCATE(B%SIGMA(6))
2132    !    ALLOCATE(B%DX(3))
2133    !    ALLOCATE(B%ORBIT(6))
2134    !    ALLOCATE(B%BBPAR,B%BEAM_BEAM,B%BBORBIT)
2135    DO I=0,N
2136       NULLIFY(B%POS(i)%NODE)
2137    ENDDO
2138    !   ALLOCATE(B%CHARGE)
2139    !   ALLOCATE(B%TIME_INSTEAD_OF_S)
2140
2141    B%X  = 0.0_dp
2142    B%U  = .FALSE.
2143    !    B%CHARGE=1
2144    !    B%TIME_INSTEAD_OF_S=.FALSE.
2145
2146    !    B%SIGMA=ZERO
2147    !    B%DX=ZERO
2148    !    B%BBPAR=ZERO
2149    !    B%ORBIT=ZERO
2150    !    B%BEAM_BEAM=MY_FALSE
2151    !    B%BBORBIT=MY_FALSE
2152  END SUBROUTINE ALLOCATE_BEAM
2153
2154  SUBROUTINE KILL_BEAM(B)
2155    IMPLICIT NONE
2156    TYPE(BEAM) , INTENT (INOUT) :: B
2157    !    IF(ASSOCIATED(B%Y)) THEN
2158    !       CALL KILL(B%Y)
2159    !       DEALLOCATE(B%Y)
2160    !    ENDIF
2161    IF(ASSOCIATED(B%N)) THEN
2162       DEALLOCATE(B%N,B%LOST,B%X,B%U,B%POS)
2163       !       DEALLOCATE(B%N,B%LOST,B%X,B%U,B%POS,B%CHARGE,B%TIME_INSTEAD_OF_S)
2164       !      DEALLOCATE(B%SIGMA,B%DX,B%BBPAR,B%ORBIT,B%BEAM_BEAM,B%BBORBIT)
2165    ENDIF
2166  END SUBROUTINE KILL_BEAM
2167
2168  SUBROUTINE KILL_BEAMS(B)
2169    IMPLICIT NONE
2170    TYPE(BEAM) , INTENT (INOUT) :: B(:)
2171    INTEGER I
2172    DO I=1,SIZE(B)
2173       CALL KILL_BEAM(B(i))
2174    ENDDO
2175  END SUBROUTINE KILL_BEAMS
2176
2177
2178  FUNCTION BEAM_IN_X(B,I)
2179    IMPLICIT NONE
2180    REAL(DP) BEAM_IN_X(6)
2181    TYPE(BEAM), INTENT(INOUT) ::B
2182    INTEGER, INTENT(IN) :: I
2183
2184    BEAM_IN_X=B%X(I,1:6)
2185
2186  END  FUNCTION BEAM_IN_X
2187
2188  SUBROUTINE X_IN_BEAM(B,X,I,DL,T)
2189    IMPLICIT NONE
2190    REAL(DP),OPTIONAL:: X(6)
2191    REAL(DP),OPTIONAL:: DL
2192    TYPE(BEAM), INTENT(INOUT) ::B
2193    TYPE(INTEGRATION_NODE),OPTIONAL,POINTER :: T
2194    INTEGER, INTENT(IN) :: I
2195
2196    if(PRESENT(X)) B%X(I,1:6)=X(1:6)
2197    IF(PRESENT(DL)) B%X(I,7)=DL
2198    IF(PRESENT(T)) B%POS(I)%NODE=>T
2199    if(.not.CHECK_STABLE) then
2200       !       write(6,*) "unstable "
2201       CALL RESET_APERTURE_FLAG
2202       b%u(I)=.true.
2203       B%LOST=B%LOST+1
2204    endif
2205
2206  END  SUBROUTINE X_IN_BEAM
2207
2208  !  Beam Beam stuff
2209
2210  subroutine s_locate_beam_beam(my_ering,sc,pos,tl,b_b)
2211    implicit none
2212    type(layout), target :: my_ering
2213    type(integration_node), pointer :: tl
2214    real(dp) sc
2215    integer pos,j
2216    logical(lp) b_b
2217
2218
2219    IF(.NOT.ASSOCIATED(my_ering%T)) THEN
2220       CALL MAKE_NODE_LAYOUT(my_ering)
2221    ENDIF
2222    ! s(1) total ld
2223    ! s(2) local integration distance
2224    !          SC=MOD(SC,MY_RING%T%END%S(1))
2225    b_b=.false.
2226    TL=>my_ering%T%START
2227    DO j=1,my_ering%T%N
2228       if(pos<1) then
2229          IF(TL%S(1)<=SC.AND.TL%NEXT%S(1)>SC) then
2230             b_b=.true.
2231             exit
2232          endif
2233       else
2234          if(j==pos) then
2235             b_b=.true.
2236             exit
2237          endif
2238       endif
2239       TL=>TL%NEXT
2240    ENDDO
2241    if(b_b.and.(tl%cas==case0.or.tl%cas==caset)) then
2242       write(6,*) " Beam-Beam position at ",tl%parent_fibre%mag%name
2243       if(.not.associated(tl%BB)) call alloc(tl%BB)
2244       write(6,*) tl%pos,tl%parent_fibre%mag%name,' created'
2245       b_b=my_true
2246    else
2247       b_b=my_false
2248       write(6,*) " Beam-Beam position not found "
2249    endif
2250  end  subroutine s_locate_beam_beam
2251
2252  subroutine locate_beam_beam(my_ering,a,ent,tl,b_b)
2253    implicit none
2254    type(layout), target :: my_ering
2255    type(integration_node), pointer :: tl,tmin,tp
2256    real(dp) dmin,a(3),ent(3,3),d,cos
2257    integer pos,j
2258    logical(lp) b_b
2259
2260    dmin=mybig
2261
2262    IF(.NOT.ASSOCIATED(my_ering%T)) THEN
2263       CALL MAKE_NODE_LAYOUT(my_ering)
2264       call FILL_SURVEY_DATA_IN_NODE_LAYOUT(my_ering)
2265    ENDIF
2266    IF(.NOT.ASSOCIATED(my_ering%T%start%B)) THEN
2267       call FILL_SURVEY_DATA_IN_NODE_LAYOUT(my_ering)
2268    ENDIF
2269
2270    b_b=.false.
2271    TL=>my_ering%T%START
2272    tmin=>tl
2273    tp=>tl
2274    DO j=1,my_ering%T%N
2275       if(tl%cas==case0.or.tl%cas==caset) then
2276          d=sqrt((a(1)-tl%a(1))**2+(a(2)-tl%a(2))**2+(a(3)-tl%a(3))**2)
2277          if(d<=dmin) then
2278             dmin=d
2279             tmin=>tl
2280             tp=>tl%parent_fibre%previous%t2%previous%previous
2281          endif
2282       endif
2283       TL=>TL%NEXT
2284    ENDDO
2285    if((tmin%cas==case0.or.tmin%cas==caset)) then
2286
2287       write(6,*) " Tentative Beam-Beam position at ",tl%parent_fibre%mag%name
2288       write(6,*) tmin%pos,tmin%parent_fibre%mag%name,' created'
2289       b_b=my_true
2290
2291       cos=0.0_dp
2292       do j=1,3
2293          cos=cos+ (a(j)-tmin%a(j))*tmin%ent(1,j)
2294       enddo
2295       tl=>tp
2296       if(cos<0.0_dp) then
2297          write(6,*) " Beam-Beam position replaced at ",tl%parent_fibre%mag%name,tl%cas
2298          write(6,*) tl%pos,tl%parent_fibre%mag%name,' created'
2299       endif
2300    else
2301       b_b=my_false
2302       write(6,*) " Beam-Beam position not found "
2303    endif
2304    if(b_b.and.(.not.associated(tl%BB))) call alloc(tl%BB)
2305
2306  end  subroutine locate_beam_beam
2307
2308end module ptc_multiparticle
Note: See TracBrowser for help on using the repository browser.