source: PSPA/madxPSPA/src/madx_ptc_setcavs.f90 @ 447

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

import madx-5.01.00
File size: 16.3 KB
Line 
1module madx_ptc_setcavs_module
2  USE madx_keywords
3  USE madx_ptc_intstate_module
4  implicit none
5  public
6
7  ! flag for debugging ranges from 0 (no debug printout) to 10 (the most detailed)
8  logical(lp), public                     :: cavsareset   = .false.
9  ! flag that indicates if cavities were already set for the current setup
10
11  !********************************************************************************************
12  !********************************************************************************************
13  !********************************************************************************************
14
15contains
16  subroutine setcavities(my_ring, maxaccel)
17    implicit                none
18    type(layout),target  :: my_ring
19    type(internal_state) :: localis ! internal state to be use in this routine = intstate+totalpath+time
20    logical(lp)          :: maxaccel
21    real(dp)             :: charge    ! charge of an accelerated particle
22    !      use madx_keywords
23    integer              :: i,j!,currentelement=1       !iterators
24    integer              :: mf1,mf2
25    type(fibre), pointer :: p     ! skowron: temporary variable: current fibre
26    type(work)           :: startfen
27    type(work)           :: nfen      ! New Fibre ENergy
28    integer, pointer     :: poscav(:) !array keeping indexes of cavities
29    real(dp),allocatable :: phasecav(:) !array keeping phases of cavities
30    real(dp)             :: patchprecision=c_1d_8
31    logical(lp)          :: patchenergy=.true.
32    logical(lp)          :: patchnext=.true.
33    real(dp)             :: prevbeta0  !just a temporary real variable
34    real (dp)            :: x(1:6)   ! track vector -
35    ! here we always use closed orbit track, that is all its relative coordinates are 0
36    real(dp)             :: position=zero !synchronous particle position
37    real(kind(1d0))     :: get_value
38    !------------------------------------------------------
39100 format (a20, f10.4, a10, f10.4, a10, f10.4)
40110 format (8f10.4, l2, i3)
41120 format (i16, f16.3, f16.4, e16.3, e16.3, f16.3)
42130 format (a12, i5, a6, a30, a6, f10.4, a20, f10.4)
43    !------------------------------------------------------
44
45
46    !Below we enforce that x(6) is cT, and it is time of flight from the start
47    !we use time T=x(6)/ctime to find the time of arrival to a cavity so we can adjust its phase optimally
48   
49    !This is not needed, and even to contrary, it creates a bug:
50    !if the tracking is done without totalpath, cavities should be tuned to such x(6)=0 gives max acceleration
51
52    localis = getintstate()
53   
54!    localis = localis - nocavity0 + totalpath0
55    localis = localis - nocavity0
56    if (getdebug() > 1) then
57       print *, "I am in setcavities "
58       call print(localis,6)
59    endif
60
61    patchnext=.true.
62
63    charge = get_value('beam ', "charge ")
64
65    if (getdebug() > 1) then
66       call kanalnummer(mf1)
67       open(unit=mf1,file='sychrpart.txt')
68       call kanalnummer(mf2)
69       open(unit=mf2,file='twcavsettings.txt')
70       write(mf2,'(6a16)') "!ElNo     ","Ref.Momentum","Phase","Frequency [Hz]","Voltage","DeltaE"
71    endif
72
73    nfen = 0
74    startfen = 0
75    x(:)=zero
76
77    call locate_all_twcav(my_ring,poscav)
78    if ( getdebug() > 2 ) write(6,*) "There are ", size(poscav), " Cavities in the line."
79    if ( size(poscav) == 0) then
80       if (getdebug() > 1) then
81          close(mf1);close(mf2);
82       endif
83       return
84    endif
85
86    allocate(phasecav(size(poscav)))
87
88    !  Here is tracking element by element
89    p=>my_ring%start
90
91    startfen=p  !setting up start energy for record
92    nfen=p      ! current fibre energy
93
94    if ( getdebug() > 1 ) then
95        print *, 'c_%feed_p0c = ', c_%feed_p0c
96    endif
97
98    if ( getdebug() > 2 ) write (*,*) 'START TRACKING TILL THE FIRST CAVITY'
99
100    i = 1
101
102    do j=1,size(poscav)
103
104       if ( getdebug() > 2 ) then
105          write (*,*) 'Current cavity no is j=',j
106          write (*,*) 'Setting beam momentum AND tracking ', nfen%p0c,' till the cavity (',poscav(j),')'
107       endif
108
109
110       do i=i,poscav(j)-1 !from the current element i to the current cavity j
111
112          p = nfen   ! set current reference energy
113          call track(my_ring,x,i,i+1,localis)
114
115          if ( .not. c_%stable_da) then
116             call fort_warn('setcavities: ','DA got unstable')
117             call seterrorflag(10,"setcavities ","DA got unstable");
118
119             deallocate(poscav);
120             deallocate(phasecav);
121             if (getdebug() > 1) then
122                close(mf1);close(mf2);
123             endif
124             return
125          endif
126
127          if ( getdebug()>1 ) then
128             write (6,*) ' i=',i,' name=',p%mag%name, &
129                  ' beta0 ', nfen%beta0, &
130                  ' newpos ', x(6), &
131                  ' Current energy ',nfen%energy
132             write(6,'(6f8.4)') x
133          endif
134
135          if (getdebug() > 1) then
136             write (mf1,*) ' '
137             write (mf1,130) 'i=',i,' name=',p%mag%name,' p0c=',p%mag%p%p0c, ' Current energy ',nfen%energy
138             write (mf1,'(6f8.4)') x
139          endif
140
141          p=>p%next
142       enddo
143
144       if (associated(p%next)) then
145          if ( (p%next%mag%kind==kind21) .or. (p%next%mag%kind==kind4) ) then
146             write(6,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
147             write(6,*) "!!!                                     !!!"
148             write(6,*) "!!!    CONESCUTIVE        CAVITIES      !!!"
149             write(6,*) "!!!                                     !!!"
150             write(6,*) "!!!     Currently it is forbiden        !!!"
151             write(6,*) "!!! to place one cavity after another.  !!!"
152             write(6,*) "!!! Plese insert a marker between them. !!!"
153             write(6,*) "!!!                                     !!!"
154             write(6,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
155             stop
156          endif
157       endif
158
159       !AT THIS POINT WE HAVE ARRIVED TO A CAVITY
160       ! p point to this cavity
161
162       ! set the reference energy in this cavity
163
164       p = nfen
165
166       if ( getdebug() > 1 ) then
167          write (6,130) 'i=',i,' name=',p%mag%name,' p0c=',p%mag%p%p0c, ' Current energy ',nfen%energy
168          write (6,'(6f8.4)') x
169       endif
170
171       !TUNE CAVITY
172       call setcavity(p,x,phasecav(j),charge,maxaccel)
173
174       if (getdebug() > 1) then
175          write(mf2,120) poscav(j), p%mag%p%p0c, p%mag%phas*c_360/twopi, p%mag%freq, p%mag%volt, p%mag%delta_e
176       endif
177
178       !TRACK CAVITY
179       call track(my_ring,x,poscav(j),poscav(j)+1,localis)
180
181       if ( .not. c_%stable_da) then
182          call fort_warn('setcavities: ','DA got unstable')
183          call seterrorflag(10,"setcavities ","DA got unstable");
184
185          deallocate(poscav);
186          deallocate(phasecav);
187          if (getdebug() > 1) then
188             close(mf1);close(mf2);
189          endif
190          return
191       endif
192
193       if (getdebug() > 1) then
194          write (mf1,*) ' '
195          write (mf1,130) 'poscav(j)=',poscav(j),' name=',p%mag%name,' p0c=',p%mag%p%p0c, ' Current energy ',nfen%energy
196          write (mf1,'(6f8.4)') x
197       endif
198
199       if ( getdebug() > 2 ) then
200          write(6,'(a, 6f12.8)') ' Track parameters after cavity ',x
201          write(*,100) 'Old Fibre: energy=',nfen%energy,' momentum=',nfen%p0c,' kinetic=',nfen%kinetic
202       endif
203
204       ! GET NEW ENERGY AFTER THE CAVITY
205       prevbeta0 = nfen%beta0
206
207       nfen= x(5)*nfen%p0c
208
209       if ( getdebug() > 2 ) then
210          write(6,100) 'New Fibre: energy=',nfen%energy,' momentum=',nfen%p0c,' kinetic=',nfen%kinetic
211          write(6,110) nfen
212          write(6,'(a10, f8.4)') 'Relative E increase', (nfen%energy-startfen%energy)/startfen%kinetic
213       endif
214
215       if (getdebug()>1) write (6,*) 'beta0 ', prevbeta0, ' oldpos ', position, ' newpos ', x(6), ' Current energy ',nfen%energy
216       position = x(6)
217
218       !PATCH THE NEXT ELEMENT ON ENTRANCE
219       p%next = nfen
220
221       if ( getdebug() > 2 ) write (*,*) 'Finding patch for j=',j,' ',p%mag%name
222       call find_patch(p,next=patchnext,ENERGY_PATCH=patchenergy,PREC=patchprecision)
223
224       i=poscav(j)+1
225       p=>p%next
226
227       !from this point on we do not need to calculate TOF cause there is no further cavs to set
228    enddo
229
230    if ( getdebug() > 2 ) then
231       write (*,*) 'Loop over cavities done'
232       write (*,*) 'Current element is ', p%mag%name
233       write (*,*) 'Doing loop from the first element after the last cavity to the END'
234    endif
235
236    do i=i,my_ring%n !setting beam energies to the end of line
237       p = nfen
238       call track(my_ring,x,i,i+1,localis)
239
240       if ( .not. c_%stable_da) then
241          call fort_warn('setcavities: ','DA got unstable')
242          call seterrorflag(10,"setcavities ","DA got unstable");
243
244          deallocate(poscav);
245          deallocate(phasecav);
246          if (getdebug() > 1) then
247             close(mf1);close(mf2);
248          endif
249          return
250       endif
251
252       if (getdebug() > 1) then
253          write (mf1,*) ' '
254          write (mf1,130) 'i=',i,' name=',p%mag%name,' p0c=',p%mag%p%p0c, ' Current energy ',nfen%energy
255          write (mf1,'(6f8.4)') x
256       endif
257
258       if ( getdebug() > 1 ) then
259          write(6,*) ' i=',i,' name=',p%mag%name, &
260               ' beta0 ', nfen%beta0, &
261               ' newpos ', x(6), &
262               ' Current energy ',nfen%energy
263          write(6,'(6f8.4)') x
264       endif
265
266       p=>p%next
267    enddo
268
269
270    if (getdebug() > 1) then
271       write (mf1,*) ' '
272       write (mf1,*) 'END'
273       write (mf1,'(6f8.4)') x
274
275       write(6,*) 'PARAMETERS AT THE END OF LINE:'
276       write(6,'(a, 6f8.4)') ' Track parameters ',x
277       write(*,100) 'START energy=',startfen%energy,' momentum=',startfen%p0c,' kinetic=',startfen%kinetic
278       write(6,100) 'END energy=',nfen%energy,' momentum=',nfen%p0c,' kinetic=',nfen%kinetic
279       write(6,110) nfen
280       write(6,'(a10, f8.4)') 'Relative E increase', (nfen%energy-startfen%energy)/startfen%kinetic
281    endif
282
283    ! FINISHED HERE
284
285    patchnext=.false.
286    p=>my_ring%start
287
288    do i=1,my_ring%n
289       if ( associated(p%mag) .eqv. .false.) then
290          if (getdebug() > 1 ) then
291              print *, 'Fibre no. ',i,' has no mag assigned to it'
292          endif
293          cycle
294       endif
295       if ( getdebug() > 2 ) then
296          write(6,*) 'Name: ', p%mag%name, ' Kind: ', p%mag%kind
297       endif
298
299       if( (p%mag%kind == kind21) .or. (p%mag%kind == kind4)) then
300
301          if(p%next%patch%energy==1) then
302             p%patch%energy=2
303             p%next%patch%energy=0
304          endif
305
306          if ( getdebug() > 1 ) then
307             write (6,*) 'Cavity ',i,' name ',p%mag%name,' phase ', p%mag%phas,' Volt ',p%mag%volt, &
308                  & ' length ', p%mag%l
309             write(6,*) 'DELTAE ', p%mag%DELTA_E
310          endif
311
312
313       endif
314       p=>p%next
315    enddo
316
317    cavsareset = .true. !module field indicating that cavities were set appriopriately
318    deallocate(poscav);
319    deallocate(phasecav);
320
321    if (getdebug() > 1) then
322       close(mf1);close(mf2);
323    endif
324
325    !****************************************************************************************
326    !*********  E N D   O F   PTC_TRACKCAVS  ************************************************
327    !****************************************************************************************
328    !________________________________________________________________________________________
329
330  contains  ! what follows are internal subroutines of ptc_trackline
331    !____________________________________________________________________________________________
332
333    subroutine setcavity(f, x, phase_rel, charge, ene)
334      implicit none
335      type(fibre),intent(inout):: f         ! fiber -> here must be a cavity, i.e. kind21 (tw) or kind4 (rf)
336      real(dp)                 :: x(6)      ! reference particle coordinates (closed orbit for a circular machine)
337      real(dp)                 :: phase_rel ! final relative phase
338      real(dp)                 :: phase_old ! for printing old phase
339      real(dp)                 :: charge    ! charge of an particle
340      logical(lp)              :: ene       ! switches if cavity should always maximally accelerate
341      ! the reference track; lag is calculated
342      !      logical(lp)              :: givendene = .false. ! makes cavity always accelerate about a given value;
343      !      integer(4)               :: tmp
344      ! volt is calculated; lag and freq is preserved
345      real(dp)                 :: de_mev ! delta energy
346      real(dp)                 :: arrivtime !time of arrival
347!      type(internal_state)     :: globalis ! internal state to be use in the tracking
348
349      arrivtime = x(6)/clight
350      if (getdebug()>2) then
351          print *, 'arrivtime = ', arrivtime
352      endif
353
354      if( (f%mag%kind/=kind21) .and. (f%mag%kind/=kind4) ) then
355         write(6,*) " fatal error: not a Cavity "
356         stop
357      endif
358
359      if ( f%mag%kind==kind21) then
360         if(f%mag%cav21%psi/=zero) then
361            write(6,*) " warning: backwards wave present ",f%mag%cav21%psi
362            f%mag%cav21%psi=zero   ! removing backward waves
363            f%magp%cav21%psi=zero   ! removing backward waves
364         endif
365      endif
366
367!      globalis = getintstate()
368!      print*, "Total path is ", globalis%totalpath
369!      chargesign = charge/abs(charge)
370     
371      phase_old = f%mag%phas
372     
373      if(ene) then
374
375         if ( getdebug() > 2 ) then
376            de_mev=f%mag%volt*f%mag%l
377            write(*,*) '   Max Energy to gain: ', de_mev, ' MeV, x(6)', x(6)
378         endif
379
380         f%mag%phas = pi/two - twopi*f%mag%freq*arrivtime - f%mag%lag ! here we tune to be on the crest and then we add the lag
381         f%magp%phas= f%mag%phas
382         phase_rel=f%mag%phas
383
384      else
385
386         f%mag%phas = - f%mag%lag
387         f%magp%phas= f%mag%phas
388         phase_rel=f%mag%phas
389
390      endif
391
392
393      !    write (*,*) 'energy (t/f)? :',ene, 'charge: ', charge
394      if ( getdebug() > 1 ) then
395         write(6,*) 'Cavity settings:'
396         write(6,*)                    '    Name   ', f%mag%name
397         write(6,'(a12,f12.5,a10,l1)') '    Charge ', charge,' max ene? : ',ene
398         write(6,'(a12,f12.5,a10)')    '    Volt ',   f%mag%volt,' MV '
399         write(6,'(a12,f12.5,a10)')    '    DELTAE ', f%mag%delta_e, ' GeV '
400         write(6,'(a12,f12.5,a10)')    '    Length ', f%mag%l,' m'
401         write(6,'(a12,f12.3,a10)')    '    Phase ',  f%mag%phas, ' rad '
402         write(6,'(a12,f12.3,a10)')    '    Old Ph ', phase_old, ' rad '
403         write(6,'(a12,f12.0,a10)')    '    Freq ',   f%mag%freq, ' Hz '
404         write(6,'(a12,f12.5,a10,f12.4,a10)') '    Lag ',    f%mag%lag/twopi*c_360,' deg ', f%mag%lag,' rad '
405         write(6,'(a12,f12.5,a10)')    '    P0c ',    f%mag%p%p0c, 'GeV/c'
406      endif
407
408    end subroutine setcavity
409    !____________________________________________________________________________________________
410
411    subroutine locate_all_twcav(r,pos)
412      implicit none
413      type(layout), target, intent(inout) :: r
414      type(fibre), pointer:: p
415      integer, pointer ::  pos(:)
416      integer i,ic
417      ic=0
418      p=>r%start
419      do i=1,r%n
420         if( (p%mag%kind==kind21) .or. (p%mag%kind==kind4) ) then
421            if(p%mag%freq/=zero) then
422               ic=ic+1
423            endif
424         endif
425         p=>p%next
426      enddo
427      allocate(pos(ic))
428      pos=0
429      ic=0
430      p=>r%start
431      do i=1,r%n
432         if( (p%mag%kind==kind21) .or. (p%mag%kind==kind4) ) then
433            if(p%mag%freq/=zero) then
434               ic=ic+1
435               pos(ic)=i
436            endif
437         endif
438         p=>p%next
439      enddo
440
441    end subroutine locate_all_twcav
442
443  end subroutine setcavities
444
445end module madx_ptc_setcavs_module
446
447
448!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
449!
450!         if (f%mag%delta_e > 0) then
451!            kf= -f%mag%delta_e/f%mag%l/f%dir/charge
452!            if(givendene) then
453!               !push about given energy, calculates the needed voltage keeping lag and freq
454!               f%mag%volt=kf
455!               f%magp%volt=kf
456!            else
457!               !uses simply the peak voltage and phase
458!               f%mag%volt=sign(one,kf*f%mag%volt) * f%mag%volt
459!               f%magp%volt=f%mag%volt
460!            endif
461!            f%mag%phas=pi/two-twopi*f%mag%freq*arrivtime-c_%phase0
462!            f%magp%phas=f%mag%phas
463!         endif
464!         phase_rel=f%mag%phas+twopi*f%mag%freq*arrivtime
465
466!       tmp = f%mag%phas/twopi
467!       print *, tmp, f%mag%phas/twopi
468!       f%mag%phas = f%mag%phas - twopi*dble(tmp)
469!       f%magp%phas=f%mag%phas
470!       print *, f%mag%phas
Note: See TracBrowser for help on using the repository browser.