subroutine lieinit(no1,nv1,nd1,ndpt1,iref1,nis)
implicit none
integer i,iref1,nd1,ndc1,ndim,ndpt1,nis,no1,nv1
double precision ang,ra,st
!! Lieinit initializes AD Package and Lielib
parameter (ndim=3)
dimension st(ndim),ang(ndim),ra(ndim)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer iref
common /resfile/iref
integer itu
common /tunedef/itu
integer lienot
common /dano/lienot
integer idpr
common /printing/ idpr
double precision xintex
common /integratedex/ xintex(0:20)
integer idao,is,iscrri
double precision rs
common/dascr/is(100),rs(100),iscrri(100),idao
integer nplane
double precision epsplane,xplane
common /choice/ xplane(ndim),epsplane,nplane(ndim)
!+CA DASCR
call daexter
do i=1,ndim
nplane(i)=2*i-1
ang(i)=0.d0
ra(i)=0.d0
st(i)=1.d0
enddo
no=no1
nv=nv1
nd=nd1
nd2=2*nd1
do i=1,100
is(i)=0
enddo
write(*, *) "lieinit: calling daini"
call daini(no,nv,0)
if(nis.gt.0)call etallnom(is,nis,'$$IS ')
if(ndpt1.eq.0) then
ndpt=0
ndt=0
ndc1=0
else
ndpt=ndpt1
ndc1=1
if(ndpt.eq.nd2) then
ndt=nd2-1
else
ndt=nd2
if(ndpt.ne.nd2-1) then
write(6,*) ' LETHAL ERROR IN LIEINIT'
stop
endif
endif
endif
ndc=ndc1
ndc2=2*ndc1
iref=0
call initpert(st,ang,ra)
iref=iref1
if(iref1.eq.0) then
itu=0
else
itu=1
endif
if(iref1.eq.0) iref=-1
if(idpr.eq.1)write(6,*) ' NO = ',no,' IN DA-CALCULATIONS '
do i=0,20
xintex(i)=0.d0
enddo
xintex( 0)= 1.000000000000000
xintex( 1)= 5.000000000000000e-001
xintex( 2)= 8.333333333333334e-002
xintex( 3)= 0.000000000000000e+000
xintex( 4)= -1.388888888888898e-003
xintex( 5)= 0.000000000000000e+000
xintex( 6)= 3.306878306878064e-005
xintex( 7)= 0.d0
xintex( 8)= -8.267195767165669e-007
xintex( 9)= 0.d0
xintex( 10)= 4.592886537931051e-008
return
end
subroutine flowpara(ifl,jtu)
implicit none
integer iflow,jtune
common /vecflow/ iflow,jtune
integer ifl,jtu
iflow=ifl
jtune=jtu
return
end
subroutine pertpeek(st,ang,ra)
implicit none
integer i,ndim,ndim2,nreso,ntt
double precision ang,ra,st
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
parameter (nreso=20)
dimension st(ndim),ang(ndim),ra(ndim)
double precision angle,dsta,rad,sta
common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer mx,nres
common /reson/mx(ndim,nreso),nres
integer iref
common /resfile/iref
do i=1,nd
st(i)=sta(i)
ang(i)=angle(i)
ra(i)=rad(i)
enddo
return
end
subroutine inputres(mx1,nres1)
implicit none
integer i,j,ndim,ndim2,nreso,ntt
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
parameter (nreso=20)
integer mx1(ndim,nreso),nres1
integer mx,nres
common /reson/mx(ndim,nreso),nres
nres=nres1
do i=1,nreso
do j=1,ndim
mx(j,i)=0
enddo
enddo
do i=1,nres
do j=1,ndim
mx(j,i)=mx1(j,i)
enddo
enddo
return
end
subroutine respoke(mres,nre,ire)
implicit none
integer i,ire,j,ndim,ndim2,nre,nreso,ntt
double precision ang,ra,st
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
parameter (nreso=20)
integer mres(ndim,nreso)
double precision angle,dsta,rad,sta
common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer mx,nres
common /reson/mx(ndim,nreso),nres
integer iref
common /resfile/iref
dimension ang(ndim),ra(ndim),st(ndim)
iref=ire
nres=nre
do j=1,nreso
do i=1,nd
mx(i,j)=mres(i,j)
enddo
enddo
call initpert(st,ang,ra)
return
end
subroutine liepeek(iia,icoast)
implicit none
integer ndim,ndim2,nreso,ntt
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
parameter (nreso=20)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer iia(*),icoast(*)
iia(1)=no
iia(2)=nv
iia(3)=nd
iia(4)=nd2
icoast(1)=ndc
icoast(2)=ndc2
icoast(3)=ndt
icoast(4)=ndpt
return
end
subroutine lienot(not)
implicit none
integer no,not
call danot(not)
no=not
return
end
subroutine etallnom(x,n,nom)
implicit none
integer i,n,nd2
! CREATES A AD-VARIABLE WHICH CAN BE DESTROYED BY DADAL
! allocates vector of n polynomials and give it the name NOM=A10
integer x(*),i1(4),i2(4)
character*10 nom
do i=1,iabs(n)
x(i)=0
enddo
call daallno(x,iabs(n),nom)
if(n.lt.0) then
call liepeek(i1,i2)
nd2=i1(4)
do i=nd2+1,-n
call davar(x(i),0.d0,i)
enddo
endif
return
end
subroutine etall(x,n)
implicit none
integer i,n,nd2
! allocates vector of n polynomials
integer x(*),i1(4),i2(4)
do i=1,iabs(n)
x(i)=0
enddo
call daallno(x,iabs(n),'ETALL ')
if(n.lt.0) then
call liepeek(i1,i2)
nd2=i1(4)
do i=nd2+1,-n
call davar(x(i),0.d0,i)
enddo
endif
return
end
subroutine etall1(x)
implicit none
integer x
call daallno(x,1,'ETALL ')
return
end
subroutine dadal1(x)
implicit none
integer x
call dadal(x,1)
return
end
subroutine etppulnv(x,xi,xff)
implicit none
integer i,ndim,ndim2,ntt
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*)
double precision xi(*),xff(*),xf(ntt),xii(ntt)
do i=1,nv
xii(i)=xi(i)
enddo
do i=nv+1,ntt
xii(i)=0.d0
enddo
call ppush(x,nv,xii,xf)
do i=1,nv
xff(i)=xf(i)
enddo
return
end
subroutine etmtree(y,x)
implicit none
integer i,ie,iv,ndim,ndim2,nt,ntt
! ROUTINES USING THE MAP IN AD-FORM
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
dimension ie(ntt),iv(ntt)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*),y(*)
nt=nv-nd2
if(nt.gt.0) then
call etallnom(ie,nt,'IE ')
do i=nd2+1,nv
call davar(ie(i-nd2),0.d0,i)
enddo
do i=nd2+1,nv
iv(i)=ie(i-nd2)
enddo
endif
do i=1,nd2
iv(i)=y(i)
enddo
call mtree(iv,nv,x,nv)
if(nt.gt.0) then
call dadal(ie,nt)
endif
return
end
subroutine etppush(x,xi)
implicit none
integer i,ndim,ndim2,ntt
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*)
double precision xi(*),xf(ntt),xii(ntt)
do i=1,nd2
xii(i)=xi(i)
enddo
call ppush(x,nv,xii,xf)
do i=1,nd2
xi(i)=xf(i)
enddo
return
end
subroutine etppush2(x,xi,xff)
implicit none
integer i,ndim,ndim2,ntt
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*)
double precision xi(*),xff(*),xf(ntt),xii(ntt)
do i=1,nd2
xii(i)=xi(i)
enddo
call ppush(x,nv,xii,xf)
do i=1,nd2
xff(i)=xf(i)
enddo
return
end
subroutine ppushlnv(x,xi,xff,nd1)
implicit none
integer i,nd1,ndim,ndim2,ntt
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*)
double precision xi(*),xff(*),xf(ntt),xii(ntt)
do i=1,nd1
xii(i)=xi(i)
enddo
do i=nd1+1,ntt
xii(i)=0.d0
enddo
call ppush(x,nv,xii,xf)
do i=1,nd1
xff(i)=xf(i)
enddo
return
end
subroutine etcct(x,y,z)
implicit none
integer i,ie,iv,ndim,ndim2,nt,ntt
! Z=XoY
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
dimension ie(ntt),iv(ntt)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*),y(*),z(*)
nt=nv-nd2
if(nt.gt.0) then
call etallnom(ie,nt,'IE ')
do i=nd2+1,nv
call davar(ie(i-nd2),0.d0,i)
enddo
do i=nd2+1,nv
iv(i)=ie(i-nd2)
enddo
endif
do i=1,nd2
iv(i)=y(i)
enddo
call dacct(x,nd2,iv,nv,z,nd2)
if(nt.gt.0) then
call dadal(ie,nt)
endif
return
end
subroutine trx(h,rh,y)
implicit none
integer i,ie,iv,ndim,ndim2,nt,ntt
! :RH: = Y :H: Y^-1 = :HoY:
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
dimension ie(ntt),iv(ntt)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h,rh
integer y(*)
!
nt=nv-nd2
if(nt.gt.0) then
call etallnom(ie,nt,'IE ')
do i=nd2+1,nv
call davar(ie(i-nd2),0.d0,i)
enddo
do i=nd2+1,nv
iv(i)=ie(i-nd2)
enddo
endif
do i=1,nd2
iv(i)=y(i)
enddo
call dacct(h,1,iv,nv,rh,1)
if(nt.gt.0) then
call dadal(ie,nt)
endif
return
end
subroutine trxflo(h,rh,y)
implicit none
integer j,k,ndim,ndim2,ntt
! *RH* = Y *H* Y^-1 CHANGE OF A VECTOR FLOW OPERATOR
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h(*),rh(*),y(*)
integer yi(ndim2),ht(ndim2),b1,b2
!
!
call etallnom(yi,nd2 ,'YI ')
call etallnom(ht,nd2 ,'HT ')
call etallnom(b1,1,'B1 ')
call etallnom(b2,1,'B2 ')
call etinv(y,yi)
!----- HT= H o Y
call etcct(h,y,ht)
!----
call daclrd(rh)
do j=1,nd2
do k=1,nd2
call dader(k,yi(j),b1)
call trx(b1,b2,y)
call damul(b2,ht(k),b1)
call daadd(b1,rh(j),b2)
call dacop(b2,rh(j))
enddo
enddo
call dadal(b2,1)
call dadal(b1,1)
call dadal(ht,nd2)
call dadal(yi,nd2)
return
end
subroutine simil(a,x,ai,y)
implicit none
integer ndim,ndim2,ntt
! Y= AoXoAI
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*),y(*),a(*),ai(*)
integer w(ndim2),v(ndim2)
!
call etallnom(w,nd2 ,'W ')
call etallnom(v,nd2 ,'V ')
call etcct(a,x,w)
call etcct(w,ai,v)
call dacopd(v,y)
call dadal(v,nd2)
call dadal(w,nd2)
return
end
subroutine etini(x)
implicit none
integer i,ndim,ndim2,ntt
! X=IDENTITY
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*)
!*DAEXT(NO,NV) X(NDIM2)
do i=1,nd2
call davar(x(i),0.d0,i)
enddo
return
end
subroutine etinv(x,y)
implicit none
integer i,ie1,ie2,iv1,iv2,ndim,ndim2,nt,ntt
! Y=X^-1
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
dimension ie1(ntt),ie2(ntt),iv1(ntt),iv2(ntt)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*),y(*)
nt=nv-nd2
if(nt.gt.0) then
do i=1,nt
ie1(i)=0
ie2(i)=0
enddo
call etallnom(ie1,nt,'IE1 ')
call etallnom(ie2,nt,'IE2 ')
do i=nd2+1,nv
call davar(ie1(i-nd2),0.d0,i)
enddo
do i=nd2+1,nv
iv1(i)=ie1(i-nd2)
iv2(i)=ie2(i-nd2)
enddo
endif
do i=1,nd2
iv1(i)=x(i)
iv2(i)=y(i)
enddo
call dainv(iv1,nv,iv2,nv)
if(nt.gt.0) then
call dadal(ie2,nt)
call dadal(ie1,nt)
endif
return
end
subroutine etpin(x,y,jj)
implicit none
integer i,ie1,ie2,iv1,iv2,jj,ndim,ndim2,nt,ntt
! Y=PARTIAL INVERSION OF X SEE BERZ'S PACKAGE
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
dimension ie1(ntt),ie2(ntt),iv1(ntt),iv2(ntt),jj(*)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*),y(*)
nt=nv-nd2
if(nt.gt.0) then
do i=1,nt
ie1(i)=0
ie2(i)=0
enddo
call etallnom(ie1,nt,'IE1 ')
call etallnom(ie2,nt,'IE2 ')
do i=nd2+1,nv
call davar(ie1(i-nd2),0.d0,i)
enddo
do i=nd2+1,nv
iv1(i)=ie1(i-nd2)
iv2(i)=ie2(i-nd2)
enddo
endif
do i=1,nd2
iv1(i)=x(i)
iv2(i)=y(i)
enddo
call dapin(iv1,nv,iv2,nv,jj)
if(nt.gt.0) then
call dadal(ie2,nt)
call dadal(ie1,nt)
endif
return
end
subroutine dapek0(v,x,jj)
implicit none
integer i,jj,ndim2,ntt
double precision x
!- MORE EXTENSIONS OF BASIC BERZ'S PACKAGE
parameter (ndim2=6)
parameter (ntt=40)
integer v(*),jd(ntt)
dimension x(*)
do i=1,ntt
jd(i)=0
enddo
do i=1,jj
call dapek(v(i),jd,x(i))
enddo
return
end
subroutine dapok0(v,x,jj)
implicit none
integer i,jj,ndim2,ntt
double precision x
parameter (ndim2=6)
parameter (ntt=40)
integer v(*),jd(ntt)
dimension x(*)
do i=1,ntt
jd(i)=0
enddo
do i=1,jj
call dapok(v(i),jd,x(i))
enddo
return
end
subroutine dapokzer(v,jj)
implicit none
integer i,jj,ndim2,ntt
parameter (ndim2=6)
parameter (ntt=40)
integer v(*),jd(ntt)
do i=1,ntt
jd(i)=0
enddo
do i=1,jj
call dapok(v(i),jd,0.d0)
enddo
return
end
subroutine davar0(v,x,jj)
implicit none
integer i,jj,ndim2,ntt
double precision x
parameter (ndim2=6)
parameter (ntt=40)
integer v(*)
dimension x(*)
do i=1,jj
call davar(v(i),x(i),i)
enddo
return
end
subroutine comcfu(b,f1,f2,c)
implicit none
double precision f1,f2
external f1,f2
! Complex dacfu
integer b(*),c(*),t(4)
call etall(t,4)
call dacfu(b(1),f1,t(1))
call dacfu(b(1),f2,t(2))
call dacfu(b(2),f1,t(3))
call dacfu(b(2),f2,t(4))
call dasub(t(1),t(4),c(1))
call daadd(t(2),t(3),c(2))
call dadal(t,4)
return
end
subroutine take(h,m,ht)
implicit none
integer i,m,ndim,ntt
double precision r
! HT= H_M (TAKES M^th DEGREE PIECE ALL VARIABLES INCLUDED)
parameter (ndim=3)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h,ht,j(ntt)
integer b1,b2,b3
!
!
call etallnom(b1,1,'B1 ')
call etallnom(b2,1,'B2 ')
call etallnom(b3,1,'B3 ')
if(no.ge.2) then
if(m.eq.0) then
do i=1,ntt
j(i)=0
enddo
call dapek(h,j,r)
call dacon(ht,r)
else
call danot(m)
call dacop(h,b1)
call danot(m-1)
call dacop(b1,b2)
call danot(no)
call dasub(b1,b2,b3)
call dacop(b3,ht)
endif
else
do i=1,ntt
j(i)=0
enddo
if(m.eq.0) then
call dapek(h,j,r)
call dacon(ht,r)
elseif(m.eq.1) then
do i=1,nv
j(i)=1
call dapek(h,j,r)
call dapok(b3,j,r)
j(i)=0
enddo
call dacop(b3,ht)
else
call daclr(ht)
endif
endif
call dadal(b3,1)
call dadal(b2,1)
call dadal(b1,1)
return
end
subroutine taked(h,m,ht)
implicit none
integer i,m,ndim2,ntt
! \VEC{HT}= \VEC{H_M} (TAKES M^th DEGREE PIECE ALL VARIABLES INCLUDED)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h(*),ht(*),j(ntt)
integer b1,b2,x(ndim2)
!
call etallnom(b1,1,'B1 ')
call etallnom(b2,1,'B2 ')
call etallnom(x,nd2 ,'X ')
do i=1,ntt
j(i)=0
enddo
do i=1,nd2
call take(h(i),m,ht(i))
enddo
call dadal(x,nd2)
call dadal(b2,1)
call dadal(b1,1)
return
end
subroutine daclrd(h)
implicit none
integer i,ndim2,ntt
! clear a map : a vector of nd2 polynomials
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h(*)
do i=1,nd2
call daclr(h(i))
enddo
return
end
subroutine dacopd(h,ht)
implicit none
integer i,ndim2,ntt
! H goes into HT (nd2 array)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h(*),ht(*)
do i=1,nd2
call dacop(h(i),ht(i))
enddo
return
end
subroutine dacmud(h,sca,ht)
implicit none
integer i,ndim2,ntt
double precision sca
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h(*),ht(*)
do i=1,nd2
call dacmu(h(i),sca,ht(i))
enddo
return
end
subroutine dalind(h,rh,ht,rt,hr)
implicit none
integer i,ndim2
double precision rh,rt
parameter (ndim2=6)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h(*),ht(*),hr(*)
integer b(ndim2)
!
call etallnom(b,nd2 ,'B ')
do i=1,nd2
call dalin(h(i),rh,ht(i),rt,b(i))
enddo
call dacopd(b,hr)
call dadal(b,nd2)
return
end
subroutine daread(h,nd1,mfile,xipo)
implicit none
integer i,mfile,nd1,ndim2,ntt
double precision rx,xipo
! read a map
parameter (ndim2=6)
parameter (ntt=40)
integer h(*),j(ntt)
do i=1,ntt
j(i)=0
enddo
do i=1,nd1
call darea(h(i),mfile)
call dapek(h(i),j,rx)
rx=rx*xipo
call dapok(h(i),j,rx)
enddo
return
end
subroutine daprid(h,n1,n2,mfile)
implicit none
integer i,mfile,n1,n2,ndim2,ntt
! print a map
parameter (ndim2=6)
parameter (ntt=40)
integer h(*)
if(mfile.le.0) return
do i=n1,n2
call dapri(h(i),mfile)
enddo
return
end
subroutine prresflo(h,eps,mfile)
implicit none
integer i,mfile,ndim2,ntt
double precision deps,eps,filtres
! print a map in resonance basis for human consumption (useless)
parameter (ndim2=6)
parameter (ntt=40)
integer h(*) ,b(ndim2) ,c(ndim2)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ifilt
common /filtr/ ifilt
external filtres
call etall(b,nd2)
call etall(c,nd2)
call dacopd(h,c)
do i=1,nd2
ifilt=(-1)**i
call dacfu(c(i),filtres,h(i))
enddo
deps=-1.d0
call daeps(deps)
call daeps(eps)
call dacopd(c,h)
call daeps(deps)
call dadal(c,nd2)
call dadal(b,nd2)
return
end
double precision function filtres(j)
implicit none
integer i,ic,ndim
parameter (ndim=3)
! PARAMETER (NTT=40)
! INTEGER J(NTT)
integer j(*)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer ifilt
common /filtr/ ifilt
filtres=1.d0
ic=0
do i=1,(nd2-ndc2)
ic=ic+j(i)*(-1)**(i+1)
enddo
ic=ic+ifilt
if(ic.lt.0) filtres=0.d0
if(ic.eq.0.and.ifilt.eq.1) then
filtres=0.0d0
endif
return
end
subroutine daflo(h,x,y)
implicit none
integer i,ndim,ndim2,ntt
! LIE EXPONENT ROUTINES WITH FLOW OPERATORS
! \VEC{H}.GRAD X =Y
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h(*),x,y
integer b1,b2,b3
!
call etallnom(b1,1,'B1 ')
call etallnom(b2,1,'B2 ')
call etallnom(b3,1,'B3 ')
call daclr(b1)
call daclr(b2)
do i=1,nd2
call dader(i,x,b2)
call damul(b2,h(i),b3)
call daadd(b3,b1,b2)
call dacop(b2,b1)
enddo
call dacop(b1,y)
call dadal(b3,1)
call dadal(b2,1)
call dadal(b1,1)
return
end
subroutine daflod(h,x,y)
implicit none
integer i,ndim,ndim2,ntt
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h(*),x(*),y(*)
integer b1(ndim2),b2(ndim2)
!
call etall(b1,nd2)
call etall(b2,nd2)
call dacopd(h,b1)
call dacopd(x,b2)
do i=1,nd2
call daflo(b1,b2(i),y(i))
enddo
call dadal(b1,nd2)
call dadal(b2,nd2)
return
end
subroutine intd(v,h,sca)
implicit none
integer i,ndim,ndim2,ntt
double precision dlie,sca
! IF SCA=-1.D0
! \VEC{V}.GRAD = J GRAD H . GRAD = :H:
! IF SCA=1.D0
! \VEC{V}.GRAD = GRAD H . GRAD
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
external dlie
integer v(*),h
integer b1,b2,b3,b4,x(ndim2)
!
!
call etallnom(b1,1,'B1 ')
call etallnom(b2,1,'B2 ')
call etallnom(b3,1,'B3 ')
call etallnom(b4,1,'B4 ')
call etallnom(x,nd2 ,'X ')
call daclr(b4)
call daclr(h)
call etini(x)
do i=1,nd
call dacfu(v(2*i-1),dlie,b3)
call dacfu(v(2*i),dlie,b1)
call damul(b1,x(2*i-1),b2)
call damul(b3,x(2*i),b1)
call dalin(b2,1.d0,b1,sca,b3)
call daadd(b3,b4,b2)
call dacop(b2,b4)
enddo
call dacop(b4,h)
call dadal(x,nd2)
call dadal(b4,1)
call dadal(b3,1)
call dadal(b2,1)
call dadal(b1,1)
return
end
subroutine difd(h1,v,sca)
implicit none
integer i,ndim,ndim2,ntt
double precision sca
! INVERSE OF INTD ROUTINE
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer v(*),h1
integer b1,h
call etall(b1,1)
call etall(h,1)
call dacop(h1,h)
do i=1,nd
call dader(2*i-1,h,v(2*i))
call dader(2*i,h,b1)
call dacmu(b1,sca,v(2*i-1))
enddo
call dadal(h,1)
call dadal(b1,1)
return
end
subroutine expflo(h,x,y,eps,nrmax)
implicit none
integer i,ndim,ndim2,nrmax,ntt
double precision coe,eps,r,rbefore
! DOES EXP( \VEC{H} ) X = Y
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer idpr
common /printing/ idpr
integer h(*),x,y
integer b1,b2,b3,b4
logical more
!
!
call etallnom(b1,1,'B1 ')
call etallnom(b2,1,'B2 ')
call etallnom(b3,1,'B3 ')
call etallnom(b4,1,'B4 ')
call dacop(x,b4)
call dacop(x,b1)
more=.true.
rbefore=1.d30
do i=1,nrmax
coe=1.d0/dble(i)
call dacmu(b1,coe,b2)
call daflo(h,b2,b1)
call daadd(b4,b1,b3)
call daabs(b1,r)
if(more) then
if(r.gt.eps) then
rbefore=r
goto 100
else
rbefore=r
more=.false.
endif
else
if(r.ge.rbefore) then
call dacop(b3,y)
call dadal(b4,1)
call dadal(b3,1)
call dadal(b2,1)
call dadal(b1,1)
return
endif
rbefore=r
endif
100 continue
call dacop(b3,b4)
enddo
if(idpr.ge.0) then
write(6,*) ' NORM ',eps,' NEVER REACHED IN EXPFLO '
write(6,*) 'NEW IDPR '
read(5,*) idpr
endif
call dacop(b3,y)
call dadal(b4,1)
call dadal(b3,1)
call dadal(b2,1)
call dadal(b1,1)
return
end
subroutine expflod(h,x,w,eps,nrmax)
implicit none
integer j,ndim,ndim2,nrmax,ntt
double precision eps
! DOES EXP( \VEC{H} ) \VEC{X} = \VEC{Y}
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*),w(*),h(*)
integer b0,v(ndim2)
!
!
call etallnom(b0,1,'B0 ')
call etallnom(v,nd2 ,'V ')
call dacopd(x,v)
do j=1,nd2
call expflo(h,v(j),b0,eps,nrmax)
call dacop(b0,v(j))
enddo
call dacopd(v,w)
call dadal(v,nd2)
call dadal(b0,1)
return
end
subroutine facflo(h,x,w,nrmin,nrmax,sca,ifac)
implicit none
integer i,ifac,ndim,ndim2,nmax,nrmax,nrmin,ntt
double precision eps,sca
! IFAC=1
! DOES EXP(SCA \VEC{H}_MRMIN ) ... EXP(SCA \VEC{H}_NRMAX ) X= Y
! IFAC=-1
! DOES EXP(SCA \VEC{H}_NRMAX ) ... EXP(SCA \VEC{H}_MRMIN ) X= Y
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x,w,h(*)
integer bm(ndim2),b0(ndim2),v
!
call etallnom(bm,nd2 ,'BM ')
call etallnom(b0,nd2 ,'B0 ')
call etallnom(v,1 ,'V ')
call dacop(x,v)
eps=-1.d0
call daeps(eps)
nmax=100
!
! IFAC =1 ---> V = EXP(:SCA*H(NRMAX):)...EXP(:SCA*H(NRMIN):)X
if(ifac.eq.1) then
do i=nrmax,nrmin,-1
call taked(h,i,b0)
call dacmud(b0,sca,bm)
call expflo(bm,v,b0(1),eps,nmax)
call dacop(b0(1),v)
enddo
else
! IFAC =-1 ---> V = EXP(:SCA*H(NRMIN):)...EXP(:SCA*H(NRMAX):)X
do i=nrmin,nrmax
call taked(h,i,b0)
call dacmud(b0,sca,bm)
call expflo(bm,v,b0(1),eps,nmax)
call dacop(b0(1),v)
enddo
endif
call dacop(v,w)
call dadal(v,1)
call dadal(b0,nd2)
call dadal(bm,nd2)
return
end
subroutine facflod(h,x,w,nrmin,nrmax,sca,ifac)
implicit none
integer i,ifac,ndim,ndim2,nrmax,nrmin,ntt
double precision sca
! IFAC=1
! DOES EXP(SCA \VEC{H}_MRMIN ) ... EXP(SCA \VEC{H}_NRMAX ) \VEC{X}= \VEC{Y}
! IFAC=-1
! DOES EXP(SCA \VEC{H}_NRMAX ) ... EXP(SCA \VEC{H}_MRMIN ) \VEC{X}= \VEC{Y}
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*),w(*),h(*)
do i=1,nd2
call facflo(h,x(i),w(i),nrmin,nrmax,sca,ifac)
enddo
return
end
subroutine fexpo(h,x,w,nrmin,nrmax,sca,ifac)
implicit none
integer ifac,ndim,ndim2,nrma,nrmax,nrmi,nrmin,ntt
double precision sca
! WRAPPED ROUTINES FOR THE OPERATOR \VEC{H}=:H:
! WRAPPING FACFLOD
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*),w(*),h
integer v(ndim2)
nrmi=nrmin-1
nrma=nrmax-1
call etall(v,nd2)
call difd(h,v,-1.d0)
call facflod(v,x,w,nrmi,nrma,sca,ifac)
call dadal(v,nd2)
return
end
subroutine etcom(x,y,h)
implicit none
integer i,j,ndim,ndim2,ntt
! ETCOM TAKES THE BRACKET OF TWO VECTOR FIELDS.
parameter (ndim2=6)
parameter (ndim=3)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h(*),x(*),y(*),t1,t2,t3(ndim2)
call etall(t1,1)
call etall(t2,1)
call etall(t3,nd2)
do j=1,nd2
do i=1,nd2
call dader(i,x(j),t1)
call dader(i,y(j),t2)
call damul(x(i),t2,t2)
call damul(y(i),t1,t1)
call dalin(t2,1.d0,t1,-1.d0,t1)
call daadd(t1,t3(j),t3(j))
enddo
enddo
call dacopd(t3,h)
call dadal(t1,1)
call dadal(t2,1)
call dadal(t3,nd2)
return
end
subroutine etpoi(x,y,h)
implicit none
integer i,ndim,ndim2,ntt
! ETPOI TAKES THE POISSON BRACKET OF TWO FUNCTIONS
parameter (ndim2=6)
parameter (ndim=3)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h,x,y,t1,t2,t3
call etall(t1,1)
call etall(t2,1)
call etall(t3,1)
do i=1,nd
call dader(2*i-1,x,t1)
call dader(2*i,y,t2)
call damul(t1,t2,t1)
call dalin(t1,1.d0,t3,1.d0,t3)
call dader(2*i-1,y,t1)
call dader(2*i,x,t2)
call damul(t1,t2,t1)
call dalin(t1,-1.d0,t3,1.d0,t3)
enddo
call dacop(t3,h)
call dadal(t1,1)
call dadal(t2,1)
call dadal(t3,1)
return
end
subroutine exp1d(h,x,y,eps,non)
implicit none
integer ndim,ndim2,non,ntt
double precision eps
! WRAPPING EXPFLO
parameter (ndim2=6)
parameter (ndim=3)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer idpr
common /printing/ idpr
integer h,x,y
integer v(ndim2)
call etall(v,nd2)
call difd(h,v,-1.d0)
call expflo(v,x,y,eps,non)
call dadal(v,nd2)
return
end
subroutine expnd2(h,x,w,eps,nrmax)
implicit none
integer j,ndim,ndim2,nrmax,ntt
double precision eps
! WRAPPING EXPFLOD USING EXP1D
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*),w(*),h
integer b0,v(ndim2)
!
!
call etallnom(b0,1,'B0 ')
call etallnom(v,nd2 ,'V ')
call dacopd(x,v)
do j=1,nd2
call exp1d(h,v(j),b0,eps,nrmax)
call dacop(b0,v(j))
enddo
call dacopd(v,w)
call dadal(v,nd2)
call dadal(b0,1)
return
end
subroutine flofacg(xy,h,epsone)
implicit none
integer i,k,kk,ndim,ndim2,nrmax,ntt
double precision eps,epsone,r,xn,xnbefore,xnorm,xnorm1,xx
! GENERAL ONE EXPONENT FACTORIZATION
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer idpr
common /printing/ idpr
logical more
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
double precision xintex
common /integratedex/ xintex(0:20)
integer xy(*),x(ndim2),h(*)
integer v(ndim2),w(ndim2),t(ndim2), z(ndim2)
integer jj(ntt)
jj(1)=1
!
call etallnom(v,nd2 ,'V ')
call etallnom(w,nd2 ,'W ')
call etallnom(t,nd2 ,'T ')
call etallnom(x,nd2 ,'Z ')
call etallnom(z,nd2 ,'Z ')
call etini(v)
call daclrd(w)
xnorm1=0.d0
do i=1,nd2
call daabs(xy(i),r)
xnorm1=xnorm1+r
enddo
xnbefore=1.d36
more=.false.
eps=1.e-9
nrmax=1000
xn=10000.d0
do k=1,nrmax
call dacmud(h,-1.d0,t)
call expflod(t,xy,x,eps,nrmax)
call dalind(x,1.d0,v,-1.d0,t)
! write(20,*) "$$$$$$$$$$$$$$",k,"$$$$$$$$$$$$$$$$$$$$"
! call daprid(t,1,1,20)
if(xn.lt.epsone) then
if(idpr.ge.0) write(6,*) "xn quadratic",xn
call daflod(t,t,w)
call dalind(t,1.d0,w,-0.5d0,t)
call dacopd(t,z)
call dacopd(t,w)
! second order in W
call etcom(h,w,x)
call etcom(x,w,x)
! END OF order in W
do kk=1,10
call etcom(h,w,w)
call dalind(z,1.d0,w,xintex(kk),z)
enddo
call dacopd(z,t)
xx=1.d0/12.d0
call dalind(x,xx,h,1.d0,h)
endif
call dalind(t,1.d0,h,1.d0,h)
xnorm=0.d0
do i=1,nd2
call daabs(t(i),r)
xnorm=xnorm+r
enddo
xn=xnorm/xnorm1
if(xn.ge.epsone.and.(idpr.ge.0)) write(6,*)" xn linear ",xn
if(xn.lt.eps.or.more) then
more=.true.
if(xn.ge.xnbefore) goto 1000
xnbefore=xn
endif
enddo
1000 if(idpr.ge.0) write(6,*) " iteration " , k
call dadal(x,nd2)
call dadal(w,nd2)
call dadal(v,nd2)
call dadal(t,nd2)
call dadal(z,nd2)
return
end
subroutine flofac(xy,x,h)
implicit none
integer k,ndim,ndim2,ntt
! GENERAL DRAGT-FINN FACTORIZATION
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer xy(*),x(*),h(*)
integer v(ndim2),w(ndim2)
!
call etallnom(v,nd2 ,'V ')
call etallnom(w,nd2 ,'W ')
call dacopd(xy,x)
call dacopd(x,v)
call daclrd(w)
call danot(1)
call etinv(v,w)
call danot(no)
call etcct(x,w,v)
call danot(1)
call dacopd(xy,x)
call danot(no)
call dacopd(v,w)
call daclrd(h)
do k=2,no
call taked(w,k,v)
call dalind(v,1.d0,h,1.d0,h)
call facflod(h,w,v,k,k,-1.d0,-1)
call dacopd(v,w)
enddo
call dadal(w,nd2)
call dadal(v,nd2)
return
end
subroutine liefact(xy,x,h)
implicit none
integer ndim,ndim2,ntt
! SYMPLECTIC DRAGT-FINN FACTORIZATION WRAPPING FLOFAC
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer xy(*),x(*),h
integer v(ndim2)
call etall(v,nd2)
call flofac(xy,x,v)
call intd(v,h,-1.d0)
!
call dadal(v,nd2)
return
end
logical*1 function mapnorm(x,ft,a2,a1,xy,h,nord)
implicit none
integer isi,ndim,ndim2,nord,ntt
!--NORMALIZATION ROUTINES OF LIELIB
!- WRAPPING MAPNORMF
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*),a1(*),a2(*),ft,xy(*),h,hf(ndim2),ftf(ndim2)
logical*1 mapnormf
call etall(ftf,nd2)
call etall(hf,nd2)
isi=0
mapnorm = mapnormf(x,ftf,a2,a1,xy,hf,nord,isi)
call intd(hf,h,-1.d0)
call intd(ftf,ft,-1.d0)
call dadal(ftf,nd2)
call dadal(hf,nd2)
return
end
subroutine gettura(psq,radsq)
implicit none
integer ik,ndim,ndim2,ntt
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
double precision psq(ndim),radsq(ndim)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
double precision ps,rads
common /tunerad/ ps(ndim),rads(ndim)
do ik=1,nd
psq(ik)=ps(ik)
radsq(ik)=rads(ik)
enddo
return
end
subroutine setidpr(idprint,nplan)
implicit none
integer idprint,ik,ndim,ndim2,nplan
parameter (ndim=3)
parameter (ndim2=6)
dimension nplan(ndim)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer idpr
common /printing/ idpr
integer nplane
double precision epsplane,xplane
common /choice/ xplane(ndim),epsplane,nplane(ndim)
do ik=1,nd
nplane(ik)=nplan(ik)
enddo
idpr=idprint
return
end
subroutine idprset(idprint)
implicit none
integer idprint,ndim,ndim2
parameter (ndim=3)
parameter (ndim2=6)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer idpr
common /printing/ idpr
integer nplane
double precision epsplane,xplane
common /choice/ xplane(ndim),epsplane,nplane(ndim)
idpr=idprint
return
end
logical*1 function mapnormf(x,ft,a2,a1,xy,h,nord,isi)
implicit none
integer ij,isi,ndim,ndim2,nord,ntt
double precision angle,p,rad,st,x2pi,x2pii
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
dimension angle(ndim),st(ndim),p(ndim),rad(ndim)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer x(*),a1(*),a2(*),ft(*),xy(*),h(*)
integer itu
common /tunedef/itu
integer idpr
common /printing/ idpr
integer iflow,jtune
common /vecflow/ iflow,jtune
double precision ps,rads
common /tunerad/ ps(ndim),rads(ndim)
integer a1i(ndim2),a2i(ndim2)
logical*1 midbflo
!
call etallnom(a1i,nd2 ,'A1I ')
call etallnom(a2i,nd2 ,'A2I ')
! frank/etienne
do itu=1,ndim
angle(itu)=0.d0
p(itu)=0.d0
st(itu)=0.d0
rad(itu)=0.d0
ps(itu)=0.d0
rads(itu)=0.d0
enddo
jtune=isi
x2pii=1.d0/datan(1.d0)/8.d0
x2pi=datan(1.d0)*8.d0
call dacopd(x,xy)
! go to fix point in the parameters + pt to order nord>=1
call gofix(xy,a1,a1i,nord)
call simil(a1i,xy,a1,xy)
! linear part
mapnormf = midbflo(xy,a2,a2i,angle,rad,st)
do ij=1,nd-ndc
p(ij)=angle(ij)*(st(ij)*(x2pii-1.d0)+1.d0)
enddo
if(ndc.eq.1) p(nd)=angle(nd)
if(idpr.ge.0) then
write(6,*) 'tune ',(p(ij),ij=1,nd)
write(6,*) 'damping ', (rad(ij),ij=1,nd)
endif
do ij=1,nd ! -ndc Frank
ps(ij)=p(ij)
rads(ij)=rad(ij)
enddo
call initpert(st,angle,rad)
call simil(a2i,xy,a2,xy)
call dacopd(xy,a2i)
! write(6,*) 'Entering orderflo'
call orderflo(h,ft,xy,angle,rad)
do ij=1,nd-ndc
p(ij)=angle(ij)
if(angle(ij).gt.x2pi/2.d0.and.st(ij).gt.0.d0.and.itu.eq.1)then
p(ij)=angle(ij)-x2pi
write(6,*) ij,' TH TUNE MODIFIED IN H2 TO ',p(ij)/x2pi
endif
enddo
call h2pluflo(h,p,rad)
! CALL TAKED(A2I,1,XY)
call taked(a2i,1,a1i)
call etcct(xy,a1i,xy)
call dadal(a2i,nd2)
call dadal(a1i,nd2)
return
end
subroutine gofix(xy,a1,a1i,nord)
implicit none
integer i,ndim,ndim2,nord,ntt
double precision xic
! GETTING TO THE FIXED POINT AND CHANGING TIME APPROPRIATELY IN THE
! COASTING BEAM CASE
!****************************************************************
! X = A1 XY A1I WHERE X IS TO THE FIXED POINT TO ORDER NORD
! for ndpt not zero, works in all cases. (coasting beam: eigenvalue
!1 in Jordan form)
!****************************************************************
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer xy(*),a1(*),a1i(*)
integer x(ndim2),w(ndim2),v(ndim2),rel(ndim2)
!
call etallnom(x,nd2 , 'X ')
call etallnom(w,nd2 , 'W ')
call etallnom(v,nd2 , 'V ')
call etallnom(rel,nd2 ,'REL ')
! COMPUTATION OF A1 AND A1I USING DAINV
call etini(rel)
call danot(nord)
call etini(v)
do i=1,nd2-ndc2
call dacop(xy(i),x(i))
call dalin(x(i),1.d0,rel(i),-1.d0,v(i))
enddo
call etinv(v,w)
call daclrd(x)
if(ndc.eq.1) then
call davar(x(ndpt),0.d0,ndpt)
endif
call etcct(w,x,v)
if(ndc.eq.1) then
call daclr(v(nd2))
call daclr(v(nd2-ndc))
endif
call dalind(rel,1.d0,v,1.d0,a1)
call dalind(rel,1.d0,v,-1.d0,a1i)
if(ndpt.ne.0) then
! CORRECTIONS
call daclrd(w)
call daclrd(v)
call daclrd(x)
do i=1,nd2-ndc2
call dalin(a1(i),1.d0,rel(i),-1.d0,w(i))
enddo
! COMPUTE Deta/Ddelta
call dacopd(w,a1)
do i=1,nd2-ndc2
call dader(ndpt,w(i),w(i))
enddo
! COMPUTE J*Deta/dDELTA
do i=1,nd-ndc
call dacmu(w(2*i),1.d0,v(2*i-1) )
call dacmu(w(2*i-1),-1.d0,v(2*i) )
enddo
xic=(-1)**(ndt)
do i=1,nd2-ndc2
call damul(v(i),rel(i),x(1))
call daadd(x(1),w(ndt),w(ndt))
call dacop(a1(i),w(i))
enddo
call dacmu(w(ndt),xic,w(ndt))
call expflod(w,rel,a1,1.d-7,10000)
! END OF CORRECTIONS
call etinv(a1,a1i)
endif
call danot(no)
call dadal(rel,nd2)
call dadal(v,nd2)
call dadal(w,nd2)
call dadal(x,nd2)
return
end
double precision function transver(j)
implicit none
integer i,ic,ndim
! USED IN A DACFU CALL OF GOFIX
parameter (ndim=3)
! PARAMETER (NTT=40)
! INTEGER J(NTT)
integer j(*)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
transver=1.d0
ic=0
do i=1,nd2-ndc2
ic=ic+j(i)
enddo
if(ic.ne.1) transver=0.d0
return
end
subroutine orderflo(h,ft,x,ang,ra)
implicit none
integer k,ndim,ndim2,ntt
double precision ang,ra
!- NONLINEAR NORMALIZATION PIECE OF MAPNORMF
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer idpr
common /printing/ idpr
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
dimension ang(ndim),ra(ndim)
integer x(*),h(*),ft(*)
integer w(ndim2),v(ndim2),rel(ndim2)
integer roi(ndim2)
integer b1(ndim2),b5(ndim2),b6(ndim2),b9(ndim2)
!
call etallnom(w,nd2 ,'W ')
call etallnom(v,nd2 ,'V ')
call etallnom(rel,nd2 ,'REL ')
call etallnom(roi,nd2 ,'ROI ')
call etallnom(b1,nd2 ,'B1 ')
call etallnom(b5,nd2 ,'B5 ')
call etallnom(b6,nd2 ,'B6 ')
call etallnom(b9,nd2 ,'B9 ')
call rotiflo(roi,ang,ra)
call etini(rel)
call daclrd(h)
call daclrd(ft)
call etcct(x,roi,x)
do k=2,no
! IF K>2 V = H(K)^-1 X(K)
call facflod(h,x,v,2,k-1,-1.d0,-1)
! EXTRACTING K TH DEGREE OF V ----> W
call taked(v,k,w)
! write(16,*) "$$$$$$$$ K $$$$$$$$$$", k
! W = EXP(B5) + ...
call dacopd(w,b5)
! CALL INTD(W,B5,-1.D0)
! B5 ON EXIT IS THE NEW CONTRIBUTION TO H
! B6 IS THE NEW CONTRIBUTION TO FT
call nuanaflo(b5,b6)
call dalind(b5,1.d0,h,1.d0,b1)
call dacopd(b1,h)
! EXP(B9) = EXP( : ROTI B6 :)
call trxflo(b6,b9,roi)
! V = EXP(-B6) REL
call facflod(b6,rel,v,k,k,-1.d0,1)
! W = V o X
call etcct(v,x,w)
if(idpr.ge.0) then
write(6,*) ' ORDERFLO K= ', k
endif
! X = EXP(B9) W
call facflod(b9,w,x,k,k,1.d0,1)
! B6 IS THE NEW CONTRIBUTION TO FT
call dalind(b6,1.d0,ft,1.d0,b1)
call dacopd(b1,ft)
enddo
call dadal(b9,nd2)
call dadal(b6,nd2)
call dadal(b5,nd2)
call dadal(b1,nd2)
call dadal(roi,nd2)
call dadal(rel,nd2)
call dadal(v,nd2)
call dadal(w,nd2)
return
end
subroutine nuanaflo(h,ft)
implicit none
integer i,ndim,ndim2,ntt
double precision dfilt,filt,xgam,xgbm
! RESONANCE DENOMINATOR OPERATOR (1-R^-1)^-1
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer iflow,jtune
common /vecflow/ iflow,jtune
external xgam,xgbm,dfilt,filt
integer h(*),ft(*),br(ndim2),bi(ndim2),c(ndim2),ci(ndim2)
integer t1(2),t2(2)
call etall(br,nd2)
call etall(bi,nd2)
call etall(c,nd2)
call etall(ci,nd2)
call ctorflo(h,br,bi)
! FILTERING RESONANCES AND TUNE SHIFTS
! ASSUMING REALITY I.E. B(2*I-1)=CMPCJG(B(2*I))
do i=1,nd2
iflow=i
call dacfu(br(i),filt,c(i))
call dacfu(bi(i),filt,ci(i))
enddo
call rtocflo(c,ci,h)
do i=1,nd2
iflow=i
call dacfu(br(i),dfilt,br(i))
call dacfu(bi(i),dfilt,bi(i))
enddo
! NOW WE MUST REORDER C AND CI TO SEPARATE THE REAL AND IMAGINARY PART
! THIS IS NOT NECESSARY WITH :H: OPERATORS
do i=1,nd2
t1(1)=br(i)
t1(2)=bi(i)
t2(1)=c(i)
t2(2)=ci(i)
iflow=i
call comcfu(t1,xgam,xgbm,t2)
enddo
call rtocflo(c,ci,ft)
call dadal(br,nd2)
call dadal(bi,nd2)
call dadal(c,nd2)
call dadal(ci,nd2)
return
end
double precision function xgam(j)
implicit none
integer i,ic,ij,ik,ndim,ndim2
double precision ad,ans,as,ex,exh
! XGAM AND XGBM ARE THE EIGENVALUES OF THE OPERATOR NEWANAFLO
parameter (ndim=3)
parameter (ndim2=6)
! PARAMETER (NTT=40)
integer iflow,jtune
common /vecflow/ iflow,jtune
double precision angle,dsta,rad,sta
common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
! INTEGER J(NTT),JJ(NDIM),JP(NDIM)
integer j(*),jj(ndim),jp(ndim)
xgam=0.d0
ad=0.d0
as=0.d0
ic=0
do i=1,nd-ndc
ik=2*i-1
ij=2*i
jp(i)=j(ik)+j(ij)
jj(i)=j(ik)-j(ij)
if(ik.eq.iflow.or.ij.eq.iflow) then
jj(i)=jj(i)+(-1)**iflow
jp(i)=jp(i)-1
endif
ic=ic+iabs(jj(i))
enddo
do i=1,nd-ndc
ad=dsta(i)*dble(jj(i))*angle(i)-dble(jp(i))*rad(i)+ad
as=sta(i)*dble(jj(i))*angle(i)+as
enddo
exh=dexp(ad/2.d0)
ex=exh**2
ans=4.d0*ex*(dsinh(ad/2.d0)**2+dsin(as/2.d0)**2)
xgam=2.d0*(-exh*dsinh(ad/2.d0)+ex*dsin(as/2.d0)**2)/ans
return
end
double precision function xgbm(j)
implicit none
integer i,ic,ij,ik,ndim,ndim2
double precision ad,ans,as,ex,exh
parameter (ndim=3)
parameter (ndim2=6)
! PARAMETER (NTT=40)
integer iflow,jtune
common /vecflow/ iflow,jtune
double precision angle,dsta,rad,sta
common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
! INTEGER J(NTT),JJ(NDIM),JP(NDIM)
integer j(*),jj(ndim),jp(ndim)
xgbm=0.d0
ad=0.d0
as=0.d0
ic=0
do i=1,nd-ndc
ik=2*i-1
ij=2*i
jp(i)=j(ik)+j(ij)
jj(i)=j(ik)-j(ij)
if(ik.eq.iflow.or.ij.eq.iflow) then
jj(i)=jj(i)+(-1)**iflow
jp(i)=jp(i)-1
endif
ic=ic+iabs(jj(i))
enddo
do i=1,nd-ndc
ad=dsta(i)*dble(jj(i))*angle(i)-dble(jp(i))*rad(i)+ad
as=sta(i)*dble(jj(i))*angle(i)+as
enddo
exh=dexp(ad/2.d0)
ex=exh**2
ans=4.d0*ex*(dsinh(ad/2.d0)**2+dsin(as/2.d0)**2)
xgbm=dsin(as)*ex/ans
return
end
double precision function filt(j)
implicit none
integer i,ic,ic1,ic2,ij,ik,ji,ndim,ndim2,nreso
! PROJECTION FUNCTIONS ON THE KERNEL ANMD RANGE OF (1-R^-1)
!- THE KERNEL OF (1-R^-1)
parameter (ndim=3)
parameter (ndim2=6)
! PARAMETER (NTT=40)
parameter (nreso=20)
double precision angle,dsta,rad,sta
common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer mx,nres
common /reson/mx(ndim,nreso),nres
integer iflow,jtune
common /vecflow/ iflow,jtune
! INTEGER J(NTT),JJ(NDIM)
integer j(*),jj(ndim)
filt=1.d0
ic=0
do i=1,nd-ndc
ik=2*i-1
ij=2*i
jj(i)=j(ik)-j(ij)
if(ik.eq.iflow.or.ij.eq.iflow) then
jj(i)=jj(i)+(-1)**iflow
endif
ic=ic+iabs(jj(i))
enddo
if(ic.eq.0.and.jtune.eq.0) return
do i=1,nres
ic1=1
ic2=1
do ji=1,nd-ndc
if(mx(ji,i).ne.jj(ji)) ic1=0
if(mx(ji,i).ne.-jj(ji)) ic2=0
if(ic1.eq.0.and.ic2.eq.0) goto 3
enddo
return
3 continue
enddo
filt=0.d0
return
end
double precision function dfilt(j)
implicit none
integer ndim,ndim2,nreso
double precision fil,filt
!- THE RANGE OF (1-R^-1)^1
!- CALLS FILT AND EXCHANGES 1 INTO 0 AND 0 INTO 1.
parameter (ndim=3)
parameter (ndim2=6)
! PARAMETER (NTT=40)
parameter (nreso=20)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer mx,nres
common /reson/mx(ndim,nreso),nres
external filt
! INTEGER J(NTT)
integer j(*)
fil=filt(j)
if(fil.gt.0.5d0) then
dfilt=0.d0
else
dfilt=1.d0
endif
return
end
subroutine dhdjflo(h,t)
implicit none
integer i,ndim,ndim2,ntt
double precision coe,x2pi
! CONVENIENT TUNE SHIFT FINDED FOR SYMPLECTIC CASE (NU,DL)(H)=T
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer h(*),t(*)
integer b1(ndim2),b2(ndim2),bb1,bb2
!
call etall(b1,nd2)
call etall(b2,nd2)
call etall(bb1,1)
call etall(bb2,1)
x2pi=datan(1.d0)*8.d0
call ctorflo(h,b1,b2)
coe=1.d0/x2pi
do i=1,nd-ndc
call datra(2*i,b2(2*i),bb1)
call dacmu(bb1,coe,t(i+nd))
call dacop(t(i+nd),bb1)
call daclr(bb2)
call rtoc(bb1,bb2,bb1)
call dacop(bb1,t(i))
enddo
if(ndpt.ne.0) then
call dacop(h(ndt),t(nd))
call dacop(b1(ndt),t(nd2))
endif
call dadal(bb2,1)
call dadal(bb1,1)
call dadal(b2,nd2)
call dadal(b1,nd2)
return
end
subroutine dhdj(h,t)
implicit none
integer i,ndim,ndim2,ntt
double precision coe,x2pi
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer h,t(*)
integer b1,b2,bb1,bb2
!
call etallnom(b1,1,'B1 ')
call etallnom(b2,1,'B2 ')
call etallnom(bb1,1,'BB1 ')
call etallnom(bb2,1,'BB2 ')
x2pi=datan(1.d0)*8.d0
call ctor(h,b1,b2)
coe=-2.d0/x2pi
do i=1,nd-ndc
call dader(2*i-1,b1,b2)
call datra(2*i,b2,bb2)
call dacmu(bb2,coe,t(i+nd))
call dacop(t(i+nd),bb2)
call daclr(b2)
call rtoc(bb2,b2,bb1)
call dacop(bb1,t(i))
enddo
if(ndpt.eq.nd2) then
call dader(ndpt,h,t(nd))
call dader(ndpt,b1,t(nd2))
call dacmu(t(nd),-1.d0,t(nd))
call dacmu(t(nd2),-1.d0,t(nd2))
endif
if(ndt.eq.nd2) then
call dader(ndpt,h,t(nd))
call dader(ndpt,b1,t(nd2))
endif
call dadal(bb2,1)
call dadal(bb1,1)
call dadal(b2,1)
call dadal(b1,1)
return
end
subroutine h2pluflo(h,ang,ra)
implicit none
integer i,j,ndim,ndim2,ntt
double precision ang,r1,r2,ra,st
! POKES IN \VEC{H} ANGLES AND DAMPING COEFFFICIENTS
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
double precision angle,dsta,rad,sta
common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
dimension ang(ndim),st(ndim),ra(ndim),j(ntt)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer h(*)
!*DAEXT(NO,NV) H
do i=1,nd
st(i)=2.d0*sta(i)-1.d0
enddo
do i=1,ntt
j(i)=0
enddo
do i=1,nd-ndc
j(2*i-1)=1
r1=-ang(i)
!-----
call dapok(h(2*i),j,r1)
r2=ra(i)
call dapok(h(2*i-1),j,r2)
j(2*i-1)=0
j(2*i)=1
r1=ang(i)*st(i)
call dapok(h(2*i-1),j,r1)
call dapok(h(2*i),j,r2)
j(2*i)=0
enddo
if(ndpt.eq.nd2-1) then
j(ndpt)=1
call dapok(h(ndt),j,ang(nd))
elseif(ndpt.eq.nd2) then
j(ndpt)=1
call dapok(h(ndt),j,-ang(nd))
endif
return
end
subroutine rotflo(ro,ang,ra)
implicit none
integer i,ndim,ndim2,ntt
double precision ang,ch,co,ra,sh,si,sim,xx
! CREATES R AND R^-1 USING THE EXISTING ANGLES AND DAMPING
! COULD BE REPLACED BY A CALL H2PLUFLO FOLLOWED BY EXPFLOD
! CREATES R
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
dimension co(ndim),si(ndim),ang(ndim),ra(ndim)
integer j(ntt)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer ro(*)
call daclrd(ro)
do i=1,nd-ndc
xx=dexp(ra(i))
if(ista(i).eq.0) then
call hyper(ang(i),ch,sh)
co(i)=ch*xx
si(i)=-sh*xx
else
co(i)=dcos(ang(i))*xx
si(i)=dsin(ang(i))*xx
endif
enddo
do i=1,nd-ndc
if(ista(i).eq.0)then
sim=si(i)
else
sim=-si(i)
endif
j(2*i-1)=1
call dapok(ro(2*i-1),j,co(i))
call dapok(ro(2*i),j,sim)
j(2*i-1)=0
j(2*i)=1
call dapok(ro(2*i),j,co(i))
call dapok(ro(2*i-1),j,si(i))
j(2*i)=0
enddo
if(ndc.eq.1) then
j(ndt)=1
call dapok(ro(ndt),j,1.d0)
call dapok(ro(ndpt),j,0.d0)
j(ndt)=0
j(ndpt)=1
call dapok(ro(ndt),j,ang(nd))
call dapok(ro(ndpt),j,1.d0)
j(ndpt)=0
endif
return
end
subroutine rotiflo(roi,ang,ra)
implicit none
integer i,ndim,ndim2,ntt
double precision ang,ch,co,ra,sh,si,sim,simv,xx
! CREATES R^-1
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
dimension co(ndim),si(ndim),ang(ndim),ra(ndim)
integer j(ntt)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer roi(*)
do i=1,10
j(i)=0
enddo
call daclrd(roi)
do i=1,nd-ndc
xx=dexp(-ra(i))
if(ista(i).eq.0) then
call hyper(ang(i),ch,sh)
co(i)=ch*xx
si(i)=-sh*xx
else
co(i)=dcos(ang(i))*xx
si(i)=dsin(ang(i))*xx
endif
enddo
do i=1,nd-ndc
if(ista(i).eq.0)then
sim=si(i)
else
sim=-si(i)
endif
j(2*i-1)=1
call dapok(roi(2*i-1),j,co(i))
simv=-sim
call dapok(roi(2*i),j,simv)
j(2*i-1)=0
j(2*i)=1
simv=-si(i)
call dapok(roi(2*i),j,co(i))
call dapok(roi(2*i-1),j,simv)
j(2*i)=0
enddo
if(ndc.eq.1) then
j(ndt)=1
call dapok(roi(ndt),j,1.d0)
call dapok(roi(ndpt),j,0.d0)
j(ndt)=0
j(ndpt)=1
call dapok(roi(ndt),j,-ang(nd))
call dapok(roi(ndpt),j,1.d0)
j(ndpt)=0
endif
return
end
subroutine hyper(a,ch,sh)
implicit none
double precision a,ch,sh,x,xi
! USED IN ROTIFLO AND ROTFLO
x=dexp(a)
xi=1.d0/x
ch=(x+xi)/2.d0
sh=(x-xi)/2.d0
return
end
subroutine ctor(c1,r2,i2)
implicit none
integer ndim2,ntt
! CHANGES OF BASIS
! C1------> R2+I R1
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer c1,r2,i2
integer b1,b2,x(ndim2)
!
!
call etallnom(b1,1,'B1 ')
call etallnom(b2,1,'B2 ')
call etallnom(x,nd2 ,'X ')
call ctoi(c1,b1)
call etcjg(x)
call trx(b1,b2,x)
call dalin(b1,.5d0,b2,.5d0,r2)
call dalin(b1,.5d0,b2,-.5d0,i2)
call dadal(x,nd2)
call dadal(b2,1)
call dadal(b1,1)
return
end
subroutine rtoc(r1,i1,c2)
implicit none
integer ndim2,ntt
! INVERSE OF CTOR
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer c2,r1,i1
integer b1
!
call etallnom(b1,1,'B1 ')
call daadd(r1,i1,b1)
call itoc(b1,c2)
call dadal(b1,1)
return
end
subroutine ctorflo(c,dr,di)
implicit none
integer ndim,ndim2,ntt
! FLOW CTOR
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer dr(*),di(*),c(*)
call ctord(c,dr,di)
call resvec(dr,di,dr,di)
return
end
subroutine rtocflo(dr,di,c)
implicit none
integer ndim,ndim2,ntt
! FLOW RTOC
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer dr(*),di(*),c(*),er(ndim2),ei(ndim2)
call etall(er,nd2)
call etall(ei,nd2)
call reelflo(dr,di,er,ei)
call rtocd(er,ei,c)
call dadal(er,nd2)
call dadal(ei,nd2)
return
end
subroutine ctord(c,cr,ci)
implicit none
integer i,ndim,ndim2,ntt
! ROUTINES USED IN THE INTERMEDIATE STEPS OF CTORFLO AND RTOCFLO
! SAME AS CTOR OVER ARRAYS CONTAINING ND2 COMPONENTS
! ROUTINE USEFUL IN INTERMEDIATE FLOW CHANGE OF BASIS
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer c(*),ci(*),cr(*)
do i=1,nd2
call ctor(c(i),cr(i),ci(i))
enddo
return
end
subroutine rtocd(cr,ci,c)
implicit none
integer i,ndim,ndim2,ntt
! INVERSE OF CTORD
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer c(*),ci(*),cr(*)
do i=1,nd2
call rtoc(cr(i),ci(i),c(i))
enddo
return
end
subroutine resvec(cr,ci,dr,di)
implicit none
integer i,ndim,ndim2,ntt
! DOES THE SPINOR PART IN CTORFLO
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
double precision angle,dsta,rad,sta
common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer dr(*),di(*),ci(*),cr(*),tr(2),ti(2)
call etall(tr,2)
call etall(ti,2)
do i=1,nd-ndc
if(ista(i).eq.1) then
call dasub(cr(2*i-1),ci(2*i),tr(1))
call daadd(ci(2*i-1),cr(2*i),ti(1))
call daadd(cr(2*i-1),ci(2*i),tr(2))
call dasub(ci(2*i-1),cr(2*i),ti(2))
call dacop(tr(1),dr(2*i-1))
call dacop(tr(2),dr(2*i))
call dacop(ti(1),di(2*i-1))
call dacop(ti(2),di(2*i))
else
call daadd(cr(2*i-1),cr(2*i),tr(1))
call daadd(ci(2*i-1),ci(2*i),ti(1))
call dasub(cr(2*i-1),cr(2*i),tr(2))
call dasub(ci(2*i-1),ci(2*i),ti(2))
call dacop(tr(1),dr(2*i-1))
call dacop(tr(2),dr(2*i))
call dacop(ti(1),di(2*i-1))
call dacop(ti(2),di(2*i))
endif
enddo
do i=nd2-ndc2+1,nd2
call dacop(cr(i),dr(i))
call dacop(ci(i),di(i))
enddo
call dadal(tr,2)
call dadal(ti,2)
return
end
subroutine reelflo(c,ci,f,fi)
implicit none
integer i,ndim,ndim2,ntt
! DOES THE SPINOR PART IN RTOCFLO
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
double precision angle,dsta,rad,sta
common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer c(*),ci(*),f(*),fi(*),e(ndim2),ei(ndim2)
call etall(e,nd2)
call etall(ei,nd2)
do i=1,nd-ndc
call dalin(c(2*i-1),0.5d0,c(2*i),0.5d0,e(2*i-1))
call dalin(ci(2*i-1),0.5d0,ci(2*i),0.5d0,ei(2*i-1))
if(ista(i).eq.1) then
call dalin(ci(2*i-1),0.5d0,ci(2*i),-0.5d0,e(2*i))
call dalin(c(2*i-1),-0.5d0,c(2*i),0.5d0,ei(2*i))
else
call dalin(ci(2*i-1),0.5d0,ci(2*i),-0.5d0,ei(2*i))
call dalin(c(2*i-1),0.5d0,c(2*i),-0.5d0,e(2*i))
endif
enddo
do i=nd2-ndc2+1,nd2
call dacop(c(i),e(i))
call dacop(ci(i),ei(i))
enddo
call dacopd(e,f)
call dacopd(ei,fi)
call dadal(e,nd2)
call dadal(ei,nd2)
return
end
subroutine compcjg(cr,ci,dr,di)
implicit none
integer ndim,ndim2,ntt
! TAKES THE COMPLEX CONJUGATE IN RESONANCE BASIS OF A POLYNOMIAL
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer dr,di,ci,cr,x(ndim2)
call etall(x,nd2)
call etcjg(x)
call trx(cr,dr,x)
call trx(ci,di,x)
call dacmu(di,-1.d0,di)
call dadal(x,nd2)
return
end
logical*1 function midbflo(c,a2,a2i,q,a,st)
implicit none
integer i,j,ndim,ndim2,ntt
double precision a,ch,cm,cr,q,r,sa,sai,shm, &
&st,x2pi
! LINEAR EXACT NORMALIZATION USING EIGENVALUE PACKAGE OF NERI
parameter (ntt=40)
parameter (ndim2=6)
parameter (ndim=3)
integer jx(ntt)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
dimension cr(ndim2,ndim2),st(ndim),q(ndim),a(ndim)
dimension sa(ndim2,ndim2),sai(ndim2,ndim2),cm(ndim2,ndim2)
integer c(*),a2(*),a2i(*)
logical*1 mapflol
!*DAEXT(NO,NV) C(NDIM2),A2(NDIM2),A2I(NDIM2)
x2pi=datan(1.d0)*8.d0
do i=1,ntt
jx(i)=0
enddo
! frank/etienne
do i=1,ndim
st(i)=0d0
q(i)=0d0
a(i)=0d0
enddo
! frank/etienne
do i=1,ndim2
! frank/etienne
do j=1,ndim2
sai(i,j)=0.d0
sa(i,j)=0.d0
cm(i,j)=0.d0
cr(i,j)=0.d0
enddo
enddo
do i=1,nd2
do j=1,nd2
jx(j)=1
call dapek(c(i),jx,r)
jx(j)=0
cm(i,j)=r
enddo
enddo
midbflo = mapflol(sa,sai,cr,cm,st)
do i=1,nd-ndc
if(st(i)+0.001.gt.1.d0) then
a(i)=dsqrt(cr(2*i-1,2*i-1)**2+cr(2*i-1,2*i)**2)
q(i)=dacos(cr(2*i-1,2*i-1)/a(i))
a(i)=dlog(a(i))
if(cr(2*i-1,2*i).lt.0.d0) q(i)=x2pi-q(i)
else
a(i)=dsqrt(cr(2*i-1,2*i-1)**2-cr(2*i-1,2*i)**2)
ch=cr(2*i-1,2*i-1)/a(i)
shm=cr(2*i-1,2*i)/a(i)
! CH=CH+DSQRT(CH**2-1.D0)
! q(i)=DLOG(CH)
q(i)=-dlog(ch+shm)
! IF(cr(2*i-1,2*i).gt.0.d0) Q(I)=-Q(I)
a(i)=dlog(a(i))
endif
enddo
if(ndc.eq.0) then
if(st(3)+0.001.gt.1.d0.and.nd.eq.3.and.q(nd).gt.0.5d0) &
& q(3)=q(3)-x2pi
else
q(nd)=cr(ndt,ndpt)
endif
call daclrd(a2)
call daclrd(a2i)
do i=1,nd2
do j=1,nd2
jx(j)=1
r=sa(i,j)
if(r.ne.0.d0)call dapok(a2(i),jx,r)
jx(j)=1
r=sai(i,j)
if(r.ne.0.d0)call dapok(a2i(i),jx,r)
jx(j)=0
enddo
enddo
return
end
logical*1 function mapflol(sa,sai,cr,cm,st)
implicit none
integer i,ier,iunst,j,l,n,n1,ndim,ndim2
double precision ap,ax,cm,cr, &
&p,rd,rd1,ri,rr,s1,sa,sai,st,vi,vr,w,x,x2pi,xd,xj,xsu,xx
parameter (ndim2=6)
parameter (ndim=3)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
!---- FROM TRACKING CODE
integer idpr
common /printing/ idpr
integer nplane
double precision epsplane,xplane
common /choice/ xplane(ndim),epsplane,nplane(ndim)
! ---------------------
dimension cr(ndim2,ndim2),xj(ndim2,ndim2),n(ndim),x(ndim)
dimension rr(ndim2),ri(ndim2),sa(ndim2,ndim2),xx(ndim) &
&,sai(ndim2,ndim2),cm(ndim2,ndim2),w(ndim2,ndim2),st(ndim)
dimension vr(ndim2,ndim2),vi(ndim2,ndim2),s1(ndim2,ndim2),p(ndim2)
logical*1 eig6
x2pi=datan(1.d0)*8.d0
n1=0
! frank/etienne
do i=1,ndim2
do j=1,ndim2
cr(j,i)=cm(i,j)
xj(i,j)=0.d0
s1(i,j)=0.d0
enddo
enddo
! frank/etienne
do i=1,ndim
n(i)=0
xj(2*i-1,2*i)=1.d0
xj(2*i,2*i-1)=-1.d0
enddo
! frank/etienne
do i=1,ndim2
do j=1,ndim2
sai(i,j)=0.d0
w(i,j)=cm(i,j)
enddo
enddo
if(ndc.eq.1) then
s1(nd2-ndc,nd2-ndc)=1.d0
s1(nd2,nd2)=1.d0
sai(nd2-ndc,nd2-ndc)=1.d0
sai(nd2,nd2)=1.d0
endif
call mulnd2(xj,w)
call mulnd2(cr,w)
if(idpr.ge.0.or.idpr.eq.-102) then
write(6,*)'Check of the symplectic condition on the linear part'
xsu=0.d0
do i=1,nd2
write(6,'(6(2x,g23.16))') ( w(i,j), j = 1, nd2 )
do j=1,nd2
xsu=xsu+dabs(w(i,j))
enddo
enddo
write(6,*)'deviation for symplecticity ',100.d0*(xsu-nd2)/xsu, &
& ' %'
endif
mapflol = eig6(cr,rr,ri,vr,vi)
if(idpr.ge.0) then
write(6,*) ' '
write(6,*) ' Index Real Part ', &
& ' ArcSin(Imaginary Part)/2/pi'
write(6,*) ' '
do i=1,nd-ndc
rd1=dsqrt(rr(2*i-1)**2+ri(2*i-1)**2)
rd=dsqrt(rr(2*i)**2+ri(2*i)**2)
write(6,*) 2*i-1,rr(2*i-1),dasin(ri(2*i-1)/rd1)/x2pi
write(6,*) 2*i,rr(2*i),dasin(ri(2*i)/rd)/x2pi
write(6,*) ' alphas ', dlog(dsqrt(rd*rd1))
enddo
if ( idpr.ge. 0) then
write(6,*) &
& ' select ',nd-ndc, &
& ' eigenplanes (odd integers <0 real axis)'
! read(5,*) (n(i),i=1,nd-ndc)
n(1) = 1
n(2) = 3
n(3) = 5
else
n(1) = 1
n(2) = 3
n(3) = 5
endif
elseif(idpr.eq.-100) then
do i=1,nd-ndc
n(i)=nplane(i)
enddo
elseif(idpr.eq.-101.or.idpr.eq.-102) then
do i=1,nd-ndc
if(ri(2*i).ne.0.d0) then
n(i)=2*i-1
else
n(i)=-2*i+1
endif
enddo
else
do i=1,nd-ndc
n(i)=2*i-1
enddo
endif
iunst=0
do i=1,nd-ndc ! Frank NDC kept
if(n(i).lt.0) then
n(i)=-n(i)
st(i)=0.d0
iunst=1
else
st(i)=1.d0
endif
x(i)=0.d0
xx(i)=1.d0
do j=1,nd-ndc
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)
enddo
enddo
do i=1,nd-ndc
if(x(i).lt.0.d0) xx(i)=-1.d0
x(i)=dsqrt(dabs(x(i)))
enddo
do i=1,nd2-ndc2
do j=1,nd-ndc
if(st(j)+0.001.gt.1.d0) then
sai(2*j-1,i)=vr(i,n(j))*xx(j)/x(j)
sai(2*j,i)=vi(i,n(j))/x(j)
else
ax=vr(i,n(j))*xx(j)/x(j)
ap=vi(i,n(j))/x(j)
sai(2*j-1,i)=(ax+ap)/dsqrt(2.d0)
sai(2*j,i)=(ap-ax)/dsqrt(2.d0)
endif
enddo
enddo
if(idpr.eq.-101.or.idpr.eq.-102) then
call movearou(sai)
endif
! adjust sa such that sa(1,2)=0 and sa(3,4)=0. (courant-snyder-edwards-teng
! phase advances)
if(iunst.ne.1) then
do i=1,nd-ndc
p(i)=datan(-sai(2*i-1,2*i)/sai(2*i,2*i))
s1(2*i-1,2*i-1)=dcos(p(i))
s1(2*i,2*i)=dcos(p(i))
s1(2*i-1,2*i)=dsin(p(i))
s1(2*i,2*i-1)=-dsin(p(i))
enddo
call mulnd2(s1,sai)
! adjust sa to have sa(1,1)>0 and sa(3,3)>0 rotate by pi if necessary.
do i=1,nd-ndc
xd=1.d0
if(sai(2*i-1,2*i-1).lt.0.d0) xd=-1.d0
s1(2*i-1,2*i-1)=xd
s1(2*i-1,2*i)=0.d0
s1(2*i,2*i-1)=0.d0
s1(2*i,2*i)=xd
enddo
call mulnd2(s1,sai)
! sa is now uniquely and unambigeously determined.
endif
do i=1,nd2
do l=1,nd2
sa(i,l)=sai(i,l)
enddo
enddo
call matinv(sai,sa,nd2,6,ier)
call mulnd2(sai,cm)
do i=1,nd2
do j=1,nd2
cr(i,j)=sa(i,j)
enddo
enddo
call mulnd2(cm,cr)
return
end
subroutine mulnd2(rt,r)
implicit none
integer i,ia,j,ndim,ndim2
double precision r,rt,rtt
parameter (ndim2=6)
parameter (ndim=3)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
dimension rt(ndim2,ndim2),r(ndim2,ndim2),rtt(ndim2,ndim2)
do i=1,nd2
do j=1,nd2
rtt(i,j)=0.d0
enddo
enddo
do i=1,nd2
do j=1,nd2
do ia=1,nd2
rtt(i,ia)=rt(i,j)*r(j,ia)+rtt(i,ia)
enddo
enddo
enddo
do i=1,nd2
do j=1,nd2
r(i,j)=rtt(i,j)
enddo
enddo
return
end
subroutine movearou(rt)
implicit none
integer ipause, mypause
integer i,ic,j,ndim,ndim2
double precision rt,rto,s,xr,xrold,xy,xyz,xz,xzy,yz
parameter (ndim2=6)
parameter (ndim=3)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer idpr
common /printing/ idpr
dimension rt(ndim2,ndim2),rto(ndim2,ndim2)
dimension xy(ndim2,ndim2),xz(ndim2,ndim2),yz(ndim2,ndim2)
dimension xyz(ndim2,ndim2),xzy(ndim2,ndim2)
dimension s(ndim2,ndim2)
do i=1,nd2
do j=1,nd2
s(i,j)=0.d0
s(i,i)=1.d0
xy(i,j)=0.d0
xz(i,j)=0.d0
yz(i,j)=0.d0
xyz(i,j)=0.d0
xzy(i,j)=0.d0
enddo
enddo
xy(1,3)=1.d0
xy(3,1)=1.d0
xy(2,4)=1.d0
xy(4,2)=1.d0
xy(5,5)=1.d0
xy(6,6)=1.d0
xz(1,5)=1.d0
xz(5,1)=1.d0
xz(2,6)=1.d0
xz(6,2)=1.d0
xz(3,3)=1.d0
xz(4,4)=1.d0
yz(3,5)=1.d0
yz(5,3)=1.d0
yz(4,6)=1.d0
yz(6,4)=1.d0
yz(1,1)=1.d0
yz(2,2)=1.d0
xyz(1,3)=1.d0
xyz(3,5)=1.d0
xyz(5,1)=1.d0
xyz(2,4)=1.d0
xyz(4,6)=1.d0
xyz(6,2)=1.d0
xzy(1,5)=1.d0
xzy(5,3)=1.d0
xzy(3,1)=1.d0
xzy(2,6)=1.d0
xzy(6,4)=1.d0
xzy(4,2)=1.d0
ic=0
xrold=1000000000.d0
call movemul(rt,s,rto,xr)
! write(6,*) xr,xrold
! do i=1,6
! write(6,'(6(1x,1pe12.5))') (RTO(i,j),j=1,6)
! enddo
! ipause=mypause(0)
if(xr.lt.xrold) then
xrold=xr
endif
if(nd.ge.2) then
call movemul(rt,xy,rto,xr)
if(xr.lt.xrold) then
xrold=xr
ic=1
endif
endif
if(nd.eq.3) then
call movemul(rt,xz,rto,xr)
if(xr.lt.xrold) then
xrold=xr
ic=2
endif
call movemul(rt,yz,rto,xr)
if(xr.lt.xrold) then
xrold=xr
ic=3
endif
call movemul(rt,xyz,rto,xr)
if(xr.lt.xrold) then
xrold=xr
ic=4
endif
call movemul(rt,xzy,rto,xr)
if(xr.lt.xrold) then
xrold=xr
ic=5
endif
endif
if(ic.eq.0) then
call movemul(rt,s,rto,xr)
if(idpr.gt.-101) write(6,*) " no exchanged"
elseif(ic.eq.1) then
call movemul(rt,xy,rto,xr)
if(idpr.gt.-101) write(6,*) " x-y exchanged"
elseif(ic.eq.2) then
call movemul(rt,xz,rto,xr)
if(idpr.gt.-101) write(6,*) " x-z exchanged"
elseif(ic.eq.3) then
call movemul(rt,yz,rto,xr)
if(idpr.gt.-101) write(6,*) " y-z exchanged"
elseif(ic.eq.4) then
call movemul(rt,xyz,rto,xr)
if(idpr.gt.-101) write(6,*) " x-y-z permuted"
elseif(ic.eq.5) then
call movemul(rt,xzy,rto,xr)
if(idpr.gt.-101) write(6,*) " x-z-y permuted"
endif
do i=1,nd2
do j=1,nd2
rt(i,j)=rto(i,j)
enddo
enddo
return
end
subroutine movemul(rt,xy,rto,xr)
implicit none
integer i,j,k,ndim,ndim2
double precision rt,rto,xr,xy
parameter (ndim2=6)
parameter (ndim=3)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
dimension rt(ndim2,ndim2)
dimension xy(ndim2,ndim2),rto(ndim2,ndim2)
do i=1,nd2
do j=1,nd2
rto(i,j)=0.d0
enddo
enddo
do i=1,nd2
do j=1,nd2
do k=1,nd2
rto(i,k)=xy(i,j)*rt(j,k)+rto(i,k)
enddo
enddo
enddo
xr=0.d0
do i=1,nd2
do j=1,nd2
xr=xr+dabs(rto(i,j))
enddo
enddo
do i=1,nd
xr=xr-dabs(rto(2*i-1,2*i-1))
xr=xr-dabs(rto(2*i-1,2*i))
xr=xr-dabs(rto(2*i,2*i))
xr=xr-dabs(rto(2*i,2*i-1))
enddo
return
end
subroutine initpert(st,ang,ra)
implicit none
integer i,ndim,ndim2,nn,nreso
double precision ang,ra,st
! X-RATED
!- SETS UP ALL THE COMMON BLOCKS RELEVANT TO NORMAL FORM AND THE BASIS
!- CHANGES INSIDE MAPNORMF
parameter (ndim=3)
parameter (ndim2=6)
parameter (nreso=20)
dimension st(ndim),ang(ndim),ra(ndim)
double precision angle,dsta,rad,sta
common /stable/sta(ndim),dsta(ndim),angle(ndim),rad(ndim)
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer mx,nres
common /reson/mx(ndim,nreso),nres
integer iref
common /resfile/iref
if(iref.gt.0) then
write(6,*) iref
read(iref,*) nres
if(nres.ge.nreso) then
write(6,*) ' NRESO IN LIELIB TOO SMALL '
stop999
endif
elseif(iref.eq.0) then
nres=0
endif
if(nres.ne.0) write(6,*)' warning resonances left in the map'
if(iref.gt.0) then
do i=1,nres
read(iref,*) (mx(nn,i),nn=1,nd-ndc)
enddo
endif
do i=nres+1,nreso
do nn=1,ndim
mx(nn,i)=0
enddo
enddo
! frank/Etienne
do i=1,ndim
angle(i)=0.d0
rad(i)=0.d0
sta(i)=0.d0
dsta(i)=1.d0-sta(i)
ista(i)=0
idsta(i)=0
enddo
do i=1,nd ! Frank -ndc
angle(i)=ang(i)
rad(i)=ra(i)
sta(i)=st(i)
dsta(i)=1.d0-sta(i)
enddo
do i=1,nd
ista(i)=idint(sta(i)+.01)
idsta(i)=idint(dsta(i)+.01)
enddo
return
end
double precision function dlie(j)
implicit none
integer i,ndim
parameter (ndim=3)
! PARAMETER (NTT=40)
! INTEGER J(NTT)
integer j(*)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
dlie=0.d0
do i=1,nd
dlie=dble(j(2*i-1)+j(2*i))+dlie
enddo
dlie=dlie+1.d0
dlie=1.d0/dlie
return
end
double precision function rext(j)
implicit none
integer i,lie,mo,ndim
parameter (ndim=3)
! PARAMETER (NTT=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
integer j(*)
lie=0
do i=1,nd-ndc
lie=ista(i)*j(2*i)+lie
enddo
mo=mod(lie,4)+1
goto(11,12,13,14),mo
11 rext = 1.d0
return
12 rext = -1.d0
return
13 rext = -1.d0
return
14 rext = 1.d0
return
end
subroutine cpart(h,ch)
implicit none
integer ndim,ntt
double precision rext
parameter (ndim=3)
parameter (ntt=40)
external rext
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer h,ch
call dacfu(h,rext,ch)
return
end
subroutine ctoi(f1,f2)
implicit none
integer ndim2,ntt
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer f1,f2
integer b1,x(ndim2)
!
!
call etallnom(b1,1,'B1 ')
call etallnom(x,nd2 ,'X ')
call cpart(f1,b1)
call etctr(x)
call trx(b1,f2,x)
call dadal(x,nd2)
call dadal(b1,1)
return
end
subroutine itoc(f1,f2)
implicit none
integer ndim2,ntt
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer f1,f2
integer b1,x(ndim2)
!
call etallnom(b1,1,'B1 ')
call etallnom(x,nd2 ,'X ')
call etrtc(x)
call trx(f1,b1,x)
call cpart(b1,f2)
call dadal(x,nd2)
call dadal(b1,1)
return
end
subroutine etrtc(x)
implicit none
integer i,ndim,ndim2,ntt
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer x(*)
integer rel(ndim2)
!
!
call etallnom(rel,nd2 ,'REL ')
call etini(rel)
call etini(x)
do i=1,nd-ndc
call daadd(rel(2*i-1),rel(2*i),x(2*i-1))
call dasub(rel(2*i-1),rel(2*i),x(2*i))
enddo
call dadal(rel,nd2)
return
end
subroutine etctr(x)
implicit none
integer i,ndim,ndim2,ntt
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer x(*)
integer rel(ndim2)
!
!
call etallnom(rel,nd2 ,'REL ')
call etini(rel)
call etini(x)
do i=1,nd-ndc
call dalin(rel(2*i-1),.5d0,rel(2*i),.5d0,x(2*i-1))
call dalin(rel(2*i-1),.5d0,rel(2*i),-.5d0,x(2*i))
enddo
call dadal(rel,nd2)
return
end
subroutine etcjg(x)
implicit none
integer i,ndim,ndim2,ntt
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer idsta,ista
common /istable/ista(ndim),idsta(ndim)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer x(*)
integer rel(ndim2)
!
!
call etallnom(rel,nd2 ,'REL ')
call etini(rel)
call etini(x)
do i=1,nd-ndc
if(ista(i).eq.1) then
call dacop(rel(2*i-1),x(2*i))
call dacop(rel(2*i),x(2*i-1))
else
call dacop(rel(2*i-1),x(2*i-1))
call dacop(rel(2*i),x(2*i))
endif
enddo
call dadal(rel,nd2)
return
end
logical*1 function eig6(fm,reval,aieval,revec,aievec)
implicit none
integer jet,ndim2
!**************************************************************************
! Diagonalization routines of NERI
!ccccccccccccccccc
!
! this routine finds the eigenvalues and eigenvectors
! of the full matrix fm.
! the eigenvectors are normalized so that the real and
! imaginary part of vectors 1, 3, and 5 have +1 antisymmetric
! product:
! revec1 J aivec1 = 1 ; revec3 J aivec3 = 1 ;
! revec5 J aivec5 = 1.
! the eigenvectors 2 ,4, and 6 have the opposite normalization.
! written by F. Neri, Feb 26 1986.
!
parameter (ndim2=6)
integer nn
integer ilo,ihi,mdim,info
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
double precision reval(ndim2),aieval(ndim2), &
&revec(ndim2,ndim2),aievec(ndim2,ndim2)
double precision fm(ndim2,ndim2),aa(ndim2,ndim2)
integer i,i1
double precision ort(ndim2),vv(ndim2,ndim2)
eig6 = .true.
! copy matrix to temporary storage (the matrix aa is destroyed)
do i=1,nd2-ndc2
do i1=1,nd2-ndc2
aa(i1,i) = fm(i1,i)
enddo
enddo
ilo = 1
ihi = nd2-ndc2
mdim = ndim2
nn = nd2-ndc2
! compute eigenvalues and eigenvectors using double
! precision Eispack routines:
call ety(mdim,nn,ilo,ihi,aa,ort)
call etyt(mdim,nn,ilo,ihi,aa,ort,vv)
call ety2(mdim,nn,ilo,ihi,aa,reval,aieval,vv,info)
if ( info .ne. 0 ) then
write(6,*) ' ERROR IN EIG6'
eig6 = .false.
return
endif
! call neigv(vv,pbkt)
do i=1,nd-ndc
do jet=1,nd2-ndc2
revec(jet,2*i-1)=vv(jet,2*i-1)
revec(jet,2*i)=vv(jet,2*i-1)
aievec(jet,2*i-1)=vv(jet,2*i)
aievec(jet,2*i)=-vv(jet,2*i)
enddo
enddo
do i=1,nd2-ndc2
if(dabs(reval(i)**2+aieval(i)**2 -1.d0).gt.1.d-10) then
write(6,*) ' EIG6: Eigenvalues off the unit circle!'
eig6 = .false.
endif
enddo
return
end
subroutine ety(nm,n,low,igh,a,ort)
implicit none
integer i,j,m,n,ii,jj,la,mp,nm,igh,kp1,low
double precision a(nm,n),ort(igh)
double precision f,g,h,scale
!
! this subroutine is a translation of the algol procedure orthes,
! num. math. 12, 349-368(1968) by martin and wilkinson.
! handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).
!
! given a real general matrix, this subroutine
! reduces a submatrix situated in rows and columns
! low through igh to upper hessenberg form by
! orthogonal similarity transformations.
!
! on input-
!
! nm must be set to the row dimension of two-dimensional
! array parameters as declared in the calling program
! dimension statement,
!
! n is the order of the matrix,
!
! low and igh are integers determined by the balancing
! subroutine balanc. if balanc has not been used,
! set low=1, igh=n,
!
! a contains the input matrix.
!
! on output-
!
! a contains the hessenberg matrix. information about
! the orthogonal transformations used in the reduction
! is stored in the remaining triangle under the
! hessenberg matrix,
!
! ort contains further information about the transformations.
! only elements low through igh are used.
!
! fortran routine by b. s. garbow
! modified by filippo neri.
!
!
la = igh - 1
kp1 = low + 1
if (la .lt. kp1) go to 200
!
do m = kp1, la
h = 0.0
ort(m) = 0.0
scale = 0.0
! ********** scale column (algol tol then not needed) **********
do i = m, igh
scale = scale + dabs(a(i,m-1))
enddo
!
if (scale .eq. 0.0) go to 180
mp = m + igh
! ********** for i=igh step -1 until m do -- **********
do ii = m, igh
i = mp - ii
ort(i) = a(i,m-1) / scale
h = h + ort(i) * ort(i)
enddo
!
g = -dsign(dsqrt(h),ort(m))
h = h - ort(m) * g
ort(m) = ort(m) - g
! ********** form (i-(u*ut)/h) * a **********
do j = m, n
f = 0.0
! ********** for i=igh step -1 until m do -- **********
do ii = m, igh
i = mp - ii
f = f + ort(i) * a(i,j)
enddo
!
f = f / h
!
do i = m, igh
a(i,j) = a(i,j) - f * ort(i)
enddo
!
enddo
! ********** form (i-(u*ut)/h)*a*(i-(u*ut)/h) **********
do i = 1, igh
f = 0.0
! ********** for j=igh step -1 until m do -- **********
do jj = m, igh
j = mp - jj
f = f + ort(j) * a(i,j)
enddo
!
f = f / h
!
do j = m, igh
a(i,j) = a(i,j) - f * ort(j)
enddo
!
enddo
!
ort(m) = scale * ort(m)
a(m,m-1) = scale * g
180 continue
enddo
!
200 return
! ********** last card of ety **********
end
subroutine etyt(nm,n,low,igh,a,ort,z)
implicit none
integer i,j,n,kl,mm,mp,nm,igh,low,mp1
double precision a(nm,igh),ort(igh),z(nm,n)
double precision g
!
! this subroutine is a translation of the algol procedure ortrans,
! num. math. 16, 181-204(1970) by peters and wilkinson.
! handbook for auto. comp., vol.ii-linear algebra, 372-395(1971).
!
! this subroutine accumulates the orthogonal similarity
! transformations used in the reduction of a real general
! matrix to upper hessenberg form by ety.
!
! on input-
!
! nm must be set to the row dimension of two-dimensional
! array parameters as declared in the calling program
! dimension statement,
!
! n is the order of the matrix,
!
! low and igh are integers determined by the balancing
! subroutine balanc. if balanc has not been used,
! set low=1, igh=n,
!
! a contains information about the orthogonal trans-
! formations used in the reduction by orthes
! in its strict lower triangle,
!
! ort contains further information about the trans-
! formations used in the reduction by ety.
! only elements low through igh are used.
!
! on output-
!
! z contains the transformation matrix produced in the
! reduction by ety,
!
! ort has been altered.
!
! fortran routine by b. s. garbow.
! modified by f. neri.
!
!
! ********** initialize z to identity matrix **********
do i = 1, n
!
do j = 1, n
z(i,j) = 0.0
enddo
!
z(i,i) = 1.0
enddo
!
kl = igh - low - 1
if (kl .lt. 1) go to 200
! ********** for mp=igh-1 step -1 until low+1 do -- **********
do mm = 1, kl
mp = igh - mm
if (a(mp,mp-1) .eq. 0.0) go to 140
mp1 = mp + 1
!
do i = mp1, igh
ort(i) = a(i,mp-1)
enddo
!
do j = mp, igh
g = 0.0
!
do i = mp, igh
g = g + ort(i) * z(i,j)
enddo
! ********** divisor below is negative of h formed in orthes.
! double division avoids possible underflow **********
g = (g / ort(mp)) / a(mp,mp-1)
!
do i = mp, igh
z(i,j) = z(i,j) + g * ort(i)
enddo
!
enddo
!
140 continue
enddo
!
200 return
! ********** last card of etyt **********
end
subroutine ety2(nm,n,low,igh,h,wr,wi,z,ierr)
implicit none
integer i,j,k,l,m,n,en,ii,jj,ll,mm,na,nm,nn, &
&igh,its,low,mp2,enm2,ierr
double precision h(nm,n),wr(n),wi(n),z(nm,n)
double precision p,q,r,s,t,w,x,y,ra,sa,vi,vr,zz,norm,machep
logical notlas
double precision z3r,z3i
!
!
!
! this subroutine is a translation of the algol procedure hqr2,
! num. math. 16, 181-204(1970) by peters and wilkinson.
! handbook for auto. comp., vol.ii-linear algebra, 372-395(1971).
!
! this subroutine finds the eigenvalues and eigenvectors
! of a real upper hessenberg matrix by the qr method. the
! eigenvectors of a real general matrix can also be found
! if elmhes and eltran or orthes and ortran have
! been used to reduce this general matrix to hessenberg form
! and to accumulate the similarity transformations.
!
! on input-
!
! nm must be set to the row dimension of two-dimensional
! array parameters as declared in the calling program
! dimension statement,
!
! n is the order of the matrix,
!
! low and igh are integers determined by the balancing
! subroutine balanc. if balanc has not been used,
! set low=1, igh=n,
!
! h contains the upper hessenberg matrix,
!
! z contains the transformation matrix produced by eltran
! after the reduction by elmhes, or by ortran after the
! reduction by orthes, if performed. if the eigenvectors
! of the hessenberg matrix are desired, z must contain the
! identity matrix.
!
! on output-
!
! h has been destroyed,
!
! wr and wi contain the real and imaginary parts,
! respectively, of the eigenvalues. the eigenvalues
! are unordered except that complex conjugate pairs
! of values appear consecutively with the eigenvalue
! having the positive imaginary part first. if an
! error exit is made, the eigenvalues should be correct
! for indices ierr+1,...,n,
!
! z contains the real and imaginary parts of the eigenvectors.
! if the i-th eigenvalue is real, the i-th column of z
! contains its eigenvector. if the i-th eigenvalue is complex
! with positive imaginary part, the i-th and (i+1)-th
! columns of z contain the real and imaginary parts of its
! eigenvector. the eigenvectors are unnormalized. if an
! error exit is made, none of the eigenvectors has been found,
!
! ierr is set to
! zero for normal return,
! j if the j-th eigenvalue has not been
! determined after 200 iterations.
!
! arithmetic is double precision. complex division
! is simulated by routin etdiv.
!
! fortran routine by b. s. garbow.
! modified by f. neri.
!
!
! ********** machep is a machine dependent parameter specifying
! the relative precision of floating point arithmetic.
!
! **********
machep = 1.d-17
! machep = r1mach(4)
!
ierr = 0
norm = 0.0
k = 1
! ********** store roots isolated by balanc
! and compute matrix norm **********
do i = 1, n
!
do j = k, n
norm = norm + dabs(h(i,j))
enddo
!
k = i
if (i .ge. low .and. i .le. igh) go to 50
wr(i) = h(i,i)
wi(i) = 0.0
50 continue
enddo
!
en = igh
t = 0.0
! ********** search for next eigenvalues **********
60 if (en .lt. low) go to 340
its = 0
na = en - 1
enm2 = na - 1
! ********** look for single small sub-diagonal element
! for l=en step -1 until low do -- **********
70 do ll = low, en
l = en + low - ll
if (l .eq. low) go to 100
s = dabs(h(l-1,l-1)) + dabs(h(l,l))
if (s .eq. 0.0) s = norm
if (dabs(h(l,l-1)) .le. machep * s) go to 100
enddo
! ********** form shift **********
100 x = h(en,en)
if (l .eq. en) go to 270
y = h(na,na)
w = h(en,na) * h(na,en)
if (l .eq. na) go to 280
if (its .eq. 200) go to 1000
if (its .ne. 10 .and. its .ne. 20) go to 130
! ********** form exceptional shift **********
t = t + x
!
do i = low, en
h(i,i) = h(i,i) - x
enddo
!
s = dabs(h(en,na)) + dabs(h(na,enm2))
x = 0.75 * s
y = x
w = -0.4375 * s * s
130 its = its + 1
! ********** look for two consecutive small
! sub-diagonal elements.
! for m=en-2 step -1 until l do -- **********
do mm = l, enm2
m = enm2 + l - mm
zz = h(m,m)
r = x - zz
s = y - zz
p = (r * s - w) / h(m+1,m) + h(m,m+1)
q = h(m+1,m+1) - zz - r - s
r = h(m+2,m+1)
s = dabs(p) + dabs(q) + dabs(r)
p = p / s
q = q / s
r = r / s
if (m .eq. l) go to 150
if (dabs(h(m,m-1)) * (dabs(q) + dabs(r)) .le. machep * dabs(p) &
& * (dabs(h(m-1,m-1)) + dabs(zz) + dabs(h(m+1,m+1)))) go to 150
enddo
!
150 mp2 = m + 2
!
do i = mp2, en
h(i,i-2) = 0.0
if (i .eq. mp2) go to 160
h(i,i-3) = 0.0
160 continue
enddo
! ********** double qr step involving rows l to en and
! columns m to en **********
do k = m, na
notlas = k .ne. na
if (k .eq. m) go to 170
p = h(k,k-1)
q = h(k+1,k-1)
r = 0.0
if (notlas) r = h(k+2,k-1)
x = dabs(p) + dabs(q) + dabs(r)
if (x .eq. 0.0) go to 260
p = p / x
q = q / x
r = r / x
170 s = dsign(dsqrt(p*p+q*q+r*r),p)
if (k .eq. m) go to 180
h(k,k-1) = -s * x
go to 190
180 if (l .ne. m) h(k,k-1) = -h(k,k-1)
190 p = p + s
x = p / s
y = q / s
zz = r / s
q = q / p
r = r / p
! ********** row modification **********
do j = k, n
p = h(k,j) + q * h(k+1,j)
if (.not. notlas) go to 200
p = p + r * h(k+2,j)
h(k+2,j) = h(k+2,j) - p * zz
200 h(k+1,j) = h(k+1,j) - p * y
h(k,j) = h(k,j) - p * x
enddo
!
j = min0(en,k+3)
! ********** column modification **********
do i = 1, j
p = x * h(i,k) + y * h(i,k+1)
if (.not. notlas) go to 220
p = p + zz * h(i,k+2)
h(i,k+2) = h(i,k+2) - p * r
220 h(i,k+1) = h(i,k+1) - p * q
h(i,k) = h(i,k) - p
enddo
! ********** accumulate transformations **********
do i = low, igh
p = x * z(i,k) + y * z(i,k+1)
if (.not. notlas) go to 240
p = p + zz * z(i,k+2)
z(i,k+2) = z(i,k+2) - p * r
240 z(i,k+1) = z(i,k+1) - p * q
z(i,k) = z(i,k) - p
enddo
!
260 continue
enddo
!
go to 70
! ********** one root found **********
270 h(en,en) = x + t
wr(en) = h(en,en)
wi(en) = 0.0
en = na
go to 60
! ********** two roots found **********
280 p = (y - x) / 2.0
q = p * p + w
zz = dsqrt(dabs(q))
h(en,en) = x + t
x = h(en,en)
h(na,na) = y + t
if (q .lt. 0.0) go to 320
! ********** real pair **********
zz = p + dsign(zz,p)
wr(na) = x + zz
wr(en) = wr(na)
if (zz .ne. 0.0) wr(en) = x - w / zz
wi(na) = 0.0
wi(en) = 0.0
x = h(en,na)
s = dabs(x) + dabs(zz)
p = x / s
q = zz / s
r = dsqrt(p*p+q*q)
p = p / r
q = q / r
! ********** row modification **********
do j = na, n
zz = h(na,j)
h(na,j) = q * zz + p * h(en,j)
h(en,j) = q * h(en,j) - p * zz
enddo
! ********** column modification **********
do i = 1, en
zz = h(i,na)
h(i,na) = q * zz + p * h(i,en)
h(i,en) = q * h(i,en) - p * zz
enddo
! ********** accumulate transformations **********
do i = low, igh
zz = z(i,na)
z(i,na) = q * zz + p * z(i,en)
z(i,en) = q * z(i,en) - p * zz
enddo
!
go to 330
! ********** complex pair **********
320 wr(na) = x + p
wr(en) = x + p
wi(na) = zz
wi(en) = -zz
330 en = enm2
go to 60
! ********** all roots found. backsubstitute to find
! vectors of upper triangular form **********
340 if (norm .eq. 0.0) go to 1001
! ********** for en=n step -1 until 1 do -- **********
do nn = 1, n
en = n + 1 - nn
p = wr(en)
q = wi(en)
na = en - 1
if (q.lt.0) goto 710
if (q.eq.0) goto 600
if (q.gt.0) goto 800
! ********** real vector **********
600 m = en
h(en,en) = 1.0
if (na .eq. 0) go to 800
! ********** for i=en-1 step -1 until 1 do -- **********
do ii = 1, na
i = en - ii
w = h(i,i) - p
r = h(i,en)
if (m .gt. na) go to 620
!
do j = m, na
r = r + h(i,j) * h(j,en)
enddo
!
620 if (wi(i) .ge. 0.0) go to 630
zz = w
s = r
go to 700
630 m = i
if (wi(i) .ne. 0.0) go to 640
t = w
if (w .eq. 0.0) t = machep * norm
h(i,en) = -r / t
go to 700
! ********** solve real equations **********
640 x = h(i,i+1)
y = h(i+1,i)
q = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i)
t = (x * s - zz * r) / q
h(i,en) = t
if (dabs(x) .le. dabs(zz)) go to 650
h(i+1,en) = (-r - w * t) / x
go to 700
650 h(i+1,en) = (-s - y * t) / zz
700 continue
enddo
! ********** end real vector **********
go to 800
! ********** complex vector **********
710 m = na
! ********** last vector component chosen imaginary so that
! eigenvector matrix is triangular **********
if (dabs(h(en,na)) .le. dabs(h(na,en))) go to 720
h(na,na) = q / h(en,na)
h(na,en) = -(h(en,en) - p) / h(en,na)
go to 730
! 720 z3 = cmplx(0.0,-h(na,en)) / cmplx(h(na,na)-p,q)
! h(na,na) = real(z3)
! h(na,en) = aimag(z3)
720 call etdiv(z3r,z3i,0.d0,-h(na,en),h(na,na)-p,q)
h(na,na) = z3r
h(na,en) = z3i
730 h(en,na) = 0.0
h(en,en) = 1.0
enm2 = na - 1
if (enm2 .eq. 0) go to 800
! ********** for i=en-2 step -1 until 1 do -- **********
do ii = 1, enm2
i = na - ii
w = h(i,i) - p
ra = 0.0
sa = h(i,en)
!
do j = m, na
ra = ra + h(i,j) * h(j,na)
sa = sa + h(i,j) * h(j,en)
enddo
!
if (wi(i) .ge. 0.0) go to 770
zz = w
r = ra
s = sa
go to 790
770 m = i
if (wi(i) .ne. 0.0) go to 780
! z3 = cmplx(-ra,-sa) / cmplx(w,q)
! h(i,na) = real(z3)
! h(i,en) = aimag(z3)
call etdiv(z3r,z3i,-ra,-sa,w,q)
h(i,na) = z3r
h(i,en) = z3i
go to 790
! ********** solve complex equations **********
780 x = h(i,i+1)
y = h(i+1,i)
vr = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) - q * q
vi = (wr(i) - p) * 2.0 * q
if (vr .eq. 0.0 .and. vi .eq. 0.0) vr = machep * norm &
& * (dabs(w) + dabs(q) + dabs(x) + dabs(y) + dabs(zz))
! z3 = cmplx(x*r-zz*ra+q*sa,x*s-zz*sa-q*ra) / cmplx(vr,vi)
! h(i,na) = real(z3)
! h(i,en) = aimag(z3)
call etdiv(z3r,z3i,x*r-zz*ra+q*sa,x*s-zz*sa-q*ra,vr,vi)
h(i,na) = z3r
h(i,en) = z3i
if (dabs(x) .le. dabs(zz) + dabs(q)) go to 785
h(i+1,na) = (-ra - w * h(i,na) + q * h(i,en)) / x
h(i+1,en) = (-sa - w * h(i,en) - q * h(i,na)) / x
go to 790
! 785 z3 = cmplx(-r-y*h(i,na),-s-y*h(i,en)) / cmplx(zz,q)
! h(i+1,na) = real(z3)
! h(i+1,en) = aimag(z3)
785 call etdiv(z3r,z3i,-r-y*h(i,na),-s-y*h(i,en),zz,q)
h(i+1,na) = z3r
h(i+1,en) = z3i
790 continue
enddo
! ********** end complex vector **********
800 continue
enddo
! ********** end back substitution.
! vectors of isolated roots **********
do i = 1, n
if (i .ge. low .and. i .le. igh) go to 840
!
do j = i, n
z(i,j) = h(i,j)
enddo
!
840 continue
enddo
! ********** multiply by transformation matrix to give
! vectors of original full matrix.
! for j=n step -1 until low do -- **********
do jj = low, n
j = n + low - jj
m = min0(j,igh)
!
do i = low, igh
zz = 0.0
!
do k = low, m
zz = zz + z(i,k) * h(k,j)
enddo
!
z(i,j) = zz
enddo
enddo
!
go to 1001
! ********** set error -- no convergence to an
! eigenvalue after 200 iterations **********
1000 ierr = en
1001 return
! ********** last card of ety2 **********
end
subroutine etdiv(a,b,c,d,e,f)
implicit none
! computes the complex division
! a + ib = (c + id)/(e + if)
! very slow, but tries to be as accurate as
! possible by changing the order of the
! operations, so to avoid under(over)flow
! problems.
! Written by F. Neri Feb. 12 1986
!
double precision a,b,c,d,e,f
double precision s,t
double precision cc,dd,ee,ff
double precision temp
integer flip
flip = 0
cc = c
dd = d
ee = e
ff = f
if( dabs(f).ge.dabs(e) ) then
ee = f
ff = e
cc = d
dd = c
flip = 1
endif
s = 1.d0/ee
t = 1.d0/(ee+ ff*(ff*s))
if ( dabs(ff) .ge. dabs(s) ) then
temp = ff
ff = s
s = temp
endif
if( dabs(dd) .ge. dabs(s) ) then
a = t*(cc + s*(dd*ff))
else if ( dabs(dd) .ge. dabs(ff) ) then
a = t*(cc + dd*(s*ff))
else
a = t*(cc + ff*(s*dd))
endif
if ( dabs(cc) .ge. dabs(s)) then
b = t*(dd - s*(cc*ff))
else if ( dabs(cc) .ge. dabs(ff)) then
b = t*(dd - cc*(s*ff))
else
b = t*(dd - ff*(s*cc))
endif
if (flip.ne.0 ) then
b = -b
endif
return
end
subroutine sympl3(m)
!**********************************************************
!
! SYMPL3
!
!
! On return ,the matrix m(*,*), supposed to be almost
! symplectic on entry is made exactly symplectic by
! using a non iterative, constructive method.
!
!**********************************************************
!
! Written by F. Neri Feb 7 1986
!
implicit none
integer n
parameter ( n = 3 )
integer kp,kq,lp,lq,jp,jq,i
double precision m(2*n,2*n)
double precision qq,pq,qp,pp
!
do kp=2,2*n,2
kq = kp-1
do lp=2,kp-2,2
lq = lp-1
qq = 0.d0
pq = 0.d0
qp = 0.d0
pp = 0.d0
do jp=2,2*n,2
jq = jp-1
qq = qq + m(lq,jq)*m(kq,jp) - m(lq,jp)*m(kq,jq)
pq = pq + m(lp,jq)*m(kq,jp) - m(lp,jp)*m(kq,jq)
qp = qp + m(lq,jq)*m(kp,jp) - m(lq,jp)*m(kp,jq)
pp = pp + m(lp,jq)*m(kp,jp) - m(lp,jp)*m(kp,jq)
enddo
! write(6,*) qq,pq,qp,pp
do i=1,2*n
m(kq,i) = m(kq,i) - qq*m(lp,i) + pq*m(lq,i)
m(kp,i) = m(kp,i) - qp*m(lp,i) + pp*m(lq,i)
enddo
enddo
qp = 0.d0
do jp=2,2*n,2
jq = jp-1
qp = qp + m(kq,jq)*m(kp,jp) - m(kq,jp)*m(kp,jq)
enddo
! write(6,*) qp
do i=1,2*n
m(kp,i) = m(kp,i)/qp
enddo
!
! Maybe the following is a better idea ( uses sqrt and is slower )
! sign = 1.d0
! if ( qp.lt.0.d0 ) sign = -1.d0
! OR, BETTER:
! if ( qp.le.0.d0 ) then complain
! qp = dabs(qp)
! qp = dsqrt(qp)
! do 600 i=1,2*n
! m(kq,i) = m(kq,i)/qp
! m(kp,i) = sign*m(kp,i)/qp
! 600 continue
enddo
return
end
logical*1 function averaged(f,a,flag,fave)
implicit none
integer isi,ndim,ndim2,nord,ntt
double precision avepol
! TAKES THE AVERAGE OF A FUNCTION F
! FLAG TRUE A=ONE TURN MAP
! FALSE A=A_SCRIPT
!
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer idpr
common /printing/ idpr
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer f,fave,a(*)
integer cosi,sine
logical flag
external avepol
integer a1(ndim2),a2(ndim2),xy(ndim2),hf(ndim2),ftf(ndim2)
logical*1 mapnormf
if(.not.flag) then
call etall(cosi,1)
call etall(sine,1)
call trx(f,f,a)
call ctor(f,cosi,sine)
call dacfu(cosi,avepol,fave)
call dadal(cosi,1)
call dadal(sine,1)
else
call etall(cosi,1)
call etall(sine,1)
call etall(ftf,nd2)
call etall(hf,nd2)
call etall(a2,nd2)
call etall(a1,nd2)
call etall(xy,nd2)
isi=0
nord=1
averaged = mapnormf(a,ftf,a2,a1,xy,hf,nord,isi)
nord=no
call etcct(a1,a2,xy)
call facflod(ftf,xy,a1,2,nord,1.d0,-1)
call trx(f,f,a1)
call ctor(f,cosi,sine)
call dacfu(cosi,avepol,fave)
call dadal(cosi,1)
call dadal(sine,1)
call dadal(ftf,nd2)
call dadal(hf,nd2)
call dadal(a2,nd2)
call dadal(a1,nd2)
call dadal(xy,nd2)
endif
return
end
double precision function avepol(j)
implicit none
integer i,ndim
parameter (ndim=3)
! PARAMETER (NTT=40)
! INTEGER J(NTT)
integer j(*)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
avepol=1.d0
do i=1,(nd-ndc)
if(j(2*i).ne.j(2*i-1)) then
avepol=0.d0
return
endif
enddo
return
end
logical function couplean(map1,tune,map2,oneturn)
implicit none
integer i,ndim,ndim2,no1,nord,ntt
double precision crazy,tpi
! map1 ascript a1 not there
! tune 2 or 3 tunes
! map2 ascript with a couple parameter in nv
! oneturn map created with tunes and map2
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer map1(*),oneturn(*),map2(*),ftf,hf
integer xy(ndim2),m1(ndim2),m2(ndim2),a2(ndim2),a1(ndim2)
integer cs,h
double precision killnonl,planar,psq(ndim),radsq(ndim)
double precision tune(ndim)
external killnonl,planar
logical*1 mapnorm
call etall(ftf,1)
call etall(hf,1)
call etall(a1,nd2)
call etall(a2,nd2)
call etall(m1,nd2)
call etall(m2,nd2)
call etall(xy,nd2)
call etall(cs,1)
call etall(h,1)
! map1 is an a-script, the last nv entry should be empty
! this a-script should around the fixed point to all orders
! one order is lost because I use PB-field
tpi=datan(1.d0)*8.d0
do i=1,nd2
call dacfu(map1(i),killnonl,m1(i))
enddo
call etini(xy)
call daclr(cs)
do i=1,nd-ndc
call dasqr(xy(2*i-1),a2(2*i-1))
call dasqr(xy(2*i),a2(2*i))
call daadd(a2(2*i-1),a2(2*i),ftf)
crazy=-tune(i)*tpi/2.d0
call dacmu(ftf,crazy,ftf)
call daadd(ftf,cs,cs)
enddo
call etinv(m1,m2)
call trx(cs,h,m2)
call dacfu(h,planar,cs)
call dasub(h,cs,h)
call davar(a2(1),0.d0,nv)
call damul(a2(1),h,h)
call daadd(cs,h,h)
call expnd2(h,xy,xy,1.d-9,1000)
call dacopd(xy,oneturn)
nord=1
couplean = mapnorm(xy,ftf,a2,a1,m2,hf,nord)
call gettura(psq,radsq)
write(6,*) (psq(i),i=1,nd)
call etini(xy)
no1=no
call fexpo(ftf,xy,xy,3,no1,1.d0,-1)
call etcct(a2,xy,map2)
call etcct(a1,map2,map2)
call dadal(ftf,1)
call dadal(hf,1)
call dadal(a1,nd2)
call dadal(a2,nd2)
call dadal(m1,nd2)
call dadal(m2,nd2)
call dadal(xy,nd2)
call dadal(cs,1)
call dadal(h,1)
return
end
double precision function planar(j)
implicit none
integer i,ndim
parameter (ndim=3)
! PARAMETER (NTT=40)
! INTEGER J(NTT)
integer j(*)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
planar=0.d0
do i=1,(nd-ndc)
if(j(2*i).eq.j(2*i-1)) then
planar=1.d0
return
endif
if(j(2*i).eq.2) then
planar=1.d0
return
endif
if(j(2*i-1).eq.2) then
planar=1.d0
return
endif
enddo
return
end
double precision function killnonl(j)
implicit none
integer i,ic,ndim
parameter (ndim=3)
! PARAMETER (NTT=40)
! INTEGER J(NTT)
integer j(*)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer ndc,ndc2,ndpt,ndt
common /coast/ndc,ndc2,ndt,ndpt
killnonl=1.d0
ic=0
do i=1,nd2-ndc2
ic=ic+j(i)
enddo
if(ic.gt.1) killnonl=0.d0
if(j(nv).ne.0) killnonl=0.d0
return
end
subroutine fexpo1(h,x,w,nrmin,nrmax,sca,ifac)
implicit none
integer ifac,ndim,ndim2,nrma,nrmax,nrmi,nrmin,ntt
double precision sca
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x,w,h
integer v(ndim2)
nrmi=nrmin-1
nrma=nrmax-1
call etall(v,nd2)
call difd(h,v,-1.d0)
call facflo(v,x,w,nrmi,nrma,sca,ifac)
call dadal(v,nd2)
return
end
subroutine etcctpar(x,ix,xj,z)
implicit none
integer i,ie,ix,ndim,ndim2,ntt
double precision xj
parameter (ndim=3)
parameter (ndim2=6)
parameter (ntt=40)
dimension xj(*)
dimension ie(ntt)
integer nd,nd2,no,nv
common /ii/no,nv,nd,nd2
integer x(*),z(*)
call etallnom(ie,nv,'IE ')
do i=1,nd2
call davar(ie(i),0.d0,i)
enddo
do i=nd2+1,nv
call dacon(ie(i),xj(i-nd2))
enddo
call dacct(x,ix,ie,nv,z,ix)
call dadal(ie,nv)
return
end