source: PSPA/madxPSPA/src/madx_ptc_knobs.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: 57.9 KB
Line 
1module madx_ptc_knobs_module
2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
3  ! madx_ptc_knobs module
4  ! Piotr K. Skowronski (CERN)
5  !
6  ! This module contains the routines for PTC knobs manipulation (aka pol_block's)
7  ! MAD-X command ptc_knob configures the knobs.
8  ! - addknob sets a pol_block(s) in this module
9  ! - setknobs sets the configured previously knobs on a layout
10  use madx_keywords
11  use madx_ptc_intstate_module, only : getdebug
12  use twiss0fi
13  implicit none
14
15  include "madx_ptc_knobs.inc"
16  save
17  private
18  !============================================================================================
19  !  PUBLIC INTERFACE
20  public                         :: putusertable
21  public                         :: addpush
22  public                         :: cleartables
23
24  public                         :: addknob   !adds knobs (called from c)
25  public                         :: addknobi
26  public                         :: resetknobs ! removes knobs
27  public                         :: setknobs ! sets added knobs on a layout
28  public                         :: getnknobsm ! returns number of knobs that are magnets properties (set via pol_block)
29  public                         :: getnknobis ! returns number if knobs that are initial condition at the beginning of a line
30  public                         :: getnknobsall ! returns number of all knobs getnknobsm+getnknobis
31  public                         :: resultswithknobs ! called by ptc_twiss to dump the requested results
32  public                         :: parametrictwiss  ! returns lattice functions with dependence on knobs
33  public                         :: getknobinicond
34  public                         :: getnpushes ! returns number of selected parameters to be stored
35
36  public                         :: writeparresults
37  public                         :: killparresult
38  public                         :: finishknobs
39  public                         :: twissfctn
40  public                         :: getknobsnames
41  public                         :: getfctnsnames
42  public                         :: getfunctionat
43  public                         :: getfunctionvalueat
44  public                         :: getlengthat
45  public                         :: setparvalue
46  public                         :: setknobvalue
47  public                         :: filltables
48  public                         :: setnameintbles
49
50  !============================================================================================
51  !  PRIVATE
52  !    data structures
53
54  integer, parameter                    ::  maxnpushes=100
55  type (tablepush_poly), private        ::  pushes(maxnpushes)
56  integer                               ::  npushes = 0
57  character(20), private                ::  tables(maxnpushes)  !tables names of existing pushes - each is listed only ones
58  integer,  private                     ::  ntables = 0 !number of distinctive tables
59  real(dp), private, allocatable        ::  Dismom(:,:)    ! <xnormal_(2*i-1)**(2j)>= dismon(i,j)*I_i**j
60
61
62  integer, parameter                    ::  maxnpolblocks=20
63  type (pol_block), target              ::  polblocks(maxnpolblocks)
64  integer                               ::  npolblocks=0
65  integer                               ::  nknobs=0
66  type(pol_block_inicond)               ::  knobi
67  integer                               ::  nknobi
68
69
70  integer, parameter                    ::  maxpar = 100
71  integer                               ::  nmapels = 0
72  !  type(mapelresult),target              ::  mapels(maxpar)
73
74  real(dp), allocatable                 ::  spos(:)
75  real(dp), allocatable                 ::  deltaes(:) !array with energy increase for each element with respect to the beginning
76  real(dp), allocatable                 ::  parvals(:) ! temp array with parameter values, to find out a function value for some parameters
77  !  type(taylor), allocatable             ::  results(:,:)
78  type(universal_taylor), allocatable, target   ::  results(:,:)
79
80  character(48)                         ::  twisstablename
81
82  integer, parameter                    ::  textbufferlength = 100000
83  character(textbufferlength)           ::  textbuffer
84
85
86  integer                               ::  currentrow = 0
87
88  type(real_8), private, dimension(3)   ::  testold
89  type(real_8), private, dimension(3)   ::  phase
90
91  integer, private, allocatable         ::  E(:)  !array to pick taylor coefficients with .sub. and .par.
92  type(taylor), private                 ::  ave(6,6,3)
93  type(real_8), private                 ::  tx, ty
94  type(real_8), private                 ::  dph
95  type(real_8), private                 ::  test
96  integer, private                      ::  taylorsallocated = 0
97
98  !  integer, private, dimension(6)        :: j1 = (/1,0,0,0,0,0/)
99  !  integer, private, dimension(6)        :: j2 = (/0,1,0,0,0,0/)
100  !  integer, private, dimension(6)        :: j3 = (/0,0,1,0,0,0/)
101  !  integer, private, dimension(6)        :: j4 = (/0,0,0,1,0,0/)
102  integer, private, dimension(6)        :: j5 = (/0,0,0,0,1,0/)
103  integer, private, dimension(6)        :: j6 = (/0,0,0,0,0,1/)
104  integer, private, dimension(6,6)       :: fo = &
105       reshape(          (/1,0,0,0,0,0,&
106       0,1,0,0,0,0,&
107       0,0,1,0,0,0,&
108       0,0,0,1,0,0,&
109       0,0,0,0,1,0,&
110       0,0,0,0,0,1 /), &
111       (/6,6/) )
112
113  !============================================================================================
114  !  PRIVATE
115  !    routines
116  private                               :: augment_counts
117  private                               :: issuchtableexist
118  private                               :: putnameintables
119  private                               :: allocparresult
120
121
122
123
124
125
126contains
127  !____________________________________________________________________________________________
128
129  function getknobinicond()
130    implicit none
131    type(pol_block_inicond) :: getknobinicond
132
133    getknobinicond%beta = knobi%beta
134    getknobinicond%alfa = knobi%alfa
135    getknobinicond%dispersion = knobi%dispersion
136
137  end function getknobinicond
138
139  !____________________________________________________________________________________________
140
141
142  subroutine nullify(inpbic)
143    implicit none
144    type(pol_block_inicond)  :: inpbic
145
146    inpbic%alfa = 0
147    inpbic%beta = 0
148    inpbic%dispersion = 0
149
150  end subroutine nullify
151
152  subroutine putusertable(n,name,s,denergy,y,A)!here y is transfermap and A is a_ (in twiss it is called y)
153    !puts the coefficients in tables as defined in array pushes
154    implicit none
155    integer              :: n !fibre number
156    character(*)         :: name !fibre name
157    real(kind(1d0))      :: s  !position along the orbit
158    real(dp)             :: denergy
159    type(real_8),target  :: y(6)!input 6 dimensional function (polynomial) : Full MAP: A*YC*A_1
160    type(real_8),target  :: A(6)!input 6 dimensional function (polynomial) : Full MAP: A*YC*A_1
161    type(real_8),pointer :: e !element in array
162    real(kind(1d0))      :: coeff
163    integer              :: i !iterator
164    logical              :: pblockson
165
166    if (getdebug() > 3) then
167       print *,"madx_ptc_tablepush :putusertable n ",n," name ", name
168    endif
169
170
171    if ((getnknobsall() > 0) .and. (currentrow > 0)) then
172       !if there are any knobs and knobs where initialized (i.e. currentrow > 0)
173       ! otherwise the results table is not allocated
174       pblockson = .true.  !it means there are parameters on
175       spos(currentrow) = s
176       deltaes(currentrow) = denergy
177       call parametric_coord(A)
178       call parametrictwiss(A)
179    else
180       pblockson = .false.
181    endif
182   
183    if (name /= '$end$ ') then
184      call putnameintables()
185    else
186      call setnameintbles('$end$ ')!twiss gets names automatically from C element iterator
187      if (getdebug()>3) then
188        print*,"PTC_NORMAL mode"
189      endif
190    endif 
191
192    if (c_%npara == 6) then
193       call sixdmode()
194    else
195       do i=1,npushes
196
197          e => y(pushes(i)%element)
198
199          if (pushes(i)%pushtab) then
200             coeff  = e.sub.(pushes(i)%monomial)
201             if (getdebug()>3) then
202                write(6,'(a13, a10, a3, f9.6, a10, i1, 5(a13), i3)') &
203                     &        "4or5D putting",pushes(i)%monomial,"=",coeff," arr_row ", pushes(i)%element,&
204                     &        " in table ", pushes(i)%tabname," at column ", pushes(i)%colname, &
205                     &        " for fibre no ",n
206             endif
207
208             call double_to_table_curr(pushes(i)%tabname, pushes(i)%colname, coeff);
209          endif
210
211          if ( pblockson .and. (pushes(i)%index > 0) ) then
212
213             results(currentrow,pushes(i)%index) = e.par.(pushes(i)%monomial)
214
215             if (getdebug()>3) then
216                write(6,*) &
217                     &        "4or5D putting ",pushes(i)%monomial," arr_row ", pushes(i)%element,&
218                     &        " at named ", pushes(i)%colname, &
219                     &        " for fibre no ",n
220                write(6,*) "currentrow is ", currentrow," index ",pushes(i)%index
221                call print(results(currentrow,pushes(i)%index),6)
222             endif
223          endif
224
225       enddo
226    endif
227
228    call augment_counts()
229
230    if (currentrow > 0) then
231       !if currnet row == 0 it means that knobs were not initialized
232       currentrow = currentrow + 1
233    endif
234
235    !____________________________________________________________________________________________
236  contains
237    !____________________________________________________________________________________________
238    subroutine sixdmode()
239      implicit none
240      character bufchar
241      character*(6) monstr
242
243      do i=1,npushes
244
245         if (pushes(i)%element < 5) then
246            e => y(pushes(i)%element)
247         else
248            if (pushes(i)%element == 5) then
249               e => y(6) !6th coordinate  is d(cT) or cT
250            else
251               e => y(5) !5th coordinate is dp/p
252            endif
253         endif
254
255         monstr = pushes(i)%monomial
256         bufchar = monstr(5:5)
257         monstr(5:5) = monstr(6:6)
258         monstr(6:6) = bufchar
259
260         if (pushes(i)%pushtab) then
261            coeff = e.sub.monstr
262
263            if (getdebug()>3) then
264               write(6,'(a21, a10, a3, f9.6, a10, i1, 5(a13), i3)') &
265                    &        "6D mode Put in table ",pushes(i)%monomial,"=",coeff," arr_row ", pushes(i)%element,&
266                    &        " in table ", pushes(i)%tabname," at column ", pushes(i)%colname, &
267                    &        " for fibre no ",n
268            endif
269
270            call double_to_table_curr(pushes(i)%tabname, pushes(i)%colname, coeff);
271         endif
272
273         if ( pblockson .and. (pushes(i)%index > 0) ) then
274
275            results(currentrow,pushes(i)%index) = e.par.monstr
276
277
278            if (getdebug()>3) then
279               write(6,*) &
280                    &        "6Dmode parametric ",pushes(i)%monomial," arr_row ", pushes(i)%element,&
281                    &        " at named ", pushes(i)%colname, &
282                   &        " for fibre no ",n
283               write(6,*) "currentrow is ", currentrow," index ",pushes(i)%index
284               call print(results(currentrow,pushes(i)%index),6)
285            endif
286
287         endif
288
289      enddo
290    end subroutine sixdmode
291
292    !____________________________________________________________________________________________
293
294  end subroutine putusertable
295  !____________________________________________________________________________________________
296
297  subroutine fillusertables()
298    implicit none
299    integer                         :: i,e
300    integer                         :: nelems
301    real(kind(1d0))                 :: coeff
302    type(universal_taylor), pointer :: t
303
304    call cleartables()
305
306    nelems = ubound(results,1)
307    do e=1, nelems
308       do i=1,npushes
309          if ( (pushes(i)%pushtab) .and. (pushes(i)%index > 0) ) then
310             t => results(e,pushes(i)%index)
311             coeff = gettaylorvalue(t)
312             call double_to_table_curr(pushes(i)%tabname, pushes(i)%colname, coeff);
313          endif
314       enddo
315
316       call augment_counts()
317
318    enddo
319  end subroutine fillusertables
320  !____________________________________________________________________________________________
321
322  subroutine filltwisstable()
323    !puts the coefficients in tables as defined in array pushes
324    implicit none
325    integer              :: i,ii !iterator
326    integer              :: nelems !iterator
327    integer, parameter   :: n = 6 !counts
328    integer, parameter   :: fillntwisses  = disp4 - beta11 + 1
329    integer, parameter   :: ntwissesde = gama33 - beta11 + 1
330    real(kind(1d0))      :: opt_fun(fundim)
331    type(universal_taylor), pointer :: t
332
333    if (.not. ALLOCATED(results)) then
334       return
335    endif
336
337    call reset_count(twisstablename)
338
339    nelems = ubound(results,1)
340
341    if (nelems .lt. (currentrow-1)) then
342       call fort_warn("filltwisstable",&
343            "It seems the last ptc_twiss has failed")
344       nelems = currentrow - 1
345    endif
346
347
348
349    do i=1, nelems
350
351       do ii=beta11,ntwisses
352          t => results(i,ii)
353          opt_fun(ii)=  gettaylorvalue(t)
354       enddo
355
356       do ii=beta11,ntwissesde
357          opt_fun(ii) = opt_fun(ii)*deltaes(i)
358       enddo
359
360       call vector_to_table_curr(twisstablename, 'beta11 ', opt_fun(beta11), fillntwisses)
361       call vector_to_table_curr(twisstablename, 'x ', opt_fun(kn_x), n)
362
363       call augmentcountonly(twisstablename)
364
365    enddo
366
367  end subroutine filltwisstable
368
369  !____________________________________________________________________________________________
370
371  subroutine inittables()
372    implicit none
373
374    !     print*, " no=",c_%no
375    !     print*, " nv=",c_%nv
376    !     print*, " nd=",c_%nd
377    !     print*, " nd2=",c_%nd2
378    !     print*, " ndpt=",c_%ndpt
379    !     print*, " npara_fpp=",c_%npara_fpp
380    !     print*, " npara=",c_%npara
381    !     print*, " np_pol=",c_%np_pol
382    !     print*, " setknob=",c_%setknob
383    !
384
385    deallocate(dismom)
386    allocate(dismom(c_%nd,0:c_%no/2))
387
388  end subroutine inittables
389  !____________________________________________________________________________________________
390
391  subroutine cleartables()
392    implicit none
393    integer  :: i ! iterator
394    ! we do not deal here with the main twiss table, hence we do not clear it
395    do i=1,ntables
396       if (getdebug()>3) then
397           print *,"Clearing ",tables(i)
398       end if
399       call reset_count(tables(i))
400    enddo
401
402  end subroutine cleartables
403  !____________________________________________________________________________________________
404
405  subroutine augment_counts()
406    implicit none
407    integer  :: i ! iterator
408
409    do i=1,ntables
410       if (getdebug()>3) then
411           print *,"Augmenting ",tables(i)
412       end if
413       call augmentcountonly(tables(i)) !we need to use special augement,
414       !cause the regular one looks for for
415       !variables names like columns to fill the table
416    enddo
417
418  end subroutine augment_counts
419  !____________________________________________________________________________________________
420  subroutine setnameintbles(name)
421    implicit none
422    character(*)         :: name !fibre name
423    integer       :: i ! iterator
424   
425    do i=1,ntables
426       if (getdebug()>2) then
427           print *,"Putting name in ",tables(i)
428       end if
429       call string_to_table_curr(tables(i),"name ",name)
430    enddo
431
432  end subroutine setnameintbles
433  !____________________________________________________________________________________________
434  subroutine putnameintables()
435    implicit none
436    integer       :: i ! iterator
437   
438    do i=1,ntables
439       if (getdebug()>2) then
440           print *,"Putting name in ",tables(i)
441       end if
442       !puts automatically name of the current element
443       !does not work for ptc_normal
444       call string_to_table_curr(tables(i),"name ","name ")
445    enddo
446
447  end subroutine putnameintables
448  !____________________________________________________________________________________________
449
450
451  subroutine addpush(table,column,element,monomial)
452    use twissafi
453    implicit none
454    integer   table(*)
455    integer   column(*)
456    integer   element
457    integer   monomial(*)
458    logical   addtable
459    logical   parametric
460    real(kind(1d0))            :: get_value
461    character(48) charconv
462
463
464    parametric = get_value('ptc_select ','parametric ') .ne. 0
465    if ( (parametric .eqv. .false.) .and. (table(1) == 0) ) then
466       call fort_warn("addpush","Neither table specified neither parametric value. Ignoring")
467       return
468    endif
469
470    npushes = npushes + 1
471
472    pushes(npushes)%tabname = charconv(table)
473    pushes(npushes)%colname = charconv(column)
474
475    pushes(npushes)%element = element
476    pushes(npushes)%monomial = charconv(monomial)
477    !imput "int string"  has the length of the string at the first plave
478    pushes(npushes)%tabname(table(1)+1:table(1)+1)=achar(0)
479    pushes(npushes)%colname(column(1)+1:column(1)+1)=achar(0)
480    pushes(npushes)%monomial(monomial(1)+1:monomial(1)+1)=achar(0)
481
482    if (column(1) == 0) then
483       write(pushes(npushes)%colname,'(i1,a1,a6)') element,"_",pushes(npushes)%monomial
484    endif
485
486
487    if (table(1) > 0) then
488       pushes(npushes)%pushtab = .true.
489       addtable = .not. issuchtableexist(pushes(npushes)%tabname) !add to table to the list onl
490
491    else
492       pushes(npushes)%pushtab = .false.
493       addtable = .false. !no table name
494    endif
495
496    if (parametric) then
497       nmapels = nmapels + 1
498       pushes(npushes)%index = ntwisses+nmapels
499    else
500       pushes(npushes)%index = 0
501    endif
502
503
504    if (getdebug()>3) then
505       print  *,"madx_ptc_tablepush : addpush(<",&
506            &          pushes(npushes)%element,">,<",pushes(npushes)%monomial,">)"
507       print  *,"madx_ptc_tablepush : colname <",pushes(npushes)%colname,">"
508       print  *,"madx_ptc_tablepush : parametric results index ", pushes(npushes)%index
509       if (pushes(npushes)%pushtab) then
510          print  *,"madx_ptc_tablepush : table <",pushes(npushes)%tabname,">"
511       else
512          print  *,"madx_ptc_tablepush : not pushing to table"
513       endif
514    endif
515
516
517    if ( addtable) then
518       ntables = ntables + 1
519       tables(ntables) = pushes(npushes)%tabname
520       if (getdebug()>3)  then
521          print *,"Table has been added to the tables list ", tables(ntables)
522       endif
523
524    endif
525
526  end subroutine addpush
527  !____________________________________________________________________________________________
528
529
530  logical(lp) function issuchtableexist(tname)
531    implicit none
532    character(20) :: tname !name of the table to be checked if already is listed in table names array
533    integer       :: i! iterator
534
535    issuchtableexist = .false.
536
537    do i=1, ntables
538       if (tables(i) == tname) then
539          issuchtableexist = .true.
540          return
541       endif
542    enddo
543
544  end function issuchtableexist
545
546
547  !____________________________________________________________________________________________
548  subroutine addknobi(nameIA)
549    use twissafi
550    implicit none
551    integer     :: nameIA(*)
552    character(48)              :: name
553    character(48) charconv
554
555    if (nknobi >= maxnpolblocks) then
556       call fort_warn("addknob","Can not add more knobs, array with initial knobs if full")
557       return
558    endif
559
560    name = charconv(nameIA)
561    select case(name(1:6))
562    case ('beta11')
563       knobi%beta(1) = nknobi+1
564    case ('beta22')
565       knobi%beta(2) = nknobi+1
566    case ('beta33')
567       knobi%beta(3) = nknobi+1
568    case ('alfa11')
569       knobi%alfa(1) = nknobi+1
570    case ('alfa22')
571       knobi%alfa(2) = nknobi+1
572    case ('alfa33')
573       knobi%alfa(3) = nknobi+1
574    case ('disp1')
575       knobi%dispersion(1) = nknobi+1
576    case ('disp2')
577       knobi%dispersion(2) = nknobi+1
578    case ('disp3')
579       knobi%dispersion(3) = nknobi+1
580    case ('disp4')
581       knobi%dispersion(4) = nknobi+1
582
583    case default
584       print*, "Variable not recognized"
585       print*,"parameter ",name,"not recognized"
586       return
587    end select
588
589    nknobi = nknobi + 1
590
591
592  end subroutine addknobi
593  !____________________________________________________________________________________________
594  subroutine addknob(fibrenameIA)
595    use twissafi
596    implicit none
597    integer     :: fibrenameIA(*)
598    character(48)              :: fibrename
599    integer     :: nint, ndble, k, int_arr(nmax), char_l(nmax), i
600    real(kind(1d0))            :: d_arr(nmax)
601    character(400)             :: char_a
602    type (pol_block), pointer  :: pb
603    logical(lp) :: exactmatch
604    real(kind(1d0))            :: get_value
605    character(48) charconv
606
607
608    if (npolblocks >= maxnpolblocks) then
609       call fort_warn("addknob","Can not add more knobs, array with pol_blocks if full")
610       return
611    endif
612
613    npolblocks = npolblocks + 1
614    pb => polblocks(npolblocks)
615    pb = 0
616
617    fibrename = charconv(fibrenameIA)
618
619    k = index(fibrename,':')
620    if (k > 0) then
621       pb%name = fibrename(1:k-1)
622    else
623       pb%name = fibrename
624    endif
625
626    if (getdebug() > 1) then
627        print *,"addknob: pb%name is ", pb%name," npolblocks=",npolblocks
628    end if
629
630
631    exactmatch = get_value('ptc_knob ','exactmatch ') .ne. 0
632    if (exactmatch) then
633       if (getdebug() > 1) then
634           print*,"addknob: Using Exact name match: ", fibrename
635       end if
636       pb%vorname = fibrename
637    else
638       pb%n_name =  len_trim(pb%name)
639       if (getdebug() > 1) then
640          print*,"addknob: Using Not Exact name match:"
641          print*,"    all elements starting with ", pb%name
642          print*,"    number of first letters    ", pb%n_name
643       endif
644    endif
645
646    call comm_para('kn ', nint, ndble, k, int_arr, d_arr, char_a, char_l)
647    if (getdebug()>2) then
648        print*, "there is ",nint, " kn's: ", int_arr(1:nint)
649    end if
650    do i = 1, nint
651       if (int_arr(i) < 0) then
652          exit
653       endif
654       int_arr(i) = int_arr(i) + 1 !madx numerates from 0
655       nknobs = nknobs + 1
656       pb%ibn(int_arr(i)) = nknobs
657       if (getdebug()>0) then
658           print*, "Set normal mulitpole component ", int_arr(i),"as ", nknobs, "parameter of PTC"
659       end if
660    enddo
661
662    call comm_para('ks ', nint, ndble, k, int_arr, d_arr, char_a, char_l)
663    if (getdebug()>2) then
664        print*, "there is ",nint, " ks's: ", int_arr(1:nint)
665    end if
666    do i = 1, nint
667       if (int_arr(i) < 0) then
668          exit
669       endif
670       int_arr(i) = int_arr(i) + 1 !madx numerates from 0
671       nknobs = nknobs + 1
672       pb%ian(int_arr(i)) = nknobs
673       if (getdebug()>0) then
674           print*, "Set skew mulitpole component ", int_arr(i)," as ", nknobs, "parameter of PTC"
675       end if
676    enddo
677
678  end subroutine addknob
679  !____________________________________________________________________________________________
680
681  subroutine setknobs(alayout)
682    use twissafi
683    implicit none
684    type(layout),target      :: alayout
685    type (pol_block), pointer  :: pb
686    integer     :: i,j,k
687    character(48) charconv
688
689!    if (.not. associated(alayout)) then
690!       call fort_warn("setknobs","Passed pointer to a layout is not associated")
691!       return
692!    endif
693
694    twisstablename = table_name !defined in twissa.fi
695
696    if (getdebug() > 2 ) then
697       print*, "setknobs: There are ", npolblocks, "pol_blocks"
698    endif
699
700    if ((npolblocks == 0) .and. (nknobs == 0) .and. (nknobi == 0)) then
701       if (ALLOCATED(results)) call killparresult()
702       currentrow = -1
703       return
704    endif
705
706    do i=1, npolblocks
707       pb => polblocks(i)
708       if (getdebug() > 2 ) then
709          print*, "setknobs: Setting pol_block ", i," for ",pb%name
710       endif
711       alayout = pb
712    enddo
713
714    if (ALLOCATED(results)) call killparresult()
715
716    call allocparresult(alayout%n)
717
718    do i=1,c_%nd2
719       do j=1,c_%nd2
720          do k=1,c_%nd
721             call alloc(ave(i,j,k))
722          enddo
723       enddo
724    enddo
725
726    call alloc(phase)
727    call alloc(testold)
728    call alloc(test)
729    call alloc(dph)
730    call alloc(tx, ty)
731
732    taylorsallocated = 1
733
734    currentrow = 1
735
736    !    print*, "setknobs: All pol_blocks set"
737
738  end subroutine setknobs
739  !_________________________________________________________________________________
740  subroutine allocparresult(n)
741    implicit none
742    integer       :: n
743    integer       :: nr, np
744    integer       :: i,j
745
746    np = c_%nv - c_%npara
747    if (np <= 0) then
748       call fort_warn("addpush","Number of parameters is 0")
749       currentrow = -1  ! this is the signal that initialization has failed
750       return;
751    endif
752
753    allocate(parvals(np))
754    parvals(:) = zero
755
756    nr = ntwisses+nmapels
757    allocate(results(n,nr))
758    do i=1,n
759       do j=1,nr
760          results(i,j)=0
761       enddo
762    enddo
763    allocate(spos(n))
764    allocate(deltaes(n))
765    allocate(e(c_%npara_fpp))
766
767
768  end subroutine allocparresult
769  !_________________________________________________________________________________
770
771  subroutine killparresult
772    implicit none
773    integer       :: i,j
774
775    if (.not. ALLOCATED(results)) then
776       return
777    endif
778
779    !    remnestant of the taylor type
780    if(getdebug() > 2) then
781       print*, "killparresult: Shape of the current array: "
782       print*, "1D",lbound(results,1),ubound(results,1)
783       print*, "2D",lbound(results,2),ubound(results,2)
784    endif
785
786    do i=lbound(results,1),ubound(results,1)
787       do j=lbound(results,2),ubound(results,2)
788          !          print*,i,j
789          results(i,j) = -1
790          !           call kill(results(i,j))
791       enddo
792    enddo
793
794    deallocate(spos)
795    deallocate(deltaes)
796    deallocate(results)
797    deallocate(parvals)
798    deallocate(e)
799
800    currentrow = 0
801
802  end subroutine killparresult
803  !____________________________________________________________________________________________
804
805  subroutine resetknobs()
806    implicit none
807    integer     :: i
808
809    call nullify(knobi)
810    nknobi = 0
811
812    do i=1, npolblocks
813       polblocks(i) = 0
814    enddo
815    nknobs = 0
816    npolblocks = 0
817  end subroutine resetknobs
818  !____________________________________________________________________________________________
819
820  function getnknobsm()
821    implicit none
822    integer     :: getnknobsm
823
824    getnknobsm = nknobs
825
826  end function getnknobsm
827  !____________________________________________________________________________________________
828
829  function getnknobis()
830    implicit none
831    integer     :: getnknobis
832
833    getnknobis = nknobi
834
835  end function getnknobis
836  !____________________________________________________________________________________________
837
838  function getnknobsall()
839    implicit none
840    integer     :: getnknobsall
841
842    getnknobsall = nknobs + nknobi
843
844  end function getnknobsall
845  !____________________________________________________________________________________________
846
847  function getnpushes()
848    implicit none
849    integer     :: getnpushes
850
851    getnpushes = npushes
852
853  end function getnpushes
854  !____________________________________________________________________________________________
855
856  subroutine getfctnsnames
857    implicit none
858    character*(100) ::name
859    character*(8) ::twname
860    integer i,j
861    integer last
862
863    do i=1, ntwisses
864       twname = twissnames(i) ! twissname now known to be of fixed length 7
865       twname(8:8) = achar(0)
866       call madxv_setfctnname(i,twname)
867    enddo
868
869    do i=ntwisses, ntwisses+nmapels
870       do j=1,npushes
871          if (pushes(j)%index /= i ) cycle
872
873          name = pushes(j)%colname
874
875          last = len_trim(name) + 1
876          if (last > 100) last = 100
877          name(last:last) = achar(0)
878
879          call madxv_setfctnname(i, name)
880
881       enddo
882
883    enddo
884
885
886  end subroutine getfctnsnames
887  !____________________________________________________________________________________________
888
889  subroutine getknobsnames()
890    implicit none
891    character*(100) ::name
892    character*(48)  ::pname
893    integer last
894    integer i,j,n
895
896    n = 0
897
898    !magnets properties
899    do i=1, npolblocks
900       do j=1, nmax
901          if (polblocks(i)%ian(j) /= 0 ) then
902             n = n + 1
903
904             if (polblocks(i)%vorname == ' ') then
905                write(name,'(A1,i2.2,A2,A,A,I2)') "p",polblocks(i)%ian(j),": ",& !convertion to madx counting
906                     polblocks(i)%name(1:len_trim(polblocks(i)%name)), "* skew ",j-1
907             else
908                write(name,'(A1,i2.2,A2,A,A,I2)') "p",polblocks(i)%ian(j),": ",& !convertion to madx counting
909                     polblocks(i)%vorname(1:len_trim(polblocks(i)%vorname)), " skew ",j-1
910             endif
911             print*, "n=",n, name
912
913             last = len_trim(name) + 1
914             if (last > 100) last = 100
915             name(last:last) = achar(0)
916
917             call madxv_setknobname(n,name)
918          endif
919
920          if (polblocks(i)%ibn(j) /= 0 ) then
921             n = n + 1
922             if (polblocks(i)%vorname == ' ') then
923                write(name,'(A1,i2.2,A2,A,A,I2)') "p",polblocks(i)%ibn(j),": ",& !convertion to madx counting
924                     polblocks(i)%name(1:len_trim(polblocks(i)%name)) , "* normal ",j-1
925             else
926                write(name,'(A1,i2.2,A2,A,A,I2)') "p",polblocks(i)%ibn(j),": ",& !convertion to madx counting
927                     polblocks(i)%vorname(1:len_trim(polblocks(i)%vorname)) , " normal ",j-1
928             endif
929             print*, "n=",n, name
930
931             last = len_trim(name) + 1
932             if (last > 100) last = 100
933             name(last:last) = achar(0)
934             call madxv_setknobname(n,name)
935          endif
936
937       enddo
938    enddo
939
940    !initial conditions
941    do i=1, nknobi
942       n = n + 1
943       pname = 'NULL'
944       if (knobi%beta(1) == i) pname='initial beta11'
945       if (knobi%beta(2) == i) pname='initial beta22'
946       if (knobi%beta(3) == i) pname='initial beta33'
947       if (knobi%alfa(1) == i) pname='initial alfa11'
948       if (knobi%alfa(2) == i) pname='initial alfa22'
949       if (knobi%alfa(3) == i) pname='initial alfa33'
950       if (knobi%dispersion(1) == i) pname='initial disp1'
951       if (knobi%dispersion(2) == i) pname='initial disp2'
952       if (knobi%dispersion(3) == i) pname='initial disp3'
953       if (knobi%dispersion(4) == i) pname='initial disp3'
954
955       if (pname(1:4) == 'NULL') then
956          print*,"Knob of Initial Condition ", i," can not be recognized"
957          cycle
958       endif
959
960       write(name,'(A1,i2.2,A2,A15)') "p",n,": ",pname
961       last = len_trim(name) + 1
962       if (last > 48) last = 48
963       name(last:last) = achar(0)
964       call madxv_setknobname(n,name)
965    enddo
966
967
968  end subroutine getknobsnames
969  !____________________________________________________________________________________________
970
971  subroutine getfunctionat( el, n)
972    implicit none
973    integer n, el
974    integer eqlength
975
976    if (.not. ALLOCATED(results)) then
977       return
978    endif
979
980
981    if ( (el < lbound(results,1)) .or. (el > ubound(results,1)) ) then
982       return
983    endif
984
985    if ( (n < lbound(results,2)) .or. (n > ubound(results,2)) ) then
986       return
987    endif
988
989    call getpareq(results(el,n),textbuffer)
990
991    eqlength = ( LEN_TRIM(textbuffer) + 1 )
992
993    if (eqlength > textbufferlength) eqlength = textbufferlength
994    textbuffer(eqlength:eqlength) = achar(0)
995    call madxv_setfunctionat(el, n, textbuffer)
996
997
998  end subroutine getfunctionat
999  !____________________________________________________________________________________________
1000
1001
1002
1003  !____________________________________________________________________________________________
1004
1005  function getfunctionvalueat( n, el)
1006    implicit none
1007    real(dp)                        :: getfunctionvalueat
1008    integer                         :: n, el
1009    type(universal_taylor), pointer :: t
1010
1011    if (.not. ALLOCATED(results)) then
1012       return
1013    endif
1014
1015    if ( (el < lbound(results,1)) .or. (el > ubound(results,1)) ) then
1016       return
1017    endif
1018
1019    if ( (n < lbound(results,2)) .or. (n > ubound(results,2)) ) then
1020       return
1021    endif
1022
1023    getfunctionvalueat = zero
1024
1025    t => results(el,n)
1026
1027    if (.not. associated(t)) then
1028       print*,"Getting function ",n," at el ",el
1029       getfunctionvalueat = zero
1030    else
1031       getfunctionvalueat = gettaylorvalue(t)
1032    endif
1033
1034
1035  end function getfunctionvalueat
1036  !____________________________________________________________________________________________
1037
1038  function gettaylorvalue(ut)
1039    implicit none
1040    real(dp)              :: gettaylorvalue
1041    type(universal_taylor), pointer :: ut
1042    integer i,ii,np
1043    real(dp)   ::  c,p
1044
1045
1046    if (.not. associated(ut)) then
1047       call fort_warn("gettaylorvalue",&
1048            "provided taylor is not associated")
1049
1050       gettaylorvalue = zero
1051       return
1052    endif
1053
1054    if (.not. associated(ut%n)) then
1055       call fort_warn("gettaylorvalue",&
1056            "provided taylor is not associated")
1057
1058       gettaylorvalue = zero
1059       return
1060    endif
1061
1062    if (.not. allocated(parvals)) then
1063       call fort_warn("gettaylorvalue",&
1064            "array with parameter values is not allocated")
1065
1066       gettaylorvalue = zero
1067       return
1068    endif
1069
1070
1071    if(ut%n == 0) then
1072       if (getdebug() > 3) then
1073          print*,"no coefficients in the taylor"
1074       endif
1075
1076       gettaylorvalue = zero
1077       return
1078    endif
1079
1080    gettaylorvalue = zero
1081
1082    do i = 1,ut%n
1083
1084       c = ut%c(i)
1085       !       print*, "coef",i," is ",c
1086
1087       do ii = c_%npara + 1, ut%nv
1088
1089          if (ut%j(i,ii) /= 0) then
1090
1091             np =  ii - c_%npara
1092             p = parvals(np)
1093             !           print*, "par",np," is ",p
1094             if (ut%j(i,ii) /= 1) then
1095                p = p**ut%j(i,ii)
1096                !             print*, "par",np," to power",ut%j(i,ii), " is ",p
1097             endif
1098
1099             c = c*p
1100             !           print*,"coef",i, "X par",np," to power",ut%j(i,ii), " is ",c
1101
1102          endif
1103
1104       enddo
1105       !       print*,"COEF",i, " is ",c
1106       gettaylorvalue = gettaylorvalue + c
1107       !       print*,"FUNCTION after ", i," is ",gettaylorvalue
1108
1109    enddo
1110
1111
1112
1113  end function gettaylorvalue
1114
1115  !____________________________________________________________________________________________
1116  !____________________________________________________________________________________________
1117  !____________________________________________________________________________________________
1118
1119  subroutine twissfctn(y,ave)   !  Computes <x_i x_j> assuming linearity and no parameters
1120    implicit none
1121    type(real_8) y(6)
1122    real(dp) ave(6,6,3)
1123    integer e(6)
1124    integer i,j,k
1125
1126    print*,"c_%nd2 is ", c_%nd2
1127    print*,"c_%nd is ", c_%nd
1128
1129    e=0
1130    ave=zero
1131    do i=1,c_%nd2
1132       do j=i,c_%nd2
1133          do k=1,c_%nd
1134             e(k*2-1)=1
1135             ave(i,j,k)=ave(i,j,k)+(y(i).sub.e)*(y(j).sub.e)
1136             e(k*2-1)=0
1137             e(k*2)=1
1138             ave(i,j,k)=ave(i,j,k)+(y(i).sub.e)*(y(j).sub.e)
1139             e(2*k)=0
1140             ave(j,i,k)=ave(i,j,k)
1141
1142             !              if (ave(i,j,k) /= zero) then
1143             !                write(6,'(a3,i1,i1,i1,a3,f16.8)') 'ave', i,j,k,' = ',ave(i,j,k)
1144             !              endif
1145          enddo
1146       enddo
1147    enddo
1148
1149
1150    !    print*, "Beta11", ave(1,1,1)
1151    !    print*, "Alfa11", ave(1,1,2)
1152    !    print*, "Gama11", ave(1,1,3)
1153    !
1154    !    print*, "Beta12", ave(3,3,1)
1155    !    print*, "Beta12", ave(3,3,2)
1156    !    print*, "Beta12", ave(3,3,3)
1157    !
1158    !    print*, "Beta13", ave(5,5,3)
1159    !
1160    !    print*, "Beta21", ave(1,3,1)
1161    !    print*, "Beta22", ave(1,5,1)
1162    !    print*, "Beta23", ave(3,5,1)
1163    !
1164    !    print*, "Beta31", ave(3,1,1)
1165    !    print*, "Beta32", ave(5,1,1)
1166    !    print*, "Beta33", ave(5,5,1)
1167
1168
1169  end subroutine twissfctn
1170
1171  !_________________________________________________________________________________
1172  subroutine parametric_coord(y)   !  Computes <x_i x_j> assuming linearity and with parameters
1173    implicit none
1174    type(real_8) y(6)
1175
1176    e(:)=0
1177
1178    results(currentrow, kn_x)  = y(1).par.e
1179    results(currentrow, kn_px) = y(2).par.e
1180
1181    results(currentrow, kn_y)  = y(3).par.e
1182    results(currentrow, kn_py) = y(4).par.e
1183
1184    if (c_%npara_fpp == 5) then
1185       results(currentrow, kn_dp) = y(5).par.e
1186       results(currentrow, kn_t)  = zero
1187    else if(c_%npara_fpp == 6) then
1188       results(currentrow, kn_dp) = y(5).par.e
1189       results(currentrow, kn_t)  = y(6).par.e
1190    else
1191       results(currentrow, kn_dp) = zero
1192       results(currentrow, kn_t)  = zero
1193    endif
1194
1195  end subroutine parametric_coord
1196
1197  !_________________________________________________________________________________
1198  subroutine parametrictwiss(y)   !  Computes <x_i x_j> assuming linearity and with parameters
1199    implicit none
1200    type(real_8) y(6)
1201    real(dp) :: epsil=1e-12
1202    integer i,j,k
1203
1204    e=0
1205    do i=1,c_%nd2
1206       do j=i,c_%nd2
1207          do k=1,c_%nd
1208             ave(i,j,k)=zero
1209             e(k*2-1)=1
1210             ave(i,j,k)= ave(i,j,k) + (y(i)%t.par.e)*(y(j)%t.par.e)
1211             e(k*2-1)=0
1212             e(k*2)=1
1213             ave(i,j,k)=morph(ave(i,j,k))+ (y(i).par.e)*(y(j).par.e) !line * does the same, here taylor is morphed to polimorph,
1214             e(2*k)=0                                                !and above taylor component of polimorph is used explicitely
1215             ave(j,i,k)=ave(i,j,k)
1216
1217          enddo
1218       enddo
1219    enddo
1220   
1221    ave(1,2,1) = -ave(1,2,1)
1222    results(currentrow, alfa11) = ave(1,2,1)
1223    results(currentrow, beta11) = ave(1,1,1)
1224    results(currentrow, gama11) = ave(2,2,1)
1225
1226    ave(3,4,2) = -ave(3,4,2)
1227    results(currentrow, alfa22) = ave(3,4,2)
1228    results(currentrow, beta22) = ave(3,3,2)
1229    results(currentrow, gama22) = ave(4,4,2)
1230
1231    !----------------
1232
1233    ave(1,2,2) = -ave(1,2,2)
1234    results(currentrow, alfa12) = ave(1,2,2)
1235    results(currentrow, beta12) = ave(1,1,2)
1236    results(currentrow, gama12) = ave(2,2,2)
1237
1238    ave(3,4,1) = -ave(3,4,1)
1239    results(currentrow, alfa21) = ave(3,4,1)
1240    results(currentrow, beta21) = ave(3,3,1)
1241    results(currentrow, gama21) = ave(4,4,1)
1242
1243
1244    !----------------
1245
1246    if ( c_%nd == 3) then
1247       ave(1,2,3) = -ave(1,2,3)
1248       results(currentrow, alfa13) = ave(1,2,3) !-
1249       results(currentrow, beta13) = ave(1,1,3)  !-
1250       results(currentrow, gama13) = ave(2,2,3)  !-
1251
1252       ave(3,4,3) =-ave(3,4,3)
1253       results(currentrow, alfa23) = ave(3,4,3) !-
1254       results(currentrow, beta23) = ave(3,3,3)  !-
1255       results(currentrow, gama23) = ave(4,4,3)  !-
1256
1257       ave(5,6,3) = -ave(5,6,3)
1258       results(currentrow, alfa33) = ave(5,6,3) !-
1259       results(currentrow, beta33) = ave(6,6,3) !-
1260       results(currentrow, gama33) = ave(5,5,3)
1261
1262       ave(5,6,2) =-ave(5,6,2)
1263       results(currentrow, alfa32) = ave(5,6,2) !
1264       results(currentrow, beta32) = ave(6,6,2)  !-
1265       results(currentrow, gama32) = ave(5,5,2)  !-
1266
1267       ave(5,6,1) = -ave(5,6,1)
1268       results(currentrow, alfa31) = ave(5,6,1) !-
1269       results(currentrow, beta31) = ave(6,6,1)  !-
1270       results(currentrow, gama31) = ave(5,5,1)  !-
1271
1272    else
1273
1274       results(currentrow, alfa13) = zero
1275       results(currentrow, beta13) = zero
1276       results(currentrow, gama13) = zero
1277       results(currentrow, alfa23) = zero
1278       results(currentrow, beta23) = zero
1279       results(currentrow, gama23) = zero
1280       results(currentrow, alfa33) = zero
1281       results(currentrow, beta33) = zero
1282       results(currentrow, gama33) = zero
1283       results(currentrow, alfa32) = zero
1284       results(currentrow, beta32) = zero
1285       results(currentrow, gama32) = zero  !-
1286       results(currentrow, alfa31) = zero
1287       results(currentrow, beta31) = zero
1288       results(currentrow, gama31) = zero
1289
1290    endif
1291
1292    if( (c_%npara==5)       .or.  (c_%ndpt/=0) ) then
1293       !when there is no cavity it gives us dispersions
1294       do i=1,4
1295          ave(i,1,1)=(Y(i)%t.par.J5)
1296       enddo
1297    elseif (c_%nd2 == 6) then
1298       do i=1,4
1299          ave(i,1,1) =              (Y(i)%t.par.J5)*(Y(6)%t.par.J6)
1300          ave(i,1,1) = ave(i,1,1) + (Y(i)%t.par.J6)*(Y(5)%t.par.J5)
1301       enddo
1302    else
1303       do i=1,4
1304          ave(i,1,1)=zero
1305       enddo
1306    endif
1307
1308    results(currentrow, disp1) = ave(1,1,1)
1309    results(currentrow, disp2) = ave(2,1,1)
1310    results(currentrow, disp3) = ave(3,1,1)
1311    results(currentrow, disp4) = ave(4,1,1)
1312
1313
1314!!!!!!!!!!!!!!!!
1315    ! phase advance!
1316!!!!!!!!!!!!!!!!
1317
1318
1319    k = 2
1320    if(c_%nd2==6.and.c_%ndpt==0) k = 3
1321
1322    do i=1, k
1323       tx = Y(2*i -1).PAR.fo(2*i,:)
1324       ty = Y(2*i-1).PAR.fo(2*i-1,:)
1325       TEST=ATAN2(tx,ty) /TWOPI
1326
1327       IF(test<zero.AND.abs(test)>EPSIL)  TEST=TEST+one
1328       DPH=TEST-TESTOLD(i)
1329
1330       IF(dph<one.AND.abs(dph)>EPSIL) DPH=DPH+one
1331       IF(dph>half) DPH=DPH-one
1332
1333       PHASE(i)=PHASE(i)+DPH
1334       TESTOLD(i)=TEST
1335
1336    enddo
1337
1338
1339    results(currentrow, mu1) = phase(1)
1340    results(currentrow, mu2) = phase(2)
1341    results(currentrow, mu3) = phase(3)
1342
1343
1344
1345
1346    results(currentrow, beta11p) = zero
1347    results(currentrow, beta12p) = zero
1348    results(currentrow, beta13p) = zero
1349    results(currentrow, beta21p) = zero
1350    results(currentrow, beta22p) = zero
1351    results(currentrow, beta23p) = zero
1352    results(currentrow, beta31p) = zero
1353    results(currentrow, beta32p) = zero
1354    results(currentrow, beta33p) = zero
1355    results(currentrow, alfa11p) = zero
1356    results(currentrow, alfa12p) = zero
1357    results(currentrow, alfa13p) = zero
1358    results(currentrow, alfa21p) = zero
1359    results(currentrow, alfa22p) = zero
1360    results(currentrow, alfa23p) = zero
1361    results(currentrow, alfa31p) = zero
1362    results(currentrow, alfa32p) = zero
1363    results(currentrow, alfa33p) = zero
1364    results(currentrow, gama11p) = zero
1365    results(currentrow, gama12p) = zero
1366    results(currentrow, gama13p) = zero
1367    results(currentrow, gama21p) = zero
1368    results(currentrow, gama22p) = zero
1369    results(currentrow, gama23p) = zero
1370    results(currentrow, gama31p) = zero
1371    results(currentrow, gama32p) = zero
1372    results(currentrow, gama33p) = zero
1373
1374    results(currentrow, disp1p) = zero
1375    results(currentrow, disp2p) = zero
1376    results(currentrow, disp3p) = zero
1377    results(currentrow, disp4p) = zero
1378
1379    results(currentrow, disp1p2) = zero
1380    results(currentrow, disp2p2) = zero
1381    results(currentrow, disp3p2) = zero
1382    results(currentrow, disp4p2) = zero
1383
1384    results(currentrow, disp1p3) = zero
1385    results(currentrow, disp2p3) = zero
1386    results(currentrow, disp3p3) = zero
1387    results(currentrow, disp4p3) = zero
1388
1389   
1390!    print*, "parametrictwiss(",currentrow,",bet11):", beta12
1391!    call print(results(currentrow,beta12),6)
1392   
1393    !     print*, "Beta X"
1394    !     call daprint(ave(1,1,1),6)
1395    !
1396    !     print*, "Beta Y"
1397    !     call daprint(ave(3,3,2),6)
1398    !
1399    !     print*, "Beta Z"
1400    !     call daprint(ave(5,5,3),6)
1401
1402
1403
1404  end subroutine parametrictwiss
1405  !_________________________________________________________________________________
1406
1407  subroutine finishknobs
1408    implicit none
1409    integer i,j,k
1410
1411
1412    if (taylorsallocated == 0) then
1413       return
1414    endif
1415
1416    call kill(phase)
1417    call kill(testold)
1418    call kill(test)
1419    call kill(dph)
1420    call kill(tx, ty)
1421
1422    do i=1,c_%nd2
1423       do j=1,c_%nd2
1424          do k=1,c_%nd
1425             call kill(ave(i,j,k))
1426          enddo
1427       enddo
1428    enddo
1429
1430    taylorsallocated = 0
1431
1432  end subroutine finishknobs
1433  !_________________________________________________________________________________
1434
1435  subroutine resultswithknobs(n,name,y)
1436    implicit none
1437    integer              :: n !fibre number
1438    character(*)         :: name !fibre name
1439    type(real_8),target  :: y(6)!input 6 dimensional function (polynomial)
1440    !    real(kind(1d0))      :: coeff
1441    !    integer              :: i,ii,j,k !iterator
1442    !    type(taylor)         :: ave(6,6,3)
1443
1444    print*," resultswithknobs not yet implemented ",n,name
1445    call print(y(1),6)
1446
1447  end subroutine resultswithknobs
1448
1449  !_________________________________________________________________________________
1450  subroutine writeparresults(filenameIA)
1451    use twissafi
1452    implicit none
1453    integer   filenameIA(*)
1454    integer       :: mf !macro file descriptor
1455    character(48) :: filename
1456    character(48) :: fmat
1457    integer       :: i,j, nel, slen
1458    logical       :: fmt_ptc, fmt_tex
1459    integer, parameter         :: length=16
1460    character(length)          :: name
1461    integer                    :: get_string
1462    integer                    :: restart_sequ,advance_node
1463    character(48) charconv
1464   
1465    if(getdebug()>1) then
1466     print*," writeparresults "
1467     print*,"nv=",c_%nv
1468     print*,"nd2=",c_%nd2
1469     print*,"np=",c_%np
1470     print*,"ndpt=",c_%ndpt
1471     print*,"=>",c_%nv-c_%nd2-c_%np
1472    endif
1473   
1474    if (.not. ALLOCATED(results)) then
1475       call fort_warn("writeparresults","Array with parametric results is not present.")
1476       print*, "writeparresults tip: it might have been erased the ptc_end command."
1477       return
1478    endif
1479
1480    if (filenameIA(1) > 0) then
1481       filename = charconv(filenameIA)
1482       call kanalnummer(mf)
1483       open(unit=mf,file=filename)
1484    else
1485       mf = 6
1486    endif
1487
1488    fmt_ptc = .false.
1489    fmt_tex = .false.
1490
1491    i = get_string('ptc_printparametric ','format ',fmat)
1492    if (i > 0) then
1493       print*,"ptc_printparametric: format is ", fmat(1:i)
1494
1495       select case(fmat(1:i))
1496       case ('ptc')
1497          print*, "ptc_printparametric: Recognized PTC native format"
1498          fmt_ptc = .true.
1499       case ('tex')
1500          print*, "ptc_printparametric: Recognized LaTeX native format"
1501          fmt_tex = .true.
1502       case default
1503          print*, "ptc_printparametric: Format not recognized - using default PTC"
1504          fmt_ptc = .true.
1505       end select
1506    else
1507       print*, "Got empty Format - using default PTC"
1508       fmt_ptc = .true.
1509    endif
1510
1511
1512    print*,"ptc_printparametric : currentrow is ", currentrow
1513
1514    nel = restart_sequ()
1515
1516    do i=1,currentrow-1
1517
1518       call node_name(name,length)
1519
1520       write(mf,*) "Magnet ", i," ",name(1:length)
1521
1522
1523       do j=1,ntwisses
1524          !        print*, "Writing i,j",i,j
1525          write(mf,*) twissnames(j)
1526          if (fmt_ptc)  call print(results(i,j),mf)
1527          if (fmt_tex)  call printpareq(results(i,j),mf)
1528          write(mf,*) " "
1529       enddo
1530
1531
1532       do j=1,npushes
1533          if (pushes(j)%index < 1 ) cycle
1534          !        print*, "Writing i, j->index",i,j,pushes(j)%index
1535          slen = len_trim(pushes(j)%colname) - 1
1536          write(mf,*) pushes(j)%colname(1:slen)
1537          if (fmt_ptc)  call print(results(i,pushes(j)%index),mf)
1538          if (fmt_tex)  call printpareq(results(i,pushes(j)%index),mf)
1539          write(mf,*) " "
1540       enddo
1541       write(mf,*) "======================================="
1542
1543       nel = advance_node()
1544
1545    enddo
1546
1547    if (mf /= 6) close(mf)
1548
1549  end subroutine writeparresults
1550  !_________________________________________________________________________________
1551  function getparname(n)
1552    implicit none
1553    character*(3)  :: getparname
1554    integer        :: n
1555
1556    write(getparname,'(A1,i2.2)') "p",n - c_%npara
1557
1558  end function getparname
1559  !_________________________________________________________________________________
1560
1561  subroutine setparvalue(n,v)
1562    implicit none
1563    integer :: n
1564    real :: v
1565
1566    if (.not. allocated(parvals)) then
1567       call fort_warn("setparvalue",&
1568            "array with parameter values is not allocated")
1569       return
1570    endif
1571
1572    if ( (n < 1) .and. (n > ubound(parvals,1)) ) then
1573       call fort_warn("setparvalue","Array index out of range")
1574    endif
1575
1576    parvals(n) = v
1577
1578  end subroutine setparvalue
1579  !____________________________________________________________________________________________
1580
1581  subroutine setknobvalue(fibrenameIA)
1582    use twissafi
1583    implicit none
1584    integer                    :: fibrenameIA(*)
1585    character(48)              :: fibrename
1586    integer                    :: i
1587    integer                    :: kn, ks, par
1588    real                       :: v
1589    real(kind(1d0))            :: get_value
1590    logical(lp)                :: refreshtables
1591    character(48) charconv
1592    par = -1
1593    fibrename = charconv(fibrenameIA)
1594
1595    !    print *,"setknobvalue: fibrename is ", fibrename
1596
1597
1598    kn = get_value('ptc_setknobvalue ','kn ') + 1 !madx numerates from 0
1599    ks = get_value('ptc_setknobvalue ','ks ') + 1
1600
1601    if ( (kn>0) .and. (ks>0) ) then
1602       call fort_warn("setknobvalue","Both kn and ks can not be specified together");
1603       return
1604    endif
1605
1606    if ( (kn<=0) .and. (ks<=0) ) then
1607
1608       select case(fibrename(1:6))
1609       case ('BETA11')
1610          par = knobi%beta(1)
1611       case ('BETA22')
1612          par = knobi%beta(2)
1613       case ('BETA33')
1614          par = knobi%beta(3)
1615       case ('ALFA11')
1616          par = knobi%ALFA(1)
1617       case ('ALFA22')
1618          par = knobi%ALFA(2)
1619       case ('ALFA33')
1620          par = knobi%alfa(3)
1621       case ('DISP1 ')
1622          par = knobi%DISPersion(1)
1623       case ('DISP2 ')
1624          par = knobi%DISPersion(2)
1625       case ('DISP3 ')
1626          par = knobi%DISPersion(3)
1627       case ('DISP4 ')
1628          par = knobi%dispersion(4)
1629
1630       case default
1631          print*, "Name of initial condition parameter ",fibrename," not recognized"
1632          return
1633       end select
1634
1635       par = par + getnknobsm() !init cond starts as parameters after field components
1636    else
1637       do i=1, npolblocks
1638          if ( polblocks(i)%name == fibrename(1:nlp)) then
1639
1640             if ( kn>0 ) then
1641                par = polblocks(i)%ibn(kn)
1642             elseif ( ks>0 ) then
1643                par = polblocks(i)%ian(ks)
1644             endif
1645
1646             exit
1647          endif
1648       enddo
1649
1650       if (par < 0) then
1651          call fort_warn("setknobvalue","There is no knob defined on such element");
1652          return
1653       endif
1654
1655    endif
1656
1657
1658    v = get_value('ptc_setknobvalue ','value ')
1659
1660    if (getdebug() > 1) then
1661       print*, "Setting parameter ",par,"(el=",fibrename(1:16),", kn=",kn,", ks=",ks," ) to ", v
1662    endif
1663
1664    if (getdebug() > 1) then
1665        print*, "Setting par ", par, " to ", v, fibrename(1:16)
1666    end if
1667    call setparvalue(par, v)
1668
1669    refreshtables = get_value('ptc_setknobvalue ','refreshtables ') .ne. 0
1670
1671    if (refreshtables) then
1672       call filltables()
1673    endif
1674
1675
1676  end subroutine setknobvalue
1677  !____________________________________________________________________________________________
1678
1679  subroutine filltables
1680    implicit none
1681    call filltwisstable()
1682    !    call cleartables()
1683    call fillusertables()
1684  end subroutine filltables
1685  !____________________________________________________________________________________________
1686
1687  function getlengthat(n)
1688    implicit none
1689    real(dp)         :: getlengthat
1690    integer          :: n
1691
1692    print*, "getlengthat, n is ", n
1693
1694    if (.not. ALLOCATED(spos)) then
1695       return
1696    endif
1697
1698    if ( (n < 1) .and. (n > ubound(spos,1)) ) then
1699       call fort_warn("getlengthat","position out of range")
1700       getlengthat = -one
1701    endif
1702
1703    print*, "getlengthat, spos at n is ", spos(n)
1704
1705    getlengthat = spos(n)
1706
1707  end function getlengthat
1708  !_________________________________________________________________________________
1709
1710
1711  subroutine printpareq(ut,iunit)
1712    implicit none
1713    type(universal_taylor)  :: ut
1714    integer                 :: iunit, eqlen
1715     
1716    if( .not. associated (ut%n)) then
1717       call fort_warn("printpareq", "this universal taylor is void")
1718       write(iunit,'(A)') "this universal taylor is void"
1719       return
1720    endif
1721   
1722    if(ut%nv /= c_%nv) then
1723       call fort_warn("printpareq",&
1724            "number of variables of this universal taylor is different from currnet TPSA")
1725       call print(ut,6)
1726       print*,"nv=",c_%nv
1727       print*,"nd2=",c_%nd2
1728       print*,"np=",c_%np
1729       print*,"ndpt=",c_%ndpt
1730       print*,"=>",c_%nv-c_%nd2-c_%np
1731       return
1732    endif
1733
1734
1735    call getpareq(ut,textbuffer)
1736    eqlen = len_trim(textbuffer) + 1;
1737    write(iunit,'(A)')  textbuffer(1:eqlen)
1738   
1739
1740  end subroutine printpareq
1741  !_________________________________________________________________________________
1742
1743  subroutine getpareq(ut,string)
1744    implicit none
1745    type(universal_taylor) :: ut
1746    integer                :: i,ii
1747    integer                :: cpos, last
1748    character              :: sign
1749    character*(20)         :: parname
1750    character(*)           :: string
1751
1752    !    print*,"getpareq"
1753
1754    if (.not. associated(ut%n)) then
1755       call fort_warn("getpareq",&
1756            "provided taylor is not allocated")
1757       write(string,'(A)') ' 0 '
1758       return
1759    endif
1760
1761
1762    if(ut%nv /= c_%nv) then
1763       call fort_warn("getpareq",&
1764            "number of variables of this universal taylor is different from currnet TPSA")
1765       return
1766    endif
1767
1768
1769    if(ut%n == 0) then
1770       write(string,'(A)') ' 0 '
1771       if (getdebug() > 3) then
1772          print*,"no coefficients in the taylor"
1773       endif
1774       return
1775    endif
1776
1777    cpos = 1
1778    last = len(string)
1779    sign = ' '
1780
1781    if (getdebug() > 3) then
1782       print*,"There is ", ut%n, " coefficients "
1783    endif
1784
1785    do i = 1,ut%n
1786
1787       if ( ut%c(i) < zero ) then
1788          sign = ' '
1789       endif
1790
1791       write(string(cpos:last),'(A1,A1,G21.14)') " ",sign,ut%c(i);
1792       cpos = len_trim(string) + 1;
1793
1794       if ( ( cpos + 200) > last) then
1795          call fort_warn("routine madx_ptc_knobs.f90::getpareq",&
1796               "Buffer for taylor equation is too small! Bailing out!")
1797          stop;
1798       endif
1799
1800       do ii = 1, ut%nv
1801          !         print*, "cpos=",cpos," last=",last
1802          if (ut%j(i,ii) /= 0) then
1803             write(string(cpos:last),'(A1)') "*"; cpos = len_trim(string) + 1;
1804             parname = getparname(ii)
1805             write(string(cpos:last),'(a)') parname(1:LEN_TRIM(parname)); cpos = len_trim(string) + 1;
1806
1807             if (ut%j(i,ii) /= 1) then
1808                write(string(cpos:last),'(a1)')  "^"; cpos = len_trim(string) + 1;
1809                write(parname,'(i3)')  ut%j(i,ii)
1810                parname = ADJUSTL(parname)
1811                write(string(cpos:last),'(a)')  parname(1:LEN_TRIM(parname)); cpos = len_trim(string) + 1;
1812             endif
1813          endif
1814
1815       enddo
1816
1817       write(string(cpos:last),'(A1)') " "; cpos = len_trim(string) + 1;
1818       sign = '+'
1819    enddo
1820
1821
1822    !    print*, string
1823  end subroutine getpareq
1824
1825
1826end module madx_ptc_knobs_module
1827
1828
1829
1830! backup code for the results verification
1831!     call twissfctn(scv,ave)
1832!     print*, "==============================================================="
1833!     print*, "==============================================================="
1834!     print*, "==============================================================="
1835!     print*, "==============================================================="
1836!     do ii=1,c_%nd2
1837!       do jj=ii,c_%nd2
1838!         do kk=1,c_%nd
1839!
1840!          if (ave(ii,jj,kk) /= zero) then
1841!            write(6,'(a3,i1,i1,i1,a3,f16.8)') 'ave', ii,jj,kk,' = ',ave(ii,jj,kk)
1842!            do ll=1,3
1843!              do mm=1,3
1844!
1845!                v = abs(ave(ii,jj,kk) - tw%beta(ll,mm))
1846!                if ( v < 1e-16_dp) then
1847!                  print*, "v(",ll,",",mm,")=",v
1848!                  write(6,'(a3,i1,i1,i1,a7,i1,i1)') 'ave', ii,jj,kk,' = beta',ll,mm
1849!                endif
1850!
1851!                v = abs(ave(ii,jj,kk) + tw%alfa(ll,mm))
1852!                if ( v < 1e-16_dp) then
1853!                  print*, "v(",ll,",",mm,")=",v
1854!                  write(6,'(a3,i1,i1,i1,a7,i1,i1)') 'ave', ii,jj,kk,' = alfa',ll,mm
1855!                endif
1856!
1857!
1858!                v = abs(ave(ii,jj,kk) - tw%gama(ll,mm))
1859!                if ( v < 1e-16_dp) then
1860!                  print*, "v(",ll,",",mm,")=",v
1861!                  write(6,'(a3,i1,i1,i1,a7,i1,i1)') 'ave', ii,jj,kk,' = gama',ll,mm
1862!                endif
1863!
1864!              enddo
1865!
1866!              v = abs(ave(ii,jj,kk) - tw%mu(ll))
1867!              if ( v < 1e-16_dp) then
1868!                print*, "v(",ll,",",mm,")=",v
1869!                write(6,'(a3,i1,i1,i1,a7,i1)') 'ave', ii,jj,kk,' = mu',ll
1870!              endif
1871!
1872!              v = abs(ave(ii,jj,kk) - tw%disp(ll))
1873!              if ( v < 1e-16_dp) then
1874!                print*, "v(",ll,",",mm,")=",v
1875!                write(6,'(a3,i1,i1,i1,a7,i1)') 'ave', ii,jj,kk,' = disp',ll
1876!              endif
1877!
1878!              v = abs(ave(ii,jj,kk) - tw%disp(ll+1))
1879!              if ( v < 1e-16_dp) then
1880!                print*, "v(",ll,",",mm,")=",v
1881!                write(6,'(a3,i1,i1,i1,a7,i1)') 'ave', ii,jj,kk,' = disp',ll+1
1882!              endif
1883!
1884!
1885!            enddo
1886!
1887!          endif
1888!
1889!       enddo
1890!    enddo
1891! enddo
1892!
1893!     print*, "==============================================================="
1894!     print*, "==============================================================="
1895!     print*, "==============================================================="
1896!     print*, "==============================================================="
1897
1898
1899!____________________________________________________________________________________________
1900!____________________________________________________________________________________________
1901!____________________________________________________________________________________________
1902!____________________________________________________________________________________________
1903
1904
1905function w_ptc_getnknobs()
1906  use madx_ptc_knobs_module
1907  implicit none
1908  integer w_ptc_getnknobs
1909
1910  w_ptc_getnknobs = getnknobsall()
1911
1912end function w_ptc_getnknobs
1913!____________________________________________________________________________________________
1914
1915
1916function w_ptc_getlengthat(n)
1917  use madx_ptc_knobs_module
1918  implicit none
1919  real(kind(1d0)) :: w_ptc_getlengthat
1920  integer         :: n
1921
1922  w_ptc_getlengthat = getlengthat(n)
1923
1924end function w_ptc_getlengthat
1925!____________________________________________________________________________________________
1926
1927
1928
1929subroutine w_ptc_getfctnsnames()
1930  use madx_ptc_knobs_module
1931  implicit none
1932
1933  call getfctnsnames()
1934
1935end subroutine w_ptc_getfctnsnames
1936!____________________________________________________________________________________________
1937
1938subroutine w_ptc_getknobsnames()
1939  use madx_ptc_knobs_module
1940  implicit none
1941
1942  call getknobsnames()
1943
1944end subroutine w_ptc_getknobsnames
1945!____________________________________________________________________________________________
1946
1947
1948subroutine w_ptc_getfunctionat(e,n)
1949  use madx_ptc_knobs_module
1950  implicit none
1951  integer e,n
1952
1953  call getfunctionat(e,n)
1954
1955end subroutine w_ptc_getfunctionat
1956
1957!____________________________________________________________________________________________
1958
1959
1960function w_ptc_getfunctionvalueat(n,e)
1961  use madx_ptc_knobs_module
1962  implicit none
1963  real(kind(1d0)) :: w_ptc_getfunctionvalueat
1964  integer e,n
1965
1966
1967  w_ptc_getfunctionvalueat = getfunctionvalueat(n,e)
1968
1969end function w_ptc_getfunctionvalueat
1970
1971!____________________________________________________________________________________________
1972
1973
1974subroutine w_ptc_rviewer()
1975  implicit none
1976  integer :: rviewer
1977  integer :: res
1978  res =  rviewer()
1979end subroutine w_ptc_rviewer
1980
1981!____________________________________________________________________________________________
1982
1983
1984subroutine w_ptc_setparvalue(n,v)
1985  use madx_ptc_knobs_module
1986  implicit none
1987  integer :: n
1988  real    :: v
1989
1990  call setparvalue(n,v)
1991
1992end subroutine w_ptc_setparvalue
Note: See TracBrowser for help on using the repository browser.