Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

madx_ptc_normal.f90 @ 478

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

Line
1	module madx_ptc_normal_module
2	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
3	! madx_ptc_normal module
4	! Frank Schmidt (CERN)
5	!
6	use madx_ptc_module
7	!_________________________________________________________________
8	implicit none
9
10	public
11
12
13
14	private double_from_normal_t1,display_table_results
15
16	contains
17	!____________________________________________________________________________________________
18
19
20	SUBROUTINE ptc_normal()
21	USE ptc_results
22	USE madx_ptc_intstate_module
23	implicit none
24	logical(lp) closed_orbit,normal,maptable
25	integer no,mynd2,npara,nda,icase,flag_index,why(9)
26	integer i,ii,j1,k,l,starti
27	integer n_rows,row,n_haml,n_gnfu,nres,mynres,n1,n2
28	integer,external :: select_ptc_idx, minimum_acceptable_order, &
29	string_from_table_row, double_from_table_row, double_to_table_row
30	real(dp) x(6),deltap0,deltap,dt
31	integer :: indexa(4)
32	integer :: row_haml(101)
33	integer :: index1(1000,2)
34	real(kind(1d0)) get_value,val_ptc,tmp
35	character(len = 5) name_var
36	type(real_8) y(6)
37	type(real_8) :: theAscript(6) ! used here to compute dispersion's derivatives
38	type(normalform) :: n_t2 ! normal from type 2 for hamiltonian terms
39
40	!------------------------------------------------------------------------------
41
42	if(universe.le.0.or.EXCEPTION.ne.0) then
43	call fort_warn('return from ptc_normal: ',' no universe created')
44	return
45	endif
46	if(index_mad.le.0.or.EXCEPTION.ne.0) then
47	call fort_warn('return from ptc_normal: ',' no layout created')
48	return
49	endif
50
51	! do the check early so it leaves immediately and does not leave allocated objects
52	no = get_value('ptc_normal ','no ')
53
54	if (no < minimum_acceptable_order()) then
55	print*, "minimum acceptable order: ", minimum_acceptable_order()
56	call seterrorflag(11,"ptc_normal ","Order of calculation is not sufficient to calculate required parameters")
57	call fort_warn('ptc_normal: ',&
58	'Order of calculation (parameter no) is not sufficient to calculate required parameters')
59	return
60	endif
61
62	call cleartables() !defined in madx_ptc_knobs
63
64	nda=0
65
66	icase = get_value('ptc_normal ','icase ')
67	deltap0 = get_value('ptc_normal ','deltap ')
68
69	!
70	deltap = zero
71	call my_state(icase,deltap,deltap0)
72	CALL UPDATE_STATES
73
74	x(:)=zero
75	if(mytime) then
76	call Convert_dp_to_dt (deltap, dt)
77	else
78	dt=deltap
79	endif
80	if(icase.eq.5) x(5)=dt
81	closed_orbit = get_value('ptc_normal ','closed_orbit ') .ne. 0
82	if(closed_orbit) then
83	! pass starting point for closed orbit search
84	x(1)=get_value('ptc_normal ','x ')
85	x(2)=get_value('ptc_normal ','px ')
86	x(3)=get_value('ptc_normal ','y ')
87	x(4)=get_value('ptc_normal ','py ')
88	x(6)=get_value('ptc_normal ','t ')
89	x(5)=x(5)+get_value('ptc_normal ','pt ')
90
91	call find_orbit(my_ring,x,1,default,c_1d_7)
92	CALL write_closed_orbit(icase,x)
93	endif
94
95
96	call init(default,no,nda,BERZ,mynd2,npara)
97
98
99	call alloc(y)
100	y=npara
101	Y=X
102
103	c_%watch_user=.true.
104	call track(my_ring,y,1,default)
105	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
106	call fort_warn('ptc_normal: ','DA got unstable')
107	call seterrorflag(10,"ptc_normal ","DA got unstable ");
108	return
109	endif
110	call PRODUCE_APERTURE_FLAG(flag_index)
111	if(flag_index/=0) then
112	call ANALYSE_APERTURE_FLAG(flag_index,why)
113	Write(6,*) "ptc_normal unstable (map production)-programs continues "
114	Write(6,*) why ! See produce aperture flag routine in sd_frame
115	CALL kill(y)
116	c_%watch_user=.false.
117	return
118	endif
119	c_%watch_user=.false.
120	!if (getdebug()>1)
121	print77=.false.
122	read77 =.false.
123
124	call daprint(y,18)
125
126	maptable = get_value('ptc_normal ','maptable ') .ne. 0
127	if(maptable) then
128	call makemaptable(y)
129	endif
130
131
132	normal = get_value('ptc_normal ','normal ') .ne. 0
133	if(normal) then
134
135	!------ Find the number of occurences of the attribute 'haml'
136
137	n_rows = select_ptc_idx()
138
139	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
140	!! Q (+chroma+anharmon) and DX
141	!here we do normal type 1
142	call alloc(n)
143	if (getdebug() > 0) print*,"Normal Form Type 1"
144	n=y
145
146	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
147	call fort_warn('ptc_normal: ','Fatal Error: DA in NormalForm got unstable')
148	stop
149	endif
150
151	if (getdebug() > 1) then
152	write(19,'(/a/)') 'Dispersion, First and Higher Orders'
153	call daprint(n%A1,19)
154	endif
155
156
157	if (getdebug() > 1) then
158	write(19,'(/a/)') 'Tunes, Chromaticities and Anharmonicities'
159	! call daprint(n%A_t,19)
160	! call daprint(n%A,19)
161	call daprint(n%dhdj,19) ! orig one
162	! call daprint(pbrh,19)
163	endif
164
165
166	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
167	!! HAML
168	!!here we do normal type 2
169	!! It has resonanses left in the normal form
170	!! which are defined with normal%M
171
172	call alloc(n_t2)
173
174	nres = 0
175	n_haml = 0
176	n_gnfu = 0
177
178	if (n_rows > 0) then
179	do row = 1,n_rows
180	name_var=" "
181	k = string_from_table_row('normal_results ', 'name ', row, name_var)
182
183
184	if (name_var(:4) .eq. 'gnfu') then
185	if (getdebug() > 2) print*,"ptc_normal: adding gnfu"
186	n_gnfu = n_gnfu + 1
187	endif
188	if (name_var(:4) .eq. 'haml') then
189	n_haml = n_haml + 1
190	row_haml(n_haml) = row
191	endif
192	enddo
193
194	if (n_haml > 0) then
195	call alloc(pbrh)
196	do j1 =1,n_haml
197	row = row_haml(j1)
198	k = double_from_table_row("normal_results ", "value ", row, doublenum)
199	mynres = int(doublenum)
200	row = row_haml(j1) - 3*mynres + 2
201	starti = 1
202	if (j1 .eq. 1) then
203	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
204	indexa(1) = int(doublenum)
205	k = double_from_table_row("normal_results ", "order2 ", row, doublenum)
206	indexa(2) = int(doublenum)
207	k = double_from_table_row("normal_results ", "order3 ", row, doublenum)
208	indexa(3) = int(doublenum)
209	k = double_from_table_row("normal_results ", "order4 ", row, doublenum)
210	indexa(4) = int(doublenum)
211	index1(1,1) = indexa(1) - indexa(2)
212	index1(1,2) = indexa(3) - indexa(4)
213	n_t2%m(1,1)= index1(1,1)
214	n_t2%m(2,1)= index1(1,2)
215	if(c_%nd2.eq.6) n_t2%m(3,1)= indexa(3)
216	nres = 1
217	starti = 2
218	endif
219	!============ nres is the number of resonances to be set
220
221	if (mynres .ge. starti) then
222	do i = starti,mynres
223	ii = row + 3*(i-1)
224	k = double_from_table_row("normal_results ", "order1 ", ii, doublenum)
225	indexa(1) = int(doublenum)
226	k = double_from_table_row("normal_results ", "order2 ", ii, doublenum)
227	indexa(2) = int(doublenum)
228	k = double_from_table_row("normal_results ", "order3 ", ii, doublenum)
229	indexa(3) = int(doublenum)
230	k = double_from_table_row("normal_results ", "order4 ", ii, doublenum)
231	indexa(4) = int(doublenum)
232	n1 = indexa(1) - indexa(2)
233	n2 = indexa(3) - indexa(4)
234	do l = 1,nres
235	if (n1 .eq. index1(l,1) .and. n2 .eq. index1(l,2)) goto 100
236	enddo
237	nres = nres + 1
238	index1(nres,1) = n1
239	index1(nres,2) = n2
240	n_t2%m(1,nres)= n1
241	n_t2%m(2,nres)= n2
242	if(c_%nd2.eq.6) n_t2%m(3,nres)= indexa(3)
243	100 continue
244	enddo
245	endif
246	enddo
247	n_t2%nres = nres
248	endif
249
250	endif
251	!------------------------------------------------------------------------
252
253	if (n_haml > 0) then
254
255	if (getdebug() > 0) print*,"Normal Form Type 2 (for Hamiltonian Terms)"
256	n_t2=y
257
258	if (( .not. check_stable ) .or. ( .not. c_%stable_da )) then
259	call fort_warn('ptc_normal: ','Fatal Error: DA in NormalForm got unstable')
260	stop
261	endif
262
263
264
265	pbrh = n_t2%normal%pb
266
267	endif
268
269	if (n_gnfu > 0) then
270	call alloc(pbrg)
271	pbrg = n%a%pb
272	endif
273
274	!------ get values and store them in the table 'normal_results' ---------
275	if (n_rows > 0) then
276	do row = 1,n_rows
277	name_var=" "
278	k = string_from_table_row("normal_results ", "name ", row, name_var)
279
280	if (name_var(:3) .eq. 'ham') then
281	val_ptc = double_from_normal_t2(name_var, row, icase) !!! HERE IS THE RETRIVAL
282	else
283	val_ptc = double_from_normal_t1(name_var, row, icase) !!! HERE IS THE RETRIVAL
284	endif
285
286	if (name_var(:4) .ne. 'haml' .and. name_var(:4) .ne. 'gnfu') then
287	k = double_to_table_row("normal_results ", "value ", row, val_ptc)
288	endif
289
290	enddo
291	endif
292
293	call alloc(theAscript)
294	theAscript = X+n%A_t;
295
296	call putusertable(1,'$end$ ',dt ,dt,y,theAscript)
297	call kill(theAscript)
298
299
300
301	if (n_gnfu > 0) call kill(pbrg)
302	if (n_haml > 0) call kill(pbrh)
303
304
305	call kill(n)
306	call kill(n_t2)
307
308	endif
309
310
311
312	CALL kill(y)
313
314	close(18)
315
316	if (getdebug() > 1) then
317	close(19)
318	endif
319
320	! f90flush is not portable, and useless...
321	! call f90flush(18,my_false)
322	! call f90flush(19,my_false)
323
324	END subroutine ptc_normal
325	!________________________________________________________________________________
326
327	FUNCTION double_from_normal_t1(name_var,row,icase)
328	USE ptc_results
329	implicit none
330	logical(lp) name_l
331	integer,intent(IN) :: row,icase
332	real(dp) double_from_normal_t1, d_val, d_val1, d_val2
333	integer ii,i1,i2,jj
334	integer j,k,ind(6)
335	integer double_from_table_row
336	character(len = 4) name_var
337	character(len = 2) name_var1
338	character(len = 3) name_var2
339
340	name_l = .false.
341	double_from_normal_t1 = zero
342
343	name_var1 = name_var
344	SELECT CASE (name_var1)
345	CASE ("dx")
346	ind(:)=0
347	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
348	ind(5) = int(doublenum)
349	if (ind(5) == 0) ind(5) = 1
350	ind(6) = 0
351	d_val = n%A1%V(1).sub.ind
352	CASE ('dy')
353	ind(:)=0
354	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
355	ind(5) = int(doublenum)
356	if (ind(5) == 0) ind(5) = 1
357	ind(6) = 0
358	d_val = n%A1%V(3).sub.ind
359	CASE ('q1')
360	ind(:)=0
361	d_val = n%dhdj%V(c_%nd+1).sub.ind
362	CASE ('q2')
363	ind(:)=0
364	d_val = n%dhdj%V(c_%nd+2).sub.ind
365	CASE DEFAULT
366	name_l = .true.
367	END SELECT
368	if (name_l) then
369	name_l = .false.
370	name_var2 = name_var
371	SELECT CASE (name_var2)
372	CASE ('dpx')
373	ind(:)=0
374	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
375	ind(5) = int(doublenum)
376	if (ind(5) == 0) ind(5) = 1
377	ind(6) = 0
378	d_val = n%A1%V(2).sub.ind
379	CASE ('dpy')
380	ind(:)=0
381	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
382	ind(5) = int(doublenum)
383	if (ind(5) == 0) ind(5) = 1
384	ind(6) = 0
385	d_val = n%A1%V(4).sub.ind
386	CASE ('dq1')
387	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
388	ind(:)=0
389	ind(5) = int(doublenum)
390	if (ind(5) == 0) ind(5) = 1
391	ind(6) = 0
392	d_val = n%dhdj%V(c_%nd+1).sub.ind ! LD: was V(3)
393	CASE ('dq2')
394	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
395	ind(:)=0
396	ind(5) = int(doublenum)
397	if (ind(5) == 0) ind(5) = 1
398	ind(6) = 0
399	d_val = n%dhdj%V(c_%nd+2).sub.ind ! LD: was V(4)
400	CASE DEFAULT
401	name_l = .true.
402	END SELECT
403	endif
404	if (name_l) then
405	SELECT CASE (name_var)
406	CASE ('anhx')
407	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
408	do j = 1,2
409	ind(j) = int(doublenum)
410	enddo
411	k = double_from_table_row("normal_results ", "order2 ", row, doublenum)
412	do j = 3,4
413	ind(j) = int(doublenum)
414	enddo
415	k = double_from_table_row("normal_results ", "order3 ", row, doublenum)
416	ind(5) = int(doublenum)
417	ind(6) = 0
418	d_val = n%dhdj%V(c_%nd+1).sub.ind ! LD: was V(3)
419	CASE ('anhy')
420	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
421	do j = 1,2
422	ind(j) = int(doublenum)
423	enddo
424	k = double_from_table_row("normal_results ", "order2 ", row, doublenum)
425	do j = 3,4
426	ind(j) = int(doublenum)
427	enddo
428	k = double_from_table_row("normal_results ", "order3 ", row, doublenum)
429	ind(5) = int(doublenum)
430	ind(6) = 0
431	d_val = n%dhdj%V(c_%nd+2).sub.ind ! LD: was V(4)
432	CASE ('eign')
433	ii=(icase/2)*2
434	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
435	i1 = int(doublenum)
436	if(i1.gt.ii) call aafail('return from double_from_normal_t1: ',' wrong # of eigenvectors')
437	jj=0
438	if(i1.eq.5) jj=6
439	if(i1.eq.6) i1=5
440	if(jj.eq.6) i1=jj
441	k = double_from_table_row("normal_results ", "order2 ", row, doublenum)
442	i2 = int(doublenum)
443	if(i2.gt.ii) call aafail('return from double_from_normal_t1: ',' eigenvectors too many components')
444	jj=0
445	if(i2.eq.5) jj=6
446	if(i2.eq.6) i2=5
447	if(jj.eq.6) i2=jj
448	ind(:)=0
449	ind(i2)=1
450	double_from_normal_t1 = n%A_t%V(i1).sub.ind
451	if(mytime.and.i2.eq.6) double_from_normal_t1 = -double_from_normal_t1
452	RETURN
453	CASE ('gnfc')
454	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
455	ind(1) = int(doublenum)
456	k = double_from_table_row("normal_results ", "order2 ", row, doublenum)
457	ind(2) = int(doublenum)
458	k = double_from_table_row("normal_results ", "order3 ", row, doublenum)
459	ind(3) = int(doublenum)
460	k = double_from_table_row("normal_results ", "order4 ", row, doublenum)
461	ind(4) = int(doublenum)
462	ind(5) = 0
463	ind(6) = 0
464	d_val = pbrg%cos%h.sub.ind
465	double_from_normal_t1 = d_val
466	RETURN
467	CASE ('gnfs')
468	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
469	ind(1) = int(doublenum)
470	k = double_from_table_row("normal_results ", "order2 ", row, doublenum)
471	ind(2) = int(doublenum)
472	k = double_from_table_row("normal_results ", "order3 ", row, doublenum)
473	ind(3) = int(doublenum)
474	k = double_from_table_row("normal_results ", "order4 ", row, doublenum)
475	ind(4) = int(doublenum)
476	ind(5) = 0
477	ind(6) = 0
478	d_val = pbrg%sin%h.sub.ind
479	double_from_normal_t1 = d_val
480	RETURN
481	CASE ('gnfa')
482	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
483	ind(1) = int(doublenum)
484	k = double_from_table_row("normal_results ", "order2 ", row, doublenum)
485	ind(2) = int(doublenum)
486	k = double_from_table_row("normal_results ", "order3 ", row, doublenum)
487	ind(3) = int(doublenum)
488	k = double_from_table_row("normal_results ", "order4 ", row, doublenum)
489	ind(4) = int(doublenum)
490	ind(5) = 0
491	ind(6) = 0
492	d_val1 = pbrg%cos%h.sub.ind
493	d_val2 = pbrg%sin%h.sub.ind
494	double_from_normal_t1 = SQRT(d_val12 + d_val22)
495	RETURN
496	CASE ('gnfu')
497	double_from_normal_t1 = zero ! should never arrive here
498	RETURN
499	CASE DEFAULT
500	print *,"--Error in the table normal_results-- Unknown input: ",name_var
501	END SELECT
502	endif
503	double_from_normal_t1 = d_val(factorial(ind(1))factorial(ind(3))*factorial(ind(5)))
504	END FUNCTION double_from_normal_t1
505	!________________________________________________
506	!Extraction of normal type 2 variables: hamiltonian and generating functions
507
508	FUNCTION double_from_normal_t2(name_var,row,icase)
509	USE ptc_results
510	implicit none
511	logical(lp) name_l
512	integer,intent(IN) :: row,icase
513	real(dp) double_from_normal_t2, d_val, d_val1, d_val2
514	integer ii,i1,i2,jj
515	integer j,k,ind(6)
516	integer double_from_table_row
517	character(len = 4) name_var
518	character(len = 2) name_var1
519	character(len = 3) name_var2
520
521	double_from_normal_t2 = zero
522
523	name_var1 = name_var
524
525	! if (LEN_TRIM(name_var))
526
527	SELECT CASE (name_var)
528	CASE ('hamc')
529	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
530	ind(1) = int(doublenum)
531	k = double_from_table_row("normal_results ", "order2 ", row, doublenum)
532	ind(2) = int(doublenum)
533	k = double_from_table_row("normal_results ", "order3 ", row, doublenum)
534	ind(3) = int(doublenum)
535	k = double_from_table_row("normal_results ", "order4 ", row, doublenum)
536	ind(4) = int(doublenum)
537	ind(5) = 0
538	ind(6) = 0
539	d_val = pbrh%cos%h.sub.ind
540	double_from_normal_t2 = d_val
541	RETURN
542	CASE ('hams')
543	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
544	ind(1) = int(doublenum)
545	k = double_from_table_row("normal_results ", "order2 ", row, doublenum)
546	ind(2) = int(doublenum)
547	k = double_from_table_row("normal_results ", "order3 ", row, doublenum)
548	ind(3) = int(doublenum)
549	k = double_from_table_row("normal_results ", "order4 ", row, doublenum)
550	ind(4) = int(doublenum)
551	ind(5) = 0
552	ind(6) = 0
553	d_val = pbrh%sin%h.sub.ind
554	double_from_normal_t2 = d_val
555	RETURN
556	CASE ('hama')
557	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
558	ind(1) = int(doublenum)
559	k = double_from_table_row("normal_results ", "order2 ", row, doublenum)
560	ind(2) = int(doublenum)
561	k = double_from_table_row("normal_results ", "order3 ", row, doublenum)
562	ind(3) = int(doublenum)
563	k = double_from_table_row("normal_results ", "order4 ", row, doublenum)
564	ind(4) = int(doublenum)
565	ind(5) = 0
566	ind(6) = 0
567	d_val1 = pbrh%cos%h.sub.ind
568	d_val2 = pbrh%sin%h.sub.ind
569	double_from_normal_t2 = SQRT(d_val12 + d_val22)
570	RETURN
571	CASE ('haml')
572	double_from_normal_t2 = zero
573	RETURN
574	CASE DEFAULT
575	print *,"--Error in the table normal_results-- Unknown input: ",name_var
576	END SELECT
577
578	double_from_normal_t2 = d_val(factorial(ind(1))factorial(ind(3))*factorial(ind(5)))
579	END FUNCTION double_from_normal_t2
580	!_________________________________________________________________
581
582	SUBROUTINE display_table_results()
583	implicit none
584	integer, external :: select_ptc_idx, string_from_table_row, double_from_table_row
585	character(len = 4) name_var
586	integer row,k
587	integer :: ord(3)
588
589	print *,"Variable name Order 1 order 2 order 3 Value "
590	do row = 1 , select_ptc_idx()
591	name_var=" "
592	k = string_from_table_row("normal_results ", "name ", row, name_var)
593	k = double_from_table_row("normal_results ", "order1 ", row, doublenum)
594	ord(1) = int(doublenum)
595	k = double_from_table_row("normal_results ", "order2 ", row, doublenum)
596	ord(2) = int(doublenum)
597	k = double_from_table_row("normal_results ", "order3 ", row, doublenum)
598	ord(3) = int(doublenum)
599	k = double_from_table_row("normal_results ", "value ", row, doublenum)
600	WRITE(*,100) name_var,ord(1),ord(2),ord(3),doublenum
601	enddo
602	100 FORMAT(3X,A4,14X,I1,8X,I1,8X,I1,5X,f25.18)
603	END SUBROUTINE display_table_results
604	END MODULE madx_ptc_normal_module

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

Context Navigation

source: PSPA/madxPSPA/src/madx_ptc_normal.f90 @ 478

Download in other formats: