Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

madx_ptc_distrib.f90 @ 430

Last change on this file since 430 was 430, checked in by touze, 11 years ago
import madx-5.01.00
File size: 21.2 KB

Line
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: AYCA_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)*(2k))
362	endif
363	v = v + (function_to_average.sub.e) * (sigmas(5)*(k)) !was to 2k
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=filternormmoments(i,e(2i-1),e(2*i))
589	if (getdebug() > 4) then
590	print, "normmoments(",i, e(2i-1), e(2i),")=", normmoments(i,e(2i-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*ip**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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: PSPA/madxPSPA/src/madx_ptc_distrib.f90 @ 430

Download in other formats: