source: PSPA/madxPSPA/libs/ptc/src/d_lielib.f90 @ 430

Last change on this file since 430 was 430, checked in by touze, 11 years ago

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