Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

LieLib.f @ 3

Last change on this file since 3 was 3, checked in by zhangj, 12 years ago
Initiale import
Property svn:executable set to ``*
File size: 109.4 KB

Line
1	subroutine lieinit(no1,nv1,nd1,ndpt1,iref1,nis)
2	implicit none
3	integer i,iref1,nd1,ndc1,ndim,ndpt1,nis,no1,nv1
4	double precision ang,ra,st
5	!! Lieinit initializes AD Package and Lielib
6	parameter (ndim=3)
7	dimension st(ndim),ang(ndim),ra(ndim)
8	integer nd,nd2,no,nv
9	common /ii/no,nv,nd,nd2
10	integer ndc,ndc2,ndpt,ndt
11	common /coast/ndc,ndc2,ndt,ndpt
12	integer iref
13	common /resfile/iref
14	integer itu
15	common /tunedef/itu
16	integer lienot
17	common /dano/lienot
18	integer idpr
19	common /printing/ idpr
20	double precision xintex
21	common /integratedex/ xintex(0:20)
22	integer idao,is,iscrri
23	double precision rs
24	common/dascr/is(100),rs(100),iscrri(100),idao
25	integer nplane
26	double precision epsplane,xplane
27	common /choice/ xplane(ndim),epsplane,nplane(ndim)
28	!+CA DASCR
29	call daexter
30	do i=1,ndim
31	nplane(i)=2*i-1
32	ang(i)=0.d0
33	ra(i)=0.d0
34	st(i)=1.d0
35	enddo
36	no=no1
37	nv=nv1
38	nd=nd1
39	nd2=2*nd1
40	do i=1,100
41	is(i)=0
42	enddo
43	write(, ) "lieinit: calling daini"
44	call daini(no,nv,0)
45	if(nis.gt.0)call etallnom(is,nis,'$$IS ')
46	if(ndpt1.eq.0) then
47	ndpt=0
48	ndt=0
49	ndc1=0
50	else
51	ndpt=ndpt1
52	ndc1=1
53	if(ndpt.eq.nd2) then
54	ndt=nd2-1
55	else
56	ndt=nd2
57	if(ndpt.ne.nd2-1) then
58	write(6,*) ' LETHAL ERROR IN LIEINIT'
59	stop
60	endif
61	endif
62	endif
63	ndc=ndc1
64	ndc2=2*ndc1
65	iref=0
66	call initpert(st,ang,ra)
67	iref=iref1
68	if(iref1.eq.0) then
69	itu=0
70	else
71	itu=1
72	endif
73	if(iref1.eq.0) iref=-1
74
75	if(idpr.eq.1)write(6,*) ' NO = ',no,' IN DA-CALCULATIONS '
76
77	do i=0,20
78	xintex(i)=0.d0
79	enddo
80	xintex( 0)= 1.000000000000000
81	xintex( 1)= 5.000000000000000e-001
82	xintex( 2)= 8.333333333333334e-002
83	xintex( 3)= 0.000000000000000e+000
84	xintex( 4)= -1.388888888888898e-003
85	xintex( 5)= 0.000000000000000e+000
86	xintex( 6)= 3.306878306878064e-005
87	xintex( 7)= 0.d0
88	xintex( 8)= -8.267195767165669e-007
89	xintex( 9)= 0.d0
90	xintex( 10)= 4.592886537931051e-008
91
92	return
93	end
94	subroutine flowpara(ifl,jtu)
95	implicit none
96	integer iflow,jtune
97	common /vecflow/ iflow,jtune
98	integer ifl,jtu
99	iflow=ifl
100	jtune=jtu
101	return
102	end
103	subroutine pertpeek(st,ang,ra)
104	implicit none
105	integer i,ndim,ndim2,nreso,ntt
106	double precision ang,ra,st
107	parameter (ndim=3)
108	parameter (ndim2=6)
109	parameter (ntt=40)
110	parameter (nreso=20)
111	dimension st(ndim),ang(ndim),ra(ndim)
112	double precision angle,dsta,rad,sta
113	common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
114	integer idsta,ista
115	common /istable/ista(ndim),idsta(ndim)
116	integer nd,nd2,no,nv
117	common /ii/no,nv,nd,nd2
118	integer ndc,ndc2,ndpt,ndt
119	common /coast/ndc,ndc2,ndt,ndpt
120	integer mx,nres
121	common /reson/mx(ndim,nreso),nres
122	integer iref
123	common /resfile/iref
124	do i=1,nd
125	st(i)=sta(i)
126	ang(i)=angle(i)
127	ra(i)=rad(i)
128	enddo
129	return
130	end
131	subroutine inputres(mx1,nres1)
132	implicit none
133	integer i,j,ndim,ndim2,nreso,ntt
134	parameter (ndim=3)
135	parameter (ndim2=6)
136	parameter (ntt=40)
137	parameter (nreso=20)
138	integer mx1(ndim,nreso),nres1
139	integer mx,nres
140	common /reson/mx(ndim,nreso),nres
141
142	nres=nres1
143	do i=1,nreso
144	do j=1,ndim
145	mx(j,i)=0
146	enddo
147	enddo
148
149	do i=1,nres
150	do j=1,ndim
151	mx(j,i)=mx1(j,i)
152	enddo
153	enddo
154	return
155	end
156	subroutine respoke(mres,nre,ire)
157	implicit none
158	integer i,ire,j,ndim,ndim2,nre,nreso,ntt
159	double precision ang,ra,st
160	parameter (ndim=3)
161	parameter (ndim2=6)
162	parameter (ntt=40)
163	parameter (nreso=20)
164	integer mres(ndim,nreso)
165	double precision angle,dsta,rad,sta
166	common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
167	integer idsta,ista
168	common /istable/ista(ndim),idsta(ndim)
169	integer nd,nd2,no,nv
170	common /ii/no,nv,nd,nd2
171	integer ndc,ndc2,ndpt,ndt
172	common /coast/ndc,ndc2,ndt,ndpt
173	integer mx,nres
174	common /reson/mx(ndim,nreso),nres
175	integer iref
176	common /resfile/iref
177	dimension ang(ndim),ra(ndim),st(ndim)
178	iref=ire
179	nres=nre
180	do j=1,nreso
181	do i=1,nd
182	mx(i,j)=mres(i,j)
183	enddo
184	enddo
185	call initpert(st,ang,ra)
186	return
187	end
188	subroutine liepeek(iia,icoast)
189	implicit none
190	integer ndim,ndim2,nreso,ntt
191	parameter (ndim=3)
192	parameter (ndim2=6)
193	parameter (ntt=40)
194	parameter (nreso=20)
195	integer nd,nd2,no,nv
196	common /ii/no,nv,nd,nd2
197	integer ndc,ndc2,ndpt,ndt
198	common /coast/ndc,ndc2,ndt,ndpt
199	integer iia(),icoast()
200
201	iia(1)=no
202	iia(2)=nv
203	iia(3)=nd
204	iia(4)=nd2
205	icoast(1)=ndc
206	icoast(2)=ndc2
207	icoast(3)=ndt
208	icoast(4)=ndpt
209
210	return
211	end
212	subroutine lienot(not)
213	implicit none
214	integer no,not
215
216	call danot(not)
217	no=not
218
219	return
220	end
221	subroutine etallnom(x,n,nom)
222	implicit none
223	integer i,n,nd2
224	! CREATES A AD-VARIABLE WHICH CAN BE DESTROYED BY DADAL
225	! allocates vector of n polynomials and give it the name NOM=A10
226	integer x(*),i1(4),i2(4)
227	character*10 nom
228	do i=1,iabs(n)
229	x(i)=0
230	enddo
231	call daallno(x,iabs(n),nom)
232	if(n.lt.0) then
233	call liepeek(i1,i2)
234	nd2=i1(4)
235	do i=nd2+1,-n
236	call davar(x(i),0.d0,i)
237	enddo
238	endif
239	return
240	end
241	subroutine etall(x,n)
242	implicit none
243	integer i,n,nd2
244	! allocates vector of n polynomials
245	integer x(*),i1(4),i2(4)
246	do i=1,iabs(n)
247	x(i)=0
248	enddo
249	call daallno(x,iabs(n),'ETALL ')
250	if(n.lt.0) then
251	call liepeek(i1,i2)
252	nd2=i1(4)
253	do i=nd2+1,-n
254	call davar(x(i),0.d0,i)
255	enddo
256	endif
257	return
258	end
259	subroutine etall1(x)
260	implicit none
261	integer x
262	call daallno(x,1,'ETALL ')
263	return
264	end
265	subroutine dadal1(x)
266	implicit none
267	integer x
268	call dadal(x,1)
269	return
270	end
271	subroutine etppulnv(x,xi,xff)
272	implicit none
273	integer i,ndim,ndim2,ntt
274	parameter (ndim=3)
275	parameter (ndim2=6)
276	parameter (ntt=40)
277	integer nd,nd2,no,nv
278	common /ii/no,nv,nd,nd2
279	integer x(*)
280	double precision xi(),xff(),xf(ntt),xii(ntt)
281
282	do i=1,nv
283	xii(i)=xi(i)
284	enddo
285	do i=nv+1,ntt
286	xii(i)=0.d0
287	enddo
288
289	call ppush(x,nv,xii,xf)
290
291	do i=1,nv
292	xff(i)=xf(i)
293	enddo
294
295	return
296	end
297	subroutine etmtree(y,x)
298	implicit none
299	integer i,ie,iv,ndim,ndim2,nt,ntt
300	! ROUTINES USING THE MAP IN AD-FORM
301	parameter (ndim=3)
302	parameter (ndim2=6)
303	parameter (ntt=40)
304	dimension ie(ntt),iv(ntt)
305	integer nd,nd2,no,nv
306	common /ii/no,nv,nd,nd2
307	integer x(),y()
308
309	nt=nv-nd2
310	if(nt.gt.0) then
311	call etallnom(ie,nt,'IE ')
312	do i=nd2+1,nv
313	call davar(ie(i-nd2),0.d0,i)
314	enddo
315	do i=nd2+1,nv
316	iv(i)=ie(i-nd2)
317	enddo
318	endif
319	do i=1,nd2
320	iv(i)=y(i)
321	enddo
322	call mtree(iv,nv,x,nv)
323	if(nt.gt.0) then
324	call dadal(ie,nt)
325	endif
326	return
327	end
328	subroutine etppush(x,xi)
329	implicit none
330	integer i,ndim,ndim2,ntt
331	parameter (ndim=3)
332	parameter (ndim2=6)
333	parameter (ntt=40)
334	integer nd,nd2,no,nv
335	common /ii/no,nv,nd,nd2
336	integer x(*)
337	double precision xi(*),xf(ntt),xii(ntt)
338
339	do i=1,nd2
340	xii(i)=xi(i)
341	enddo
342
343	call ppush(x,nv,xii,xf)
344
345	do i=1,nd2
346	xi(i)=xf(i)
347	enddo
348
349	return
350	end
351	subroutine etppush2(x,xi,xff)
352	implicit none
353	integer i,ndim,ndim2,ntt
354	parameter (ndim=3)
355	parameter (ndim2=6)
356	parameter (ntt=40)
357	integer nd,nd2,no,nv
358	common /ii/no,nv,nd,nd2
359	integer x(*)
360	double precision xi(),xff(),xf(ntt),xii(ntt)
361
362	do i=1,nd2
363	xii(i)=xi(i)
364	enddo
365
366	call ppush(x,nv,xii,xf)
367
368	do i=1,nd2
369	xff(i)=xf(i)
370	enddo
371
372	return
373	end
374	subroutine ppushlnv(x,xi,xff,nd1)
375	implicit none
376	integer i,nd1,ndim,ndim2,ntt
377	parameter (ndim=3)
378	parameter (ndim2=6)
379	parameter (ntt=40)
380	integer nd,nd2,no,nv
381	common /ii/no,nv,nd,nd2
382	integer x(*)
383	double precision xi(),xff(),xf(ntt),xii(ntt)
384
385	do i=1,nd1
386	xii(i)=xi(i)
387	enddo
388	do i=nd1+1,ntt
389	xii(i)=0.d0
390	enddo
391
392	call ppush(x,nv,xii,xf)
393
394	do i=1,nd1
395	xff(i)=xf(i)
396	enddo
397
398	return
399	end
400	subroutine etcct(x,y,z)
401	implicit none
402	integer i,ie,iv,ndim,ndim2,nt,ntt
403	! Z=XoY
404	parameter (ndim=3)
405	parameter (ndim2=6)
406	parameter (ntt=40)
407	dimension ie(ntt),iv(ntt)
408	integer nd,nd2,no,nv
409	common /ii/no,nv,nd,nd2
410	integer x(),y(),z(*)
411
412	nt=nv-nd2
413	if(nt.gt.0) then
414	call etallnom(ie,nt,'IE ')
415	do i=nd2+1,nv
416	call davar(ie(i-nd2),0.d0,i)
417	enddo
418	do i=nd2+1,nv
419	iv(i)=ie(i-nd2)
420	enddo
421	endif
422	do i=1,nd2
423	iv(i)=y(i)
424	enddo
425	call dacct(x,nd2,iv,nv,z,nd2)
426	if(nt.gt.0) then
427	call dadal(ie,nt)
428	endif
429	return
430	end
431	subroutine trx(h,rh,y)
432	implicit none
433	integer i,ie,iv,ndim,ndim2,nt,ntt
434	! :RH: = Y :H: Y^-1 = :HoY:
435	parameter (ndim=3)
436	parameter (ndim2=6)
437	parameter (ntt=40)
438	dimension ie(ntt),iv(ntt)
439	integer nd,nd2,no,nv
440	common /ii/no,nv,nd,nd2
441	integer h,rh
442	integer y(*)
443	!
444	nt=nv-nd2
445	if(nt.gt.0) then
446	call etallnom(ie,nt,'IE ')
447	do i=nd2+1,nv
448	call davar(ie(i-nd2),0.d0,i)
449	enddo
450	do i=nd2+1,nv
451	iv(i)=ie(i-nd2)
452	enddo
453	endif
454	do i=1,nd2
455	iv(i)=y(i)
456	enddo
457	call dacct(h,1,iv,nv,rh,1)
458	if(nt.gt.0) then
459	call dadal(ie,nt)
460	endif
461	return
462	end
463	subroutine trxflo(h,rh,y)
464	implicit none
465	integer j,k,ndim,ndim2,ntt
466	! RH = Y H Y^-1 CHANGE OF A VECTOR FLOW OPERATOR
467	parameter (ndim=3)
468	parameter (ndim2=6)
469	parameter (ntt=40)
470	integer nd,nd2,no,nv
471	common /ii/no,nv,nd,nd2
472	integer h(),rh(),y(*)
473	integer yi(ndim2),ht(ndim2),b1,b2
474	!
475	!
476	call etallnom(yi,nd2 ,'YI ')
477	call etallnom(ht,nd2 ,'HT ')
478	call etallnom(b1,1,'B1 ')
479	call etallnom(b2,1,'B2 ')
480
481	call etinv(y,yi)
482	!----- HT= H o Y
483	call etcct(h,y,ht)
484	!----
485	call daclrd(rh)
486	do j=1,nd2
487	do k=1,nd2
488	call dader(k,yi(j),b1)
489	call trx(b1,b2,y)
490	call damul(b2,ht(k),b1)
491	call daadd(b1,rh(j),b2)
492	call dacop(b2,rh(j))
493	enddo
494	enddo
495
496	call dadal(b2,1)
497	call dadal(b1,1)
498	call dadal(ht,nd2)
499	call dadal(yi,nd2)
500	return
501	end
502	subroutine simil(a,x,ai,y)
503	implicit none
504	integer ndim,ndim2,ntt
505	! Y= AoXoAI
506	parameter (ndim=3)
507	parameter (ndim2=6)
508	parameter (ntt=40)
509	integer nd,nd2,no,nv
510	common /ii/no,nv,nd,nd2
511	integer x(),y(),a(),ai()
512
513	integer w(ndim2),v(ndim2)
514	!
515	call etallnom(w,nd2 ,'W ')
516	call etallnom(v,nd2 ,'V ')
517
518	call etcct(a,x,w)
519	call etcct(w,ai,v)
520
521	call dacopd(v,y)
522
523	call dadal(v,nd2)
524	call dadal(w,nd2)
525	return
526	end
527	subroutine etini(x)
528	implicit none
529	integer i,ndim,ndim2,ntt
530	! X=IDENTITY
531	parameter (ndim=3)
532	parameter (ndim2=6)
533	parameter (ntt=40)
534	integer nd,nd2,no,nv
535	common /ii/no,nv,nd,nd2
536	integer x(*)
537	!*DAEXT(NO,NV) X(NDIM2)
538	do i=1,nd2
539	call davar(x(i),0.d0,i)
540	enddo
541	return
542	end
543	subroutine etinv(x,y)
544	implicit none
545	integer i,ie1,ie2,iv1,iv2,ndim,ndim2,nt,ntt
546	! Y=X^-1
547	parameter (ndim=3)
548	parameter (ndim2=6)
549	parameter (ntt=40)
550	dimension ie1(ntt),ie2(ntt),iv1(ntt),iv2(ntt)
551	integer nd,nd2,no,nv
552	common /ii/no,nv,nd,nd2
553
554	integer x(),y()
555
556	nt=nv-nd2
557	if(nt.gt.0) then
558	do i=1,nt
559	ie1(i)=0
560	ie2(i)=0
561	enddo
562	call etallnom(ie1,nt,'IE1 ')
563	call etallnom(ie2,nt,'IE2 ')
564	do i=nd2+1,nv
565	call davar(ie1(i-nd2),0.d0,i)
566	enddo
567	do i=nd2+1,nv
568	iv1(i)=ie1(i-nd2)
569	iv2(i)=ie2(i-nd2)
570	enddo
571	endif
572	do i=1,nd2
573	iv1(i)=x(i)
574	iv2(i)=y(i)
575	enddo
576
577	call dainv(iv1,nv,iv2,nv)
578	if(nt.gt.0) then
579	call dadal(ie2,nt)
580	call dadal(ie1,nt)
581	endif
582	return
583	end
584	subroutine etpin(x,y,jj)
585	implicit none
586	integer i,ie1,ie2,iv1,iv2,jj,ndim,ndim2,nt,ntt
587	! Y=PARTIAL INVERSION OF X SEE BERZ'S PACKAGE
588	parameter (ndim=3)
589	parameter (ndim2=6)
590	parameter (ntt=40)
591	dimension ie1(ntt),ie2(ntt),iv1(ntt),iv2(ntt),jj(*)
592	integer nd,nd2,no,nv
593	common /ii/no,nv,nd,nd2
594
595	integer x(),y()
596
597	nt=nv-nd2
598	if(nt.gt.0) then
599	do i=1,nt
600	ie1(i)=0
601	ie2(i)=0
602	enddo
603	call etallnom(ie1,nt,'IE1 ')
604	call etallnom(ie2,nt,'IE2 ')
605	do i=nd2+1,nv
606	call davar(ie1(i-nd2),0.d0,i)
607	enddo
608	do i=nd2+1,nv
609	iv1(i)=ie1(i-nd2)
610	iv2(i)=ie2(i-nd2)
611	enddo
612	endif
613	do i=1,nd2
614	iv1(i)=x(i)
615	iv2(i)=y(i)
616	enddo
617
618	call dapin(iv1,nv,iv2,nv,jj)
619	if(nt.gt.0) then
620	call dadal(ie2,nt)
621	call dadal(ie1,nt)
622	endif
623	return
624	end
625	subroutine dapek0(v,x,jj)
626	implicit none
627	integer i,jj,ndim2,ntt
628	double precision x
629	!- MORE EXTENSIONS OF BASIC BERZ'S PACKAGE
630	parameter (ndim2=6)
631	parameter (ntt=40)
632	integer v(*),jd(ntt)
633	dimension x(*)
634	do i=1,ntt
635	jd(i)=0
636	enddo
637	do i=1,jj
638	call dapek(v(i),jd,x(i))
639	enddo
640	return
641	end
642	subroutine dapok0(v,x,jj)
643	implicit none
644	integer i,jj,ndim2,ntt
645	double precision x
646	parameter (ndim2=6)
647	parameter (ntt=40)
648	integer v(*),jd(ntt)
649	dimension x(*)
650	do i=1,ntt
651	jd(i)=0
652	enddo
653	do i=1,jj
654	call dapok(v(i),jd,x(i))
655	enddo
656	return
657	end
658	subroutine dapokzer(v,jj)
659	implicit none
660	integer i,jj,ndim2,ntt
661	parameter (ndim2=6)
662	parameter (ntt=40)
663	integer v(*),jd(ntt)
664	do i=1,ntt
665	jd(i)=0
666	enddo
667	do i=1,jj
668	call dapok(v(i),jd,0.d0)
669	enddo
670	return
671	end
672	subroutine davar0(v,x,jj)
673	implicit none
674	integer i,jj,ndim2,ntt
675	double precision x
676	parameter (ndim2=6)
677	parameter (ntt=40)
678	integer v(*)
679	dimension x(*)
680	do i=1,jj
681	call davar(v(i),x(i),i)
682	enddo
683	return
684	end
685	subroutine comcfu(b,f1,f2,c)
686	implicit none
687	double precision f1,f2
688	external f1,f2
689	! Complex dacfu
690	integer b(),c(),t(4)
691	call etall(t,4)
692
693	call dacfu(b(1),f1,t(1))
694	call dacfu(b(1),f2,t(2))
695	call dacfu(b(2),f1,t(3))
696	call dacfu(b(2),f2,t(4))
697
698	call dasub(t(1),t(4),c(1))
699	call daadd(t(2),t(3),c(2))
700	call dadal(t,4)
701	return
702	end
703	subroutine take(h,m,ht)
704	implicit none
705	integer i,m,ndim,ntt
706	double precision r
707	! HT= H_M (TAKES M^th DEGREE PIECE ALL VARIABLES INCLUDED)
708	parameter (ndim=3)
709	parameter (ntt=40)
710	integer nd,nd2,no,nv
711	common /ii/no,nv,nd,nd2
712	integer h,ht,j(ntt)
713
714	integer b1,b2,b3
715	!
716	!
717	call etallnom(b1,1,'B1 ')
718	call etallnom(b2,1,'B2 ')
719	call etallnom(b3,1,'B3 ')
720
721	if(no.ge.2) then
722	if(m.eq.0) then
723	do i=1,ntt
724	j(i)=0
725	enddo
726	call dapek(h,j,r)
727	call dacon(ht,r)
728	else
729	call danot(m)
730	call dacop(h,b1)
731	call danot(m-1)
732	call dacop(b1,b2)
733	call danot(no)
734	call dasub(b1,b2,b3)
735	call dacop(b3,ht)
736	endif
737	else
738	do i=1,ntt
739	j(i)=0
740	enddo
741	if(m.eq.0) then
742	call dapek(h,j,r)
743	call dacon(ht,r)
744	elseif(m.eq.1) then
745	do i=1,nv
746	j(i)=1
747	call dapek(h,j,r)
748	call dapok(b3,j,r)
749	j(i)=0
750	enddo
751	call dacop(b3,ht)
752	else
753	call daclr(ht)
754	endif
755	endif
756
757	call dadal(b3,1)
758	call dadal(b2,1)
759	call dadal(b1,1)
760	return
761	end
762	subroutine taked(h,m,ht)
763	implicit none
764	integer i,m,ndim2,ntt
765	! \VEC{HT}= \VEC{H_M} (TAKES M^th DEGREE PIECE ALL VARIABLES INCLUDED)
766	parameter (ndim2=6)
767	parameter (ntt=40)
768	integer nd,nd2,no,nv
769	common /ii/no,nv,nd,nd2
770	integer h(),ht(),j(ntt)
771
772	integer b1,b2,x(ndim2)
773	!
774	call etallnom(b1,1,'B1 ')
775	call etallnom(b2,1,'B2 ')
776	call etallnom(x,nd2 ,'X ')
777
778
779	do i=1,ntt
780	j(i)=0
781	enddo
782
783	do i=1,nd2
784	call take(h(i),m,ht(i))
785	enddo
786	call dadal(x,nd2)
787	call dadal(b2,1)
788	call dadal(b1,1)
789	return
790	end
791	subroutine daclrd(h)
792	implicit none
793	integer i,ndim2,ntt
794	! clear a map : a vector of nd2 polynomials
795	parameter (ndim2=6)
796	parameter (ntt=40)
797	integer nd,nd2,no,nv
798	common /ii/no,nv,nd,nd2
799	integer h(*)
800	do i=1,nd2
801	call daclr(h(i))
802	enddo
803	return
804	end
805	subroutine dacopd(h,ht)
806	implicit none
807	integer i,ndim2,ntt
808	! H goes into HT (nd2 array)
809	parameter (ndim2=6)
810	parameter (ntt=40)
811	integer nd,nd2,no,nv
812	common /ii/no,nv,nd,nd2
813	integer h(),ht()
814	do i=1,nd2
815	call dacop(h(i),ht(i))
816	enddo
817	return
818	end
819	subroutine dacmud(h,sca,ht)
820	implicit none
821	integer i,ndim2,ntt
822	double precision sca
823	parameter (ndim2=6)
824	parameter (ntt=40)
825	integer nd,nd2,no,nv
826	common /ii/no,nv,nd,nd2
827	integer h(),ht()
828	do i=1,nd2
829	call dacmu(h(i),sca,ht(i))
830	enddo
831	return
832	end
833	subroutine dalind(h,rh,ht,rt,hr)
834	implicit none
835	integer i,ndim2
836	double precision rh,rt
837	parameter (ndim2=6)
838	integer nd,nd2,no,nv
839	common /ii/no,nv,nd,nd2
840	integer h(),ht(),hr(*)
841
842	integer b(ndim2)
843	!
844	call etallnom(b,nd2 ,'B ')
845
846	do i=1,nd2
847	call dalin(h(i),rh,ht(i),rt,b(i))
848	enddo
849	call dacopd(b,hr)
850	call dadal(b,nd2)
851	return
852	end
853	subroutine daread(h,nd1,mfile,xipo)
854	implicit none
855	integer i,mfile,nd1,ndim2,ntt
856	double precision rx,xipo
857	! read a map
858	parameter (ndim2=6)
859	parameter (ntt=40)
860	integer h(*),j(ntt)
861	do i=1,ntt
862	j(i)=0
863	enddo
864	do i=1,nd1
865	call darea(h(i),mfile)
866	call dapek(h(i),j,rx)
867	rx=rx*xipo
868	call dapok(h(i),j,rx)
869	enddo
870	return
871	end
872	subroutine daprid(h,n1,n2,mfile)
873	implicit none
874	integer i,mfile,n1,n2,ndim2,ntt
875	! print a map
876	parameter (ndim2=6)
877	parameter (ntt=40)
878	integer h(*)
879	if(mfile.le.0) return
880	do i=n1,n2
881	call dapri(h(i),mfile)
882	enddo
883	return
884	end
885	subroutine prresflo(h,eps,mfile)
886	implicit none
887	integer i,mfile,ndim2,ntt
888	double precision deps,eps,filtres
889	! print a map in resonance basis for human consumption (useless)
890	parameter (ndim2=6)
891	parameter (ntt=40)
892	integer h(*) ,b(ndim2) ,c(ndim2)
893	integer nd,nd2,no,nv
894	common /ii/no,nv,nd,nd2
895	integer ifilt
896	common /filtr/ ifilt
897	external filtres
898	call etall(b,nd2)
899	call etall(c,nd2)
900	call dacopd(h,c)
901	do i=1,nd2
902	ifilt=(-1)**i
903	call dacfu(c(i),filtres,h(i))
904	enddo
905
906	deps=-1.d0
907	call daeps(deps)
908	call daeps(eps)
909
910	call dacopd(c,h)
911	call daeps(deps)
912	call dadal(c,nd2)
913	call dadal(b,nd2)
914	return
915	end
916	double precision function filtres(j)
917	implicit none
918	integer i,ic,ndim
919	parameter (ndim=3)
920	! PARAMETER (NTT=40)
921	! INTEGER J(NTT)
922	integer j(*)
923	integer nd,nd2,no,nv
924	common /ii/no,nv,nd,nd2
925	integer ndc,ndc2,ndpt,ndt
926	common /coast/ndc,ndc2,ndt,ndpt
927	integer ifilt
928	common /filtr/ ifilt
929	filtres=1.d0
930	ic=0
931	do i=1,(nd2-ndc2)
932	ic=ic+j(i)(-1)*(i+1)
933	enddo
934	ic=ic+ifilt
935	if(ic.lt.0) filtres=0.d0
936	if(ic.eq.0.and.ifilt.eq.1) then
937	filtres=0.0d0
938	endif
939	return
940	end
941	subroutine daflo(h,x,y)
942	implicit none
943	integer i,ndim,ndim2,ntt
944	! LIE EXPONENT ROUTINES WITH FLOW OPERATORS
945
946	! \VEC{H}.GRAD X =Y
947	parameter (ndim=3)
948	parameter (ndim2=6)
949	parameter (ntt=40)
950	integer nd,nd2,no,nv
951	common /ii/no,nv,nd,nd2
952	integer h(*),x,y
953	integer b1,b2,b3
954	!
955	call etallnom(b1,1,'B1 ')
956	call etallnom(b2,1,'B2 ')
957	call etallnom(b3,1,'B3 ')
958
959	call daclr(b1)
960	call daclr(b2)
961	do i=1,nd2
962	call dader(i,x,b2)
963	call damul(b2,h(i),b3)
964	call daadd(b3,b1,b2)
965	call dacop(b2,b1)
966	enddo
967	call dacop(b1,y)
968	call dadal(b3,1)
969	call dadal(b2,1)
970	call dadal(b1,1)
971	return
972	end
973	subroutine daflod(h,x,y)
974	implicit none
975	integer i,ndim,ndim2,ntt
976	parameter (ndim=3)
977	parameter (ndim2=6)
978	parameter (ntt=40)
979	integer nd,nd2,no,nv
980	common /ii/no,nv,nd,nd2
981	integer h(),x(),y(*)
982	integer b1(ndim2),b2(ndim2)
983	!
984	call etall(b1,nd2)
985	call etall(b2,nd2)
986
987	call dacopd(h,b1)
988	call dacopd(x,b2)
989
990	do i=1,nd2
991	call daflo(b1,b2(i),y(i))
992	enddo
993
994	call dadal(b1,nd2)
995	call dadal(b2,nd2)
996	return
997	end
998	subroutine intd(v,h,sca)
999	implicit none
1000	integer i,ndim,ndim2,ntt
1001	double precision dlie,sca
1002	! IF SCA=-1.D0
1003	! \VEC{V}.GRAD = J GRAD H . GRAD = :H:
1004
1005	! IF SCA=1.D0
1006	! \VEC{V}.GRAD = GRAD H . GRAD
1007	parameter (ndim=3)
1008	parameter (ndim2=6)
1009	parameter (ntt=40)
1010	integer nd,nd2,no,nv
1011	common /ii/no,nv,nd,nd2
1012	external dlie
1013	integer v(*),h
1014
1015	integer b1,b2,b3,b4,x(ndim2)
1016	!
1017	!
1018	call etallnom(b1,1,'B1 ')
1019	call etallnom(b2,1,'B2 ')
1020	call etallnom(b3,1,'B3 ')
1021	call etallnom(b4,1,'B4 ')
1022	call etallnom(x,nd2 ,'X ')
1023
1024	call daclr(b4)
1025	call daclr(h)
1026	call etini(x)
1027	do i=1,nd
1028	call dacfu(v(2*i-1),dlie,b3)
1029	call dacfu(v(2*i),dlie,b1)
1030	call damul(b1,x(2*i-1),b2)
1031	call damul(b3,x(2*i),b1)
1032	call dalin(b2,1.d0,b1,sca,b3)
1033	call daadd(b3,b4,b2)
1034	call dacop(b2,b4)
1035	enddo
1036	call dacop(b4,h)
1037	call dadal(x,nd2)
1038	call dadal(b4,1)
1039	call dadal(b3,1)
1040	call dadal(b2,1)
1041	call dadal(b1,1)
1042	return
1043	end
1044	subroutine difd(h1,v,sca)
1045	implicit none
1046	integer i,ndim,ndim2,ntt
1047	double precision sca
1048	! INVERSE OF INTD ROUTINE
1049	parameter (ndim=3)
1050	parameter (ndim2=6)
1051	parameter (ntt=40)
1052	integer nd,nd2,no,nv
1053	common /ii/no,nv,nd,nd2
1054	integer v(*),h1
1055	integer b1,h
1056	call etall(b1,1)
1057	call etall(h,1)
1058	call dacop(h1,h)
1059	do i=1,nd
1060	call dader(2i-1,h,v(2i))
1061	call dader(2*i,h,b1)
1062	call dacmu(b1,sca,v(2*i-1))
1063	enddo
1064	call dadal(h,1)
1065	call dadal(b1,1)
1066	return
1067	end
1068	subroutine expflo(h,x,y,eps,nrmax)
1069	implicit none
1070	integer i,ndim,ndim2,nrmax,ntt
1071	double precision coe,eps,r,rbefore
1072	! DOES EXP( \VEC{H} ) X = Y
1073	parameter (ndim=3)
1074	parameter (ndim2=6)
1075	parameter (ntt=40)
1076	integer idpr
1077	common /printing/ idpr
1078	integer h(*),x,y
1079	integer b1,b2,b3,b4
1080	logical more
1081	!
1082	!
1083	call etallnom(b1,1,'B1 ')
1084	call etallnom(b2,1,'B2 ')
1085	call etallnom(b3,1,'B3 ')
1086	call etallnom(b4,1,'B4 ')
1087
1088	call dacop(x,b4)
1089	call dacop(x,b1)
1090	more=.true.
1091	rbefore=1.d30
1092	do i=1,nrmax
1093	coe=1.d0/dble(i)
1094	call dacmu(b1,coe,b2)
1095	call daflo(h,b2,b1)
1096	call daadd(b4,b1,b3)
1097	call daabs(b1,r)
1098	if(more) then
1099	if(r.gt.eps) then
1100	rbefore=r
1101	goto 100
1102	else
1103	rbefore=r
1104	more=.false.
1105	endif
1106	else
1107	if(r.ge.rbefore) then
1108	call dacop(b3,y)
1109	call dadal(b4,1)
1110	call dadal(b3,1)
1111	call dadal(b2,1)
1112	call dadal(b1,1)
1113	return
1114	endif
1115	rbefore=r
1116	endif
1117	100 continue
1118	call dacop(b3,b4)
1119	enddo
1120	if(idpr.ge.0) then
1121	write(6,*) ' NORM ',eps,' NEVER REACHED IN EXPFLO '
1122	write(6,*) 'NEW IDPR '
1123	read(5,*) idpr
1124	endif
1125	call dacop(b3,y)
1126	call dadal(b4,1)
1127	call dadal(b3,1)
1128	call dadal(b2,1)
1129	call dadal(b1,1)
1130	return
1131	end
1132	subroutine expflod(h,x,w,eps,nrmax)
1133	implicit none
1134	integer j,ndim,ndim2,nrmax,ntt
1135	double precision eps
1136	! DOES EXP( \VEC{H} ) \VEC{X} = \VEC{Y}
1137	parameter (ndim=3)
1138	parameter (ndim2=6)
1139	parameter (ntt=40)
1140	integer nd,nd2,no,nv
1141	common /ii/no,nv,nd,nd2
1142	integer x(),w(),h(*)
1143	integer b0,v(ndim2)
1144	!
1145	!
1146	call etallnom(b0,1,'B0 ')
1147	call etallnom(v,nd2 ,'V ')
1148
1149	call dacopd(x,v)
1150	do j=1,nd2
1151	call expflo(h,v(j),b0,eps,nrmax)
1152	call dacop(b0,v(j))
1153	enddo
1154	call dacopd(v,w)
1155	call dadal(v,nd2)
1156	call dadal(b0,1)
1157	return
1158	end
1159	subroutine facflo(h,x,w,nrmin,nrmax,sca,ifac)
1160	implicit none
1161	integer i,ifac,ndim,ndim2,nmax,nrmax,nrmin,ntt
1162	double precision eps,sca
1163	! IFAC=1
1164	! DOES EXP(SCA \VEC{H}_MRMIN ) ... EXP(SCA \VEC{H}_NRMAX ) X= Y
1165	! IFAC=-1
1166	! DOES EXP(SCA \VEC{H}_NRMAX ) ... EXP(SCA \VEC{H}_MRMIN ) X= Y
1167	parameter (ndim=3)
1168	parameter (ndim2=6)
1169	parameter (ntt=40)
1170	integer nd,nd2,no,nv
1171	common /ii/no,nv,nd,nd2
1172	integer x,w,h(*)
1173	integer bm(ndim2),b0(ndim2),v
1174	!
1175	call etallnom(bm,nd2 ,'BM ')
1176	call etallnom(b0,nd2 ,'B0 ')
1177	call etallnom(v,1 ,'V ')
1178
1179	call dacop(x,v)
1180
1181	eps=-1.d0
1182	call daeps(eps)
1183	nmax=100
1184	!
1185	! IFAC =1 ---> V = EXP(:SCAH(NRMAX):)...EXP(:SCAH(NRMIN):)X
1186	if(ifac.eq.1) then
1187	do i=nrmax,nrmin,-1
1188	call taked(h,i,b0)
1189	call dacmud(b0,sca,bm)
1190
1191	call expflo(bm,v,b0(1),eps,nmax)
1192	call dacop(b0(1),v)
1193	enddo
1194	else
1195	! IFAC =-1 ---> V = EXP(:SCAH(NRMIN):)...EXP(:SCAH(NRMAX):)X
1196	do i=nrmin,nrmax
1197	call taked(h,i,b0)
1198	call dacmud(b0,sca,bm)
1199
1200	call expflo(bm,v,b0(1),eps,nmax)
1201	call dacop(b0(1),v)
1202	enddo
1203	endif
1204	call dacop(v,w)
1205	call dadal(v,1)
1206	call dadal(b0,nd2)
1207	call dadal(bm,nd2)
1208	return
1209	end
1210	subroutine facflod(h,x,w,nrmin,nrmax,sca,ifac)
1211	implicit none
1212	integer i,ifac,ndim,ndim2,nrmax,nrmin,ntt
1213	double precision sca
1214	! IFAC=1
1215	! DOES EXP(SCA \VEC{H}_MRMIN ) ... EXP(SCA \VEC{H}_NRMAX ) \VEC{X}= \VEC{Y}
1216	! IFAC=-1
1217	! DOES EXP(SCA \VEC{H}_NRMAX ) ... EXP(SCA \VEC{H}_MRMIN ) \VEC{X}= \VEC{Y}
1218	parameter (ndim=3)
1219	parameter (ndim2=6)
1220	parameter (ntt=40)
1221	integer nd,nd2,no,nv
1222	common /ii/no,nv,nd,nd2
1223	integer x(),w(),h(*)
1224
1225	do i=1,nd2
1226	call facflo(h,x(i),w(i),nrmin,nrmax,sca,ifac)
1227	enddo
1228
1229	return
1230	end
1231	subroutine fexpo(h,x,w,nrmin,nrmax,sca,ifac)
1232	implicit none
1233	integer ifac,ndim,ndim2,nrma,nrmax,nrmi,nrmin,ntt
1234	double precision sca
1235	! WRAPPED ROUTINES FOR THE OPERATOR \VEC{H}=:H:
1236	! WRAPPING FACFLOD
1237	parameter (ndim=3)
1238	parameter (ndim2=6)
1239	parameter (ntt=40)
1240	integer nd,nd2,no,nv
1241	common /ii/no,nv,nd,nd2
1242	integer x(),w(),h
1243
1244	integer v(ndim2)
1245
1246	nrmi=nrmin-1
1247	nrma=nrmax-1
1248	call etall(v,nd2)
1249	call difd(h,v,-1.d0)
1250	call facflod(v,x,w,nrmi,nrma,sca,ifac)
1251
1252	call dadal(v,nd2)
1253
1254	return
1255	end
1256	subroutine etcom(x,y,h)
1257	implicit none
1258	integer i,j,ndim,ndim2,ntt
1259	! ETCOM TAKES THE BRACKET OF TWO VECTOR FIELDS.
1260	parameter (ndim2=6)
1261	parameter (ndim=3)
1262	parameter (ntt=40)
1263	integer nd,nd2,no,nv
1264	common /ii/no,nv,nd,nd2
1265	integer h(),x(),y(*),t1,t2,t3(ndim2)
1266
1267	call etall(t1,1)
1268	call etall(t2,1)
1269	call etall(t3,nd2)
1270
1271	do j=1,nd2
1272	do i=1,nd2
1273
1274	call dader(i,x(j),t1)
1275	call dader(i,y(j),t2)
1276	call damul(x(i),t2,t2)
1277	call damul(y(i),t1,t1)
1278	call dalin(t2,1.d0,t1,-1.d0,t1)
1279	call daadd(t1,t3(j),t3(j))
1280
1281	enddo
1282	enddo
1283
1284	call dacopd(t3,h)
1285
1286	call dadal(t1,1)
1287	call dadal(t2,1)
1288	call dadal(t3,nd2)
1289	return
1290	end
1291	subroutine etpoi(x,y,h)
1292	implicit none
1293	integer i,ndim,ndim2,ntt
1294	! ETPOI TAKES THE POISSON BRACKET OF TWO FUNCTIONS
1295	parameter (ndim2=6)
1296	parameter (ndim=3)
1297	parameter (ntt=40)
1298	integer nd,nd2,no,nv
1299	common /ii/no,nv,nd,nd2
1300	integer h,x,y,t1,t2,t3
1301
1302	call etall(t1,1)
1303	call etall(t2,1)
1304	call etall(t3,1)
1305
1306	do i=1,nd
1307
1308	call dader(2*i-1,x,t1)
1309	call dader(2*i,y,t2)
1310	call damul(t1,t2,t1)
1311
1312	call dalin(t1,1.d0,t3,1.d0,t3)
1313	call dader(2*i-1,y,t1)
1314	call dader(2*i,x,t2)
1315	call damul(t1,t2,t1)
1316
1317	call dalin(t1,-1.d0,t3,1.d0,t3)
1318
1319	enddo
1320
1321	call dacop(t3,h)
1322
1323	call dadal(t1,1)
1324	call dadal(t2,1)
1325	call dadal(t3,1)
1326	return
1327	end
1328	subroutine exp1d(h,x,y,eps,non)
1329	implicit none
1330	integer ndim,ndim2,non,ntt
1331	double precision eps
1332	! WRAPPING EXPFLO
1333	parameter (ndim2=6)
1334	parameter (ndim=3)
1335	parameter (ntt=40)
1336	integer nd,nd2,no,nv
1337	common /ii/no,nv,nd,nd2
1338	integer idpr
1339	common /printing/ idpr
1340	integer h,x,y
1341
1342	integer v(ndim2)
1343
1344	call etall(v,nd2)
1345	call difd(h,v,-1.d0)
1346	call expflo(v,x,y,eps,non)
1347
1348	call dadal(v,nd2)
1349
1350	return
1351	end
1352	subroutine expnd2(h,x,w,eps,nrmax)
1353	implicit none
1354	integer j,ndim,ndim2,nrmax,ntt
1355	double precision eps
1356	! WRAPPING EXPFLOD USING EXP1D
1357	parameter (ndim=3)
1358	parameter (ndim2=6)
1359	parameter (ntt=40)
1360	integer nd,nd2,no,nv
1361	common /ii/no,nv,nd,nd2
1362	integer x(),w(),h
1363
1364	integer b0,v(ndim2)
1365	!
1366	!
1367	call etallnom(b0,1,'B0 ')
1368	call etallnom(v,nd2 ,'V ')
1369
1370	call dacopd(x,v)
1371	do j=1,nd2
1372	call exp1d(h,v(j),b0,eps,nrmax)
1373	call dacop(b0,v(j))
1374	enddo
1375	call dacopd(v,w)
1376	call dadal(v,nd2)
1377	call dadal(b0,1)
1378	return
1379	end
1380	subroutine flofacg(xy,h,epsone)
1381	implicit none
1382	integer i,k,kk,ndim,ndim2,nrmax,ntt
1383	double precision eps,epsone,r,xn,xnbefore,xnorm,xnorm1,xx
1384	! GENERAL ONE EXPONENT FACTORIZATION
1385	parameter (ndim=3)
1386	parameter (ndim2=6)
1387	parameter (ntt=40)
1388	integer idpr
1389	common /printing/ idpr
1390	logical more
1391	integer nd,nd2,no,nv
1392	common /ii/no,nv,nd,nd2
1393	double precision xintex
1394	common /integratedex/ xintex(0:20)
1395	integer xy(),x(ndim2),h()
1396	integer v(ndim2),w(ndim2),t(ndim2), z(ndim2)
1397	integer jj(ntt)
1398	jj(1)=1
1399	!
1400	call etallnom(v,nd2 ,'V ')
1401	call etallnom(w,nd2 ,'W ')
1402	call etallnom(t,nd2 ,'T ')
1403	call etallnom(x,nd2 ,'Z ')
1404	call etallnom(z,nd2 ,'Z ')
1405
1406	call etini(v)
1407	call daclrd(w)
1408	xnorm1=0.d0
1409	do i=1,nd2
1410	call daabs(xy(i),r)
1411	xnorm1=xnorm1+r
1412	enddo
1413	xnbefore=1.d36
1414	more=.false.
1415	eps=1.e-9
1416	nrmax=1000
1417	xn=10000.d0
1418	do k=1,nrmax
1419	call dacmud(h,-1.d0,t)
1420	call expflod(t,xy,x,eps,nrmax)
1421	call dalind(x,1.d0,v,-1.d0,t)
1422	! write(20,*) "$$$$$$$$$$$$$$",k,"$$$$$$$$$$$$$$$$$$$$"
1423	! call daprid(t,1,1,20)
1424	if(xn.lt.epsone) then
1425	if(idpr.ge.0) write(6,*) "xn quadratic",xn
1426	call daflod(t,t,w)
1427	call dalind(t,1.d0,w,-0.5d0,t)
1428	call dacopd(t,z)
1429	call dacopd(t,w)
1430	! second order in W
1431	call etcom(h,w,x)
1432	call etcom(x,w,x)
1433	! END OF order in W
1434
1435	do kk=1,10
1436	call etcom(h,w,w)
1437	call dalind(z,1.d0,w,xintex(kk),z)
1438	enddo
1439	call dacopd(z,t)
1440	xx=1.d0/12.d0
1441	call dalind(x,xx,h,1.d0,h)
1442	endif
1443
1444	call dalind(t,1.d0,h,1.d0,h)
1445	xnorm=0.d0
1446	do i=1,nd2
1447	call daabs(t(i),r)
1448	xnorm=xnorm+r
1449	enddo
1450	xn=xnorm/xnorm1
1451	if(xn.ge.epsone.and.(idpr.ge.0)) write(6,*)" xn linear ",xn
1452	if(xn.lt.eps.or.more) then
1453	more=.true.
1454	if(xn.ge.xnbefore) goto 1000
1455	xnbefore=xn
1456	endif
1457	enddo
1458	1000 if(idpr.ge.0) write(6,*) " iteration " , k
1459	call dadal(x,nd2)
1460	call dadal(w,nd2)
1461	call dadal(v,nd2)
1462	call dadal(t,nd2)
1463	call dadal(z,nd2)
1464	return
1465	end
1466	subroutine flofac(xy,x,h)
1467	implicit none
1468	integer k,ndim,ndim2,ntt
1469	! GENERAL DRAGT-FINN FACTORIZATION
1470	parameter (ndim=3)
1471	parameter (ndim2=6)
1472	parameter (ntt=40)
1473	integer nd,nd2,no,nv
1474	common /ii/no,nv,nd,nd2
1475	integer xy(),x(),h(*)
1476	integer v(ndim2),w(ndim2)
1477	!
1478	call etallnom(v,nd2 ,'V ')
1479	call etallnom(w,nd2 ,'W ')
1480
1481	call dacopd(xy,x)
1482	call dacopd(x,v)
1483	call daclrd(w)
1484	call danot(1)
1485	call etinv(v,w)
1486	call danot(no)
1487	call etcct(x,w,v)
1488	call danot(1)
1489	call dacopd(xy,x)
1490	call danot(no)
1491	call dacopd(v,w)
1492	call daclrd(h)
1493	do k=2,no
1494	call taked(w,k,v)
1495	call dalind(v,1.d0,h,1.d0,h)
1496	call facflod(h,w,v,k,k,-1.d0,-1)
1497	call dacopd(v,w)
1498	enddo
1499	call dadal(w,nd2)
1500	call dadal(v,nd2)
1501	return
1502	end
1503	subroutine liefact(xy,x,h)
1504	implicit none
1505	integer ndim,ndim2,ntt
1506	! SYMPLECTIC DRAGT-FINN FACTORIZATION WRAPPING FLOFAC
1507	parameter (ndim=3)
1508	parameter (ndim2=6)
1509	parameter (ntt=40)
1510	integer nd,nd2,no,nv
1511	common /ii/no,nv,nd,nd2
1512	integer xy(),x(),h
1513
1514	integer v(ndim2)
1515
1516	call etall(v,nd2)
1517
1518	call flofac(xy,x,v)
1519	call intd(v,h,-1.d0)
1520	!
1521	call dadal(v,nd2)
1522
1523	return
1524	end
1525
1526
1527	logical*1 function mapnorm(x,ft,a2,a1,xy,h,nord)
1528	implicit none
1529	integer isi,ndim,ndim2,nord,ntt
1530	!--NORMALIZATION ROUTINES OF LIELIB
1531	!- WRAPPING MAPNORMF
1532	parameter (ndim=3)
1533	parameter (ndim2=6)
1534	parameter (ntt=40)
1535	integer nd,nd2,no,nv
1536	common /ii/no,nv,nd,nd2
1537	integer x(),a1(),a2(),ft,xy(),h,hf(ndim2),ftf(ndim2)
1538	logical*1 mapnormf
1539
1540	call etall(ftf,nd2)
1541	call etall(hf,nd2)
1542	isi=0
1543	mapnorm = mapnormf(x,ftf,a2,a1,xy,hf,nord,isi)
1544	call intd(hf,h,-1.d0)
1545	call intd(ftf,ft,-1.d0)
1546	call dadal(ftf,nd2)
1547	call dadal(hf,nd2)
1548
1549	return
1550	end
1551	subroutine gettura(psq,radsq)
1552	implicit none
1553	integer ik,ndim,ndim2,ntt
1554	parameter (ndim=3)
1555	parameter (ndim2=6)
1556	parameter (ntt=40)
1557	double precision psq(ndim),radsq(ndim)
1558	integer nd,nd2,no,nv
1559	common /ii/no,nv,nd,nd2
1560	double precision ps,rads
1561	common /tunerad/ ps(ndim),rads(ndim)
1562
1563	do ik=1,nd
1564	psq(ik)=ps(ik)
1565	radsq(ik)=rads(ik)
1566	enddo
1567
1568	return
1569	end
1570	subroutine setidpr(idprint,nplan)
1571	implicit none
1572	integer idprint,ik,ndim,ndim2,nplan
1573	parameter (ndim=3)
1574	parameter (ndim2=6)
1575	dimension nplan(ndim)
1576	integer nd,nd2,no,nv
1577	common /ii/no,nv,nd,nd2
1578	integer idpr
1579	common /printing/ idpr
1580	integer nplane
1581	double precision epsplane,xplane
1582	common /choice/ xplane(ndim),epsplane,nplane(ndim)
1583
1584	do ik=1,nd
1585	nplane(ik)=nplan(ik)
1586	enddo
1587	idpr=idprint
1588
1589	return
1590	end
1591	subroutine idprset(idprint)
1592	implicit none
1593	integer idprint,ndim,ndim2
1594	parameter (ndim=3)
1595	parameter (ndim2=6)
1596	integer nd,nd2,no,nv
1597	common /ii/no,nv,nd,nd2
1598	integer idpr
1599	common /printing/ idpr
1600	integer nplane
1601	double precision epsplane,xplane
1602	common /choice/ xplane(ndim),epsplane,nplane(ndim)
1603
1604	idpr=idprint
1605
1606	return
1607	end
1608
1609
1610	logical*1 function mapnormf(x,ft,a2,a1,xy,h,nord,isi)
1611	implicit none
1612	integer ij,isi,ndim,ndim2,nord,ntt
1613	double precision angle,p,rad,st,x2pi,x2pii
1614	parameter (ndim=3)
1615	parameter (ndim2=6)
1616	parameter (ntt=40)
1617	dimension angle(ndim),st(ndim),p(ndim),rad(ndim)
1618	integer nd,nd2,no,nv
1619	common /ii/no,nv,nd,nd2
1620	integer ndc,ndc2,ndpt,ndt
1621	common /coast/ndc,ndc2,ndt,ndpt
1622	integer x(),a1(),a2(),ft(),xy(),h()
1623	integer itu
1624	common /tunedef/itu
1625	integer idpr
1626	common /printing/ idpr
1627	integer iflow,jtune
1628	common /vecflow/ iflow,jtune
1629	double precision ps,rads
1630	common /tunerad/ ps(ndim),rads(ndim)
1631	integer a1i(ndim2),a2i(ndim2)
1632	logical*1 midbflo
1633	!
1634	call etallnom(a1i,nd2 ,'A1I ')
1635	call etallnom(a2i,nd2 ,'A2I ')
1636	! frank/etienne
1637	do itu=1,ndim
1638	angle(itu)=0.d0
1639	p(itu)=0.d0
1640	st(itu)=0.d0
1641	rad(itu)=0.d0
1642	ps(itu)=0.d0
1643	rads(itu)=0.d0
1644	enddo
1645	jtune=isi
1646	x2pii=1.d0/datan(1.d0)/8.d0
1647	x2pi=datan(1.d0)*8.d0
1648	call dacopd(x,xy)
1649	! go to fix point in the parameters + pt to order nord>=1
1650	call gofix(xy,a1,a1i,nord)
1651	call simil(a1i,xy,a1,xy)
1652	! linear part
1653	mapnormf = midbflo(xy,a2,a2i,angle,rad,st)
1654	do ij=1,nd-ndc
1655	p(ij)=angle(ij)(st(ij)(x2pii-1.d0)+1.d0)
1656	enddo
1657	if(ndc.eq.1) p(nd)=angle(nd)
1658	if(idpr.ge.0) then
1659	write(6,*) 'tune ',(p(ij),ij=1,nd)
1660	write(6,*) 'damping ', (rad(ij),ij=1,nd)
1661	endif
1662	do ij=1,nd ! -ndc Frank
1663	ps(ij)=p(ij)
1664	rads(ij)=rad(ij)
1665	enddo
1666	call initpert(st,angle,rad)
1667	call simil(a2i,xy,a2,xy)
1668	call dacopd(xy,a2i)
1669	! write(6,*) 'Entering orderflo'
1670	call orderflo(h,ft,xy,angle,rad)
1671	do ij=1,nd-ndc
1672	p(ij)=angle(ij)
1673	if(angle(ij).gt.x2pi/2.d0.and.st(ij).gt.0.d0.and.itu.eq.1)then
1674	p(ij)=angle(ij)-x2pi
1675	write(6,*) ij,' TH TUNE MODIFIED IN H2 TO ',p(ij)/x2pi
1676	endif
1677	enddo
1678	call h2pluflo(h,p,rad)
1679	! CALL TAKED(A2I,1,XY)
1680	call taked(a2i,1,a1i)
1681	call etcct(xy,a1i,xy)
1682
1683	call dadal(a2i,nd2)
1684	call dadal(a1i,nd2)
1685	return
1686	end
1687	subroutine gofix(xy,a1,a1i,nord)
1688	implicit none
1689	integer i,ndim,ndim2,nord,ntt
1690	double precision xic
1691	! GETTING TO THE FIXED POINT AND CHANGING TIME APPROPRIATELY IN THE
1692	! COASTING BEAM CASE
1693
1694	!****************************************************************
1695	! X = A1 XY A1I WHERE X IS TO THE FIXED POINT TO ORDER NORD
1696	! for ndpt not zero, works in all cases. (coasting beam: eigenvalue
1697	!1 in Jordan form)
1698	!****************************************************************
1699	parameter (ndim=3)
1700	parameter (ndim2=6)
1701	parameter (ntt=40)
1702	integer nd,nd2,no,nv
1703	common /ii/no,nv,nd,nd2
1704	integer ndc,ndc2,ndpt,ndt
1705	common /coast/ndc,ndc2,ndt,ndpt
1706	integer xy(),a1(),a1i(*)
1707
1708	integer x(ndim2),w(ndim2),v(ndim2),rel(ndim2)
1709	!
1710	call etallnom(x,nd2 , 'X ')
1711	call etallnom(w,nd2 , 'W ')
1712	call etallnom(v,nd2 , 'V ')
1713	call etallnom(rel,nd2 ,'REL ')
1714
1715	! COMPUTATION OF A1 AND A1I USING DAINV
1716	call etini(rel)
1717
1718	call danot(nord)
1719
1720	call etini(v)
1721
1722	do i=1,nd2-ndc2
1723	call dacop(xy(i),x(i))
1724	call dalin(x(i),1.d0,rel(i),-1.d0,v(i))
1725	enddo
1726	call etinv(v,w)
1727	call daclrd(x)
1728	if(ndc.eq.1) then
1729	call davar(x(ndpt),0.d0,ndpt)
1730	endif
1731	call etcct(w,x,v)
1732	if(ndc.eq.1) then
1733	call daclr(v(nd2))
1734	call daclr(v(nd2-ndc))
1735	endif
1736	call dalind(rel,1.d0,v,1.d0,a1)
1737	call dalind(rel,1.d0,v,-1.d0,a1i)
1738
1739	if(ndpt.ne.0) then
1740
1741	! CORRECTIONS
1742	call daclrd(w)
1743	call daclrd(v)
1744	call daclrd(x)
1745
1746	do i=1,nd2-ndc2
1747	call dalin(a1(i),1.d0,rel(i),-1.d0,w(i))
1748	enddo
1749
1750	! COMPUTE Deta/Ddelta
1751	call dacopd(w,a1)
1752
1753	do i=1,nd2-ndc2
1754	call dader(ndpt,w(i),w(i))
1755	enddo
1756	! COMPUTE J*Deta/dDELTA
1757
1758	do i=1,nd-ndc
1759	call dacmu(w(2i),1.d0,v(2i-1) )
1760	call dacmu(w(2i-1),-1.d0,v(2i) )
1761	enddo
1762
1763	xic=(-1)**(ndt)
1764
1765	do i=1,nd2-ndc2
1766	call damul(v(i),rel(i),x(1))
1767	call daadd(x(1),w(ndt),w(ndt))
1768	call dacop(a1(i),w(i))
1769	enddo
1770	call dacmu(w(ndt),xic,w(ndt))
1771
1772	call expflod(w,rel,a1,1.d-7,10000)
1773	! END OF CORRECTIONS
1774
1775	call etinv(a1,a1i)
1776
1777	endif
1778
1779	call danot(no)
1780
1781	call dadal(rel,nd2)
1782	call dadal(v,nd2)
1783	call dadal(w,nd2)
1784	call dadal(x,nd2)
1785	return
1786	end
1787	double precision function transver(j)
1788	implicit none
1789	integer i,ic,ndim
1790	! USED IN A DACFU CALL OF GOFIX
1791	parameter (ndim=3)
1792	! PARAMETER (NTT=40)
1793	! INTEGER J(NTT)
1794	integer j(*)
1795	integer nd,nd2,no,nv
1796	common /ii/no,nv,nd,nd2
1797	integer ndc,ndc2,ndpt,ndt
1798	common /coast/ndc,ndc2,ndt,ndpt
1799
1800	transver=1.d0
1801	ic=0
1802	do i=1,nd2-ndc2
1803	ic=ic+j(i)
1804	enddo
1805	if(ic.ne.1) transver=0.d0
1806	return
1807	end
1808	subroutine orderflo(h,ft,x,ang,ra)
1809	implicit none
1810	integer k,ndim,ndim2,ntt
1811	double precision ang,ra
1812	!- NONLINEAR NORMALIZATION PIECE OF MAPNORMF
1813	parameter (ndim=3)
1814	parameter (ndim2=6)
1815	parameter (ntt=40)
1816	integer idpr
1817	common /printing/ idpr
1818	integer nd,nd2,no,nv
1819	common /ii/no,nv,nd,nd2
1820	dimension ang(ndim),ra(ndim)
1821	integer x(),h(),ft(*)
1822	integer w(ndim2),v(ndim2),rel(ndim2)
1823	integer roi(ndim2)
1824	integer b1(ndim2),b5(ndim2),b6(ndim2),b9(ndim2)
1825	!
1826	call etallnom(w,nd2 ,'W ')
1827	call etallnom(v,nd2 ,'V ')
1828	call etallnom(rel,nd2 ,'REL ')
1829	call etallnom(roi,nd2 ,'ROI ')
1830	call etallnom(b1,nd2 ,'B1 ')
1831	call etallnom(b5,nd2 ,'B5 ')
1832	call etallnom(b6,nd2 ,'B6 ')
1833	call etallnom(b9,nd2 ,'B9 ')
1834	call rotiflo(roi,ang,ra)
1835	call etini(rel)
1836	call daclrd(h)
1837	call daclrd(ft)
1838	call etcct(x,roi,x)
1839	do k=2,no
1840	! IF K>2 V = H(K)^-1 X(K)
1841	call facflod(h,x,v,2,k-1,-1.d0,-1)
1842	! EXTRACTING K TH DEGREE OF V ----> W
1843	call taked(v,k,w)
1844	! write(16,*) "$$$$$$$$ K $$$$$$$$$$", k
1845	! W = EXP(B5) + ...
1846	call dacopd(w,b5)
1847	! CALL INTD(W,B5,-1.D0)
1848	! B5 ON EXIT IS THE NEW CONTRIBUTION TO H
1849	! B6 IS THE NEW CONTRIBUTION TO FT
1850	call nuanaflo(b5,b6)
1851	call dalind(b5,1.d0,h,1.d0,b1)
1852	call dacopd(b1,h)
1853	! EXP(B9) = EXP( : ROTI B6 :)
1854	call trxflo(b6,b9,roi)
1855
1856	! V = EXP(-B6) REL
1857	call facflod(b6,rel,v,k,k,-1.d0,1)
1858	! W = V o X
1859	call etcct(v,x,w)
1860	if(idpr.ge.0) then
1861	write(6,*) ' ORDERFLO K= ', k
1862	endif
1863	! X = EXP(B9) W
1864	call facflod(b9,w,x,k,k,1.d0,1)
1865	! B6 IS THE NEW CONTRIBUTION TO FT
1866	call dalind(b6,1.d0,ft,1.d0,b1)
1867	call dacopd(b1,ft)
1868	enddo
1869	call dadal(b9,nd2)
1870	call dadal(b6,nd2)
1871	call dadal(b5,nd2)
1872	call dadal(b1,nd2)
1873	call dadal(roi,nd2)
1874	call dadal(rel,nd2)
1875	call dadal(v,nd2)
1876	call dadal(w,nd2)
1877	return
1878	end
1879	subroutine nuanaflo(h,ft)
1880	implicit none
1881	integer i,ndim,ndim2,ntt
1882	double precision dfilt,filt,xgam,xgbm
1883	! RESONANCE DENOMINATOR OPERATOR (1-R^-1)^-1
1884	parameter (ndim=3)
1885	parameter (ndim2=6)
1886	parameter (ntt=40)
1887	integer nd,nd2,no,nv
1888	common /ii/no,nv,nd,nd2
1889	integer iflow,jtune
1890	common /vecflow/ iflow,jtune
1891	external xgam,xgbm,dfilt,filt
1892	integer h(),ft(),br(ndim2),bi(ndim2),c(ndim2),ci(ndim2)
1893	integer t1(2),t2(2)
1894
1895	call etall(br,nd2)
1896	call etall(bi,nd2)
1897	call etall(c,nd2)
1898	call etall(ci,nd2)
1899
1900	call ctorflo(h,br,bi)
1901
1902	! FILTERING RESONANCES AND TUNE SHIFTS
1903	! ASSUMING REALITY I.E. B(2I-1)=CMPCJG(B(2I))
1904
1905	do i=1,nd2
1906	iflow=i
1907	call dacfu(br(i),filt,c(i))
1908	call dacfu(bi(i),filt,ci(i))
1909	enddo
1910	call rtocflo(c,ci,h)
1911
1912	do i=1,nd2
1913
1914	iflow=i
1915	call dacfu(br(i),dfilt,br(i))
1916	call dacfu(bi(i),dfilt,bi(i))
1917	enddo
1918	! NOW WE MUST REORDER C AND CI TO SEPARATE THE REAL AND IMAGINARY PART
1919	! THIS IS NOT NECESSARY WITH :H: OPERATORS
1920
1921	do i=1,nd2
1922	t1(1)=br(i)
1923	t1(2)=bi(i)
1924	t2(1)=c(i)
1925	t2(2)=ci(i)
1926	iflow=i
1927	call comcfu(t1,xgam,xgbm,t2)
1928	enddo
1929
1930	call rtocflo(c,ci,ft)
1931
1932	call dadal(br,nd2)
1933	call dadal(bi,nd2)
1934	call dadal(c,nd2)
1935	call dadal(ci,nd2)
1936
1937	return
1938	end
1939	double precision function xgam(j)
1940	implicit none
1941	integer i,ic,ij,ik,ndim,ndim2
1942	double precision ad,ans,as,ex,exh
1943	! XGAM AND XGBM ARE THE EIGENVALUES OF THE OPERATOR NEWANAFLO
1944	parameter (ndim=3)
1945	parameter (ndim2=6)
1946	! PARAMETER (NTT=40)
1947	integer iflow,jtune
1948	common /vecflow/ iflow,jtune
1949	double precision angle,dsta,rad,sta
1950	common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
1951	integer idsta,ista
1952	common /istable/ista(ndim),idsta(ndim)
1953	integer nd,nd2,no,nv
1954	common /ii/no,nv,nd,nd2
1955	integer ndc,ndc2,ndpt,ndt
1956	common /coast/ndc,ndc2,ndt,ndpt
1957	! INTEGER J(NTT),JJ(NDIM),JP(NDIM)
1958	integer j(*),jj(ndim),jp(ndim)
1959	xgam=0.d0
1960	ad=0.d0
1961	as=0.d0
1962	ic=0
1963	do i=1,nd-ndc
1964	ik=2*i-1
1965	ij=2*i
1966	jp(i)=j(ik)+j(ij)
1967	jj(i)=j(ik)-j(ij)
1968	if(ik.eq.iflow.or.ij.eq.iflow) then
1969	jj(i)=jj(i)+(-1)**iflow
1970	jp(i)=jp(i)-1
1971	endif
1972	ic=ic+iabs(jj(i))
1973	enddo
1974
1975	do i=1,nd-ndc
1976	ad=dsta(i)dble(jj(i))angle(i)-dble(jp(i))*rad(i)+ad
1977	as=sta(i)dble(jj(i))angle(i)+as
1978	enddo
1979
1980	exh=dexp(ad/2.d0)
1981	ex=exh**2
1982	ans=4.d0ex(dsinh(ad/2.d0)2+dsin(as/2.d0)2)
1983	xgam=2.d0(-exhdsinh(ad/2.d0)+exdsin(as/2.d0)*2)/ans
1984
1985	return
1986	end
1987	double precision function xgbm(j)
1988	implicit none
1989	integer i,ic,ij,ik,ndim,ndim2
1990	double precision ad,ans,as,ex,exh
1991	parameter (ndim=3)
1992	parameter (ndim2=6)
1993	! PARAMETER (NTT=40)
1994	integer iflow,jtune
1995	common /vecflow/ iflow,jtune
1996	double precision angle,dsta,rad,sta
1997	common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
1998	integer idsta,ista
1999	common /istable/ista(ndim),idsta(ndim)
2000	integer nd,nd2,no,nv
2001	common /ii/no,nv,nd,nd2
2002	integer ndc,ndc2,ndpt,ndt
2003	common /coast/ndc,ndc2,ndt,ndpt
2004	! INTEGER J(NTT),JJ(NDIM),JP(NDIM)
2005	integer j(*),jj(ndim),jp(ndim)
2006	xgbm=0.d0
2007	ad=0.d0
2008	as=0.d0
2009	ic=0
2010	do i=1,nd-ndc
2011	ik=2*i-1
2012	ij=2*i
2013	jp(i)=j(ik)+j(ij)
2014	jj(i)=j(ik)-j(ij)
2015	if(ik.eq.iflow.or.ij.eq.iflow) then
2016	jj(i)=jj(i)+(-1)**iflow
2017	jp(i)=jp(i)-1
2018	endif
2019	ic=ic+iabs(jj(i))
2020	enddo
2021
2022	do i=1,nd-ndc
2023	ad=dsta(i)dble(jj(i))angle(i)-dble(jp(i))*rad(i)+ad
2024	as=sta(i)dble(jj(i))angle(i)+as
2025	enddo
2026
2027	exh=dexp(ad/2.d0)
2028	ex=exh**2
2029	ans=4.d0ex(dsinh(ad/2.d0)2+dsin(as/2.d0)2)
2030	xgbm=dsin(as)*ex/ans
2031
2032	return
2033	end
2034	double precision function filt(j)
2035	implicit none
2036	integer i,ic,ic1,ic2,ij,ik,ji,ndim,ndim2,nreso
2037	! PROJECTION FUNCTIONS ON THE KERNEL ANMD RANGE OF (1-R^-1)
2038	!- THE KERNEL OF (1-R^-1)
2039	parameter (ndim=3)
2040	parameter (ndim2=6)
2041	! PARAMETER (NTT=40)
2042	parameter (nreso=20)
2043	double precision angle,dsta,rad,sta
2044	common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
2045	integer idsta,ista
2046	common /istable/ista(ndim),idsta(ndim)
2047	integer nd,nd2,no,nv
2048	common /ii/no,nv,nd,nd2
2049	integer ndc,ndc2,ndpt,ndt
2050	common /coast/ndc,ndc2,ndt,ndpt
2051	integer mx,nres
2052	common /reson/mx(ndim,nreso),nres
2053	integer iflow,jtune
2054	common /vecflow/ iflow,jtune
2055	! INTEGER J(NTT),JJ(NDIM)
2056	integer j(*),jj(ndim)
2057
2058	filt=1.d0
2059
2060	ic=0
2061	do i=1,nd-ndc
2062	ik=2*i-1
2063	ij=2*i
2064	jj(i)=j(ik)-j(ij)
2065	if(ik.eq.iflow.or.ij.eq.iflow) then
2066	jj(i)=jj(i)+(-1)**iflow
2067	endif
2068	ic=ic+iabs(jj(i))
2069	enddo
2070
2071	if(ic.eq.0.and.jtune.eq.0) return
2072
2073	do i=1,nres
2074	ic1=1
2075	ic2=1
2076	do ji=1,nd-ndc
2077	if(mx(ji,i).ne.jj(ji)) ic1=0
2078	if(mx(ji,i).ne.-jj(ji)) ic2=0
2079	if(ic1.eq.0.and.ic2.eq.0) goto 3
2080	enddo
2081	return
2082	3 continue
2083	enddo
2084
2085	filt=0.d0
2086	return
2087	end
2088	double precision function dfilt(j)
2089	implicit none
2090	integer ndim,ndim2,nreso
2091	double precision fil,filt
2092	!- THE RANGE OF (1-R^-1)^1
2093	!- CALLS FILT AND EXCHANGES 1 INTO 0 AND 0 INTO 1.
2094	parameter (ndim=3)
2095	parameter (ndim2=6)
2096	! PARAMETER (NTT=40)
2097	parameter (nreso=20)
2098	integer nd,nd2,no,nv
2099	common /ii/no,nv,nd,nd2
2100	integer ndc,ndc2,ndpt,ndt
2101	common /coast/ndc,ndc2,ndt,ndpt
2102	integer mx,nres
2103	common /reson/mx(ndim,nreso),nres
2104	external filt
2105	! INTEGER J(NTT)
2106	integer j(*)
2107
2108	fil=filt(j)
2109	if(fil.gt.0.5d0) then
2110	dfilt=0.d0
2111	else
2112	dfilt=1.d0
2113	endif
2114	return
2115	end
2116	subroutine dhdjflo(h,t)
2117	implicit none
2118	integer i,ndim,ndim2,ntt
2119	double precision coe,x2pi
2120	! CONVENIENT TUNE SHIFT FINDED FOR SYMPLECTIC CASE (NU,DL)(H)=T
2121	parameter (ndim=3)
2122	parameter (ndim2=6)
2123	parameter (ntt=40)
2124	integer nd,nd2,no,nv
2125	common /ii/no,nv,nd,nd2
2126	integer ndc,ndc2,ndpt,ndt
2127	common /coast/ndc,ndc2,ndt,ndpt
2128	integer h(),t()
2129
2130	integer b1(ndim2),b2(ndim2),bb1,bb2
2131	!
2132	call etall(b1,nd2)
2133	call etall(b2,nd2)
2134	call etall(bb1,1)
2135	call etall(bb2,1)
2136
2137	x2pi=datan(1.d0)*8.d0
2138	call ctorflo(h,b1,b2)
2139	coe=1.d0/x2pi
2140
2141	do i=1,nd-ndc
2142	call datra(2i,b2(2i),bb1)
2143	call dacmu(bb1,coe,t(i+nd))
2144	call dacop(t(i+nd),bb1)
2145	call daclr(bb2)
2146	call rtoc(bb1,bb2,bb1)
2147	call dacop(bb1,t(i))
2148	enddo
2149
2150	if(ndpt.ne.0) then
2151	call dacop(h(ndt),t(nd))
2152	call dacop(b1(ndt),t(nd2))
2153	endif
2154
2155	call dadal(bb2,1)
2156	call dadal(bb1,1)
2157	call dadal(b2,nd2)
2158	call dadal(b1,nd2)
2159	return
2160	end
2161	subroutine dhdj(h,t)
2162	implicit none
2163	integer i,ndim,ndim2,ntt
2164	double precision coe,x2pi
2165	parameter (ndim=3)
2166	parameter (ndim2=6)
2167	parameter (ntt=40)
2168	integer nd,nd2,no,nv
2169	common /ii/no,nv,nd,nd2
2170	integer ndc,ndc2,ndpt,ndt
2171	common /coast/ndc,ndc2,ndt,ndpt
2172	integer h,t(*)
2173
2174	integer b1,b2,bb1,bb2
2175	!
2176	call etallnom(b1,1,'B1 ')
2177	call etallnom(b2,1,'B2 ')
2178	call etallnom(bb1,1,'BB1 ')
2179	call etallnom(bb2,1,'BB2 ')
2180
2181	x2pi=datan(1.d0)*8.d0
2182	call ctor(h,b1,b2)
2183	coe=-2.d0/x2pi
2184	do i=1,nd-ndc
2185	call dader(2*i-1,b1,b2)
2186	call datra(2*i,b2,bb2)
2187	call dacmu(bb2,coe,t(i+nd))
2188	call dacop(t(i+nd),bb2)
2189	call daclr(b2)
2190	call rtoc(bb2,b2,bb1)
2191	call dacop(bb1,t(i))
2192	enddo
2193
2194	if(ndpt.eq.nd2) then
2195	call dader(ndpt,h,t(nd))
2196	call dader(ndpt,b1,t(nd2))
2197	call dacmu(t(nd),-1.d0,t(nd))
2198	call dacmu(t(nd2),-1.d0,t(nd2))
2199	endif
2200	if(ndt.eq.nd2) then
2201	call dader(ndpt,h,t(nd))
2202	call dader(ndpt,b1,t(nd2))
2203	endif
2204	call dadal(bb2,1)
2205	call dadal(bb1,1)
2206	call dadal(b2,1)
2207	call dadal(b1,1)
2208	return
2209	end
2210	subroutine h2pluflo(h,ang,ra)
2211	implicit none
2212	integer i,j,ndim,ndim2,ntt
2213	double precision ang,r1,r2,ra,st
2214	! POKES IN \VEC{H} ANGLES AND DAMPING COEFFFICIENTS
2215	parameter (ndim=3)
2216	parameter (ndim2=6)
2217	parameter (ntt=40)
2218	double precision angle,dsta,rad,sta
2219	common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
2220	dimension ang(ndim),st(ndim),ra(ndim),j(ntt)
2221	integer nd,nd2,no,nv
2222	common /ii/no,nv,nd,nd2
2223	integer ndc,ndc2,ndpt,ndt
2224	common /coast/ndc,ndc2,ndt,ndpt
2225	integer h(*)
2226	!*DAEXT(NO,NV) H
2227
2228	do i=1,nd
2229	st(i)=2.d0*sta(i)-1.d0
2230	enddo
2231
2232	do i=1,ntt
2233	j(i)=0
2234	enddo
2235
2236	do i=1,nd-ndc
2237	j(2*i-1)=1
2238	r1=-ang(i)
2239	!-----
2240	call dapok(h(2*i),j,r1)
2241
2242	r2=ra(i)
2243	call dapok(h(2*i-1),j,r2)
2244	j(2*i-1)=0
2245
2246	j(2*i)=1
2247	r1=ang(i)*st(i)
2248	call dapok(h(2*i-1),j,r1)
2249	call dapok(h(2*i),j,r2)
2250	j(2*i)=0
2251
2252	enddo
2253
2254	if(ndpt.eq.nd2-1) then
2255	j(ndpt)=1
2256	call dapok(h(ndt),j,ang(nd))
2257	elseif(ndpt.eq.nd2) then
2258	j(ndpt)=1
2259	call dapok(h(ndt),j,-ang(nd))
2260	endif
2261	return
2262	end
2263	subroutine rotflo(ro,ang,ra)
2264	implicit none
2265	integer i,ndim,ndim2,ntt
2266	double precision ang,ch,co,ra,sh,si,sim,xx
2267	! CREATES R AND R^-1 USING THE EXISTING ANGLES AND DAMPING
2268	! COULD BE REPLACED BY A CALL H2PLUFLO FOLLOWED BY EXPFLOD
2269	! CREATES R
2270	parameter (ndim=3)
2271	parameter (ndim2=6)
2272	parameter (ntt=40)
2273	integer idsta,ista
2274	common /istable/ista(ndim),idsta(ndim)
2275	dimension co(ndim),si(ndim),ang(ndim),ra(ndim)
2276	integer j(ntt)
2277	integer nd,nd2,no,nv
2278	common /ii/no,nv,nd,nd2
2279	integer ndc,ndc2,ndpt,ndt
2280	common /coast/ndc,ndc2,ndt,ndpt
2281	integer ro(*)
2282	call daclrd(ro)
2283	do i=1,nd-ndc
2284	xx=dexp(ra(i))
2285	if(ista(i).eq.0) then
2286	call hyper(ang(i),ch,sh)
2287	co(i)=ch*xx
2288	si(i)=-sh*xx
2289	else
2290	co(i)=dcos(ang(i))*xx
2291	si(i)=dsin(ang(i))*xx
2292	endif
2293	enddo
2294	do i=1,nd-ndc
2295	if(ista(i).eq.0)then
2296	sim=si(i)
2297	else
2298	sim=-si(i)
2299	endif
2300	j(2*i-1)=1
2301	call dapok(ro(2*i-1),j,co(i))
2302	call dapok(ro(2*i),j,sim)
2303	j(2*i-1)=0
2304	j(2*i)=1
2305	call dapok(ro(2*i),j,co(i))
2306	call dapok(ro(2*i-1),j,si(i))
2307	j(2*i)=0
2308	enddo
2309
2310	if(ndc.eq.1) then
2311	j(ndt)=1
2312	call dapok(ro(ndt),j,1.d0)
2313	call dapok(ro(ndpt),j,0.d0)
2314	j(ndt)=0
2315	j(ndpt)=1
2316	call dapok(ro(ndt),j,ang(nd))
2317	call dapok(ro(ndpt),j,1.d0)
2318	j(ndpt)=0
2319	endif
2320
2321	return
2322	end
2323	subroutine rotiflo(roi,ang,ra)
2324	implicit none
2325	integer i,ndim,ndim2,ntt
2326	double precision ang,ch,co,ra,sh,si,sim,simv,xx
2327	! CREATES R^-1
2328	parameter (ndim=3)
2329	parameter (ndim2=6)
2330	parameter (ntt=40)
2331	integer idsta,ista
2332	common /istable/ista(ndim),idsta(ndim)
2333	dimension co(ndim),si(ndim),ang(ndim),ra(ndim)
2334	integer j(ntt)
2335	integer nd,nd2,no,nv
2336	common /ii/no,nv,nd,nd2
2337	integer ndc,ndc2,ndpt,ndt
2338	common /coast/ndc,ndc2,ndt,ndpt
2339	integer roi(*)
2340
2341	do i=1,10
2342	j(i)=0
2343	enddo
2344
2345	call daclrd(roi)
2346	do i=1,nd-ndc
2347	xx=dexp(-ra(i))
2348	if(ista(i).eq.0) then
2349	call hyper(ang(i),ch,sh)
2350	co(i)=ch*xx
2351	si(i)=-sh*xx
2352	else
2353	co(i)=dcos(ang(i))*xx
2354	si(i)=dsin(ang(i))*xx
2355	endif
2356	enddo
2357	do i=1,nd-ndc
2358	if(ista(i).eq.0)then
2359	sim=si(i)
2360	else
2361	sim=-si(i)
2362	endif
2363	j(2*i-1)=1
2364	call dapok(roi(2*i-1),j,co(i))
2365	simv=-sim
2366	call dapok(roi(2*i),j,simv)
2367	j(2*i-1)=0
2368	j(2*i)=1
2369	simv=-si(i)
2370	call dapok(roi(2*i),j,co(i))
2371	call dapok(roi(2*i-1),j,simv)
2372	j(2*i)=0
2373	enddo
2374
2375	if(ndc.eq.1) then
2376	j(ndt)=1
2377	call dapok(roi(ndt),j,1.d0)
2378	call dapok(roi(ndpt),j,0.d0)
2379	j(ndt)=0
2380	j(ndpt)=1
2381	call dapok(roi(ndt),j,-ang(nd))
2382	call dapok(roi(ndpt),j,1.d0)
2383	j(ndpt)=0
2384	endif
2385
2386	return
2387	end
2388	subroutine hyper(a,ch,sh)
2389	implicit none
2390	double precision a,ch,sh,x,xi
2391	! USED IN ROTIFLO AND ROTFLO
2392	x=dexp(a)
2393	xi=1.d0/x
2394	ch=(x+xi)/2.d0
2395	sh=(x-xi)/2.d0
2396	return
2397	end
2398	subroutine ctor(c1,r2,i2)
2399	implicit none
2400	integer ndim2,ntt
2401	! CHANGES OF BASIS
2402	! C1------> R2+I R1
2403	parameter (ndim2=6)
2404	parameter (ntt=40)
2405	integer nd,nd2,no,nv
2406	common /ii/no,nv,nd,nd2
2407	integer c1,r2,i2
2408	integer b1,b2,x(ndim2)
2409	!
2410	!
2411	call etallnom(b1,1,'B1 ')
2412	call etallnom(b2,1,'B2 ')
2413	call etallnom(x,nd2 ,'X ')
2414
2415	call ctoi(c1,b1)
2416	call etcjg(x)
2417	call trx(b1,b2,x)
2418	call dalin(b1,.5d0,b2,.5d0,r2)
2419	call dalin(b1,.5d0,b2,-.5d0,i2)
2420	call dadal(x,nd2)
2421	call dadal(b2,1)
2422	call dadal(b1,1)
2423	return
2424	end
2425	subroutine rtoc(r1,i1,c2)
2426	implicit none
2427	integer ndim2,ntt
2428	! INVERSE OF CTOR
2429	parameter (ndim2=6)
2430	parameter (ntt=40)
2431	integer nd,nd2,no,nv
2432	common /ii/no,nv,nd,nd2
2433	integer c2,r1,i1
2434
2435	integer b1
2436	!
2437	call etallnom(b1,1,'B1 ')
2438
2439	call daadd(r1,i1,b1)
2440	call itoc(b1,c2)
2441	call dadal(b1,1)
2442	return
2443	end
2444	subroutine ctorflo(c,dr,di)
2445	implicit none
2446	integer ndim,ndim2,ntt
2447	! FLOW CTOR
2448	parameter (ndim=3)
2449	parameter (ndim2=6)
2450	parameter (ntt=40)
2451	integer dr(),di(),c(*)
2452
2453	call ctord(c,dr,di)
2454	call resvec(dr,di,dr,di)
2455
2456	return
2457	end
2458	subroutine rtocflo(dr,di,c)
2459	implicit none
2460	integer ndim,ndim2,ntt
2461	! FLOW RTOC
2462	parameter (ndim=3)
2463	parameter (ndim2=6)
2464	parameter (ntt=40)
2465	integer nd,nd2,no,nv
2466	common /ii/no,nv,nd,nd2
2467	integer dr(),di(),c(*),er(ndim2),ei(ndim2)
2468
2469	call etall(er,nd2)
2470	call etall(ei,nd2)
2471
2472	call reelflo(dr,di,er,ei)
2473	call rtocd(er,ei,c)
2474
2475	call dadal(er,nd2)
2476	call dadal(ei,nd2)
2477
2478	return
2479	end
2480	subroutine ctord(c,cr,ci)
2481	implicit none
2482	integer i,ndim,ndim2,ntt
2483	! ROUTINES USED IN THE INTERMEDIATE STEPS OF CTORFLO AND RTOCFLO
2484	! SAME AS CTOR OVER ARRAYS CONTAINING ND2 COMPONENTS
2485	! ROUTINE USEFUL IN INTERMEDIATE FLOW CHANGE OF BASIS
2486	parameter (ndim=3)
2487	parameter (ndim2=6)
2488	parameter (ntt=40)
2489	integer nd,nd2,no,nv
2490	common /ii/no,nv,nd,nd2
2491	integer c(),ci(),cr(*)
2492	do i=1,nd2
2493	call ctor(c(i),cr(i),ci(i))
2494	enddo
2495	return
2496	end
2497	subroutine rtocd(cr,ci,c)
2498	implicit none
2499	integer i,ndim,ndim2,ntt
2500	! INVERSE OF CTORD
2501	parameter (ndim=3)
2502	parameter (ndim2=6)
2503	parameter (ntt=40)
2504	integer nd,nd2,no,nv
2505	common /ii/no,nv,nd,nd2
2506	integer c(),ci(),cr(*)
2507	do i=1,nd2
2508	call rtoc(cr(i),ci(i),c(i))
2509	enddo
2510	return
2511	end
2512	subroutine resvec(cr,ci,dr,di)
2513	implicit none
2514	integer i,ndim,ndim2,ntt
2515	! DOES THE SPINOR PART IN CTORFLO
2516	parameter (ndim=3)
2517	parameter (ndim2=6)
2518	parameter (ntt=40)
2519	double precision angle,dsta,rad,sta
2520	common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
2521	integer idsta,ista
2522	common /istable/ista(ndim),idsta(ndim)
2523	integer ndc,ndc2,ndpt,ndt
2524	common /coast/ndc,ndc2,ndt,ndpt
2525	integer nd,nd2,no,nv
2526	common /ii/no,nv,nd,nd2
2527	integer dr(),di(),ci(),cr(),tr(2),ti(2)
2528
2529	call etall(tr,2)
2530	call etall(ti,2)
2531
2532	do i=1,nd-ndc
2533	if(ista(i).eq.1) then
2534	call dasub(cr(2i-1),ci(2i),tr(1))
2535	call daadd(ci(2i-1),cr(2i),ti(1))
2536	call daadd(cr(2i-1),ci(2i),tr(2))
2537	call dasub(ci(2i-1),cr(2i),ti(2))
2538	call dacop(tr(1),dr(2*i-1))
2539	call dacop(tr(2),dr(2*i))
2540	call dacop(ti(1),di(2*i-1))
2541	call dacop(ti(2),di(2*i))
2542	else
2543	call daadd(cr(2i-1),cr(2i),tr(1))
2544	call daadd(ci(2i-1),ci(2i),ti(1))
2545	call dasub(cr(2i-1),cr(2i),tr(2))
2546	call dasub(ci(2i-1),ci(2i),ti(2))
2547	call dacop(tr(1),dr(2*i-1))
2548	call dacop(tr(2),dr(2*i))
2549	call dacop(ti(1),di(2*i-1))
2550	call dacop(ti(2),di(2*i))
2551	endif
2552	enddo
2553
2554	do i=nd2-ndc2+1,nd2
2555	call dacop(cr(i),dr(i))
2556	call dacop(ci(i),di(i))
2557	enddo
2558
2559	call dadal(tr,2)
2560	call dadal(ti,2)
2561	return
2562	end
2563	subroutine reelflo(c,ci,f,fi)
2564	implicit none
2565	integer i,ndim,ndim2,ntt
2566	! DOES THE SPINOR PART IN RTOCFLO
2567	parameter (ndim=3)
2568	parameter (ndim2=6)
2569	parameter (ntt=40)
2570	double precision angle,dsta,rad,sta
2571	common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
2572	integer idsta,ista
2573	common /istable/ista(ndim),idsta(ndim)
2574	integer ndc,ndc2,ndpt,ndt
2575	common /coast/ndc,ndc2,ndt,ndpt
2576	integer nd,nd2,no,nv
2577	common /ii/no,nv,nd,nd2
2578	integer c(),ci(),f(),fi(),e(ndim2),ei(ndim2)
2579
2580	call etall(e,nd2)
2581	call etall(ei,nd2)
2582
2583	do i=1,nd-ndc
2584	call dalin(c(2i-1),0.5d0,c(2i),0.5d0,e(2*i-1))
2585	call dalin(ci(2i-1),0.5d0,ci(2i),0.5d0,ei(2*i-1))
2586	if(ista(i).eq.1) then
2587	call dalin(ci(2i-1),0.5d0,ci(2i),-0.5d0,e(2*i))
2588	call dalin(c(2i-1),-0.5d0,c(2i),0.5d0,ei(2*i))
2589	else
2590	call dalin(ci(2i-1),0.5d0,ci(2i),-0.5d0,ei(2*i))
2591	call dalin(c(2i-1),0.5d0,c(2i),-0.5d0,e(2*i))
2592	endif
2593	enddo
2594
2595	do i=nd2-ndc2+1,nd2
2596	call dacop(c(i),e(i))
2597	call dacop(ci(i),ei(i))
2598	enddo
2599
2600	call dacopd(e,f)
2601	call dacopd(ei,fi)
2602
2603	call dadal(e,nd2)
2604	call dadal(ei,nd2)
2605	return
2606	end
2607	subroutine compcjg(cr,ci,dr,di)
2608	implicit none
2609	integer ndim,ndim2,ntt
2610	! TAKES THE COMPLEX CONJUGATE IN RESONANCE BASIS OF A POLYNOMIAL
2611	parameter (ndim=3)
2612	parameter (ndim2=6)
2613	parameter (ntt=40)
2614	integer nd,nd2,no,nv
2615	common /ii/no,nv,nd,nd2
2616	integer dr,di,ci,cr,x(ndim2)
2617
2618	call etall(x,nd2)
2619
2620	call etcjg(x)
2621	call trx(cr,dr,x)
2622	call trx(ci,di,x)
2623	call dacmu(di,-1.d0,di)
2624
2625	call dadal(x,nd2)
2626	return
2627	end
2628
2629
2630	logical*1 function midbflo(c,a2,a2i,q,a,st)
2631	implicit none
2632	integer i,j,ndim,ndim2,ntt
2633	double precision a,ch,cm,cr,q,r,sa,sai,shm, &
2634	&st,x2pi
2635	! LINEAR EXACT NORMALIZATION USING EIGENVALUE PACKAGE OF NERI
2636	parameter (ntt=40)
2637	parameter (ndim2=6)
2638	parameter (ndim=3)
2639	integer jx(ntt)
2640	integer nd,nd2,no,nv
2641	common /ii/no,nv,nd,nd2
2642	integer ndc,ndc2,ndpt,ndt
2643	common /coast/ndc,ndc2,ndt,ndpt
2644	dimension cr(ndim2,ndim2),st(ndim),q(ndim),a(ndim)
2645	dimension sa(ndim2,ndim2),sai(ndim2,ndim2),cm(ndim2,ndim2)
2646	integer c(),a2(),a2i(*)
2647	logical*1 mapflol
2648
2649	!*DAEXT(NO,NV) C(NDIM2),A2(NDIM2),A2I(NDIM2)
2650	x2pi=datan(1.d0)*8.d0
2651
2652	do i=1,ntt
2653	jx(i)=0
2654	enddo
2655
2656	! frank/etienne
2657	do i=1,ndim
2658	st(i)=0d0
2659	q(i)=0d0
2660	a(i)=0d0
2661	enddo
2662	! frank/etienne
2663	do i=1,ndim2
2664	! frank/etienne
2665	do j=1,ndim2
2666	sai(i,j)=0.d0
2667	sa(i,j)=0.d0
2668	cm(i,j)=0.d0
2669	cr(i,j)=0.d0
2670	enddo
2671	enddo
2672
2673	do i=1,nd2
2674	do j=1,nd2
2675	jx(j)=1
2676	call dapek(c(i),jx,r)
2677	jx(j)=0
2678	cm(i,j)=r
2679	enddo
2680	enddo
2681
2682	midbflo = mapflol(sa,sai,cr,cm,st)
2683	do i=1,nd-ndc
2684	if(st(i)+0.001.gt.1.d0) then
2685	a(i)=dsqrt(cr(2i-1,2i-1)*2+cr(2i-1,2i)*2)
2686	q(i)=dacos(cr(2i-1,2i-1)/a(i))
2687	a(i)=dlog(a(i))
2688	if(cr(2i-1,2i).lt.0.d0) q(i)=x2pi-q(i)
2689	else
2690	a(i)=dsqrt(cr(2i-1,2i-1)*2-cr(2i-1,2i)*2)
2691	ch=cr(2i-1,2i-1)/a(i)
2692	shm=cr(2i-1,2i)/a(i)
2693	! CH=CH+DSQRT(CH**2-1.D0)
2694	! q(i)=DLOG(CH)
2695	q(i)=-dlog(ch+shm)
2696	! IF(cr(2i-1,2i).gt.0.d0) Q(I)=-Q(I)
2697	a(i)=dlog(a(i))
2698	endif
2699	enddo
2700
2701	if(ndc.eq.0) then
2702	if(st(3)+0.001.gt.1.d0.and.nd.eq.3.and.q(nd).gt.0.5d0) &
2703	& q(3)=q(3)-x2pi
2704	else
2705	q(nd)=cr(ndt,ndpt)
2706	endif
2707
2708	call daclrd(a2)
2709	call daclrd(a2i)
2710
2711	do i=1,nd2
2712	do j=1,nd2
2713	jx(j)=1
2714	r=sa(i,j)
2715	if(r.ne.0.d0)call dapok(a2(i),jx,r)
2716	jx(j)=1
2717	r=sai(i,j)
2718	if(r.ne.0.d0)call dapok(a2i(i),jx,r)
2719	jx(j)=0
2720	enddo
2721	enddo
2722
2723	return
2724	end
2725
2726
2727	logical*1 function mapflol(sa,sai,cr,cm,st)
2728	implicit none
2729	integer i,ier,iunst,j,l,n,n1,ndim,ndim2
2730	double precision ap,ax,cm,cr, &
2731	&p,rd,rd1,ri,rr,s1,sa,sai,st,vi,vr,w,x,x2pi,xd,xj,xsu,xx
2732	parameter (ndim2=6)
2733	parameter (ndim=3)
2734	integer nd,nd2,no,nv
2735	common /ii/no,nv,nd,nd2
2736	integer ndc,ndc2,ndpt,ndt
2737	common /coast/ndc,ndc2,ndt,ndpt
2738	!---- FROM TRACKING CODE
2739	integer idpr
2740	common /printing/ idpr
2741	integer nplane
2742	double precision epsplane,xplane
2743	common /choice/ xplane(ndim),epsplane,nplane(ndim)
2744	! ---------------------
2745	dimension cr(ndim2,ndim2),xj(ndim2,ndim2),n(ndim),x(ndim)
2746	dimension rr(ndim2),ri(ndim2),sa(ndim2,ndim2),xx(ndim) &
2747	&,sai(ndim2,ndim2),cm(ndim2,ndim2),w(ndim2,ndim2),st(ndim)
2748	dimension vr(ndim2,ndim2),vi(ndim2,ndim2),s1(ndim2,ndim2),p(ndim2)
2749	logical*1 eig6
2750
2751	x2pi=datan(1.d0)*8.d0
2752	n1=0
2753	! frank/etienne
2754	do i=1,ndim2
2755	do j=1,ndim2
2756	cr(j,i)=cm(i,j)
2757	xj(i,j)=0.d0
2758	s1(i,j)=0.d0
2759	enddo
2760	enddo
2761	! frank/etienne
2762	do i=1,ndim
2763	n(i)=0
2764	xj(2i-1,2i)=1.d0
2765	xj(2i,2i-1)=-1.d0
2766	enddo
2767	! frank/etienne
2768	do i=1,ndim2
2769	do j=1,ndim2
2770	sai(i,j)=0.d0
2771	w(i,j)=cm(i,j)
2772	enddo
2773	enddo
2774	if(ndc.eq.1) then
2775	s1(nd2-ndc,nd2-ndc)=1.d0
2776	s1(nd2,nd2)=1.d0
2777	sai(nd2-ndc,nd2-ndc)=1.d0
2778	sai(nd2,nd2)=1.d0
2779	endif
2780	call mulnd2(xj,w)
2781	call mulnd2(cr,w)
2782	if(idpr.ge.0.or.idpr.eq.-102) then
2783	write(6,*)'Check of the symplectic condition on the linear part'
2784	xsu=0.d0
2785	do i=1,nd2
2786	write(6,'(6(2x,g23.16))') ( w(i,j), j = 1, nd2 )
2787	do j=1,nd2
2788	xsu=xsu+dabs(w(i,j))
2789	enddo
2790	enddo
2791	write(6,)'deviation for symplecticity ',100.d0(xsu-nd2)/xsu, &
2792	& ' %'
2793	endif
2794	mapflol = eig6(cr,rr,ri,vr,vi)
2795	if(idpr.ge.0) then
2796	write(6,*) ' '
2797	write(6,*) ' Index Real Part ', &
2798	& ' ArcSin(Imaginary Part)/2/pi'
2799	write(6,*) ' '
2800	do i=1,nd-ndc
2801	rd1=dsqrt(rr(2i-1)2+ri(2i-1)**2)
2802	rd=dsqrt(rr(2i)2+ri(2i)**2)
2803	write(6,) 2i-1,rr(2i-1),dasin(ri(2i-1)/rd1)/x2pi
2804	write(6,) 2i,rr(2i),dasin(ri(2i)/rd)/x2pi
2805	write(6,) ' alphas ', dlog(dsqrt(rdrd1))
2806	enddo
2807	if ( idpr.ge. 0) then
2808	write(6,*) &
2809	& ' select ',nd-ndc, &
2810	& ' eigenplanes (odd integers <0 real axis)'
2811	! read(5,*) (n(i),i=1,nd-ndc)
2812	n(1) = 1
2813	n(2) = 3
2814	n(3) = 5
2815	else
2816	n(1) = 1
2817	n(2) = 3
2818	n(3) = 5
2819	endif
2820	elseif(idpr.eq.-100) then
2821	do i=1,nd-ndc
2822	n(i)=nplane(i)
2823	enddo
2824	elseif(idpr.eq.-101.or.idpr.eq.-102) then
2825	do i=1,nd-ndc
2826	if(ri(2*i).ne.0.d0) then
2827	n(i)=2*i-1
2828	else
2829	n(i)=-2*i+1
2830	endif
2831	enddo
2832	else
2833	do i=1,nd-ndc
2834	n(i)=2*i-1
2835	enddo
2836	endif
2837	iunst=0
2838	do i=1,nd-ndc ! Frank NDC kept
2839	if(n(i).lt.0) then
2840	n(i)=-n(i)
2841	st(i)=0.d0
2842	iunst=1
2843	else
2844	st(i)=1.d0
2845	endif
2846	x(i)=0.d0
2847	xx(i)=1.d0
2848	do j=1,nd-ndc
2849	x(i)=vr(2j-1,n(i))vi(2j,n(i))-vr(2j,n(i))vi(2j-1,n(i))+ &
2850	& x(i)
2851	enddo
2852	enddo
2853
2854	do i=1,nd-ndc
2855	if(x(i).lt.0.d0) xx(i)=-1.d0
2856	x(i)=dsqrt(dabs(x(i)))
2857	enddo
2858	do i=1,nd2-ndc2
2859	do j=1,nd-ndc
2860	if(st(j)+0.001.gt.1.d0) then
2861	sai(2j-1,i)=vr(i,n(j))xx(j)/x(j)
2862	sai(2*j,i)=vi(i,n(j))/x(j)
2863	else
2864	ax=vr(i,n(j))*xx(j)/x(j)
2865	ap=vi(i,n(j))/x(j)
2866	sai(2*j-1,i)=(ax+ap)/dsqrt(2.d0)
2867	sai(2*j,i)=(ap-ax)/dsqrt(2.d0)
2868	endif
2869	enddo
2870	enddo
2871	if(idpr.eq.-101.or.idpr.eq.-102) then
2872	call movearou(sai)
2873	endif
2874	! adjust sa such that sa(1,2)=0 and sa(3,4)=0. (courant-snyder-edwards-teng
2875	! phase advances)
2876	if(iunst.ne.1) then
2877	do i=1,nd-ndc
2878	p(i)=datan(-sai(2i-1,2i)/sai(2i,2i))
2879	s1(2i-1,2i-1)=dcos(p(i))
2880	s1(2i,2i)=dcos(p(i))
2881	s1(2i-1,2i)=dsin(p(i))
2882	s1(2i,2i-1)=-dsin(p(i))
2883	enddo
2884	call mulnd2(s1,sai)
2885	! adjust sa to have sa(1,1)>0 and sa(3,3)>0 rotate by pi if necessary.
2886	do i=1,nd-ndc
2887	xd=1.d0
2888	if(sai(2i-1,2i-1).lt.0.d0) xd=-1.d0
2889	s1(2i-1,2i-1)=xd
2890	s1(2i-1,2i)=0.d0
2891	s1(2i,2i-1)=0.d0
2892	s1(2i,2i)=xd
2893	enddo
2894	call mulnd2(s1,sai)
2895	! sa is now uniquely and unambigeously determined.
2896	endif
2897	do i=1,nd2
2898	do l=1,nd2
2899	sa(i,l)=sai(i,l)
2900	enddo
2901	enddo
2902	call matinv(sai,sa,nd2,6,ier)
2903
2904	call mulnd2(sai,cm)
2905	do i=1,nd2
2906	do j=1,nd2
2907	cr(i,j)=sa(i,j)
2908	enddo
2909	enddo
2910
2911	call mulnd2(cm,cr)
2912
2913	return
2914	end
2915	subroutine mulnd2(rt,r)
2916	implicit none
2917	integer i,ia,j,ndim,ndim2
2918	double precision r,rt,rtt
2919	parameter (ndim2=6)
2920	parameter (ndim=3)
2921	integer nd,nd2,no,nv
2922	common /ii/no,nv,nd,nd2
2923	dimension rt(ndim2,ndim2),r(ndim2,ndim2),rtt(ndim2,ndim2)
2924	do i=1,nd2
2925	do j=1,nd2
2926	rtt(i,j)=0.d0
2927	enddo
2928	enddo
2929	do i=1,nd2
2930	do j=1,nd2
2931	do ia=1,nd2
2932	rtt(i,ia)=rt(i,j)*r(j,ia)+rtt(i,ia)
2933	enddo
2934	enddo
2935	enddo
2936
2937	do i=1,nd2
2938	do j=1,nd2
2939	r(i,j)=rtt(i,j)
2940	enddo
2941	enddo
2942	return
2943	end
2944	subroutine movearou(rt)
2945	implicit none
2946	integer ipause, mypause
2947	integer i,ic,j,ndim,ndim2
2948	double precision rt,rto,s,xr,xrold,xy,xyz,xz,xzy,yz
2949	parameter (ndim2=6)
2950	parameter (ndim=3)
2951	integer nd,nd2,no,nv
2952	common /ii/no,nv,nd,nd2
2953	integer idpr
2954	common /printing/ idpr
2955	dimension rt(ndim2,ndim2),rto(ndim2,ndim2)
2956	dimension xy(ndim2,ndim2),xz(ndim2,ndim2),yz(ndim2,ndim2)
2957	dimension xyz(ndim2,ndim2),xzy(ndim2,ndim2)
2958	dimension s(ndim2,ndim2)
2959	do i=1,nd2
2960	do j=1,nd2
2961	s(i,j)=0.d0
2962	s(i,i)=1.d0
2963	xy(i,j)=0.d0
2964	xz(i,j)=0.d0
2965	yz(i,j)=0.d0
2966	xyz(i,j)=0.d0
2967	xzy(i,j)=0.d0
2968	enddo
2969	enddo
2970
2971	xy(1,3)=1.d0
2972	xy(3,1)=1.d0
2973	xy(2,4)=1.d0
2974	xy(4,2)=1.d0
2975	xy(5,5)=1.d0
2976	xy(6,6)=1.d0
2977
2978	xz(1,5)=1.d0
2979	xz(5,1)=1.d0
2980	xz(2,6)=1.d0
2981	xz(6,2)=1.d0
2982	xz(3,3)=1.d0
2983	xz(4,4)=1.d0
2984
2985	yz(3,5)=1.d0
2986	yz(5,3)=1.d0
2987	yz(4,6)=1.d0
2988	yz(6,4)=1.d0
2989	yz(1,1)=1.d0
2990	yz(2,2)=1.d0
2991
2992	xyz(1,3)=1.d0
2993	xyz(3,5)=1.d0
2994	xyz(5,1)=1.d0
2995	xyz(2,4)=1.d0
2996	xyz(4,6)=1.d0
2997	xyz(6,2)=1.d0
2998
2999	xzy(1,5)=1.d0
3000	xzy(5,3)=1.d0
3001	xzy(3,1)=1.d0
3002	xzy(2,6)=1.d0
3003	xzy(6,4)=1.d0
3004	xzy(4,2)=1.d0
3005
3006	ic=0
3007	xrold=1000000000.d0
3008	call movemul(rt,s,rto,xr)
3009	! write(6,*) xr,xrold
3010	! do i=1,6
3011	! write(6,'(6(1x,1pe12.5))') (RTO(i,j),j=1,6)
3012	! enddo
3013	! ipause=mypause(0)
3014	if(xr.lt.xrold) then
3015	xrold=xr
3016	endif
3017
3018	if(nd.ge.2) then
3019	call movemul(rt,xy,rto,xr)
3020	if(xr.lt.xrold) then
3021	xrold=xr
3022	ic=1
3023	endif
3024	endif
3025
3026	if(nd.eq.3) then
3027	call movemul(rt,xz,rto,xr)
3028	if(xr.lt.xrold) then
3029	xrold=xr
3030	ic=2
3031	endif
3032	call movemul(rt,yz,rto,xr)
3033	if(xr.lt.xrold) then
3034	xrold=xr
3035	ic=3
3036	endif
3037	call movemul(rt,xyz,rto,xr)
3038	if(xr.lt.xrold) then
3039	xrold=xr
3040	ic=4
3041	endif
3042	call movemul(rt,xzy,rto,xr)
3043	if(xr.lt.xrold) then
3044	xrold=xr
3045	ic=5
3046	endif
3047	endif
3048
3049	if(ic.eq.0) then
3050	call movemul(rt,s,rto,xr)
3051	if(idpr.gt.-101) write(6,*) " no exchanged"
3052	elseif(ic.eq.1) then
3053	call movemul(rt,xy,rto,xr)
3054	if(idpr.gt.-101) write(6,*) " x-y exchanged"
3055	elseif(ic.eq.2) then
3056	call movemul(rt,xz,rto,xr)
3057	if(idpr.gt.-101) write(6,*) " x-z exchanged"
3058	elseif(ic.eq.3) then
3059	call movemul(rt,yz,rto,xr)
3060	if(idpr.gt.-101) write(6,*) " y-z exchanged"
3061	elseif(ic.eq.4) then
3062	call movemul(rt,xyz,rto,xr)
3063	if(idpr.gt.-101) write(6,*) " x-y-z permuted"
3064	elseif(ic.eq.5) then
3065	call movemul(rt,xzy,rto,xr)
3066	if(idpr.gt.-101) write(6,*) " x-z-y permuted"
3067	endif
3068
3069	do i=1,nd2
3070	do j=1,nd2
3071	rt(i,j)=rto(i,j)
3072	enddo
3073	enddo
3074
3075	return
3076	end
3077	subroutine movemul(rt,xy,rto,xr)
3078	implicit none
3079	integer i,j,k,ndim,ndim2
3080	double precision rt,rto,xr,xy
3081	parameter (ndim2=6)
3082	parameter (ndim=3)
3083	integer nd,nd2,no,nv
3084	common /ii/no,nv,nd,nd2
3085	dimension rt(ndim2,ndim2)
3086	dimension xy(ndim2,ndim2),rto(ndim2,ndim2)
3087
3088	do i=1,nd2
3089	do j=1,nd2
3090	rto(i,j)=0.d0
3091	enddo
3092	enddo
3093
3094	do i=1,nd2
3095	do j=1,nd2
3096	do k=1,nd2
3097	rto(i,k)=xy(i,j)*rt(j,k)+rto(i,k)
3098	enddo
3099	enddo
3100	enddo
3101
3102	xr=0.d0
3103	do i=1,nd2
3104	do j=1,nd2
3105	xr=xr+dabs(rto(i,j))
3106	enddo
3107	enddo
3108	do i=1,nd
3109	xr=xr-dabs(rto(2i-1,2i-1))
3110	xr=xr-dabs(rto(2i-1,2i))
3111	xr=xr-dabs(rto(2i,2i))
3112	xr=xr-dabs(rto(2i,2i-1))
3113	enddo
3114	return
3115	end
3116	subroutine initpert(st,ang,ra)
3117	implicit none
3118	integer i,ndim,ndim2,nn,nreso
3119	double precision ang,ra,st
3120	! X-RATED
3121	!- SETS UP ALL THE COMMON BLOCKS RELEVANT TO NORMAL FORM AND THE BASIS
3122	!- CHANGES INSIDE MAPNORMF
3123	parameter (ndim=3)
3124	parameter (ndim2=6)
3125	parameter (nreso=20)
3126	dimension st(ndim),ang(ndim),ra(ndim)
3127	double precision angle,dsta,rad,sta
3128	common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
3129	integer idsta,ista
3130	common /istable/ista(ndim),idsta(ndim)
3131	integer nd,nd2,no,nv
3132	common /ii/no,nv,nd,nd2
3133	integer ndc,ndc2,ndpt,ndt
3134	common /coast/ndc,ndc2,ndt,ndpt
3135	integer mx,nres
3136	common /reson/mx(ndim,nreso),nres
3137	integer iref
3138	common /resfile/iref
3139
3140	if(iref.gt.0) then
3141	write(6,*) iref
3142	read(iref,*) nres
3143	if(nres.ge.nreso) then
3144	write(6,*) ' NRESO IN LIELIB TOO SMALL '
3145	stop999
3146	endif
3147	elseif(iref.eq.0) then
3148	nres=0
3149	endif
3150	if(nres.ne.0) write(6,*)' warning resonances left in the map'
3151	if(iref.gt.0) then
3152	do i=1,nres
3153	read(iref,*) (mx(nn,i),nn=1,nd-ndc)
3154	enddo
3155	endif
3156	do i=nres+1,nreso
3157	do nn=1,ndim
3158	mx(nn,i)=0
3159	enddo
3160	enddo
3161	! frank/Etienne
3162	do i=1,ndim
3163	angle(i)=0.d0
3164	rad(i)=0.d0
3165	sta(i)=0.d0
3166	dsta(i)=1.d0-sta(i)
3167	ista(i)=0
3168	idsta(i)=0
3169	enddo
3170	do i=1,nd ! Frank -ndc
3171	angle(i)=ang(i)
3172	rad(i)=ra(i)
3173	sta(i)=st(i)
3174	dsta(i)=1.d0-sta(i)
3175	enddo
3176	do i=1,nd
3177	ista(i)=idint(sta(i)+.01)
3178	idsta(i)=idint(dsta(i)+.01)
3179	enddo
3180	return
3181	end
3182	double precision function dlie(j)
3183	implicit none
3184	integer i,ndim
3185	parameter (ndim=3)
3186	! PARAMETER (NTT=40)
3187	! INTEGER J(NTT)
3188	integer j(*)
3189	integer nd,nd2,no,nv
3190	common /ii/no,nv,nd,nd2
3191	dlie=0.d0
3192	do i=1,nd
3193	dlie=dble(j(2i-1)+j(2i))+dlie
3194	enddo
3195	dlie=dlie+1.d0
3196	dlie=1.d0/dlie
3197	return
3198	end
3199	double precision function rext(j)
3200	implicit none
3201	integer i,lie,mo,ndim
3202	parameter (ndim=3)
3203	! PARAMETER (NTT=40)
3204	integer nd,nd2,no,nv
3205	common /ii/no,nv,nd,nd2
3206	integer ndc,ndc2,ndpt,ndt
3207	common /coast/ndc,ndc2,ndt,ndpt
3208	integer idsta,ista
3209	common /istable/ista(ndim),idsta(ndim)
3210	integer j(*)
3211	lie=0
3212	do i=1,nd-ndc
3213	lie=ista(i)j(2i)+lie
3214	enddo
3215	mo=mod(lie,4)+1
3216	goto(11,12,13,14),mo
3217	11 rext = 1.d0
3218	return
3219	12 rext = -1.d0
3220	return
3221	13 rext = -1.d0
3222	return
3223	14 rext = 1.d0
3224	return
3225	end
3226	subroutine cpart(h,ch)
3227	implicit none
3228	integer ndim,ntt
3229	double precision rext
3230	parameter (ndim=3)
3231	parameter (ntt=40)
3232	external rext
3233	integer nd,nd2,no,nv
3234	common /ii/no,nv,nd,nd2
3235	integer h,ch
3236	call dacfu(h,rext,ch)
3237	return
3238	end
3239	subroutine ctoi(f1,f2)
3240	implicit none
3241	integer ndim2,ntt
3242	parameter (ndim2=6)
3243	parameter (ntt=40)
3244	integer nd,nd2,no,nv
3245	common /ii/no,nv,nd,nd2
3246	integer f1,f2
3247	integer b1,x(ndim2)
3248	!
3249	!
3250	call etallnom(b1,1,'B1 ')
3251	call etallnom(x,nd2 ,'X ')
3252
3253	call cpart(f1,b1)
3254	call etctr(x)
3255	call trx(b1,f2,x)
3256	call dadal(x,nd2)
3257	call dadal(b1,1)
3258	return
3259	end
3260	subroutine itoc(f1,f2)
3261	implicit none
3262	integer ndim2,ntt
3263	parameter (ndim2=6)
3264	parameter (ntt=40)
3265	integer nd,nd2,no,nv
3266	common /ii/no,nv,nd,nd2
3267	integer f1,f2
3268	integer b1,x(ndim2)
3269	!
3270	call etallnom(b1,1,'B1 ')
3271	call etallnom(x,nd2 ,'X ')
3272
3273	call etrtc(x)
3274	call trx(f1,b1,x)
3275	call cpart(b1,f2)
3276	call dadal(x,nd2)
3277	call dadal(b1,1)
3278	return
3279	end
3280	subroutine etrtc(x)
3281	implicit none
3282	integer i,ndim,ndim2,ntt
3283	parameter (ndim=3)
3284	parameter (ndim2=6)
3285	parameter (ntt=40)
3286	integer nd,nd2,no,nv
3287	common /ii/no,nv,nd,nd2
3288	integer ndc,ndc2,ndpt,ndt
3289	common /coast/ndc,ndc2,ndt,ndpt
3290	integer x(*)
3291
3292	integer rel(ndim2)
3293	!
3294	!
3295	call etallnom(rel,nd2 ,'REL ')
3296
3297	call etini(rel)
3298	call etini(x)
3299	do i=1,nd-ndc
3300	call daadd(rel(2i-1),rel(2i),x(2*i-1))
3301	call dasub(rel(2i-1),rel(2i),x(2*i))
3302	enddo
3303	call dadal(rel,nd2)
3304	return
3305	end
3306	subroutine etctr(x)
3307	implicit none
3308	integer i,ndim,ndim2,ntt
3309	parameter (ndim=3)
3310	parameter (ndim2=6)
3311	parameter (ntt=40)
3312	integer nd,nd2,no,nv
3313	common /ii/no,nv,nd,nd2
3314	integer ndc,ndc2,ndpt,ndt
3315	common /coast/ndc,ndc2,ndt,ndpt
3316	integer x(*)
3317	integer rel(ndim2)
3318	!
3319	!
3320	call etallnom(rel,nd2 ,'REL ')
3321
3322	call etini(rel)
3323	call etini(x)
3324	do i=1,nd-ndc
3325	call dalin(rel(2i-1),.5d0,rel(2i),.5d0,x(2*i-1))
3326	call dalin(rel(2i-1),.5d0,rel(2i),-.5d0,x(2*i))
3327	enddo
3328	call dadal(rel,nd2)
3329	return
3330	end
3331	subroutine etcjg(x)
3332	implicit none
3333	integer i,ndim,ndim2,ntt
3334	parameter (ndim=3)
3335	parameter (ndim2=6)
3336	parameter (ntt=40)
3337	integer idsta,ista
3338	common /istable/ista(ndim),idsta(ndim)
3339	integer nd,nd2,no,nv
3340	common /ii/no,nv,nd,nd2
3341	integer ndc,ndc2,ndpt,ndt
3342	common /coast/ndc,ndc2,ndt,ndpt
3343	integer x(*)
3344
3345	integer rel(ndim2)
3346	!
3347	!
3348	call etallnom(rel,nd2 ,'REL ')
3349
3350	call etini(rel)
3351	call etini(x)
3352	do i=1,nd-ndc
3353	if(ista(i).eq.1) then
3354	call dacop(rel(2i-1),x(2i))
3355	call dacop(rel(2i),x(2i-1))
3356	else
3357	call dacop(rel(2i-1),x(2i-1))
3358	call dacop(rel(2i),x(2i))
3359	endif
3360	enddo
3361	call dadal(rel,nd2)
3362	return
3363	end
3364
3365
3366	logical*1 function eig6(fm,reval,aieval,revec,aievec)
3367	implicit none
3368	integer jet,ndim2
3369	!**************************************************************************
3370
3371	! Diagonalization routines of NERI
3372
3373	!ccccccccccccccccc
3374	!
3375	! this routine finds the eigenvalues and eigenvectors
3376	! of the full matrix fm.
3377	! the eigenvectors are normalized so that the real and
3378	! imaginary part of vectors 1, 3, and 5 have +1 antisymmetric
3379	! product:
3380	! revec1 J aivec1 = 1 ; revec3 J aivec3 = 1 ;
3381	! revec5 J aivec5 = 1.
3382	! the eigenvectors 2 ,4, and 6 have the opposite normalization.
3383	! written by F. Neri, Feb 26 1986.
3384	!
3385	parameter (ndim2=6)
3386	integer nn
3387	integer ilo,ihi,mdim,info
3388	integer nd,nd2,no,nv
3389	common /ii/no,nv,nd,nd2
3390	integer ndc,ndc2,ndpt,ndt
3391	common /coast/ndc,ndc2,ndt,ndpt
3392	double precision reval(ndim2),aieval(ndim2), &
3393	&revec(ndim2,ndim2),aievec(ndim2,ndim2)
3394	double precision fm(ndim2,ndim2),aa(ndim2,ndim2)
3395	integer i,i1
3396	double precision ort(ndim2),vv(ndim2,ndim2)
3397
3398	eig6 = .true.
3399	! copy matrix to temporary storage (the matrix aa is destroyed)
3400	do i=1,nd2-ndc2
3401	do i1=1,nd2-ndc2
3402	aa(i1,i) = fm(i1,i)
3403	enddo
3404	enddo
3405	ilo = 1
3406	ihi = nd2-ndc2
3407	mdim = ndim2
3408	nn = nd2-ndc2
3409	! compute eigenvalues and eigenvectors using double
3410	! precision Eispack routines:
3411	call ety(mdim,nn,ilo,ihi,aa,ort)
3412	call etyt(mdim,nn,ilo,ihi,aa,ort,vv)
3413	call ety2(mdim,nn,ilo,ihi,aa,reval,aieval,vv,info)
3414	if ( info .ne. 0 ) then
3415	write(6,*) ' ERROR IN EIG6'
3416	eig6 = .false.
3417	return
3418	endif
3419	! call neigv(vv,pbkt)
3420	do i=1,nd-ndc
3421	do jet=1,nd2-ndc2
3422	revec(jet,2i-1)=vv(jet,2i-1)
3423	revec(jet,2i)=vv(jet,2i-1)
3424	aievec(jet,2i-1)=vv(jet,2i)
3425	aievec(jet,2i)=-vv(jet,2i)
3426	enddo
3427	enddo
3428	do i=1,nd2-ndc2
3429	if(dabs(reval(i)2+aieval(i)2 -1.d0).gt.1.d-10) then
3430	write(6,*) ' EIG6: Eigenvalues off the unit circle!'
3431	eig6 = .false.
3432	endif
3433	enddo
3434	return
3435	end
3436
3437
3438	subroutine ety(nm,n,low,igh,a,ort)
3439	implicit none
3440	integer i,j,m,n,ii,jj,la,mp,nm,igh,kp1,low
3441	double precision a(nm,n),ort(igh)
3442	double precision f,g,h,scale
3443	!
3444	! this subroutine is a translation of the algol procedure orthes,
3445	! num. math. 12, 349-368(1968) by martin and wilkinson.
3446	! handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).
3447	!
3448	! given a real general matrix, this subroutine
3449	! reduces a submatrix situated in rows and columns
3450	! low through igh to upper hessenberg form by
3451	! orthogonal similarity transformations.
3452	!
3453	! on input-
3454	!
3455	! nm must be set to the row dimension of two-dimensional
3456	! array parameters as declared in the calling program
3457	! dimension statement,
3458	!
3459	! n is the order of the matrix,
3460	!
3461	! low and igh are integers determined by the balancing
3462	! subroutine balanc. if balanc has not been used,
3463	! set low=1, igh=n,
3464	!
3465	! a contains the input matrix.
3466	!
3467	! on output-
3468	!
3469	! a contains the hessenberg matrix. information about
3470	! the orthogonal transformations used in the reduction
3471	! is stored in the remaining triangle under the
3472	! hessenberg matrix,
3473	!
3474	! ort contains further information about the transformations.
3475	! only elements low through igh are used.
3476	!
3477	! fortran routine by b. s. garbow
3478	! modified by filippo neri.
3479	!
3480	!
3481	la = igh - 1
3482	kp1 = low + 1
3483	if (la .lt. kp1) go to 200
3484	!
3485	do m = kp1, la
3486	h = 0.0
3487	ort(m) = 0.0
3488	scale = 0.0
3489	! ******** scale column (algol tol then not needed) ********
3490	do i = m, igh
3491	scale = scale + dabs(a(i,m-1))
3492	enddo
3493	!
3494	if (scale .eq. 0.0) go to 180
3495	mp = m + igh
3496	! ******** for i=igh step -1 until m do -- ********
3497	do ii = m, igh
3498	i = mp - ii
3499	ort(i) = a(i,m-1) / scale
3500	h = h + ort(i) * ort(i)
3501	enddo
3502	!
3503	g = -dsign(dsqrt(h),ort(m))
3504	h = h - ort(m) * g
3505	ort(m) = ort(m) - g
3506	! ********** form (i-(uut)/h) a **********
3507	do j = m, n
3508	f = 0.0
3509	! ******** for i=igh step -1 until m do -- ********
3510	do ii = m, igh
3511	i = mp - ii
3512	f = f + ort(i) * a(i,j)
3513	enddo
3514	!
3515	f = f / h
3516	!
3517	do i = m, igh
3518	a(i,j) = a(i,j) - f * ort(i)
3519	enddo
3520	!
3521	enddo
3522	! ********** form (i-(uut)/h)a(i-(uut)/h) **********
3523	do i = 1, igh
3524	f = 0.0
3525	! ******** for j=igh step -1 until m do -- ********
3526	do jj = m, igh
3527	j = mp - jj
3528	f = f + ort(j) * a(i,j)
3529	enddo
3530	!
3531	f = f / h
3532	!
3533	do j = m, igh
3534	a(i,j) = a(i,j) - f * ort(j)
3535	enddo
3536	!
3537	enddo
3538	!
3539	ort(m) = scale * ort(m)
3540	a(m,m-1) = scale * g
3541	180 continue
3542	enddo
3543	!
3544	200 return
3545	! ******** last card of ety ********
3546	end
3547	subroutine etyt(nm,n,low,igh,a,ort,z)
3548	implicit none
3549	integer i,j,n,kl,mm,mp,nm,igh,low,mp1
3550	double precision a(nm,igh),ort(igh),z(nm,n)
3551	double precision g
3552	!
3553	! this subroutine is a translation of the algol procedure ortrans,
3554	! num. math. 16, 181-204(1970) by peters and wilkinson.
3555	! handbook for auto. comp., vol.ii-linear algebra, 372-395(1971).
3556	!
3557	! this subroutine accumulates the orthogonal similarity
3558	! transformations used in the reduction of a real general
3559	! matrix to upper hessenberg form by ety.
3560	!
3561	! on input-
3562	!
3563	! nm must be set to the row dimension of two-dimensional
3564	! array parameters as declared in the calling program
3565	! dimension statement,
3566	!
3567	! n is the order of the matrix,
3568	!
3569	! low and igh are integers determined by the balancing
3570	! subroutine balanc. if balanc has not been used,
3571	! set low=1, igh=n,
3572	!
3573	! a contains information about the orthogonal trans-
3574	! formations used in the reduction by orthes
3575	! in its strict lower triangle,
3576	!
3577	! ort contains further information about the trans-
3578	! formations used in the reduction by ety.
3579	! only elements low through igh are used.
3580	!
3581	! on output-
3582	!
3583	! z contains the transformation matrix produced in the
3584	! reduction by ety,
3585	!
3586	! ort has been altered.
3587	!
3588	! fortran routine by b. s. garbow.
3589	! modified by f. neri.
3590	!
3591	!
3592	! ******** initialize z to identity matrix ********
3593	do i = 1, n
3594	!
3595	do j = 1, n
3596	z(i,j) = 0.0
3597	enddo
3598	!
3599	z(i,i) = 1.0
3600	enddo
3601	!
3602	kl = igh - low - 1
3603	if (kl .lt. 1) go to 200
3604	! ******** for mp=igh-1 step -1 until low+1 do -- ********
3605	do mm = 1, kl
3606	mp = igh - mm
3607	if (a(mp,mp-1) .eq. 0.0) go to 140
3608	mp1 = mp + 1
3609	!
3610	do i = mp1, igh
3611	ort(i) = a(i,mp-1)
3612	enddo
3613	!
3614	do j = mp, igh
3615	g = 0.0
3616	!
3617	do i = mp, igh
3618	g = g + ort(i) * z(i,j)
3619	enddo
3620	! ********** divisor below is negative of h formed in orthes.
3621	! double division avoids possible underflow **********
3622	g = (g / ort(mp)) / a(mp,mp-1)
3623	!
3624	do i = mp, igh
3625	z(i,j) = z(i,j) + g * ort(i)
3626	enddo
3627	!
3628	enddo
3629	!
3630	140 continue
3631	enddo
3632	!
3633	200 return
3634	! ******** last card of etyt ********
3635	end
3636	subroutine ety2(nm,n,low,igh,h,wr,wi,z,ierr)
3637	implicit none
3638	integer i,j,k,l,m,n,en,ii,jj,ll,mm,na,nm,nn, &
3639	&igh,its,low,mp2,enm2,ierr
3640	double precision h(nm,n),wr(n),wi(n),z(nm,n)
3641	double precision p,q,r,s,t,w,x,y,ra,sa,vi,vr,zz,norm,machep
3642	logical notlas
3643	double precision z3r,z3i
3644	!
3645	!
3646	!
3647	! this subroutine is a translation of the algol procedure hqr2,
3648	! num. math. 16, 181-204(1970) by peters and wilkinson.
3649	! handbook for auto. comp., vol.ii-linear algebra, 372-395(1971).
3650	!
3651	! this subroutine finds the eigenvalues and eigenvectors
3652	! of a real upper hessenberg matrix by the qr method. the
3653	! eigenvectors of a real general matrix can also be found
3654	! if elmhes and eltran or orthes and ortran have
3655	! been used to reduce this general matrix to hessenberg form
3656	! and to accumulate the similarity transformations.
3657	!
3658	! on input-
3659	!
3660	! nm must be set to the row dimension of two-dimensional
3661	! array parameters as declared in the calling program
3662	! dimension statement,
3663	!
3664	! n is the order of the matrix,
3665	!
3666	! low and igh are integers determined by the balancing
3667	! subroutine balanc. if balanc has not been used,
3668	! set low=1, igh=n,
3669	!
3670	! h contains the upper hessenberg matrix,
3671	!
3672	! z contains the transformation matrix produced by eltran
3673	! after the reduction by elmhes, or by ortran after the
3674	! reduction by orthes, if performed. if the eigenvectors
3675	! of the hessenberg matrix are desired, z must contain the
3676	! identity matrix.
3677	!
3678	! on output-
3679	!
3680	! h has been destroyed,
3681	!
3682	! wr and wi contain the real and imaginary parts,
3683	! respectively, of the eigenvalues. the eigenvalues
3684	! are unordered except that complex conjugate pairs
3685	! of values appear consecutively with the eigenvalue
3686	! having the positive imaginary part first. if an
3687	! error exit is made, the eigenvalues should be correct
3688	! for indices ierr+1,...,n,
3689	!
3690	! z contains the real and imaginary parts of the eigenvectors.
3691	! if the i-th eigenvalue is real, the i-th column of z
3692	! contains its eigenvector. if the i-th eigenvalue is complex
3693	! with positive imaginary part, the i-th and (i+1)-th
3694	! columns of z contain the real and imaginary parts of its
3695	! eigenvector. the eigenvectors are unnormalized. if an
3696	! error exit is made, none of the eigenvectors has been found,
3697	!
3698	! ierr is set to
3699	! zero for normal return,
3700	! j if the j-th eigenvalue has not been
3701	! determined after 200 iterations.
3702	!
3703	! arithmetic is double precision. complex division
3704	! is simulated by routin etdiv.
3705	!
3706	! fortran routine by b. s. garbow.
3707	! modified by f. neri.
3708	!
3709	!
3710	! ********** machep is a machine dependent parameter specifying
3711	! the relative precision of floating point arithmetic.
3712	!
3713	! **********
3714	machep = 1.d-17
3715	! machep = r1mach(4)
3716	!
3717	ierr = 0
3718	norm = 0.0
3719	k = 1
3720	! ********** store roots isolated by balanc
3721	! and compute matrix norm **********
3722	do i = 1, n
3723	!
3724	do j = k, n
3725	norm = norm + dabs(h(i,j))
3726	enddo
3727	!
3728	k = i
3729	if (i .ge. low .and. i .le. igh) go to 50
3730	wr(i) = h(i,i)
3731	wi(i) = 0.0
3732	50 continue
3733	enddo
3734	!
3735	en = igh
3736	t = 0.0
3737	! ******** search for next eigenvalues ********
3738	60 if (en .lt. low) go to 340
3739	its = 0
3740	na = en - 1
3741	enm2 = na - 1
3742	! ********** look for single small sub-diagonal element
3743	! for l=en step -1 until low do -- **********
3744	70 do ll = low, en
3745	l = en + low - ll
3746	if (l .eq. low) go to 100
3747	s = dabs(h(l-1,l-1)) + dabs(h(l,l))
3748	if (s .eq. 0.0) s = norm
3749	if (dabs(h(l,l-1)) .le. machep * s) go to 100
3750	enddo
3751	! ******** form shift ********
3752	100 x = h(en,en)
3753	if (l .eq. en) go to 270
3754	y = h(na,na)
3755	w = h(en,na) * h(na,en)
3756	if (l .eq. na) go to 280
3757	if (its .eq. 200) go to 1000
3758	if (its .ne. 10 .and. its .ne. 20) go to 130
3759	! ******** form exceptional shift ********
3760	t = t + x
3761	!
3762	do i = low, en
3763	h(i,i) = h(i,i) - x
3764	enddo
3765	!
3766	s = dabs(h(en,na)) + dabs(h(na,enm2))
3767	x = 0.75 * s
3768	y = x
3769	w = -0.4375 * s * s
3770	130 its = its + 1
3771	! ********** look for two consecutive small
3772	! sub-diagonal elements.
3773	! for m=en-2 step -1 until l do -- **********
3774	do mm = l, enm2
3775	m = enm2 + l - mm
3776	zz = h(m,m)
3777	r = x - zz
3778	s = y - zz
3779	p = (r * s - w) / h(m+1,m) + h(m,m+1)
3780	q = h(m+1,m+1) - zz - r - s
3781	r = h(m+2,m+1)
3782	s = dabs(p) + dabs(q) + dabs(r)
3783	p = p / s
3784	q = q / s
3785	r = r / s
3786	if (m .eq. l) go to 150
3787	if (dabs(h(m,m-1)) * (dabs(q) + dabs(r)) .le. machep * dabs(p) &
3788	& * (dabs(h(m-1,m-1)) + dabs(zz) + dabs(h(m+1,m+1)))) go to 150
3789	enddo
3790	!
3791	150 mp2 = m + 2
3792	!
3793	do i = mp2, en
3794	h(i,i-2) = 0.0
3795	if (i .eq. mp2) go to 160
3796	h(i,i-3) = 0.0
3797	160 continue
3798	enddo
3799	! ********** double qr step involving rows l to en and
3800	! columns m to en **********
3801	do k = m, na
3802	notlas = k .ne. na
3803	if (k .eq. m) go to 170
3804	p = h(k,k-1)
3805	q = h(k+1,k-1)
3806	r = 0.0
3807	if (notlas) r = h(k+2,k-1)
3808	x = dabs(p) + dabs(q) + dabs(r)
3809	if (x .eq. 0.0) go to 260
3810	p = p / x
3811	q = q / x
3812	r = r / x
3813	170 s = dsign(dsqrt(pp+qq+r*r),p)
3814	if (k .eq. m) go to 180
3815	h(k,k-1) = -s * x
3816	go to 190
3817	180 if (l .ne. m) h(k,k-1) = -h(k,k-1)
3818	190 p = p + s
3819	x = p / s
3820	y = q / s
3821	zz = r / s
3822	q = q / p
3823	r = r / p
3824	! ******** row modification ********
3825	do j = k, n
3826	p = h(k,j) + q * h(k+1,j)
3827	if (.not. notlas) go to 200
3828	p = p + r * h(k+2,j)
3829	h(k+2,j) = h(k+2,j) - p * zz
3830	200 h(k+1,j) = h(k+1,j) - p * y
3831	h(k,j) = h(k,j) - p * x
3832	enddo
3833	!
3834	j = min0(en,k+3)
3835	! ******** column modification ********
3836	do i = 1, j
3837	p = x * h(i,k) + y * h(i,k+1)
3838	if (.not. notlas) go to 220
3839	p = p + zz * h(i,k+2)
3840	h(i,k+2) = h(i,k+2) - p * r
3841	220 h(i,k+1) = h(i,k+1) - p * q
3842	h(i,k) = h(i,k) - p
3843	enddo
3844	! ******** accumulate transformations ********
3845	do i = low, igh
3846	p = x * z(i,k) + y * z(i,k+1)
3847	if (.not. notlas) go to 240
3848	p = p + zz * z(i,k+2)
3849	z(i,k+2) = z(i,k+2) - p * r
3850	240 z(i,k+1) = z(i,k+1) - p * q
3851	z(i,k) = z(i,k) - p
3852	enddo
3853	!
3854	260 continue
3855	enddo
3856	!
3857	go to 70
3858	! ******** one root found ********
3859	270 h(en,en) = x + t
3860	wr(en) = h(en,en)
3861	wi(en) = 0.0
3862	en = na
3863	go to 60
3864	! ******** two roots found ********
3865	280 p = (y - x) / 2.0
3866	q = p * p + w
3867	zz = dsqrt(dabs(q))
3868	h(en,en) = x + t
3869	x = h(en,en)
3870	h(na,na) = y + t
3871	if (q .lt. 0.0) go to 320
3872	! ******** real pair ********
3873	zz = p + dsign(zz,p)
3874	wr(na) = x + zz
3875	wr(en) = wr(na)
3876	if (zz .ne. 0.0) wr(en) = x - w / zz
3877	wi(na) = 0.0
3878	wi(en) = 0.0
3879	x = h(en,na)
3880	s = dabs(x) + dabs(zz)
3881	p = x / s
3882	q = zz / s
3883	r = dsqrt(pp+qq)
3884	p = p / r
3885	q = q / r
3886	! ******** row modification ********
3887	do j = na, n
3888	zz = h(na,j)
3889	h(na,j) = q * zz + p * h(en,j)
3890	h(en,j) = q * h(en,j) - p * zz
3891	enddo
3892	! ******** column modification ********
3893	do i = 1, en
3894	zz = h(i,na)
3895	h(i,na) = q * zz + p * h(i,en)
3896	h(i,en) = q * h(i,en) - p * zz
3897	enddo
3898	! ******** accumulate transformations ********
3899	do i = low, igh
3900	zz = z(i,na)
3901	z(i,na) = q * zz + p * z(i,en)
3902	z(i,en) = q * z(i,en) - p * zz
3903	enddo
3904	!
3905	go to 330
3906	! ******** complex pair ********
3907	320 wr(na) = x + p
3908	wr(en) = x + p
3909	wi(na) = zz
3910	wi(en) = -zz
3911	330 en = enm2
3912	go to 60
3913	! ********** all roots found. backsubstitute to find
3914	! vectors of upper triangular form **********
3915	340 if (norm .eq. 0.0) go to 1001
3916	! ******** for en=n step -1 until 1 do -- ********
3917	do nn = 1, n
3918	en = n + 1 - nn
3919	p = wr(en)
3920	q = wi(en)
3921	na = en - 1
3922	if (q.lt.0) goto 710
3923	if (q.eq.0) goto 600
3924	if (q.gt.0) goto 800
3925	! ******** real vector ********
3926	600 m = en
3927	h(en,en) = 1.0
3928	if (na .eq. 0) go to 800
3929	! ******** for i=en-1 step -1 until 1 do -- ********
3930	do ii = 1, na
3931	i = en - ii
3932	w = h(i,i) - p
3933	r = h(i,en)
3934	if (m .gt. na) go to 620
3935	!
3936	do j = m, na
3937	r = r + h(i,j) * h(j,en)
3938	enddo
3939	!
3940	620 if (wi(i) .ge. 0.0) go to 630
3941	zz = w
3942	s = r
3943	go to 700
3944	630 m = i
3945	if (wi(i) .ne. 0.0) go to 640
3946	t = w
3947	if (w .eq. 0.0) t = machep * norm
3948	h(i,en) = -r / t
3949	go to 700
3950	! ******** solve real equations ********
3951	640 x = h(i,i+1)
3952	y = h(i+1,i)
3953	q = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i)
3954	t = (x * s - zz * r) / q
3955	h(i,en) = t
3956	if (dabs(x) .le. dabs(zz)) go to 650
3957	h(i+1,en) = (-r - w * t) / x
3958	go to 700
3959	650 h(i+1,en) = (-s - y * t) / zz
3960	700 continue
3961	enddo
3962	! ******** end real vector ********
3963	go to 800
3964	! ******** complex vector ********
3965	710 m = na
3966	! ********** last vector component chosen imaginary so that
3967	! eigenvector matrix is triangular **********
3968	if (dabs(h(en,na)) .le. dabs(h(na,en))) go to 720
3969	h(na,na) = q / h(en,na)
3970	h(na,en) = -(h(en,en) - p) / h(en,na)
3971	go to 730
3972	! 720 z3 = cmplx(0.0,-h(na,en)) / cmplx(h(na,na)-p,q)
3973	! h(na,na) = real(z3)
3974	! h(na,en) = aimag(z3)
3975	720 call etdiv(z3r,z3i,0.d0,-h(na,en),h(na,na)-p,q)
3976	h(na,na) = z3r
3977	h(na,en) = z3i
3978	730 h(en,na) = 0.0
3979	h(en,en) = 1.0
3980	enm2 = na - 1
3981	if (enm2 .eq. 0) go to 800
3982	! ******** for i=en-2 step -1 until 1 do -- ********
3983	do ii = 1, enm2
3984	i = na - ii
3985	w = h(i,i) - p
3986	ra = 0.0
3987	sa = h(i,en)
3988	!
3989	do j = m, na
3990	ra = ra + h(i,j) * h(j,na)
3991	sa = sa + h(i,j) * h(j,en)
3992	enddo
3993	!
3994	if (wi(i) .ge. 0.0) go to 770
3995	zz = w
3996	r = ra
3997	s = sa
3998	go to 790
3999	770 m = i
4000	if (wi(i) .ne. 0.0) go to 780
4001	! z3 = cmplx(-ra,-sa) / cmplx(w,q)
4002	! h(i,na) = real(z3)
4003	! h(i,en) = aimag(z3)
4004	call etdiv(z3r,z3i,-ra,-sa,w,q)
4005	h(i,na) = z3r
4006	h(i,en) = z3i
4007	go to 790
4008	! ******** solve complex equations ********
4009	780 x = h(i,i+1)
4010	y = h(i+1,i)
4011	vr = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) - q * q
4012	vi = (wr(i) - p) * 2.0 * q
4013	if (vr .eq. 0.0 .and. vi .eq. 0.0) vr = machep * norm &
4014	& * (dabs(w) + dabs(q) + dabs(x) + dabs(y) + dabs(zz))
4015	! z3 = cmplx(xr-zzra+qsa,xs-zzsa-qra) / cmplx(vr,vi)
4016	! h(i,na) = real(z3)
4017	! h(i,en) = aimag(z3)
4018	call etdiv(z3r,z3i,xr-zzra+qsa,xs-zzsa-qra,vr,vi)
4019	h(i,na) = z3r
4020	h(i,en) = z3i
4021	if (dabs(x) .le. dabs(zz) + dabs(q)) go to 785
4022	h(i+1,na) = (-ra - w * h(i,na) + q * h(i,en)) / x
4023	h(i+1,en) = (-sa - w * h(i,en) - q * h(i,na)) / x
4024	go to 790
4025	! 785 z3 = cmplx(-r-yh(i,na),-s-yh(i,en)) / cmplx(zz,q)
4026	! h(i+1,na) = real(z3)
4027	! h(i+1,en) = aimag(z3)
4028	785 call etdiv(z3r,z3i,-r-yh(i,na),-s-yh(i,en),zz,q)
4029	h(i+1,na) = z3r
4030	h(i+1,en) = z3i
4031	790 continue
4032	enddo
4033	! ******** end complex vector ********
4034	800 continue
4035	enddo
4036	! ********** end back substitution.
4037	! vectors of isolated roots **********
4038	do i = 1, n
4039	if (i .ge. low .and. i .le. igh) go to 840
4040	!
4041	do j = i, n
4042	z(i,j) = h(i,j)
4043	enddo
4044	!
4045	840 continue
4046	enddo
4047	! ********** multiply by transformation matrix to give
4048	! vectors of original full matrix.
4049	! for j=n step -1 until low do -- **********
4050	do jj = low, n
4051	j = n + low - jj
4052	m = min0(j,igh)
4053	!
4054	do i = low, igh
4055	zz = 0.0
4056	!
4057	do k = low, m
4058	zz = zz + z(i,k) * h(k,j)
4059	enddo
4060	!
4061	z(i,j) = zz
4062	enddo
4063	enddo
4064	!
4065	go to 1001
4066	! ********** set error -- no convergence to an
4067	! eigenvalue after 200 iterations **********
4068	1000 ierr = en
4069	1001 return
4070	! ******** last card of ety2 ********
4071	end
4072	subroutine etdiv(a,b,c,d,e,f)
4073	implicit none
4074	! computes the complex division
4075	! a + ib = (c + id)/(e + if)
4076	! very slow, but tries to be as accurate as
4077	! possible by changing the order of the
4078	! operations, so to avoid under(over)flow
4079	! problems.
4080	! Written by F. Neri Feb. 12 1986
4081	!
4082	double precision a,b,c,d,e,f
4083	double precision s,t
4084	double precision cc,dd,ee,ff
4085	double precision temp
4086	integer flip
4087	flip = 0
4088	cc = c
4089	dd = d
4090	ee = e
4091	ff = f
4092	if( dabs(f).ge.dabs(e) ) then
4093	ee = f
4094	ff = e
4095	cc = d
4096	dd = c
4097	flip = 1
4098	endif
4099	s = 1.d0/ee
4100	t = 1.d0/(ee+ ff(ffs))
4101	if ( dabs(ff) .ge. dabs(s) ) then
4102	temp = ff
4103	ff = s
4104	s = temp
4105	endif
4106	if( dabs(dd) .ge. dabs(s) ) then
4107	a = t(cc + s(dd*ff))
4108	else if ( dabs(dd) .ge. dabs(ff) ) then
4109	a = t(cc + dd(s*ff))
4110	else
4111	a = t(cc + ff(s*dd))
4112	endif
4113	if ( dabs(cc) .ge. dabs(s)) then
4114	b = t(dd - s(cc*ff))
4115	else if ( dabs(cc) .ge. dabs(ff)) then
4116	b = t(dd - cc(s*ff))
4117	else
4118	b = t(dd - ff(s*cc))
4119	endif
4120	if (flip.ne.0 ) then
4121	b = -b
4122	endif
4123	return
4124	end
4125	subroutine sympl3(m)
4126	!**********************************************************
4127	!
4128	! SYMPL3
4129	!
4130	!
4131	! On return ,the matrix m(,), supposed to be almost
4132	! symplectic on entry is made exactly symplectic by
4133	! using a non iterative, constructive method.
4134	!
4135	!**********************************************************
4136	!
4137	! Written by F. Neri Feb 7 1986
4138	!
4139	implicit none
4140	integer n
4141	parameter ( n = 3 )
4142	integer kp,kq,lp,lq,jp,jq,i
4143	double precision m(2n,2n)
4144	double precision qq,pq,qp,pp
4145	!
4146	do kp=2,2*n,2
4147	kq = kp-1
4148	do lp=2,kp-2,2
4149	lq = lp-1
4150	qq = 0.d0
4151	pq = 0.d0
4152	qp = 0.d0
4153	pp = 0.d0
4154	do jp=2,2*n,2
4155	jq = jp-1
4156	qq = qq + m(lq,jq)m(kq,jp) - m(lq,jp)m(kq,jq)
4157	pq = pq + m(lp,jq)m(kq,jp) - m(lp,jp)m(kq,jq)
4158	qp = qp + m(lq,jq)m(kp,jp) - m(lq,jp)m(kp,jq)
4159	pp = pp + m(lp,jq)m(kp,jp) - m(lp,jp)m(kp,jq)
4160	enddo
4161	! write(6,*) qq,pq,qp,pp
4162	do i=1,2*n
4163	m(kq,i) = m(kq,i) - qqm(lp,i) + pqm(lq,i)
4164	m(kp,i) = m(kp,i) - qpm(lp,i) + ppm(lq,i)
4165	enddo
4166	enddo
4167	qp = 0.d0
4168	do jp=2,2*n,2
4169	jq = jp-1
4170	qp = qp + m(kq,jq)m(kp,jp) - m(kq,jp)m(kp,jq)
4171	enddo
4172	! write(6,*) qp
4173	do i=1,2*n
4174	m(kp,i) = m(kp,i)/qp
4175	enddo
4176	!
4177	! Maybe the following is a better idea ( uses sqrt and is slower )
4178	! sign = 1.d0
4179	! if ( qp.lt.0.d0 ) sign = -1.d0
4180	! OR, BETTER:
4181	! if ( qp.le.0.d0 ) then complain
4182	! qp = dabs(qp)
4183	! qp = dsqrt(qp)
4184	! do 600 i=1,2*n
4185	! m(kq,i) = m(kq,i)/qp
4186	! m(kp,i) = sign*m(kp,i)/qp
4187	! 600 continue
4188	enddo
4189	return
4190	end
4191
4192
4193	logical*1 function averaged(f,a,flag,fave)
4194	implicit none
4195	integer isi,ndim,ndim2,nord,ntt
4196	double precision avepol
4197	! TAKES THE AVERAGE OF A FUNCTION F
4198	! FLAG TRUE A=ONE TURN MAP
4199	! FALSE A=A_SCRIPT
4200	!
4201	parameter (ndim=3)
4202	parameter (ndim2=6)
4203	parameter (ntt=40)
4204	integer idpr
4205	common /printing/ idpr
4206	integer nd,nd2,no,nv
4207	common /ii/no,nv,nd,nd2
4208	integer f,fave,a(*)
4209	integer cosi,sine
4210	logical flag
4211	external avepol
4212
4213	integer a1(ndim2),a2(ndim2),xy(ndim2),hf(ndim2),ftf(ndim2)
4214
4215	logical*1 mapnormf
4216
4217	if(.not.flag) then
4218	call etall(cosi,1)
4219	call etall(sine,1)
4220	call trx(f,f,a)
4221	call ctor(f,cosi,sine)
4222	call dacfu(cosi,avepol,fave)
4223	call dadal(cosi,1)
4224	call dadal(sine,1)
4225	else
4226
4227	call etall(cosi,1)
4228	call etall(sine,1)
4229	call etall(ftf,nd2)
4230	call etall(hf,nd2)
4231	call etall(a2,nd2)
4232	call etall(a1,nd2)
4233	call etall(xy,nd2)
4234
4235	isi=0
4236	nord=1
4237	averaged = mapnormf(a,ftf,a2,a1,xy,hf,nord,isi)
4238	nord=no
4239	call etcct(a1,a2,xy)
4240	call facflod(ftf,xy,a1,2,nord,1.d0,-1)
4241	call trx(f,f,a1)
4242	call ctor(f,cosi,sine)
4243	call dacfu(cosi,avepol,fave)
4244
4245	call dadal(cosi,1)
4246	call dadal(sine,1)
4247	call dadal(ftf,nd2)
4248	call dadal(hf,nd2)
4249	call dadal(a2,nd2)
4250	call dadal(a1,nd2)
4251	call dadal(xy,nd2)
4252
4253	endif
4254
4255	return
4256	end
4257	double precision function avepol(j)
4258	implicit none
4259	integer i,ndim
4260	parameter (ndim=3)
4261	! PARAMETER (NTT=40)
4262	! INTEGER J(NTT)
4263	integer j(*)
4264	integer nd,nd2,no,nv
4265	common /ii/no,nv,nd,nd2
4266	integer ndc,ndc2,ndpt,ndt
4267	common /coast/ndc,ndc2,ndt,ndpt
4268	avepol=1.d0
4269	do i=1,(nd-ndc)
4270	if(j(2i).ne.j(2i-1)) then
4271	avepol=0.d0
4272	return
4273	endif
4274	enddo
4275
4276	return
4277	end
4278
4279
4280	logical function couplean(map1,tune,map2,oneturn)
4281	implicit none
4282	integer i,ndim,ndim2,no1,nord,ntt
4283	double precision crazy,tpi
4284	! map1 ascript a1 not there
4285	! tune 2 or 3 tunes
4286
4287	! map2 ascript with a couple parameter in nv
4288	! oneturn map created with tunes and map2
4289
4290	parameter (ndim=3)
4291	parameter (ndim2=6)
4292	parameter (ntt=40)
4293	integer ndc,ndc2,ndpt,ndt
4294	common /coast/ndc,ndc2,ndt,ndpt
4295	integer nd,nd2,no,nv
4296	common /ii/no,nv,nd,nd2
4297	integer map1(),oneturn(),map2(*),ftf,hf
4298	integer xy(ndim2),m1(ndim2),m2(ndim2),a2(ndim2),a1(ndim2)
4299	integer cs,h
4300
4301	double precision killnonl,planar,psq(ndim),radsq(ndim)
4302	double precision tune(ndim)
4303	external killnonl,planar
4304
4305	logical*1 mapnorm
4306
4307	call etall(ftf,1)
4308	call etall(hf,1)
4309	call etall(a1,nd2)
4310	call etall(a2,nd2)
4311	call etall(m1,nd2)
4312	call etall(m2,nd2)
4313	call etall(xy,nd2)
4314	call etall(cs,1)
4315	call etall(h,1)
4316
4317	! map1 is an a-script, the last nv entry should be empty
4318	! this a-script should around the fixed point to all orders
4319	! one order is lost because I use PB-field
4320
4321	tpi=datan(1.d0)*8.d0
4322	do i=1,nd2
4323	call dacfu(map1(i),killnonl,m1(i))
4324	enddo
4325
4326	call etini(xy)
4327	call daclr(cs)
4328
4329	do i=1,nd-ndc
4330	call dasqr(xy(2i-1),a2(2i-1))
4331	call dasqr(xy(2i),a2(2i))
4332	call daadd(a2(2i-1),a2(2i),ftf)
4333	crazy=-tune(i)*tpi/2.d0
4334	call dacmu(ftf,crazy,ftf)
4335	call daadd(ftf,cs,cs)
4336	enddo
4337
4338	call etinv(m1,m2)
4339	call trx(cs,h,m2)
4340
4341	call dacfu(h,planar,cs)
4342	call dasub(h,cs,h)
4343	call davar(a2(1),0.d0,nv)
4344
4345	call damul(a2(1),h,h)
4346	call daadd(cs,h,h)
4347	call expnd2(h,xy,xy,1.d-9,1000)
4348
4349	call dacopd(xy,oneturn)
4350
4351	nord=1
4352	couplean = mapnorm(xy,ftf,a2,a1,m2,hf,nord)
4353
4354	call gettura(psq,radsq)
4355	write(6,*) (psq(i),i=1,nd)
4356
4357	call etini(xy)
4358	no1=no
4359	call fexpo(ftf,xy,xy,3,no1,1.d0,-1)
4360	call etcct(a2,xy,map2)
4361	call etcct(a1,map2,map2)
4362
4363	call dadal(ftf,1)
4364	call dadal(hf,1)
4365	call dadal(a1,nd2)
4366	call dadal(a2,nd2)
4367	call dadal(m1,nd2)
4368	call dadal(m2,nd2)
4369	call dadal(xy,nd2)
4370	call dadal(cs,1)
4371	call dadal(h,1)
4372
4373	return
4374	end
4375	double precision function planar(j)
4376	implicit none
4377	integer i,ndim
4378	parameter (ndim=3)
4379	! PARAMETER (NTT=40)
4380	! INTEGER J(NTT)
4381	integer j(*)
4382	integer nd,nd2,no,nv
4383	common /ii/no,nv,nd,nd2
4384	integer ndc,ndc2,ndpt,ndt
4385	common /coast/ndc,ndc2,ndt,ndpt
4386	planar=0.d0
4387	do i=1,(nd-ndc)
4388	if(j(2i).eq.j(2i-1)) then
4389	planar=1.d0
4390	return
4391	endif
4392	if(j(2*i).eq.2) then
4393	planar=1.d0
4394	return
4395	endif
4396	if(j(2*i-1).eq.2) then
4397	planar=1.d0
4398	return
4399	endif
4400	enddo
4401
4402	return
4403	end
4404	double precision function killnonl(j)
4405	implicit none
4406	integer i,ic,ndim
4407	parameter (ndim=3)
4408	! PARAMETER (NTT=40)
4409	! INTEGER J(NTT)
4410	integer j(*)
4411	integer nd,nd2,no,nv
4412	common /ii/no,nv,nd,nd2
4413	integer ndc,ndc2,ndpt,ndt
4414	common /coast/ndc,ndc2,ndt,ndpt
4415
4416	killnonl=1.d0
4417
4418	ic=0
4419	do i=1,nd2-ndc2
4420	ic=ic+j(i)
4421	enddo
4422	if(ic.gt.1) killnonl=0.d0
4423	if(j(nv).ne.0) killnonl=0.d0
4424
4425	return
4426	end
4427	subroutine fexpo1(h,x,w,nrmin,nrmax,sca,ifac)
4428	implicit none
4429	integer ifac,ndim,ndim2,nrma,nrmax,nrmi,nrmin,ntt
4430	double precision sca
4431	parameter (ndim=3)
4432	parameter (ndim2=6)
4433	parameter (ntt=40)
4434	integer nd,nd2,no,nv
4435	common /ii/no,nv,nd,nd2
4436	integer x,w,h
4437
4438	integer v(ndim2)
4439
4440	nrmi=nrmin-1
4441	nrma=nrmax-1
4442	call etall(v,nd2)
4443	call difd(h,v,-1.d0)
4444	call facflo(v,x,w,nrmi,nrma,sca,ifac)
4445	call dadal(v,nd2)
4446
4447	return
4448	end
4449	subroutine etcctpar(x,ix,xj,z)
4450	implicit none
4451	integer i,ie,ix,ndim,ndim2,ntt
4452	double precision xj
4453	parameter (ndim=3)
4454	parameter (ndim2=6)
4455	parameter (ntt=40)
4456	dimension xj(*)
4457	dimension ie(ntt)
4458	integer nd,nd2,no,nv
4459	common /ii/no,nv,nd,nd2
4460	integer x(),z()
4461
4462	call etallnom(ie,nv,'IE ')
4463	do i=1,nd2
4464	call davar(ie(i),0.d0,i)
4465	enddo
4466	do i=nd2+1,nv
4467	call dacon(ie(i),xj(i-nd2))
4468	enddo
4469
4470	call dacct(x,ix,ie,nv,z,ix)
4471
4472	call dadal(ie,nv)
4473	return
4474	end

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

Context Navigation

source: TRACY3/trunk/tracy/TPSA/LieLib.f @ 3

Download in other formats: