Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

madx_ptc_twiss.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: 99.1 KB

Line
1	module madx_ptc_twiss_module
2	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3	! madx_ptc_distrib module
4	! Piotr K. Skowronski , Frank Schmidt (CERN)
5	!
6	! This module contains service for twiss distributions with PTC
7	use madx_ptc_module
8	use madx_ptc_intstate_module, only : getdebug
9	USE madx_ptc_setcavs_module
10	USE madx_ptc_knobs_module
11	USE madx_ptc_distrib_module
12
13	implicit none
14
15	save
16	private
17
18	!============================================================================================
19	! PUBLIC INTERFACE
20	public :: twiss,ptc_twiss
21
22
23
24	!============================================================================================
25	! PRIVATE
26	! data structures
27
28	!PSk 2011.01.05 goes global to the modules so the slice tracking produces it for the summ table
29	type(real_8) :: theTransferMap(6)
30	type(universal_taylor) :: unimap(6)
31
32	type twiss
33
34	logical(lp) nf
35	real(dp), dimension(3,3) :: beta,alfa,gama
36	real(dp), dimension(3,3) :: beta_p,alfa_p,gama_p ! derivatives of the above w.r.t delta_p
37	real(dp), dimension(3) :: mu
38	real(dp), dimension(6) :: disp
39	real(dp), dimension(6) :: disp_p ! derivatives of the dispersion w.r.t delta_p
40	real(dp), dimension(6) :: disp_p2 ! second derivatives of the dispersion w.r.t delta_p
41	real(dp), dimension(6) :: disp_p3 ! third order derivatives of dispersion w.r.t delta_p
42	real(dp), dimension(3) :: tune
43	real(dp), dimension(6,6) :: eigen
44	end type twiss
45
46	interface assignment (=)
47	module procedure equaltwiss
48	module procedure zerotwiss
49	module procedure normalform_normalform
50	end interface assignment (=)
51
52	interface alloc
53	module procedure alloctwiss
54	end interface alloc
55
56	interface kill
57	module procedure killtwiss
58	end interface kill
59
60	type(normalform) :: Normal
61	real(dp), private, dimension(2,ndim2) :: angp
62	real(dp), private, dimension(ndim2) :: dicu
63	integer, private, parameter :: ndd=ndim2
64	integer, private, dimension(4) :: iia,icoast
65	integer, private :: np
66
67	!new lattice function
68	real(dp), private, dimension(3) :: testold
69	real(dp), private, dimension(3) :: phase
70
71	character(len=5), private, dimension(5), parameter :: str5 = (/'10000','01000','00100','00010','00001'/)
72
73	integer, private, allocatable :: J(:)
74	integer, private, dimension(6) :: j5 = (/0,0,0,0,1,0/)
75	integer, private, dimension(6) :: j6 = (/0,0,0,0,0,1/)
76	integer, private, dimension(6,6) :: fo = &
77	reshape( (/1,0,0,0,0,0,&
78	0,1,0,0,0,0,&
79	0,0,1,0,0,0,&
80	0,0,0,1,0,0,&
81	0,0,0,0,1,0,&
82	0,0,0,0,0,1 /), &
83	(/6,6/) )
84
85	logical :: slice_magnets, center_magnets, deltap_dependency, isRing
86
87	logical :: resetBetaExtrema(3,3);
88
89	real(dp) :: minBeta(3,3) ! jluc: to store extremas of Twiss functions (show-up in summary table
90	real(dp) :: maxBeta(3,3) ! jluc: to store extremas of Twiss functions (show-up in summary table)
91
92	!============================================================================================
93	! PRIVATE
94	! routines
95	private zerotwiss,equaltwiss,killtwiss
96
97
98
99	contains
100	!____________________________________________________________________________________________
101
102	subroutine equaltwiss(s1,Y)
103	implicit none
104	type(twiss), intent(inout)::s1
105	type(real_8), intent(in)::Y(6)
106	integer jj,i,k, ndel, n
107	real(dp) :: lat(0:6,6,3)
108	real(dp) :: test, dph
109	real(dp) :: epsil=1e-12 !
110	integer :: J(lnv)
111
112
113	lat = zero
114
115	n=3 ! 1 2 3 are tunes
116
117	ndel=0
118	if(c_%ndpt/=0) then
119	! print*, "We are at the mode 6D + nocav"
120	ndel=1 !this is 6D without cavity (MADX icase=56)
121	endif
122
123	J=0;
124	do i=1,c_%nd2-2*ndel
125	do jj=i,c_%nd2-2*ndel
126	do k=1,c_%nd-ndel
127	n=n+1
128	J(2*k-1)=1
129	lat(i,jj,k)= (Y(i)%t.sub.J)*(Y(jj)%t.sub.J)
130	J(2*k-1)=0
131	J(2*k)=1
132	lat(i,jj,k)=lat(i,jj,k) + (Y(i)%t.sub.J)*(Y(jj)%t.sub.J)
133	lat(jj,i,k)=lat(i,jj,k)
134	J(2*k)=0
135	enddo
136	enddo
137	enddo
138
139	J=0
140	!here ND2=4 and delta is present nd2=6 and delta is a constant
141	! print*,"nv",c_%nv,"nd2",c_%nd2,"np",c_%np,"ndpt",c_%ndpt ,"=>",c_%nv-c_%nd2-c_%np
142	if( (c_%npara==5) .or. (c_%ndpt/=0) ) then
143	!when there is no cavity it gives us dispersions
144	do i=1,4
145	lat(0,i,1)=(Y(i)%t.sub.J5)
146	enddo
147	elseif (c_%nd2 == 6) then
148	do i=1,4
149	lat(0,i,1) = (Y(i)%t.sub.J5)*(Y(6)%t.sub.J6)
150	lat(0,i,1) = lat(0,i,1) + (Y(i)%t.sub.J6)*(Y(5)%t.sub.J5)
151	enddo
152	else
153	do i=1,4
154	lat(0,i,1)=zero
155	enddo
156	endif
157
158
159	!!!!!!!!!!!!!!!!
160	! phase advance!
161	!!!!!!!!!!!!!!!!
162
163	k = 2
164	if(c_%nd2==6.and.c_%ndpt==0) k = 3
165
166	j=0
167	do i=1, k
168	jj=2*i -1
169	TEST=ATAN2((Y(2i -1).SUB.fo(2i,:)),(Y(2i-1).SUB.fo(2i-1,:)))/TWOPI
170	if (i == 3) then
171	TEST = ATAN2((Y(6).SUB.fo(5,:)),(Y(6).SUB.fo(6,:)))/TWOPI
172	endif
173
174
175	IF(TEST<zero.AND.abs(TEST)>EPSIL)TEST=TEST+one
176	DPH=TEST-TESTOLD(i)
177	IF(DPH<zero.AND.abs(DPH)>EPSIL) DPH=DPH+one
178	IF(DPH>half) DPH=DPH-one
179
180	PHASE(i)=PHASE(i)+DPH
181	TESTOLD(i)=TEST
182
183	enddo
184
185
186	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
187	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
189	do i=1,c_%nd
190	do k=1,c_%nd
191	s1%beta(k,i)= lat(2k-1,2k-1,i)
192	s1%alfa(k,i)=-lat(2k-1,2k ,i)
193	s1%gama(k,i)= lat(2k ,2k ,i)
194	enddo
195	enddo
196
197	! --- derivatives of the Twiss parameters w.r.t delta_p
198	if (deltap_dependency) then
199	if( (c_%npara==5) .or. (c_%ndpt/=0) ) then ! condition to be checked
200	call computeDeltapDependency(y,s1)
201	endif
202	endif
203	! ---
204
205
206	!when there is no cavity it gives us dispersions
207	do i=1,c_%nd2-2*ndel
208	s1%disp(i)=lat(0,i,1)
209	enddo
210
211	if (c_%nd == 3) then
212	do i=1,c_%nd
213	test = s1%beta(3,i)
214	s1%beta(3,i) = s1%gama(3,i)
215	s1%gama(3,i) = test
216	enddo
217	endif
218
219	!--- track the Twiss functions' extremas
220	call trackBetaExtrema(1,1,s1%beta(1,1))
221	call trackBetaExtrema(2,2,s1%beta(2,2))
222	!---
223
224	s1%mu=phase
225
226	do k=1,3
227	do i=1,6
228	s1%eigen(k2-1,i) = Y(k2-1).sub.fo(i,:)
229	s1%eigen(k2 ,i) = Y(k2 ).sub.fo(i,:)
230	enddo
231	enddo
232
233	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
234	call fort_warn('ptc_twiss: ','DA in twiss got unstable')
235	call seterrorflag(10,"ptc_twiss ","DA got unstable in twiss ");
236	return
237	endif
238
239	end subroutine equaltwiss
240
241
242
243	subroutine alloctwiss(s1)
244	implicit none
245	type (twiss),intent(inout)::s1
246
247	s1=0
248	end subroutine alloctwiss
249	!_________________________________________________________________
250
251	subroutine killtwiss(s1)
252	implicit none
253	type (twiss),intent(inout)::s1
254
255	s1%nf=.false.
256	s1%beta(:,:)=zero
257	s1%beta_p(:,:)=zero
258	s1%alfa(:,:)=zero
259	s1%alfa_p(:,:)=zero
260	s1%gama(:,:)=zero
261	s1%gama_p(:,:)=zero
262	s1%mu(:)=zero
263	s1%disp(:)=zero
264	s1%disp_p(:)=zero
265	s1%disp_p2(:)=zero
266	s1%disp_p3(:)=zero
267	s1%tune(:)=zero
268	s1%eigen(:,:)=zero
269	end subroutine killtwiss
270	!_________________________________________________________________
271
272	!________________________________________________________________________________
273
274	subroutine zerotwiss(s1,i)
275	implicit none
276	type(twiss), intent(inout)::s1
277	integer, intent(in)::i
278
279	if(i==0) then
280
281	call liepeek(iia,icoast)
282	NP=iia(2)-c_%nd2
283
284	s1%nf=.false.
285	s1%beta(:,:)=zero
286	s1%beta_p(:,:)=zero
287	s1%alfa(:,:)=zero
288	s1%alfa_p(:,:)=zero
289	s1%gama(:,:)=zero
290	s1%gama_p(:,:)=zero
291	s1%mu(:)=zero
292	s1%disp(:)=zero
293	s1%disp_p(:)=zero
294	s1%disp_p2(:)=zero
295	s1%disp_p3(:)=zero
296	s1%tune(:)=zero
297	s1%eigen(:,:)=zero
298	dicu(:)=zero
299	angp(:,:)=zero
300	endif
301
302	end subroutine zerotwiss
303
304
305
306
307	!_________________________________________________________________
308
309	subroutine ptc_twiss(tab_name,summary_tab_name)
310	use twissafi
311	implicit none
312	logical(lp) :: closed_orbit,beta_flg, slice, goslice
313	integer :: k,i,ii
314	integer :: no,mynd2,npara,nda,icase,flag_index,why(9),my_nv,nv_min
315	character(200) :: whymsg
316	integer :: ioptfun,iii,restart_sequ,advance_node,mf1,mf2
317	integer :: tab_name(*)
318	integer :: summary_tab_name(*)
319	real(dp) :: x(6)
320	real(dp) :: deltap0,deltap ,d_val
321	real(kind(1d0)) :: get_value,suml,s
322	integer :: posstart, posnow
323	integer :: geterrorflag !C function that returns errorflag value
324	type(real_8) :: y(6)
325	type(twiss) :: tw
326	type(fibre), POINTER :: current
327	type(integration_node), pointer :: nodePtr, stopNode
328	type(work) :: startfen !Fibre energy at the start
329	real(dp) :: r,re(ndim2,ndim2),dt
330	logical(lp) :: initial_matrix_manual, initial_matrix_table, initial_map_manual
331	logical(lp) :: initial_distrib_manual, initial_ascript_manual, writetmap
332	integer :: row, rmatrix
333	real(dp) :: emi(3)
334	logical(lp) :: isputdata
335	integer :: countSkipped
336	character(48) :: summary_table_name
337	character(12) :: tmfile='transfer.map'
338	character(48) charconv !routine
339
340	if(universe.le.0.or.EXCEPTION.ne.0) then
341	call fort_warn('return from ptc_twiss: ',' no universe created')
342	call seterrorflag(1,"ptc_twiss ","no universe created till now");
343	return
344	endif
345	if(index_mad.le.0.or.EXCEPTION.ne.0) then
346	call fort_warn('return from ptc_twiss: ',' no layout created')
347	call seterrorflag(2,"ptc_twiss ","no layout created till now");
348	return
349	endif
350
351	call resetBetaExtremas()
352
353	!skipnormalform = my_false
354	countSkipped = 0
355
356	!all zeroing
357	testold = zero
358	phase = zero
359	do i=1,6
360	unimap(i) = zero
361	enddo
362
363	if (getdebug() > 1) then
364	print*,"ptc_twiss"
365	endif
366	call momfirstinit()
367
368	!------------------------------------------------------------------------------
369	table_name = charconv(tab_name)
370	summary_table_name = charconv(summary_tab_name)
371
372	if (getdebug() > 1) then
373	print*,"ptc_twiss: Table name is ",table_name
374	print*,"ptc_twiss: Summary table name is", summary_table_name
375	endif
376
377	call cleartables() !defined in madx_ptc_knobs
378
379	call kill_para(my_ring) !removes all the previous parameters
380
381	nda = getnknobsall() !defined in madx_ptc_knobs
382	suml=zero
383
384	icase = get_value('ptc_twiss ','icase ')
385
386	deltap0 = get_value('ptc_twiss ','deltap ')
387
388	rmatrix = get_value('ptc_twiss ','rmatrix ')
389
390	deltap = zero
391
392	call my_state(icase,deltap,deltap0)
393
394	CALL UPDATE_STATES
395
396	! check that deltap-dependency selected if icase=5, which stands for
397	! 4 dimensions and deltap/p as an external parameter.
398	! indeed the simplified formulas we applied for derivation w.r.t deltap assume
399	! that deltap is an externally supplied constant parameter, which is no longer
400	! the case when icase = 6 as deltap then becomes a phase-space state-variable.
401	! and for icase=4, there is no dispersion.
402
403	deltap_dependency = get_value('ptc_twiss ','deltap_dependency ') .ne. 0
404
405	if (deltap_dependency .and. .not.((icase.eq.5) .or. (icase.eq.56))) then
406	call fort_warn('ptc_twiss: ','derivation w.r.t deltap assume deltap is fixed parameter')
407	call fort_warn('ptc_twiss: ','derivation w.r.t deltap assume icase=5 (formula differentiation) or icase=56 (from map)')
408	endif
409
410	x(:)=zero
411	if(mytime) then
412	call Convert_dp_to_dt (deltap, dt)
413	else
414	dt=deltap
415	endif
416	if(icase.eq.5 .or. icase.eq.56) x(5)=dt
417
418	closed_orbit = get_value('ptc_twiss ','closed_orbit ') .ne. 0
419
420	if( closed_orbit .and. (icav .gt. 0) .and. (my_ring%closed .eqv. .false.)) then
421	call fort_warn('return from ptc_twiss: ',' Closed orbit requested on not closed layout.')
422	call seterrorflag(3,"ptc_twiss ","Closed orbit requested on not closed layout.");
423	return
424	endif
425
426	if(closed_orbit) then
427
428	if ( .not. c_%stable_da) then
429	call fort_warn('ptc_twiss: ','DA got unstable even before finding orbit')
430	call seterrorflag(10,"ptc_twiss ","DA got unstable even before finding orbit");
431	stop
432	! return
433	endif
434
435	! pass starting point for closed orbit search
436	x(1)=get_value('ptc_twiss ','x ')
437	x(2)=get_value('ptc_twiss ','px ')
438	x(3)=get_value('ptc_twiss ','y ')
439	x(4)=get_value('ptc_twiss ','py ')
440	x(6)=get_value('ptc_twiss ','t ')
441	x(5)=x(5)+get_value('ptc_twiss ','pt ')
442
443
444
445	!print*, "Looking for orbit"
446	!w_p = 0
447	call find_orbit(my_ring,x,1,default,c_1d_7)
448
449	if ( .not. check_stable) then
450	call fort_warn('ptc_twiss: ','Can not find closed orbit')
451	call seterrorflag(10,"ptc_twiss ","Can not find closed orbit");
452	return
453	! return
454	endif
455
456	! print*, "From closed orbit", w_p%nc
457	! if ( w_p%nc .gt. 0) then
458	! do i=1,w_p%nc
459	! call fort_warn('ptc_twiss: ',w_p%c(i))
460	! call seterrorflag(10,"ptc_twiss ",w_p%c(i));
461	! enddo
462	! endif
463
464	CALL write_closed_orbit(icase,x)
465
466	elseif(my_ring%closed .eqv. .true.) then
467	print*, "Closed orbit specified by the user!"
468	!CALL write_closed_orbit(icase,x) at this position it isn't read
469	endif
470
471	mynd2 = 0
472	npara = 0
473	no = get_value('ptc_twiss ','no ')
474	if ( no .lt. 1 ) then
475	call fort_warn('madx_ptc_twiss.f90 <ptc_twiss>:','Order in twiss is smaller then 1')
476	print*, "Order is ", no
477	return
478	endif
479
480	writetmap = get_value('ptc_twiss ','writetmap ') .ne. 0
481
482	!this must be before initialization of the Bertz
483
484	initial_distrib_manual = get_value('ptc_twiss ','initial_moments_manual ') .ne. 0
485	if (initial_distrib_manual) then
486	if (getdebug() > 1) then
487	print*,"Initializing map with initial_moments_manual=true"
488	endif
489	call readinitialdistrib()
490	endif
491
492
493	call init(default,no,nda,BERZ,mynd2,npara)
494
495	!This must be before init map
496	call alloc(y)
497	y=npara
498	Y=X
499
500	! if (maxaccel .eqv. .false.) then
501	! cavsareset = .false.
502	! endif
503
504	if ( (cavsareset .eqv. .false.) .and. (my_ring%closed .eqv. .false.) ) then
505
506	call setcavities(my_ring,maxaccel)
507	if (geterrorflag() /= 0) then
508	return
509	endif
510	endif
511
512	call setknobs(my_ring)
513
514	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
515	! INIT Y that is tracked !
516	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
517	call initmap(dt)
518
519	if (geterrorflag() /= 0) then
520	!if arror occured then return
521	return
522	endif
523
524	call alloc(theTransferMap)
525	theTransferMap = npara
526	theTransferMap = X
527
528
529	!############################################################################
530	!############################################################################
531	!############################################################################
532
533	slice_magnets = get_value('ptc_twiss ','slice_magnets ') .ne. 0
534	center_magnets = get_value('ptc_twiss ','center_magnets ') .ne. 0
535
536	slice = slice_magnets .or. center_magnets
537
538	if ( slice) then
539	call make_node_layout(my_ring)
540	call getBeamBeam()
541	endif
542
543
544	!############################################################################
545	!############################################################################
546	!############################################################################
547
548	call alloc(tw)
549
550	!Y
551
552	!the initial twiss is needed to initialize propely calculation of some variables f.g. phase advance
553	tw=y
554	phase = zero !we have to do it after the very initial twiss params calculation above
555
556	current=>MY_RING%start
557	startfen = 0
558	startfen = current!setting up start energy for record
559	suml=zero
560	iii=restart_sequ()
561	print77=.false.
562	read77=.false.
563
564	if (getdebug() > 2) then
565	call kanalnummer(mf1)
566	open(unit=mf1,file='ptctwiss.txt')
567	print *, "ptc_twiss: internal state is:"
568	call print(default,6)
569	endif
570
571	call killsavedmaps() !delete all maps, if present
572	mapsorder = 0 !it is set at the end, so we are sure the twiss was successful
573
574	savemaps = get_value('ptc_twiss ','savemaps ') .ne. 0
575
576	if (savemaps) then
577	allocate(maps(MY_RING%n))
578	do i=1,MY_RING%n
579	do ii=1,6
580	call alloc(maps(i)%unimap(ii),0,0)
581	maps(i)%unimap(ii) = zero !this initializes and allocates the variables
582	enddo
583	enddo
584	else
585	nullify(maps) !assurance
586	endif
587
588
589
590
591
592	do i=1,MY_RING%n
593
594	if (getdebug() > 1) then
595	write(6,*) "##########################################"
596	write(6,'(i4, 1x,a, f10.6)') i,current%mag%name, suml
597	write(6,'(a1,a,a1)') ">",current%mag%vorname,"<"
598	write(6,'(a, f12.6, a)') "Ref Momentum ",current%mag%p%p0c," GeV/c"
599	! if (associated(current%mag%BN)) write(6,*) "k1=", current%mag%BN(2)
600	endif
601
602
603	! Can not do this trick because beam beam can be defined within those elements
604	! so even stupid markers will occur at least twice in the twiss table
605	! skowron 2012.07.03
606	!if (slice) then
607	goslice = .true.
608	! if (current%mag%kind==kind0) then ! this is a MARKER
609	! goslice = .false.
610	! elseif(current%mag%kind==kind1) then ! this is a DRIFT, they go in one step anyway
611	! goslice = .false.
612	! elseif(current%mag%kind==kind11) then ! this is a MONITOR, they go in one step anyway
613	! goslice = .false.
614	! elseif(current%mag%kind==kind12) then ! this is a HMONITOR, they go in one step anyway
615	! goslice = .false.
616	! elseif(current%mag%kind==kind13) then ! this is a VMONITOR, they go in one step anyway
617	! goslice = .false.
618	! elseif(current%mag%kind==kind14) then ! this is a INSTRUMENT, they go in one step anyway
619	! goslice = .false.
620	! endif
621	!endif
622
623	if (slice .and. goslice) then
624
625	if (getdebug() > 1) then
626	write(6,*) "##### SLICE MAGNETS"
627	endif
628
629	nodePtr => current%t1
630
631	posstart = nodePtr%pos
632	!posnow = posstart;
633
634	if (associated(current%next)) then
635	!write(6,*) " Next node exists"
636	stopNode => current%next%t1
637	else
638	stopNode => current%t2
639	endif
640
641	do while ( .not. (associated(nodePtr, stopNode)) )
642
643	s = nodePtr%next%s(1) ! s(1) is the total arc-length, s(3) the total integration-distance
644
645	!I do not know what JL meant here
646	if ((s .eq. 0d0) .and. (nodePtr%pos .eq. (my_ring%t%n+posstart-1))) then
647	s = nodePtr%s(1) + nodePtr%next%s(5) ! s of previous node + local offset
648	endif
649
650	if (getdebug() > 2) then
651	write(6,*) "##### SLICE MAGNETS NODE ",&
652	& nodePtr%pos," => ",nodePtr%pos+1," s=",s
653	endif
654
655	if (nda > 0) then
656	call track_probe_x(my_ring,y,+default, & ! +default in case of extra parameters !?
657	& node1=nodePtr%pos,node2=nodePtr%pos+1)
658	call track_probe_x(my_ring,theTransferMap,+default, & ! +default in case of extra parameters !?
659	& node1=nodePtr%pos,node2=nodePtr%pos+1)
660	else
661	call track_probe_x(my_ring,y,default, &
662	& node1=nodePtr%pos,node2=nodePtr%pos+1)
663	call track_probe_x(my_ring,theTransferMap,default, &
664	& node1=nodePtr%pos,node2=nodePtr%pos+1)
665	endif
666
667	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
668	call fort_warn('ptc_twiss: ','DA got unstable')
669	call seterrorflag(10,"ptc_twiss ","DA got unstable ");
670	if (getdebug() > 2) close(mf1)
671	return
672	endif
673
674	call PRODUCE_APERTURE_FLAG(flag_index)
675	if(flag_index/=0) then
676	call ANALYSE_APERTURE_FLAG(flag_index,why)
677	Write(6,*) "ptc_twiss unstable (Twiss parameters) element: ",i," name: ",current%MAG%name,"-programs continues "
678	write(whymsg,*) 'APERTURE error: ',why
679	call fort_warn('ptc_twiss: ',whymsg)
680	call seterrorflag(10,"ptc_twiss: ",whymsg);
681	! Write(6,*) why ! See produce aperture flag routine in sd_frame
682	goto 100
683	endif
684
685	isputdata = .false.
686
687	!!!!!!!!!!!!!!!
688
689	if (center_magnets ) then
690	if ( associated(nodePtr,current%tm) ) then
691	if (mod(current%mag%p%nst,2)/=0) then !checking if number od slices is even
692	!here it is odd (we should be in the middle)
693	if (current%mag%L==0) then
694	! this is a zero-length marker or monitor or equivalent keep it
695	isputdata = .true.
696	elseif (current%mag%kind==31) then ! this is a DRIFT
697	! write(0,*) 'SKIP a drift of length',fibrePtr%mag%L,'was cut into an odd number of slices: ',fibrePtr%mag%p%nst
698	countSkipped = countSkipped + 1
699	else ! neither a zero-length element, nor a drift
700	! write(0,*) 'SKIP element of name',fibrePtr%mag%name,'and kind',fibrePtr%mag%kind
701	countSkipped = countSkipped + 1
702	endif
703	else ! this element has an even number of slices
704	isputdata = .true.
705	endif
706	endif
707	else
708	if (nodePtr%next%cas==case0) then
709	!not center_magnets, take every reasonable node
710	! this an inner integration node i.e. neither an extremity nor a fringe node, both to be discarded
711	isputdata = .true.
712	else
713	if (getdebug() > 2) then
714	write(6,*) " Not Saving data CASE=",nodePtr%next%cas
715	endif
716	endif
717	endif
718
719	if (isputdata ) then
720
721	if (getdebug() > 2) then
722	write(6,*) " Saving data CASE=",nodePtr%next%cas
723	endif
724
725	tw = y ! set the twiss parameters, with y being equal to the A_ phase advance
726	suml = s;
727
728	call puttwisstable(theTransferMap)
729	call putusertable(i,current%mag%name,suml,getdeltae(),theTransferMap, y)
730
731	!else
732	! write(6,*) " NOT Saving data"
733	endif
734
735	nodePtr => nodePtr%next
736	enddo
737
738	if (isputdata .eqv. .false.) then ! always save the last point if it was not yet saved
739	!write(6,*) " END OF ELEMENT Saving data"
740
741	if (getdebug() > 2) then
742	write(6,*) " Saving any way, it is the last node"
743	endif
744
745	tw = y ! set the twiss parameters, with y being equal to the A_ phase advance
746	if (s > suml) then !work around against last element having s=0
747	suml = s;
748	endif
749
750	call puttwisstable(theTransferMap)
751	call putusertable(i,current%mag%name,suml,getdeltae(),theTransferMap, y)
752
753	endif
754
755	else
756
757	if (nda > 0) then
758	! if (getnknobis() > 0) c_%knob = my_true
759	!print*, "parametric",i,c_%knob
760	call track(my_ring,y,i,i+1,+default)
761	call track(my_ring,theTransferMap,i,i+1,+default)
762	else
763	call track(my_ring,y,i,i+1, default)
764	call track(my_ring,theTransferMap,i,i+1,default)
765	endif
766
767
768	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
769	call fort_warn('ptc_twiss: ','DA got unstable')
770	call seterrorflag(10,"ptc_twiss ","DA got unstable ");
771	if (getdebug() > 2) close(mf1)
772	return
773	endif
774
775	call PRODUCE_APERTURE_FLAG(flag_index)
776	if(flag_index/=0) then
777	call ANALYSE_APERTURE_FLAG(flag_index,why)
778	Write(6,*) "ptc_twiss unstable (Twiss parameters) element: ",i," name: ",current%MAG%name,"-programs continues "
779	write(whymsg,*) 'APERTURE error: ',why
780	call fort_warn('ptc_twiss: ',whymsg)
781	call seterrorflag(10,"ptc_twiss: ",whymsg);
782	! Write(6,*) why ! See produce aperture flag routine in sd_frame
783	goto 100
784	endif
785
786	if (getdebug() > 2) then
787	write(mf1,*) "##########################################"
788	write(mf1,'(i4, 1x,a, f10.6)') i,current%mag%name, suml
789	call print(y,mf1)
790	endif
791
792	suml=suml+current%MAG%P%ld
793
794	if (savemaps) then
795	do ii=1,6
796	maps(i)%unimap(ii) = y(ii)
797	enddo
798	maps(i)%s = suml
799	maps(i)%name = current%mag%name
800	endif
801
802	! compute the Twiss parameters
803	tw=y
804
805	call puttwisstable(theTransferMap)
806
807	call putusertable(i,current%mag%name,suml,getdeltae(),theTransferMap,y)
808
809
810	endif
811
812	iii=advance_node()
813	current=>current%next
814	enddo
815
816	100 continue
817
818
819	call orbitRms(summary_table_name) ! this function fills the summary table and header with the closed orbits RMS and extrema
820
821	! relocated the following here to avoid side-effect
822	print77=.false.
823	read77=.false.
824
825
826	if (writetmap) then
827	!'=== TRANSFER MAP ==='
828	call kanalnummer(mf2)
829	open(unit=mf2,file=tmfile)
830	call print(theTransferMap,mf2)
831	close(mf2)
832	endif
833
834	if (getdebug() > 2) then
835
836
837	call kanalnummer(mf2)
838	! avoid file name conflict between slice.madx and center.madx
839	! which are located under same testing directory
840	if (get_value('ptc_twiss ','center_magnets ').ne.0) then
841	open(unit=mf2,file='end_center.map')
842	else
843	open(unit=mf2,file='end.map')
844	endif
845
846	call print(y,mf2)
847
848	close(mf2)
849
850	endif
851
852
853	! 26 november 2009
854	if(isRing .eqv. .true.) then
855	call oneTurnSummary(isRing, theTransferMap , x, suml)
856	else
857	print*, "Reduced SUMM Table (closed orbit not requested)"
858	call onePassSummary(theTransferMap , x, suml)
859	endif
860
861
862
863	call set_option('ptc_twiss_summary ',1)
864	! 26 november 2009: comment the following and replace by the above
865	! if ( (momentumCompactionToggle .eqv. .true.) .and. (getenforce6D() .eqv. .false.)) then
866	! ! only makes sense if the lattice is a ring (skipped for a line lattice)
867	! call oneTurnSummary()
868	! call set_option('ptc_twiss_summary ', 1)
869	! else
870	! call set_option('ptc_twiss_summary ',0) ! for time-being, do not support lines
871	! endif
872
873	if (getdebug() > 1) then
874	write(6,*) "##########################################"
875	write(6,*) "##########################################"
876	write(6,*) "### END OF PTC_TWISS #######"
877	write(6,*) "##########################################"
878	write(6,*) "##########################################"
879	! if (associated(current%mag%BN)) write(6,*) "k1=", current%mag%BN(2)
880	endif
881
882	c_%watch_user=.false.
883
884	call kill(tw)
885	CALL kill(y)
886
887	CALL kill(theTransferMap)
888
889	do i=1,6
890	call kill(unimap(i))
891	enddo
892
893	call finishknobs()
894
895	if (savemaps) then !do it at the end, so we are sure the twiss was successful
896	mapsorder = no
897	mapsicase = icase
898	if (getnmoments() > 0) call ptc_moments(no*2) !calcualate moments with the maximum available order
899	endif
900
901
902	! f90flush is not portable, and useless...
903	! call f90flush(20,my_false)
904
905	if (getdebug() > 2) close(mf1)
906
907	!****************************************************************************************
908	!******* E N D O F PTC_TWISS **********************************************
909	!****************************************************************************************
910	!________________________________________________________________________________________
911
912	contains ! what follows are internal subroutines of ptc_twiss
913	!____________________________________________________________________________________________
914
915	subroutine initmap(dt)
916	implicit none
917	integer :: double_from_table_row
918	integer :: mman, mtab, mascr, mdistr !these variable allow to check if the user did not put too many options
919	integer :: mmap
920	real(dp) :: dt
921
922	beta_flg = (get_value('ptc_twiss ','betx ').gt.0) .and. (get_value('ptc_twiss ','bety ').gt.0)
923
924	mman = get_value('ptc_twiss ','initial_matrix_manual ')
925	mtab = get_value('ptc_twiss ','initial_matrix_table ')
926	mascr = get_value('ptc_twiss ','initial_ascript_manual ')
927	mdistr = get_value('ptc_twiss ','initial_moments_manual ')
928	mmap = get_value('ptc_twiss ','initial_map_manual ')
929
930
931	isRing = .false. ! set to true in the case the map
932
933
934
935	! is calculated over a ring, about the closed orbit. Later-on in the same subroutine.
936
937
938	initial_matrix_manual = mman .ne. 0
939	initial_map_manual = mmap .ne. 0
940	initial_matrix_table = mtab .ne. 0
941	initial_ascript_manual = mascr .ne. 0
942
943
944	if ( (mman + mtab + mascr + mdistr + mmap) > 1) then
945	call seterrorflag(11,"ptc_twiss ","Ambigous command options");
946	print*, "Only one of the following switches might be on:"
947	print*, "initial_matrix_manual = ", initial_matrix_manual
948	print*, "initial_map_manual = ", initial_map_manual
949	print*, "initial_matrix_table = ", initial_matrix_table
950	print*, "initial_ascript_manual = ", initial_ascript_manual
951	print*, "initial_moments_manual = ", mdistr
952	endif
953
954	! print*, "initial_distrib_manual is ",initial_distrib_manual
955
956	if(initial_matrix_table) then
957	k = double_from_table_row("map_table ", "nv ", 1, doublenum)
958	if(k.ne.-1) then
959	call liepeek(iia,icoast)
960	my_nv=int(doublenum)
961	nv_min=min(c_%npara,my_nv)
962	else
963	initial_matrix_table=.false.
964	endif
965	endif
966
967	if(initial_matrix_table) then
968
969	if (getdebug() > 1) then
970	print*,"Initializing map with initial_matrix_table=true"
971	endif
972	call readmatrixfromtable()
973
974	elseif(initial_ascript_manual) then
975
976	if (getdebug() > 1) then
977	print*,"Initializing map with initial_ascript_manual=true"
978	endif
979	call readinitialascript()
980	if (geterrorflag() /= 0) then
981	return
982	endif
983
984	elseif(initial_map_manual) then
985	if (getdebug() > 1) then
986	print*,"Initializing map with initial_map_manual=true"
987	endif
988	call readinitialmap()
989	if (geterrorflag() /= 0) then
990	return
991	endif
992
993	elseif(initial_matrix_manual) then
994
995	if (getdebug() > 1) then
996	print*,"Initializing map with initial_matrix_manual=true"
997	endif
998	call readinitialmatrix()
999
1000	if (geterrorflag() /= 0) then
1001	return
1002	endif
1003	elseif (initial_distrib_manual) then
1004	!matrix is already prepared beforehand
1005	if (getdebug() > 1) then
1006	print*,"Initializing map with initial_moments_manual=true"
1007	print*, "Initializing map from prepared UniTaylor"
1008	endif
1009	call readreforbit() !reads x
1010
1011	do i=1, c_%nd2
1012	y(i) = unimap(i)
1013	enddo
1014
1015	if (geterrorflag() /= 0) then
1016	return
1017	endif
1018	elseif(beta_flg) then
1019
1020	if (getdebug() > 1) then
1021	print*,"Initializing map with initial twiss parameters"
1022	endif
1023
1024	call readinitialtwiss(dt)
1025
1026	if (geterrorflag() /= 0) then
1027	return
1028	endif
1029	else
1030
1031	if (getdebug() > 1) then
1032	print*,"Initializing map from one turn map: Start Map"
1033	call print(y,6)
1034	endif
1035
1036	isRing = .true. ! compute momemtum compaction factor, tunes, chromaticies for ring
1037
1038	call track(my_ring,y,1,default)
1039	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
1040	call fort_warn('ptc_twiss: ','DA got unstable (one turn map production)')
1041	call seterrorflag(10,"ptc_twiss ","DA got unstable (one turn map production)");
1042	return
1043	endif
1044
1045	call PRODUCE_APERTURE_FLAG(flag_index)
1046	if(flag_index/=0) then
1047	call ANALYSE_APERTURE_FLAG(flag_index,why)
1048
1049	write(whymsg,*) 'APERTURE unstable (one turn map production) - programs continues: ',why
1050	call fort_warn('ptc_twiss: ',whymsg)
1051	call seterrorflag(10,"ptc_twiss: ",whymsg);
1052	! Write(6,*) "ptc_twiss unstable (map production)-programs continues "
1053	! Write(6,*) why ! See produce aperture flag routine in sd_frame
1054	c_%watch_user=.false.
1055	CALL kill(y)
1056	return
1057	endif
1058
1059	if (getdebug() > 1) then
1060	print*,"Initializing map from one turn map. One Turn Map"
1061	call print(y,6)
1062	endif
1063
1064	call maptoascript()
1065
1066	call reademittance()
1067
1068
1069	endif
1070	end subroutine initmap
1071	!____________________________________________________________________________________________
1072
1073	function getdeltae()
1074	implicit none
1075	real(dp) :: getdeltae
1076	type(work) :: cfen !current fibre energy
1077
1078	cfen = 0 ! do not remove -> if it is removed energy is wrong because = adds energy to the previous value
1079
1080	if ( (associated(current%next) .eqv. .false. ) .or. (associated( current%next, my_ring%start)) ) then
1081	! if current is the last element in the sequence i.e.
1082	! p%next == NULL (LINE) OR
1083	! p%next points the first element (CIRCLE)
1084	cfen=current
1085
1086	if (getdebug()>1) then
1087	!if it is the last element in the line
1088	print *, 'It is the last element ', current%mag%name
1089	!(it is always marker, i.e element that does not change reference energy)
1090	print *, 'Its reference energy is ', cfen%p0c
1091	endif
1092
1093	!take its reference energy
1094	else
1095	cfen=current%next ! energy after passing this element
1096	endif
1097
1098	getdeltae = cfen%energy/startfen%energy
1099
1100	if (getdebug() > 2) then
1101	write(mf1,'(3(a, f10.6))') "Ref Momentum ",cfen%p0c," Energy ", cfen%energy," DeltaE ",getdeltae
1102	endif
1103
1104	end function getdeltae
1105	!____________________________________________________________________________________________
1106
1107	subroutine puttwisstable(transfermap)
1108	implicit none
1109	include "madx_ptc_knobs.inc"
1110	integer i1,i2,ii,i1a,i2a
1111	real(kind(1d0)) :: opt_fun(150),myx ! opt_fun(72) -> opt_fun(81)
1112	! increase to 150 to have extra space beyond what's needed to accomodate additional derivatives w.r.t. delta_p
1113	real(kind(1d0)) :: deltae
1114	type(real_8), target :: transfermap(6)
1115	! added on 3 November 2010 to hold Edwards & Teng parametrization
1116	real(dp) :: betx,bety,alfx,alfy,R11,R12,R21,R22
1117	! to convert between Ripken and Edwards-Teng parametrization
1118	real(dp) :: kappa,u,ax,ay,kx,ky,kxy2,usqrt,bx,by,cx,cy,cosvp,sinvp,cosvm,sinvm,cosv2,sinv2,cosv1,sinv1
1119	real(dp) :: deltaeValue
1120
1121	if (getdebug() > 2) then
1122	write(mf1,*) "##########################################"
1123	write(mf1,*) ""
1124	write(mf1,'(i4, 1x,a, f10.6)') i,current%mag%name,suml
1125	write(mf1,*) ""
1126	call print(y,mf1)
1127	endif
1128
1129
1130	deltae = getdeltae()
1131
1132	call double_to_table_curr(table_name, 's ', suml)
1133
1134	doublenum = deltae * startfen%energy
1135	call double_to_table_curr(table_name, 'energy ', doublenum)
1136
1137
1138	opt_fun(:)=zero
1139
1140	call liepeek(iia,icoast)
1141	allocate(j(c_%npara))
1142	j(:)=0
1143	do ii=1,c_%npara ! fish
1144	opt_fun(ii)=y(ii)%T.sub.j
1145	enddo
1146
1147	myx=opt_fun(6)
1148	opt_fun(6)=opt_fun(5)
1149	opt_fun(5)=myx
1150	deallocate(j)
1151
1152	ioptfun=6
1153	call vector_to_table_curr(table_name, 'x ', opt_fun(1), ioptfun)
1154
1155
1156	opt_fun(1) = transfermap(1).sub.fo(1,:)
1157	opt_fun(2) = transfermap(1).sub.fo(2,:)
1158	opt_fun(3) = transfermap(1).sub.fo(3,:)
1159	opt_fun(4) = transfermap(1).sub.fo(4,:)
1160	opt_fun(5) = transfermap(1).sub.fo(6,:)
1161	opt_fun(6) = transfermap(1).sub.fo(5,:)
1162
1163
1164	opt_fun(7) = transfermap(2).sub.fo(1,:)
1165	opt_fun(8) = transfermap(2).sub.fo(2,:)
1166	opt_fun(9) = transfermap(2).sub.fo(3,:)
1167	opt_fun(10)= transfermap(2).sub.fo(4,:)
1168	opt_fun(11)= transfermap(2).sub.fo(6,:)
1169	opt_fun(12)= transfermap(2).sub.fo(5,:)
1170
1171	opt_fun(13)= transfermap(3).sub.fo(1,:)
1172	opt_fun(14)= transfermap(3).sub.fo(2,:)
1173	opt_fun(15)= transfermap(3).sub.fo(3,:)
1174	opt_fun(16)= transfermap(3).sub.fo(4,:)
1175	opt_fun(17)= transfermap(3).sub.fo(6,:)
1176	opt_fun(18)= transfermap(3).sub.fo(5,:)
1177
1178	opt_fun(19)= transfermap(4).sub.fo(1,:)
1179	opt_fun(20)= transfermap(4).sub.fo(2,:)
1180	opt_fun(21)= transfermap(4).sub.fo(3,:)
1181	opt_fun(22)= transfermap(4).sub.fo(4,:)
1182	opt_fun(23)= transfermap(4).sub.fo(6,:)
1183	opt_fun(24)= transfermap(4).sub.fo(5,:)
1184
1185
1186	opt_fun(25)= transfermap(6).sub.fo(1,:)
1187	opt_fun(26)= transfermap(6).sub.fo(2,:)
1188	opt_fun(27)= transfermap(6).sub.fo(3,:)
1189	opt_fun(28)= transfermap(6).sub.fo(4,:)
1190	opt_fun(29)= transfermap(6).sub.fo(6,:)
1191	opt_fun(30)= transfermap(6).sub.fo(5,:)
1192
1193
1194	opt_fun(31)= transfermap(5).sub.fo(1,:)
1195	opt_fun(32)= transfermap(5).sub.fo(2,:)
1196	opt_fun(33)= transfermap(5).sub.fo(3,:)
1197	opt_fun(34)= transfermap(5).sub.fo(4,:)
1198	opt_fun(35)= transfermap(5).sub.fo(6,:)
1199	opt_fun(36)= transfermap(5).sub.fo(5,:)
1200
1201	ioptfun=36
1202	call vector_to_table_curr(table_name, 're11 ', opt_fun(1), ioptfun)
1203
1204
1205
1206	opt_fun(beta11)= tw%beta(1,1) * deltae ! beta11=1
1207	opt_fun(beta12)= tw%beta(1,2) * deltae
1208	opt_fun(beta13)= tw%beta(1,3) * deltae
1209	opt_fun(beta21)= tw%beta(2,1) * deltae
1210	opt_fun(beta22)= tw%beta(2,2) * deltae
1211	opt_fun(beta23)= tw%beta(2,3) * deltae
1212	opt_fun(beta31)= tw%beta(3,1) * deltae
1213	opt_fun(beta32)= tw%beta(3,2) * deltae
1214	opt_fun(beta33)= tw%beta(3,3) * deltae
1215
1216	opt_fun(alfa11)= tw%alfa(1,1) * deltae
1217	opt_fun(alfa12)= tw%alfa(1,2) * deltae
1218	opt_fun(alfa13)= tw%alfa(1,3) * deltae
1219	opt_fun(alfa21)= tw%alfa(2,1) * deltae
1220	opt_fun(alfa22)= tw%alfa(2,2) * deltae
1221	opt_fun(alfa23)= tw%alfa(2,3) * deltae
1222	opt_fun(alfa31)= tw%alfa(3,1) * deltae
1223	opt_fun(alfa32)= tw%alfa(3,2) * deltae
1224	opt_fun(alfa33)= tw%alfa(3,3) * deltae
1225
1226	opt_fun(gama11)= tw%gama(1,1) * deltae
1227	opt_fun(gama12)= tw%gama(1,2) * deltae
1228	opt_fun(gama13)= tw%gama(1,3) * deltae
1229	opt_fun(gama21)= tw%gama(2,1) * deltae
1230	opt_fun(gama22)= tw%gama(2,2) * deltae
1231	opt_fun(gama23)= tw%gama(2,3) * deltae
1232	opt_fun(gama31)= tw%gama(3,1) * deltae
1233	opt_fun(gama32)= tw%gama(3,2) * deltae
1234	opt_fun(gama33)= tw%gama(3,3) * deltae
1235
1236
1237	! --- derivatives of Twiss paramters w.r.t delta_p
1238	! NOW why do we need to multiply by deltae, as for the other Twiss parameters?
1239	if (deltap_dependency) then
1240	opt_fun(beta11p)= tw%beta_p(1,1) * deltae
1241	opt_fun(beta12p)= tw%beta_p(1,2) * deltae
1242	opt_fun(beta13p)= tw%beta_p(1,3) * deltae
1243	opt_fun(beta22p)= tw%beta_p(2,1) * deltae
1244	opt_fun(beta22p)= tw%beta_p(2,2) * deltae
1245	opt_fun(beta23p)= tw%beta_p(2,3) * deltae
1246	opt_fun(beta32p)= tw%beta_p(3,2) * deltae
1247	opt_fun(beta33p)= tw%beta_p(3,3) * deltae
1248
1249	opt_fun(alfa11p)= tw%alfa_p(1,1) * deltae
1250	opt_fun(alfa12p)= tw%alfa_p(1,2) * deltae
1251	opt_fun(alfa13p)= tw%alfa_p(1,3) * deltae
1252	opt_fun(alfa21p)= tw%alfa_p(2,1) * deltae
1253	opt_fun(alfa22p)= tw%alfa_p(2,2) * deltae
1254	opt_fun(alfa23p)= tw%alfa_p(2,3) * deltae
1255	opt_fun(alfa31p)= tw%alfa_p(3,1) * deltae
1256	opt_fun(alfa32p)= tw%alfa_p(3,2) * deltae
1257	opt_fun(alfa33p)= tw%alfa_p(3,3) * deltae
1258
1259	opt_fun(gama11p)= tw%gama_p(1,1) * deltae
1260	opt_fun(gama12p)= tw%gama_p(1,2) * deltae
1261	opt_fun(gama13p)= tw%gama_p(1,3) * deltae
1262	opt_fun(gama21p)= tw%gama_p(2,1) * deltae
1263	opt_fun(gama22p)= tw%gama_p(2,2) * deltae
1264	opt_fun(gama23p)= tw%gama_p(2,3) * deltae
1265	opt_fun(gama31p)= tw%gama_p(3,1) * deltae
1266	opt_fun(gama32p)= tw%gama_p(3,2) * deltae
1267	opt_fun(gama33p)= tw%gama_p(3,3) * deltae
1268	endif
1269	! --- end
1270
1271	! march 10th: do we need to multiply by deltae the following?
1272	opt_fun(mu1)=tw%mu(1) !* deltae
1273	opt_fun(mu2)=tw%mu(2) !* deltae
1274	opt_fun(mu3)=tw%mu(3) !* deltae
1275
1276	! write(0,*),"DEBUG = ", tw%mu(3) ! should give Qs
1277
1278	opt_fun(disp1)=tw%disp(1) ! was 31 instead of 57
1279	opt_fun(disp2)=tw%disp(2) ! was 32 instead of 58
1280	opt_fun(disp3)=tw%disp(3) ! was 33 instead of 59
1281	opt_fun(disp4)=tw%disp(4) ! was 34 instead of 60
1282
1283	! 9 march 2009: add 4 items
1284
1285	if (deltap_dependency) then
1286	opt_fun(disp1p) = tw%disp_p(1)
1287	opt_fun(disp2p) = tw%disp_p(2)
1288	opt_fun(disp3p) = tw%disp_p(3)
1289	opt_fun(disp4p) = tw%disp_p(4)
1290	! 20 july 2009: add 8 items for second derivatives w.r.t deltap
1291	opt_fun(disp1p2) = tw%disp_p2(1)
1292	opt_fun(disp2p2) = tw%disp_p2(2)
1293	opt_fun(disp3p2) = tw%disp_p2(3)
1294	opt_fun(disp4p2) = tw%disp_p2(4)
1295	opt_fun(disp1p3) = tw%disp_p3(1)
1296	opt_fun(disp2p3) = tw%disp_p3(2)
1297	opt_fun(disp3p3) = tw%disp_p3(3)
1298	opt_fun(disp4p3) = tw%disp_p3(4)
1299	endif
1300
1301	! JLUC TODO
1302	! opt_fun(61)=zero disp4 is now 61 in madx_ptc_knobs.inc
1303	! jluc: left the following umodified, except 36->62
1304	opt_fun(62+4+8+6)=zero ! was 36 instead of 62 => on 9 march add 4 => on 3 July 2009 add 4 => on 3 November 2010 add 6
1305	do i1=1,6
1306	if(i1.le.4) then
1307	i1a=i1
1308	elseif(i1.eq.5) then
1309	i1a=6
1310	else
1311	i1a=5
1312	endif
1313	do i2=1,6
1314	if(i2.le.4) then
1315	i2a=i2
1316	elseif(i2.eq.5) then
1317	i2a=6
1318	else
1319	i2a=5
1320	endif
1321	! idiotic counting reset (62+4+8+6) ==> 73
1322	ii=73+(i1a-1)*6+i2a ! was 36 instead of 62 => now (62+4) instead of 62 => 62+4+4
1323	opt_fun(ii)=tw%eigen(i1,i2) * deltae
1324	! where do these eigen values go? no such dedicated column defined in madx_ptc_knobs.inc
1325	if(mytime.and.i2a.eq.6) opt_fun(ii)=-opt_fun(ii)
1326	enddo
1327	enddo
1328
1329	if (getdebug() > 2) then
1330	! this part of the code is out of sync with the rest
1331	write(6,'(a,1(f9.4,1x))') current%MAG%name,suml
1332	write(6,'(a,1(f10.7,1x))') "Delta E ", deltae
1333	write(6,'(a,3(i8.0,1x))') "idxes ", beta11,beta22,beta33
1334	write(6,'(a,3(f9.4,1x))') "betas raw ", tw%beta(1,1),tw%beta(2,2),tw%beta(3,3)
1335	write(6,'(a,3(f9.4,1x))') "betas w/ener ", opt_fun(1),opt_fun(5),opt_fun(9)
1336	write(6,'(a,4(f9.4,1x))') "dispersions ", opt_fun(57),opt_fun(58),opt_fun(59),opt_fun(60)
1337	write(6,'(a,3(f9.4,1x))') "phase adv. ", tw%mu(1),tw%mu(2),tw%mu(3)
1338	endif
1339
1340	! the following works : we see the list of all elements in sequence - what about twiss_ptc_line & twiss_ptc_ring?
1341	! jluc debug - begin
1342	!write(28,'(a,1(f8.4,1x))') current%MAG%name,suml
1343	! jluc debug - end
1344
1345	! on July 3rd 2009, add another 8 for second/third derivatives of dispersions w.r.t. deltap
1346	ioptfun=81+4+8+6 !72->81 to accomodate additional derivatives w.r.t. delta_p => should one add 4 to this one, as above?
1347	! actually 3*21+6 (???) elements from beta11 to include up to disp6p
1348	! on november 3rd 2010, added 6
1349
1350	! overwrote the above for which I am not sure where the value comes from
1351	ioptfun = 79 + 36 ! 79 as for ntwisses in madx_ptc_knobs.inc + 36 eigenvalues
1352
1353	! LD: update columns access from twiss column description in mad_gcst.c
1354	! WAS: fill contiguous data in one-go, from beta11 up to mu1, mu2, mu3
1355	! NOW: fill contiguous data in one-go, from beta11 up to the end
1356	! eigenvalue seems to be retrieved ~50 lines above...
1357	ioptfun = 18*3+16+3+36+1 ! = 110
1358	call vector_to_table_curr(table_name, 'beta11 ', opt_fun(1), ioptfun)
1359
1360	! convert between the Ripken and Edwards-Teng parametrization
1361	! according to the formulas in "BETATRON MOTION WITH COUPLING OF HORIZONTAL AND VERTICAL DEGREES OF FREEDOM"
1362	! from V. A. Lebedev and S. A. Bogacz
1363
1364	deltaeValue = deltae ! equals 1.0 unless there is a cavity
1365
1366	r11 = zero
1367	r12 = zero
1368	r21 = zero
1369	r22 = zero
1370
1371	if (tw%beta(1,2)==zero .and. tw%beta(2,1)==zero) then
1372
1373	! in case there is absolutely no coupling kx and ky will be zero and u will be NaN
1374	! and betx, bety, alfx, alfy will also evaluate as NaN if we apply the above formulae
1375	! therefore we simply copy beta11 into betx and beta22 into bety in this case, so as
1376	! to get the same values between twiss and ptc_twiss
1377	! beta11, alfa11 etc... are multiplied by deltae before output
1378	! hence we reflect this in the formula from Lebedev
1379	betx = tw%beta(1,1) * deltaeValue
1380	bety = tw%beta(2,2) * deltaeValue
1381	alfx = tw%alfa(1,1) * deltaeValue
1382	alfy = tw%alfa(2,2) * deltaeValue
1383
1384	else
1385
1386	kx=sqrt(tw%beta(1,2)/tw%beta(1,1)); ! multiplication by deltae in numerator and denominator
1387	ky=sqrt(tw%beta(2,1)/tw%beta(2,2));
1388
1389	! beta11, alfa11 etc... are multiplied by deltae before output
1390	ax=kxtw%alfa(1,1) deltaeValue -tw%alfa(1,2) * deltaeValue /kx;
1391	! hence we reflect this in the formula from Lebedev
1392	ay=kytw%alfa(2,2) deltaeValue -tw%alfa(2,1) * deltaeValue /ky;
1393	kxy2=kxkxky*ky;
1394	if((abs(kxkx-kyky).gt.TINY(ONE)).and.(abs(one-kxy2).gt.TINY(ONE))) then
1395	usqrt=kxy2(one+(axax-ayay)/(kxkx-kyky)(one-kxy2))
1396	if(usqrt.gt.TINY(ONE)) then
1397	usqrt=sqrt(usqrt)
1398	! u=(-kxy2+usqrt)/(one-kxy2)
1399	if(kxy2.le.usqrt) THEN
1400	u=(-kxy2+usqrt)/(one-kxy2)
1401	else
1402	u=(-kxy2-usqrt)/(one-kxy2)
1403	endif
1404	else
1405	u=-kxy2/(one-kxy2)
1406	endif
1407
1408	! betx, bety, alfx, alfy are the values computed by twiss with very good precision
1409	! beta11, alfa11 etc... are multiplied by deltae before output
1410	! hence we reflect this in the formula from Lebedev
1411
1412	kappa=one-u
1413
1414	betx = (tw%beta(1,1)/kappa) * deltaeValue
1415	bety = (tw%beta(2,2)/kappa) * deltaeValue
1416	alfx = (tw%alfa(1,1)/kappa) * deltaeValue
1417	alfy = (tw%alfa(2,2)/kappa) * deltaeValue
1418
1419	bx = kx*kappa+u/kx
1420	by = ky*kappa-u/ky
1421	cx = kx*kappa-u/kx
1422	cy = ky*kappa+u/ky
1423
1424	cosvp = (axay-bxcy)/(ayay+cycy)
1425	sinvp = (axcy+aybx)/(ayay+cycy)
1426	cosvm = (axay+bycx)/(axax+cxcx)
1427	sinvm = (axby-aycx)/(axax+cxcx)
1428
1429	cosv2 = sqrt((one+cosvpcosvm-sinvpsinvm)/two)
1430	sinv2 = -sqrt((one-cosvpcosvm+sinvpsinvm)/two)
1431	cosv1 = -sqrt((one+cosvpcosvm+sinvpsinvm)/two)
1432	sinv1 = -sqrt((one-cosvpcosvm-sinvpsinvm)/two)
1433
1434	r11 = sqrt(tw%beta(2,2)/tw%beta(1,2))(tw%alfa(1,2)sinv2+u*cosv2)/kappa
1435	r12 = sqrt(tw%beta(1,1)tw%beta(2,1))sinv1/kappa
1436	r21 = (cosv2(tw%alfa(1,2)kappa-tw%alfa(2,2)u)-sinv2&
1437	(ukappa+tw%alfa(1,2)tw%alfa(2,2)))/(kappasqrt(tw%beta(1,2)tw%beta(2,2)))
1438	r22 = -sqrt(tw%beta(1,1)/tw%beta(2,1))(ucosv1+tw%alfa(2,1)*sinv1)/kappa
1439	else
1440	call fort_warn("ptc_twiss","Argument of sqrt(kxy2(1+(axax-ayay)/(kxkx-kyky)(one-kxy2))) smaller than TINY")
1441	print,"Argument of sqrt(kxy2(1+(axax-ayay)/(kxkx-kyky)*(one-kxy2))) is: ",&
1442	kxy2(1+(axax-ayay)/(kxkx-kyky)(one-kxy2))
1443	betx = tw%beta(1,1) * deltaeValue
1444	bety = tw%beta(2,2) * deltaeValue
1445	alfx = tw%alfa(1,1) * deltaeValue
1446	alfy = tw%alfa(2,2) * deltaeValue
1447	endif
1448
1449	endif
1450
1451	! Edwards-Teng parameters go into betx, bety, alfx, alfy which are at the beginning of twiss_table_cols in madxl.h
1452	call double_to_table_curr(table_name, 'betx ', betx ) ! non contiguous with the above table entries
1453	call double_to_table_curr(table_name, 'bety ', bety ) ! hence we must store these values one by one
1454	call double_to_table_curr(table_name, 'alfx ', alfx )
1455	call double_to_table_curr(table_name, 'alfy ', alfy )
1456	call double_to_table_curr(table_name, 'r11 ', r11 )
1457	call double_to_table_curr(table_name, 'r12 ', r12 )
1458	call double_to_table_curr(table_name, 'r21 ', r21 )
1459	call double_to_table_curr(table_name, 'r22 ', r22 )
1460
1461	call augment_count(table_name)
1462
1463
1464	end subroutine puttwisstable
1465	!____________________________________________________________________________________________
1466
1467	subroutine readrematrix
1468	!reads covariance matrix of the initial distribution
1469	implicit none
1470
1471	re(1,1) = get_value('ptc_twiss ','re11 ')
1472	re(1,2) = get_value('ptc_twiss ','re12 ')
1473	re(1,3) = get_value('ptc_twiss ','re13 ')
1474	re(1,4) = get_value('ptc_twiss ','re14 ')
1475	re(1,5) = get_value('ptc_twiss ','re16 ')
1476	re(1,6) = get_value('ptc_twiss ','re15 ')
1477	re(2,1) = get_value('ptc_twiss ','re21 ')
1478	re(2,2) = get_value('ptc_twiss ','re22 ')
1479	re(2,3) = get_value('ptc_twiss ','re23 ')
1480	re(2,4) = get_value('ptc_twiss ','re24 ')
1481	re(2,5) = get_value('ptc_twiss ','re26 ')
1482	re(2,6) = get_value('ptc_twiss ','re25 ')
1483	re(3,1) = get_value('ptc_twiss ','re31 ')
1484	re(3,2) = get_value('ptc_twiss ','re32 ')
1485	re(3,3) = get_value('ptc_twiss ','re33 ')
1486	re(3,4) = get_value('ptc_twiss ','re34 ')
1487	re(3,5) = get_value('ptc_twiss ','re36 ')
1488	re(3,6) = get_value('ptc_twiss ','re35 ')
1489	re(4,1) = get_value('ptc_twiss ','re41 ')
1490	re(4,2) = get_value('ptc_twiss ','re42 ')
1491	re(4,3) = get_value('ptc_twiss ','re43 ')
1492	re(4,4) = get_value('ptc_twiss ','re44 ')
1493	re(4,5) = get_value('ptc_twiss ','re46 ')
1494	re(4,6) = get_value('ptc_twiss ','re45 ')
1495	re(5,1) = get_value('ptc_twiss ','re61 ')
1496	re(5,2) = get_value('ptc_twiss ','re62 ')
1497	re(5,3) = get_value('ptc_twiss ','re63 ')
1498	re(5,4) = get_value('ptc_twiss ','re64 ')
1499	re(5,5) = get_value('ptc_twiss ','re66 ')
1500	re(5,6) = get_value('ptc_twiss ','re65 ')
1501	re(6,1) = get_value('ptc_twiss ','re51 ')
1502	re(6,2) = get_value('ptc_twiss ','re52 ')
1503	re(6,3) = get_value('ptc_twiss ','re53 ')
1504	re(6,4) = get_value('ptc_twiss ','re54 ')
1505	re(6,5) = get_value('ptc_twiss ','re56 ')
1506	re(6,6) = get_value('ptc_twiss ','re55 ')
1507
1508	end subroutine readrematrix
1509	!_________________________________________________________________
1510
1511	subroutine readreforbit
1512	!reads covariance
1513	implicit none
1514	x(:)=zero
1515	x(1)=get_value('ptc_twiss ','x ')
1516	x(2)=get_value('ptc_twiss ','px ')
1517	x(3)=get_value('ptc_twiss ','y ')
1518	x(4)=get_value('ptc_twiss ','py ')
1519	x(5)=get_value('ptc_twiss ','pt ')
1520	x(6)=get_value('ptc_twiss ','t ')
1521	end subroutine readreforbit
1522	!_________________________________________________________________
1523
1524	subroutine readinitialdistrib
1525	!reads covariance matrix of the initial distribution
1526	implicit none
1527	include 'madx_ptc_distrib.inc'
1528	type(taylor) ht
1529	type(pbfield) h
1530	type(damap) id
1531	type(normalform) norm
1532	real(dp) lam
1533	integer nd,nd_m
1534	logical fake_3
1535	integer jc(6)
1536	type(real_8) yy(6)
1537	integer :: dodo = 0
1538	real(dp) x(6)
1539
1540	if(dodo==1) then
1541	x=zero
1542	call find_orbit(my_ring,x,1,default,c_1d_7)
1543	write(6,*) x
1544	call init(default,1,0,berz)
1545	call alloc(yy)
1546	call alloc(id)
1547	id=1
1548	yy=x+id
1549	call track(my_ring,yy,1,default)
1550	call print(yy,6)
1551	stop 999
1552	endif
1553
1554	jc(1)=2
1555	jc(2)=1
1556	jc(3)=4
1557	jc(4)=3
1558	jc(5)=6
1559	jc(6)=5
1560
1561	print*, "We are at initialization with moments of the distribution"
1562
1563	call readrematrix()
1564
1565	print*, re(1,:)
1566	print*, re(2,:)
1567	print*, re(3,:)
1568	print*, re(4,:)
1569	print*, re(5,:)
1570	print*, re(6,:)
1571
1572	nd=2
1573	if(icase==6) nd=3
1574	nd_m=nd
1575	fake_3=.false.
1576	if (getdistrtype(3) /= distr_gauss.and.nd==3) then
1577	!here we have flat in delta
1578	nd_m=2
1579	fake_3=.true.
1580	endif
1581
1582	call init(2,nd,0,0)
1583
1584	call alloc(ht)
1585	call alloc(h)
1586	call alloc(id)
1587	call alloc(norm)
1588
1589	do i = 1,nd_m*2
1590	do ii = 1,nd_m*2
1591	ht=ht+re(i,ii)(-1)(ii+i)(1.0_dp.mono.jc(i))*(1.0_dp.mono.jc(ii))
1592	enddo
1593	enddo
1594	ht=-ht*pi
1595
1596	lam=ten*full_abs(ht)
1597	ht=ht/lam
1598	if(fake_3) then !1959 is the yearh of birth of Etienne (just a number)
1599	ht=ht+(0.1959e0_dp.mono.'000020')+(0.1959e0_dp.mono.'000002')
1600	endif
1601	h=ht
1602	id=1
1603	id=texp(h,id)
1604	norm=id
1605	emi=zero
1606	emi(1:nd_m)=norm%tune(1:nd_m)*lam
1607	re=norm%a_t
1608
1609
1610	do i = 1,c_%nd2
1611	!call print(norm%a_t%v(i),6)
1612	unimap(i) = norm%a_t%v(i)
1613	enddo
1614
1615	! re=id
1616	call setemittances(emi(1),emi(2),emi(3))
1617
1618	call kill(h)
1619	call kill(ht)
1620	call kill(id)
1621	call kill(norm)
1622
1623
1624	end subroutine readinitialdistrib
1625	!_________________________________________________________________
1626
1627	subroutine readmatrixfromtable ! 26 april 2010: changed this routine
1628	! because the format of the map_table has now changed to contain all terms,
1629	! and not only the zeroth and first order ones...
1630	implicit none
1631	integer :: double_from_table_row, table_length
1632	! following added 26 april 2010
1633	integer :: order, nrows
1634	integer :: nx, nxp, ny, nyp, ndeltap, nt, index
1635	real(dp):: coeff
1636	!character(6) :: selector
1637	integer, dimension(6) :: jj ! 3 may 2010
1638
1639	order = get_value('ptc_twiss ','no ')
1640
1641	row = 1 ! starts at one
1642
1643	nrows = table_length("map_table ")
1644
1645	do while(row .le. nrows) ! k=0 when read okay. k=-3 when the table has no row
1646	k = double_from_table_row("map_table ","coef ", row,doublenum)
1647	!write(0,*) 'k=',k
1648	!write(0,*) 'coef=',doublenum
1649	coeff=doublenum
1650	k = double_from_table_row("map_table ","n_vector ",row,doublenum)
1651	index = int(doublenum)
1652	k = double_from_table_row("map_table ","nx ",row,doublenum)
1653	nx = int(doublenum)
1654	!write(0,*) 'index=',index
1655	!write(0,*) 'nx=',nx
1656	k = double_from_table_row("map_table ","nxp ",row,doublenum)
1657	nxp = int(doublenum)
1658	!write(0,*) 'nxp=',nxp
1659	k = double_from_table_row("map_table ","ny ",row,doublenum)
1660	ny = int(doublenum)
1661	!write(0,*) 'ny=',ny
1662	k = double_from_table_row("map_table ","nyp ",row,doublenum)
1663	nyp = int(doublenum)
1664	!write(0,*) 'nyp=',nyp
1665	k = double_from_table_row("map_table ","ndeltap ",row,doublenum)
1666	ndeltap = int(doublenum)
1667	!write(0,*) 'ndeltap=',ndeltap
1668	k = double_from_table_row("map_table ","nt ",row,doublenum)
1669	nt = int(doublenum)
1670	!write(0,*) 'nt=',nt
1671
1672	! y(index)%T.sub.j=coeff ! are we able to do this? NO
1673
1674	! code proposed by piotr to achieve the equivalent of y(1).sub.'12345'=4.0
1675	!oldv = Y(1).sub.'12345'
1676	!newtoset = (4 - oldv).mono.'12345'
1677	!Y(1) = Y(1) + newtoset
1678
1679	! code proposed by etienne to achieve the same
1680	!call pok(y(1),j,4.d0)
1681
1682	if (k.eq.0) then
1683	jj(1)=nx
1684	jj(2)=nxp
1685	jj(3)=ny
1686	jj(4)=nyp
1687	jj(5)=ndeltap
1688	jj(6)=nt
1689	call pok(y(index)%T,jj,coeff)
1690	! the following gives the same result as the above
1691	!oldv = y(index).sub.jj
1692	!newtoset = (coeff - oldv).mono.jj ! mono for monomial
1693	!y(index)%t=y(index)%t+newtoset
1694	! failed to compile the following work in one line
1695	!Y(index) = Y(index) + ((coeff-(Y(index).sub.jj)).mono.jj)
1696	endif
1697
1698
1699
1700
1701	row = row+1
1702
1703	enddo
1704
1705	! call daprint(y,28) ! to be compared with fort.18 created by ptc_normal
1706
1707	call maptoascript()
1708
1709
1710	end subroutine readmatrixfromtable
1711
1712	!_________________________________________________________________
1713	subroutine readinitialmap ! from fort.18 file
1714	implicit none
1715	type(damap) :: map
1716	! call readMapFromFort18(y)
1717	call alloc(map)
1718	call dainput(map,18)
1719	y = map
1720	call kill(map)
1721	call maptoascript()
1722	call reademittance()
1723	end subroutine readinitialmap
1724	!_________________________________________________________________
1725
1726	!_________________________________________________________________
1727
1728	subroutine readinitialmatrix
1729	!reads initial map elements from MAD-X ptc_twiss command parameters
1730	implicit none
1731
1732	call readrematrix() !reads re
1733	call readreforbit() !reads x
1734	call initmapfrommatrix()
1735	call maptoascript()
1736	call reademittance()
1737
1738	end subroutine readinitialmatrix
1739	!_________________________________________________________________
1740
1741	subroutine readinitialascript
1742	!reads initial map elements from MAD-X ptc_twiss command parameters
1743	implicit none
1744	type(damap) :: map
1745	call alloc(map)
1746	call dainput(map,19)
1747	y = x+map
1748	call kill(map)
1749	call reademittance()
1750
1751	end subroutine readinitialascript
1752	!_________________________________________________________________
1753
1754	subroutine reademittance
1755	!initializes Y(6) from re(6,6)
1756	implicit none
1757	real(dp) :: emix,emiy,emiz
1758
1759	emix = get_value('probe ','ex ')
1760	emiy = get_value('probe ','ey ')
1761	emiz = get_value('probe ','et ')
1762
1763	call setemittances(emix,emiy,emiz)
1764
1765	end subroutine reademittance
1766	!_________________________________________________________________
1767
1768	subroutine initmapfrommatrix
1769	!initializes Y(6) from re(6,6)
1770	implicit none
1771	real(dp),dimension(ndim2)::reval,aieval
1772	real(dp),dimension(ndim2,ndim2)::revec,aievec
1773
1774	call liepeek(iia,icoast)
1775	if (getdebug() > 1) then
1776	write (6,'(8a8)') "no","nv","nd","nd2","ndc","ndc2","ndt","ndpt"
1777	write (6,'(8i8)') iia(1),iia(2),iia(3),iia(4),icoast(1),icoast(2),icoast(3),icoast(4)
1778	print*, "c_%npara is ", c_%npara
1779	endif
1780
1781	allocate(j(c_%nv))
1782	j(:)=0
1783	do i = 1,c_%npara
1784	do ii = 1,c_%npara
1785	j(ii)=1
1786	r=re(i,ii)-(y(i)%T.sub.j)
1787	y(i)%T=y(i)%T+(r.mono.j)
1788	j(ii)=0
1789	enddo
1790	enddo
1791	deallocate(j)
1792
1793	call eig6(re,reval,aieval,revec,aievec)
1794	do i=1,iia(4)-icoast(2)
1795	if(abs(reval(i)2+aieval(i)2 -one).gt.c_1d_10) then
1796	call fort_warn("ptc_twiss","Fatal Error: Eigenvalues off the unit circle!")
1797	stop
1798	endif
1799	enddo
1800	end subroutine initmapfrommatrix
1801	!_________________________________________________________________
1802
1803	subroutine maptoascript
1804	!Performes normal form on a map, and plugs A_ in its place
1805	implicit none
1806
1807	call alloc(normal)
1808	normal = y
1809
1810	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
1811	call fort_warn('ptc_twiss: ','Error: DA in twiss got unstable during Normal Form')
1812	call seterrorflag(10,"ptc_twiss ","DA in twiss got unstable during Normal Form");
1813	return
1814	endif
1815
1816	y = x + normal%a_t
1817	call kill(normal)
1818
1819	end subroutine maptoascript
1820	!_________________________________________________________________
1821
1822	subroutine readinitialtwiss(dt)
1823	!Reads initial twiss parameters from MAD-X command
1824	implicit none
1825	real(kind(1d0)) alpha(3),beta(3),disp(4),mu(3)
1826	type(real_8) al(3),be(3),di(4)
1827	type(pol_block_inicond) :: inicondknobs
1828	integer k_system
1829	real(dp) sizept
1830	real(kind(1d0)) emiz
1831	real(dp) dt
1832
1833	beta(1) = get_value('ptc_twiss ','betx ')
1834	beta(2) = get_value('ptc_twiss ','bety ')
1835	beta(3) = get_value('ptc_twiss ','betz ')
1836	alpha(1) = get_value('ptc_twiss ','alfx ')
1837	alpha(2) = get_value('ptc_twiss ','alfy ')
1838	alpha(3) = get_value('ptc_twiss ','alfz ')
1839	disp(1) = get_value('ptc_twiss ','dx ')
1840	disp(2) = get_value('ptc_twiss ','dpx ')
1841	disp(3) = get_value('ptc_twiss ','dy ')
1842	disp(4) = get_value('ptc_twiss ','dpy ')
1843	mu(1) = get_value('ptc_twiss ','mux ')
1844	mu(2) = get_value('ptc_twiss ','muy ')
1845	mu(3) = get_value('ptc_twiss ','muz ')
1846
1847	if (getdebug() > 1) then
1848	print*,"TW: Ini betas ",beta
1849	print*,"TW: Ini alfas ",alpha
1850	endif
1851
1852	if (c_%nd == 3) then
1853	if (beta(3) <= zero) then
1854	call fort_warn("ptc_twiss","Fatal Error: 6D requested and betz is smaller then or equal to 0!")
1855	stop
1856	endif
1857
1858	beta(3) = (one+alpha(3)**2)/beta(3)
1859	alpha(3) =-alpha(3)
1860
1861	endif
1862
1863
1864	x(:)=zero
1865	x(1)=get_value('ptc_twiss ','x ')
1866	x(2)=get_value('ptc_twiss ','px ')
1867	x(3)=get_value('ptc_twiss ','y ')
1868	x(4)=get_value('ptc_twiss ','py ')
1869	! x(5)=get_value('ptc_twiss ','t ')
1870	! x(6)=get_value('ptc_twiss ','pt ')
1871	x(6)=get_value('ptc_twiss ','t ')
1872	x(5)=get_value('ptc_twiss ','pt ')
1873	!frs plug deltap
1874	if(icase.eq.5) x(5) = x(5) + dt
1875
1876
1877	if (getdebug() > 0 ) then
1878	CALL write_closed_orbit(icase,x)
1879	endif
1880
1881	call reademittance()
1882
1883	!Here we initialize Y(6)
1884
1885	call alloc(be); call alloc(al); call alloc(di)
1886
1887	! code to power knows
1888	!
1889
1890
1891	do i=1,c_%nd
1892	be(i)=beta(i)
1893	al(i)=alpha(i)
1894	enddo
1895
1896	do i=1,4
1897	di(i)=disp(i)
1898	enddo
1899
1900	if (getnknobis() > 0) then
1901	!POWER KNOBS
1902	c_%knob = my_true
1903
1904	k_system= (c_%npara - c_%nd2) + getnknobsm()
1905	inicondknobs = getknobinicond()
1906
1907	do i=1,c_%nd
1908
1909	if(inicondknobs%beta(i)/=0) then
1910	if (getdebug() > 1) then
1911	print*,"Beta ",i," is knob no. ", inicondknobs%beta(i)
1912	endif
1913	call make_it_knob(be(i),k_system+inicondknobs%beta(i))
1914	endif
1915
1916	if(inicondknobs%alfa(i)/=0) then
1917	if (getdebug() > 1) then
1918	print*,"Alfa ",i," is knob no. ", inicondknobs%alfa(i)
1919	endif
1920	call make_it_knob(al(i),k_system+inicondknobs%alfa(i))
1921	endif
1922	enddo
1923
1924	do i=1,4
1925	if(inicondknobs%dispersion(i)/=0) then
1926	if (getdebug() > 1) then
1927	print*,"Dispersion ",i," is knob no. ", inicondknobs%dispersion(i)
1928	endif
1929	call make_it_knob(di(i),k_system+inicondknobs%dispersion(i))
1930	endif
1931	enddo
1932
1933	endif
1934
1935
1936	y=x
1937
1938	do i=1,c_%nd
1939	! print*, " ", i, beta(i)
1940	! call print(y(2*i-1),6)
1941	! call print(y(2*i ),6)
1942
1943	y(2i-1)= x(2i-1) + sqrt(be(i)) * morph((one.mono.(2*i-1)) )
1944	y(2i)= x(2i) + one/sqrt(be(i)) * &
1945	(morph( (one.mono.(2i)) )-(al(i)) morph((one.mono.(2*i-1))))
1946
1947	! call print(y(2*i-1),6)
1948	! call print(y(2*i ),6)
1949	enddo
1950
1951
1952
1953	!--moments--!
1954	if( (c_%npara==5) ) then
1955
1956	!print*, "c_%npara ",c_%npara, " c_%ndpt ", c_%ndpt
1957
1958	if ( (beta(3) .gt. zero) .and. (c_%ndpt/=0) ) then
1959
1960	!Option one: sigma(5) is sqrt of emittance as in other two dimensions
1961	! print*, "Init X5 with betaz ", beta(3)
1962
1963	if (c_%nd < 3) then !otherwise it was already done, to be cleaned cause it is ugly and bug prone
1964	beta(3) = (one+alpha(3)**2)/beta(3)
1965	alpha(3) =-alpha(3)
1966	endif
1967
1968	y(5) = x(5) + sqrt( beta(3) )*morph((one.mono.5))
1969	y(6)= x(6) + one/sqrt(beta(3)) * (morph( (one.mono.6) )-(alpha(3)) * morph(one.mono.5))
1970
1971	emiz = get_value('probe ','et ')
1972	if ( emiz .le. 0 ) then
1973	sizept = get_value('probe ','sige ')
1974	emiz = sizept/sqrt(beta(3))
1975	! print*, "Calculated Emittance ", emiz
1976	else
1977	emiz = sqrt(emiz)
1978	! print*, "Read Emittance ", emiz
1979	endif
1980
1981	call setsigma(5, emiz)
1982	call setsigma(6, emiz)
1983
1984	else
1985	!by default we have no knowledge about longitudinal phase space, so init dp/p to ident
1986	! print*, "Init X5 with ONE"
1987	!frs we need here the initial value of pt and t should not hurt
1988	y(5) = x(5) + morph((one.mono.5))
1989	y(6) = x(6) + morph((one.mono.5))
1990	call setsigma(5, get_value('probe ','sige '))
1991	call setsigma(6, get_value('probe ','sigt '))
1992	endif
1993
1994	endif
1995
1996	if ( (icase .gt. 5) .and. (get_value('probe ','et ') .le. 0) ) then !6 and 56
1997	!beta(3) is converted to gamma already (in 3rd coord the canonical planes are swapped)
1998	emiz = get_value('probe ','sige ')/sqrt(beta(3))
1999	print*, "icase=",icase," nd=",c_%nd, "sigma(5)=", emiz
2000	call setsigma(5, emiz)
2001	call setsigma(6, emiz)
2002	endif
2003
2004
2005
2006	if(icase/=4) then
2007	do i=1,4
2008	y(i)= y(i) + di(i) * morph((one.mono.5))
2009	enddo
2010	endif
2011
2012
2013	if ( getdebug() > 2) then
2014	print*," Read the following BETA0 block in module ptc_twiss"
2015	print*," Twiss parameters:"
2016	write (6,'(6a8)') "betx","alfx","bety","alfy","betz","alfz"
2017	write (6,'(6f8.4)') beta(1), alpha(1), beta(2), alpha(2), beta(3), alpha(3)
2018	write (6,'(4a8)') "dx","dpx","dy","dpy"
2019	write (6,'(4f8.4)') disp
2020	write (6,'(3a8)') "mux","muy","muz"
2021	write (6,'(3f8.4)') mu
2022
2023	print*," Track:"
2024	write(6,'(6f8.4)') x
2025	endif
2026
2027
2028	end subroutine readinitialtwiss
2029	!____________________________________________________________________________________________
2030
2031	subroutine onePassSummary(oneTurnMap,state,suml)
2032
2033	implicit none
2034	type(real_8),target :: oneTurnMap(6)
2035	real(dp), target :: state(6) ! six-dimensional phase-space state (usually referred-to as 'x')
2036	real(dp) :: suml ! cumulative length along the ring
2037	real(dp) :: rdp_mmilion ! float with zero (0)
2038	real(dp) :: deltap ! float with zero (0)
2039
2040	rdp_mmilion= -1e6;
2041
2042	call double_to_table_curr( summary_table_name, 'length ', suml ) ! total length of the machine
2043
2044
2045	call double_to_table_curr( summary_table_name, 'alpha_c ', rdp_mmilion ) ! momemtum compaction factor
2046	call double_to_table_curr( summary_table_name, 'alpha_c_p ', rdp_mmilion) ! derivative w.r.t delta-p/p
2047	call double_to_table_curr( summary_table_name, 'alpha_c_p2 ', rdp_mmilion) ! 2nd order derivative
2048	call double_to_table_curr( summary_table_name, 'alpha_c_p3 ', rdp_mmilion) ! 3rd order derivative
2049	call double_to_table_curr( summary_table_name, 'eta_c ', rdp_mmilion) ! associated phase-slip factor
2050	call double_to_table_curr( summary_table_name, 'gamma_tr ', rdp_mmilion) ! associated transition energy
2051
2052	call double_to_table_curr( summary_table_name, 'q1 ', tw%mu(1))
2053	call double_to_table_curr( summary_table_name, 'q2 ', tw%mu(2))
2054	call double_to_table_curr( summary_table_name, 'dq1 ', rdp_mmilion)
2055	call double_to_table_curr( summary_table_name, 'dq2 ', rdp_mmilion)
2056
2057	call double_to_table_curr( summary_table_name, 'qs ', rdp_mmilion)
2058	call double_to_table_curr( summary_table_name, 'beta_x_min ', rdp_mmilion)
2059	call double_to_table_curr( summary_table_name, 'beta_x_max ', rdp_mmilion)
2060	call double_to_table_curr( summary_table_name, 'beta_y_min ', rdp_mmilion)
2061	call double_to_table_curr( summary_table_name, 'beta_y_max ', rdp_mmilion)
2062
2063	deltap = get_value('ptc_twiss ','deltap ')
2064	call double_to_table_curr( summary_table_name, 'deltap ', deltap)
2065
2066
2067	call double_to_table_curr( summary_table_name,'orbit_x ', rdp_mmilion)
2068	call double_to_table_curr( summary_table_name,'orbit_px ', rdp_mmilion)
2069	call double_to_table_curr( summary_table_name,'orbit_y ', rdp_mmilion)
2070	call double_to_table_curr( summary_table_name,'orbit_py ', rdp_mmilion)
2071
2072	call double_to_table_curr( summary_table_name,'xcorms ', rdp_mmilion)
2073	call double_to_table_curr( summary_table_name,'ycorms ', rdp_mmilion)
2074	call double_to_table_curr( summary_table_name,'pxcorms ', rdp_mmilion)
2075	call double_to_table_curr( summary_table_name,'pycorms ', rdp_mmilion)
2076
2077	call double_to_table_curr( summary_table_name,'xcomax ', rdp_mmilion)
2078	call double_to_table_curr( summary_table_name,'ycomax ', rdp_mmilion)
2079	call double_to_table_curr( summary_table_name,'pxcomax ', rdp_mmilion)
2080	call double_to_table_curr( summary_table_name,'pycomax ', rdp_mmilion)
2081
2082	call double_to_table_curr( summary_table_name,'orbit_pt ', rdp_mmilion)
2083	call double_to_table_curr( summary_table_name,'orbit_-cT ', rdp_mmilion)
2084
2085	call augment_count( summary_table_name ); ! only one row actually...
2086
2087	end subroutine onePassSummary
2088	! jluc
2089	! compute momemtum-compaction factor in the same fashion it is carried-out in twiss.F
2090
2091	subroutine oneTurnSummary(isRing,oneTurnMap,state,suml)
2092
2093	implicit none
2094	logical :: isRing ! true for rings, false for lines
2095	type(real_8),target :: oneTurnMap(6)
2096	real(dp), target :: state(6) ! six-dimensional phase-space state (usually referred-to as 'x')
2097	real(dp) :: suml ! cumulative length along the ring
2098	type(fibre), pointer :: fibrePtr
2099	real(dp) :: alpha_c, eta_c ! momentum-compaction factor & phase-slip factor
2100	real(dp) :: alpha_c_p ! first order derivative w.r.t delta-p/p
2101	real(dp) :: alpha_c_p2 ! second order derivative w.r.t delta-p/p
2102	real(dp) :: alpha_c_p3 ! third order derivative w.r.t delta-p/p
2103	real(dp) :: gamma_tr ! gamma_transition, or "transition energy" above which the particles' arrival time
2104	! with respect to other particles is determined by its path length instead of by its velocity
2105	real(dp) :: deltap
2106	real(dp) :: betaRelativistic, gammaRelativistic
2107	real(dp) :: fractionalTunes(ndim)
2108	real(dp) :: chromaticities(2)
2109	integer :: i,j
2110	real(dp) :: sd ! as in twiss.F
2111	type(real_8) :: theAscript(6) ! used here to compute dispersion's derivatives
2112	type(damap) :: yy ! added on November 6th to retreive momemtum compaction without differentiating the formula
2113	integer, dimension(6,6) :: coeffSelector = &
2114	reshape( (/1,0,0,0,0,0, &
2115	0,1,0,0,0,0, &
2116	0,0,1,0,0,0, &
2117	0,0,0,1,0,0, &
2118	0,0,0,0,1,0, &
2119	0,0,0,0,0,1/), (/6,6/))
2120	type(normalform) theNormalForm
2121	real(dp) :: dispersion(4)
2122	real(dp) :: dispersion_p(4) ! derivative of the dispersion w.r.t delta-p/p
2123	integer :: debugFiles
2124	integer :: icase
2125	integer :: order
2126	real(dp) :: rdp_mmilion
2127	rdp_mmilion= -1e6;
2128
2129	order = get_value('ptc_twiss ', 'no ')
2130
2131	! should end-up gracefully here in case the topology of the lattice is not those of a closed-ring
2132
2133	debugFiles = 0 ! set it to one and fort.21, fort.22 and fort.23 are created
2134
2135
2136	! 2. retreive the relativistic parameters beta and gamma
2137	! (beta=v/c, gamma=E/mc^2 and gamma=1/sqrt(1-beta^2))
2138	betaRelativistic = get_value('probe ','beta ');
2139	gammaRelativistic = get_value('probe ','gamma ');
2140	icase = get_value('ptc_twiss ','icase ') ! mind the trailing space
2141
2142	! now retrieve the one-turn map's coefficients
2143	if (debugFiles .eq. 1) then
2144	do i=1,6
2145	do j=1,6
2146	write(21,*) "r(",i,j,")=",oneTurnMap(i).sub.coeffSelector(j,:)
2147	enddo
2148	enddo
2149	endif
2150
2151	! 4. retreive the dispersion coefficients
2152	! (may be the coefficient of delta of the map?)
2153	! decompose the map via a normal form to get the dispersion...
2154	call alloc(theNormalForm)
2155	theNormalForm = oneTurnMap
2156
2157
2158	if (debugFiles .eq. 1) then
2159	call daprint(theNormalForm%A1,23) ! supposed to print dispersion's first and higher orders
2160	! according to h_definition.f90: type normalform contains A1 as dispersion
2161	! (would need to go through DHDJ to get the tune...)
2162	endif
2163	! first order dispersions !?
2164	! (at least checked that these values match those computed in twiss.F)
2165	dispersion(1) = theNormalForm%A1%v(1).sub.'000010'
2166	dispersion(2) = theNormalForm%A1%v(2).sub.'000010'
2167	dispersion(3) = theNormalForm%A1%v(3).sub.'000010'
2168	dispersion(4) = theNormalForm%A1%v(4).sub.'000010'
2169
2170	if (debugFiles .eq. 1) then
2171	do i=1,4
2172	write(21,*) "dispersion(",i,")=", dispersion(i)
2173	enddo
2174	endif
2175
2176
2177	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2178	!!!!!!!!!!!! !!!!!!!!!!!!!!!!
2179	!!!!!!!!!!!! COMPACTION FACTOR, ITS HIGHER ORDERS !!!!!!!!!!!!!!!!
2180	!!!!!!!!!!!! GAMMA TRANSITION, SLIP FACTOR !!!!!!!!!!!!!!!!
2181	!!!!!!!!!!!! !!!!!!!!!!!!!!!!
2182	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2183
2184
2185
2186	! in 4D we have no knowledge of longitudinal dynamics, it is not calculated
2187	! the same 5D with time=false
2188
2189	eta_c = rdp_mmilion
2190	alpha_c = rdp_mmilion
2191	gamma_tr = rdp_mmilion
2192
2193	alpha_c_p = rdp_mmilion
2194	alpha_c_p2 = rdp_mmilion
2195	alpha_c_p3 = rdp_mmilion
2196
2197	!call daprint(theNormalForm%dhdj,6)
2198
2199
2200	if ( (icase.gt.4) .and. (default%time) ) then
2201
2202	! Here R56 is dT/ddelta
2203	! 5. apply formulas from twiss.F:
2204	!sd = r(5,6)+r(5,1)disp(1)+...+r(5,4)disp(4)
2205	!print*,"ALPHA_C, GAMMA TR : TIME ON"
2206
2207	sd = - oneTurnMap(6).sub.coeffSelector(5,:) ! 5/6 swap MADX/PTC
2208	!print*,'sd(0)',sd
2209	do i=1,4
2210	!print*, 'Disp',i,'=', dispersion(i)
2211	sd = sd - (oneTurnMap(6).sub.coeffSelector(i,:))*dispersion(i)
2212	enddo
2213	!print*,'sd(f)',sd
2214
2215	eta_c = -sd * betaRelativistic**2 / suml
2216	alpha_c = one / gammaRelativistic**2 + eta_c
2217	gamma_tr = one / sqrt(alpha_c)
2218
2219	elseif( (icase.eq.5) .and. (default%time .eqv. .false.) ) then
2220
2221	! Here R56 is dL/ddelta
2222	! so we get alpha_c first from transfer matrix
2223
2224	sd = +1.0*(oneTurnMap(6).sub.coeffSelector(5,:)) ! 5/6 swap MADX/PTC
2225	!print*,'sd(0)',sd
2226	do i=1,4
2227	!print*, 'Disp',i,'=', dispersion(i)
2228	sd = sd + (oneTurnMap(6).sub.coeffSelector(i,:))*dispersion(i)
2229	enddo
2230	!print*,'sd(f)',sd
2231	alpha_c = sd/suml
2232	eta_c = alpha_c - one / gammaRelativistic**2
2233	gamma_tr = one / sqrt(alpha_c)
2234
2235	elseif( (icase.eq.56) .and. (default%time .eqv. .false.) ) then
2236
2237	!print*,"ALPHA_C, GAMMA TR : 56D TIME OFF"
2238
2239
2240	call alloc(yy)
2241	do i=1,c_%nd2 ! c_%nd2 is 6 when icase is 56 or 6 (but 4 when icase=5)
2242	yy%v(i) = oneTurnMap(i)%t
2243	enddo
2244	yy = theNormalForm%A1*(-1)yy*theNormalForm%A1 ! takes away all dispersion dependency
2245	!write(0,*) 'for yy, c_%nd2 is ',c_%nd2 ! 0 is stderr
2246	alpha_c = (yy%v(6).sub.'000010')/suml
2247	gamma_tr = one / sqrt(alpha_c)! overwrite the value obtained from the Twiss formula
2248	eta_c = alpha_c - one / gammaRelativistic**2
2249
2250	alpha_c_p = 2.0*(yy%v(6).sub.'000020')/suml
2251
2252	if (order.ge.3) then
2253	alpha_c_p2 = 3.02.0(yy%v(6).sub.'000030')/suml
2254	endif
2255
2256	if (order.ge.4) then
2257	alpha_c_p3 = 4.03.02.0*(yy%v(6).sub.'000040')/suml
2258	endif
2259
2260
2261	call kill(yy)
2262
2263	endif
2264
2265	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2266	!!!! HIGHER ORDERS IN dP/P !!!!
2267	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2268
2269	if ( (icase.eq.5) .and. (default%time .eqv. .true.) ) then
2270
2271	!print*,"ALPHA_C dp/dp derivatives : 5D TIME ON"
2272
2273	if (getdebug() > 2) then
2274	call kanalnummer(mf1)
2275	open(unit=mf1,file='oneTurnMap.5dt.txt')
2276	call daprint(oneTurnMap,mf1)
2277	close(mf1)
2278	endif
2279
2280	! compute delta-p/p dependency of alpha_c
2281	! first order derivatives of the dispersions
2282	! assuming icase=5
2283	!call daprint(oneTurnMap,88)
2284
2285	! always assume time=true, so that the fifth phase-space variable is deltap instead of pt!
2286	! otherwise should issue a warning or an error!
2287
2288	call alloc(theAscript)
2289	! proceed slightly differently to retreive the dispersion's derivatives
2290	! (so far only managed to get the dispersion on the normal form, but not its derivatives)
2291	theAscript = state + theNormalForm%a_t
2292
2293	! note: should reuse information computed previously than redo this...
2294	! if icase=5, Taylor series expansion disp = disp + delta_p
2295	do i=1,4
2296	! apparently, we are up by a factor two
2297	dispersion_p(i) = 2.0*(theAscript(i)%t.sub.'000020') ! as usual in this file
2298	enddo
2299
2300	! compute derivative of the formula for eta_c
2301
2302	sd = 0.0
2303	! first dr65/ddeltap
2304	sd = sd - 1.0(oneTurnMap(6).sub.'100010')dispersion(1) &
2305	& - 1.0(oneTurnMap(6).sub.'010010')dispersion(2) &
2306	& - 1.0(oneTurnMap(6).sub.'001010')dispersion(3) &
2307	& - 1.0(oneTurnMap(6).sub.'000110')dispersion(4) &
2308	& - 2.0*(oneTurnMap(6).sub.'000020') &
2309	! new
2310	& - 1.0(oneTurnMap(6).sub.'000011')(oneTurnMap(6).sub.'000001') &
2311	! then terms dr6i/ddeltap*dispersion(i)
2312	& - 2.0(oneTurnMap(6).sub.'200000')dispersion(1)**2 &
2313	& - (oneTurnMap(6).sub.'110000')dispersion(1)dispersion(2) &
2314	& - (oneTurnMap(6).sub.'101000')dispersion(1)dispersion(3) &
2315	& - (oneTurnMap(6).sub.'100100')dispersion(1)dispersion(4) &
2316	! new
2317	& - (oneTurnMap(6).sub.'100010')*dispersion(1) &
2318	! dr62/ddeltap*disperion(2)
2319	& - (oneTurnMap(6).sub.'110000')dispersion(1)dispersion(2) &
2320	& - 2.0(oneTurnMap(6).sub.'020000')dispersion(2)**2 &
2321	& - (oneTurnMap(6).sub.'011000')dispersion(2)dispersion(3) &
2322	& - (oneTurnMap(6).sub.'010100')dispersion(2)dispersion(4) &
2323	! new
2324	& - (oneTurnMap(6).sub.'010010')*dispersion(2) &
2325	! dr63/ddeltap*dispersion(3)
2326	& - (oneTurnMap(6).sub.'101000')dispersion(1)dispersion(3) &
2327	& - (oneTurnMap(6).sub.'011000')dispersion(2)dispersion(3) &
2328	& - 2.0(oneTurnMap(6).sub.'002000')dispersion(3)**2 &
2329	& - (oneTurnMap(6).sub.'001100')dispersion(3)dispersion(4) &
2330	! new
2331	& - 1(oneTurnMap(6).sub.'001010')dispersion(3) &
2332	! dr64/ddeltap*dispersion(4)
2333	& - (oneTurnMap(6).sub.'100100')dispersion(1)dispersion(4) &
2334	& - (oneTurnMap(6).sub.'010100')dispersion(2)dispersion(4) &
2335	& - (oneTurnMap(6).sub.'001100')dispersion(3)dispersion(4) &
2336	& - 2.0(oneTurnMap(6).sub.'000200')dispersion(4)**2 &
2337	! new
2338	& - 1(oneTurnMap(6).sub.'000110')dispersion(4)
2339
2340	! terms involving derivatives of the dispersions
2341	do i=1,4
2342	! print*, 'sd=',sd,' disp=',dispersion_p(i)
2343	sd = sd - (oneTurnMap(6).sub.coeffSelector(i,:))*dispersion_p(i)
2344	enddo
2345
2346	! print*, 'sd=',sd
2347	alpha_c_p = sd * (betaRelativistic**2) / suml
2348
2349	! eventually, one could differentiate the above formula to obtain alpha_c_p2
2350	! but for the time-being, expect icase to be 56 to compute alpha_c_p2 and alpha_c_p3.
2351
2352	alpha_c_p2 = rdp_mmilion
2353	alpha_c_p3 = rdp_mmilion
2354
2355
2356	call kill(theAscript)
2357
2358
2359
2360	elseif ( (icase.eq.56) .and. (default%time .eqv. .true.) ) then ! here one may obtain the pathlength derivatives from the map
2361
2362	!print*,"ALPHA_C dp/dp derivatives : 56D TIME ON"
2363	if (getdebug() > 2) then
2364	call kanalnummer(mf1)
2365	open(unit=mf1,file='oneTurnMap.56dt.txt')
2366	call daprint(oneTurnMap,mf1)
2367	close(mf1)
2368	endif
2369
2370	call alloc(yy)
2371	do i=1,c_%nd2 ! c_%nd2 is 6 when icase is 56 or 6 (but 4 when icase=5)
2372	yy%v(i) = oneTurnMap(i)%t
2373	enddo
2374	yy = theNormalForm%A1*(-1)yy*theNormalForm%A1 ! takes away all dispersion dependency
2375
2376	alpha_c_p = -2.0*(yy%v(6).sub.'000020')/suml
2377
2378	if (order.ge.3) then
2379	alpha_c_p2 = -3.02.0(yy%v(6).sub.'000030')/suml
2380	endif
2381
2382	if (order.ge.4) then
2383	alpha_c_p3 = -4.03.02.0*(yy%v(6).sub.'000040')/suml
2384	endif
2385
2386	call kill(yy)
2387
2388	endif
2389
2390	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2391	!!!!!!!!!!!! !!!!!!!!!!!!!!!!
2392	!!!!!!!!!!!! END OF COMPACTION FACTOR !!!!!!!!!!!!!!!!
2393	!!!!!!!!!!!! !!!!!!!!!!!!!!!!
2394	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2395
2396
2397
2398	! also output the tune ...
2399	fractionalTunes = theNormalForm%tune
2400	! the above is exactly equivalent to the following two lines (i.e. returns frac.tune)
2401	fractionalTunes(1) = theNormalForm%DHDJ%v(1).sub.'0000' ! as in So_fitting.f90
2402	fractionalTunes(2) = theNormalForm%DHDJ%v(2).sub.'0000' ! as in So_fitting.f90
2403	! 26 november 2009
2404	if (icase.eq.6) then
2405	fractionalTunes(3) = - theNormalForm%DHDJ%v(3).sub.'0000' ! to enter here icase must be 6 but not only:
2406	! there must be a cavity otherwise icase is set internally to 56 by my_state in ptc_module.f90
2407	!write(0,*) 'fractionTunes(3)=', fractionalTunes(3)
2408	! in the above, inserted minus sign to match the 'phase' or 'tw%mu(3)'
2409	! computed as atan2(ascript(6).sub.'000010',ascript(6).sub.'000001')/2*pi
2410	else
2411	fractionalTunes(3) = 0.0
2412	!write(0,*) 'nullify fractionalTunes(3), icase=',icase
2413	endif
2414	! Q: is it possible to get the actual total tune, as returned by twiss.F?
2415	! => no, not with a map...
2416
2417	! ... as well as the chromaticities
2418	if (icase.eq.5 .or. icase.eq.56) then
2419	chromaticities(1) = theNormalForm%DHDJ%v(1).sub.'00001' ! as in So_fitting.f90
2420	chromaticities(2) = theNormalForm%DHDJ%v(2).sub.'00001' ! as in So_fitting.f90
2421	! to get chromaticities, went to higher order with above "call init_default(default,2,0)"
2422	else
2423	! if icase = 6, delta_p is a phase-space variable and not an external parameter hence we can't compute chromaticies
2424	chromaticities(1) = 0.0
2425	chromaticities(2) = 0.0
2426	endif
2427
2428	! for debug: check the values by printing the map
2429	if (debugFiles .eq. 1) then
2430	call daprint(oneTurnMap,25) ! prints the one-turn map on file 25
2431	call daprint(theNormalForm%dhdj,26) ! print tunes, chromaticities and anharmonicities
2432	! as done in madx_ptc_normal.f90
2433	endif
2434
2435	call kill(theNormalForm)
2436
2437	call kill(oneTurnMap)
2438
2439
2440	!write(6,*)
2441	!write(6,*) "Momentum compaction (momentum compaction factor & phase-slip factor):"
2442	!write(6,*) "alpha=", alpha_c, "eta=", eta_c
2443	!write(6,*) "fractional tunes=", fractionalTunes
2444	!write(6,*) "chromaticities=", chromaticities
2445
2446	deltap = get_value('ptc_twiss ','deltap ')
2447
2448	! write the data into the ptc_twiss_summary table
2449	! warning: must preserve the order EXACTLY
2450	call double_to_table_curr( summary_table_name, 'length ', suml ) ! total length of the machine
2451	call double_to_table_curr( summary_table_name, 'alpha_c ', alpha_c ) ! momemtum compaction factor
2452	call double_to_table_curr( summary_table_name, 'alpha_c_p ', alpha_c_p) ! derivative w.r.t delta-p/p
2453	call double_to_table_curr( summary_table_name, 'alpha_c_p2 ', alpha_c_p2) ! 2nd order derivative
2454	call double_to_table_curr( summary_table_name, 'alpha_c_p3 ', alpha_c_p3) ! 3rd order derivative
2455	call double_to_table_curr( summary_table_name, 'eta_c ', eta_c ) ! associated phase-slip factor
2456	call double_to_table_curr( summary_table_name, 'gamma_tr ', gamma_tr) ! associated transition energy
2457	call double_to_table_curr( summary_table_name, 'q1 ', fractionalTunes(1))
2458	call double_to_table_curr( summary_table_name, 'q2 ', fractionalTunes(2))
2459	call double_to_table_curr( summary_table_name, 'dq1 ', chromaticities(1))
2460	call double_to_table_curr( summary_table_name, 'dq2 ', chromaticities(2))
2461	! 26 november 2009
2462	call double_to_table_curr( summary_table_name, 'qs ', fractionalTunes(3))
2463	! write the extremas of the Twiss functions
2464	! for the time-being, do not bother about the coupling terms
2465	call double_to_table_curr( summary_table_name, 'beta_x_min ', minBeta(1,1))
2466	call double_to_table_curr( summary_table_name, 'beta_x_max ', maxBeta(1,1))
2467	call double_to_table_curr( summary_table_name, 'beta_y_min ', minBeta(2,2))
2468	call double_to_table_curr( summary_table_name, 'beta_y_max ', maxBeta(2,2))
2469	call double_to_table_curr( summary_table_name, 'deltap ', deltap)
2470	! the 6-d closed orbit
2471	call double_to_table_curr( summary_table_name,'orbit_x ',state(1))
2472	call double_to_table_curr( summary_table_name,'orbit_px ', state(2))
2473	call double_to_table_curr( summary_table_name,'orbit_y ', state(3))
2474	call double_to_table_curr( summary_table_name,'orbit_py ', state(4))
2475	! warning: if 'time=false', the last two phase-space state-variables
2476	! should be deltap/p and path-length respectively
2477	call double_to_table_curr( summary_table_name,'orbit_pt ', state(5))
2478	call double_to_table_curr( summary_table_name,'orbit_-cT ', state(6))
2479
2480	call augment_count( summary_table_name ); ! only one row actually...
2481
2482	end subroutine oneTurnSummary
2483
2484
2485	END subroutine ptc_twiss
2486	!____________________________________________________________________________________________
2487
2488	subroutine computeDeltapDependency(y,s1)
2489	implicit none
2490	type(real_8), intent(in) :: y(6)
2491	type(twiss), intent(inout) :: s1
2492	integer :: k,i
2493	integer :: J(lnv) ! the map's coefficient selector, as usual
2494	integer :: Jderiv(lnv) ! to store the map's coefficient selector of the derivative w.r.t deltap
2495	real(kind(1d0)) :: get_value ! C-function
2496	integer :: no ! order must be at equal to 2 to be able to get terms of the form x*deltap
2497	! required to evaluate the derivatives of Twiss parameters w.r.t deltap
2498	integer :: ndel ! as in subroutine 'equaltwiss'...
2499
2500	! in order to avoid this message, should prevent entering this subroutine
2501	! in case none of the Twiss derivatives is selected...
2502
2503	no = get_value('ptc_twiss ','no ')
2504	if ( no .lt. 2 ) then
2505	call fort_warn('madx_ptc_twiss.f90 <ptc_twiss>:','Order in computeDeltapDependency() is smaller then 2')
2506	print*, "Order is ", no
2507	return
2508	endif
2509
2510	J=0
2511	do i=1,c_%nd ! c_%nd is global variable
2512	do k=1,c_%nd ! the same
2513	J(2*i-1)=1
2514	Jderiv = J ! vector copy
2515	Jderiv(5)=1 ! the delta_p coefficient
2516	s1%beta_p(k,i)= (Y(2k-1).sub.Jderiv)(Y(2k-1).sub.J) + (Y(2k-1).sub.J)(Y(2k-1).sub.Jderiv)
2517	s1%alfa_p(k,i)= -(Y(2k-1).sub.Jderiv)(Y(2k).sub.J) - (Y(2k-1).sub.J)(Y(2k).sub.Jderiv)
2518	s1%gama_p(k,i)= (Y(2k).sub.Jderiv)(Y(2k).sub.J) + (Y(2k).sub.J)(Y(2k).sub.Jderiv)
2519	J(2*i-1)=0
2520	J(2*i)=1
2521	Jderiv = J ! vector copy
2522	Jderiv(5)=1 ! the delta_p coefficient
2523	s1%beta_p(k,i) = s1%beta_p(k,i) &
2524	+ (Y(2k-1).sub.Jderiv)(Y(2k-1).sub.J) + (Y(2k-1).sub.J)(Y(2k-1).sub.Jderiv)
2525	s1%alfa_p(k,i)= s1%alfa_p(k,i) &
2526	- (Y(2k-1).sub.Jderiv)(Y(2k).sub.J) - (Y(2k-1).sub.J)(Y(2k).sub.Jderiv)
2527	s1%gama_p(k,i)= s1%gama_p(k,i) &
2528	+ (Y(2k).sub.Jderiv)(Y(2k).sub.J)+(Y(2k).sub.J)(Y(2k).sub.Jderiv)
2529	J(2*i)=0
2530	enddo
2531	enddo
2532
2533	! the computations above match the following formulas, obtained by derivation of the Twiss parameters using the chain-rule
2534	! beta derivatives w.r.t delta_p
2535	! beta11 = two * (y(1)%t.sub.'100000')(y(1)%t.sub.'100010') + two (y(1)%t.sub.'010000')*(y(1)%t.sub.'010010')
2536	! beta12 = two * (y(1)%t.sub.'001000')(y(1)%t.sub.'001010') + two (y(1)%t.sub.'000100')*(y(1)%t.sub.'000110')
2537	! beta21 = two * (y(3)%t.sub.'100000')(y(3)%t.sub.'100010') + two (y(3)%t.sub.'010000')*(y(3)%t.sub.'010010')
2538	! beta22 = two * (y(3)%t.sub.'001000')(y(3)%t.sub.'001010') + two (y(3)%t.sub.'000100')*(y(3)%t.sub.'000110')
2539	! alpha derivatives w.r.t delta_p
2540	! alfa11 = -((y(1)%t.sub.'100010')(y(2)%t.sub.'100000')+(y(1)%t.sub.'100000')(y(2).sub.'100010')+&
2541	! (y(1)%t.sub.'010010')(y(2)%t.sub.'010000')+(y(1)%t.sub.'010000')(y(2)%t.sub.'010010'))
2542	! alfa12 = -((y(1)%t.sub.'001010')(y(2)%t.sub.'001000')+(y(1)%t.sub.'001000')(y(2)%t.sub.'001010')+&
2543	! (y(1)%t.sub.'000110')(y(2)%t.sub.'000100')+(y(1)%t.sub.'000100')(y(2)%t.sub.'000110'))
2544	! alfa21 = -((y(3)%t.sub.'100010')(y(4)%t.sub.'100000')+(y(3)%t.sub.'100000')(y(4)%t.sub.'100010')+&
2545	! (y(3)%t.sub.'010010')(y(4)%t.sub.'010000')+(y(3)%t.sub.'010000')(y(4)%t.sub.'010010'))
2546	! alfa22 = -((y(3)%t.sub.'001010')(y(4)%t.sub.'001000')+(y(3)%t.sub.'001000')(y(4)%t.sub.'001010')+&
2547	! (y(3)%t.sub.'000110')(y(4)%t.sub.'000100')+(y(3)%t.sub.'000100')(y(4)%t.sub.'000110'))
2548	! gamma derivatives w.r.t delta_p
2549	! gama11 = (y(2)%t.sub.'100010')(y(2)%t.sub.'100000')+(y(2)%t.sub.'100000')(y(2)%t.sub.'100010')+&
2550	! (y(2)%t.sub.'010010')(y(2)%t.sub.'010000')+(y(2)%t.sub.'010000')(y(2)%t.sub.'010010')
2551	! gama12 = (y(2)%t.sub.'001010')(y(2)%t.sub.'001000')+(y(2)%t.sub.'001000')(y(2)%t.sub.'001010')+&
2552	! (y(2)%t.sub.'000110')(y(2)%t.sub.'000100')+(y(2)%t.sub.'000100')(y(2)%t.sub.'000110')
2553	! gama21 = (y(4)%t.sub.'100010')(y(4)%t.sub.'100000')+(y(4)%t.sub.'100000')(y(4)%t.sub.'100010')+&
2554	! (y(4)%t.sub.'010010')(y(4)%t.sub.'010000')+(y(4)%t.sub.'010000')(y(4)%t.sub.'010010')
2555	! gama22 = (y(4)%t.sub.'001010')(y(4)%t.sub.'001000')+(y(4)%t.sub.'001000')(y(4)%t.sub.'001010')+&
2556	! (y(4)%t.sub.'000110')(y(4)%t.sub.'000100')+(y(4)%t.sub.'000100')(y(4)%t.sub.'000110')
2557
2558	! now compute deltap dependencies of the dispersion
2559
2560	! code to be differentiated, as in subroutine 'equaltwiss'
2561	! J=0
2562	! !here ND2=4 and delta is present nd2=6 and delta is a constant
2563	! ! print*,"nv",c_%nv,"nd2",c_%nd2,"np",c_%np,"ndpt",c_%ndpt ,"=>",c_%nv-c_%nd2-c_%np
2564	! if( (c_%npara==5) .or. (c_%ndpt/=0) ) then
2565	! !when there is no cavity it gives us dispersions
2566	! do i=1,4
2567	! lat(0,i,1)=(Y(i)%t.sub.J5)
2568	! enddo
2569	! elseif (c_%nd2 == 6) then
2570	! do i=1,4
2571	! lat(0,i,1) = (Y(i)%t.sub.J5)*(Y(6)%t.sub.J6)
2572	! lat(0,i,1) = lat(0,i,1) + (Y(i)%t.sub.J6)*(Y(5)%t.sub.J5)
2573	! enddo
2574	! else
2575	! do i=1,4
2576	! lat(0,i,1)=zero
2577	! enddo
2578	! endif
2579	!
2580	! ...
2581	!
2582	! !when there is no cavity it gives us dispersions
2583	! do i=1,c_%nd2-2*ndel
2584	! s1%disp(i)=lat(0,i,1)
2585	! enddo
2586
2587
2588
2589	! differentiating the code that computes the dispersion...
2590	ndel=0
2591	! as in subroutine 'equaltwiss'...
2592	if (c_%ndpt/=0) then
2593	ndel=1
2594	endif
2595	if ((c_%npara ==5) .or. (c_%ndpt /= 0)) then
2596	do i=1,c_%nd2-2*ndel ! should it be 1 to 4?
2597	! in 4D+deltap case, coefficients are those of the Taylor series development
2598	! of deltap (factorial)
2599	! disp = disp0 + sigma (1/n!) disp_pn
2600	s1%disp_p(i) = 2.0*(y(i)%t.sub.'000020')
2601	if (no.gt.2) then
2602	s1%disp_p2(i) = 3.02.0(y(i)%t.sub.'000030') ! assume at least order 3 for the map
2603	else
2604	call fort_warn("ptc_twiss ","assume no>=3 for dispersion's 2nd order derivatives w.r.t delta-p")
2605	endif
2606	if (no.gt.3) then
2607	s1%disp_p3(i) = 4.03.02.0*(y(i)%t.sub.'000040')
2608	! assume at least order 4 for the map
2609	else
2610	call fort_warn("ptc_twiss ","assume no>=4 for dispersion's 3rd order derivatives w.r.t delta-p")
2611	endif
2612	enddo
2613	elseif (c_%nd2 ==6) then
2614	do i=1,c_%nd2-2*ndel ! should it be 1 to 4?
2615	! finally never enter here: these differentiation formulas turned-out not to apply to the case
2616	! where deltap is a phase-space state-variable rather than an externally supplied
2617	! parameter.
2618	! u'v+uv'
2619	s1%disp_p(i) = &
2620	2.0(y(i)%t.sub.'000020')(y(6)%t.sub.'000001')+&
2621	(y(i)%t.sub.'000010')*(y(6)%t.sub.'000011')+&
2622	(y(i)%t.sub.'000011')*(y(5)%t.sub.'000010')+&
2623	(y(i)%t.sub.'000001')2.0(y(5)%t.sub.'000020')
2624	enddo
2625	else
2626	do i=1,c_%nd2-2*ndel ! should it be 1 to 4?
2627	s1%disp_p(i) = 0
2628	enddo
2629	endif
2630	! write(6,*) "dispersion derivative (1) = ",s1%disp_p(1)
2631	! checked the above yields non-zero value...
2632
2633	end subroutine computeDeltapDependency
2634	!____________________________________________________________________________________________
2635
2636	! --- a set of routines to track extremas of Twiss functions
2637	subroutine resetBetaExtremas()
2638	integer i,j
2639	do i=1,3
2640	do j=1,3
2641	resetBetaExtrema(i,j) = .true.
2642	enddo
2643	enddo
2644	end subroutine resetBetaExtremas
2645	subroutine trackBetaExtrema(i,j,value)
2646	implicit none
2647	integer :: i,j
2648	real(dp) :: value
2649	if (resetBetaExtrema(i,j)) then
2650	resetBetaExtrema(i,j) = .false.
2651	minBeta(i,j) = value
2652	maxBeta(i,j) = value
2653	! write(80,*) 'first time in trackBetaExtrema for ',i,j
2654	else
2655	if (minBeta(i,j) .gt. value) then
2656	minBeta(i,j) = value
2657	elseif (maxBeta(i,j) .lt. value) then
2658	maxBeta(i,j) = value
2659	endif
2660	endif
2661	end subroutine trackBetaExtrema
2662	! --- end of set of routines
2663
2664
2665	subroutine orbitRms(summary_table_name)
2666	implicit none
2667	character(48) :: summary_table_name
2668
2669	real(dp) :: state(6) ! 6 dimensional state space usually referred to as 'x'
2670	real(dp) :: x(6) ! the 6 dimensional state space
2671	real(kind(1d0)) :: xrms(6)
2672	real(kind(1d0)) :: xcomax, pxcomax, ycomax, pycomax
2673	integer :: i, j
2674
2675	real(kind(1d0)) :: get_value
2676
2677	if (get_value('ptc_twiss ','closed_orbit ').eq.0) then
2678	! for a line or if we don't mention the closed_orbit, xcorms makes no sense
2679	call double_to_table_curr(summary_table_name,'xcorms ',0d0)
2680	call double_to_table_curr(summary_table_name,'pxcorms ',0d0)
2681	call double_to_table_curr(summary_table_name,'ycorms ',0d0)
2682	call double_to_table_curr(summary_table_name,'pycorms ',0d0)
2683	call double_to_table_curr(summary_table_name,'xcomax ',0d0)
2684	call double_to_table_curr(summary_table_name,'pxcomax ',0d0)
2685	call double_to_table_curr(summary_table_name,'ycomax ',0d0)
2686	call double_to_table_curr(summary_table_name,'pycomax ',0d0)
2687	else
2688
2689
2690	call make_node_layout(my_ring) ! essential: the way to look inside the magnets
2691	state = zero
2692	call find_orbit(my_ring,state,1,default,c_1d_7) ! 1 for the first element
2693
2694	xcomax = state(1)
2695	pxcomax = state(2)
2696	ycomax = state(3)
2697	pycomax = state(4)
2698
2699	x=state
2700	xrms = zero
2701	do i=1,my_ring%n
2702	!call find_orbit(my_ring,state,i,default,1.d-5) ! i for the ith element?
2703	call track(my_ring,x,i,i+1,default) ! track x directly!
2704	! write(0,*) "xco(find_orbit)=",state(1),"xco(tracked)=",x(1)
2705	do j=1,6
2706	xrms(j) = xrms(j) + x(j)*x(j)
2707	enddo
2708	if (x(1)>xcomax) then
2709	xcomax = x(1)
2710	endif
2711	if (x(2)>pxcomax) then
2712	pxcomax = x(2)
2713	endif
2714	if (x(3)>ycomax) then
2715	ycomax = x(3)
2716	endif
2717	if (x(4)>pycomax) then
2718	pycomax = x(4)
2719	endif
2720	enddo
2721
2722	xrms = sqrt(xrms / my_ring%n)
2723
2724	call double_to_table_curr(summary_table_name,'xcorms ',xrms(1))
2725	call double_to_table_curr(summary_table_name,'pxcorms ',xrms(2))
2726	call double_to_table_curr(summary_table_name,'ycorms ',xrms(3))
2727	call double_to_table_curr(summary_table_name,'pycorms ',xrms(4))
2728	call double_to_table_curr(summary_table_name,'xcomax ',xcomax)
2729	call double_to_table_curr(summary_table_name,'pxcomax ',pxcomax)
2730	call double_to_table_curr(summary_table_name,'ycomax ',ycomax)
2731	call double_to_table_curr(summary_table_name,'pycomax ',pycomax)
2732
2733
2734	endif
2735
2736	end subroutine orbitRms
2737	!____________________________________________________________________________________________
2738
2739
2740	subroutine getBeamBeam()
2741	implicit none
2742	integer :: i,e,elcode
2743	integer, external :: restart_sequ, & ! restart beamline and return number of beamline node
2744	advance_node ! advance to the next node in expanded sequence
2745	! =0 (end of range), =1 (else)
2746	REAL(KIND(1d0)), external :: node_value !/returns value for parameter par of current element /
2747	type(fibre), pointer :: p
2748	real (dp) :: fk
2749
2750	TYPE(INTEGRATION_NODE),POINTER :: CURR_SLICE
2751
2752	e=restart_sequ()
2753	p=>my_ring%start
2754	do e=1, my_ring%n !in slices e goes to nelem + 1 because the last slice is the fist one.
2755
2756	elcode=node_value('mad8_type ')
2757
2758	if (elcode .eq. 22) then
2759
2760	if (getdebug() > 1 ) then
2761	write(6,*) " Beam-Beam position at element named >>",p%mag%name,"<<"
2762	endif
2763
2764	CURR_SLICE => p%t1
2765
2766	do while (.not. (CURR_SLICE%cas==case0.or.CURR_SLICE%cas==caset) )
2767	if (associated(CURR_SLICE,p%t2)) exit
2768	CURR_SLICE => CURR_SLICE%next
2769	!print*, CURR_SLICE%cas
2770	enddo
2771
2772	!print *, 'BB Node Case NO: ',CURR_SLICE%cas
2773
2774	if(((CURR_SLICE%cas==case0).or.(CURR_SLICE%cas==caset))) then !must be 0 or 3
2775
2776	if(.not.associated(CURR_SLICE%BB)) call alloc(CURR_SLICE%BB)
2777
2778	call getfk(fk)
2779	CURR_SLICE%bb%fk = fk
2780	CURR_SLICE%bb%sx = node_value('sigx ')
2781	CURR_SLICE%bb%sy = node_value('sigy ')
2782	CURR_SLICE%bb%xm = node_value('xma ')
2783	CURR_SLICE%bb%ym = node_value('yma ')
2784	CURR_SLICE%bb%PATCH=.true.
2785	if (getdebug() > 2 ) then
2786	print*, "BB fk=",CURR_SLICE%bb%fk
2787	print*, "BB sx=",CURR_SLICE%bb%sx
2788	print*, "BB sy=",CURR_SLICE%bb%sy
2789	print*, "BB xm=",CURR_SLICE%bb%xm
2790	print*, "BB ym=",CURR_SLICE%bb%ym
2791	endif
2792
2793	do_beam_beam = .true.
2794
2795	else
2796	call fort_warn('getBeamBeam: ','Bad node case for BeamBeam')
2797	endif
2798
2799	endif
2800
2801	i=advance_node() ! c-code go to the next node -> the passed value is never used, just to shut up a compiler
2802	p=>p%next
2803	enddo
2804
2805	end subroutine getBeamBeam
2806
2807
2808	end module madx_ptc_twiss_module

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

Context Navigation

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

Download in other formats: