Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sodd.f90 @ 445

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

Line
1	subroutine soddin(ierr)
2	!---------------------------------------------------------------------
3	!---------------------------------------------------------------------
4	! The name SODD stands for "Second Order Detuning and Distortion"
5	! ---- - - - -
6	!---------------------------------------------------------------------
7	!---------------------------------------------------------------------
8	!
9	! Programm SODD consists of 3 parts:
10	!
11	! 1.) Subroutine Detune calculates the detuning function terms in
12	! first and second order in the strength of the multipoles.
13	! 2.) Subroutine Distort1 calculates the distortion function and
14	! Hamiltonian terms in first order in the strength of the
15	! multipoles.
16	! 3.) Subroutine Distort2 calculates the distortion function and
17	! Hamiltonian terms in second order in the strength of the
18	! multipoles.
19	!
20	!
21	! Theory: Analytical Calculation of Smear and Tune Shift SSC-232
22	! Authors: J. Bengtsson and J. Irwin
23	! Date: February 1990
24	!
25	! Program Author: Frank Schmidt, CERN SL-Group
26	! Date: November 1998 - January 1999
27	! Maximum Order: parameter mmul=11 (22-Pole)
28	! Comments to: Frank.Schmidt@cern.ch
29	!
30	!---------------------------------------------------------------------
31	!
32	! Requirements:
33	! -------------
34	!
35	! The program expects a file "fc.34" which contains 8 columns:
36	! - Position of Multipole [m]
37	! - Name of Multipole
38	! - Multipole Type (>0 erect Multipole; <0 skew Multipole);
39	! eg "3" is an erect sextupole, "-5" is a skew decapole
40	! - Single particle strength `a la SIXTRACK
41	! - Horizontal Betafunction
42	! - Vertical Betafunction
43	! - Horizontal Phase-advance
44	! - Vertical Phase-advance
45	!
46	! The last line in this file is special: All entries are given at
47	! the end of the accelerator with the name "END" and the artificial
48	! multipole type "100".
49	!
50	! This file can be automatically produced with SIXTRACK when the
51	! linear optics parameter are calculated with the "LINE" block
52	! (see web locaction: http://wwwslap.cern.ch/~frs/). It
53	! is also produced directly from MAD with the DOOM program by
54	! invoking the MAD-SIXTRACK convertor using the `special' flag (see
55	! web locaction: http://wwwslap.cern.ch/~hansg/doom/doom_six.html
56	! for details).
57	!
58	! The program needs roughly 50Mbytes in particular for the second
59	! order distortion function due to pairs of multipoles up to 11th
60	! order.
61	!
62	!---------------------------------------------------------------------
63	!
64	! Program Manual:
65	! ---------------
66	!
67	! a.) Program selection:
68	! Specify the program:
69	! - iprog=1,3,5,7 => Run Program Detune
70	! - iprog=2,3,6,7 => Run Program Distort1
71	! - iprog=4,5,6,7 => Run Program Distort2
72	! b.) Order selection:
73	! Specify low (n1) and high (n2) limit of orders to be studied
74	! erect and skew elements are denoted with positive and negative
75	! values respectively; eg.:
76	! n1=-4, n2=10 => the orders (-4,-3,-2 and 3-10) will be treated
77	! c.) Position Range:
78	! Specify a position range (etl1,etl2) in [m], the multipoles
79	! located between this range will be used for the calculation
80	! d.) Print_out Switch:
81	! Output files meant for spreadsheets, i.e. no comments. Human
82	! readable output is found in fort.6
83	! iu_on = 0 => No printout
84	! iu_on = 1 => Output at the end of the position range
85	! iu_on = 2 => At each Multipole in the position range
86	! In Detail:
87	!
88	! Program Detune
89	!
90	! Multipole
91	! Strength
92	! Unit iu_on Order Contents in each column
93	! ----------------------------------------------------------------
94	! 70 1 1 multipole order
95	! (hor.,ver.) plane => (1,2)
96	! hor. or ver. detuning
97	! order of horizontal Emittance
98	! order of vertical Emittance
99	! ----------------------------------------------------------------
100	! 71 2 1 multipole order
101	! (hor.,ver.) plane => (1,2)
102	! hor. or ver. detuning
103	! order of horizontal Emittance
104	! order of vertical Emittance
105	! ----------------------------------------------------------------
106	! 72 1 2 first multipole order
107	! second multipole order
108	! horizontal detuning
109	! order of horizontal Emittance
110	! order of vertical Emittance
111	! ----------------------------------------------------------------
112	! 73 1 2 first multipole order
113	! second multipole order
114	! vertical detuning
115	! order of horizontal Emittance
116	! order of vertical Emittance
117	! ---------------------------------------------------------------
118	!
119	! Program Distort1
120	!
121	! Multipole
122	! Strength
123	! Unit iu_on Order Contents in each column
124	! ----------------------------------------------------------------
125	! 74 1 1 multipole order
126	! cosine part of distortion
127	! sine part of distortion
128	! amplitude of distortion
129	! j
130	! k
131	! l
132	! m
133	! ----------------------------------------------------------------
134	! 75 1 1 multipole order
135	! cosine part of Hamiltonian
136	! sine part of Hamiltonian
137	! amplitude of Hamiltonian
138	! j
139	! k
140	! l
141	! m
142	! ----------------------------------------------------------------
143	! 76 2 1 multipole order
144	! appearance number in position range
145	! number of resonance
146	! position
147	! cosine part of Hamiltonian
148	! sine part of Hamiltonian
149	! amplitude of Hamiltonian
150	! j
151	! k
152	! l
153	! m
154	! ----------------------------------------------------------------
155	! 77 2 1 multipole order
156	! appearance number in position range
157	! number of resonance
158	! position
159	! cosine part of distortion
160	! sine part of distortion
161	! amplitude of distortion
162	! j
163	! k
164	! l
165	! m
166	! ----------------------------------------------------------------
167	!
168	! Program Distort2
169	!
170	! Multipole
171	! Strength
172	! Unit iu_on Order Contents in each column
173	! ----------------------------------------------------------------
174	! 78 1 2 first multipole order
175	! second multipole order
176	! cosine part of distortion
177	! sine part of distortion
178	! amplitude of distortion
179	! j
180	! k
181	! l
182	! m
183	! ----------------------------------------------------------------
184	! 79 1 2 first multipole order
185	! second multipole order
186	! cosine part of Hamiltonian
187	! sine part of Hamiltonian
188	! amplitude of Hamiltonian
189	! j
190	! k
191	! l
192	! m
193	! ----------------------------------------------------------------
194	!
195	!---------------------------------------------------------------------
196
197	implicit none
198
199	integer ierr
200
201	integer ier,iprog,mh,mmul,mmul2,mmult,mmultf,mmultw,mmultx, &
202	&n3,nblz
203	double precision etl20,pieni
204	integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, &
205	&iu_on,n1,n2
206	double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, &
207	&fac2,fac4,facd2,fact,fact0,factb,four,ham,one, &
208	&phi,pi,pi2,pihi,qx,qy,sbeta,sign,two,zero
209	real tlim,time0,time1,time2,time3,time4
210	real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
211	integer j70,j71,j72,j73,j74,j75,j76,j77,j78,j79
212	integer table_size_70,table_size_71,table_size_72,table_size_73, &
213	&table_size_74,table_size_75,table_size_76,table_size_77, &
214	&table_size_78,table_size_79
215	character*16 strn
216	character*18 comment
217	parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100)
218	parameter(mmul2=mmul+1,mmult=8mmul*2)
219	parameter(mmultw=2mmul,mmultx=mmult/7,mmultf=2mmultx*mmult)
220
221	!--- szcompar is the size of the arrays
222	!--- returned by the routine comm_para
223	integer szcompar
224	parameter (szcompar = 100)
225
226	!--- szchara is the size of the character strings char_a
227	!--- returned by the routine comm_para
228	integer szchara
229	parameter (szchara = 400)
230
231	common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul)
232	common/c1/ tlim,time0,time1,time2,time3,time4
233	common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
234	common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz)
235	common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz)
236	common/c4q/ phi(2,-mmul:mmul,nblz)
237	common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul)
238	common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2
239	common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult)
240	common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx)
241	common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul)
242	common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf)
243	common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul)
244	common/c12/ ihv(mmultx),iplane(mmultx,2)
245	common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh)
246	common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh)
247	common/indeces/j70,j71,j72,j73,j74,j75,j76,j77,j78,j79
248	common/table_sizes/table_size_70,table_size_71,table_size_72, &
249	&table_size_73,table_size_74,table_size_75, &
250	&table_size_76,table_size_77,table_size_78, &
251	&table_size_79
252
253	integer detun,disto1,disto2,prend,prblup
254	integer nint, ndble, k, int_arr(szcompar), char_l(szcompar)
255	double precision d_arr(szcompar)
256	character * (szchara) char_a
257
258	!---------------------------------------------------------------------
259
260	ierr = 0
261	j70 = 0
262	j71 = 0
263	j72 = 0
264	j73 = 0
265	j74 = 0
266	j75 = 0
267	j76 = 0
268	j77 = 0
269	j78 = 0
270	j79 = 0
271	table_size_70 = 2
272	table_size_71 = 2
273	table_size_72 = 2
274	table_size_73 = 2
275	table_size_74 = 2
276	table_size_75 = 2
277	table_size_76 = 2
278	table_size_77 = 2
279	table_size_78 = 2
280	table_size_79 = 2
281
282	!--- get detun
283
284	detun = 0
285	call comm_para('detune ', nint, ndble, k, int_arr, d_arr, &
286	&char_a, char_l)
287	if (nint .gt. 0) detun = int_arr(1)
288
289	!--- get disto1
290
291	disto1 = 0
292	call comm_para('distort1 ', nint, ndble, k, int_arr, d_arr, &
293	&char_a, char_l)
294	if (nint .gt. 0) disto1 = int_arr(1)
295
296	!--- get disto2
297
298	disto2 = 0
299	call comm_para('distort2 ', nint, ndble, k, int_arr, d_arr, &
300	&char_a, char_l)
301	if (nint .gt. 0) disto2 = int_arr(1)
302
303	iprog = detun + disto12 + disto24
304	if(iprog.le.0.or.iprog.gt.7) call prror(0,14,iprog)
305	if(iprog.eq.1.or.iprog.eq.3.or.iprog.eq.5.or.iprog.eq.7) &
306	&write(6,*) ' Program <DETUNE> will be executed '
307	if(iprog.eq.2.or.iprog.eq.3.or.iprog.eq.6.or.iprog.eq.7) &
308	&write(6,*) ' Program <DISTORT1> will be executed '
309	if(iprog.eq.4.or.iprog.eq.5.or.iprog.eq.6.or.iprog.eq.7) &
310	&write(6,*) ' Program <DISTORT2> will be executed '
311
312
313	char_a = ' '
314	call comm_para('multipole_order_range ', nint, ndble, k, int_arr, &
315	&d_arr,char_a, char_l)
316	n1 = int_arr(1)
317	n2 = int_arr(2)
318	if(n1.gt.n2) then
319	n3=n1
320	n1=n2
321	n2=n3
322	endif
323	if(abs(n1).gt.mmul.or.abs(n2).gt.mmul) call prror(0,15,mmul)
324
325	!--- get etl1 & etl2 (start position & stop position)
326
327	etl1 = 0.
328	etl2 = 0.
329	call comm_para('start_stop ', nint, ndble, k, int_arr, d_arr, &
330	&char_a, char_l)
331	if (ndble .gt. 0) then
332	etl1 = d_arr(1)
333	etl2 = d_arr(2)
334	endif
335
336	print *,"Start and End Position in [m]: ", etl1,etl2
337
338	if(etl1.lt.0) etl1=-etl1
339	if(etl2.lt.0) etl2=-etl2
340	if(etl2.lt.etl1) then
341	etl20=etl2
342	etl1=etl2
343	etl2=etl20
344	endif
345
346	!--- get no print
347
348	iu_on = 0
349	call comm_para('noprint ', nint, ndble, k, int_arr, d_arr, &
350	&char_a, char_l)
351
352	if (int_arr(1) .eq. 0) then
353
354	!--- get prend
355
356	prend = 0
357	call comm_para('print_at_end ', nint, ndble, k, int_arr, d_arr, &
358	&char_a, char_l)
359	if (nint .gt. 0) prend = int_arr(1)
360
361	!--- get prblup
362
363	prblup = 0
364	call comm_para('print_all ', nint, ndble, k, int_arr, d_arr, &
365	&char_a, char_l)
366	if (nint .gt. 0) prblup = int_arr(1)
367
368	iu_on = prend + prblup*2
369	endif
370
371	if(iu_on.lt.0.or.iu_on.gt.3) iu_on=0
372	if(iu_on.eq.0) write(6,*) ' No print_out requested '
373	if(iu_on.eq.1.or.iu_on.eq.3) then
374	write(6,*) ' Print_out at end of structure '
375	if(detun.eq.1) then
376	open(70,file='detune_1_end',form='formatted',status='unknown',&
377	&iostat=ier)
378	if(ier.ne.0) call prror(0,16,70)
379	write(70,10070)
380	open(72,file='detune_2_hor',form='formatted',status='unknown',&
381	&iostat=ier)
382	if(ier.ne.0) call prror(0,16,72)
383	write(72,10072)
384	open(73,file='detune_2_ver',form='formatted', &
385	&status='unknown',iostat=ier)
386	if(ier.ne.0) call prror(0,16,73)
387	write(73,10073)
388	endif
389	if(disto1.eq.1) then
390	open(74,file='distort_1_f_end',form='formatted', &
391	&status='unknown',iostat=ier)
392	if(ier.ne.0) call prror(0,16,74)
393	write(74,10074)
394	open(75,file='distort_1_h_end',form='formatted', &
395	&status='unknown',iostat=ier)
396	if(ier.ne.0) call prror(0,16,75)
397	write(75,10074)
398	endif
399	if(disto2.eq.1) then
400	open(78,file='distort_2_f_end',form='formatted', &
401	&status='unknown',iostat=ier)
402	if(ier.ne.0) call prror(0,16,78)
403	write(78,10078)
404	open(79,file='distort_2_h_end',form='formatted', &
405	&status='unknown',iostat=ier)
406	if(ier.ne.0) call prror(0,16,79)
407	write(79,10078)
408	endif
409	endif
410	if(iu_on.eq.2.or.iu_on.eq.3) then
411	write(6,*) ' Print_out at each multipole '
412	if(detun.eq.1) then
413	open(71,file='detune_1_all',form='formatted', &
414	&status='unknown',iostat=ier)
415	if(ier.ne.0) call prror(0,16,71)
416	write(71,10070)
417	endif
418	if(disto1.eq.1) then
419	open(76,file='distort_1_f_all',form='formatted', &
420	&status='unknown',iostat=ier)
421	if(ier.ne.0) call prror(0,16,76)
422	write(76,10076)
423	open(77,file='distort_1_h_all',form='formatted', &
424	&status='unknown',iostat=ier)
425	if(ier.ne.0) call prror(0,16,77)
426	write(77,10076)
427	endif
428	endif
429	open(34,file='fc.34',form='formatted',status='unknown', &
430	&recl=8192,iostat=ier)
431	if(ier.ne.0) call prror(0,16,34)
432	!-----------------------------------------------------------------------
433	! Initialisation
434	!-----------------------------------------------------------------------
435	call init
436	call timest(tlim)
437	call timex(tim1)
438	!-----------------------------------------------------------------------
439	! Read Data
440	!-----------------------------------------------------------------------
441	call readdat
442	!-----------------------------------------------------------------------
443
444	call timex(tim2)
445	write(6,10000) tim2-tim1
446	if(iprog.eq.1.or.iprog.eq.3.or.iprog.eq.5.or.iprog.eq.7) &
447	&call detune
448	call timex(tim2)
449	if(iprog.eq.2.or.iprog.eq.3.or.iprog.eq.6.or.iprog.eq.7) &
450	&call distort1
451	call timex(tim2)
452	if(iprog.eq.4.or.iprog.eq.5.or.iprog.eq.6.or.iprog.eq.7) &
453	&call distort2
454	call timex(tim6)
455	write(6,10010) tim6-tim1
456
457	continue
458	if(detun.eq.1) then
459	close(70,status='keep')
460	close(71,status='keep')
461	close(72,status='keep')
462	close(73,status='keep')
463	endif
464	if(disto1.eq.1) then
465	close(74,status='keep')
466	close(75,status='keep')
467	close(76,status='keep')
468	close(77,status='keep')
469	endif
470	if(disto2.eq.1) then
471	close(78,status='keep')
472	close(79,status='keep')
473	endif
474
475	10000 format(/80('-')//'Reading Data took ',f10.3,' second(s)', &
476	&' of Execution Time'//80('-')//)
477	10010 format(/80('-')//'Program <SODD> used a total of ',f10.3, &
478	&' second(s) of Execution Time'//80('-')//)
479	10070 format('MulOrd ',' Plane',4x,'Hor/Vert detuning',3x,'Hinv Vinv')
480	10072 format('MulOrd1 ','MulOrd2 ',3x,'Horizontal detuning', &
481	&2x,'Hinv Vinv')
482	10073 format('MulOrd1 ','MulOrd2 ',3x,'Vertical detuning ', &
483	&2x,'Hinv Vinv')
484	10074 format('MulOrd',11x,'Cosine',17x,'Sine',19x,'Amplitude', &
485	&8x,'j',3x,'k',3x,'l',3x,'m',1x)
486	10076 format('MulOrd',3x,'Loc',2x,'Res',6x,'Position[m]',14x,'Cosine', &
487	&17x,'Sine',19x,'Amplitude', &
488	&8x,'j',3x,'k',3x,'l',3x,'m',1x)
489	10078 format('MulOrd1',' MulOrd2',10x,'Cosine',17x,'Sine',19x, &
490	&'Amplitude',8x,'j',3x,'k',3x,'l',3x,'m',1x)
491	return
492	end subroutine soddin
493
494	subroutine detune
495	!-----------------------------------------------------------------------
496	!--Program Detune
497	!-----------------------------------------------------------------------
498	implicit none
499	integer i,ia,iaa,iah,iamin,ib,ic,id,ih1,ii,imu,imuty,j,jj,jstep,k,&
500	&k1,k2,k3,l,l0,l1,l2,l3,m,mh,mmul,mmul2,mmult,mmultf,mmultw, &
501	&mmultx,nblz
502	double precision betdx,betdxi,fac,facd0,facd1,phix,phiy,pieni, &
503	&sigde,sinar,tij
504	integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, &
505	&iu_on,n1,n2
506	double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, &
507	&fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, &
508	&qy,sbeta,sign,two,zero
509	real tlim,time0,time1,time2,time3,time4
510	real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
511	character*16 strn
512	character*18 comment
513	parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100)
514	parameter(mmul2=mmul+1,mmult=8mmul*2)
515	parameter(mmultw=2mmul,mmultx=mmult/7,mmultf=2mmultx*mmult)
516	common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul)
517	common/c1/ tlim,time0,time1,time2,time3,time4
518	common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
519	common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz)
520	common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz)
521	common/c4q/ phi(2,-mmul:mmul,nblz)
522	common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul)
523	common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2
524	common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult)
525	common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx)
526	common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul)
527	common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf)
528	common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul)
529	common/c12/ ihv(mmultx),iplane(mmultx,2)
530	common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh)
531	common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh)
532	dimension tij(0:mmul2,-mmul2:mmul2)
533	dimension betdxi(mmultx),betdx(mmultx)
534	dimension imu(2)
535	!---------------------------------------------------------------------
536	iamin=2
537	jstep=2
538	!-----------------------------------------------------------------------
539	! Read Data
540	!-----------------------------------------------------------------------
541	do i=0,mmul2
542	do j=-mmul2,mmul2
543	tij(i,j)=zero
544	enddo
545	enddo
546	!---------------------------------------------------------------------
547	!--
548	!--Detuning First Order in Multipole Strength
549	!--
550	!---------------------------------------------------------------------
551	do ia=2,mmul,jstep
552	imu(1)=ia
553	imu(2)=0
554	iah=ia/2
555	if(ityc(ia).gt.0) then
556	write(6,10000) 'Structure has: ',ityc(ia),' single ', &
557	&comment(ia)
558	ih1=ia/2-1
559	do i=1,ityc(ia)
560	fac=-two/dble(ia)*bstr(ia,i)
561	sigde=-one
562	do j=1,iah
563	jj=j-1
564	sigde=-sigde
565	det1(1,j)=sigdefaccosav(jj)cosav(iah-jj)dble(iah-jj)* &
566	&fact(ia,1+jj2)beta(1,ia,i)*(iah-jj)beta(2,ia,i)**jj
567	det(1,ia,iah-j,jj)=det(1,ia,iah-j,jj)+det1(1,j)
568	if(j.eq.iah) then
569	sigde=-sigde
570	det1(2,iah)=sigdefacdble(iah)cosav(iah) &
571	&beta(2,ia,i)**iah
572	det(2,ia,0,iah-1)=det(2,ia,0,iah-1)+det1(2,iah)
573	endif
574	enddo
575	if(iu_on.gt.0) call detwri(0,imu,ih1)
576	enddo
577	!---------------------------------------------------------------------
578	!--
579	!--Printing First Order Terms
580	!--
581	!---------------------------------------------------------------------
582	call detwri(1,imu,ih1)
583	endif
584	enddo
585	call timex(tim3)
586	write(6,10030) tim3-tim2
587	!---------------------------------------------------------------------
588	!--
589	!--Calculation of Detuning Poisson Brackets
590	!--
591	!-----Needed for Second Order Terms
592	!--
593	!---------------------------------------------------------------------
594	do imuty=1,2
595	call timex(tim3)
596	!---------------------------------------------------------------------
597	!--
598	!--Detuning Second Order in Multipole Strength
599	!--
600	!---------------------------------------------------------------------
601	do ia=iamin,mmul
602	do ib=ia,mmul,jstep
603	if(imuty.eq.1) ic=ia
604	if(imuty.eq.1) id=ib
605	if(imuty.eq.2) ic=-ia
606	if(imuty.eq.2) id=-ib
607	if(ityc(ic).gt.0.and.ityc(id).gt.0) then
608	imu(1)=ic
609	imu(2)=id
610	do k=1,mmultx
611	ihv(k)=0
612	iplane(k,1)=0
613	iplane(k,2)=0
614	betexp(k,1)=zero
615	betexp(k,2)=zero
616	betexp(k,3)=zero
617	betexp(k,4)=zero
618	fac2(k)=zero
619	do m=1,mmult
620	itij(k,m,1)=0
621	itij(k,m,2)=0
622	itij(k,m,3)=0
623	fac4(k,m)=zero
624	enddo
625	enddo
626	call caldet2(ia,ib,imuty)
627	do ii=1,ityc(ic)
628	do k=1,icc
629	betdxi(k)=(beta(1,ic,ii)*betexp(k,1)) &
630	&(beta(2,ic,ii)**betexp(k,2))
631	enddo
632	do jj=1,ityc(id)
633	if(ii.eq.1.and.jj.eq.1) then
634	if(ia.ne.ib) then
635	write(6,10050) 'Pairs of: ',ityc(ic),comment(ic), &
636	&' and: ',ityc(id),comment(id)
637	else
638	write(6,10060) 'Quadratic Effect of: ', &
639	&ityc(ic),comment(ic)
640	endif
641	endif
642	!--Phase Differences
643	phix=abs(phi(1,id,jj)-phi(1,ic,ii))-qx
644	phiy=abs(phi(2,id,jj)-phi(2,ic,ii))-qy
645	!--Phase Factor
646	if(imuty.eq.1) iaa=0
647	if(imuty.eq.2) iaa=1
648	do i=ia+iaa,0,-2
649	do j=i-ia,ia-i,2
650	if(i.ne.0.or.j.gt.0) then
651	sinar=sin(dble(i)qx+dble(j)qy)
652	if(abs(sinar).gt.pieni) then
653	tij(i,j)=cos(dble(i)phix+dble(j)phiy)/sinar
654	else
655	write(6,*) ' Warning: Results for detuning ',&
656	&'probably wrong; sinar= ',sinar
657	tij(i,j)=zero
658	endif
659	endif
660	enddo
661	enddo
662	!--Detuning Calculation
663	facd0=bstr(ic,ii)*bstr(id,jj)
664	do k=1,icc
665	k1=iplane(k,1)
666	k2=iplane(k,2)
667	k3=ihv(k)
668	if(k3.eq.1.or.(k3.eq.2.and.k1.eq.0)) then
669	betdx(k)=betdxi(k)* &
670	&(beta(1,id,jj)*betexp(k,3)) &
671	&(beta(2,id,jj)**betexp(k,4))
672	facd1=zero
673	l0=icd(k)
674	do l=1,l0
675	l1=itij(k,l,1)
676	l2=itij(k,l,2)
677	l3=itij(k,l,3)
678	facd1=facd1+fac4(k,l)*(tij(l2,l3)+ &
679	&dble(l1)*tij(l2,-l3))
680	enddo
681	det(k3,ic,k1,k2)=det(k3,ic,k1,k2)+facd0* &
682	&fac2(k)facd1betdx(k)
683	endif
684	enddo
685	do i=ia+iaa,0,-2
686	do j=i-ia,ia-i,2
687	if(i.ne.0.or.j.gt.0) then
688	tij(i,j)=zero
689	endif
690	enddo
691	enddo
692	enddo
693	enddo
694	!---------------------------------------------------------------------
695	!--
696	!--Printing Second Order Terms
697	!--
698	!---------------------------------------------------------------------
699	ih1=(ia+ib)/2-2
700	call detwri(2,imu,ih1)
701	endif
702	enddo
703	enddo
704	call timex(tim4)
705	write(6,10090) tim4-tim3
706	enddo
707	call timex(tim5)
708	write(6,10100) tim5-tim2
709	!---------------------------------------------------------------------
710	10000 format(//80('-')//5x,a,i6,2a)
711	10030 format(/'First Order Detuning Calculation took ', &
712	&f10.3,' second(s)',' of Execution Time'//80('-')//)
713	10050 format(//80('-')//10x,a,i6,2a,i6,a)
714	10060 format(//80('-')//10x,a,i6,a)
715	10090 format(/80('-')/'Second Order Detuning Calculation took ', &
716	&f10.3,' second(s)',' of Execution Time'//80('-')//)
717	10100 format(/80('-')//'Program <DETUNE> used ',f10.3, &
718	&' second(s) of Execution Time'//80('-')//)
719	return
720	end subroutine detune
721	subroutine distort1
722	!-----------------------------------------------------------------------
723	!--Program Distort1 (first order)
724	!-----------------------------------------------------------------------
725	implicit none
726	integer i,ia,iadd,ib,ic,id,id1,id11,id110,id12,id120,idiff1, &
727	&idiff2,ie,if,ii,iia,iia2,ivar,j,mh,mmul,mmul2,mmult,mmultf,mmultw,&
728	&mmultx,nblz
729	double precision distc,dists,dstc,dsts,facc,facs,fact1,factb1, &
730	&factb2,fpre,hamc,hams,hmc,hms,pieni,signii,sinar,tc,ts
731	integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, &
732	&iu_on,n1,n2
733	integer j70,j71,j72,j73,j74,j75,j76,j77,j78,j79
734	integer table_size_70,table_size_71,table_size_72,table_size_73, &
735	&table_size_74,table_size_75,table_size_76,table_size_77, &
736	&table_size_78,table_size_79
737	integer data_size
738	double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, &
739	&fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, &
740	&qy,sbeta,sign,two,zero,amp
741	real tlim,time0,time1,time2,time3,time4
742	real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
743	character*16 strn,table_name
744	character*18 comment
745	parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100)
746	parameter(mmul2=mmul+1,mmult=8mmul*2)
747	parameter(mmultw=2mmul,mmultx=mmult/7,mmultf=2mmultx*mmult)
748	common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul)
749	common/c1/ tlim,time0,time1,time2,time3,time4
750	common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
751	common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz)
752	common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz)
753	common/c4q/ phi(2,-mmul:mmul,nblz)
754	common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul)
755	common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2
756	common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult)
757	common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx)
758	common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul)
759	common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf)
760	common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul)
761	common/c12/ ihv(mmultx),iplane(mmultx,2)
762	common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh)
763	common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh)
764	common/indeces/j70,j71,j72,j73,j74,j75,j76,j77,j78,j79
765	common/table_sizes/table_size_70,table_size_71,table_size_72, &
766	&table_size_73,table_size_74,table_size_75, &
767	&table_size_76,table_size_77,table_size_78, &
768	&table_size_79
769	integer int_to_write(11)
770	double precision double_to_write(11)
771
772	dimension distc(-mmul:mmul,mh),dists(-mmul:mmul,mh)
773	dimension dstc(mh),dsts(mh),ivar(-mmul:mmul),fact1(mh)
774	dimension hamc(-mmul:mmul,mh),hams(-mmul:mmul,mh)
775	!---------------------------------------------------------------------
776	do i=1,mh
777	fact1(i)=zero
778	do j=-mmul,mmul
779	ivar(j)=0
780	ifacta(1,j,i)=0
781	ifacta(2,j,i)=0
782	ifact4(1,j,i)=0
783	ifact4(2,j,i)=0
784	ifact4(3,j,i)=0
785	ifact4(4,j,i)=0
786	factb(1,j,i)=zero
787	factb(2,j,i)=zero
788	fact0(j,i)=zero
789	distc(j,i)=zero
790	dists(j,i)=zero
791	hamc(j,i)=zero
792	hams(j,i)=zero
793	enddo
794	enddo
795	iadd=2
796	idiff1=1
797	idiff2=2
798	ivar(1)=1
799	ivar(-1)=1
800	!--Number of Resonances
801	do i=2,mmul
802	ivar(i)=iadd
803	ivar(-i)=iadd
804	if(mod(i,2).eq.0) then
805	idiff1=idiff1+idiff2
806	idiff2=idiff2+1
807	endif
808	iadd=iadd+idiff1
809	enddo
810	!--Calculation of the various Factors for Multipoles -mmul -> +mmul
811	do ia=-mmul,mmul
812	if(ityc(ia).gt.0) then
813	iia=iabs(ia)
814	fpre=one/dble((iia2*(iia+1)))
815	ib=0
816	signii=-one
817	if(ia.lt.0) then
818	iia2=iia-1
819	else
820	iia2=iia
821	endif
822	do ic=iia2,0,-2
823	signii=signii*(-one)
824	ib=ib+1
825	id=iia-ic
826	id1=id+1
827	ifacta(1,ia,ib)=ic
828	ifacta(2,ia,ib)=id
829	ifact4(1,ia,ib)=ic
830	ifact4(2,ia,ib)=0
831	ifact4(3,ia,ib)=id
832	ifact4(4,ia,ib)=0
833	factb1=abs(ifacta(1,ia,ib))/two
834	factb2=abs(ifacta(2,ia,ib))/two
835	factb(1,ia,ib)=factb1
836	factb(2,ia,ib)=factb2
837	fact0(ia,ib)=signiifprefact(iia,id1)
838	if(id.gt.0.and.ic.gt.0) then
839	ib=ib+1
840	ifacta(1,ia,ib)=ic
841	ifacta(2,ia,ib)=-id
842	ifact4(1,ia,ib)=ic
843	ifact4(2,ia,ib)=0
844	ifact4(3,ia,ib)=0
845	ifact4(4,ia,ib)=id
846	factb(1,ia,ib)=factb1
847	factb(2,ia,ib)=factb2
848	fact0(ia,ib)=signiifprefact(iia,id1)
849	endif
850	id11=1
851	id110=0
852	do ie=ic,0,-2
853	if(id.ne.0.or.ie.ne.0) then
854	if(ie.ne.ic) then
855	id110=id110+1
856	ib=ib+1
857	id11=id11+1
858	ifacta(1,ia,ib)=ie
859	ifacta(2,ia,ib)=id
860	ifact4(1,ia,ib)=ie+id110
861	ifact4(2,ia,ib)=id110
862	ifact4(3,ia,ib)=id
863	ifact4(4,ia,ib)=0
864	factb(1,ia,ib)=factb1
865	factb(2,ia,ib)=factb2
866	fact0(ia,ib)=signiifprefact(iia,id1)* &
867	&fact(ic,id11)
868	if(id.gt.0.and.ie.gt.0) then
869	ib=ib+1
870	ifacta(1,ia,ib)=ie
871	ifacta(2,ia,ib)=-id
872	ifact4(1,ia,ib)=ie+id110
873	ifact4(2,ia,ib)=id110
874	ifact4(3,ia,ib)=0
875	ifact4(4,ia,ib)=id
876	factb(1,ia,ib)=factb1
877	factb(2,ia,ib)=factb2
878	fact0(ia,ib)=signiifprefact(iia,id1)* &
879	&fact(ic,id11)
880	endif
881	endif
882	id12=1
883	id120=0
884	do if=id-2,0,-2
885	if(if.ne.0.or.ie.ne.0) then
886	id120=id120+1
887	ib=ib+1
888	id12=id12+1
889	ifacta(1,ia,ib)=ie
890	ifacta(2,ia,ib)=if
891	ifact4(1,ia,ib)=ie+id110
892	ifact4(2,ia,ib)=id110
893	ifact4(3,ia,ib)=if+id120
894	ifact4(4,ia,ib)=id120
895	factb(1,ia,ib)=factb1
896	factb(2,ia,ib)=factb2
897	fact0(ia,ib)=signiifprefact(iia,id1)* &
898	&fact(ic,id11)*fact(id,id12)
899	if(if.gt.0.and.ie.gt.0) then
900	ib=ib+1
901	ifacta(1,ia,ib)=ie
902	ifacta(2,ia,ib)=-if
903	ifact4(1,ia,ib)=ie+id110
904	ifact4(2,ia,ib)=id110
905	ifact4(3,ia,ib)=id120
906	ifact4(4,ia,ib)=if+id120
907	factb(1,ia,ib)=factb1
908	factb(2,ia,ib)=factb2
909	fact0(ia,ib)=signiifprefact(iia,id1)* &
910	&fact(ic,id11)*fact(id,id12)
911	endif
912	endif
913	enddo
914	id120=0
915	endif
916	enddo
917	id110=0
918	enddo
919	call sortres(2,ia,ib)
920	endif
921	enddo
922	!--Calculation of I Order Distortion Function
923	do iia=2,mmul
924	do ii=1,2
925	if(ii.eq.1) ia=iia
926	if(ii.eq.2) ia=-iia
927	if(ityc(ia).gt.0) then
928	write(6,10100) 'First Order Contribution of: ', &
929	&ityc(ia),comment(ia)
930	do i=1,ityc(ia)
931	do j=1,ivar(ia)
932	!--Factor of Machine Parameters
933	fact1(j)=bstr(ia,i)fact0(ia,j) &
934	&beta(1,ia,i)*factb(1,ia,j)beta(2,ia,i)**factb(2,ia,j)
935	!--Phase Factor
936	tc=sin(ifacta(1,ia,j)*(phi(1,ia,i)-qx)+ &
937	&ifacta(2,ia,j)*(phi(2,ia,i)-qy))
938	ts=cos(ifacta(1,ia,j)*(phi(1,ia,i)-qx)+ &
939	&ifacta(2,ia,j)*(phi(2,ia,i)-qy))
940	!--Generating Function Calculation
941	sinar=sin(ifacta(1,ia,j)qx+ifacta(2,ia,j)qy)
942	if(abs(sinar).gt.pieni) then
943	dstc(j)=fact1(j)*tc/sinar
944	dsts(j)=fact1(j)*ts/sinar
945	else
946	write(6,*) ' Warning: Results for distortion ', &
947	&'function probably wrong; sinar= ',sinar
948	dstc(j)=zero
949	dsts(j)=zero
950	endif
951	distc(ia,j)=distc(ia,j)+dstc(j)
952	dists(ia,j)=dists(ia,j)+dsts(j)
953	enddo
954	!--Write Generating Function Calculation of each Element
955	if(iu_on.eq.2.or.iu_on.eq.3) then
956	table_name ='distort_1_f_all '
957	data_size = ivar(ia)
958	call fit_table(table_name,data_size,table_size_76)
959	do j=1,data_size
960	amp = dsqrt(dstc(j)dstc(j)+dsts(j)dsts(j))
961	write(76,10000) ia,i,j,etl(ia,i),dstc(j),dsts(j),amp, &
962	&ifact4(1,ia,j),ifact4(2,ia,j), &
963	&ifact4(3,ia,j),ifact4(4,ia,j)
964	data_size = ivar(ia)
965	int_to_write(1) = ia
966	int_to_write(2) = i
967	int_to_write(3) = j
968	double_to_write(4) = etl(ia,i)
969	double_to_write(5) = dstc(j)
970	double_to_write(6) = dsts(j)
971	double_to_write(7) = amp
972	int_to_write(8) = ifact4(1,ia,j)
973	int_to_write(9) = ifact4(2,ia,j)
974	int_to_write(10) = ifact4(3,ia,j)
975	int_to_write(11) = ifact4(4,ia,j)
976	call write_table(table_name,11,int_to_write, &
977	&double_to_write)
978	j76 = j76 + 1
979	call augment_count(table_name)
980	enddo
981	endif
982	!--Calculate and write Hamiltonian of each Element
983	if(iu_on.eq.2.or.iu_on.eq.3) then
984	table_name ='distort_1_h_all '
985	data_size = ivar(ia)
986	call fit_table(table_name,data_size,table_size_77)
987	do j=1,data_size
988	facc=one-cos(two(ifacta(1,ia,j)qx+ &
989	&ifacta(2,ia,j)*qy))
990	facs=sin(two(ifacta(1,ia,j)qx+ifacta(2,ia,j)*qy))
991	hmc=faccdstc(j)+facsdsts(j)
992	hms=faccdsts(j)-facsdstc(j)
993	double_to_write(5)=hmc
994	double_to_write(6)=hms
995	amp = dsqrt(hmchmc+hmshms)
996	write(77,10000) ia,i,j,etl(ia,i),double_to_write(5), &
997	&double_to_write(6),amp,ifact4(1,ia,j),ifact4(2,ia,j), &
998	&ifact4(3,ia,j),ifact4(4,ia,j)
999	int_to_write(1) = ia
1000	int_to_write(2) = i
1001	int_to_write(3) = j
1002	double_to_write(4) = etl(ia,i)
1003	double_to_write(7) = amp
1004	int_to_write(8) = ifact4(1,ia,j)
1005	int_to_write(9) = ifact4(2,ia,j)
1006	int_to_write(10) = ifact4(3,ia,j)
1007	int_to_write(11) = ifact4(4,ia,j)
1008	call write_table(table_name,11,int_to_write, &
1009	&double_to_write)
1010	j77 = j77 + 1
1011	call augment_count(table_name)
1012	enddo
1013	endif
1014	enddo
1015	!--Write total Generating Function
1016	if(iu_on.eq.1.or.iu_on.eq.3) then
1017	write(6,10020)
1018	table_name ='distort_1_f_end '
1019	data_size = ivar(ia)
1020	call fit_table(table_name,data_size,table_size_74)
1021	do j=1,data_size
1022	amp = sqrt(distc(ia,j)distc(ia,j)+dists(ia,j) &
1023	&dists(ia,j))
1024	write(6,10040) distc(ia,j),dists(ia,j),amp, &
1025	&ifact4(1,ia,j),ifact4(2,ia,j),ifact4(3,ia,j), &
1026	&ifact4(4,ia,j)
1027	write(74,10010) ia,distc(ia,j),dists(ia,j),amp, &
1028	&ifact4(1,ia,j),ifact4(2,ia,j),ifact4(3,ia,j), &
1029	&ifact4(4,ia,j)
1030	int_to_write(1) = ia
1031	double_to_write(2) = distc(ia,j)
1032	double_to_write(3) = dists(ia,j)
1033	double_to_write(4) = amp
1034	int_to_write(5) = ifact4(1,ia,j)
1035	int_to_write(6) = ifact4(2,ia,j)
1036	int_to_write(7) = ifact4(3,ia,j)
1037	int_to_write(8) = ifact4(4,ia,j)
1038	call write_table(table_name,8,int_to_write, &
1039	&double_to_write)
1040	j74 = j74 + 1
1041	call augment_count(table_name)
1042	enddo
1043	!--Calculate and write total Hamiltonian
1044	write(6,10030)
1045	table_name ='distort_1_h_end '
1046	data_size = ivar(ia)
1047	call fit_table(table_name,data_size,table_size_75)
1048	do j=1,data_size
1049	facc=one-cos(two(ifacta(1,ia,j)qx+ifacta(2,ia,j)*qy))
1050	facs=sin(two(ifacta(1,ia,j)qx+ifacta(2,ia,j)*qy))
1051	hamc(ia,j)=faccdistc(ia,j)+facsdists(ia,j)
1052	hams(ia,j)=faccdists(ia,j)-facsdistc(ia,j)
1053	double_to_write(2)=hamc(ia,j)
1054	double_to_write(3)=hams(ia,j)
1055	amp = sqrt(hamc(ia,j)hamc(ia,j)+hams(ia,j)hams(ia,j))
1056	write(6,10040) double_to_write(2),double_to_write(3), &
1057	&amp,ifact4(1,ia,j),ifact4(2,ia,j), &
1058	&ifact4(3,ia,j),ifact4(4,ia,j)
1059	write(75,10010) ia,double_to_write(2), &
1060	&double_to_write(3),amp,ifact4(1,ia,j),ifact4(2,ia,j), &
1061	&ifact4(3,ia,j),ifact4(4,ia,j)
1062	int_to_write(1) = ia
1063	double_to_write(4) = amp
1064	int_to_write(5) = ifact4(1,ia,j)
1065	int_to_write(6) = ifact4(2,ia,j)
1066	int_to_write(7) = ifact4(3,ia,j)
1067	int_to_write(8) = ifact4(4,ia,j)
1068	call write_table(table_name,8,int_to_write, &
1069	&double_to_write)
1070	j75 = j75 + 1
1071	call augment_count(table_name)
1072	enddo
1073	endif
1074	endif
1075	enddo
1076	enddo
1077	call timex(tim5)
1078	write(6,10110) tim5-tim2
1079	!---------------------------------------------------------------------
1080	10000 format(i4,i7,i5,4(1pe23.15),4i4)
1081	10010 format(i4,4x,3(1pe23.15),6i4)
1082	10020 format(/80('-')/10x,'Distortionfunction Terms (first order)'/ &
1083	&80('-')/4x,'Cosine-Term',8x,'Sine-Term',8x,'Amplitude',7x,'J', &
1084	&5x,'K',5x,'L',5x,'M'/80('-'))
1085	10030 format(/80('-')/10x,'Hamiltonian Terms (first order)'/80('-')/ &
1086	&4x,'Cosine-Term',8x,'Sine-Term',8x,'Amplitude',7x,'J', &
1087	&5x,'K',5x,'L',5x,'M'/80('-'))
1088	10040 format(3(1pe17.9),4i6)
1089	10100 format(//80('-')//10x,a,i8,a)
1090	10110 format(/80('-')//'Program <DISTORT1> used ',f10.3, &
1091	&' second(s) of Execution Time'//80('-')//)
1092	return
1093	end subroutine distort1
1094	subroutine distort2
1095	!-----------------------------------------------------------------------
1096	!--Program Distort2 (second order)
1097	!-----------------------------------------------------------------------
1098	implicit none
1099	integer i,ia,ia0,ib,ib0,ic,ic0,ic00,ica,ica2,id,id0,id00,ide, &
1100	&ii,imax,imuty,itij0,jj,k,k1,k2,l,l0,l3,l4,ll3,ll4,m,mde,mh,mmul, &
1101	&mmul2,mmult,mmultf,mmultw,mmultx,nblz,num
1102	double precision arg,argij,betdx,betdxi,dl12q,dl34q,efact,ephix1, &
1103	&ephiy1,esin34,etrgij,facc,facd0,facd1,facd10,facp,facs,ffff,phix1,&
1104	&phix2,phix20,phiy1,phiy2,phiy20,pieni,sin34,spij,trgij
1105	integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, &
1106	&iu_on,n1,n2
1107	integer j70,j71,j72,j73,j74,j75,j76,j77,j78,j79
1108	integer table_size_70,table_size_71,table_size_72,table_size_73, &
1109	&table_size_74,table_size_75,table_size_76,table_size_77, &
1110	&table_size_78,table_size_79
1111	integer data_size
1112	double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, &
1113	&fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, &
1114	&qy,sbeta,sign,two,zero
1115	real tlim,time0,time1,time2,time3,time4
1116	real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
1117	character*16 strn,table_name
1118	character*18 comment
1119	parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100)
1120	parameter(mmul2=mmul+1,mmult=8mmul*2)
1121	parameter(mmultw=2mmul,mmultx=mmult/7,mmultf=2mmultx*mmult)
1122	common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul)
1123	common/c1/ tlim,time0,time1,time2,time3,time4
1124	common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
1125	common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz)
1126	common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz)
1127	common/c4q/ phi(2,-mmul:mmul,nblz)
1128	common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul)
1129	common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2
1130	common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult)
1131	common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx)
1132	common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul)
1133	common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf)
1134	common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul)
1135	common/c12/ ihv(mmultx),iplane(mmultx,2)
1136	common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh)
1137	common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh)
1138	common/indeces/j70,j71,j72,j73,j74,j75,j76,j77,j78,j79
1139	common/table_sizes/table_size_70,table_size_71,table_size_72, &
1140	&table_size_73,table_size_74,table_size_75, &
1141	&table_size_76,table_size_77,table_size_78, &
1142	&table_size_79
1143	dimension efact(mmultx,mmult),ll3(mmultx,mmult),ll4(mmultx,mmult)
1144	dimension betdxi(mmultx),betdx(mmultx)
1145	dimension facd1(mmultx,mmult,2)
1146	dimension itij0(mmultx,mmult,4),facd10(mmultx,mmult,2)
1147	dimension mde(mmul)
1148	integer int_to_write(11)
1149	double precision double_to_write(11)
1150	!---------------------------------------------------------------------
1151	!--Set the rest of the Variables to zero
1152	do i=1,mmult
1153	do ii=1,mmultx
1154	facd1(ii,i,1)=zero
1155	facd1(ii,i,2)=zero
1156	enddo
1157	enddo
1158	do i=1,mmul
1159	mde(i)=0
1160	enddo
1161	ii=0
1162	do i=4,mmul,2
1163	ii=ii+1
1164	mde(ii)=i
1165	enddo
1166	!---------------------------------------------------------------------
1167	!--
1168	!--Calculation of Driving-term Poisson Brackets
1169	!--
1170	!-----Needed for Second Order Terms
1171	!--
1172	!---------------------------------------------------------------------
1173	!---------------------------------------------------------------------
1174	!--
1175	!--Driving & Hamiltonian Terms Second Order in Multipole Strength
1176	!--
1177	!---------------------------------------------------------------------
1178	do ic00=n1,n2
1179	ic=ic00
1180	do id00=-iabs(ic),iabs(ic)
1181	id=id00
1182	ia=iabs(ic)
1183	ib=iabs(id)
1184	!--Difference between erect and skew
1185	!--imuty=1 ==> "erect"; imuty=2 ==> "skew";
1186	!--imuty=3 ==> "erect-skew"; imuty=4 ==> "skew-erect"
1187	if(ic.gt.0.and.id.gt.0) imuty=1
1188	if(ic.lt.0.and.id.lt.0) imuty=2
1189	if(ic.gt.0.and.id.lt.0) imuty=3
1190	if(ic.lt.0.and.id.gt.0) imuty=4
1191	call timex(tim3)
1192	if(ityc(ic).gt.0.and.ityc(id).gt.0) then
1193	icc=0
1194	imax=ia+ib-2
1195	do k=1,mmultx
1196	betexp(k,1)=zero
1197	betexp(k,2)=zero
1198	betexp(k,3)=zero
1199	betexp(k,4)=zero
1200	icd(k)=0
1201	fac2(k)=zero
1202	do m=1,mmult
1203	isig(k,m)=0
1204	itij(k,m,1)=0
1205	itij(k,m,2)=0
1206	itij(k,m,3)=0
1207	itij(k,m,4)=0
1208	itij(k,m,5)=0
1209	itij(k,m,6)=0
1210	itij0(k,m,1)=0
1211	itij0(k,m,2)=0
1212	itij0(k,m,3)=0
1213	itij0(k,m,4)=0
1214	fac4(k,m)=zero
1215	enddo
1216	enddo
1217	!--
1218	!--The ide(n0=0,yes=1) switch checks if the average(detuning)terms
1219	!--create contributions to the distortion function in higher order
1220	!--A full loop is needed
1221	!--
1222	ide=0
1223	if(imuty.eq.1.and.ia.ne.ib) then
1224	do k=1,mmul
1225	if(ib.eq.mde(k)) ide=1
1226	enddo
1227	endif
1228	if(ide.eq.1) then
1229	!--
1230	!--Extra (ide=1) Loop starts here
1231	!--
1232	ia0=ia
1233	ib0=ib
1234	ic0=ic
1235	id0=id
1236	ia=ib0
1237	ib=ia0
1238	ic=id0
1239	id=ic0
1240	call caldt2(ia,ib,imuty)
1241	do k=1,mmultx
1242	do m=1,mmult
1243	itij0(k,m,1)=itij(k,m,3)
1244	itij0(k,m,2)=itij(k,m,4)
1245	itij0(k,m,3)=itij(k,m,5)
1246	itij0(k,m,4)=itij(k,m,6)
1247	enddo
1248	enddo
1249	do ii=1,ityc(ic)
1250	!--Beta term of first set of multipoles
1251	do k=1,icc
1252	betdxi(k)=(beta(1,ic,ii)*betexp(k,1)) &
1253	&(beta(2,ic,ii)**betexp(k,2))
1254	enddo
1255	do jj=1,ityc(id)
1256	facd0=twobstr(ic,ii)bstr(id,jj)
1257	!--Phase Differences (phix2,phiy2) of secondary resonances change
1258	!--sign for reversed order of the position of multipoles
1259	spij=one
1260	if(phi(1,id,jj).lt.phi(1,ic,ii)) spij=-one
1261	!--Drivingterm Calculation
1262	do k=1,icc
1263	!--Beta term of second set of multipoles multiplied to first
1264	betdx(k)=betdxi(k)* &
1265	&(beta(1,id,jj)*betexp(k,3)) &
1266	&(beta(2,id,jj)**betexp(k,4))
1267	l0=icd(k)
1268	!--Cosine and sine part of driving term
1269	do l=1,l0
1270	if(isig(k,l).eq.1) then
1271	dl12q=itij(k,l,1)qx+itij(k,l,2)qy
1272	if(abs(dl12q).lt.pieni) dl12q=pieni
1273	l3=itij(k,l,3)-itij(k,l,4)
1274	l4=itij(k,l,5)-itij(k,l,6)
1275	dl34q=l3qx+l4qy
1276	if(abs(dl34q).lt.pieni) dl34q=pieni
1277	sin34=sin(dl34q)
1278	phix1=-qx
1279	phiy1=-qy
1280	trgij=sin(dl12q)
1281	if(phi(1,id,jj).eq.phi(1,ic,ii)) then
1282	phix2=phi(1,id,jj)+qx
1283	phiy2=phi(2,id,jj)+qy
1284	else
1285	phix2=phi(1,id,jj)-spij*qx
1286	phiy2=phi(2,id,jj)-spij*qy
1287	endif
1288	arg=l3phix1+l4phiy1+itij(k,l,1)* &
1289	&phix2+itij(k,l,2)*phiy2
1290	facp=fac4(k,l)facd0fac2(k)* &
1291	&betdx(k)/trgij/sin34
1292	facd10(k,l,1)=facd10(k,l,1)+facp*sin(arg)
1293	facd10(k,l,2)=facd10(k,l,2)+facp*cos(arg)
1294	endif
1295	enddo
1296	enddo
1297	enddo
1298	enddo
1299	ia=ia0
1300	ib=ib0
1301	ic=ic0
1302	id=id0
1303	endif
1304	!---------------------------------------------------------------------
1305	!--
1306	!--Main loop
1307	!--
1308	!---------------------------------------------------------------------
1309	icc=0
1310	do k=1,mmultx
1311	betexp(k,1)=zero
1312	betexp(k,2)=zero
1313	betexp(k,3)=zero
1314	betexp(k,4)=zero
1315	icd(k)=0
1316	fac2(k)=zero
1317	do m=1,mmult
1318	isig(k,m)=0
1319	itij(k,m,1)=0
1320	itij(k,m,2)=0
1321	itij(k,m,3)=0
1322	itij(k,m,4)=0
1323	itij(k,m,5)=0
1324	itij(k,m,6)=0
1325	fac4(k,m)=zero
1326	enddo
1327	enddo
1328	call caldt2(ia,ib,imuty)
1329	do k=1,icc
1330	l0=icd(k)
1331	do l=1,l0
1332	dl12q=itij(k,l,1)qx+itij(k,l,2)qy
1333	if(abs(dl12q).lt.pieni) dl12q=pieni
1334	etrgij=sin(dl12q)
1335	ll3(k,l)=itij(k,l,3)-itij(k,l,4)
1336	ll4(k,l)=itij(k,l,5)-itij(k,l,6)
1337	dl34q=ll3(k,l)qx+ll4(k,l)qy
1338	if(abs(dl34q).lt.pieni) dl34q=pieni
1339	esin34=sin(dl34q)
1340	efact(k,l)=fac4(k,l)/etrgij/esin34
1341	enddo
1342	enddo
1343	do ii=1,ityc(ic)
1344	!--Beta term of first set of multipoles
1345	do k=1,icc
1346	betdxi(k)=(beta(1,ic,ii)*betexp(k,1)) &
1347	&(beta(2,ic,ii)**betexp(k,2))
1348	enddo
1349	do jj=1,ityc(id)
1350	facd0=bstr(ic,ii)*bstr(id,jj)
1351	if(ic.ne.id.and.ic.ne.-id) facd0=facd0*two
1352	if(ii.eq.1.and.jj.eq.1) then
1353	if(ic.ne.id) then
1354	write(6,10090) 'Pairs of: ',ityc(ic),comment(ic), &
1355	&' and: ',ityc(id),comment(id)
1356	else
1357	write(6,10100) 'Quadratic Effect of: ', &
1358	&ityc(ic),comment(ic)
1359	endif
1360	endif
1361	!--Phase Differences (phix2,phiy2) of secondary resonances change
1362	!--sign for reversed order of the position of multipoles
1363	spij=one
1364	if(phi(1,id,jj).lt.phi(1,ic,ii)) spij=-one
1365	phix20=phi(1,id,jj)-phi(1,ic,ii)-spij*qx
1366	phiy20=phi(2,id,jj)-phi(2,ic,ii)-spij*qy
1367	ephix1=phi(1,ic,ii)-qx
1368	ephiy1=phi(2,ic,ii)-qy
1369	!--Drivingterm Calculation
1370	do k=1,icc
1371	!--Beta term of second set of multipoles multiplied to first
1372	betdx(k)=betdxi(k)* &
1373	&(beta(1,id,jj)*betexp(k,3)) &
1374	&(beta(2,id,jj)**betexp(k,4))
1375	ffff=facd0fac2(k)betdx(k)
1376	l0=icd(k)
1377	num=0
1378	do l=1,l0
1379	num=num+isig(k,l)
1380	enddo
1381	if (num.eq.0) then
1382	!--Cosine and sine part of driving term
1383	do l=1,l0
1384	arg=ll3(k,l)ephix1+ll4(k,l)ephiy1+ &
1385	&itij(k,l,1)phix20+itij(k,l,2)phiy20
1386	facp=ffff*efact(k,l)
1387	facd1(k,l,1)=facd1(k,l,1)+facp*sin(arg)
1388	facd1(k,l,2)=facd1(k,l,2)+facp*cos(arg)
1389	enddo
1390	else
1391	do l=1,l0
1392	if(isig(k,l).eq.0) then
1393	arg=ll3(k,l)ephix1+ll4(k,l)ephiy1+ &
1394	&itij(k,l,1)phix20+itij(k,l,2)phiy20
1395	facp=ffff*efact(k,l)
1396	facd1(k,l,1)=facd1(k,l,1)+facp*sin(arg)
1397	facd1(k,l,2)=facd1(k,l,2)+facp*cos(arg)
1398	else
1399	dl12q=itij(k,l,1)qx+itij(k,l,2)qy
1400	if(abs(dl12q).lt.pieni) dl12q=pieni
1401	l3=itij(k,l,3)-itij(k,l,4)
1402	l4=itij(k,l,5)-itij(k,l,6)
1403	dl34q=l3qx+l4qy
1404	if(abs(dl34q).lt.pieni) dl34q=pieni
1405	sin34=sin(dl34q)
1406	phix1=-qx
1407	phiy1=-qy
1408	trgij=sin(dl12q)
1409	phix2=phi(1,id,jj)
1410	phiy2=phi(2,id,jj)
1411	if(ic.eq.id.and.phi(1,id,jj).eq.phi(1,ic,ii)) &
1412	&then
1413	trgij=trgij/cos(dl12q)
1414	else
1415	phix2=phix2-spij*qx
1416	phiy2=phiy2-spij*qy
1417	endif
1418	arg=l3phix1+l4phiy1+itij(k,l,1)* &
1419	&phix2+itij(k,l,2)*phiy2
1420	facp=fac4(k,l)facd0fac2(k)* &
1421	&betdx(k)/trgij/sin34
1422	!--Factor 2 comes from trigonometric functions
1423	if(ic.eq.id) facp=two*facp
1424	facd1(k,l,1)=facd1(k,l,1)+facp*sin(arg)
1425	facd1(k,l,2)=facd1(k,l,2)+facp*cos(arg)
1426	endif
1427	enddo
1428	endif
1429	enddo
1430	enddo
1431	enddo
1432	!--Collect a total of ica non-zero terms
1433	ica=0
1434	do k1=1,mmultx
1435	do k2=1,mmult
1436	if(abs(facd1(k1,k2,1)).gt.pieni.or. &
1437	&abs(facd1(k1,k2,2)).gt.pieni) then
1438	ica=ica+1
1439	if(ica.gt.mmultf) call prror(3,10,mmultf)
1440	facd2(1,ica)=facd1(k1,k2,1)
1441	facd2(2,ica)=facd1(k1,k2,2)
1442	facd1(k1,k2,1)=zero
1443	facd1(k1,k2,2)=zero
1444	do ii=1,4
1445	ifacd2(ii,ica)=itij(k1,k2,ii+2)
1446	enddo
1447	endif
1448	enddo
1449	enddo
1450	!--Collect also the additional detuning induced terms
1451	if(ide.eq.1) then
1452	do k1=1,mmultx
1453	do k2=1,mmult
1454	if(abs(facd10(k1,k2,1)).gt.pieni.or. &
1455	&abs(facd10(k1,k2,2)).gt.pieni) then
1456	ica=ica+1
1457	if(ica.gt.mmultf) call prror(3,10,mmultf)
1458	facd2(1,ica)=facd10(k1,k2,1)
1459	facd2(2,ica)=facd10(k1,k2,2)
1460	facd10(k1,k2,1)=zero
1461	facd10(k1,k2,2)=zero
1462	do ii=1,4
1463	ifacd2(ii,ica)=itij0(k1,k2,ii)
1464	enddo
1465	endif
1466	enddo
1467	enddo
1468	endif
1469	!--Add up all contributions related to a specific resonance
1470	do k1=1,ica
1471	do k2=k1+1,ica
1472	if(ifacd2(1,k2).eq.ifacd2(1,k1).and. &
1473	&ifacd2(2,k2).eq.ifacd2(2,k1).and. &
1474	&ifacd2(3,k2).eq.ifacd2(3,k1).and. &
1475	&ifacd2(4,k2).eq.ifacd2(4,k1)) then
1476	facd2(1,k1)=facd2(1,k1)+facd2(1,k2)
1477	facd2(2,k1)=facd2(2,k1)+facd2(2,k2)
1478	facd2(1,k2)=zero
1479	facd2(2,k2)=zero
1480	do ii=1,4
1481	ifacd2(ii,k2)=0
1482	enddo
1483	endif
1484	enddo
1485	enddo
1486	ica2=0
1487	!--Collect the final total of ica2 non-zero resonance terms
1488	do k1=1,ica
1489	if(abs(facd2(1,k1)).gt.pieni.or. &
1490	&abs(facd2(2,k1)).gt.pieni) then
1491	ica2=ica2+1
1492	if(ica2.ne.k1) then
1493	facd2(1,ica2)=facd2(1,k1)
1494	facd2(2,ica2)=facd2(2,k1)
1495	facd2(1,k1)=zero
1496	facd2(2,k1)=zero
1497	do ii=1,4
1498	ifacd2(ii,ica2)=ifacd2(ii,k1)
1499	ifacd2(ii,k1)=0
1500	enddo
1501	endif
1502	endif
1503	enddo
1504	!--Order the resonance terms
1505	call sortres(3,imax,ica2)
1506	write(6,10060)
1507	table_name ='distort_2_f_end '
1508	data_size = ica2
1509	call fit_table(table_name,data_size,table_size_78)
1510	do k1=1,data_size
1511	argij=two((ifacd2(1,k1)-ifacd2(2,k1))qx+ &
1512	&(ifacd2(3,k1)-ifacd2(4,k1))*qy)
1513	facc=one-cos(argij)
1514	facs=sin(argij)
1515	ham(1,k1)=faccfacd2(1,k1)+facsfacd2(2,k1)
1516	ham(2,k1)=faccfacd2(2,k1)-facsfacd2(1,k1)
1517	double_to_write(5)=dsqrt(facd2(1,k1)*facd2(1,k1)+ &
1518	&facd2(2,k1)*facd2(2,k1))
1519	write(6,10080) facd2(1,k1),facd2(2,k1),double_to_write(5),&
1520	&ifacd2(1,k1),ifacd2(2,k1),ifacd2(3,k1),ifacd2(4,k1)
1521	if(iu_on.eq.1.or.iu_on.eq.3) &
1522	&write(78,10110) ic,id,facd2(1,k1),facd2(2,k1), &
1523	&double_to_write(5),ifacd2(1,k1),ifacd2(2,k1), &
1524	&ifacd2(3,k1),ifacd2(4,k1)
1525	int_to_write(1) = ic
1526	int_to_write(2) = id
1527	double_to_write(3) = facd2(1,k1)
1528	double_to_write(4) = facd2(2,k1)
1529	int_to_write(6) = ifacd2(1,k1)
1530	int_to_write(7) = ifacd2(2,k1)
1531	int_to_write(8) = ifacd2(3,k1)
1532	int_to_write(9) = ifacd2(4,k1)
1533	call write_table(table_name,9,int_to_write, &
1534	&double_to_write)
1535	j78 = j78 + 1
1536	call augment_count(table_name)
1537	facd2(1,k1)=zero
1538	facd2(2,k1)=zero
1539	enddo
1540	write(6,10070)
1541	table_name ='distort_2_h_end '
1542	data_size = ica2
1543	call fit_table(table_name,data_size,table_size_79)
1544	do k1=1,data_size
1545	double_to_write(5)=dsqrt(ham(1,k1)*ham(1,k1)+ &
1546	&ham(2,k1)*ham(2,k1))
1547	write(6,10080) ham(1,k1),ham(2,k1),double_to_write(5), &
1548	&ifacd2(1,k1),ifacd2(2,k1),ifacd2(3,k1),ifacd2(4,k1)
1549	if(iu_on.eq.1.or.iu_on.eq.3) &
1550	&write(79,10110) ic,id,ham(1,k1),ham(2,k1), &
1551	&double_to_write(5),ifacd2(1,k1),ifacd2(2,k1), &
1552	&ifacd2(3,k1),ifacd2(4,k1)
1553	int_to_write(1) = ic
1554	int_to_write(2) = id
1555	double_to_write(3) = ham(1,k1)
1556	double_to_write(4) = ham(2,k1)
1557	int_to_write(6) = ifacd2(1,k1)
1558	int_to_write(7) = ifacd2(2,k1)
1559	int_to_write(8) = ifacd2(3,k1)
1560	int_to_write(9) = ifacd2(4,k1)
1561	call write_table(table_name,9,int_to_write, &
1562	&double_to_write)
1563	j79 = j79 + 1
1564	call augment_count(table_name)
1565	do ii=1,4
1566	ifacd2(ii,k1)=0
1567	ham(1,k1)=zero
1568	ham(2,k1)=zero
1569	enddo
1570	enddo
1571	call timex(tim4)
1572	write(6,10030) tim4-tim3
1573	write(6,10020) tim4-tim2
1574	endif
1575	enddo
1576	enddo
1577	!---------------------------------------------------------------------
1578	call timex(tim5)
1579	write(6,10040) tim5-tim2
1580	return
1581	!---------------------------------------------------------------------
1582	10020 format(/80('-')//'Intermediate Time in <DISTORT2>: ', &
1583	&f10.3,' second(s)',' of Execution Time' &
1584	&//80('-')//)
1585	10030 format(/80('-')//'Second Order Drivingterm Calculation took ', &
1586	&f10.3,' second(s)',' of Execution Time'//80('-')//)
1587	10040 format(/80('-')//'Program <DISTORT2> used ',f10.3, &
1588	&' second(s) of Execution Time'//80('-')//)
1589	10060 format(/80('-')/10x,'Distortionfunction Terms (second order)'/ &
1590	&80('-')/4x,'Cosine-Term',8x,'Sine-Term',8x,'Amplitude',7x, &
1591	&'J',5x,'K',5x,'L',5x,'M'/80('-'))
1592	10070 format(/80('-')/10x,'Hamiltonian Terms (second order)'/80('-')/ &
1593	&4x,'Cosine-Term',8x,'Sine-Term',8x,'Amplitude',7x, &
1594	&'J',5x,'K',5x,'L',5x,'M'/80('-'))
1595	10080 format(3(1pe17.9),4i6)
1596	10090 format(//80('-')//10x,a,i8,2a,i8,a)
1597	10100 format(//80('-')//10x,a,i8,a)
1598	10110 format(i4,4x,i4,4x,3(1pe23.15),4i4)
1599	end subroutine distort2
1600	subroutine caldet2(ia,ib,imuty)
1601	!---------------------------------------------------------------------
1602	!--
1603	!--Calculation of the poisson bracket to derive the detuning function
1604	!--Second order in the strength of the multipoles
1605	!--
1606	!---------------------------------------------------------------------
1607	implicit none
1608	integer ia,ib,ic,id,ido,if2,ihoff,ihw,ii,ii1,ii11,iiend,iii,iii1, &
1609	&ijnu,ijnu1,imuty,ipl,iti,ivoff,jj,jj1,jj11,jjend,jjj,jjj1,k,k1,kk,&
1610	&l,l1,m,mend,mexp,mh,mij,mm,mmul,mmul2,mmult,mmultf,mmultw,mmultx, &
1611	&n,nblz,nend,nexp,nij,nn
1612	double precision betl,betxp,fac,fac1,fac3,fcc2,fcc4,pieni,signii, &
1613	&signjj
1614	integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, &
1615	&iu_on,n1,n2
1616	double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, &
1617	&fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, &
1618	&qy,sbeta,sign,two,zero
1619	real tlim,time0,time1,time2,time3,time4
1620	real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
1621	character*16 strn
1622	character*18 comment
1623	parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100)
1624	parameter(mmul2=mmul+1,mmult=8mmul*2)
1625	parameter(mmultw=2mmul,mmultx=mmult/7,mmultf=2mmultx*mmult)
1626	common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul)
1627	common/c1/ tlim,time0,time1,time2,time3,time4
1628	common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
1629	common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz)
1630	common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz)
1631	common/c4q/ phi(2,-mmul:mmul,nblz)
1632	common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul)
1633	common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2
1634	common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult)
1635	common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx)
1636	common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul)
1637	common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf)
1638	common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul)
1639	common/c12/ ihv(mmultx),iplane(mmultx,2)
1640	common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh)
1641	common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh)
1642	!---------------------------------------------------------------------
1643	dimension betl(mmultx,4),betxp(2,mmultx,4)
1644	dimension ihw(2,mmult),ipl(2,mmult,2),iti(2,mmultx,mmult,6)
1645	dimension fcc2(2,mmultx),fcc4(2,mmultx,mmult)
1646	do k=1,mmultx
1647	fcc2(1,k)=zero
1648	fcc2(2,k)=zero
1649	do m=1,mmult
1650	do nn=1,2
1651	ihw(nn,k)=0
1652	ipl(nn,k,1)=0
1653	ipl(nn,k,2)=0
1654	iti(nn,k,m,1)=0
1655	iti(nn,k,m,2)=0
1656	iti(nn,k,m,3)=0
1657	betl(k,1)=zero
1658	betl(k,2)=zero
1659	betl(k,3)=zero
1660	betl(k,4)=zero
1661	betxp(nn,k,1)=zero
1662	betxp(nn,k,2)=zero
1663	betxp(nn,k,3)=zero
1664	betxp(nn,k,4)=zero
1665	fcc4(nn,k,m)=zero
1666	enddo
1667	enddo
1668	enddo
1669	fac=-two/dble(42((ia+ib)/2)ia*ib)
1670	if(ia.ne.ib) fac=fac*two
1671	iiend=(ia+2)/2
1672	jjend=(ib+2)/2
1673	signii=-one
1674	if(imuty.eq.2.and.2*iiend.eq.ia+2) iiend=iiend-1
1675	do ii=1,iiend
1676	ii1=(ii-1)*2
1677	if(imuty.eq.2) ii1=ii1+1
1678	signii=signii*(-one)
1679	signjj=-one
1680	if(imuty.eq.2.and.2*jjend.eq.ib+2) jjend=jjend-1
1681	do jj=1,jjend
1682	if2=0
1683	jj1=(jj-1)*2
1684	if(imuty.eq.2) jj1=jj1+1
1685	signjj=signjj*(-one)
1686	fac1=facsigniisignjjfact(ia,1+ii1)fact(ib,1+jj1)
1687	iii=ia-ii1
1688	jjj=ib-jj1
1689	ijnu=(ii-1)*jjend+jj
1690	betl(ijnu,1)=dble(iii)/two
1691	betl(ijnu,2)=dble(ii1)/two
1692	betl(ijnu,3)=dble(jjj)/two
1693	betl(ijnu,4)=dble(jj1)/two
1694	ihoff=0
1695	ivoff=0
1696	if(iii.eq.0.and.jjj.ne.0) then
1697	fac1=fac1*fact(jjj,jjj/2+1)
1698	ihoff=1
1699	endif
1700	if(iii.ne.0.and.jjj.eq.0) then
1701	fac1=fac1*fact(iii,iii/2+1)
1702	ihoff=1
1703	endif
1704	if(ii1.eq.0.and.jj1.ne.0) then
1705	fac1=fac1*fact(jj1,jj1/2+1)
1706	ivoff=1
1707	endif
1708	if(ii1.ne.0.and.jj1.eq.0) then
1709	fac1=fac1*fact(ii1,ii1/2+1)
1710	ivoff=1
1711	endif
1712	k=(iii+jjj)/2
1713	k1=k-1
1714	l=(ii1+jj1)/2
1715	l1=l-1
1716	mij=min(iii,jjj)
1717	nij=min(ii1,jj1)
1718	mend=(mij+2)/2
1719	nend=(nij+2)/2
1720	iii1=0
1721	jjj1=0
1722	ii11=0
1723	jj11=0
1724	if(iii.gt.mij) iii1=(iii+2)/2-mend
1725	if(jjj.gt.mij) jjj1=(jjj+2)/2-mend
1726	if(ii1.gt.nij) ii11=(ii1+2)/2-nend
1727	if(jj1.gt.nij) jj11=(jj1+2)/2-nend
1728	!---------------------------------------------------------------------
1729	!--
1730	!--Deriving the Derivatives
1731	!--
1732	!---------------------------------------------------------------------
1733	if(k.gt.0) then
1734	if(l.gt.0.and.ivoff.eq.0.and.ihoff.eq.0) then
1735	!--
1736	!--Qx Term first Derivative by Iy for the Detuning by Ix
1737	!--
1738	ido=1
1739	!--Finding Double Appearance
1740	do ijnu1=1,ijnu-1
1741	if( &
1742	&betl(ijnu1,1).eq.betl(ijnu,3).and. &
1743	&betl(ijnu1,2).eq.betl(ijnu,4).and. &
1744	&betl(ijnu1,3).eq.betl(ijnu,1).and. &
1745	&betl(ijnu1,4).eq.betl(ijnu,2)) then
1746	ido=0
1747	if2=if2+1
1748	if(if2.gt.2) call prror(1,13,if2)
1749	fcc2(if2,ijnu1)=two*fcc2(if2,ijnu1)
1750	endif
1751	enddo
1752	if(ido.eq.1) then
1753	if2=if2+1
1754	if(if2.gt.2) call prror(1,13,if2)
1755	ihw(if2,ijnu)=1
1756	ipl(if2,ijnu,1)=k1
1757	ipl(if2,ijnu,2)=l1
1758	betxp(if2,ijnu,1)=betl(ijnu,1)
1759	betxp(if2,ijnu,2)=betl(ijnu,2)
1760	betxp(if2,ijnu,3)=betl(ijnu,3)
1761	betxp(if2,ijnu,4)=betl(ijnu,4)
1762	fcc2(if2,ijnu)=fac1dble(kii1jj1)/dble(2*(k1+l1))
1763	do m=1,mend
1764	mexp=mij-(m-1)*2
1765	fac3=fact(iii,iii1+m)*fact(jjj,jjj1+m)
1766	do n=1,nend
1767	nexp=nij-(n-1)*2
1768	nn=(m-1)*nend+n
1769	fcc4(if2,ijnu,nn)= &
1770	&fac3fact(ii1,ii11+n)fact(jj1,jj11+n)* &
1771	&dble(nexp)*(one/dble(ii1)+one/dble(jj1))/two
1772	if(mexp.ne.0) then
1773	iti(if2,ijnu,nn,1)=-1
1774	iti(if2,ijnu,nn,2)=mexp
1775	iti(if2,ijnu,nn,3)=nexp
1776	else
1777	iti(if2,ijnu,nn,1)=0
1778	iti(if2,ijnu,nn,2)=0
1779	iti(if2,ijnu,nn,3)=nexp
1780	endif
1781	enddo
1782	enddo
1783	endif
1784	endif
1785	if(k1.gt.0.and.ihoff.eq.0) then
1786	!--
1787	!--Qx Term first Derivative by Ix for the Detuning by Ix
1788	!--
1789	ido=1
1790	!--Finding Double Appearance
1791	do ijnu1=1,ijnu-1
1792	if( &
1793	&betl(ijnu1,1).eq.betl(ijnu,3).and. &
1794	&betl(ijnu1,2).eq.betl(ijnu,4).and. &
1795	&betl(ijnu1,3).eq.betl(ijnu,1).and. &
1796	&betl(ijnu1,4).eq.betl(ijnu,2)) then
1797	ido=0
1798	if2=if2+1
1799	if(if2.gt.2) call prror(1,13,if2)
1800	fcc2(if2,ijnu1)=two*fcc2(if2,ijnu1)
1801	endif
1802	enddo
1803	if(ido.eq.1) then
1804	if2=if2+1
1805	if(if2.gt.2) call prror(1,13,if2)
1806	ihw(if2,ijnu)=1
1807	ipl(if2,ijnu,1)=k1-1
1808	ipl(if2,ijnu,2)=l
1809	betxp(if2,ijnu,1)=betl(ijnu,1)
1810	betxp(if2,ijnu,2)=betl(ijnu,2)
1811	betxp(if2,ijnu,3)=betl(ijnu,3)
1812	betxp(if2,ijnu,4)=betl(ijnu,4)
1813	fcc2(if2,ijnu)=fac1dble(k1iii*jjj)/ &
1814	&dble(2**(k1+l1))
1815	!--Terms Iy**n where n>0
1816	if(ii1.ne.0.and.jj1.ne.0) then
1817	do n=1,nend
1818	nexp=nij-(n-1)*2
1819	fac3=fact(ii1,ii11+n)*fact(jj1,jj11+n)
1820	do m=1,mend
1821	mexp=mij-(m-1)*2
1822	mm=(n-1)*mend+m
1823	fcc4(if2,ijnu,mm)= &
1824	&fac3fact(iii,iii1+m)fact(jjj,jjj1+m)* &
1825	&dble(mexp)*(one/dble(iii)+one/dble(jjj))/two
1826	if(abs(fcc4(if2,ijnu,mm)).gt.pieni) then
1827	if(nexp.ne.0) then
1828	iti(if2,ijnu,mm,1)=1
1829	iti(if2,ijnu,mm,2)=mexp
1830	iti(if2,ijnu,mm,3)=nexp
1831	else if(mexp.ne.0) then
1832	iti(if2,ijnu,mm,1)=0
1833	iti(if2,ijnu,mm,2)=mexp
1834	iti(if2,ijnu,mm,3)=0
1835	endif
1836	endif
1837	enddo
1838	enddo
1839	!--Terms Iy**0=1
1840	else
1841	do m=1,mend
1842	mexp=mij-(m-1)*2
1843	fcc4(if2,ijnu,m)= &
1844	&fact(iii,iii1+m)fact(jjj,jjj1+m) &
1845	&dble(mexp)*(one/dble(iii)+one/dble(jjj))/two
1846	if(abs(fcc4(if2,ijnu,m)).gt.pieni) then
1847	iti(if2,ijnu,m,1)=0
1848	iti(if2,ijnu,m,2)=mexp
1849	iti(if2,ijnu,m,3)=0
1850	endif
1851	enddo
1852	endif
1853	endif
1854	else if(ihoff.eq.1.and.ivoff.eq.0) then
1855	!--
1856	!--Horizontal Term only one x, first Derivative by Iy and then Ix
1857	!--
1858	ido=1
1859	!--Finding Double Appearance
1860	do ijnu1=1,ijnu-1
1861	if( &
1862	&betl(ijnu1,1).eq.betl(ijnu,3).and. &
1863	&betl(ijnu1,2).eq.betl(ijnu,4).and. &
1864	&betl(ijnu1,3).eq.betl(ijnu,1).and. &
1865	&betl(ijnu1,4).eq.betl(ijnu,2)) then
1866	ido=0
1867	if2=if2+1
1868	if(if2.gt.2) call prror(1,13,if2)
1869	fcc2(if2,ijnu1)=two*fcc2(if2,ijnu1)
1870	endif
1871	enddo
1872	if(ido.eq.1) then
1873	if2=if2+1
1874	if(if2.gt.2) call prror(1,13,if2)
1875	ihw(if2,ijnu)=1
1876	ipl(if2,ijnu,1)=k1
1877	ipl(if2,ijnu,2)=l1
1878	betxp(if2,ijnu,1)=betl(ijnu,1)
1879	betxp(if2,ijnu,2)=betl(ijnu,2)
1880	betxp(if2,ijnu,3)=betl(ijnu,3)
1881	betxp(if2,ijnu,4)=betl(ijnu,4)
1882	fcc2(if2,ijnu)=fac1dble(kii1jj1)/dble(2*(k1+l1))
1883	do n=1,nend
1884	nexp=nij-(n-1)*2
1885	fcc4(if2,ijnu,n)= &
1886	&fact(ii1,ii11+n)fact(jj1,jj11+n) &
1887	&dble(nexp)*(one/dble(ii1)+one/dble(jj1))/two
1888	iti(if2,ijnu,n,1)=0
1889	iti(if2,ijnu,n,2)=0
1890	iti(if2,ijnu,n,3)=nexp
1891	enddo
1892	endif
1893	else if(iii.eq.1.and.jjj.eq.1) then
1894	!--
1895	!--Vertical Term for uneven Multipoles
1896	!--
1897	if2=if2+1
1898	if(if2.gt.2) call prror(1,13,if2)
1899	ihw(if2,ijnu)=2
1900	ipl(if2,ijnu,1)=0
1901	ipl(if2,ijnu,2)=l-1
1902	betxp(if2,ijnu,1)=betl(ijnu,1)
1903	betxp(if2,ijnu,2)=betl(ijnu,2)
1904	betxp(if2,ijnu,3)=betl(ijnu,3)
1905	betxp(if2,ijnu,4)=betl(ijnu,4)
1906	fcc2(if2,ijnu)=fac1dble(liiijjj)/dble(2*(l-1))
1907	do n=1,nend
1908	fac3=fact(ii1,ii11+n)*fact(jj1,jj11+n)
1909	nexp=nij-(n-1)*2
1910	fcc4(if2,ijnu,n)=fac3* &
1911	&(one/dble(iii)+one/dble(jjj))/two
1912	if(nexp.ne.0) then
1913	iti(if2,ijnu,n,1)=1
1914	iti(if2,ijnu,n,2)=1
1915	iti(if2,ijnu,n,3)=nexp
1916	else
1917	iti(if2,ijnu,n,1)=0
1918	iti(if2,ijnu,n,2)=1
1919	iti(if2,ijnu,n,3)=nexp
1920	endif
1921	enddo
1922	endif
1923	!--
1924	!--Vertical Term for even Multipoles
1925	!--
1926	else if(ivoff.eq.0) then
1927	if2=if2+1
1928	if(if2.gt.2) call prror(1,13,if2)
1929	ihw(if2,ijnu)=2
1930	ipl(if2,ijnu,1)=0
1931	ipl(if2,ijnu,2)=l1-1
1932	betxp(if2,ijnu,1)=betl(ijnu,1)
1933	betxp(if2,ijnu,2)=betl(ijnu,2)
1934	betxp(if2,ijnu,3)=betl(ijnu,3)
1935	betxp(if2,ijnu,4)=betl(ijnu,4)
1936	fcc2(if2,ijnu)=fac1dble(l1ii1jj1)/dble(2*(l-2))
1937	do n=1,nend
1938	nexp=nij-(n-1)*2
1939	fcc4(if2,ijnu,n)= &
1940	&fact(ii1,ii11+n)fact(jj1,jj11+n) &
1941	&dble(nexp)*(one/dble(ii1)+one/dble(jj1))/two
1942	iti(if2,ijnu,n,1)=0
1943	iti(if2,ijnu,n,2)=0
1944	iti(if2,ijnu,n,3)=nexp
1945	enddo
1946	endif
1947	enddo
1948	enddo
1949	!--Write out Parameters properly
1950	ic=0
1951	do nn=1,2
1952	do kk=1,mmultx
1953	if(abs(fcc2(nn,kk)).gt.pieni) then
1954	id=0
1955	ic=ic+1
1956	if(ic.gt.mmultx) call prror(1,6,mmultx)
1957	ihv(ic)=ihw(nn,kk)
1958	iplane(ic,1)=ipl(nn,kk,1)
1959	iplane(ic,2)=ipl(nn,kk,2)
1960	betexp(ic,1)=betxp(nn,kk,1)
1961	betexp(ic,2)=betxp(nn,kk,2)
1962	betexp(ic,3)=betxp(nn,kk,3)
1963	betexp(ic,4)=betxp(nn,kk,4)
1964	fac2(ic)=fcc2(nn,kk)
1965	do mm=1,mmult
1966	if(abs(fcc4(nn,kk,mm)).gt.pieni) then
1967	id=id+1
1968	if(id.gt.mmult) call prror(1,7,mmult)
1969	itij(ic,id,1)=iti(nn,kk,mm,1)
1970	itij(ic,id,2)=iti(nn,kk,mm,2)
1971	itij(ic,id,3)=iti(nn,kk,mm,3)
1972	fac4(ic,id)=fcc4(nn,kk,mm)
1973	endif
1974	enddo
1975	icd(ic)=id
1976	endif
1977	enddo
1978	icc=ic
1979	enddo
1980	return
1981	end subroutine caldet2
1982	subroutine caldt2(ia,ib,imuty)
1983	!---------------------------------------------------------------------
1984	!--
1985	!--Calculation of the poisson bracket to derive the distortion function
1986	!--Second order in the strength of the multipoles
1987	!--
1988	!---------------------------------------------------------------------
1989	implicit none
1990	integer ia,ib,ic,ica,icol1,icol2,icol3,icol4,icol5,icol6,id,ido1, &
1991	&ido2,ido3,ido4,ihoff,ihoff1,ihoff2,ihvall,ii,ii1,iiend,iii,ijnu, &
1992	&imuty,isiga,iti,ivoff,ivoff1,ivoff2,jj,jj1,jjend,jjj,k,k1,k2,k3, &
1993	&k4,k41,k412s,k412s1,k412s2,k41s,k42,k42s,k43,k434s,k434s1,k434s2, &
1994	&k44,k5,kk,l1,m,mh,mm,mmul,mmul2,mmult,mmultf,mmultw,mmultx,nblz, &
1995	&nn,nnanf,nnend,nsig
1996	double precision betl,diab,fac11,fac12,fac21,fac22,faca,facc, &
1997	&facc1,facc2,fcc2,fcc4,pieni,signii,signjj
1998	integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, &
1999	&iu_on,n1,n2
2000	double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, &
2001	&fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, &
2002	&qy,sbeta,sign,two,zero
2003	real tlim,time0,time1,time2,time3,time4
2004	real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
2005	character*16 strn
2006	character*18 comment
2007	parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100)
2008	parameter(mmul2=mmul+1,mmult=8mmul*2)
2009	parameter(mmultw=2mmul,mmultx=mmult/7,mmultf=2mmultx*mmult)
2010	common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul)
2011	common/c1/ tlim,time0,time1,time2,time3,time4
2012	common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
2013	common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz)
2014	common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz)
2015	common/c4q/ phi(2,-mmul:mmul,nblz)
2016	common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul)
2017	common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2
2018	common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult)
2019	common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx)
2020	common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul)
2021	common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf)
2022	common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul)
2023	common/c12/ ihv(mmultx),iplane(mmultx,2)
2024	common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh)
2025	common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh)
2026	dimension iti(2,mmultx,mmult,6),isiga(2,mmultx,mmult)
2027	dimension fcc2(2,mmultx),fcc4(2,mmultx,mmult)
2028	dimension betl(mmultx,4)
2029	!---------------------------------------------------------------------
2030	do k=1,mmultx
2031	fcc2(1,k)=zero
2032	fcc2(2,k)=zero
2033	do m=1,mmult
2034	do nn=1,2
2035	iti(nn,k,m,1)=0
2036	iti(nn,k,m,2)=0
2037	iti(nn,k,m,3)=0
2038	iti(nn,k,m,4)=0
2039	iti(nn,k,m,5)=0
2040	iti(nn,k,m,6)=0
2041	isiga(nn,k,m)=0
2042	betl(k,1)=zero
2043	betl(k,2)=zero
2044	betl(k,3)=zero
2045	betl(k,4)=zero
2046	fcc4(nn,k,m)=zero
2047	enddo
2048	enddo
2049	enddo
2050	diab=dble(ia+ib)/two
2051	facc=one/four/(two*diab)/dble(iaib)
2052	!--Prefactor I (General factors)
2053	iiend=(ia+2)/2
2054	jjend=(ib+2)/2
2055	!--(iiend,jjend)=>Total number of terms in (x+iy)**n
2056	!--(i,j) is on left or right side of the Poisson bracket respectively
2057	signii=-one
2058	if((imuty.eq.4.or.imuty.eq.2).and.2*iiend.eq.ia+2) iiend=iiend-1
2059	do ii=1,iiend
2060	ii1=(ii-1)*2
2061	if(imuty.eq.4.or.imuty.eq.2) ii1=ii1+1
2062	signii=signii*(-one)
2063	signjj=-one
2064	if((imuty.eq.3.or.imuty.eq.2).and.2*jjend.eq.ib+2) jjend=jjend-1
2065	do jj=1,jjend
2066	jj1=(jj-1)*2
2067	!--(ii1,jj1)=>Increment of powers (0->{(iiend,jjend)-1}*2)
2068	!-- goes in steps of 2
2069	!--These are the exponents of the second variable y
2070	if(imuty.eq.3.or.imuty.eq.2) jj1=jj1+1
2071	signjj=signjj*(-one)
2072	facc1=faccsigniisignjjfact(ia,1+ii1)fact(ib,1+jj1)
2073	!--Prefactor II (signs, binominal coefficients)
2074	iii=ia-ii1
2075	jjj=ib-jj1
2076	!--(iii,jjj)=>Exponent of the first variable x
2077	ijnu=(ii-1)*jjend+jj
2078	if(ijnu.gt.mmultx) call prror(3,6,mmultx)
2079	!--inju=># of current case treated of a total of (iiend*jjend)
2080	betl(ijnu,1)=dble(iii)/two
2081	betl(ijnu,2)=dble(ii1)/two
2082	betl(ijnu,3)=dble(jjj)/two
2083	betl(ijnu,4)=dble(jj1)/two
2084	ihoff=0
2085	ihoff1=0
2086	ihoff2=0
2087	ivoff=0
2088	ivoff1=0
2089	ivoff2=0
2090	if(iii.eq.0.and.jjj.eq.0) ihoff=1
2091	if(iii.eq.0.and.jjj.ne.0.and.ihoff.eq.0) ihoff1=1
2092	if(iii.ne.0.and.jjj.eq.0.and.ihoff.eq.0) ihoff2=1
2093	if(ii1.eq.0.and.jj1.eq.0) ivoff=1
2094	if(ii1.eq.0.and.jj1.ne.0.and.ivoff.eq.0) ivoff1=1
2095	if(ii1.ne.0.and.jj1.eq.0.and.ivoff.eq.0) ivoff2=1
2096	ihvall=ihoff+ihoff1+ihoff2+ivoff+ivoff1+ivoff2
2097	!--(ihoff,ivoff)=0/1 (x,y) on (at least on one side/on no sides) of
2098	!-- Poisson bracket
2099	!--(ihoff1,ivoff1)=0/1 (x,y) on (both sides/on left side) of
2100	!-- Poisson bracket
2101	!--(ihoff2,ivoff2)=0/1 (x,y) on (both sides/on right side) of
2102	!-- Poisson bracket
2103	facc2=-facc1/(two**(diab-one))/four
2104	!--Prefactor III (Invariant->emittance) and -1/4d0 from evaluation of
2105	!-- trigonometric functions
2106	!---------------------------------------------------------------------
2107	!--
2108	!--Deriving the Derivatives
2109	!--
2110	!---------------------------------------------------------------------
2111	nnanf=1
2112	nnend=2
2113	if(ihoff.eq.1.or.ihoff1.eq.1.or.ihoff2.eq.1) then
2114	nnanf=2
2115	nnend=2
2116	endif
2117	if(ivoff.eq.1.or.ivoff1.eq.1.or.ivoff2.eq.1) nnend=1
2118	if(ihvall.eq.0) then
2119	nnanf=1
2120	nnend=2
2121	endif
2122	!--Differentiation nn=1 => horizontal nn=2 => vertical
2123	!-- ==> exchange planes (ido1,ido2 => ido3,ido4)
2124	!-- ==> exchange terms of secondary resonances (icol1 => icol2)
2125	!-- ==> exchange primary resonance terms (icol3,icol4 => icol5,icol6)
2126	do nn=nnanf,nnend
2127	if(nn.eq.1) then
2128	ido1=iii
2129	ido2=jjj
2130	ido3=ii1
2131	ido4=jj1
2132	icol1=1
2133	icol2=2
2134	icol3=3
2135	icol4=4
2136	icol5=5
2137	icol6=6
2138	!--nsig see below
2139	nsig=3
2140	endif
2141	if(nn.eq.2) then
2142	ido1=ii1
2143	ido2=jj1
2144	ido3=iii
2145	ido4=jjj
2146	icol1=2
2147	icol2=1
2148	icol3=5
2149	icol4=6
2150	icol5=3
2151	icol6=4
2152	!--nsig is set to 4 for the vertical differentiation as the
2153	!--horizontal stays always positive
2154	nsig=4
2155	endif
2156	fcc2(nn,ijnu)=facc2dble(ido1ido2)
2157	ica=0
2158	!--(k1,k2) always differentiation nn=1 horizontal nn=2 vertical
2159	do k1=0,ido1,2
2160	k41=ido1-k1
2161	!--k41 see k44
2162	fac11=dble(k41)/dble(ido1)*fact(ido1,1+k1/2)
2163	fac12=fact(ido1,1+k1/2)
2164	do k2=0,ido2,2
2165	k42=ido2-k2
2166	!--k42 see k44
2167	fac21=fact(ido2,1+k2/2)
2168	fac22=dble(k42)/dble(ido2)*fact(ido2,1+k2/2)
2169	!--(l1) sums over the 3-4 possible sign combinations
2170	do l1=1,nsig
2171	faca=sign(l1,1)fac11fac21- &
2172	&sign(l1,2)fac12fac22
2173	if(k41.eq.0) faca=-sign(l1,2)fac12fac22
2174	if(k42.eq.0) faca=sign(l1,1)fac11fac21
2175	!--Condition 1 & 5 clear
2176	!--2 Condition avoids duplication of l1 by l2
2177	!--3 Condition avoids duplication of l1 by l3
2178	!--4 Condition avoids duplication of l1 with any other
2179	if((abs(faca).gt.pieni).and. &
2180	&(l1.ne.2.or.k41.ne.0).and. &
2181	&(l1.ne.3.or.k42.ne.0).and. &
2182	&(l1.ne.4.or.(k41.ne.0.and.k42.ne.0)).and. &
2183	&((k41+k42).ge.1)) then
2184	!--(k3,k4) plane with no differentiation nn=1 vertical nn=2 horizontal
2185	!--((ido3+3),(ido4+3)) allow for constant term in the
2186	!--undifferentiated part
2187	do k3=0,2*ido3,2
2188	do k4=0,2*ido4,2
2189	!--k43,k44 control the secondary resonance
2190	!--i.e. the same that appear in first order
2191	k43=ido3-k3
2192	k44=ido4-k4
2193	!--(k41s,k42s) takes care of the sign of k41 and k42
2194	k41s=sign(l1,1)*k41
2195	k42s=sign(l1,2)*k42
2196	!--(k412s) sum of the 2 terms with differentiation
2197	k412s=k41s+k42s
2198	k412s1=(ido1+ido2+k412s)/2-1
2199	k412s2=k412s1-k412s
2200	!--(k434s) sum of the 2 terms without differentiation
2201	k434s=k43+k44
2202	k434s1=(ido3+ido4+k434s)/2
2203	k434s2=k434s1-k434s
2204	!--1. (0,0) Secondary Resonance has to be excluded
2205	!--2. Select resonances in the desired nomenclature
2206	!--at least the sum must be larger than 0
2207	!--and the resonance should be written in descending order starting
2208	!--with the horizontal one
2209	if( &
2210	&(k42s.ne.0.or.k44.ne.0).and. &
2211	&(k412s.gt.0.or.k434s.gt.0).and. &
2212	&((nn.eq.1.and.k412s1.ge.k412s2).or. &
2213	&(nn.eq.2.and.k434s1.ge.k434s2))) then
2214	!--Program check: finds duplicated cases
2215	do k5=1,ica
2216	if( &
2217	&iti(nn,ijnu,k5,icol1).eq.k42s.and. &
2218	&iti(nn,ijnu,k5,icol2).eq.k44.and. &
2219	&iti(nn,ijnu,k5,icol3).eq.k412s1.and. &
2220	&iti(nn,ijnu,k5,icol4).eq.k412s2.and. &
2221	&iti(nn,ijnu,k5,icol5).eq.k434s1.and. &
2222	&iti(nn,ijnu,k5,icol6).eq.k434s2 &
2223	&) call prror(3,9,2)
2224	enddo
2225	ica=ica+1
2226	if(ica.gt.mmult) call prror(3,5,mmult)
2227	fcc4(nn,ijnu,ica)=faca* &
2228	&fact(ido3,1+k3/2)*fact(ido4,1+k4/2)
2229	iti(nn,ijnu,ica,icol1)=k42s
2230	iti(nn,ijnu,ica,icol2)=k44
2231	iti(nn,ijnu,ica,icol3)=k412s1
2232	iti(nn,ijnu,ica,icol4)=k412s2
2233	iti(nn,ijnu,ica,icol5)=k434s1
2234	iti(nn,ijnu,ica,icol6)=k434s2
2235	if(k41.eq.0.and.k43.eq.0) isiga(nn,ijnu,ica)=1
2236	endif
2237	enddo
2238	enddo
2239	endif
2240	enddo
2241	enddo
2242	enddo
2243	enddo
2244	enddo
2245	enddo
2246	!--Write out Parameters properly
2247	ic=0
2248	do nn=1,2
2249	do kk=1,mmultx
2250	if(abs(fcc2(nn,kk)).gt.pieni) then
2251	id=0
2252	ic=ic+1
2253	if(ic.gt.mmultx) call prror(3,6,mmultx)
2254	betexp(ic,1)=betl(kk,1)
2255	betexp(ic,2)=betl(kk,2)
2256	betexp(ic,3)=betl(kk,3)
2257	betexp(ic,4)=betl(kk,4)
2258	fac2(ic)=fcc2(nn,kk)
2259	do mm=1,mmult
2260	if(abs(fcc4(nn,kk,mm)).gt.pieni) then
2261	id=id+1
2262	if(id.gt.mmult) call prror(3,7,mmult)
2263	if(isiga(nn,kk,mm).eq.1) then
2264	isig(ic,id)=1
2265	isiga(nn,kk,mm)=0
2266	endif
2267	itij(ic,id,1)=iti(nn,kk,mm,1)
2268	itij(ic,id,2)=iti(nn,kk,mm,2)
2269	itij(ic,id,3)=iti(nn,kk,mm,3)
2270	itij(ic,id,4)=iti(nn,kk,mm,4)
2271	itij(ic,id,5)=iti(nn,kk,mm,5)
2272	itij(ic,id,6)=iti(nn,kk,mm,6)
2273	fac4(ic,id)=fcc4(nn,kk,mm)
2274	endif
2275	enddo
2276	icd(ic)=id
2277	endif
2278	enddo
2279	icc=ic
2280	enddo
2281	return
2282	end subroutine caldt2
2283	subroutine init
2284	!---------------------------------------------------------------------
2285	!--Initilisation
2286	!---------------------------------------------------------------------
2287	implicit none
2288	integer i,ii,j,k,mh,mmul,mmul2,mmult,mmultf,mmultw,mmultx,nblz
2289	double precision pieni
2290	integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, &
2291	&iu_on,n1,n2
2292	double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, &
2293	&fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, &
2294	&qy,sbeta,sign,two,zero
2295	real tlim,time0,time1,time2,time3,time4
2296	real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
2297	character*16 strn
2298	character*18 comment
2299	parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100)
2300	parameter(mmul2=mmul+1,mmult=8mmul*2)
2301	parameter(mmultw=2mmul,mmultx=mmult/7,mmultf=2mmultx*mmult)
2302	common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul)
2303	common/c1/ tlim,time0,time1,time2,time3,time4
2304	common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
2305	common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz)
2306	common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz)
2307	common/c4q/ phi(2,-mmul:mmul,nblz)
2308	common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul)
2309	common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2
2310	common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult)
2311	common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx)
2312	common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul)
2313	common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf)
2314	common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul)
2315	common/c12/ ihv(mmultx),iplane(mmultx,2)
2316	common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh)
2317	common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh)
2318	!---------------------------------------------------------------------
2319	tlim=1e7
2320	zero=dble(0d0)
2321	one=1d0
2322	two=2d0
2323	four=4d0
2324	pi=atan(one)*four
2325	pi2=pi*two
2326	pihi=two/pi
2327	comment(2)=' Erect Quadrupoles'
2328	comment(-2)=' Skew Quadrupoles '
2329	comment(3)=' Erect Sextupoles '
2330	comment(-3)=' Skew Sextupoles '
2331	comment(4)=' Erect Octupoles '
2332	comment(-4)=' Skew Octupoles '
2333	comment(5)=' Erect Decapoles '
2334	comment(-5)=' Skew Decapoles '
2335	comment(6)=' Erect Dodecapoles'
2336	comment(-6)=' Skew Dodecapoles '
2337	comment(7)=' Erect 14th-Pole '
2338	comment(-7)=' Skew 14th-Pole '
2339	comment(8)=' Erect 16th-Pole '
2340	comment(-8)=' Skew 16th-Pole '
2341	comment(9)=' Erect 18th-Pole '
2342	comment(-9)=' Skew 18th-Pole '
2343	comment(10)=' Erect 20th-Pole '
2344	comment(-10)=' Skew 20th-Pole '
2345	comment(11)=' Erect 22th-Pole '
2346	comment(-11)=' Skew 22th-Pole '
2347	!--
2348	!--Binominals
2349	!--
2350	cosav(0)=one
2351	do i=2,mmul,2
2352	ii=i/2
2353	cosav(ii)=cosav(ii-1)*dble(i-1)/dble(i)
2354	enddo
2355	do j=0,mmul
2356	do i=0,mmul2
2357	fact(j,i)=one
2358	enddo
2359	enddo
2360	do i=2,mmul
2361	do j=2,mmul
2362	if(j.ge.i) then
2363	fact(j,i)=fact(j-1,i)+fact(j-1,i-1)
2364	endif
2365	enddo
2366	enddo
2367	do i=-mmul,mmul
2368	do j=0,mmul
2369	do k=0,mmul
2370	det(1,i,j,k)=zero
2371	det(2,i,j,k)=zero
2372	enddo
2373	enddo
2374	enddo
2375	do i=0,mmul
2376	det1(1,i)=zero
2377	det1(2,i)=zero
2378	enddo
2379	!--
2380	!--sign array
2381	!--
2382	sign(1,1)=one
2383	sign(1,2)=one
2384	sign(2,1)=-one
2385	sign(2,2)=one
2386	sign(3,1)=one
2387	sign(3,2)=-one
2388	sign(4,1)=-one
2389	sign(4,2)=-one
2390	return
2391	end subroutine init
2392	subroutine readdat
2393	!---------------------------------------------------------------------
2394	!--
2395	!--Reading Lattice Information from File # 34
2396	!--
2397	!---------------------------------------------------------------------
2398	implicit none
2399	integer i,ich1,ich2,ii,isize,itype,j,jj,mh,mmul,mmul2,mmult, &
2400	&mmultf,mmultw,mmultx,nblz,ndum,num
2401	double precision betx,bety,et,phix0,phiy0,pieni,str
2402	character*16 name
2403	character*200 ch
2404	character*205 ch1
2405	integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, &
2406	&iu_on,n1,n2
2407	double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, &
2408	&fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, &
2409	&qy,sbeta,sign,two,zero, factor
2410	real tlim,time0,time1,time2,time3,time4
2411	real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
2412	character*16 strn
2413	character*18 comment
2414	parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100)
2415	parameter(mmul2=mmul+1,mmult=8mmul*2)
2416	parameter(mmultw=2mmul,mmultx=mmult/7,mmultf=2mmultx*mmult)
2417	common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul)
2418	common/c1/ tlim,time0,time1,time2,time3,time4
2419	common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
2420	common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz)
2421	common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz)
2422	common/c4q/ phi(2,-mmul:mmul,nblz)
2423	common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul)
2424	common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2
2425	common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult)
2426	common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx)
2427	common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul)
2428	common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf)
2429	common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul)
2430	common/c12/ ihv(mmultx),iplane(mmultx,2)
2431	common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh)
2432	common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh)
2433	!---------------------------------------------------------------------
2434
2435	do i=1,nblz
2436	do j=-mmul,mmul
2437	ityc(j)=0
2438	bstr(j,i)=zero
2439	strn(j,i)=" "
2440	beta(1,j,i)=zero
2441	beta(2,j,i)=zero
2442	sbeta(1,j,i)=zero
2443	sbeta(2,j,i)=zero
2444	phi(1,j,i)=zero
2445	phi(2,j,i)=zero
2446	etl(j,i)=zero
2447	enddo
2448	enddo
2449	do ii=1,10000000
2450	read(34,'(A150)',end=900) ch
2451	ich1=0
2452	ich2=0
2453	do jj=1,200
2454	if(ich1.eq.0.and.ch(jj:jj).ne.' ') ich1=1
2455	if(ich1.eq.1.and.ch(jj:jj).eq.' ') ich1=2
2456	if(ich1.eq.2.and.ch(jj:jj).ne.' ') then
2457	ich1=3
2458	ich2=jj
2459	endif
2460	if(ich1.eq.3.and.ch(jj:jj).eq.' ') then
2461	ch1(1:200)=ch(1:ich2-1)//''''//ch(ich2:jj-1)//''''// &
2462	&ch(jj:200)//' / '
2463	goto 210
2464	endif
2465	enddo
2466	210 continue
2467	read(ch1,*) et,name,itype,str,betx,bety,phix0,phiy0
2468	if(itype.eq.100) then
2469	qx=phix0*pi
2470	qy=phiy0*pi
2471	goto 5
2472	endif
2473	if(et.ge.etl1.and.et.le.etl2) then
2474	if(abs(str).ge.pieni.and.itype.ge.n1.and.itype.le.n2) then
2475	ityc(itype)=ityc(itype)+1
2476	isize=ityc(itype)
2477	if(isize.gt.nblz) call prror(3,3,nblz)
2478	etl(itype,isize)=et
2479	strn(itype,isize)=name
2480	factor = 10*(3(iabs(itype)-2))
2481	bstr(itype,isize)=str*factor
2482	beta(1,itype,isize)=betx
2483	beta(2,itype,isize)=bety
2484	sbeta(1,itype,isize)=dsqrt(betx)
2485	sbeta(2,itype,isize)=dsqrt(bety)
2486	phi(1,itype,isize)=phix0twopi
2487	phi(2,itype,isize)=phiy0twopi
2488	endif
2489	endif
2490	enddo
2491	call prror(3,4,num)
2492	5 continue
2493	write(6,*) ' Qx= ',qx/pi,' Qy= ',qy/pi
2494	close(34)
2495	return
2496	900 close(34)
2497	call prror(3,1,ndum)
2498	end subroutine readdat
2499	subroutine detwri(icase,imu,ih1)
2500	!---------------------------------------------------------------------
2501	!--Organize Writing for Detune1
2502	!---------------------------------------------------------------------
2503	implicit none
2504	integer i,ic,icase,id,ih1,imu,j,k,mh,mmul,mmul2,mmult,mmultf, &
2505	&mmultw,mmultx,nblz
2506	double precision pieni
2507	integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, &
2508	&iu_on,n1,n2
2509	double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, &
2510	&fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2, &
2511	&pihi,qx,qy,sbeta,sign,two,zero
2512	real tlim,time0,time1,time2,time3,time4
2513	real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
2514	integer j70,j71,j72,j73,j74,j75,j76,j77,j78,j79
2515	integer table_size_70,table_size_71,table_size_72,table_size_73, &
2516	&table_size_74,table_size_75,table_size_76,table_size_77, &
2517	&table_size_78,table_size_79
2518	character*16 strn,table_name
2519	character*18 comment
2520	parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100)
2521	parameter(mmul2=mmul+1,mmult=8mmul*2)
2522	parameter(mmultw=2mmul,mmultx=mmult/7,mmultf=2mmultx*mmult)
2523	common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul)
2524	common/c1/ tlim,time0,time1,time2,time3,time4
2525	common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
2526	common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz)
2527	common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz)
2528	common/c4q/ phi(2,-mmul:mmul,nblz)
2529	common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul)
2530	common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2
2531	common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult)
2532	common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx)
2533	common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul)
2534	common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf)
2535	common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul)
2536	common/c12/ ihv(mmultx),iplane(mmultx,2)
2537	common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh)
2538	common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh)
2539	common/indeces/j70,j71,j72,j73,j74,j75,j76,j77,j78,j79
2540	common/table_sizes/table_size_70,table_size_71,table_size_72, &
2541	&table_size_73,table_size_74,table_size_75, &
2542	&table_size_76,table_size_77,table_size_78, &
2543	&table_size_79
2544	integer int_to_write(11)
2545	double precision double_to_write(11)
2546	dimension imu(2)
2547	integer data_size
2548	!---------------------------------------------------------------------
2549	ic=imu(1)
2550	id=imu(2)
2551	if(icase.eq.0) then
2552
2553	!--- Zero order multiple strength order
2554
2555	!--- handle unit 71 (detune_1_all)
2556
2557	if(iu_on.eq.2.or.iu_on.eq.3) then
2558	table_name = 'detune_1_all '
2559	data_size = 2*ih1+2+j71
2560	call fit_table(table_name,data_size,table_size_71)
2561	do i=ih1,0,-1
2562	write(71,10040) ic,1,det(1,ic,i,ih1-i)/pi2,i,ih1-i
2563	int_to_write(1) = ic
2564	int_to_write(2) = 1
2565	double_to_write(3) = det(1,ic,i,ih1-i)/pi2
2566	int_to_write(4) = i
2567	int_to_write(5) = ih1-i
2568	call write_table(table_name,5,int_to_write,double_to_write)
2569	j71 = j71 + 1
2570	call augment_count(table_name)
2571	enddo
2572	do i=ih1-1,0,-1
2573	write(71,10040) ic,2, &
2574	&det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)/pi2,i+1,ih1-i-1
2575	int_to_write(2) = 2
2576	double_to_write(3) = det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)&
2577	&/pi2
2578	int_to_write(4) = i+1
2579	int_to_write(5) = ih1-i-1
2580	call write_table(table_name,5,int_to_write,double_to_write)
2581	j71 = j71 + 1
2582	call augment_count(table_name)
2583	enddo
2584	write(71,10040) ic,2,det(2,ic,0,ih1)/pi2,0,ih1
2585	double_to_write(3) = det(2,ic,0,ih1)/pi2
2586	int_to_write(4) = 0
2587	int_to_write(5) = ih1
2588	call write_table(table_name,5,int_to_write,double_to_write)
2589	j71 = j71 + 1
2590	call augment_count(table_name)
2591	endif
2592
2593	!--- First order multiple strength order
2594
2595	!--- handle unit 70 (detune_1_end)
2596
2597	else if(icase.eq.1) then
2598	write(6,10000)
2599	table_name = 'detune_1_end '
2600	data_size = 2*ih1+2+j70
2601	call fit_table(table_name,data_size,table_size_70)
2602	do i=ih1,0,-1
2603	write(6,10010) det(1,ic,i,ih1-i)/pi2,i,ih1-i
2604	if(iu_on.eq.1.or.iu_on.eq.3) then
2605	write(70,10040) ic,1,det(1,ic,i,ih1-i)/pi2,i,ih1-i
2606	int_to_write(1) = ic
2607	int_to_write(2) = 1
2608	double_to_write(3) = det(1,ic,i,ih1-i)/pi2
2609	int_to_write(4) = i
2610	int_to_write(5) = ih1-i
2611	call write_table(table_name,5,int_to_write,double_to_write)
2612	j70 = j70 + 1
2613	call augment_count(table_name)
2614	endif
2615	enddo
2616	write(6,10020)
2617	do i=ih1-1,0,-1
2618	write(6,10010) det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)/pi2, &
2619	&i+1,ih1-i-1
2620	if(iu_on.eq.1.or.iu_on.eq.3) then
2621	write(70,10040) ic,2,det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)&
2622	&/pi2,i+1,ih1-i-1
2623	int_to_write(2) = 2
2624	double_to_write(3) = det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)&
2625	&/pi2
2626	int_to_write(4) = i+1
2627	int_to_write(5) = ih1-i-1
2628	call write_table(table_name,5,int_to_write,double_to_write)
2629	j70 = j70 + 1
2630	call augment_count(table_name)
2631	endif
2632	enddo
2633	write(6,10010) det(2,ic,0,ih1)/pi2,0,ih1
2634	if(iu_on.eq.1.or.iu_on.eq.3) then
2635	write(70,10040) ic,2,det(2,ic,0,ih1)/pi2,0,ih1
2636	double_to_write(3) = det(2,ic,0,ih1)/pi2
2637	int_to_write(4) = 0
2638	int_to_write(5) = ih1
2639	call write_table(table_name,5,int_to_write,double_to_write)
2640	j70 = j70 + 1
2641	call augment_count(table_name)
2642	endif
2643	else if(icase.eq.2) then
2644
2645	!--- Second order multiple strength order
2646
2647	!--- handle unit 72 (detune_2_hor)
2648
2649	if(iu_on.eq.1.or.iu_on.eq.3) then
2650	write(6,10000)
2651	table_name = 'detune_2_hor '
2652	data_size = ih1+1+j72
2653	call fit_table(table_name,data_size,table_size_72)
2654	endif
2655	do i=ih1,0,-1
2656	write(6,10010) det(1,ic,i,ih1-i)/pi2,i,ih1-i
2657	if(iu_on.eq.1.or.iu_on.eq.3) then
2658	write(72,10050) ic,id,det(1,ic,i,ih1-i)/pi2,i,ih1-i
2659	int_to_write(1) = ic
2660	int_to_write(2) = id
2661	double_to_write(3) = det(1,ic,i,ih1-i)/pi2
2662	int_to_write(4) = i
2663	int_to_write(5) = ih1-i
2664	call write_table(table_name,5,int_to_write,double_to_write)
2665	j72 = j72 + 1
2666	call augment_count(table_name)
2667	endif
2668	enddo
2669
2670	!--- handle unit 73 (detune_2_ver)
2671
2672	if(iu_on.eq.1.or.iu_on.eq.3) then
2673	write(6,10020)
2674	table_name = 'detune_2_ver '
2675	data_size = ih1+2+j73
2676	call fit_table(table_name,data_size,table_size_73)
2677	endif
2678	do i=ih1-1,0,-1
2679	write(6,10010) det(1,ic,i,ih1-i)*dble(ih1-i)/ &
2680	&dble(i+1)/pi2,i+1,ih1-i-1
2681	if(iu_on.eq.1.or.iu_on.eq.3) then
2682	write(73,10050) ic,id, &
2683	&det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)/pi2,i+1,ih1-i-1
2684	int_to_write(1) = ic
2685	int_to_write(2) = id
2686	double_to_write(3) = det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)&
2687	&/pi2
2688	int_to_write(4) = i+1
2689	int_to_write(5) = ih1-i-1
2690	call write_table(table_name,5,int_to_write,double_to_write)
2691	j73 = j73 + 1
2692	call augment_count(table_name)
2693	endif
2694	enddo
2695	write(6,10010) det(2,ic,0,ih1)/pi2,0,ih1
2696	if(iu_on.eq.1.or.iu_on.eq.3) then
2697	write(73,10050) ic,id, &
2698	&det(2,ic,0,ih1)/pi2,0,ih1
2699	double_to_write(3) = det(2,ic,0,ih1)/pi2
2700	int_to_write(4) = 0
2701	int_to_write(5) = ih1
2702	call write_table(table_name,5,int_to_write,double_to_write)
2703	j73 = j73 + 1
2704	call augment_count(table_name)
2705	endif
2706	do i=-mmul,mmul
2707	do j=0,mmul
2708	do k=0,mmul
2709	det(1,i,j,k)=zero
2710	det(2,i,j,k)=zero
2711	enddo
2712	enddo
2713	enddo
2714	do i=0,mmul
2715	det1(1,i)=zero
2716	det1(2,i)=zero
2717	enddo
2718	endif
2719	return
2720	!---------------------------------------------------------------------
2721	10000 format(/80('-')/8x,'Horizontal Coefficient',4x,'E-x',3x,'E-y'/ &
2722	&80('-'))
2723	10010 format(10x,1pe20.12,2i6)
2724	10020 format(/80('-')/8x,' Vertical Coefficient',4x,'E-x',3x,'E-y'/ &
2725	&80('-'))
2726	10040 format(i4,3x,i4,2x,1pe23.15,i4,1x,i4)
2727	10050 format(i4,4x,i4,3x,1pe23.15,1x,i4,1x,i4)
2728
2729	end subroutine detwri
2730	subroutine sortres(icase,im1,ic)
2731	!---------------------------------------------------------------------
2732	!--
2733	!--Order the resonance terms
2734	!--
2735	!---------------------------------------------------------------------
2736	implicit none
2737	integer ic,icase,idum1,idum2,idum3,idum4,idum5,idum6,im,im1,is, &
2738	&is0,k0,k1,k2,k3,k4,l2,l4,m,mh,mmul,mmul2,mmult,mmultf,mmultw, &
2739	&mmultx,nblz
2740	double precision dum1,dum2,dum3,dum4,pieni
2741	integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, &
2742	&iu_on,n1,n2
2743	double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, &
2744	&fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, &
2745	&qy,sbeta,sign,two,zero
2746	real tlim,time0,time1,time2,time3,time4
2747	real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
2748	character*16 strn
2749	character*18 comment
2750	parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100)
2751	parameter(mmul2=mmul+1,mmult=8mmul*2)
2752	parameter(mmultw=2mmul,mmultx=mmult/7,mmultf=2mmultx*mmult)
2753	common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul)
2754	common/c1/ tlim,time0,time1,time2,time3,time4
2755	common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8
2756	common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz)
2757	common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz)
2758	common/c4q/ phi(2,-mmul:mmul,nblz)
2759	common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul)
2760	common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2
2761	common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult)
2762	common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx)
2763	common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul)
2764	common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf)
2765	common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul)
2766	common/c12/ ihv(mmultx),iplane(mmultx,2)
2767	common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh)
2768	common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh)
2769	!---------------------------------------------------------------------
2770	im=iabs(im1)
2771	is0=0
2772	is=0
2773	do k0=im,0,-1
2774	do k1=k0,0,-1
2775	do k2=k1,0,-1
2776	l2=k1-k2
2777	do k3=k0-k1,k0-k1,-1
2778	do k4=k3,0,-1
2779	l4=k3-k4
2780	do m=1,ic
2781	if(icase.eq.2.and. &
2782	&k2.eq.ifact4(1,im1,m).and.l2.eq.ifact4(2,im1,m) &
2783	&.and.k4.eq.ifact4(3,im1,m).and. &
2784	&l4.eq.ifact4(4,im1,m)) then
2785	is0=is0+1
2786	dum1=fact0(im1,is0)
2787	dum2=factb(1,im1,is0)
2788	dum3=factb(2,im1,is0)
2789	idum1=ifacta(1,im1,is0)
2790	idum2=ifacta(2,im1,is0)
2791	idum3=ifact4(1,im1,is0)
2792	idum4=ifact4(2,im1,is0)
2793	idum5=ifact4(3,im1,is0)
2794	idum6=ifact4(4,im1,is0)
2795	fact0(im1,is0)=fact0(im1,m)
2796	factb(1,im1,is0)=factb(1,im1,m)
2797	factb(2,im1,is0)=factb(2,im1,m)
2798	ifacta(1,im1,is0)=ifacta(1,im1,m)
2799	ifacta(2,im1,is0)=ifacta(2,im1,m)
2800	ifact4(1,im1,is0)=ifact4(1,im1,m)
2801	ifact4(2,im1,is0)=ifact4(2,im1,m)
2802	ifact4(3,im1,is0)=ifact4(3,im1,m)
2803	ifact4(4,im1,is0)=ifact4(4,im1,m)
2804	fact0(im1,m)=dum1
2805	factb(1,im1,m)=dum2
2806	factb(2,im1,m)=dum3
2807	ifacta(1,im1,m)=idum1
2808	ifacta(2,im1,m)=idum2
2809	ifact4(1,im1,m)=idum3
2810	ifact4(2,im1,m)=idum4
2811	ifact4(3,im1,m)=idum5
2812	ifact4(4,im1,m)=idum6
2813	endif
2814	if(icase.eq.3.and. &
2815	&k2.eq.ifacd2(1,m).and.l2.eq.ifacd2(2,m) &
2816	&.and.k4.eq.ifacd2(3,m).and.l4.eq.ifacd2(4,m)) then
2817	is=is+1
2818	dum1=facd2(1,is)
2819	dum2=facd2(2,is)
2820	dum3=ham(1,is)
2821	dum4=ham(2,is)
2822	idum1=ifacd2(1,is)
2823	idum2=ifacd2(2,is)
2824	idum3=ifacd2(3,is)
2825	idum4=ifacd2(4,is)
2826	facd2(1,is)=facd2(1,m)
2827	facd2(2,is)=facd2(2,m)
2828	ham(1,is)=ham(1,m)
2829	ham(2,is)=ham(2,m)
2830	ifacd2(1,is)=ifacd2(1,m)
2831	ifacd2(2,is)=ifacd2(2,m)
2832	ifacd2(3,is)=ifacd2(3,m)
2833	ifacd2(4,is)=ifacd2(4,m)
2834	facd2(1,m)=dum1
2835	facd2(2,m)=dum2
2836	ham(1,m)=dum3
2837	ham(2,m)=dum4
2838	ifacd2(1,m)=idum1
2839	ifacd2(2,m)=idum2
2840	ifacd2(3,m)=idum3
2841	ifacd2(4,m)=idum4
2842	endif
2843	enddo
2844	enddo
2845	enddo
2846	enddo
2847	enddo
2848	enddo
2849	return
2850	end subroutine sortres
2851	subroutine prror(ipro,ier,num)
2852	implicit none
2853	integer ier,ipro,num
2854	!-----------------------------------------------------------------------
2855	!--Error Output
2856	!-----------------------------------------------------------------------
2857
2858	if(ipro.eq.0) write(6,*) ' Error occurred in Reading User Input '
2859	if(ipro.eq.1) write(6,*) ' Error occurred in Program Detune, ', &
2860	&' which calculates the first and second Order Detuning.'
2861	if(ipro.eq.2) write(6,*) ' Error occurred in Program Distort1, ',&
2862	&' which calculates the Distortion Function in first Order.'
2863	if(ipro.eq.3) write(6,*) ' Error occurred in Program Distort2, ',&
2864	&' which calculates the Distortion Function in second Order.'
2865	goto(10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160),ier
2866	10 write(6,*) ' File 34 empty or corrupted - Program stops'
2867	stop 1
2868	20 write(6,*) ' File 34 cannot be read - Program stops'
2869	stop 2
2870	30 write(6,*) ' Number of Elements nblz: ',num,' too small!'
2871	stop 3
2872	40 write(6,*) ' Input File fc.34 corrupted '
2873	stop 4
2874	50 write(6,*) ' Too many Cases: in Derivation mmult:',num
2875	stop 5
2876	60 write(6,*) ' Too many Cases: ic terms, increase mmult: ',num
2877	stop 6
2878	70 write(6,*) ' Too many Cases: id terms, increase mmultx: ',num
2879	stop 7
2880	80 write(6,*) ' Number of Resonance Terms exceeded, increase ', &
2881	&'mmultx: ',num
2882	stop 8
2883	90 write(6,*) ' Program Error unphysical Terms in Case: ',num
2884	stop 9
2885	100 write(6,*) ' Maximum Number of Resonances: ',num,' too large ', &
2886	&'==> increase mmultf'
2887	stop 10
2888	110 write(6,*) ' Inconsistent Number: ',num,' of Resonances for ', &
2889	&'the mixed Case'
2890	stop 11
2891	120 write(6,*) ' Final Ordering of Resonances is not possible as ', &
2892	&'their Number: ',num,' is too large.'
2893	stop 12
2894	130 write(6,*) ' There are only 2 derivatives possible - program', &
2895	&' error'
2896	stop 13
2897	140 write(6,*) ' Program Specifier iprog: ',num,' must lie between', &
2898	&' 1 and 7 '
2899	stop 14
2900	150 write(6,*) ' Analysis is restricted to order: ',num
2901	stop 15
2902	160 write(6,*) ' Unit: ',num,' could not be opened '
2903	stop 16
2904	end subroutine prror
2905
2906	subroutine fit_table(table_name,data_size,tab_size)
2907	implicit none
2908	integer data_size,tab_size,ll
2909	character table_name()
2910	character*16 tab_name
2911
2912	tab_name = table_name
2913	ll = len(tab_name)
2914	tab_name(ll:ll) = ' '
2915	10 if(data_size .gt. tab_size) then
2916	call double_table(tab_name)
2917	tab_size = tab_size*2
2918	goto 10
2919	endif
2920
2921	return
2922	end subroutine fit_table
2923
2924	subroutine ertab(k,table_name,column,row)
2925	implicit none
2926	integer k,row
2927	character table_name(),column()
2928	if (k .eq.0) return
2929	if (k .eq. -1) print *,"table ",table_name," does not exist"
2930	if (k .eq. -2) print *," in table ",table_name,"column ",column, &
2931	&" does not exist"
2932	if (k .eq. -3) print *,"in table ",table_name,"row",row, &
2933	&" does not exist"
2934
2935	return
2936	end subroutine ertab
2937
2938	subroutine write_table(table_name,table_type,int_to_write, &
2939	&double_to_write)
2940	implicit none
2941	integer table_type,int_to_write(11)
2942	integer k,table_type_index(11)
2943	integer j70,j71,j72,j73,j74,j75,j76,j77,j78,j79
2944	integer tab_types_5(5),tab_types_8(8),tab_types_9(9), &
2945	&tab_types_11(11)
2946	double precision double_to_write(11)
2947	character table_name()
2948	character*16 name_5(5),name_8(8),name_9(9),name_11(11)
2949	common/indeces/j70,j71,j72,j73,j74,j75,j76,j77,j78,j79
2950	data table_type_index/0,0,0,0,1,0,0,2,3,0,4/
2951	data tab_types_5/1,1,2,1,1/
2952	data tab_types_8/1,2,2,2,1,1,1,1/
2953	data tab_types_9/1,1,2,2,2,1,1,1,1/
2954	data tab_types_11/1,1,1,2,2,2,2,1,1,1,1/
2955	data name_5/'multipoleorder ','plane ','detuning ','h_inv_order ',&
2956	&'v_inv_order '/
2957	data name_8/'multipoleorder ','cosine ','sine ','amplitude ', &
2958	&'j ','k ','l ','m '/
2959	data name_9/'multipoleorder1 ','multipoleorder2 ','cosine ', &
2960	&'sine ','amplitude ','j ','k ','l ','m '/
2961	data name_11/'multipoleorder ','location ','resonance ', &
2962	&'position[m] ','cosine ','sine ','amplitude ','j ','k ','l ','m '/
2963
2964	goto(50,80,90,110) table_type_index(table_type)
2965	50 do k = 1,5
2966	if(tab_types_5(k) .eq. 2) then
2967	call double_to_table_curr(table_name,name_5(k),double_to_write(k))
2968	else
2969	call double_to_table_curr(table_name,name_5(k), &
2970	&dble(int_to_write(k)))
2971	endif
2972	enddo
2973	go to 1000
2974
2975	80 do k = 1,8
2976	if(tab_types_8(k) .eq. 2) then
2977	call double_to_table_curr(table_name,name_8(k),double_to_write(k))
2978	else
2979	call double_to_table_curr(table_name,name_8(k), &
2980	&dble(int_to_write(k)))
2981	endif
2982	enddo
2983	go to 1000
2984
2985	90 do k = 1,9
2986	if(tab_types_9(k) .eq. 2) then
2987	call double_to_table_curr(table_name,name_9(k),double_to_write(k))
2988	else
2989	call double_to_table_curr(table_name,name_9(k), &
2990	&dble(int_to_write(k)))
2991	endif
2992	enddo
2993	goto 1000
2994
2995	110 do k = 1,11
2996	if(tab_types_11(k) .eq. 2) then
2997	call double_to_table_curr(table_name,name_11(k),double_to_write(k))
2998	else
2999	call double_to_table_curr(table_name,name_11(k), &
3000	&dble(int_to_write(k)))
3001	endif
3002	enddo
3003
3004	1000 return
3005	end subroutine write_table

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

Context Navigation

source: PSPA/madxPSPA/src/sodd.f90 @ 445

Download in other formats: