Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

madx_ptc_module.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: 89.2 KB

Line
1	MODULE ptc_results
2	USE madx_keywords
3	implicit none
4	public
5	integer :: number_variables = 6
6	integer :: order = 20
7	character(len = 2), dimension(6) :: ptc_variables = (/'x ','xp','y ','yp','z ','dp'/)
8	character(len = 2) :: ptc_var
9	type(normalform) n
10	type (pbresonance) pbrg,pbrh
11	END MODULE ptc_results
12
13	MODULE madx_ptc_module
14	use S_fitting_new
15	! USE madx_keywords
16	USE madx_ptc_setcavs_module
17	USE madx_ptc_knobs_module
18	use madx_ptc_intstate_module, only : getdebug
19
20	implicit none
21	public
22	logical(lp) mytime
23	integer icav
24
25	integer :: universe=0,index_mad=0,EXCEPTION=0
26	integer ipause
27	integer,external :: mypause
28	real(kind(1d0)) get_value,node_value
29	type(layout),pointer :: MY_RING
30	type(mad_universe),target :: m_u
31	integer, private, parameter :: mynreso=20
32	integer, private, dimension(4) :: iia,icoast
33	real(dp) :: mux_default=c_0_28, muy_default=c_0_31, muz_default=c_1d_3
34	integer, private, allocatable :: J(:)
35	logical(lp) :: savemaps=.false.
36	logical(lp) :: resplit,even
37	real(dp) my_thin,my_xbend
38
39	type mapbuffer
40	type(universal_taylor) :: unimap(6)
41	real(dp) :: s
42	character(nlp+1) :: name
43	end type mapbuffer
44
45	type(mapbuffer), pointer :: maps(:) !buffered maps from the last twiss
46	integer :: mapsorder = 0 !order of the buffered maps, if 0 maps no maps buffered
47	integer :: mapsicase = 0
48
49	CONTAINS
50
51	subroutine ptc_create_universe()
52	implicit none
53	real(kind(1d0)) get_value
54
55	print77=.false.
56	read77 =.false.
57	lingyun_yang=get_value('ptc_create_universe ','ntpsa ').ne.0
58	lielib_print(6)=get_value('ptc_create_universe ','symprint ')
59
60	nullify(maps)
61
62	if (getdebug()==0) global_verbose = .false.
63	if (getdebug()>0) then
64	print*,"Now PTC"
65	endif
66	sector_nmul_max = get_value('ptc_create_universe ','sector_nmul_max ')
67
68	! print*,">>ss1<< old sector_nmul",sector_nmul
69
70	sector_nmul = get_value('ptc_create_universe ','sector_nmul ')
71
72	! print*,">>ss1<< new sector_nmul",sector_nmul
73
74	if(sector_nmul_max.lt.sector_nmul) then
75	call aafail('sector_nmul_max must be larger than sector_nmul: ',&
76	'check your ptc_create_universe input')
77	endif
78	call set_up_universe(m_u)
79	universe=universe+1
80
81
82	end subroutine ptc_create_universe
83	!_________________________________________________________________
84
85	subroutine ptc_create_layout()
86	implicit none
87	real(kind(1d0)) get_value
88
89	if(universe.le.0.or.EXCEPTION.ne.0) then
90	call fort_warn('return from ptc_create_layout: ',' no universe created')
91	return
92	endif
93
94	call append_empty_layout(m_u)
95	index_mad=index_mad+1
96	my_ring=>m_u%end
97
98	call ptc_input()
99
100	if(EXCEPTION.eq.1) then
101	call fort_warn('wrong magnet type KINDI which must be: ','1, 2, 3')
102	return
103	endif
104
105	cavsareset = .false.
106	mytime=get_value('ptc_create_layout ','time ').ne.0
107
108	if(mytime) then
109	default=getintstate()
110	default=default+time
111	call setintstate(default)
112	endif
113
114	end subroutine ptc_create_layout
115	!_________________________________________________________________
116
117	subroutine ptc_move_to_layout()
118	implicit none
119	real(kind(1d0)) get_value
120	integer my_index
121
122	if(universe.le.0.or.EXCEPTION.ne.0) then
123	call fort_warn('return from ptc_move_to_layout: ',' no universe created')
124	return
125	endif
126
127	my_index = get_value('ptc_move_to_layout ','index ')
128
129	if(my_index.gt.index_mad.or.my_index.le.0) then
130	call fort_warn('return from ptc_move_to_layout: ',' layout outside allowed range')
131	print*," Allowed range 0 < ",index_mad
132	return
133	endif
134
135	call move_to_layout_i(m_u,my_ring,my_index)
136
137	end subroutine ptc_move_to_layout
138	!_________________________________________________________________
139
140	subroutine ptc_input()
141	use twtrrfi
142	use twiss0fi
143	use name_lenfi
144	implicit none
145	logical(lp) particle,doneit,isclosedlayout
146	integer i,j,k,code,nt,icount,nn,ns,nd,mg,get_string
147	! integer get_option
148	integer double_from_table_row,table_column_exists
149	integer restart_sequ,advance_node,n_ferr,node_fd_errors
150	integer, parameter :: nt0=20000
151	real(dp) l,l_machine,energy,kin,brho,beta0,p0c,pma,e0f,lrad,charge
152	real(dp) f_errors(0:maxferr),aperture(maxnaper),normal(0:maxmul)
153	real(dp) patch_ang(3),patch_trans(3)
154	real(dp) skew(0:maxmul),field(2,0:maxmul),fieldk(2),myfield(2*maxmul+2)
155	real(dp) gamma,gamma2,gammatr2,freq,offset_deltap
156	real(dp) fint,fintx,div,muonfactor,edge,rhoi,hgap,corr,tanedg,secedg,psip
157	real(dp) sk1,sk1s,sk2,sk2s,sk3,sk3s,tilt,dum1,dum2
158	REAL(dp) :: normal_0123(0:3), skew_0123(0:3) ! <= knl(1), ksl(1)
159	real(dp) gammatr,ks,ksi
160	real(kind(1d0)) get_value,node_value
161	character(name_len) name
162	character(name_len) aptype
163	type(keywords) key
164	character(20) keymod0,keymod1
165	character(name_len) magnet_name
166	logical(lp) exact0,no_cavity_totalpath
167	integer exact1
168	integer sector_nmul_max0,sector_nmul0
169	integer model
170	integer method,method0,method1
171	integer nst0,nst1,ord_max,kk
172	REAL (dp) :: tempdp,bvk
173	logical(lp):: ptcrbend,truerbend,errors_out
174	! Etienne helical
175	character(nlp) heli(100)
176	integer mheli,helit,ihelit
177	type(fibre), pointer :: p
178	!---------------------------------------------------------------
179	!---------------------------------------------------------------
180	if (getdebug() > 1) then
181	print *, '--------------------------------------------------------------'
182	print *, '--------------------------------------------------------------'
183	print *, '------ E X E C U T I N G P T C I N P U T --------'
184	print *, '--------------------------------------------------------------'
185	print *, '--------------------------------------------------------------'
186	endif
187
188	energy=get_value('probe ','energy ')
189	pma=get_value('probe ','mass ')
190	charge=get_value('probe ','charge ')
191	bvk=get_value('probe ','bv ')
192
193	e0f=sqrt(ENERGY2-pma2)
194
195	if (getdebug() > 0) then
196	print *, 'MAD-X Beam Parameters'
197	print '(a26, e13.6)', ' Energy :',energy
198	print '(a26, e13.6)', ' Kinetic Energy :',energy-pma
199	print '(a26, e13.6)', ' Particle Rest Mass :',pma
200	print '(a26, e13.6)', ' Momentum :',e0f
201	endif
202
203
204
205	beta0=e0f/ENERGY
206
207
208	if(abs(pma-pmae)/pmae<c_0_002) then
209	if (getdebug() > 1) then
210	print *,'Executing MAKE_STATES(TRUE), i.e. ELECTRON beam'
211	endif
212	particle=.true.
213	CALL MAKE_STATES(PARTICLE)
214	elseif(abs(pma-pmap)/pmap<c_0_002) then
215	if (getdebug() > 1) then
216	print *,'Executing MAKE_STATES(FALSE), i.e. PROTON beam'
217	endif
218	particle=.false.
219	CALL MAKE_STATES(PARTICLE)
220	else
221	if (getdebug() > 1) then
222	print '(a, f8.4, a)','Executing MAKE_STATES(',pma/pmae,'), i.e. PROTON beam'
223	endif
224	muonfactor=pma/pmae
225	CALL MAKE_STATES(muonfactor)
226	endif
227
228	!the state is cleared at this stage
229	call setintstate(default)
230	!valid October 2002: oldscheme=.false.
231	!!valid October 2002: oldscheme=.true.
232
233	if (getdebug()==0) global_verbose = .false.
234
235	! with_external_frame=.false.
236	! with_internal_frame=.false.
237	! with_chart=.false.
238	! with_patch=.false.
239
240	! Global Keywords
241
242	if (getdebug() > 1) then
243	print *, '=============================================================='
244	print *, 'INPUT PARAMETERS ARE:'
245	endif
246
247	sector_nmul_max0 = sector_nmul_max
248	if (getdebug() > 1) then
249	print*,' Global max sector_nmul: ',sector_nmul_max0
250	endif
251
252	sector_nmul0 = sector_nmul
253	if (getdebug() > 1) then
254	print*,' Global sector_nmul: ',sector_nmul0
255	endif
256
257
258	model = get_value('ptc_create_layout ','model ')
259	if (getdebug() > 1) then
260	print*,' Global Model code is : ',model
261	endif
262
263	!*****************************
264	! MODEL Settings
265	!*****************************
266	select case(model)
267	CASE(1)
268	keymod0 = "DRIFT_KICK "
269	CASE(2)
270	keymod0 = "MATRIX_KICK "
271	CASE(3)
272	keymod0 = "DELTA_MATRIX_KICK"
273	CASE DEFAULT
274	PRINT *, 'EXCEPTION occured: Can not recognize model type ',model
275	EXCEPTION=1
276	index_mad=-1
277	ipause=mypause(444)
278	RETURN
279	END SELECT
280
281	if (getdebug() > 1) then
282	print*,' Global Model name (keymod0) is : ',keymod0
283	endif
284
285	method = get_value('ptc_create_layout ','method ')
286	if (getdebug() > 1) then
287	print*,' Global method is: ',method
288	endif
289
290	!*****************************
291	! METHOD Settings
292	!*****************************
293	select case(method)
294	CASE(2)
295	method0 = method
296	CASE(4)
297	method0 = method
298	CASE(6)
299	method0 = method
300	CASE DEFAULT
301	PRINT *, 'EXCEPTION occured: Can not recognize method order ',method
302	EXCEPTION=1
303	index_mad=-1
304	ipause=mypause(444)
305	RETURN
306	END SELECT
307
308	exact0 = get_value('ptc_create_layout ','exact ') .ne. 0
309	if (getdebug() > 1) then
310	print*,' Global exact is: ',exact0
311	endif
312
313	nst0 = get_value('ptc_create_layout ','nst ')
314	if (getdebug() > 1) then
315	print*,' Global Number of Integration Steps (nst) is: ',nst0
316	endif
317
318	! MAD-X specials
319	! madlength = get_option('rbarc ') .eq. 0
320	madlength = .false.
321	if (getdebug() > 1) then
322	print*,' global rbend_length: ',madlength
323	endif
324
325	mad = get_value('ptc_create_layout ','mad_mult ') .ne. 0
326	if (getdebug() > 1) then
327	print*,' global mad_mult as in mad8: ',mad
328	endif
329
330	mad8 = get_value('ptc_create_layout ','mad8 ') .ne. 0
331	if (getdebug() > 1) then
332	print*,' rbend as in mad8 (only global): ',mad8
333	endif
334
335	gamma = get_value('probe ','gamma ')
336	if (getdebug() > 1) then
337	print*,' gamma: ',gamma
338	endif
339
340	if (table_column_exists('summ ', 'gammatr ', 1, gammatr).ne.0) then
341	k = double_from_table_row('summ ','gammatr ',1,gammatr)
342	else
343	gammatr = 0
344	endif
345
346	if (getdebug() > 1) then
347	print*,' gammatr: ',gammatr
348	endif
349
350	gamma2 = gamma**2
351	gammatr2 = gammatr**2
352
353	if (getdebug() > 1) then
354	print *, '=============================================================='
355	print *, ''
356	endif
357
358	! call Set_Up(MY_RING)
359
360	if (getdebug() > 0) then
361	print *, 'Setting MADx with '
362	print *, ' energy ',energy
363	print *, ' method ',method0
364	print *, ' Num. of steps ',nst0
365	print *, ' charge ',charge
366	endif
367	! etienne helical
368	helit=0
369	call kanalnummer(mheli)
370	open(unit=mheli,file='helical.txt',status='OLD',err=1001)
371	read(mheli,*) helit
372	if(helit>100) then
373	write(6,*) " too many helical dipole ",helit
374	stop 99
375	endif
376	do ihelit=1,helit
377	read(mheli,*) heli(ihelit)
378	CALL CONTEXT(heli(ihelit))
379	enddo
380	close(mheli)
381	1001 continue
382	helit=0
383	call kanalnummer(mheli)
384	open(unit=mheli,file='sixtrack_compatible.txt',status='OLD',err=1002)
385	read(mheli,*) sixtrack_compatible
386	close(mheli)
387	1002 continue
388	! end of etienne helical
389
390	! preliminary setting
391	! my_ring%charge=1
392	initial_charge=1
393	CALL SET_MADx(energy=energy,METHOD=method0,STEP=nst0)
394	if (getdebug() > 1) then
395	print *, 'MADx is set'
396	endif
397
398	icav=0
399	nt=0
400	j=restart_sequ()
401	j=0
402	l_machine=zero
403
404	errors_out = get_value('ptc_create_layout ','errors_out ').ne.0
405	magnet_name=" "
406	if(errors_out) mg = get_string('ptc_create_layout ','magnet_name ',magnet_name)
407
408	10 continue
409	nst1=node_value("nst ")
410	if(nst1.gt.0) then
411	nstd = nst1
412	else
413	nstd = nst0
414	endif
415
416	call zero_key(key)
417
418	j=j+1
419	nt=nt+1
420	if(nt==nt0) then
421	call fort_warn("Potential problem for very large structure: ","More than 20'000 elements found")
422	endif
423	icount=0
424	l=zero
425	l=node_value('l ')
426	key%list%l=l
427	l_machine=l_machine+l
428	code=node_value('mad8_type ')
429	if(code.eq.39) code=15
430	if(code.eq.38) code=24
431	call element_name(name,name_len)
432	key%list%name=name
433
434	call node_name(name,name_len)
435	key%list%vorname=name
436
437	!frs&piotr 18 Dec 2007: sector_nmul must stay global for the time being
438	!local, if present, superseed global at current node
439
440
441	!*****************************
442	! MODEL Settings
443	!*****************************
444
445	model = node_value('model ')
446	keymod1 = " "
447	select case(model)
448	CASE(1)
449	keymod1 = "DRIFT_KICK "
450	CASE(2)
451	keymod1 = "MATRIX_KICK "
452	CASE(3)
453	keymod1 = "DELTA_MATRIX_KICK"
454	END SELECT
455
456
457	if(keymod1.ne." ") then
458	key%model=keymod1
459	else
460	key%model=keymod0
461	endif
462	method1=node_value("method ")
463	if(method1.eq.2.or.method1.eq.4.or.method1.eq.6) then
464	metd = method1
465	else
466	metd = method0
467	endif
468
469	exact1=node_value("exact ")
470
471	if(exact1.eq.0.or.exact1.eq.1) then
472	EXACT_MODEL = exact1 .ne. 0
473	else
474	EXACT_MODEL = exact0
475	endif
476
477	!special node keys
478	key%list%permfringe=node_value("permfringe ") .ne. zero
479	key%list%kill_ent_fringe=node_value("kill_ent_fringe ") .ne. zero
480	key%list%kill_exi_fringe=node_value("kill_exi_fringe ") .ne. zero
481	key%list%bend_fringe=node_value("bend_fringe ") .ne. zero
482
483	nn=name_len
484	call node_string('apertype ',aptype,nn)
485	call dzero(aperture,maxnaper)
486	call get_node_vector('aperture ',nn,aperture)
487	if(.not.((aptype.eq."circle".and.aperture(1).eq.zero).or.aptype.eq." ")) then
488	c_%APERTURE_FLAG=.true.
489	select case(aptype)
490	case("circle")
491	key%list%aperture_on=.true.
492	key%list%aperture_kind=1
493	key%list%aperture_r(1)=aperture(1)
494	key%list%aperture_r(2)=aperture(1)
495	case("ellipse")
496	key%list%aperture_on=.true.
497	key%list%aperture_kind=1
498	key%list%aperture_r(1)=aperture(1)
499	key%list%aperture_r(2)=aperture(2)
500	case("rectangle")
501	key%list%aperture_on=.true.
502	key%list%aperture_kind=2
503	key%list%aperture_x=aperture(1)
504	key%list%aperture_y=aperture(2)
505	! case("lhcscreen")
506	case("rectellipse")
507	key%list%aperture_on=.true.
508	key%list%aperture_kind=3
509	key%list%aperture_x=aperture(1)
510	key%list%aperture_y=aperture(2)
511	key%list%aperture_r(1)=aperture(3)
512	key%list%aperture_r(2)=aperture(4)
513	case("marguerite")
514	key%list%aperture_on=.true.
515	key%list%aperture_kind=4
516	key%list%aperture_r(1)=aperture(1)
517	key%list%aperture_r(2)=aperture(2)
518	case("racetrack")
519	key%list%aperture_on=.true.
520	key%list%aperture_kind=5
521	key%list%aperture_x=aperture(1)
522	key%list%aperture_y=aperture(2)
523	key%list%aperture_r(1)=aperture(3)
524	case("general")
525	key%list%aperture_kind=6
526	print*,"General aperture not implemented"
527	stop
528	end select
529	endif
530	call append_empty(my_ring)
531
532	!print *,'Element code is ',code
533
534	select case(code)
535	case(0,4,25)
536	key%magnet="marker"
537	case(22)
538	call fort_warn('ptc_input: ','Element Beam-Beam, must use slice tracking to get effect')
539	key%magnet="marker"
540	case(1,11,20,21)
541	key%magnet="drift"
542	CALL CONTEXT(key%list%name)
543
544	do ihelit=1,helit
545	IF(index(key%list%name,heli(ihelit)(1:len_trim(heli(ihelit))))/=0) then
546	key%magnet="helicaldipole"
547	write(6,*) " drift ",key%list%name, " became helical dipole in PTC "
548	endif
549	enddo
550	! end etienne Helical
551	case(2) ! PTC accepts mults
552	if(l.eq.zero) then
553	key%magnet="marker"
554	goto 100
555	endif
556	key%magnet="rbend"
557	!VK
558	CALL SUMM_MULTIPOLES_AND_ERRORS (l, key, normal_0123,skew_0123,ord_max)
559
560	tempdp=sqrt(normal_0123(0)normal_0123(0)+skew_0123(0)skew_0123(0))
561	key%list%b0=bvk(node_value('angle ')+tempdpl)
562
563	! print*, "RBEND: Angle: ", node_value('angle ')," tempdp ", tempdp, " l ", l
564	! print*, "RBEND: normal: ",normal_0123(0)," skew: ",skew_0123(0)
565
566	key%list%k(2)=node_value('k1 ')+ key%list%k(2)
567	key%list%k(3)=node_value('k2 ')+ key%list%k(3)
568	key%list%k(4)=node_value('k3 ')+ key%list%k(4)
569
570	key%list%ks(2)=node_value('k1s ')+ key%list%ks(2)
571	key%list%ks(3)=node_value('k2s ')+ key%list%ks(3)
572	key%list%ks(4)=node_value('k3s ')+ key%list%ks(4)
573
574	if(EXACT_MODEL.and.(node_value('angle ').eq.zero)) then
575	key%magnet="quadrupole"
576	key%tiltd=node_value('tilt ')
577	else
578
579	! Gymnastic needed since PTC expects MAD8 convention
580	key%list%t1=node_value('e1 ')
581	key%list%t2=node_value('e2 ')
582	key%list%hgap=node_value('hgap ')
583	! key%list%fint=node_value('fint ')
584	fint=node_value('fint ')
585	fintx=node_value('fintx ')
586	if((fintx.ne.fint).and.(fintx.gt.zero.and.fint.gt.zero)) then
587	print*," The fint and fintx must be the same at each end or each might be zero"
588	stop
589	endif
590	if(fint.gt.zero) then
591	key%list%fint=fint
592	if(fintx.eq.zero) key%list%kill_exi_fringe=my_true
593	else
594	if(fintx.gt.zero) then
595	key%list%fint=fintx
596	key%list%kill_ent_fringe=my_true
597	else
598	key%list%fint=zero
599	endif
600	endif
601	key%list%h1=node_value('h1 ')
602	key%list%h2=node_value('h2 ')
603	key%tiltd=node_value('tilt ')
604	if(tempdp.gt.0) key%tiltd=key%tiltd + atan2(skew_0123(0),normal_0123(0))
605	ptcrbend=node_value('ptcrbend ').ne.0
606	if(ptcrbend) then
607	call context(key%list%name)
608	truerbend=node_value('truerbend ').ne.0
609	if(truerbend) then
610	key%magnet="TRUERBEND"
611	if(key%list%t2/=zero) then
612	write(6,*) " The true parallel face bend "
613	write(6,*) " only accepts the total angle and e1 as an input "
614	write(6,*) " if e1=0, then the pipe angle to the entrance face is "
615	write(6,*) " angle/2. It is a normal rbend."
616	write(6,*) " If e1/=0, then the pipe angle to the entrance face is "
617	write(6,*) ' angle/2+e1 and the exit pipe makes an angle "angle/2-e1" '
618	write(6,*) " with the exit face."
619	write(6,*) " The offending non-zero t2 = (e2 - angle/2) is set to zero! "
620	write(6,*) " Make sure that this is what you want!!! "
621	! write(6,*) " CHANGE YOUR LATTICE FILRE."
622	! stop 666
623	key%list%t2=zero
624	endif
625	else
626	key%magnet="WEDGRBEND"
627	endif
628	endif
629	endif
630	if(errors_out) then
631	if(key%list%name(:len_trim(magnet_name)-1).eq. &
632	magnet_name(:len_trim(magnet_name)-1)) then
633	call string_to_table_curr('errors_dipole ', 'name ', key%list%name)
634	call double_to_table_curr('errors_dipole ', 'k0l ', bvk*key%list%b0)
635	call augment_count('errors_dipole ')
636	endif
637	endif
638	case(3) ! PTC accepts mults watch out sector_nmul defaulted to 4
639	if(l.eq.zero) then
640	key%magnet="marker"
641	goto 100
642	endif
643	key%magnet="sbend"
644	!VK
645	CALL SUMM_MULTIPOLES_AND_ERRORS (l, key, normal_0123,skew_0123,ord_max)
646	if(sector_nmul_max.lt.ord_max.and.EXACT_MODEL) call aafail('the order of multipoles in a sbend in exact mode cannot be ',&
647	&'larger than sector_mul_max: check your ptc_create_universe input')
648
649	tempdp=sqrt(normal_0123(0)normal_0123(0)+skew_0123(0)skew_0123(0))
650	key%list%b0=bvk(node_value('angle ')+ tempdpl)
651
652	key%list%k(2)=node_value('k1 ')+ key%list%k(2)
653	key%list%k(3)=node_value('k2 ')+ key%list%k(3)
654	key%list%k(4)=node_value('k3 ')+ key%list%k(4)
655
656	key%list%ks(2)=node_value('k1s ')+ key%list%ks(2)
657	key%list%ks(3)=node_value('k2s ')+ key%list%ks(3)
658	key%list%ks(4)=node_value('k3s ')+ key%list%ks(4)
659
660	key%list%t1=node_value('e1 ')
661	key%list%t2=node_value('e2 ')
662	key%list%hgap=node_value('hgap ')
663	! key%list%fint=node_value('fint ')
664	fint=node_value('fint ')
665	fintx=node_value('fintx ')
666	if((fintx.ne.fint).and.(fintx.gt.zero.and.fint.gt.zero)) then
667	print*," The fint and fintx must be the same at each end or each might be zero"
668	stop
669	endif
670	if(fint.gt.zero) then
671	key%list%fint=fint
672	if(fintx.eq.zero) key%list%kill_exi_fringe=my_true
673	else
674	if(fintx.gt.zero) then
675	key%list%fint=fintx
676	key%list%kill_ent_fringe=my_true
677	else
678	key%list%fint=zero
679	endif
680	endif
681	key%list%h1=node_value('h1 ')
682	key%list%h2=node_value('h2 ')
683	key%tiltd=node_value('tilt ')
684	if(tempdp.gt.0) key%tiltd=key%tiltd + atan2(skew_0123(0),normal_0123(0))
685	if(errors_out) then
686	if(key%list%name(:len_trim(magnet_name)-1).eq. &
687	magnet_name(:len_trim(magnet_name)-1)) then
688	call string_to_table_curr('errors_dipole ', 'name ', key%list%name)
689	call double_to_table_curr('errors_dipole ', 'k0l ', bvk*key%list%b0)
690	call augment_count('errors_dipole ')
691	endif
692	endif
693	case(5)
694	key%magnet="quadrupole"
695	!VK
696	CALL SUMM_MULTIPOLES_AND_ERRORS (l, key, normal_0123,skew_0123,ord_max)
697
698	! Read data & fill %k(:), %ks(:) arrays which are
699	! summs of multipoles and errors
700
701	! quadrupole components
702	sk1= node_value('k1 ')
703	sk1s=node_value('k1s ')
704	tilt=node_value('tilt ')
705	dum1=key%list%k(2)-normal_0123(1)
706	dum2=key%list%ks(2)-skew_0123(1)
707
708	if(dum1.ne.zero.or.dum2.ne.zero) then !
709	sk1= sk1 +dum1 !
710	sk1s=sk1s+dum2 !
711	endif !
712	if (sk1s .ne. zero) then !
713	tilt = -atan2(sk1s, sk1)/two + tilt !
714	sk1 = sqrt(sk12 + sk1s2) !
715	endif !
716	key%list%k(2) =sk1 !
717	key%list%ks(2)=zero ! added by VK !
718	key%tiltd=tilt !==========================================!
719
720	!================================================================
721	! dipole component not active in MAD-X proper
722	key%list%k(1)=key%list%k(1)+bvk*node_value('k0 ')
723
724	case(6)
725	key%magnet="sextupole"
726	!VK
727	CALL SUMM_MULTIPOLES_AND_ERRORS (l, key, normal_0123,skew_0123,ord_max)
728
729	! sextupole components
730	sk2= node_value('k2 ')
731	sk2s=node_value('k2s ')
732	tilt=node_value('tilt ')
733	dum1=key%list%k(3)-normal_0123(2)
734	dum2=key%list%ks(3)-skew_0123(2)
735
736	if(dum1.ne.zero.or.dum2.ne.zero) then !
737	sk2= sk2 +dum1 !
738	sk2s=sk2s+dum2 !
739	endif !
740	if (sk2s .ne. zero) then !
741	tilt = -atan2(sk2s, sk2)/three + tilt !
742	sk2 = sqrt(sk22 + sk2s2) !
743	endif !
744	key%list%k(3) =sk2 !
745	key%list%ks(3)=zero ! added by VK !
746	key%tiltd=tilt !==========================================!
747
748	!================================================================
749	if(errors_out) then
750	if(key%list%name(:len_trim(magnet_name)-1).eq. &
751	magnet_name(:len_trim(magnet_name)-1)) then
752	call string_to_table_curr('errors_total ', 'name ', key%list%name)
753	myfield(:) = zero
754	do kk=1,maxmul
755	myfield(2*kk-1) = key%list%k(kk)
756	myfield(2*kk) = key%list%ks(kk)
757	enddo
758	call vector_to_table_curr('errors_total ', 'k0l ', myfield(1), i)
759	call augment_count('errors_total ')
760	endif
761	endif
762
763	case(7)
764	key%magnet="octupole"
765	!VK
766	CALL SUMM_MULTIPOLES_AND_ERRORS (l, key, normal_0123,skew_0123,ord_max)
767
768	! octupole components
769	sk3= node_value('k3 ')
770	sk3s=node_value('k3s ')
771
772	tilt=node_value('tilt ')
773	dum1=key%list%k(4)-normal_0123(3)
774	dum2=key%list%ks(4)-skew_0123(3)
775
776	if(dum1.ne.zero.or.dum2.ne.zero) then !
777	sk3= sk3 +dum1 !
778	sk3s=sk3s+dum2 !
779	endif !
780	if (sk3s .ne. zero) then !
781	tilt = -atan2(sk3s, sk3)/four + tilt !
782	sk3 = sqrt(sk32 + sk3s2) !
783	endif !
784	key%list%k(4) =sk3 !
785	key%list%ks(4)=zero ! added by VK !
786
787	key%tiltd=tilt !==========================================!
788
789	!================================================================
790
791	case(8)
792	key%magnet="multipole"
793	!---- Multipole components.
794	call dzero(f_errors,maxferr+1)
795	n_ferr = node_fd_errors(f_errors)
796	call dzero(normal,maxmul+1)
797	call dzero(skew,maxmul+1)
798	call get_node_vector('knl ',nn,normal)
799	call get_node_vector('ksl ',ns,skew)
800	if(nn.ge.NMAX) nn=NMAX-1
801	if(ns.ge.NMAX) ns=NMAX-1
802	do i=1,NMAX
803	key%list%k(i)=zero
804	key%list%ks(i)=zero
805	enddo
806	skew(0)=-skew(0) ! frs error found 30.08.2008
807	key%list%thin_h_angle=bvk*normal(0)
808	key%list%thin_v_angle=bvk*skew(0)
809	lrad=node_value('lrad ')
810	if(lrad.gt.zero) then
811	key%list%thin_h_foc=normal(0)*normal(0)/lrad
812	key%list%thin_v_foc=skew(0)*skew(0)/lrad
813	endif
814	if(nn.gt.0) then
815	do i=1,nn
816	key%list%k(i+1)=normal(i)
817	enddo
818	endif
819	if(ns.gt.0) then
820	do i=1,ns
821	key%list%ks(i+1)=skew(i)
822	enddo
823	endif
824	call dzero(field,2*(maxmul+1))
825	if (n_ferr .gt. 0) then
826	call dcopy(f_errors,field,n_ferr)
827	endif
828	nd = max(nn, ns, n_ferr/2)
829	if(nd.ge.maxmul) nd=maxmul-1
830	if(n_ferr.gt.0) then
831	do i=0,nd
832	key%list%k(i+1)=key%list%k(i+1)+field(1,i)
833	key%list%ks(i+1)=key%list%ks(i+1)+field(2,i)
834	enddo
835	endif
836	key%tiltd=node_value('tilt ')
837	if(errors_out) then
838	if(key%list%name(:len_trim(magnet_name)-1).eq. &
839	magnet_name(:len_trim(magnet_name)-1)) then
840	call string_to_table_curr('errors_field ', 'name ', key%list%name)
841	call string_to_table_curr('errors_total ', 'name ', key%list%name)
842	i=2*maxmul+2
843	myfield(:) = zero
844	do kk=1,nd+1
845	myfield(2*kk-1) = field(1,kk-1)
846	myfield(2*kk) = field(2,kk-1)
847	enddo
848	call vector_to_table_curr('errors_field ', 'k0l ', myfield(1), i)
849	myfield(:) = zero
850	do kk=1,nd+1
851	myfield(2*kk-1) = key%list%k(kk)
852	myfield(2*kk) = key%list%ks(kk)
853	enddo
854	call vector_to_table_curr('errors_total ', 'k0l ', myfield(1), i)
855	call augment_count('errors_field ')
856	call augment_count('errors_total ')
857	endif
858	endif
859	case(9) ! PTC accepts mults
860	key%magnet="solenoid"
861	ks=node_value('ks ')
862	if(l.ne.zero) then
863	key%list%bsol=bvk*ks
864	else
865	ksi=node_value('ksi ')
866	lrad=node_value('lrad ')
867	if(lrad.eq.zero.and.ks.ne.zero) lrad=ksi/ks
868	if(ksi.eq.zero.or.lrad.eq.zero) then
869	key%magnet="marker"
870	print*,"Thin solenoid: ",name," has no strength - set to marker"
871	else
872	key%list%bsol=bvk*ksi/lrad
873	key%list%ls=lrad
874	endif
875	endif
876	!VK
877	CALL SUMM_MULTIPOLES_AND_ERRORS (l, key, normal_0123,skew_0123,ord_max)
878
879	case(10)
880	key%magnet="rfcavity"
881	key%list%volt=bvk*node_value('volt ')
882	freq=c_1d6*node_value('freq ')
883	key%list%lag=node_value('lag ')*twopi
884	offset_deltap=get_value('ptc_create_layout ','offset_deltap ')
885	if(offset_deltap.ne.zero) then
886	default = getintstate()
887	default=default+totalpath0
888	call setintstate(default)
889	freq=freq((gammatr2-gamma2)offset_deltap/gammatr2/gamma2+one)
890	endif
891	key%list%freq0=freq
892	key%list%n_bessel=node_value('n_bessel ')
893	key%list%harmon=one
894	if(key%list%volt.ne.zero.and.key%list%freq0.ne.zero) icav=1
895	! case(11)
896	! key%magnet="elseparator"
897	! key%list%volt=node_value('ex ')
898	! key%list%lag=atan2(node_value('ey '),node_value('ex '))
899	! key%tiltd=node_value('tilt ')
900	m_u%end%HARMONIC_NUMBER=node_value('harmon ') ! etienne_harmon
901	no_cavity_totalpath=node_value('no_cavity_totalpath ').ne.0
902	if(no_cavity_totalpath) then
903	key%list%cavity_totalpath=0
904	else
905	key%list%cavity_totalpath=1
906	endif
907	case(12)
908	! actually our SROT element
909	key%magnet="CHANGEREF"
910	call dzero(patch_ang,3)
911	call dzero(patch_trans,3)
912	patch_ang(3)=node_value('angle ')
913	key%list%patchg=2
914	do i=1,3
915	key%list%ang(i)=patch_ang(i)
916	key%list%t(i)=patch_trans(i)
917	enddo
918	case(13)
919	! actually our YROT element
920	key%magnet="CHANGEREF"
921	call dzero(patch_ang,3)
922	call dzero(patch_trans,3)
923	patch_ang(2)=-node_value('angle ')
924	key%list%patchg=2
925	do i=1,3
926	key%list%ang(i)=patch_ang(i)
927	key%list%t(i)=patch_trans(i)
928	enddo
929	case(14,15,16) ! PTC accepts mults
930	call dzero(f_errors,maxferr+1)
931	n_ferr = node_fd_errors(f_errors)
932	do i=1,NMAX
933	key%list%k(i)=zero
934	key%list%ks(i)=zero
935	enddo
936	do i = 1, 2
937	fieldk(i) = zero
938	enddo
939	if (n_ferr .gt. 0) call dcopy(f_errors, fieldk, min(2, n_ferr))
940	if (l .eq. zero) then
941	div = one
942	else
943	div = l
944	endif
945	if(code.eq.14) then
946	key%magnet="hkicker"
947	key%list%k(1)=(node_value('kick ')+node_value('chkick ')+fieldk(1)/div)
948	else if(code.eq.15) then
949	key%magnet="kicker"
950	key%list%k(1)=(node_value('hkick ')+node_value('chkick ')+fieldk(1)/div)
951	key%list%ks(1)=(node_value('vkick ')+node_value('cvkick ')+fieldk(2)/div)
952	else if(code.eq.16) then
953	key%magnet="vkicker"
954	key%list%ks(1)=(node_value('kick ')+node_value('cvkick ')+fieldk(2)/div)
955	else
956	key%magnet="marker"
957	endif
958	i=2*maxmul+2
959	if(errors_out) then
960	if(key%list%name(:len_trim(magnet_name)-1).eq. &
961	magnet_name(:len_trim(magnet_name)-1)) then
962	myfield(:) = zero
963	myfield(1)=-key%list%k(1)
964	myfield(2)=key%list%ks(1)
965	call string_to_table_curr('errors_total ', 'name ', key%list%name)
966	call vector_to_table_curr('errors_total ', 'k0l ', myfield(1), i)
967	call augment_count('errors_total ')
968	endif
969	endif
970	key%tiltd=node_value('tilt ')
971	case(17)
972	key%magnet="hmonitor"
973	case(18)
974	key%magnet="monitor"
975	case(19)
976	key%magnet="vmonitor"
977	! case(20)
978	! key%magnet="ecollimator"
979	! key%list%x_col=node_value('xsize ')
980	! key%list%y_col=node_value('ysize ')
981	! key%tiltd=node_value('tilt ')
982	! case(21)
983	! key%magnet="rcollimator"
984	! key%list%x_col=node_value('xsize ')
985	! key%list%y_col=node_value('ysize ')
986	! key%tiltd=node_value('tilt ')
987	case(33)
988	!---- This is the dipedge element
989	edge= node_value('e1 ')
990	hgap= node_value('hgap ')
991	rhoi= bvk * node_value('h ')
992	fint= node_value('fint ')
993	corr= 2 * rhoi * hgap * fint
994	if(rhoi .ne. zero .and. ( edge .ne. zero .or. corr .ne. zero )) then
995	key%magnet="multipole"
996	tanedg = tan(edge)
997	secedg = one / cos(edge)
998	psip = edge - corr * secedg * (one + sin(edge)**2)
999	key%list%hf= rhoi * tanedg
1000	key%list%vf= -rhoi * tan(psip)
1001	else
1002	key%magnet="marker"
1003	endif
1004	case(24)
1005	key%magnet="instrument"
1006	key%tiltd=node_value('tilt ')
1007	case(27)
1008	key%magnet="twcavity"
1009	key%list%volt=bvk*node_value('volt ')
1010	freq=c_1d6*node_value('freq ')
1011	key%list%lag=node_value('lag ')*twopi
1012	offset_deltap=get_value('ptc_create_layout ','offset_deltap ')
1013	default=default+totalpath0 !fringe field calculation vitally relies on it!!!!
1014	if(offset_deltap.ne.zero) then
1015	freq=freq((gammatr2-gamma2)offset_deltap/gammatr2/gamma2+one)
1016	endif
1017	key%list%freq0=freq
1018	key%list%dphas=node_value("delta_lag ")
1019	key%list%psi=node_value("psi ")
1020	key%list%harmon=one
1021	if(key%list%volt.ne.zero.and.key%list%freq0.ne.zero) icav=1
1022	case(35)
1023	key%magnet="CHANGEREF"
1024	call dzero(patch_ang,3)
1025	call dzero(patch_trans,3)
1026	call get_node_vector('patch_ang ',3,patch_ang)
1027	call get_node_vector('patch_trans ',3,patch_trans)
1028	key%list%patchg=2
1029	do i=1,3
1030	key%list%ang(i)=patch_ang(i)
1031	key%list%t(i)=patch_trans(i)
1032	enddo
1033	case(37)
1034	key%magnet="rfcavity"
1035	key%list%volt=zero
1036	do i=1,NMAX
1037	key%list%k(i)=zero
1038	key%list%ks(i)=zero
1039	enddo
1040	key%list%k(1)=node_value('volt ')*c_1d_3
1041	! vertical crab
1042	! maybe requires a flip of sign
1043	! key%list%ks(1)= (+/-) node_value('volt ')*c_1d_3
1044	!
1045	freq=c_1d6*node_value('freq ')
1046	key%list%lag=node_value('lag ')*twopi+pih
1047	offset_deltap=get_value('ptc_create_layout ','offset_deltap ')
1048	if(offset_deltap.ne.zero) then
1049	default = getintstate()
1050	default=default+totalpath0
1051	call setintstate(default)
1052	freq=freq((gammatr2-gamma2)offset_deltap/gammatr2/gamma2+one)
1053	endif
1054	key%list%freq0=freq
1055	key%list%n_bessel=0
1056	key%list%harmon=one
1057
1058	if(key%list%k(1).ne.zero.and.key%list%freq0.ne.zero) icav=1
1059	case default
1060	print*,"Element: ",name," not implemented in PTC"
1061	stop
1062	end select
1063	100 continue
1064	! if(code.ne.14.and.code.ne.15.and.code.ne.16) then
1065	do i=1,NMAX
1066	key%list%k(i)=bvk*key%list%k(i)
1067	key%list%ks(i)=bvk*key%list%ks(i)
1068	enddo
1069	! endif
1070	call create_fibre(my_ring%end,key,EXCEPTION)
1071
1072	if(advance_node().ne.0) goto 10
1073
1074	if (getdebug() > 0) then
1075	print*,' Length of machine: ',l_machine
1076	endif
1077
1078	CALL GET_ENERGY(ENERGY,kin,BRHO,beta0,P0C)
1079
1080	isclosedlayout=get_value('ptc_create_layout ','closed_layout ') .ne. 0
1081
1082	if (getdebug() > 0) then
1083	if ( isclosedlayout .eqv. .true. ) then
1084	print *,'The machine is a RING'
1085	else
1086	print *,'The machine is a LINE'
1087	endif
1088	endif
1089
1090	MY_RING%closed=isclosedlayout
1091
1092	doneit=.true.
1093	call ring_l(my_ring,doneit)
1094
1095	resplit=get_value('ptc_create_layout ','resplit ').ne.0
1096	if(resplit) then
1097	my_thin = get_value('ptc_create_layout ','thin ')
1098	my_xbend = get_value('ptc_create_layout ','xbend ')
1099	even = get_value('ptc_create_layout ','even ').ne.0
1100	resplit_cutting=2
1101	CALL THIN_LENS_resplit(my_ring,THIN=my_thin,even=even,xbend=my_xbend)
1102	endif
1103
1104	if (getdebug() > 0) then
1105	write(6,*) "------------------------------------ PTC Survey ------------------------------------"
1106	write(6,*) "Before start: ",my_ring%start%chart%f%a
1107	write(6,*) "Before end: ",my_ring%end%chart%f%b
1108	endif
1109
1110	call survey(my_ring)
1111
1112	if (getdebug() > 0) then
1113	write(6,*) "After start: ",my_ring%start%chart%f%a
1114	write(6,*) "After end: ",my_ring%end%chart%f%b
1115	endif
1116
1117	call setintstate(default)
1118
1119	if(my_ring%HARMONIC_NUMBER>0) then
1120	call get_length(my_ring,l)
1121	p=>my_ring%start
1122	do i=1,my_ring%n
1123	if(p%mag%kind==kind4) then
1124	if(p%mag%freq==zero) then
1125	write(6,*) " Bullshitting in MADX with Cavities ",my_ring%HARMONIC_NUMBER
1126	p%mag%freq=clightmy_ring%HARMONIC_NUMBERBETA0/l
1127	p%magp%freq=p%mag%freq
1128	endif
1129	endif
1130	p=>p%next
1131	enddo
1132	endif
1133
1134	if (getdebug() > 1) then
1135	print *, '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
1136	print *, '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
1137	print *, '^^^^^^ F I N I S H E D P T C I N P U T ^^^^^^^^'
1138	print *, '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
1139	print *, '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
1140	endif
1141
1142	return
1143
1144	END subroutine ptc_input
1145	!_________________________________________________________________
1146
1147	SUBROUTINE SUMM_MULTIPOLES_AND_ERRORS (l, key, normal_0123, skew_0123,ord_max)
1148	use twtrrfi ! integer, maxmul,maxferr,maxnaper
1149	implicit none
1150	! 1) read multipole coeff. and errors for a current thick element
1151	! 2) fill the error and multiploes arrays of data-bases
1152	REAL(dp), INTENT(IN) :: l
1153	type(keywords), INTENT(INOUT) :: key
1154	REAL(dp), INTENT(OUT) :: normal_0123(0:3), skew_0123(0:3) ! n/l;
1155	REAL(dp) :: normal(0:maxmul), skew (0:maxmul), &
1156	f_errors(0:maxferr), field(2,0:maxmul)
1157	INTEGER :: n_norm, n_skew, n_ferr ! number of terms in command line
1158	INTEGER :: node_fd_errors ! function
1159	integer :: i, i_count, n_dim_mult_err, ord_max
1160
1161	!initialization
1162	normal_0123(:)=zero
1163	skew_0123(:)=zero
1164	do i=1,NMAX
1165	key%list%k(i)=zero
1166	key%list%ks(i)=zero
1167	enddo
1168
1169	! real(dp) f_errors(0:maxferr),normal(0:maxmul),skew(0:maxmul)
1170	! Get multipole components on bench !-----------------------!
1171	call dzero(normal,maxmul+1) ! make zero "normal" !
1172	call dzero(skew,maxmul+1) ! make zero "skew" !
1173	! !
1174	! madxdict.h: "knl = [r, {0}], " !
1175	! "ksl = [r, {0}], " !
1176	! Assign values from the command line !
1177	call get_node_vector('knl ',n_norm,normal) !
1178	call get_node_vector('ksl ',n_skew,skew) !
1179	skew(0)=-skew(0) ! frs error found 30.08.2008
1180	if(n_norm.ge.maxmul) n_norm=maxmul-1 !
1181	if(n_skew.ge.maxmul) n_skew=maxmul-1 !
1182	ord_max=max(n_norm,n_skew) !
1183	! void get_node_vector(charpar,intlength,double* vector) !
1184	! /* returns vector for parameter par of current element */ !
1185	! !
1186	! get errors !
1187	call dzero(f_errors,maxferr+1) !
1188	n_ferr = node_fd_errors(f_errors) ! !
1189	! /* returns the field errors of a node */ !
1190	call dzero(field,2*(maxmul+1)) ! array to be zeroed. !
1191	if (n_ferr .gt. 0) then !
1192	call dcopy(f_errors,field,n_ferr) !
1193	! subroutine dcopy(in,out,n) !
1194	! Purpose: Copy arrays. !
1195	endif !
1196	!-----------------------------------------------------------!
1197
1198	! fill strength of ALL normal multipoles
1199	if(n_norm.gt.0) then ! ===============================!
1200	do i_count=0,n_norm !
1201	if(i_count.gt.0) then !
1202	if(l.ne.zero) then !
1203	key%list%k(i_count+1)=normal(i_count)/l !
1204	else !
1205	key%list%k(i_count+1)=normal(i_count) !
1206	endif !
1207	endif !
1208	if (i_count.le.3) then !
1209	if(l.ne.zero) then !
1210	normal_0123(i_count)=normal(i_count)/l !
1211	else !
1212	normal_0123(i_count)=normal(i_count) !
1213	endif !
1214	endif !
1215	enddo !
1216	endif !================================================!
1217
1218	! fill strength of ALL skew multipoles
1219	if(n_skew.gt.0) then ! ===============================!
1220	do i_count=0,n_skew !
1221	if(i_count.gt.0) then !
1222	if(l.ne.zero) then !
1223	key%list%ks(i_count+1)=skew(i_count)/l !
1224	else !
1225	key%list%ks(i_count+1)=skew(i_count) !
1226	endif !
1227	endif !
1228	if (i_count.le.3) then !
1229	if(l.ne.zero) then !
1230	skew_0123(i_count)=skew(i_count)/l !
1231	else !
1232	skew_0123(i_count)=skew(i_count) !
1233	endif !
1234	endif !
1235	enddo !
1236	endif !================================================!
1237
1238	n_dim_mult_err = max(n_norm, n_skew, n_ferr/2) !===========!
1239	if(n_dim_mult_err.ge.maxmul) n_dim_mult_err=maxmul-1 !
1240	if(n_ferr.gt.0) then !
1241	do i_count=0,n_dim_mult_err !
1242	if(l.ne.zero) then !
1243	key%list%k(i_count+1)=key%list%k(i_count+1)+ & !
1244	field(1,i_count)/l !
1245	key%list%ks(i_count+1)=key%list%ks(i_count+1)+ & !
1246	field(2,i_count)/l !
1247	else !
1248	key%list%k(i_count+1)=key%list%k(i_count+1)+ & !
1249	field(1,i_count) !
1250	key%list%ks(i_count+1)=key%list%ks(i_count+1)+ & !
1251	field(2,i_count) !
1252	endif !
1253	enddo !
1254	endif !====================================================!
1255
1256
1257
1258	END SUBROUTINE SUMM_MULTIPOLES_AND_ERRORS
1259	!----------------------------------------------------------------
1260	!_________________________________________________________________
1261
1262	SUBROUTINE REFRESH_MULTIPOLES(l, normal_0123, skew_0123,ord_max,normal,skew)
1263	use twtrrfi ! integer, maxmul,maxferr,maxnaper
1264	implicit none
1265	! 1) read multipole coeff. and errors for a current thick element
1266	! 2) fill the error and multiploes arrays of data-bases
1267	REAL(dp), INTENT(IN) :: l
1268	REAL(dp), INTENT(OUT) :: normal_0123(0:3), skew_0123(0:3) ! n/l;
1269	REAL(dp) :: normal(0:maxmul), skew (0:maxmul), &
1270	f_errors(0:maxferr), field(2,0:maxmul)
1271	INTEGER :: n_norm, n_skew, n_ferr ! number of terms in command line
1272	INTEGER :: node_fd_errors ! function
1273	integer :: i, i_count, n_dim_mult_err, ord_max
1274
1275	!initialization
1276	normal_0123(:)=zero
1277	skew_0123(:)=zero
1278
1279	! real(dp) f_errors(0:maxferr),normal(0:maxmul),skew(0:maxmul)
1280	! Get multipole components on bench !-----------------------!
1281	call dzero(normal,maxmul+1) ! make zero "normal" !
1282	call dzero(skew,maxmul+1) ! make zero "skew" !
1283	! !
1284	! madxdict.h: "knl = [r, {0}], " !
1285	! "ksl = [r, {0}], " !
1286	! Assign values from the command line !
1287	call get_node_vector('knl ',n_norm,normal) !
1288	call get_node_vector('ksl ',n_skew,skew) !
1289	skew(0)=-skew(0) ! frs error found 30.08.2008
1290	if(n_norm.ge.maxmul) n_norm=maxmul-1 !
1291	if(n_skew.ge.maxmul) n_skew=maxmul-1 !
1292	ord_max=max(n_norm,n_skew) !
1293	if(l.ne.zero) then !
1294	do i=0,maxmul
1295	normal(i)=normal(i)/l
1296	skew(i)=skew(i)/l
1297	enddo
1298	endif
1299	do i=0,3
1300	normal_0123(i)=normal(i)
1301	skew_0123(i)=skew(i)
1302	enddo
1303
1304	! get errors !
1305	call dzero(f_errors,maxferr+1) !
1306	n_ferr = node_fd_errors(f_errors) ! !
1307	! /* returns the field errors of a node */ !
1308	call dzero(field,2*(maxmul+1)) ! array to be zeroed. !
1309	if (n_ferr .gt. 0) then !
1310	call dcopy(f_errors,field,n_ferr) !
1311	endif !
1312	!-----------------------------------------------------------!
1313
1314	n_dim_mult_err = max(n_norm, n_skew, n_ferr/2) !===========!
1315	if(n_dim_mult_err.ge.maxmul) n_dim_mult_err=maxmul-1 !
1316	if(n_ferr.gt.0) then !
1317	do i_count=0,n_dim_mult_err !
1318	if(l.ne.zero) then !
1319	normal(i_count+1)=normal(i_count+1)+field(1,i_count)/l
1320	skew(i_count+1)=skew(i_count+1)+field(2,i_count)/l
1321	else !
1322	normal(i_count+1)=normal(i_count+1)+field(1,i_count)
1323	skew(i_count+1)=skew(i_count+1)+field(2,i_count)
1324	endif !
1325	enddo !
1326	endif !====================================================!
1327
1328
1329
1330	END SUBROUTINE REFRESH_MULTIPOLES
1331	!----------------------------------------------------------------
1332
1333	subroutine ptc_getnfieldcomp(fibreidx, ncomp, nval)
1334	implicit none
1335	real(kind(1d0)) :: nval
1336	integer :: fibreidx
1337	integer :: ncomp
1338	type(fibre), pointer :: p
1339	integer :: j
1340
1341	p=>my_ring%start
1342	do j=1, fibreidx
1343	p=>p%next
1344	enddo
1345
1346	ncomp = ncomp + 1
1347	nval = p%mag%BN(ncomp)
1348
1349	end subroutine ptc_getnfieldcomp
1350	!----------------------------------------------------------------
1351
1352	subroutine ptc_getsfieldcomp(fibreidx, ncomp, nval)
1353	implicit none
1354	real(kind(1d0)) :: nval
1355	integer :: fibreidx
1356	integer :: ncomp
1357	type(fibre), pointer :: p
1358	integer :: j
1359
1360	p=>my_ring%start
1361	do j=1, fibreidx
1362	p=>p%next
1363	enddo
1364
1365	ncomp = ncomp + 1
1366
1367	nval = p%mag%AN(ncomp)
1368	print*, "Returning AN",nval," for ",p%mag%name
1369
1370
1371	end subroutine ptc_getsfieldcomp
1372	!----------------------------------------------------------------
1373
1374	subroutine ptc_setfieldcomp(fibreidx)
1375	implicit none
1376	integer :: fibreidx
1377	type(fibre), pointer :: p
1378	integer :: j, i
1379	integer :: kn, ks
1380	real(dp) :: v
1381	real(kind(1d0)) get_value
1382
1383	if ( .not. associated(my_ring) ) then
1384	call fort_warn("ptc_setfieldcomp","No active PTC layout/period")
1385	return
1386	endif
1387
1388	if (getdebug()>2) then
1389	print*, "I am in ptc_setfieldcomp: Element index is ", fibreidx
1390	endif
1391
1392	if ( (fibreidx .lt. 1) .and. (fibreidx .gt. my_ring%n) ) then
1393	call fort_warn("ptc_setfieldcomp","element out of range of the current layout")
1394	return
1395	endif
1396
1397	p=>my_ring%start
1398	do j=1, fibreidx
1399	p=>p%next
1400	enddo
1401
1402	if (getdebug() > 1 ) then
1403	print*,"Found element no. ", fibreidx," named ", p%mag%name, &
1404	&" of kind ", p%mag%kind, mytype(p%mag%kind)
1405	print*,"Currently nmul is ", p%mag%p%nmul
1406
1407	write(6,*) "BNs",p%mag%BN
1408	write(6,*) "ANs",p%mag%AN
1409
1410	DO i=1,p%mag%p%nmul
1411	print*, "Polimorphic BN(",i,")"
1412	call print(p%mag%BN(i),6)
1413	print*, "Polimorphic AN(",i,")"
1414	call print(p%mag%AN(i),6)
1415	ENDDO
1416
1417	endif
1418
1419	kn = get_value('ptc_setfieldcomp ','kn ')
1420	v = get_value('ptc_setfieldcomp ','value ')
1421
1422	if (kn >= 0) then
1423	kn = kn + 1
1424
1425	if (getdebug() > 1) then
1426	print*,"Setting up KN ", kn, " from ", p%mag%BN(kn) ," to ", v
1427	endif
1428
1429	call add(p%mag, kn,0,v)
1430	call add(p%magp,kn,0,v)
1431
1432
1433	else
1434	ks = get_value('ptc_setfieldcomp ','ks ')
1435	if (ks < 0) then
1436	call fort_warn("ptc_setfieldcomp","neither kn nor ks specified")
1437	return
1438	endif
1439	ks = ks + 1
1440
1441	! print*,"Setting up skew field component ", ks," to ", v
1442
1443	if (getdebug() > 1) then
1444	print*,"Setting up KS ", ks, " from ", p%mag%AN(ks) ," to ", v
1445	endif
1446	call add(p%mag, -ks,0,v)
1447	call add(p%magp,-ks,0,v)
1448
1449	endif
1450
1451	if (getdebug() > 1 ) then
1452	write(6,*) "BNs",p%mag%BN
1453	write(6,*) "ANs",p%mag%AN
1454	write(6,*) ""
1455	endif
1456	end subroutine ptc_setfieldcomp
1457	!----------------------------------------------------------------
1458
1459	subroutine ptc_align()
1460	use twiss0fi
1461	implicit none
1462	integer j,n_align,node_al_errors
1463	integer restart_sequ,advance_node
1464	real(dp) al_errors(align_max)
1465	type(fibre), pointer :: f
1466	!---------------------------------------------------------------
1467
1468	j=restart_sequ()
1469	j=0
1470	f=>my_ring%start
1471	10 continue
1472
1473	j=j+1
1474	n_align = node_al_errors(al_errors)
1475	if (n_align.ne.0) then
1476	!write(6,'(6f11.8)') al_errors(1:6)
1477	call mad_misalign_fibre(f,al_errors(1:6))
1478	endif
1479	f=>f%next
1480	if(advance_node().ne.0) goto 10
1481
1482	END subroutine ptc_align
1483	!_________________________________________________________________
1484
1485	subroutine ptc_dumpmaps()
1486	!Dumps to file maps and/or matrixes (i.e. first order maps)
1487	implicit none
1488	type(fibre), pointer :: p
1489	type(damap) :: id !identity map used for calculating maps for each element
1490	type(real_8) :: y2(6) !polimorphes array used for calculating maps for each element
1491	type(real_8) :: yfull(6) !polimorphes array used for calculating maps for each element
1492	real(dp) :: xt(6)
1493	integer :: i !iterators
1494	integer mf1,mf2
1495	character(200) :: filename='ptcmaps.txt'
1496	character(200) :: filenamefull='ptcmaps'
1497	integer :: flag_index,why(9)
1498	character(200) :: whymsg
1499	real(kind(1d0)) :: suml=zero
1500	integer geterrorflag !C function that returns errorflag value
1501
1502	suml=zero
1503
1504	if (cavsareset .eqv. .false.) then
1505	call setcavities(my_ring,maxaccel)
1506	if (geterrorflag() /= 0) then
1507	return
1508	endif
1509	endif
1510
1511	if (getdebug() > 1) then
1512	print *, '<madx_ptc_module.f90 : ptc_dumpmaps> Maps are dumped to file ',filename
1513	endif
1514	call kanalnummer(mf1)
1515	open(unit=mf1,file=filename)
1516
1517	! write(filenamefull,*) filename,".",my_ring%start%mag%name,"-",my_ring%end%mag%name,".txt"
1518	filenamefull="ptcmaps.start-end.txt"
1519	print*, filenamefull
1520	call kanalnummer(mf2)
1521	open(unit=mf2,file=filenamefull)
1522
1523	print*, "no=1"," mynd2=",c_%nd2," npara=",c_%npara
1524	call init(getintstate(),1,c_%np_pol,berz)
1525
1526	call alloc(id);
1527	call alloc(y2);
1528	call alloc(yfull);
1529
1530	xt(:) = zero
1531	id = 1 ! making identity map
1532
1533	yfull = xt + id
1534
1535	p=>my_ring%start
1536	do i=1,my_ring%n
1537
1538
1539	y2=xt+id ! we track identity map from the current position
1540
1541	if( (p%mag%kind/=kind21) .and. (p%mag%kind/=kind4) ) then
1542
1543	call track(my_ring,y2,i,i+1,getintstate())
1544
1545	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
1546	call fort_warn('ptc_dumpmaps: ','DA got unstable')
1547	call seterrorflag(10,"ptc_dumpmaps ","DA got unstable ");
1548	close(mf1)
1549	close(mf2)
1550	return
1551	endif
1552
1553	call PRODUCE_APERTURE_FLAG(flag_index)
1554	if(flag_index/=0) then
1555	call ANALYSE_APERTURE_FLAG(flag_index,why)
1556
1557	Write(6,*) "ptc_dumpmaps: APERTURE error for element: ",i," name: ",p%MAG%name
1558	write(whymsg,*) 'APERTURE error: ',why
1559	call fort_warn('ptc_dumpmaps: ',whymsg)
1560	call seterrorflag(10,"ptc_dumpmaps: ",whymsg);
1561	c_%watch_user=.false.
1562	close(mf1)
1563	close(mf2)
1564	return
1565	endif
1566
1567	call track(my_ring,xt,i,i+1,getintstate())
1568	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
1569	call fort_warn('ptc_dumpmaps: ','DA got unstable')
1570	call seterrorflag(10,"ptc_dumpmaps ","DA got unstable ");
1571	close(mf1)
1572	close(mf2)
1573	return
1574	endif
1575
1576	call PRODUCE_APERTURE_FLAG(flag_index)
1577	if(flag_index/=0) then
1578	call ANALYSE_APERTURE_FLAG(flag_index,why)
1579	Write(6,*) "ptc_dumpmaps: APERTURE error for element: ",i," name: ",p%MAG%name
1580	write(whymsg,*) 'APERTURE error: ',why
1581	call fort_warn('ptc_dumpmaps: ',whymsg)
1582	call seterrorflag(10,"ptc_dumpmaps: ",whymsg);
1583	c_%watch_user=.false.
1584	close(mf1)
1585	close(mf2)
1586	return
1587	endif
1588	else
1589	if (getdebug() > 2) then
1590	print *, 'Track Cavity...'
1591	endif
1592
1593	call track(my_ring,y2,i,i+2,getintstate())
1594	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
1595	call fort_warn('ptc_dumpmaps: ','DA got unstable')
1596	call seterrorflag(10,"ptc_dumpmaps ","DA got unstable ");
1597	close(mf1)
1598	close(mf2)
1599	return
1600	endif
1601
1602	call PRODUCE_APERTURE_FLAG(flag_index)
1603	if(flag_index/=0) then
1604	call ANALYSE_APERTURE_FLAG(flag_index,why)
1605	Write(6,*) "ptc_dumpmaps: APERTURE error for element: ",i," name: ",p%MAG%name
1606	write(whymsg,*) 'APERTURE error: ',why
1607	call fort_warn('ptc_dumpmaps: ',whymsg)
1608	call seterrorflag(10,"ptc_dumpmaps: ",whymsg);
1609	c_%watch_user=.false.
1610	close(mf1)
1611	close(mf2)
1612	return
1613	endif
1614
1615	call track(my_ring,xt,i,i+2,getintstate())
1616	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
1617	call fort_warn('ptc_dumpmaps: ','DA got unstable')
1618	call seterrorflag(10,"ptc_dumpmaps ","DA got unstable ");
1619	close(mf1)
1620	close(mf2)
1621	return
1622	endif
1623
1624	call PRODUCE_APERTURE_FLAG(flag_index)
1625	if(flag_index/=0) then
1626	call ANALYSE_APERTURE_FLAG(flag_index,why)
1627	Write(6,*) "ptc_dumpmaps: APERTURE error for element: ",i," name: ",p%MAG%name
1628	write(whymsg,*) 'APERTURE error: ',why
1629	call fort_warn('ptc_dumpmaps: ',whymsg)
1630	call seterrorflag(10,"ptc_dumpmaps: ",whymsg);
1631	c_%watch_user=.false.
1632	close(mf1)
1633	close(mf2)
1634	return
1635	endif
1636	endif
1637
1638
1639	call track(my_ring,yfull,i,i+1,getintstate())
1640
1641	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
1642	call fort_warn('ptc_dumpmaps: ','DA got unstable')
1643	call seterrorflag(10,"ptc_dumpmaps ","DA got unstable ");
1644	close(mf1)
1645	close(mf2)
1646	return
1647	endif
1648
1649	call PRODUCE_APERTURE_FLAG(flag_index)
1650	if(flag_index/=0) then
1651	call ANALYSE_APERTURE_FLAG(flag_index,why)
1652
1653	Write(6,*) "ptc_dumpmaps: APERTURE error for element: ",i," name: ",p%MAG%name
1654	write(whymsg,*) 'APERTURE error: ',why
1655	call fort_warn('ptc_dumpmaps: ',whymsg)
1656	call seterrorflag(10,"ptc_dumpmaps: ",whymsg);
1657	c_%watch_user=.false.
1658	close(mf1)
1659	close(mf2)
1660	return
1661	endif
1662
1663	write(mf2,*) p%mag%name, suml,' m ==========================='
1664	call print(yfull,mf2)
1665
1666	suml=suml+p%MAG%P%ld
1667
1668	write(mf1,*) p%mag%name, suml,' m ==========================='
1669	if (c_%npara == 6) then
1670	call dump6dmap(y2, mf1)
1671	elseif (c_%npara == 5) then
1672	call dump5dmap(y2, mf1)
1673	elseif (c_%npara == 4) then
1674	call dump4dmap(y2, mf1)
1675	else
1676	call fort_warn("ptc_dumpmaps","c_%npara is neither 6,5 nor 4")
1677	endif
1678	p=>p%next
1679	enddo
1680
1681	close(mf1)
1682	call kill(y2);
1683	call kill(id);
1684
1685	!_________________________________________________________________
1686	!_________________________________________________________________
1687	!_________________________________________________________________
1688
1689	contains
1690	!_________________________________________________________________
1691	subroutine dump4dmap(y2, fun)
1692	implicit none
1693	double precision a1000,a0100,a0010,a0001
1694	type(real_8) :: y2(6) !polimorphes array used for calculating maps for each element
1695	integer :: fun !file unit number
1696	integer :: ii
1697
1698	if (getdebug() > 1) then
1699
1700	endif
1701
1702	do ii=1,4
1703	a1000=y2(ii).sub.'1000'
1704	a0100=y2(ii).sub.'0100'
1705	a0010=y2(ii).sub.'0010'
1706	a0001=y2(ii).sub.'0001'
1707	write(fun,'(6f13.8)') a1000, &
1708	& a0100, &
1709	& a0010, &
1710	& a0001
1711	enddo
1712
1713	end subroutine dump4dmap
1714	!_________________________________________________________________
1715
1716	subroutine dump5dmap(y2, fun)
1717	implicit none
1718	double precision a10000,a01000,a00100,a00010,a00001
1719	type(real_8) :: y2(6) !polimorphes array used for calculating maps for each element
1720	integer :: fun !file unit number
1721	integer :: ii
1722	do ii=1,5
1723	a10000=y2(ii).sub.'10000'
1724	a01000=y2(ii).sub.'01000'
1725	a00100=y2(ii).sub.'00100'
1726	a00010=y2(ii).sub.'00010'
1727	a00001=y2(ii).sub.'00001'
1728	write(fun,'(6f13.8)') a10000, &
1729	& a01000, &
1730	& a00100, &
1731	& a00010, &
1732	& a00001 !
1733	enddo
1734
1735	end subroutine dump5dmap
1736	!_________________________________________________________________
1737
1738	subroutine dump6dmap(y2, fun)
1739	implicit none
1740	double precision a100000,a010000,a001000,a000100,a000010,a000001
1741	type(real_8) :: y2(6) !polimorphes array used for calculating maps for each element
1742	integer :: fun !file unit number
1743	integer :: ii
1744
1745	do ii=1,4
1746	a100000=y2(ii).sub.'100000'
1747	a010000=y2(ii).sub.'010000'
1748	a001000=y2(ii).sub.'001000'
1749	a000100=y2(ii).sub.'000100'
1750	a000010=y2(ii).sub.'000010'
1751	a000001=y2(ii).sub.'000001'
1752	write(fun,'(6f13.8)') a100000, &
1753	& a010000, &
1754	& a001000, &
1755	& a000100, &
1756	& a000001, & !madx format has dp/p at the last column
1757	& a000010 !
1758	enddo
1759
1760	do ii=6,5,-1
1761	a100000=y2(ii).sub.'100000'
1762	a010000=y2(ii).sub.'010000'
1763	a001000=y2(ii).sub.'001000'
1764	a000100=y2(ii).sub.'000100'
1765	a000010=y2(ii).sub.'000010'
1766	a000001=y2(ii).sub.'000001'
1767	write(fun,'(6f13.8)') a100000, &
1768	& a010000, &
1769	& a001000, &
1770	& a000100, &
1771	& a000001, & !madx format has dp/p at the last column
1772	& a000010 !
1773	enddo
1774	end subroutine dump6dmap
1775
1776
1777	end subroutine ptc_dumpmaps
1778
1779	RECURSIVE FUNCTION FACTORIAL (N) &
1780	RESULT (FACTORIAL_RESULT)
1781	INTEGER :: N, FACTORIAL_RESULT
1782
1783	IF (N <= 0 ) THEN
1784	FACTORIAL_RESULT = 1
1785	ELSE
1786	FACTORIAL_RESULT = N * FACTORIAL (N-1)
1787	END IF
1788	END FUNCTION FACTORIAL
1789
1790	subroutine ptc_track()
1791	implicit none
1792	integer i,nint,ndble,nchar,int_arr(1),char_l,icase,turns,flag_index,why(9)
1793	integer j,next_start
1794	real(dp) x0(6),x(6),deltap0,deltap,dt
1795	real(dp) xx,pxx,yx,pyx,tx,deltaex,fxx,phixx,fyx,phiyx,ftx,phitx
1796	real(kind(1d0)) get_value
1797	logical(lp) closed_orbit
1798	character*12 char_a
1799	data char_a / ' ' /
1800	!------------------------------------------------------------------------------
1801
1802
1803	if(universe.le.0.or.EXCEPTION.ne.0) then
1804	call fort_warn('return from ptc_track: ',' no universe created')
1805	return
1806	endif
1807	if(index_mad.le.0.or.EXCEPTION.ne.0) then
1808	call fort_warn('return from ptc_track: ',' no layout created')
1809	return
1810	endif
1811
1812	icase = get_value('ptc_track ','icase ')
1813	deltap0 = get_value('ptc_track ','deltap ')
1814
1815	deltap = zero
1816	call my_state(icase,deltap,deltap0)
1817
1818	if (getdebug() > 2) then
1819	print *, "ptc_track: internal state is:"
1820	call print(default,6)
1821	endif
1822
1823	x0(:)=zero
1824	if(mytime) then
1825	call Convert_dp_to_dt (deltap, dt)
1826	else
1827	dt=deltap
1828	endif
1829	if(icase.eq.5) x0(5)=dt
1830	closed_orbit = get_value('ptc_track ','closed_orbit ') .ne. 0
1831	if(closed_orbit) then
1832	call find_orbit(my_ring,x0,1,default,c_1d_7)
1833	CALL write_closed_orbit(icase,x0)
1834	endif
1835
1836
1837	call comm_para('coord ',nint,ndble,nchar,int_arr,x,char_a,char_l)
1838
1839	j = next_start(xx,pxx,yx,pyx,tx,deltaex,fxx,phixx,fyx,phiyx,ftx,phitx)
1840	print*,"dat1",j,xx,pxx,yx,pyx,tx,deltaex,fxx,phixx,fyx,phiyx,ftx,phitx
1841	j = next_start(xx,pxx,yx,pyx,tx,deltaex,fxx,phixx,fyx,phiyx,ftx,phitx)
1842	print*,"dat2",j,xx,pxx,yx,pyx,tx,deltaex,fxx,phixx,fyx,phiyx,ftx,phitx
1843	j = next_start(xx,pxx,yx,pyx,tx,deltaex,fxx,phixx,fyx,phiyx,ftx,phitx)
1844	print*,"dat3",j,xx,pxx,yx,pyx,tx,deltaex,fxx,phixx,fyx,phiyx,ftx,phitx
1845
1846	x(:)=x(:)+x0(:)
1847	print*," Initial Coordinates: ", x
1848	turns = get_value('ptc_track ','turns ')
1849	c_%watch_user=.true.
1850	do i=1,turns
1851	call track(my_ring,x,1,default)
1852	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
1853	call fort_warn('ptc_track: ','DA got unstable')
1854	call seterrorflag(10,"ptc_track ","DA got unstable ");
1855	return
1856	endif
1857	call PRODUCE_APERTURE_FLAG(flag_index)
1858	if(flag_index/=0) then
1859	call ANALYSE_APERTURE_FLAG(flag_index,why)
1860	Write(6,*) "ptc_track unstable (tracking)-programs continues "
1861	Write(6,*) why ! See produce aperture flag routine in sd_frame
1862	goto 100
1863	endif
1864	enddo
1865	c_%watch_user=.false.
1866	print*," End Coordinates: ",x
1867	return
1868	100 continue
1869	c_%watch_user=.false.
1870	print*," Last Coordinates: ",x," after: ",i," turn(s)"
1871
1872	END subroutine ptc_track
1873
1874
1875
1876	!________________________________________________________________________________
1877
1878
1879	subroutine ptc_end()
1880	implicit none
1881	integer i
1882
1883	if(universe.le.0.or.EXCEPTION.ne.0) then
1884	call fort_warn('return from ptc_end: ',' no universe can be killed')
1885	return
1886	endif
1887
1888	call killsavedmaps() !module ptc_twiss -> kill buffered maps
1889
1890	! call killparresult()
1891	call resetknobs() !remove the knobs
1892
1893	if ( associated(m_u%n) .eqv. .false. ) then
1894	print*, "We attempt to kill not initialized universe!"
1895	endif
1896
1897
1898	call kill_universe(m_u)
1899	nullify(my_ring)
1900	call kill_tpsa
1901	do i=1,size(s_b)
1902	call nul_coef(s_b(i))
1903	enddo
1904	deallocate(s_b)
1905	firsttime_coef=.true.
1906
1907	universe=universe-1
1908	end subroutine ptc_end
1909
1910
1911	subroutine normalform_normalform(s1,s2)
1912	implicit none
1913	type (normalform),intent(inout)::s1
1914	type (normalform),intent(in)::s2
1915	integer i,j
1916
1917	s1%a_t=s2%a_t
1918	s1%a1=s2%a1
1919	s1%a%constant(:)=s2%a%constant(:)
1920	s1%a%Linear=s2%a%Linear
1921	s1%a%nonlinear=s2%a%nonlinear
1922	s1%a%pb=s2%a%pb
1923	s1%normal%constant(:)=s2%normal%constant(:)
1924	s1%normal%Linear=s2%normal%Linear
1925	s1%normal%nonlinear=s2%normal%nonlinear
1926	s1%normal%pb=s2%normal%pb
1927	s1%DHDJ=s2%DHDJ
1928	do i=1,ndim
1929	s1%TUNE(i)=s2%TUNE(i)
1930	s1%damping(i)=s2%damping(i)
1931	s1%plane(i)=s2%plane(i)
1932	do j=1,mynreso
1933	s1%m(i,j)=s2%m(i,j)
1934	enddo
1935	enddo
1936	s1%nord=s2%nord
1937	s1%jtune=s2%jtune
1938	s1%nres=s2%nres
1939	s1%AUTO=s2%AUTO
1940	end subroutine normalform_normalform
1941	!_________________________________________________________________
1942
1943
1944	SUBROUTINE set_PARAMETERS(R,nt,iorder,IFAM,inda,scale)
1945	!Strength of Multipole of order iorder as parameter
1946
1947	IMPLICIT NONE
1948	integer ipause, mypause
1949	logical(lp) ok
1950	INTEGER iorder,i,j,jj,k,lstr,IFAM,tot,nt,inda,min1
1951	INTEGER,parameter::ipara=100
1952	real(dp) scale(ipara),value
1953	character(20) str
1954	CHARACTER(3) STR1
1955	character(10),dimension(10)::multname
1956	type(layout) r
1957	type(fibre), POINTER :: current
1958	INTEGER,ALLOCATABLE,dimension(:)::DAFAM
1959	INTEGER,ALLOCATABLE,dimension(:,:)::FAM
1960	real(dp),ALLOCATABLE,dimension(:)::SFAM
1961	multname=(/"Dipole ","Quadrupole","Sextupole ","Octupole ","Decapole ",&
1962	"Dodecapole","14-Pole ","16-Pole ","18-Pole ","20-Pole "/)
1963
1964	ALLOCATE(FAM(IFAM,0:R%N),DAFAM(IFAM),SFAM(IFAM))
1965
1966	min1=0
1967	if(iorder.lt.0) then
1968	min1=1
1969	iorder=-iorder
1970	endif
1971	DO I=1,IFAM
1972	OK=.TRUE.
1973	DO WHILE(OK)
1974	TOT=0
1975	if(min1.eq.0) WRITE(6,*) " Identify ",multname(iorder)
1976	if(min1.eq.1) WRITE(6,*) " Identify ","SKEW-"//multname(iorder)
1977	READ(5,*) STR
1978	STR=TRIM(ADJUSTL(STR))
1979	LSTR=LEN_TRIM (STR)
1980	current=>r%start
1981	DO J=1,R%N
1982	IF(current%MAG%NAME==STR.and.current%MAG%P%NMUL==iorder) THEN
1983	TOT=TOT+1
1984	FAM(I,TOT)=J
1985	ENDIF
1986	current=>current%next
1987	ENDDO
1988	WRITE(6,*) TOT," Is that OK? YES or NO?"
1989	READ (5,*) STR1
1990	STR1=TRIM(ADJUSTL(STR1))
1991	IF(STR1(1:1)=='Y'.OR.STR1(1:1)=='y') THEN
1992	OK=.FALSE.
1993	inda=inda+1
1994	if(inda.gt.100) then
1995	write(6,*) " Problem: Only ",ipara," Parameters allowed"
1996	ipause=mypause(2002)
1997	endif
1998	DAFAM(I)=inda
1999	WRITE(6,*) " Give Scaling Factor, '0' uses Default"
2000	read(5,*) value
2001	if(value==0) then
2002	WRITE(6,*) " Take Default Scaling Value : ",scale(inda)
2003	SFAM(I)=scale(inda)
2004	else
2005	SFAM(I)=value
2006	endif
2007	ENDIF
2008	ENDDO
2009
2010	FAM(I,0)=TOT
2011	current=r%start
2012	DO JJ=1,FAM(I,0)
2013	J=FAM(I,JJ)
2014	! ALLOCATION GYMNASTIC IF Multipole NOT YET ALLOCATED
2015	IF(current%MAGP%P%NMUL<iorder) THEN
2016	CALL KILL(current%MAGP%BN,current%MAGP%P%NMUL)
2017	CALL KILL(current%MAGP%AN,current%MAGP%P%NMUL)
2018	current%MAGP%P%NMUL=iorder
2019	DEALLOCATE(current%MAGP%BN)
2020	DEALLOCATE(current%MAGP%AN)
2021	CALL ALLOC(current%MAGP%BN,iorder)
2022	CALL ALLOC(current%MAGP%AN,iorder)
2023	ALLOCATE(current%MAGP%BN(iorder),current%MAGP%AN(iorder))
2024	DO K=1,current%MAG%P%NMUL
2025	current%MAGP%BN(K)=current%MAG%BN(K)
2026	current%MAGP%AN(K)=current%MAG%AN(K)
2027	ENDDO
2028	DEALLOCATE(current%MAG%BN)
2029	DEALLOCATE(current%MAG%AN)
2030	ALLOCATE(current%MAG%BN(iorder),current%MAG%AN(iorder))
2031	call equal(current%MAG,current%MAGP)
2032	ENDIF
2033	if(min1.eq.0) then
2034	current%MAGP%BN(iorder)%I=NT+I
2035	current%MAGP%BN(iorder)%KIND=3
2036	else
2037	current%MAGP%AN(iorder)%I=NT+I
2038	current%MAGP%AN(iorder)%KIND=3
2039	endif
2040	current=>current%next
2041	ENDDO
2042	ENDDO
2043
2044	current=r%start
2045	DO I=1,IFAM
2046	DO JJ=1,1
2047	J=FAM(I,JJ)
2048	if(min1.eq.0) WRITE(6,*) current%MAG%NAME,' ', current%MAG%BN(iorder)
2049	if(min1.eq.1) WRITE(6,*) current%MAG%NAME,' ', current%MAG%AN(iorder)
2050	current=>current%next
2051	ENDDO
2052	ENDDO
2053
2054	DEALLOCATE(FAM,STAT=I)
2055	! WRITE(6,*) I
2056	DEALLOCATE(DAFAM,STAT=I)
2057	! WRITE(6,*) I
2058	DEALLOCATE(SFAM,STAT=I)
2059	! WRITE(6,*) I
2060
2061	end subroutine set_PARAMETERS
2062	!______________________________________________________________________
2063
2064	subroutine my_state(icase,deltap,deltap0)
2065	implicit none
2066	integer icase,i
2067	real(dp) deltap0,deltap
2068
2069	default = getintstate()
2070
2071	if (getdebug()>1) then
2072	print*, "icase=",icase," deltap=",deltap," deltap0=",deltap0
2073	endif
2074
2075	deltap = zero
2076	select case(icase)
2077	CASE(4)
2078	if (getdebug()>1) then
2079	print*, "my_state: Enforcing ONLY_4D+NOCAVITY and NO DELTA"
2080	endif
2081	default=default-delta0
2082	default=default+only_4d0+NOCAVITY0
2083	i=4
2084	CASE(5)
2085	if (getdebug()>1) then
2086	print*, "my_state: Enforcing DELTA"
2087	endif
2088	default=default+delta0
2089	deltap=deltap0
2090	i=5
2091	CASE(56)
2092	if (getdebug()>1) then
2093	print*, "my_state: Enforcing coasting beam"
2094	endif
2095	default = default - delta0 - only_4d0
2096	default = default + NOCAVITY0
2097	deltap=deltap0
2098	i=56
2099	CASE(6)
2100	i=6
2101	CASE DEFAULT
2102	default=default+only_4d0+NOCAVITY0
2103	i=4
2104	END SELECT
2105
2106	if (i==6) then
2107	if ( (icav==0) .and. my_ring%closed .and. (getenforce6D() .eqv. .false.)) then
2108	default = default - delta0 - only_4d0
2109	default=default + NOCAVITY0
2110	call fort_warn('return mystate: ',' no cavity - dimensionality reduced 6 -> 5 and 1/2')
2111	i=56
2112	else
2113	default = default - delta0 - only_4d0 - NOCAVITY0 !enforcing nocavity to false
2114	endif
2115	endif
2116
2117	call setintstate(default)
2118	CALL UPDATE_STATES
2119
2120	if (getdebug()>0) call print(default,6)
2121
2122	icase = i
2123
2124	end subroutine my_state
2125
2126	!______________________________________________________________________
2127
2128	subroutine f90flush(i,option)
2129	implicit none
2130	integer i,ios
2131	logical(lp) ostat, fexist,option
2132	logical fexist1, ostat1
2133	character*20 faction,faccess,fform,fwstat,fposition
2134	character*255 fname
2135	inquire(err=1,iostat=ios,&
2136	unit=i,opened=ostat1,exist=fexist1,write=fwstat)
2137	fexist = fexist1
2138	ostat = ostat1
2139	if (.not.ostat.or..not.fexist.or.fwstat.ne.'YES') return
2140	inquire(err=2,iostat=ios,&
2141	unit=i,action=faction,access=faccess,&
2142	form=fform,name=fname,position=fposition)
2143	close (unit=i,err=3)
2144	! write (,) 'Re-opening ',i,' ',option,ios,faction,faccess,fform,fposition,fname
2145	if (option) then
2146	open(err=4,iostat=ios,&
2147	unit=i,action=faction,access=faccess,form=fform,&
2148	file=fname,status='old',position='append')
2149	else
2150	open(err=4,iostat=ios,&
2151	unit=i,action=faction,access=faccess,form=fform,&
2152	file=fname,status='replace',position='rewind')
2153	endif
2154	return
2155	1 write (,)&
2156	' F90FLUSH 1st INQUIRE FAILED with IOSTAT ',ios,' on UNIT ',i
2157	stop
2158	2 write (,)&
2159	' F90FLUSH 2nd INQUIRE FAILED with IOSTAT ', ios,' on UNIT ',i
2160	stop
2161	3 write (,)&
2162	' F90FLUSH CLOSE FAILED with IOSTAT ',ios,' on UNIT ',i
2163	stop
2164	4 write (,)&
2165	' F90FLUSH RE-OPEN FAILED with IOSTAT ',ios,' on UNIT ',i
2166	stop
2167	end subroutine f90flush
2168
2169	SUBROUTINE write_closed_orbit(icase,x)
2170	implicit none
2171	INTEGER, INTENT(IN):: icase
2172	REAL (dp),INTENT(IN) :: x(6)
2173	if(icase.eq.4) then
2174	print*,"Closed orbit: ",x(1),x(2),x(3),x(4)
2175	elseif(icase.eq.5) then
2176	print*,"Closed orbit: ",x(1),x(2),x(3),x(4),x(5)
2177	elseif(icase.eq.6) then
2178	print*,"Closed orbit: ",x(1),x(2),x(3),x(4),-x(6),x(5)
2179	endif
2180	ENDSUBROUTINE write_closed_orbit
2181
2182	SUBROUTINE Convert_dp_to_dt(deltap, dt)
2183	implicit none
2184	! convert deltap=(p-p0)/p0 to dt=deltaE/p0c
2185	REAL(dp), INTENT(IN) :: deltap
2186	REAL(dp), INTENT(OUT) :: dt
2187
2188	! local
2189	real(dp) :: MASS_GeV, ENERGY,KINETIC,BRHO,BETA0,P0C,gamma0I,gambet
2190
2191	! to get "energy" value
2192	Call GET_ONE(MASS_GeV,ENERGY,KINETIC,BRHO,BETA0,P0C,gamma0I,gambet)
2193
2194	IF (beta0.gt.zero ) THEN
2195	dt=SQRT(deltap*(deltap+two)+one/beta0/beta0)-one/beta0
2196	ELSE ! exculde devision by 0
2197	call aafail('SUBR. Convert_dp_to_dt: ',' CALL GET_ONE => beta0.LE.0')
2198	ENDIF
2199
2200	END SUBROUTINE Convert_dp_to_dt
2201	!=============================================================================
2202
2203	subroutine makemaptable(y)
2204	implicit none
2205	type(real_8):: y(6)
2206	integer,parameter :: i_map_coor=10
2207	integer :: map_term, ja(6),i,ii,iii
2208	integer :: i1,i2,i3,i4,i5,i6
2209	integer :: order, no
2210	real(dp) :: coef
2211	real(kind(1d0)) :: map_coor(i_map_coor)
2212	real(kind(1d0)) :: get_value
2213	! type(universal_taylor) :: ut
2214
2215	! write(0,*) "MAP_TABLE"
2216
2217	map_term=42
2218	call make_map_table(map_term)
2219
2220	order = get_value("ptc_normal ","no ")
2221
2222	call liepeek(iia,icoast)
2223	allocate(j(c_%npara))
2224	ja(:) = 0
2225	j(:) = 0
2226
2227	goto 100 ! skip the code that was in place until 29 March 2010
2228
2229	do iii=1,c_%npara
2230	coef = y(iii)%T.sub.j
2231	! following works
2232	!coef = y(iii)%T.sub.mapSelector5variables(1)
2233	map_coor(1)=coef
2234	map_coor(2)=iii
2235	map_coor(3)=c_%npara
2236	map_coor(4)=0
2237	map_coor(5)=ja(1)
2238	map_coor(6)=ja(2)
2239	map_coor(7)=ja(3)
2240	map_coor(8)=ja(4)
2241	map_coor(9)=ja(5)
2242	map_coor(10)=ja(6)
2243	call vector_to_table_curr("ptc_normal ", 'coef ', map_coor(1), i_map_coor)
2244	call augment_count("ptc_normal ")
2245	enddo
2246
2247	do i = 1,c_%npara
2248
2249	do ii = 1,c_%npara
2250	j(ii) = 1
2251	ja(ii) = j(ii)
2252	coef = y(i)%T.sub.j
2253	map_coor(1)=coef
2254	map_coor(2)=i
2255	map_coor(3)=c_%npara! 29.06.2006 here was iia(2) - to be verified
2256	map_coor(4)=sum(ja(:))
2257	map_coor(5)=ja(1)
2258	map_coor(6)=ja(2)
2259	map_coor(7)=ja(3)
2260	map_coor(8)=ja(4)
2261	map_coor(9)=ja(5)
2262	map_coor(10)=ja(6)
2263	call vector_to_table_curr("ptc_normal ", 'coef ', map_coor(1), i_map_coor)
2264	call augment_count("ptc_normal ")
2265	j(:) = 0
2266	ja(ii) = j(ii)
2267	enddo
2268
2269	! ut = y(i)
2270	! do ii = 1,ut%n
2271	! map_coor(1)=ut%c(ii) !coef
2272	! map_coor(2)=i !index of taylor
2273	! map_coor(3)=c_%npara
2274	! map_coor(4)=sum(ut%j(i,:)) !order
2275	! map_coor(5)=ut%j(ii,1)
2276	! map_coor(6)=ut%j(ii,2)
2277	! map_coor(7)=ut%j(ii,3)
2278	! map_coor(8)=ut%j(ii,4)
2279	! map_coor(9)=ut%j(ii,5)
2280	! map_coor(10)=ut%j(ii,6)
2281	! enddo
2282
2283
2284	enddo
2285
2286	! note that the order in which the coefficients appear in the map_table slightly
2287	! differ from the order in which they appear in fort.18
2288	100 do i=1,c_%npara ! distribute exponents over 6 variables, knowing their sum
2289	do no=0,order
2290	if (c_%npara.eq.6) then
2291	do i1=no,0,-1
2292	do i2=no-i1,0,-1
2293	do i3=no-i1-i2,0,-1
2294	do i4=no-i1-i2-i3,0,-1
2295	do i5=no-i1-i2-i3-i4,0,-1
2296	do i6=no-i1-i2-i3-i4-i5,0,-1
2297	if (i1+i2+i3+i4+i5+i6==no) then
2298	!write(0,'(6(i4))'), i1,i2,i3,i4,i5,i6
2299	j(1)=i1
2300	j(2)=i2
2301	j(3)=i3
2302	j(4)=i4
2303	j(5)=i5
2304	j(6)=i6
2305	coef = y(i)%T.sub.j
2306	if (coef.ne.zero) then
2307	map_coor(1)=coef
2308	map_coor(2)=i
2309	map_coor(3)=c_%npara
2310	map_coor(4)=no
2311	map_coor(5)=j(1)
2312	map_coor(6)=j(2)
2313	map_coor(7)=j(3)
2314	map_coor(8)=j(4)
2315	map_coor(9)=j(5)
2316	map_coor(10)=j(6)
2317	call vector_to_table_curr("map_table ", 'coef ', map_coor(1), i_map_coor)
2318	call augment_count("map_table ")
2319	endif
2320	!write(0,*) 'write coef', coef
2321	endif
2322	enddo
2323	enddo
2324	enddo
2325	enddo
2326	enddo
2327	enddo
2328	elseif (c_%npara.eq.5) then ! distribute exponents over 5 variables, knowing their sum
2329	do i1=no,0,-1
2330	do i2=no-i1,0,-1
2331	do i3=no-i1-i2,0,-1
2332	do i4=no-i1-i2-i3,0,-1
2333	do i5=no-i1-i2-i3-i4,0,-1
2334	if (i1+i2+i3+i4+i5==no) then
2335	j(1)=i1
2336	j(2)=i2
2337	j(3)=i3
2338	j(4)=i4
2339	j(5)=i5
2340	coef = y(i)%T.sub.j
2341	if (coef.ne.zero) then
2342	map_coor(1)=coef
2343	map_coor(2)=i
2344	map_coor(3)=c_%npara
2345	map_coor(4)=no
2346	map_coor(5)=j(1)
2347	map_coor(6)=j(2)
2348	map_coor(7)=j(3)
2349	map_coor(8)=j(4)
2350	map_coor(9)=j(5)
2351	map_coor(10) = 0
2352	call vector_to_table_curr("map_table ", 'coef ', map_coor(1), i_map_coor)
2353	call augment_count("map_table ")
2354	endif
2355	endif
2356	enddo
2357	enddo
2358	enddo
2359	enddo
2360	enddo
2361	elseif (c_%npara.eq.4) then ! distribute exponents over 4 variables, knowing their sum
2362	do i1=no,0,-1
2363	do i2=no-i1,0,-1
2364	do i3=no-i1-i2,0,-1
2365	do i4=no-i1-i2-i3,0,-1
2366	if (i1+i2+i3+i4==no) then
2367	j(1)=i1
2368	j(2)=i2
2369	j(3)=i3
2370	j(4)=i4
2371	coef = y(i)%T.sub.j
2372	if (coef.ne.zero) then
2373	map_coor(1)=coef
2374	map_coor(2)=i
2375	map_coor(3)=c_%npara
2376	map_coor(4)=no
2377	map_coor(5)=j(1)
2378	map_coor(6)=j(2)
2379	map_coor(7)=j(3)
2380	map_coor(8)=j(4)
2381	map_coor(9)=0
2382	map_coor(10)=0
2383	call vector_to_table_curr("map_table ", 'coef ', map_coor(1), i_map_coor)
2384	call augment_count("map_table ")
2385	endif
2386	endif
2387	enddo
2388	enddo
2389	enddo
2390	enddo
2391	else
2392	call fort_warn('ptc_normal ','map output expects 4,5 or 6 variables')
2393	endif
2394	enddo
2395	enddo
2396
2397
2398
2399
2400
2401	deallocate(j)
2402
2403
2404
2405	end subroutine makemaptable
2406
2407
2408	!_________________________________________________________________
2409
2410	subroutine killsavedmaps
2411	implicit none
2412	integer i,ii
2413
2414	if(.not. savemaps) return
2415
2416	if (.not. associated(maps)) then
2417	return
2418	endif
2419
2420	do i=lbound(maps,1),ubound(maps,1)
2421	do ii=1,6
2422	call kill(maps(i)%unimap(ii))
2423	enddo
2424	enddo
2425	deallocate(maps)
2426	nullify(maps)
2427
2428	end subroutine killsavedmaps
2429	!_________________________________________________________________
2430
2431
2432	SUBROUTINE ptc_read_errors()
2433	use twtrrfi
2434	use name_lenfi
2435	implicit none
2436	integer i,k,pos,nfac(maxmul),flag,string_from_table_row,double_from_table_row,l
2437	real(dp) d(2*maxmul),b(maxmul),a(maxmul),tilt,ab,bvk
2438	character(name_len) name,name2
2439	type(fibre),pointer :: p
2440	logical(lp) :: overwrite
2441	real(kind(1d0)) get_value
2442	character*4 :: mag_index1(10)=(/'k0l ','k1l ','k2l ','k3l ','k4l ','k5l ','k6l ','k7l ','k8l ','k9l '/)
2443	character*5 :: mag_index2(10)=(/'k0sl ','k1sl ','k2sl ','k3sl ','k4sl ','k5sl ','k6sl ','k7sl ','k8sl ','k9sl '/)
2444	character*5 :: mag_index3(11)=(/'k10l ','k11l ','k12l ','k13l ','k14l ','k15l ','k16l ','k17l ','k18l ','k19l ','k20l '/)
2445	character*6 :: mag_index4(11)=(/'k10sl ','k11sl ','k12sl ','k13sl ','k14sl ','k15sl ','k16sl ', &
2446	'k17sl ','k18sl ','k19sl ','k20sl '/)
2447
2448	overwrite = get_value('ptc_read_errors ','overwrite ').ne.0
2449	bvk=get_value('probe ','bv ')
2450
2451	nfac(1)=1
2452	do i=2,maxmul
2453	nfac(i)=nfac(i-1)*(i-1)
2454	enddo
2455
2456	flag = string_from_table_row('errors_read ', 'name ',1,name)
2457
2458	if(flag.ne.0) call aafail('fill_errors reports: ',' The >>> errors_read <<< table is empty ')
2459	i=0
2460
2461	p=>my_ring%start
2462	do while(.true.)
2463	i=i+1
2464	b(:)=zero
2465	a(:)=zero
2466	d(:)=zero
2467	name2=" "
2468	flag = string_from_table_row('errors_read ', 'name ',i,name2)
2469	if(flag.ne.0) goto 100
2470	do k=1,maxmul
2471	if(k<=10) then
2472	flag = double_from_table_row('errors_read ',mag_index1(k),i,d(2*k-1))
2473	flag = double_from_table_row('errors_read ',mag_index2(k),i,d(2*k))
2474	else
2475	flag = double_from_table_row('errors_read ',mag_index3(k-10),i,d(2*k-1))
2476	flag = double_from_table_row('errors_read ',mag_index4(k-10),i,d(2*k))
2477	endif
2478	enddo
2479	if(flag.ne.0) goto 100
2480	do k=1,maxmul
2481	b(k)=d(2*k-1)/nfac(k)
2482	a(k)=d(2*k)/nfac(k)
2483	enddo
2484	name=" "
2485	name(:len_trim(name2)-1)=name2(:len_trim(name2)-1)
2486	call context(name)
2487	call move_to(my_ring,p,name,pos)
2488	tilt=-p%mag%p%tiltd
2489	if(pos/=0) then
2490	if(p%mag%l/=zero) then
2491	do k=1,maxmul
2492	b(k)=b(k)/p%mag%l
2493	a(k)=a(k)/p%mag%l
2494	enddo
2495	endif
2496	do k=1,maxmul
2497	b(k)=bvk*b(k)
2498	a(k)=bvk*a(k)
2499	enddo
2500	if(tilt/=zero) then
2501	do k=1,maxmul
2502	ab=b(k)
2503	b(k)=b(k)cos(tiltk)+a(k)sin(tiltk)
2504	a(k)=-absin(tiltk)+a(k)cos(tiltk)
2505	enddo
2506	endif
2507	do k=NMAX,1,-1
2508	if(b(k)/=zero) then
2509	if(overwrite) then
2510	call add(p,k,0,b(k))
2511	else
2512	call add(p,k,1,b(k))
2513	endif
2514	endif
2515	if(a(k)/=zero) then
2516	if(overwrite) then
2517	call add(p,-k,0,a(k))
2518	else
2519	call add(p,-k,1,a(k))
2520	endif
2521	endif
2522	enddo
2523	else
2524	write(6,*) " name,pos, dir of dna ",name, p%mag%parent_fibre%dir
2525	endif
2526	enddo
2527	100 continue
2528	return
2529
2530	end SUBROUTINE ptc_read_errors
2531
2532	subroutine ptc_refresh_k()
2533	use twtrrfi
2534	use name_lenfi
2535	implicit none
2536	integer j,code,k,pos,nfac(maxmul)
2537	integer restart_sequ,advance_node
2538	type(fibre),pointer :: p
2539	real(dp) sk,sks,tilt,b(maxmul),a(maxmul),bvk
2540	real(kind(1d0)) :: get_value,node_value
2541	character(name_len) name
2542	logical(lp) :: overwrite
2543	!---------------------------------------------------------------
2544
2545	overwrite = get_value('ptc_refresh_k ','overwrite ').ne.0
2546	bvk=get_value('probe ','bv ')
2547
2548	nfac(1)=1
2549	do j=2,maxmul
2550	nfac(j)=nfac(j-1)*(j-1)
2551	enddo
2552
2553	j=restart_sequ()
2554	j=0
2555	p=>my_ring%start
2556	10 continue
2557	b(:)=zero
2558	a(:)=zero
2559
2560	code=node_value('mad8_type ')
2561	if(code.ne.5.and.code.ne.6) goto 100
2562	if(code.eq.5) then
2563	! quadrupole components code = 5
2564	k=2
2565	sk= node_value('k1 ')
2566	sks=node_value('k1s ')
2567	tilt=node_value('tilt ')
2568	b(k)=sk
2569	if (sks .ne. zero) then
2570	tilt = -atan2(sks, sk)/two + tilt
2571	b(k)=sqrt(sk2+sks2)/abs(sk)*sk !
2572	endif
2573	elseif(code.eq.6) then
2574	! sextupole components code = 6
2575	k=3
2576	sk= node_value('k2 ')
2577	sks=node_value('k2s ')
2578	tilt=node_value('tilt ')
2579	b(k)=sk
2580	if (sks .ne. zero) then
2581	tilt = -atan2(sks, sk)/three + tilt
2582	b(k)=sqrt(sk2+sks2)/abs(sk)*sk !
2583	endif
2584	endif
2585
2586	call element_name(name,name_len)
2587	call context(name)
2588	call move_to(my_ring,p,name,pos)
2589	if(pos/=0) then
2590	b(k)=b(k)/nfac(k)
2591	if(tilt/=zero) then
2592	a(k)=-b(k)sin(tiltk)
2593	b(k)=b(k)cos(tiltk)
2594	endif
2595	b(k)=bvk*b(k)
2596	a(k)=bvk*a(k)
2597	do j=1,maxmul
2598	if(overwrite) then
2599	call add(p,j,0,b(j))
2600	call add(p,-j,0,a(j))
2601	else
2602	call add(p,j,1,b(j))
2603	call add(p,-j,1,a(j))
2604	endif
2605	enddo
2606	else
2607	write(6,*) " name,pos, dir of dna ",name, p%mag%parent_fibre%dir
2608	endif
2609
2610	100 continue
2611
2612	if(advance_node().ne.0) goto 10
2613
2614	return
2615
2616	END subroutine ptc_refresh_k
2617
2618	subroutine getfk(fk)
2619	!returns FK factor for Beam-Beam effect
2620	implicit none
2621	real (dp) :: fk,dpp
2622	real (dp) :: gamma0,beta0,beta_dp,ptot,b_dir,arad,totch
2623	real (dp) :: q,q_prime
2624	integer :: b_dir_int
2625	real(kind(1d0)) :: get_value
2626	real(kind(1d0)) :: get_variable
2627	integer :: get_option
2628	REAL(KIND(1d0)) :: node_value !/returns value for parameter par of current element /
2629
2630	!---- Calculate momentum deviation and according changes
2631	! of the relativistic factor beta0
2632
2633	gamma0 = get_value('probe ','gamma ')
2634	arad=get_value('probe ', 'arad ')
2635	totch=node_value('charge ') * get_value('probe ', 'npart ')
2636
2637	dpp = get_variable('track_deltap ')
2638	q = get_value('probe ','charge ')
2639	q_prime = node_value('charge ')
2640
2641	beta0 = sqrt(one-one/gamma0**2)
2642	ptot = beta0gamma0(one+dpp)
2643	beta_dp = ptot / sqrt(one + ptot**2)
2644	b_dir_int = node_value('bbdir ')
2645	b_dir=dble(b_dir_int)
2646	b_dir = b_dir/sqrt(b_dir*b_dir + 1.0d-32)
2647
2648	fk = twoaradtotch/gamma0 /beta0/(one+dpp)/q* &
2649	(one-beta0beta_dpb_dir)/(beta_dp+0.5(b_dir-one)b_dir*beta0)
2650
2651	end subroutine getfk
2652
2653
2654
2655
2656	subroutine putbeambeam()
2657	!returns FK factor for Beam-Beam effect
2658	implicit none
2659	real (dp) :: fk, xma, yma, sigx, sigy,s
2660	integer :: elno
2661	logical(lp) :: found
2662	TYPE(INTEGRATION_NODE),POINTER :: node
2663
2664	real(kind(1d0)), external :: get_value
2665	real(kind(1d0)), external :: get_variable
2666	integer, external :: get_option
2667	REAL(KIND(1d0)), external :: node_value !/returns value for parameter par of current element /
2668
2669	s = get_value('ptc_putbeambeam ','global_s ')
2670	xma = get_value('ptc_putbeambeam ','xma ')
2671	yma = get_value('ptc_putbeambeam ','yma ')
2672	sigx = get_value('ptc_putbeambeam ','sigx ')
2673	sigy = get_value('ptc_putbeambeam ','sigy ')
2674
2675	print*, 'Input xma, yma, sigx, sigy, s'
2676	print*, 'Input', xma, yma, sigx, sigy, s
2677
2678	if(.not.associated(my_ring%t)) then
2679	CALL MAKE_node_LAYOUT(my_ring)
2680	endif
2681
2682	elno = 0
2683	call s_locate_beam_beam(my_ring,s,elno,node,found)
2684
2685	if (.not. found) then
2686	print*,"could not find a node for beam-beam"
2687	return
2688	endif
2689
2690	print*, 'Name of element in PTC: ', node%PARENT_FIBRE%mag%name
2691
2692	write(6,*) 'node%a:',node%a ! node entrance position
2693	write(6,*) 'node%ent(1,1:3):',node%ent(1,1:3) ! node entrance e_1 vector
2694	write(6,*) 'node%ent(2,1:3):',node%ent(2,1:3) ! node entrance e_2 vector
2695	write(6,*) 'node%ent(3,1:3):',node%ent(3,1:3) ! node entrance e_3 vector
2696	write(6,*) " s variable of node and following node "
2697	write(6,*) node%s(1),node%next%s(1)
2698
2699	!N.B. If nothing else is done, the beam-beam kick is placed at the entrance of the node.
2700	!The call FIND_PATCH(t%a,t%ent,o ,mid,D,ANG) needs to be invoked to place the beam-beam kick
2701
2702	end subroutine putbeambeam
2703
2704	END MODULE madx_ptc_module

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

Context Navigation

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

Download in other formats: