1 | module 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 | |
---|
126 | contains |
---|
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 | |
---|
1826 | end 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 | |
---|
1905 | function w_ptc_getnknobs() |
---|
1906 | use madx_ptc_knobs_module |
---|
1907 | implicit none |
---|
1908 | integer w_ptc_getnknobs |
---|
1909 | |
---|
1910 | w_ptc_getnknobs = getnknobsall() |
---|
1911 | |
---|
1912 | end function w_ptc_getnknobs |
---|
1913 | !____________________________________________________________________________________________ |
---|
1914 | |
---|
1915 | |
---|
1916 | function 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 | |
---|
1924 | end function w_ptc_getlengthat |
---|
1925 | !____________________________________________________________________________________________ |
---|
1926 | |
---|
1927 | |
---|
1928 | |
---|
1929 | subroutine w_ptc_getfctnsnames() |
---|
1930 | use madx_ptc_knobs_module |
---|
1931 | implicit none |
---|
1932 | |
---|
1933 | call getfctnsnames() |
---|
1934 | |
---|
1935 | end subroutine w_ptc_getfctnsnames |
---|
1936 | !____________________________________________________________________________________________ |
---|
1937 | |
---|
1938 | subroutine w_ptc_getknobsnames() |
---|
1939 | use madx_ptc_knobs_module |
---|
1940 | implicit none |
---|
1941 | |
---|
1942 | call getknobsnames() |
---|
1943 | |
---|
1944 | end subroutine w_ptc_getknobsnames |
---|
1945 | !____________________________________________________________________________________________ |
---|
1946 | |
---|
1947 | |
---|
1948 | subroutine 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 | |
---|
1955 | end subroutine w_ptc_getfunctionat |
---|
1956 | |
---|
1957 | !____________________________________________________________________________________________ |
---|
1958 | |
---|
1959 | |
---|
1960 | function 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 | |
---|
1969 | end function w_ptc_getfunctionvalueat |
---|
1970 | |
---|
1971 | !____________________________________________________________________________________________ |
---|
1972 | |
---|
1973 | |
---|
1974 | subroutine w_ptc_rviewer() |
---|
1975 | implicit none |
---|
1976 | integer :: rviewer |
---|
1977 | integer :: res |
---|
1978 | res = rviewer() |
---|
1979 | end subroutine w_ptc_rviewer |
---|
1980 | |
---|
1981 | !____________________________________________________________________________________________ |
---|
1982 | |
---|
1983 | |
---|
1984 | subroutine 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 | |
---|
1992 | end subroutine w_ptc_setparvalue |
---|