1 | module madx_ptc_distrib_module |
---|
2 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1 |
---|
3 | ! madx_ptc_distrib module |
---|
4 | ! Piotr K. Skowronski (CERN) |
---|
5 | ! |
---|
6 | ! This module contains service for trackin |
---|
7 | |
---|
8 | use madx_keywords |
---|
9 | use madx_ptc_module |
---|
10 | USE madx_ptc_knobs_module |
---|
11 | use madx_ptc_intstate_module, only : getdebug |
---|
12 | implicit none |
---|
13 | |
---|
14 | include "madx_ptc_distrib.inc" |
---|
15 | save |
---|
16 | private |
---|
17 | |
---|
18 | !============================================================================================ |
---|
19 | ! PUBLIC INTERFACE |
---|
20 | public :: momfirstinit |
---|
21 | public :: aremomentson |
---|
22 | public :: ptc_moments |
---|
23 | public :: putmoments |
---|
24 | public :: initmoments |
---|
25 | public :: allocmoments |
---|
26 | public :: killmoments |
---|
27 | public :: getdistrtype |
---|
28 | public :: setemittances |
---|
29 | public :: setsigma |
---|
30 | public :: printsigmas |
---|
31 | public :: addmoment |
---|
32 | public :: getnmoments |
---|
33 | public :: getmomentstabcol |
---|
34 | |
---|
35 | !============================================================================================ |
---|
36 | ! PRIVATE |
---|
37 | ! data structures |
---|
38 | integer :: firstinitdone = 0 |
---|
39 | real(dp), pointer :: normmoments(:,:,:) |
---|
40 | real(dp) :: sigmas(6) |
---|
41 | integer :: distributiontype(3) |
---|
42 | |
---|
43 | integer, parameter :: maxnmoments=100 |
---|
44 | type (momentdef), private :: moments(maxnmoments) |
---|
45 | integer :: nmoments = 0 |
---|
46 | |
---|
47 | type(gmap) :: gmapa !da map needed to calculation initialized with sigmas |
---|
48 | type(taylor) :: function_to_average ! taylor used to calculate averages |
---|
49 | |
---|
50 | !============================================================================================ |
---|
51 | ! PRIVATE |
---|
52 | ! routines |
---|
53 | |
---|
54 | public :: makegaus |
---|
55 | public :: makeflat5 |
---|
56 | public :: makeflat56 |
---|
57 | private :: filter |
---|
58 | |
---|
59 | |
---|
60 | contains |
---|
61 | !____________________________________________________________________________________________ |
---|
62 | |
---|
63 | !____________________________________________________________________________________________ |
---|
64 | |
---|
65 | |
---|
66 | subroutine addmoment(x,px,y,py,t,dp,tableIA, columnIA, parametric ) |
---|
67 | use twissafi |
---|
68 | implicit none |
---|
69 | integer :: x,px,y,py,dp,t |
---|
70 | integer :: columnIA(*) |
---|
71 | integer :: tableIA(*) |
---|
72 | integer :: parametric |
---|
73 | integer, parameter :: zeroasciicode = IACHAR ( '0' ) |
---|
74 | character(48) charconv |
---|
75 | ! name of the column corresponds to MADX nomenclature (5 col is dp/p) |
---|
76 | ! iarray corresponds to |
---|
77 | |
---|
78 | nmoments = nmoments + 1 |
---|
79 | |
---|
80 | moments(nmoments)%iarray(1) = x |
---|
81 | moments(nmoments)%iarray(2) = px |
---|
82 | moments(nmoments)%iarray(3) = y |
---|
83 | moments(nmoments)%iarray(4) = py |
---|
84 | moments(nmoments)%iarray(5) = dp |
---|
85 | moments(nmoments)%iarray(6) = t |
---|
86 | |
---|
87 | moments(nmoments)%table = charconv(tableIA) |
---|
88 | moments(nmoments)%column = charconv(columnIA) |
---|
89 | |
---|
90 | moments(nmoments)%table(tableIA(1)+1:tableIA(1)+1)=achar(0) |
---|
91 | moments(nmoments)%column(columnIA(1)+1:columnIA(1)+1)=achar(0) |
---|
92 | |
---|
93 | if (parametric /= 0) then |
---|
94 | print*,"To be made as parametric variable" |
---|
95 | moments(nmoments)%index = nmoments !for the time being it is the same as the index |
---|
96 | else |
---|
97 | moments(nmoments)%index = 0 |
---|
98 | endif |
---|
99 | |
---|
100 | |
---|
101 | if (getdebug()>0) then |
---|
102 | print *,"addmoment : <", moments(nmoments)%iarray(1:6) ,& |
---|
103 | & ">,<", moments(nmoments)%column,& |
---|
104 | & ">,<", moments(nmoments)%table, ">)" |
---|
105 | endif |
---|
106 | |
---|
107 | |
---|
108 | |
---|
109 | end subroutine addmoment |
---|
110 | |
---|
111 | !____________________________________________________________________________________________ |
---|
112 | |
---|
113 | integer function getnmoments() |
---|
114 | implicit none |
---|
115 | getnmoments = nmoments |
---|
116 | end function getnmoments |
---|
117 | !____________________________________________________________________________________________ |
---|
118 | |
---|
119 | subroutine getmomentstabcol(n, tabn, coln) |
---|
120 | implicit none |
---|
121 | integer :: n |
---|
122 | character(20) :: tabn |
---|
123 | character(17) :: coln |
---|
124 | |
---|
125 | if ( (n < 1) .and. (n > nmoments ) ) then |
---|
126 | |
---|
127 | tabn(1:1) = achar(0) |
---|
128 | coln(1:1) = achar(0) |
---|
129 | call fort_warn("getmomentstabcol","index out of range") |
---|
130 | return |
---|
131 | endif |
---|
132 | |
---|
133 | |
---|
134 | tabn = moments(n)%table |
---|
135 | coln = moments(n)%column |
---|
136 | |
---|
137 | |
---|
138 | end subroutine getmomentstabcol |
---|
139 | !____________________________________________________________________________________________ |
---|
140 | |
---|
141 | subroutine momfirstinit() |
---|
142 | implicit none |
---|
143 | if (firstinitdone == 1) then |
---|
144 | return |
---|
145 | endif |
---|
146 | |
---|
147 | nullify(normmoments) |
---|
148 | firstinitdone = 1 |
---|
149 | |
---|
150 | end subroutine momfirstinit |
---|
151 | |
---|
152 | !____________________________________________________________________________________________ |
---|
153 | |
---|
154 | logical(lp) function aremomentson() |
---|
155 | implicit none |
---|
156 | |
---|
157 | if ( associated(normmoments) ) then |
---|
158 | aremomentson = my_true |
---|
159 | else |
---|
160 | aremomentson = my_false |
---|
161 | endif |
---|
162 | |
---|
163 | end function aremomentson |
---|
164 | !_________________________________________________________________ |
---|
165 | |
---|
166 | subroutine ptc_moments(order) |
---|
167 | implicit none |
---|
168 | integer :: order,mynd2,npara,nda |
---|
169 | integer :: i,ii,iii |
---|
170 | type(real_8) :: y(6) |
---|
171 | integer :: restart_sequ,advance_node |
---|
172 | |
---|
173 | if (nmoments < 1) then |
---|
174 | call fort_info("ptc_moments","No moments specified for calculation.") |
---|
175 | return |
---|
176 | endif |
---|
177 | |
---|
178 | if (mapsorder < 1) then |
---|
179 | call seterrorflag(1,"ptc_moments",& |
---|
180 | "Maps are not available. Did you run ptc_twiss with savemaps=true ?") |
---|
181 | return |
---|
182 | endif |
---|
183 | |
---|
184 | if (.not. associated(maps)) then |
---|
185 | return |
---|
186 | endif |
---|
187 | |
---|
188 | if (mapsicase == 5) then |
---|
189 | call fort_warn("ptc_moments","For the time being the calculation of moments is not available in 5D case.") |
---|
190 | return |
---|
191 | endif |
---|
192 | |
---|
193 | if ((mapsicase == 5) .and. (sigmas(5) .le. 0)) then |
---|
194 | call fort_warn("ptc_moments","Spread in dp/p in undefined and it won't be taken taken to the account") |
---|
195 | print*,"In 5D case you have to specify either" |
---|
196 | print*," - SIGE in the BEAM command or" |
---|
197 | print*," - ET in the BEAM command AND BETZ with in ptc_twiss" |
---|
198 | endif |
---|
199 | |
---|
200 | if ((mapsicase == 6) .and. (sigmas(5) .le. 0)) then |
---|
201 | call fort_warn("ptc_moments","Spread in dp/p in undefined and it won't be taken taken to the account") |
---|
202 | print*,"In 6D case you have to specify longitudinal emittance SIGE in the BEAM command" |
---|
203 | endif |
---|
204 | |
---|
205 | |
---|
206 | call initmoments() |
---|
207 | call makemomentstables(); |
---|
208 | |
---|
209 | nda = getnknobsall() !defined in madx_ptc_knobs |
---|
210 | |
---|
211 | !print*, "In moments order ", order |
---|
212 | mynd2 = 0 |
---|
213 | npara = 0 |
---|
214 | call init(default,order,nda,BERZ,mynd2,npara) |
---|
215 | |
---|
216 | call allocmoments() |
---|
217 | |
---|
218 | call alloc(y) |
---|
219 | |
---|
220 | iii=restart_sequ() |
---|
221 | |
---|
222 | do i=lbound(maps,1),ubound(maps,1) |
---|
223 | |
---|
224 | ! if (i == MY_RING%n) then |
---|
225 | ! call ptc_setdebuglevel(1) |
---|
226 | ! endif |
---|
227 | |
---|
228 | do ii=1,6 |
---|
229 | y(ii) = maps(i)%unimap(ii) |
---|
230 | enddo |
---|
231 | |
---|
232 | call putmoments(i,maps(i)%name,maps(i)%s,y) |
---|
233 | |
---|
234 | iii=advance_node() |
---|
235 | |
---|
236 | enddo |
---|
237 | 100 continue |
---|
238 | |
---|
239 | call ptc_setdebuglevel(0) |
---|
240 | |
---|
241 | call kill(y) |
---|
242 | call killmoments() |
---|
243 | |
---|
244 | end subroutine ptc_moments |
---|
245 | !____________________________________________________________________________________________ |
---|
246 | |
---|
247 | |
---|
248 | subroutine putmoments(n,name,s,y) |
---|
249 | implicit none |
---|
250 | integer :: n !fibre number |
---|
251 | character(*) :: name !fibre name |
---|
252 | real(dp) :: s !position along the orbit |
---|
253 | type(real_8),target :: y(6)!input 6 dimensional function (polynomial) : Full MAP: A*YC*A_1 |
---|
254 | real(kind(1d0)) :: v |
---|
255 | logical :: set |
---|
256 | integer :: i,j,k,e(6) |
---|
257 | integer :: debug |
---|
258 | ! integer :: last |
---|
259 | ! integer :: index !standarf function |
---|
260 | |
---|
261 | ! last = index(name,'$END') |
---|
262 | ! if ( last == 0) then |
---|
263 | debug = getdebug() |
---|
264 | ! else |
---|
265 | ! debug = 10; |
---|
266 | ! endif |
---|
267 | |
---|
268 | if ( .not. associated(normmoments) ) then |
---|
269 | return |
---|
270 | endif |
---|
271 | |
---|
272 | if (debug > 3) then |
---|
273 | |
---|
274 | print*, "#################################################" |
---|
275 | print*, "#################################################" |
---|
276 | print*, "#################################################" |
---|
277 | print*, "Moments for fibre ", n,name," at ", s,"m" |
---|
278 | print*, "ND2=",c_%nd2 |
---|
279 | print*, "NPARA=",c_%npara |
---|
280 | |
---|
281 | print*, "Function 1" |
---|
282 | call print(y(1),6) |
---|
283 | print*, "" |
---|
284 | print*, "Function 5" |
---|
285 | call print(y(5),6) |
---|
286 | print*, "" |
---|
287 | print*, "Function 6" |
---|
288 | call print(y(6),6) |
---|
289 | print*, "" |
---|
290 | |
---|
291 | print*, "GMap" |
---|
292 | call print(gmapa,6) |
---|
293 | print*, "" |
---|
294 | |
---|
295 | endif |
---|
296 | |
---|
297 | !--moments--! |
---|
298 | do i=1, nmoments |
---|
299 | |
---|
300 | function_to_average = zero |
---|
301 | set = .false. |
---|
302 | |
---|
303 | do j=1,c_%npara |
---|
304 | |
---|
305 | |
---|
306 | if (debug .gt. 3) write(*,'(a6,i1,a6,i1,a6,i1)',ADVANCE='NO') "nmom=",i," ndim=",j," pow=",moments(i)%iarray(j) |
---|
307 | do k = 1, moments(i)%iarray(j) |
---|
308 | if (set) then |
---|
309 | function_to_average = function_to_average*y(j)%t |
---|
310 | if (debug .gt. 3) write(*,'(a1)',ADVANCE='NO') "*" |
---|
311 | else |
---|
312 | function_to_average = y(j)%t |
---|
313 | set = .true. |
---|
314 | if (debug .gt. 3) write(*,'(a1)',ADVANCE='NO') "|" |
---|
315 | endif |
---|
316 | enddo |
---|
317 | if (debug .gt. 3) write(*,*) "->" |
---|
318 | |
---|
319 | enddo |
---|
320 | |
---|
321 | ! function_to_average=y(1)*y(1) ! just a function (taylor series) |
---|
322 | |
---|
323 | ! if (debug > 3) then |
---|
324 | ! print*, "function_to_average" |
---|
325 | ! call print(function_to_average,6) |
---|
326 | ! endif |
---|
327 | |
---|
328 | if (debug > 5) then |
---|
329 | print*, "Function to average" |
---|
330 | call print(function_to_average,6) |
---|
331 | endif |
---|
332 | |
---|
333 | |
---|
334 | call cfu(function_to_average,filter,function_to_average) !cycling i.e. put the form factor to the function |
---|
335 | |
---|
336 | if (debug > 4) then |
---|
337 | print*, "After cfu" |
---|
338 | call print(function_to_average,6) |
---|
339 | endif |
---|
340 | |
---|
341 | function_to_average=function_to_average.o.gmapa ! replaces x px y py ... by sigma1, sigma2, etc |
---|
342 | |
---|
343 | if (debug > 3) then |
---|
344 | print*, "Averaging (gmapped)" |
---|
345 | call print(function_to_average,6) |
---|
346 | |
---|
347 | endif |
---|
348 | |
---|
349 | v = function_to_average.sub.0 |
---|
350 | |
---|
351 | if (c_%npara == 5) then |
---|
352 | if (debug > 3) then |
---|
353 | print*, v |
---|
354 | endif |
---|
355 | e = 0 |
---|
356 | do k=1,c_%no |
---|
357 | e(5) = k |
---|
358 | if (debug > 3) then |
---|
359 | print*, "s^",k,"=",(sigmas(5)**(k)) |
---|
360 | print*, "f = ", (function_to_average.sub.e) |
---|
361 | print*, (function_to_average.sub.e) * (sigmas(5)**(2*k)) |
---|
362 | endif |
---|
363 | v = v + (function_to_average.sub.e) * (sigmas(5)**(k)) !was to 2*k |
---|
364 | if (debug > 9) then |
---|
365 | print*, v |
---|
366 | endif |
---|
367 | enddo |
---|
368 | endif |
---|
369 | |
---|
370 | |
---|
371 | call double_to_table_curr(moments(i)%table,moments(i)%column,v) |
---|
372 | |
---|
373 | if (debug > 2) then |
---|
374 | print*, "Final ",i," = ", v |
---|
375 | print*,moments(i)%iarray |
---|
376 | call print(function_to_average,6) |
---|
377 | print*,"######################################################################################" |
---|
378 | endif |
---|
379 | |
---|
380 | |
---|
381 | enddo |
---|
382 | |
---|
383 | call augmentcountmomtabs(s) |
---|
384 | |
---|
385 | |
---|
386 | end subroutine putmoments |
---|
387 | |
---|
388 | !_________________________________________________________________________________ |
---|
389 | |
---|
390 | subroutine setemittances(emix,emiy,emiz) |
---|
391 | implicit none |
---|
392 | real(dp) :: emix |
---|
393 | real(dp) :: emiy |
---|
394 | real(dp) :: emiz |
---|
395 | |
---|
396 | if (emix .lt. zero) then |
---|
397 | print *, "X Emittance is less then 0" |
---|
398 | return |
---|
399 | endif |
---|
400 | |
---|
401 | if (emiy .lt. zero) then |
---|
402 | print *, "Y Emittance is less then 0" |
---|
403 | return |
---|
404 | endif |
---|
405 | |
---|
406 | if (emiz .lt. zero) then |
---|
407 | print *, "Z Emittance is less then 0" |
---|
408 | return |
---|
409 | endif |
---|
410 | |
---|
411 | |
---|
412 | if (getdebug() > 1) then |
---|
413 | print*, "Setting Sigmas (Emittances)" |
---|
414 | endif |
---|
415 | |
---|
416 | sigmas(1) = sqrt(emix) |
---|
417 | sigmas(2) = sigmas(1) |
---|
418 | sigmas(3) = sqrt(emiy) |
---|
419 | sigmas(4) = sigmas(3) |
---|
420 | |
---|
421 | sigmas(5) = sqrt(emiz) |
---|
422 | sigmas(6) = sigmas(5) |
---|
423 | |
---|
424 | if (getdebug() > 1) then |
---|
425 | print *, "Current sigmas setemittances ", sigmas |
---|
426 | endif |
---|
427 | |
---|
428 | |
---|
429 | |
---|
430 | end subroutine setemittances |
---|
431 | !_________________________________________________________________________________ |
---|
432 | |
---|
433 | |
---|
434 | subroutine setsigma(ndim,sig) |
---|
435 | implicit none |
---|
436 | integer :: ndim |
---|
437 | real(kind(1d0)) :: sig |
---|
438 | |
---|
439 | if (sig .lt. zero) then |
---|
440 | print *, "X Emittance is less then 0" |
---|
441 | return |
---|
442 | endif |
---|
443 | |
---|
444 | if ( (ndim .le. 0) .or. (ndim .gt. 6)) then |
---|
445 | print *, "Unknown dimension code", ndim |
---|
446 | return |
---|
447 | endif |
---|
448 | |
---|
449 | |
---|
450 | if (getdebug() > 1) then |
---|
451 | print *, "Setting sigma for ", ndim |
---|
452 | print *, "Current sigmas (setsigma) ", sigmas |
---|
453 | endif |
---|
454 | |
---|
455 | sigmas(ndim) = sig |
---|
456 | |
---|
457 | end subroutine setsigma |
---|
458 | !_________________________________________________________________________________ |
---|
459 | |
---|
460 | subroutine printsigmas |
---|
461 | implicit none |
---|
462 | |
---|
463 | print*,"Sigmas:", sigmas |
---|
464 | end subroutine printsigmas |
---|
465 | !_________________________________________________________________________________ |
---|
466 | |
---|
467 | subroutine initmoments() |
---|
468 | implicit none |
---|
469 | integer :: no |
---|
470 | integer :: i ! dimension |
---|
471 | integer :: get_string |
---|
472 | character(len=48), dimension(3) :: disttypes |
---|
473 | character(len=48) :: cmdname |
---|
474 | integer, dimension(3) :: stringlength |
---|
475 | character(48) :: tmpstring |
---|
476 | integer :: getcurrentcmdname |
---|
477 | ! This routine must be called before any init in ptc_twiss is performed |
---|
478 | ! since it initialize Bertz for its purpose |
---|
479 | ! |
---|
480 | ! |
---|
481 | if ( associated(normmoments) ) then |
---|
482 | deallocate(normmoments) |
---|
483 | endif |
---|
484 | |
---|
485 | if (nmoments < 1) then |
---|
486 | ! call fort_warn("initmoments","No moments specified for calculation.") |
---|
487 | return |
---|
488 | endif |
---|
489 | |
---|
490 | |
---|
491 | i = getcurrentcmdname(cmdname); |
---|
492 | |
---|
493 | if (i .eq. 0 ) then |
---|
494 | call fort_warn("initmoments","Can not get the current command name.") |
---|
495 | return |
---|
496 | endif |
---|
497 | |
---|
498 | stringlength(1) = get_string(cmdname,'xdistr ',disttypes(1)) |
---|
499 | stringlength(2) = get_string(cmdname,'ydistr ',disttypes(2)) |
---|
500 | stringlength(3) = get_string(cmdname,'zdistr ',disttypes(3)) |
---|
501 | |
---|
502 | !we take what was available in the last ptc_twiss and go to maximum order we can go |
---|
503 | no = mapsorder |
---|
504 | if ( no < 1 ) then |
---|
505 | call fort_warn('madx_ptc_distrib.f90 <initmoments>:','Order in twiss is smaller then 1') |
---|
506 | return |
---|
507 | endif |
---|
508 | no = no*2 !we take what was available in the last ptc_twiss and go to maximum order we can go |
---|
509 | |
---|
510 | allocate(normmoments(3, 0:no, 0:no)) |
---|
511 | normmoments = zero |
---|
512 | |
---|
513 | do i=1,3 |
---|
514 | tmpstring = disttypes(i) |
---|
515 | select case(tmpstring(1:stringlength(i))) |
---|
516 | case ('gauss') |
---|
517 | if (getdebug() > 1) then |
---|
518 | print*, "initmoments: Gauss distribution for dimension ", i |
---|
519 | endif |
---|
520 | call makegaus(no,i) |
---|
521 | distributiontype(i) = distr_gauss |
---|
522 | case ('flat5') |
---|
523 | if (getdebug() > 1) then |
---|
524 | print*, "initmoments: Flat distribution for dimension ", i |
---|
525 | endif |
---|
526 | call makeflat5(no,i) |
---|
527 | distributiontype(i) = distr_flat5 |
---|
528 | case ('flat56') |
---|
529 | if (getdebug() > 1) then |
---|
530 | print*, "initmoments: Flat distribution for dimension ", i |
---|
531 | endif |
---|
532 | call makeflat56(no,i) |
---|
533 | distributiontype(i) = distr_flat56 |
---|
534 | case default |
---|
535 | call fort_warn("initmoments","Distribution not recognized") |
---|
536 | print*, "initmoments: Distribution ", tmpstring(1:stringlength(i)), "not recognized" |
---|
537 | print*, "initmoments: Using default Gaussian for dimension ", i |
---|
538 | call makegaus(no,i) |
---|
539 | distributiontype(i) = distr_gauss |
---|
540 | end select |
---|
541 | |
---|
542 | enddo |
---|
543 | |
---|
544 | end subroutine initmoments |
---|
545 | !_________________________________________________________________________________ |
---|
546 | subroutine allocmoments |
---|
547 | implicit none |
---|
548 | !allocates variables needed for moments calculations |
---|
549 | ! namely gmapa and fucntion_to_avarage |
---|
550 | |
---|
551 | if ( .not. associated(normmoments) ) then |
---|
552 | return |
---|
553 | endif |
---|
554 | |
---|
555 | call alloc(gmapa,c_%nv) |
---|
556 | gmapa = sigmas |
---|
557 | call alloc(function_to_average) |
---|
558 | |
---|
559 | end subroutine allocmoments |
---|
560 | |
---|
561 | !_________________________________________________________________________________ |
---|
562 | subroutine killmoments |
---|
563 | implicit none |
---|
564 | !cleans everything allocated in in initmoments and allocmoments |
---|
565 | |
---|
566 | if ( .not. associated(normmoments) ) then |
---|
567 | return |
---|
568 | endif |
---|
569 | |
---|
570 | deallocate(normmoments) |
---|
571 | |
---|
572 | call kill(gmapa) |
---|
573 | call kill(function_to_average) |
---|
574 | |
---|
575 | end subroutine killmoments |
---|
576 | !_________________________________________________________________________________ |
---|
577 | |
---|
578 | real(dp) function filter(e) |
---|
579 | implicit none |
---|
580 | integer e(:) |
---|
581 | integer i |
---|
582 | |
---|
583 | filter=one |
---|
584 | |
---|
585 | |
---|
586 | do i=1,c_%nd |
---|
587 | |
---|
588 | filter=filter*normmoments(i,e(2*i-1),e(2*i)) |
---|
589 | if (getdebug() > 4) then |
---|
590 | print*, "normmoments(",i, e(2*i-1), e(2*i),")=", normmoments(i,e(2*i-1),e(2*i)) |
---|
591 | endif |
---|
592 | enddo |
---|
593 | |
---|
594 | if (getdebug() > 3) then |
---|
595 | |
---|
596 | print*,"filter(",e(1:6),")=",filter |
---|
597 | print*,"==================" |
---|
598 | endif |
---|
599 | |
---|
600 | end function filter |
---|
601 | |
---|
602 | |
---|
603 | !_________________________________________________________________________________ |
---|
604 | !_________________________________________________________________________________ |
---|
605 | !_________________________________________________________________________________ |
---|
606 | !! |
---|
607 | subroutine makegaus(no,d) |
---|
608 | implicit none |
---|
609 | integer no !order |
---|
610 | integer d ! dimension number |
---|
611 | integer i,j,jn(2) |
---|
612 | type(taylor) x,p,f |
---|
613 | type(Taylorresonance) fr |
---|
614 | |
---|
615 | if (getdebug() > 1) then |
---|
616 | print*, "Making Gauss distributions for dimension ", d |
---|
617 | endif |
---|
618 | |
---|
619 | call init(no,1,0,0) |
---|
620 | |
---|
621 | call alloc(x,p,f) |
---|
622 | call alloc(fr) |
---|
623 | |
---|
624 | x=one.mono.1 |
---|
625 | p=one.mono.2 |
---|
626 | |
---|
627 | f=one |
---|
628 | do i=0,no |
---|
629 | do j=i,no |
---|
630 | |
---|
631 | if(mod(i,2)/=0) cycle |
---|
632 | if(mod(j,2)/=0) cycle |
---|
633 | f=x**i*p**j |
---|
634 | |
---|
635 | fr=f |
---|
636 | |
---|
637 | jn(1)=(i+j)/2; |
---|
638 | jn(2)=jn(1); |
---|
639 | normmoments(d,i,j)=(fr%cos.sub.jn) |
---|
640 | normmoments(d,i,j)=normmoments(d,i,j)*singlefac(jn(1))*two**(jn(1)) |
---|
641 | normmoments(d,j,i)=normmoments(d,i,j) |
---|
642 | |
---|
643 | if (getdebug() > 1) then |
---|
644 | print*, "mom(",i,",",j,")=",normmoments(d,i,j) |
---|
645 | endif |
---|
646 | enddo |
---|
647 | enddo |
---|
648 | call kill(x,p,f) |
---|
649 | call kill(fr) |
---|
650 | |
---|
651 | end subroutine makegaus |
---|
652 | |
---|
653 | !_________________________________________________________________________________ |
---|
654 | |
---|
655 | subroutine makeflat5(no,d) |
---|
656 | implicit none |
---|
657 | integer no |
---|
658 | integer d !dimension number |
---|
659 | integer i,j |
---|
660 | |
---|
661 | if (getdebug() > 1) then |
---|
662 | print*, "Making flat distribution " |
---|
663 | endif |
---|
664 | |
---|
665 | do i=0,no |
---|
666 | do j=i,no |
---|
667 | |
---|
668 | if(mod(i,2)/=0) cycle |
---|
669 | |
---|
670 | normmoments(d,i,j)=(three)**(i/2)/(i+1) !delta assumed flat distribution and |
---|
671 | normmoments(d,j,i)=normmoments(d,i,j) !and L is the delta function |
---|
672 | |
---|
673 | if (getdebug() > 1) then |
---|
674 | print*, "mom(",i,",",j,")=",normmoments(d,i,j) |
---|
675 | endif |
---|
676 | |
---|
677 | enddo |
---|
678 | enddo |
---|
679 | |
---|
680 | end subroutine makeflat5 |
---|
681 | |
---|
682 | !_________________________________________________________________________________ |
---|
683 | |
---|
684 | subroutine makeflat56(no,d) |
---|
685 | implicit none |
---|
686 | integer no |
---|
687 | integer d !dimension number |
---|
688 | integer i,j |
---|
689 | |
---|
690 | if (getdebug() > 1) then |
---|
691 | print*, "Making flat in delta and T distributions" |
---|
692 | endif |
---|
693 | |
---|
694 | do i=0,no,2 |
---|
695 | do j=i,no,2 |
---|
696 | |
---|
697 | normmoments(d,i,j)=(three)**(i/2)/(i+1) !delta assumed flat distribution and |
---|
698 | normmoments(d,i,j)=normmoments(d,i,j)*(three)**(j/2)/(j+1) !delta assumed flat distribution and |
---|
699 | |
---|
700 | normmoments(d,j,i)=normmoments(d,i,j) !and L is the delta function |
---|
701 | |
---|
702 | if (getdebug() > 1) then |
---|
703 | print*, "mom(",i,",",j,")=",normmoments(d,i,j) |
---|
704 | endif |
---|
705 | |
---|
706 | enddo |
---|
707 | enddo |
---|
708 | |
---|
709 | end subroutine makeflat56 |
---|
710 | !_________________________________________________________________________________ |
---|
711 | !_________________________________________________________________________________ |
---|
712 | !_________________________________________________________________________________ |
---|
713 | |
---|
714 | |
---|
715 | !_________________________________________________________________________________ |
---|
716 | |
---|
717 | real(dp) function singlefac(n) |
---|
718 | implicit none |
---|
719 | integer n,i |
---|
720 | |
---|
721 | singlefac=one |
---|
722 | do i=1,n |
---|
723 | singlefac=singlefac*i |
---|
724 | enddo |
---|
725 | |
---|
726 | end function singlefac |
---|
727 | |
---|
728 | !_________________________________________________________________________________ |
---|
729 | |
---|
730 | integer function getdistrtype(axis) |
---|
731 | implicit none |
---|
732 | integer axis |
---|
733 | |
---|
734 | getdistrtype = distributiontype(axis) |
---|
735 | end function getdistrtype |
---|
736 | |
---|
737 | |
---|
738 | !_________________________________________________________________________________ |
---|
739 | |
---|
740 | |
---|
741 | real(dp) function getsigma(axis) |
---|
742 | implicit none |
---|
743 | integer axis |
---|
744 | |
---|
745 | if ( (axis > 0) .and. (axis < 7) ) then |
---|
746 | getsigma = sigmas(axis) |
---|
747 | else |
---|
748 | getsigma = zero |
---|
749 | endif |
---|
750 | |
---|
751 | end function getsigma |
---|
752 | |
---|
753 | end module madx_ptc_distrib_module |
---|
754 | !_________________________________________________________________________________ |
---|
755 | !_________________________________________________________________________________ |
---|
756 | !_________________________________________________________________________________ |
---|
757 | |
---|
758 | |
---|
759 | integer function w_ptc_getnmoments() |
---|
760 | use madx_ptc_distrib_module |
---|
761 | implicit none |
---|
762 | w_ptc_getnmoments = getnmoments() |
---|
763 | end function w_ptc_getnmoments |
---|
764 | |
---|
765 | !_________________________________________________________________________________ |
---|
766 | |
---|
767 | subroutine w_ptc_getmomentstabcol(n, tabn, coln) |
---|
768 | use madx_ptc_distrib_module |
---|
769 | implicit none |
---|
770 | integer :: n |
---|
771 | character(20) :: tabn |
---|
772 | character(17) :: coln |
---|
773 | |
---|
774 | call getmomentstabcol(n, tabn, coln) |
---|
775 | |
---|
776 | end subroutine w_ptc_getmomentstabcol |
---|