source: PSPA/madxPSPA/src/madx_ptc_trackcavs.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: 26.5 KB
Line 
1module madx_ptc_trackline_module
2  use madx_ptc_module
3  use madx_ptc_intstate_module
4  use madx_ptc_setcavs_module
5  implicit none
6  save
7  public
8
9  public                              :: ptc_trackline       ! subroutine inside the module
10  public                              :: ptc_track_everystep
11
12  ! flag for debugging ranges from 0 (no debug printout) to 10 (the most detailed)
13  real(dp),allocatable :: Dismom(:,:)    ! <xnormal_(2*i-1)**(2j)>= dismon(i,j)*I_i**j
14  logical(lp)          :: onetable
15
16  !********************************************************************************************
17  !********************************************************************************************
18  !********************************************************************************************
19
20contains
21
22
23
24  subroutine ptc_track_everystep(nobs)
25    ! subroutine that performs tracking with acceleration
26    ! it is called as a result of ptc_trackline MAD-X command
27
28    implicit none
29    integer, intent (IN) :: nobs ! the maximum number of observation points >=1
30    INTEGER, ALLOCATABLE :: observedelements(:)
31    real(dp)             :: charge    ! charge of an accelerated particle
32    type(fibre), pointer :: p
33    real (dp)            :: x(1:6)
34    !    real (dp)            :: polarx(1:6)   ! track vector -
35    real (dp)            :: xp, yp,  p0
36    real (dp)            :: pathlegth = zero
37    integer              :: npart = 1
38    integer              :: n = 1
39    integer              :: nturns = 1
40    integer              :: t = 1
41    integer              :: elcode
42    logical(lp)          :: gcs
43    logical(lp)          :: rplot, dosave, doloop, isstart, isend
44!    real (dp)            :: gposx, gposy, gposz
45    integer              :: e, ni
46    integer              :: obspointnumber ! observation point number in c-code
47    integer              :: getnumberoftracks !function
48    type(internal_state) :: intstate
49    real(kind(1d0))      :: get_value
50    integer              :: get_option
51    integer, external    :: restart_sequ, & !  restart beamline and return number of beamline node
52         advance_node    !  advance to the next node in expanded sequence
53    !                    !  =0 (end of range), =1 (else)
54    REAL(KIND(1d0)), external :: node_value  !/*returns value for parameter par of current element */
55    TYPE(BEAM) :: TheBEAM
56    TYPE(INTEGRATION_NODE),POINTER :: CURR_SLICE,PREV_SLICE
57    integer             :: mf
58
59
60    !------------------------------------------------------
61    !initialization
62    npart = 1
63    n = 1
64    t = 1
65    !------------------------------------------------------
66
67    if(universe.le.0.or.EXCEPTION.ne.0) then
68       call fort_warn('return from ptc_trackline: ',' no universe created')
69       print*,"Max number of nobs ", nobs
70       return
71    endif
72    if(index_mad.le.0.or.EXCEPTION.ne.0) then
73       call fort_warn('return from ptc_trackline: ',' no layout created')
74       return
75    endif
76
77    nturns = get_value('ptc_trackline ','turns ')
78    if (getdebug() > 2) then
79        print *, 'ptc_trackline, nturns = ', nturns
80    endif
81
82    if ( (nturns > 1) .and. (my_ring%closed .eqv. .false.)) then
83       call fort_warn('WARNING: You can not make more than one turn in a line!', &
84            'Putting number of turns to 1!')
85       nturns = 1
86    endif
87   
88    onetable = get_option('onetable ') .ne. 0
89
90    gcs = get_value('ptc_trackline ','gcs ') .ne. 0
91
92    rplot = get_value('ptc_trackline ','rootntuple ') .ne. 0
93
94    intstate = getintstate()
95    if (gcs .and.  intstate%TOTALPATH==1) then
96       call fort_warn("ptc_trackline","Having global coordinates and totalpath for z is sensless")
97       gcs = .false.
98    endif
99
100
101    allocate(observedelements(1:my_ring%n)); observedelements(:)=0 ! zero means that this element is not an obs. point
102
103
104    charge = get_value('beam ', "charge ");
105    if (getdebug() > 3 ) then
106        print *, 'Read charge:', charge,' layout has charge ', my_ring%start%charge
107    endif
108
109    if (cavsareset .eqv. .false.) then
110       call setcavities(my_ring,maxaccel)
111    endif
112
113    if (getdebug() > 0) then
114        print *, 'reading tracks starting posiotions from table ....'
115    endif
116
117    call gettrack(1,x(1),x(2),x(3),x(4),x(6),x(5))
118    x(6) = -x(6)
119
120    if (getdebug() > 0) then
121        print *, 'reading.... Done'
122    endif
123
124    if (getdebug() > 0) then
125       print *, '###################################################'
126       print *, '###################################################'
127       print *, '######         TRACKING WITH PTC         ##########'
128       print *, '###################################################'
129       print *, '###################################################'
130    endif
131
132    if (rplot) then
133       call newrplot()
134    endif
135
136
137    if(.not.associated(my_ring%t))  then
138       CALL MAKE_node_LAYOUT(my_ring)
139    endif
140   
141!    c_%x_prime=.true.
142    ! Check observation points and install bb elements if any
143    if (getdebug() > 2 ) then
144      print*, "DO BEAM BEAM FLAG = ", do_beam_beam
145    endif
146   
147    e=restart_sequ()
148    p=>my_ring%start
149    do e=1, my_ring%n !in slices e goes to nelem + 1 because the last slice is the fist one.
150
151       obspointnumber=node_value('obs_point ')
152       ! instead enforce saving data at the beginning and the very end
153       !IF (e.eq.1) obspointnumber=1 ! node_value gives 0 for 1st (?)
154
155       if (obspointnumber .gt. 0) then
156          if (getdebug() > 0) then
157              print *,"Element ",e," is an observation point no. ",obspointnumber
158          endif
159          observedelements(e) = obspointnumber
160       endif
161       
162       elcode=node_value('mad8_type ')
163       
164       if (elcode .eq. 22) then
165         
166          if (getdebug() > 1 ) then
167            write(6,*) " Beam-Beam position at element named >>",p%mag%name,"<<"
168          endif
169           
170          CURR_SLICE => p%t1
171         
172          do while (.not. (CURR_SLICE%cas==case0.or.CURR_SLICE%cas==caset) )
173            if (associated(CURR_SLICE,p%t2)) exit
174            CURR_SLICE => CURR_SLICE%next
175            !print*, CURR_SLICE%cas
176          enddo
177         
178          !print *,  'BB Node Case NO: ',CURR_SLICE%cas
179         
180          if(((CURR_SLICE%cas==case0).or.(CURR_SLICE%cas==caset))) then !must be 0 or 3
181
182            if(.not.associated(CURR_SLICE%BB)) call alloc(CURR_SLICE%BB)
183
184            call getfk(xp)
185            CURR_SLICE%bb%fk = xp
186            CURR_SLICE%bb%sx = node_value('sigx ')
187            CURR_SLICE%bb%sy = node_value('sigy ')
188            CURR_SLICE%bb%xm = node_value('xma ')
189            CURR_SLICE%bb%ym = node_value('yma ')
190            CURR_SLICE%bb%PATCH=.true.
191            if (getdebug() > 2 ) then
192              print*, "BB fk=",CURR_SLICE%bb%fk
193              print*, "BB sx=",CURR_SLICE%bb%sx
194              print*, "BB sy=",CURR_SLICE%bb%sy
195              print*, "BB xm=",CURR_SLICE%bb%xm
196              print*, "BB ym=",CURR_SLICE%bb%ym
197            endif
198
199            do_beam_beam = .true.
200             
201          else
202            call fort_warn('ptc_trackline: ','Bad node case for BeamBeam')
203          endif
204         
205       endif
206       
207       obspointnumber=advance_node() ! c-code go to the next node -> the passed value is never used, just to shut up a compiler
208       p=>p%next
209    enddo
210   
211    if (getdebug() > 2 ) then
212      print*, "DO BEAM BEAM FLAG = ", do_beam_beam
213    endif
214
215    n=1
216    npart = getnumberoftracks()
217    if (getdebug() > 0) then
218        print *, 'There is ', npart,' tracks'
219    endif
220
221    !     IF(.NOT.ASSOCIATED(TheBeam%N)) THEN
222    CALL ALLOCATE_BEAM(TheBeam,npart)
223    !     ELSEIF(TheBeam%N/=npart) THEN
224    !        CALL KILL_BEAM(TheBeam)
225    !        CALL ALLOCATE_BEAM(TheBeam,npart)
226    !     ENDIF
227
228!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
229!!!!!!!!!    READS DATA FROM MADX         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
230!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
231
232    do n=1, npart
233
234       pathlegth = zero
235
236       if (getdebug() > 2 ) then
237         print *, 'Getting track ',n
238       endif 
239
240       call gettrack(n,TheBeam%X(n,1),TheBeam%X(n,2),TheBeam%X(n,3),TheBeam%X(n,4),TheBeam%X(n,6),TheBeam%X(n,5))
241
242       if (getdebug() > 1 ) then
243         write(6,'(a10,1x,i8,1x,6(f9.6,1x))') 'Track ',n,TheBeam%X(n,1:6)
244       endif
245       
246       TheBeam%X(n,7)=ZERO
247
248       if( associated(TheBeam%POS(n)%NODE) ) then
249          TheBeam%POS(n)%NODE=>my_ring%start%t1
250       endif
251
252    enddo !loop over tracks
253
254!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
255!!!!!!!!!      TRACKING       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
256!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
257
258    if (getdebug() > 1 ) then
259      call kanalnummer(mf)
260      open(unit=mf,file='thintracking_ptc.txt',POSITION='APPEND' , STATUS='UNKNOWN')
261    endif
262   
263    do t=1, nturns
264       if (getdebug() > 2 ) then
265       print*, "TURN NUMBER ",t
266       endif
267       p=>my_ring%start
268
269       PREV_SLICE => my_ring%start%T1
270       CURR_SLICE => prev_slice%next
271       e = 1
272
273       !       print*,"Name of the first element ", my_ring%start%mag%name
274       !       print*,"Position of the first element ", my_ring%start%T2%pos
275       !       print*,"Name of the last element ", my_ring%end%mag%name
276       !       print*,"Position of the last element ", my_ring%end%T2%pos
277
278       do ni=1, my_ring%end%T2%pos
279
280          if ( .not. associated(CURR_SLICE%PARENT_FIBRE, PREV_SLICE%PARENT_FIBRE) ) then
281             e = e + 1
282             p=>p%next
283          endif
284
285
286          !          call track_beam(my_ring,TheBeam,getintstate(), pos1=ni, pos2=ni+1)
287          call track_beam(my_ring,TheBeam,getintstate(), node1=ni, node2=ni+1)
288          pathlegth = curr_slice%s(3)
289
290          if (getdebug() > 2 ) then
291             !write(6,*) e, 'l=',pathlegth
292             n=1! print only first one for debug
293             write(6,'(a7,1x,i4,1x,a5,1x,i3,1x,a4,1x,f9.2,1x ,6(f12.8,1x))') &
294                     'Track ',n,'Elem ',e,'S =',pathlegth,TheBeam%X(n,1:6)
295          endif
296         
297          isstart= (t.eq.1).and.(e.eq.1)
298          isend  = (t.eq.nturns).and.(e.eq.my_ring%n)
299          dosave = (observedelements(e) .gt. 0) .or. isstart .or. isend
300          doloop = dosave .or. rplot
301          if (doloop) then
302            do n=1, npart
303
304               x = TheBeam%X(n,1:6)
305
306               p0 = p%mag%p%p0c
307
308                ! a simple hook to get a text file
309               if (getdebug() > 1 ) then
310                  write(mf,'(i8,1x, a16, 1x, 3i4, 1x,2f8.4, 1x, 7f12.8)' ) ni, p%mag%name, e, n, t, &
311                       pathlegth, TheBeam%X(n,7), &
312                       x(1), x(2) , x(3), x(4) , x(5), x(6) , p0
313               endif
314
315               if (rplot) then
316                  !For thin tracking I still do not know how to get global coordinates
317                  !   gcs = my_false !For thin tracking
318                  !   if (gcs) then
319                  !      !                write(6,'(a12,3f8.4)') "Magnet B ", p%mag%p%f%b(1), p%mag%p%f%b(2), p%mag%p%f%b(3)
320                  !      gposx = x(1)*p%chart%f%exi(1,1) + x(3)*p%chart%f%exi(1,2) + x(6)*p%chart%f%exi(1,3)
321                  !      gposy = x(1)*p%chart%f%exi(2,1) + x(3)*p%chart%f%exi(2,2) + x(6)*p%chart%f%exi(2,3)
322                  !      gposz = x(1)*p%chart%f%exi(3,1) + x(3)*p%chart%f%exi(3,2) + x(6)*p%chart%f%exi(3,3)
323                  !      !                write(6,'(a12,3f8.4)') " Rotated ", gposx,gposy,gposz
324                  !      gposx = gposx + p%chart%f%b(1)
325                  !      gposy = gposy + p%chart%f%b(2)
326                  !      gposz = gposz + p%chart%f%b(3)
327                  !
328                  !      write(6,'(a12, 2i6,3f8.4)') p%mag%name, n,e, gposx,gposy,gposz
329                  !      call plottrack(n, e, t, gposx, xp , gposy, yp , x(5), p0 , gposz)
330                  !   else
331                  !call plottrack(n, e, t, x(1), xp , x(3), yp , x(5), p0 , x(6))
332                   call plottrack(n, e, t, x(1), x(2) , x(3), x(4) , x(5), p0 , x(6))
333                  !   endif
334               endif
335
336               if ( dosave ) then
337                  if ( associated(CURR_SLICE, p%t2 ) ) then
338         !print*, "Sending to table", n, e, pathlegth
339         call putintracktable(n,t,observedelements(e),x(1), x(2) , x(3), x(4) ,x(6), x(5), pathlegth, p0)
340                  endif
341               endif
342               !fields in the table         "number", "turn", "x", "px", "y", "py", "t", "pt", "s", "e"
343
344              enddo
345         
346          endif !if doloop over tracks to save in table or root ntuple
347         
348          if (associated(CURR_SLICE%next)) then
349             PREV_SLICE => CURR_SLICE
350             CURR_SLICE => CURR_SLICE%next
351          else
352             exit;
353          endif
354
355       enddo !over elements
356
357    enddo !loop over turns
358
359   
360    if (getdebug() > 1 ) then
361      close(mf)
362    endif 
363
364    if (rplot) call rplotfinish()
365    call deletetrackstrarpositions()
366
367!    c_%x_prime=.false.
368
369    CALL KILL_BEAM(TheBeam)
370
371    deallocate (observedelements)
372    !==============================================================================
373  end subroutine ptc_track_everystep
374 !_________________________________________________________________________________
375
376   
377  subroutine putinstatustable (npart,turn,elno,elna,spos,stat,x,xini,e,mf)
378    use name_lenfi
379    implicit none
380 !   integer  :: npart,turn,nobs,stat,mf,elno
381    integer  :: npart,turn,stat,mf,elno
382!    real(kind(1d0)) :: tt
383!    character*36 table_puttab
384    character*36 table
385    character(name_len) elna
386    !hbu
387    real (dp)            :: x(1:6)
388    real (dp)            :: xini(1:6)
389    real(dp) :: spos,e
390!    integer :: get_option
391    !hbu
392    data table        / 'tracksumm        ' /
393     
394    write(mf,*) npart,spos,turn,elno,elna,xini, x, e,stat
395   
396    !"number", "turn", "x", "px", "y", "py", "t", "pt", "s", "e",
397
398    doublenum = npart
399    call double_to_table_curr(table, 'number ' , doublenum)
400    doublenum = turn
401    call double_to_table_curr(table, 'turn ' , doublenum)
402    doublenum = x(1)
403    call double_to_table_curr(table, 'x ' , doublenum)
404    doublenum = x(2)
405    call double_to_table_curr(table, 'px ' , doublenum)
406    doublenum = x(3)
407    call double_to_table_curr(table, 'y ' , doublenum)
408    doublenum = x(4)
409    call double_to_table_curr(table, 'py ' , doublenum)
410    doublenum = -x(6)
411    call double_to_table_curr(table, 't ' , doublenum)
412    doublenum = x(5)
413    call double_to_table_curr(table, 'pt ' , doublenum)
414    doublenum = spos
415    call double_to_table_curr(table, 's ' , doublenum)
416    doublenum = e
417    call double_to_table_curr(table, 'e ' , doublenum)
418
419    call augment_count(table)
420   
421  end subroutine putinstatustable
422
423  !_________________________________________________________________________________
424
425  subroutine putintracktable (npart,turn,nobs,x,px,y,py,t,pt,spos,e)
426    implicit none
427    !--- purpose: enter particle coordinates in table                      *
428    !    input:                                                            *
429    !    npart  (int)           particle number                            *
430    !    turn   (int)           turn number                                *
431    !    nobs   (int)           observation point number                   *
432    !----------------------------------------------------------------------*
433
434    !vvk
435    !      real(dp) :: tmp_coord_array(lnv), tmp_norm_array(lnv), tmp_norm
436    integer  :: npart,turn,nobs
437    real(kind(1d0)) :: tt
438    character*36 table_puttab
439    character*36 table
440    !hbu
441    real(dp) :: x,px,y,py,t,pt
442    real(dp) :: spos,e
443    !hbu
444    data table_puttab / 'track.obs$$$$.p$$$$' /
445    data table        / 'trackone           ' /
446
447    if ( onetable ) then
448       table='trackone'
449       table(9:9)= achar(0)
450    else
451       table=table_puttab
452       write(table(10:13), '(i4.4)') nobs
453       write(table(16:19), '(i4.4)') npart
454    endif
455
456
457    tt = turn
458
459    doublenum = nobs
460    call double_to_table_curr(table, 'number ' , doublenum)
461
462    doublenum = npart
463    call double_to_table_curr(table, 'number ' , doublenum)
464
465    call double_to_table_curr(table, 'turn ', tt)
466    doublenum = x
467    call double_to_table_curr(table, 'x ' , doublenum)
468
469    doublenum = px
470    call double_to_table_curr(table, 'px ', doublenum)
471
472    doublenum = y
473    call double_to_table_curr(table, 'y ' , doublenum)
474
475    doublenum = py
476    call double_to_table_curr(table, 'py ', doublenum)
477
478    doublenum = -t
479    call double_to_table_curr(table, 't ' , doublenum)
480
481    doublenum = pt
482    call double_to_table_curr(table, 'pt ', doublenum)
483
484    doublenum = spos
485    call double_to_table_curr(table, 's ' , doublenum)
486
487    doublenum = e
488    call double_to_table_curr(table, 'e ' , doublenum)
489    call augment_count(table)
490
491  end subroutine putintracktable
492
493  !_________________________________________________________________________________
494
495
496  subroutine ptc_trackline(nobs)
497    ! subroutine that performs tracking with acceleration
498    ! it is called as a result of ptc_trackline MAD-X command
499
500    implicit none
501    integer, intent (IN) :: nobs ! the maximum number of observation points >=1
502    INTEGER, ALLOCATABLE :: observedelements(:)
503    real(dp)             :: charge    ! charge of an accelerated particle
504    type(fibre), pointer :: p
505    real (dp)            :: x(1:6)
506    real (dp)            :: xini(1:6)
507    !    real (dp)            :: polarx(1:6)   ! track vector -
508    real (dp)            :: xp, yp, p0
509    real (dp)            :: pathlegth = zero
510    integer              :: npart = 1
511    integer              :: n = 1
512    integer              :: nturns = 1
513    integer              :: t = 1
514    logical(lp)          :: gcs
515    logical(lp)          :: rplot
516    real (dp)            :: gposx, gposy, gposz
517    integer              :: e
518    integer              :: apertflag
519    character(200)       :: whymsg
520    integer              :: why(9)
521    !    integer              :: rplotno
522    integer              :: obspointnumber ! observation point number in c-code
523    integer              :: getnumberoftracks !function
524    type(internal_state) :: intstate
525    real(kind(1d0))      :: get_value
526    integer              :: get_option
527    integer, external    :: restart_sequ, & !  restart beamline and return number of beamline node
528         advance_node    !  advance to the next node in expanded sequence
529    !                    !  =0 (end of range), =1 (else)
530    REAL(KIND(1d0)), external :: node_value  !/*returns value for parameter par of current element */
531    integer              :: mf
532    !type(work)           :: fen      ! Fibre ENergy
533    !------------------------------------------------------
534    !initialization
535    npart = 1
536    n = 1
537    t = 1
538    !------------------------------------------------------
539
540    if(universe.le.0.or.EXCEPTION.ne.0) then
541       call fort_warn('return from ptc_trackline: ',' no universe created')
542       print *, nobs
543       return
544    endif
545    if(index_mad.le.0.or.EXCEPTION.ne.0) then
546       call fort_warn('return from ptc_trackline: ',' no layout created')
547       return
548    endif
549
550    nturns = get_value('ptc_trackline ','turns ')
551    if (getdebug() > 2) then
552        print *, 'ptc_trackline, nturns = ', nturns
553    endif
554
555    if ( (nturns > 1) .and. (my_ring%closed .eqv. .false.)) then
556       call fort_warn('WARNING: You can not make more than one turn in a line!', &
557            'Putting number of turns to 1!')
558       nturns = 1
559    endif
560
561    onetable = get_option('onetable ') .ne. 0
562
563    gcs = get_value('ptc_trackline ','gcs ') .ne. 0
564
565    rplot = get_value('ptc_trackline ','rootntuple ') .ne. 0
566
567    intstate = getintstate()
568    if (gcs .and.  intstate%TOTALPATH==1) then
569       call fort_warn("ptc_trackline","Having global coordinates and totalpath for z is sensless")
570       call fort_warn("ptc_trackline","Disabling gcs")
571       gcs = .false.
572    endif
573
574
575    allocate(observedelements(1:my_ring%n)); observedelements(:)=0 ! zero means that this element is not an obs. point
576
577!    c_%x_prime=.true.
578
579    e=restart_sequ()
580    p=>my_ring%start
581    do e=1, my_ring%n
582
583       obspointnumber=node_value('obs_point ')
584       IF (e.eq.1) obspointnumber=1 ! node_value gives 0 for 1st (?)
585
586       if (obspointnumber .gt. 0) then
587          if (getdebug() > 0) then
588              print *,"Element ",e," is an observation point no. ",obspointnumber
589          endif
590          observedelements(e) = obspointnumber
591       endif
592
593       obspointnumber=advance_node() ! c-code go to the next node -> the passed value is never used, just to shut up a compiler
594       p=>p%next
595    enddo
596
597
598    charge = get_value('beam ', "charge ");
599    if (getdebug() > 3 ) then
600        print *, 'Read charge:', charge,' layout has charge ', my_ring%start%charge
601    endif
602
603    if (cavsareset .eqv. .false.) then
604       call setcavities(my_ring,maxaccel)
605    endif
606
607    if (getdebug() > 0) then
608        print *, 'reading tracks starting posiotions from table ....'
609    endif
610
611    call gettrack(1,x(1),x(2),x(3),x(4),x(6),x(5))
612
613    if (getdebug() > 0) then
614        print *, 'reading.... Done'
615    endif
616
617    if (getdebug() > 0) then
618       print *, '###################################################'
619       print *, '###################################################'
620       print *, '######         TRACKING WITH PTC         ##########'
621       print *, '###################################################'
622       print *, '###################################################'
623    endif
624
625    if (rplot) then
626       call newrplot()
627    endif
628
629!    call kanalnummer(mf)
630!    open(unit=mf,file='ptctracklinestatus.txt',POSITION='APPEND' , STATUS='UNKNOWN')
631
632    n=1
633    npart = getnumberoftracks()
634    if (getdebug() > 0) then
635        print *, 'There is ', npart,' tracks'
636    endif
637    do n=1, npart
638
639       pathlegth = zero
640
641       if (getdebug() > 3 ) then
642           print *, 'Getting track ',n
643       endif
644
645       call gettrack(n,x(1),x(2),x(3),x(4),x(6),x(5))
646       x(6) = -x(6)
647       if (getdebug() > 0 ) write(6,'(a10,1x,i8,1x,6(f9.6,1x))') 'Track ',n,x
648       xini = x
649       do t=1, nturns
650
651          p=>my_ring%start
652
653          do e=1, my_ring%n
654
655             !print*, p%mag%name, p%mag%P%KILL_ENT_FRINGE, p%mag%P%KILL_EXI_FRINGE,  &
656             !        p%mag%P%BEND_FRINGE,  p%mag%p%PERMFRINGE, p%mag%PERMFRINGE
657             
658             !print*, p%mag%name
659             !write(6,'(a10,1x,i8,1x,6(f12.9,1x))') 'Track ',n,x
660
661             call track(my_ring,x,e,e+1,getintstate())
662             if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
663                call fort_warn('ptc_trackline: ','DA got unstable')
664                call seterrorflag(10,"ptc_trackline ","DA got unstable ");
665                goto 100 !for the time being lets try next particle,
666                         !but most probably we will need to stop tracking and reinit
667             !goto 101
668               
669             endif
670             
671             pathlegth = pathlegth + p%mag%p%ld
672
673             if (getdebug() > 2 ) then
674                write(6,*) e, 'l=',pathlegth
675                write(6,'(5f8.4, f16.8)') x(1),x(2),x(3),x(4),x(5),x(6)
676             endif
677
678             p0 = p%mag%p%p0c
679
680             if (rplot) then
681                if (gcs) then
682                   !                write(6,'(a12,3f8.4)') "Magnet B ", p%mag%p%f%b(1), p%mag%p%f%b(2), p%mag%p%f%b(3)
683                   gposx = x(1)*p%chart%f%exi(1,1) + x(3)*p%chart%f%exi(1,2) + x(6)*p%chart%f%exi(1,3)
684                   gposy = x(1)*p%chart%f%exi(2,1) + x(3)*p%chart%f%exi(2,2) + x(6)*p%chart%f%exi(2,3)
685                   gposz = x(1)*p%chart%f%exi(3,1) + x(3)*p%chart%f%exi(3,2) + x(6)*p%chart%f%exi(3,3)
686                   !                write(6,'(a12,3f8.4)') " Rotated ", gposx,gposy,gposz
687                   gposx = gposx + p%chart%f%b(1)
688                   gposy = gposy + p%chart%f%b(2)
689                   gposz = gposz + p%chart%f%b(3)
690
691                   if (getdebug() > 3 ) write(6,'(a12, 2i6,3f8.4)') p%mag%name, n,e, gposx,gposy,gposz
692
693                   call plottrack(n, e, t, gposx, x(2) , gposy, x(4) , x(5), p0 , gposz)
694                else
695                   call plottrack(n, e, t, x(1), x(2) , x(3), x(4) , x(5), p0 , x(6))
696                endif
697             endif
698
699             if ( observedelements(e) .gt. 0) then
700                call putintracktable(n,t,observedelements(e),x(1), x(2) , x(3), x(4) , x(6), x(5), pathlegth, p0)
701             endif
702             !fields in the table         "number", "turn", "x", "px", "y", "py", "t", "pt", "s", "e"
703
704             call produce_aperture_flag(apertflag)
705             if (apertflag/=0) then
706                print *, 'Particle out of aperture!'
707
708                call ANALYSE_APERTURE_FLAG(apertflag,why)
709                Write(6,*) "ptc_trackline: APERTURE error for element: ",e," name: ",p%MAG%name
710                Write(6,*) "Message: ",messagelost
711                write(whymsg,*) 'APERTURE error: ',why
712                call fort_warn('ptc_twiss: ',whymsg)
713                call seterrorflag(10,"ptc_twiss: ",whymsg);
714               
715                goto 100 !take next track
716               
717             endif
718
719             if (e .lt. my_ring%n) p=>p%next
720
721          enddo !over elements
722
723          if (apertflag/=0) then
724             exit; !goes to the next particle
725          endif
726
727       enddo !loop over turns
728       !           npart,turn,elno,elna,spos,stat,x,xini,e,mf
729       
730       !fen = p;
731       t = t - 1
732       !call putinstatustable(n,t,e,p%previous%MAG%name,pathlegth,0,x,xini,fen%energy,mf)
733     
734100    continue!take next track
735       
736    enddo !loop over tracks
737   
738
739101 continue!finish the program
740   
741    if (rplot) call rplotfinish()
742    call deletetrackstrarpositions()
743
744    !close(mf)
745
746!    c_%x_prime=.false.
747
748    deallocate (observedelements)
749    !==============================================================================
750  end subroutine ptc_trackline
751
752
753
754
755  !_________________________________________________________________________________
756
757
758end module madx_ptc_trackline_module
759
760
761!              if (getdebug() > 3) then
762!                 write(6,*) p%mag%name
763!                 write(6,'(a12,3f8.4)') "Chart  B ", p%chart%f%b(1), p%chart%f%b(2), p%chart%f%b(3)
764!                 write(6,'(a12,3f8.4)') "Magnet B ", p%mag%p%f%b(1), p%mag%p%f%b(2), p%mag%p%f%b(3)
765!                 write(6,'(a12,3f8.4)') "Chart Exi1 ", p%chart%f%exi(1,1), p%chart%f%exi(1,2), p%chart%f%exi(1,3)
766!                 write(6,'(a12,3f8.4)') "Chart Exi2 ", p%chart%f%exi(2,1), p%chart%f%exi(2,2), p%chart%f%exi(2,3)
767!                 write(6,'(a12,3f8.4)') "Chart Exi2 ", p%chart%f%exi(3,1), p%chart%f%exi(3,2), p%chart%f%exi(3,3)
768!                 write(6,'(a12,3f8.4)') "mag Exi1 ", p%mag%p%f%exi(1,1), p%mag%p%f%exi(1,2), p%mag%p%f%exi(1,3)
769!                 write(6,'(a12,3f8.4)') "mag Exi2 ", p%mag%p%f%exi(2,1), p%mag%p%f%exi(2,2), p%mag%p%f%exi(2,3)
770!                 write(6,'(a12,3f8.4)') "mag Exi2 ", p%mag%p%f%exi(3,1), p%mag%p%f%exi(3,2), p%mag%p%f%exi(3,3)
771!              endif
772!
Note: See TracBrowser for help on using the repository browser.