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

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

import madx-5.01.00
File size: 78.1 KB
Line 
1module bbfi
2  implicit none
3  public
4  integer bbd_max
5  parameter(bbd_max=350)
6  integer :: bbd_loc(bbd_max)=0,bbd_cnt=0,bbd_flag=0,bbd_pos=0
7  double precision :: bb_kick(2,bbd_max)=0.d0
8end module bbfi
9module deltrafi
10  implicit none
11  public
12  logical :: dorad=.false.,dodamp=.false.,dorand=.false.,fastune=.false.
13  double precision :: deltax=0.d0
14end module deltrafi
15module dyntabfi
16  implicit none
17  public
18  double precision :: dynapfrac=0.d0,dktrturns=0.d0,xend=0.d0,pxend=0.d0,&
19       yend=0.d0,pyend=0.d0,tend=0.d0,ptend=0.d0,smear=0.d0,yapunov=0.d0
20end module dyntabfi
21module wmaxmin0fi
22  implicit none
23  public
24  double precision :: wxmax=0.d0,wxmin=0.d0,wymax=0.d0,wymin=0.d0,&
25       wxymax=0.d0,wxymin=0.d0
26end module wmaxmin0fi
27module tunesfi
28  implicit none
29  public
30  double precision :: x0=0.d0,y0=0.d0,tunx=0.d0,tuny=0.d0,dtune=0.d0
31end module tunesfi
32module twiss0fi
33  implicit none
34  public
35  integer align_max,fundim
36  parameter(align_max=14,fundim = 74)
37end module twiss0fi
38module twissafi
39  implicit none
40  public
41  character(48) :: table_name=' ',sectorTableName=' '
42  logical :: match_is_on=.false.
43end module twissafi
44module twisslfi
45  implicit none
46  public
47  logical ::  centre=.false.,centre_cptk=.false.,centre_bttk=.false.,first,&
48       rmatrix=.false.,sectormap=.false.,ripken=.false.
49end module twisslfi
50module twisscfi
51  use twiss0fi
52  implicit none
53  public
54  double precision :: opt_fun0(fundim)=0.d0,opt_fun(fundim)=0.d0,disp(6)=0.d0,&
55       ddisp(6)=0.d0,rmat(2,2)=0.d0,betx=0.d0,alfx=0.d0,amux=0.d0,bety=0.d0,&
56       alfy=0.d0,amuy=0.d0,bxmax=0.d0,dxmax=0.d0,bymax=0.d0,dymax=0.d0,&
57       xcomax=0.d0,ycomax=0.d0,sigxco=0.d0,sigyco=0.d0,sigdx=0.d0,sigdy=0.d0,&
58       wgt=0.d0,cosmux=0.d0,cosmuy=0.d0,wx=0.d0,phix=0.d0,dmux=0.d0,wy=0.d0,&
59       phiy=0.d0,dmuy=0.d0,synch_1=0.d0,synch_2=0.d0,synch_3=0.d0,synch_4=0.d0,&
60       synch_5=0.d0,suml=0.d0,circ=0.d0,eta=0.d0,alfa=0.d0,gamtr=0.d0,qx=0.d0,&
61       qy=0.d0,sinmux=0.d0,sinmuy=0.d0,xix=0.d0,xiy=0.d0,currpos=0.d0
62end module twisscfi
63module twissotmfi
64  implicit none
65  public
66  double precision :: rotm(6,6)=0.d0,rw(6,6)=0.d0,skick(6)=0.d0,sorb(6)=0.d0,&
67       srmat(6,6)=0.d0,stmat(6,6,6)=0.d0
68end module twissotmfi
69module max_iterate
70  implicit none
71  public
72  integer MAXITER
73  parameter(MAXITER=150)
74end module max_iterate
75module twiss_elpfi
76  implicit none
77  public
78  !---fixed positions for element parameters-----------------------------*
79  double precision :: g_elpar(50)=0.d0
80  !-general--------------------------------------------------------------*
81  integer    g_el, g_kmax, g_kmin, g_calib, g_polarity
82  parameter (g_el = 2, g_kmax = 3, g_kmin = 4, g_calib = 5, g_polarity = 6)
83  !-bend-----------------------------------------------------------------*
84  integer    b_angle , b_tilt , b_k0 , b_k0s ,                      &
85       b_k1 , b_k1s , b_e1 , b_e2 , b_k2 ,                    &
86       b_k2s , b_h1 , b_h2 , b_hgap ,                         &
87       b_fint , b_fintx , b_k3 , b_k3s
88  parameter (b_angle = 7, b_tilt = 8, b_k0 = 9, b_k0s = 10,         &
89       b_k1 = 11, b_k1s = 12, b_e1 = 13 , b_e2 = 14, b_k2 =15,&
90       b_k2s = 16, b_h1 = 17, b_h2 = 18, b_hgap = 19,         &
91       b_fint = 20, b_fintx = 21, b_k3 = 22, b_k3s = 23)
92  !-quad-----------------------------------------------------------------*
93  integer    q_tilt, q_k1 , q_k1s
94  parameter (q_tilt = 7, q_k1 = 8, q_k1s = 9)
95  !-sext-----------------------------------------------------------------*
96  integer    s_tilt, s_k2 , s_k2s
97  parameter (s_tilt = 7, s_k2 = 8, s_k2s = 9)
98  !-oct------------------------------------------------------------------*
99  integer    o_tilt, o_k3 , o_k3s
100  parameter (o_tilt = 7, o_k3 = 8, o_k3s = 9)
101  !-mult-----------------------------------------------------------------*
102  integer    m_tilt, m_lrad
103  parameter (m_tilt = 7, m_lrad = 8)
104  !-sol------------------------------------------------------------------*
105  integer    so_lrad, so_ks, so_ksi
106  parameter (so_lrad = 7, so_ks = 8, so_ksi = 9)
107  !-rfc------------------------------------------------------------------*
108  integer    r_volt, r_lag , r_freq
109  parameter (r_volt = 7, r_lag = 8, r_freq = 9)
110  !-elsep----------------------------------------------------------------*
111  integer    e_tilt, e_ex , e_ey
112  parameter (e_tilt = 7, e_ex = 8, e_ey = 9)
113  !-hkick----------------------------------------------------------------*
114  integer    h_tilt, h_lrad , h_kick, h_hkick, h_chkick
115  parameter (h_tilt = 7, h_lrad = 8, h_kick = 9, h_hkick = 10,       &
116       h_chkick = 11)
117  !-vkick----------------------------------------------------------------*
118  integer    v_tilt, v_lrad , v_kick, v_vkick, v_cvkick
119  parameter (v_tilt = 7, v_lrad = 8, v_kick = 9, v_vkick = 10,       &
120       v_cvkick = 11)
121  !-kick-----------------------------------------------------------------*
122  integer    k_tilt, k_lrad , k_hkick, k_vkick, k_chkick, k_cvkick
123  parameter (k_tilt = 7, k_lrad = 8, k_hkick = 9, k_vkick = 10,      &
124       k_chkick = 11, k_cvkick = 12)
125end module twiss_elpfi
126module emitfi
127  implicit none
128  public
129  double precision :: qx=0.d0,qy=0.d0,qs=0.d0,cg=0.d0,sum(3)=0.d0,sumu0=0.d0
130  save qx, qy, qs, cg,sum,sumu0
131end module emitfi
132module twtrrfi
133  implicit none
134  public
135  !---- maxmul is the maximum multipole order both in twiss and trrun
136  integer maxmul,maxferr,maxnaper
137  parameter(maxmul=20,maxferr=50,maxnaper=100)
138end module twtrrfi
139module ibsdbfi
140  implicit none
141  public
142  integer :: bunch=0
143  double precision :: circ=0.d0,clight=0.d0,arad=0.d0,freq0=0.d0,alpha=0.d0,&
144       amass=0.d0,charge=0.d0,en0=0.d0,gammas=0.d0,gamma=0.d0,ex=0.d0,ey=0.d0,&
145       et=0.d0,sigt=0.d0,sige=0.d0,betas=0.d0,beta=0.d0,parnum=0.d0,&
146       currnt=0.d0,sigx=0.d0,sigy=0.d0,alfa=0.d0
147end module ibsdbfi
148module matchfi
149  implicit none
150  public
151  integer :: icovar=0,ilevel=0
152  double precision :: edm=0.d0,fmin=0.d0
153end module matchfi
154module name_lenfi
155  implicit none
156  public
157  integer name_len
158  parameter(name_len=48)
159end module name_lenfi
160module physconsfi
161  implicit none
162  public
163  double precision :: amu0=0.d0,elamda=0.d0,emass=0.d0,eps0=0.d0,erad=0.d0,&
164       hbar=0.d0,plamda=0.d0,pmass=0.d0,prad=0.d0,qelect=0.d0,mumass=0.d0
165end module physconsfi
166module touschekfi
167  implicit none
168  public
169  integer :: bunch=0
170  double precision :: circ=0.d0,clight=0.d0,arad=0.d0,freq0=0.d0,amass=0.d0,&
171       charge=0.d0,en0=0.d0,gammas=0.d0,gamma=0.d0,ex=0.d0,ey=0.d0,et=0.d0,&
172       sigt=0.d0,sige=0.d0,betas=0.d0,beta=0.d0,parnum=0.d0,currnt=0.d0,&
173       alfa=0.d0,um1=0.d0,deltap=0.d0,fb1=0.d0,fb2=0.d0
174end module touschekfi
175module trackfi
176  implicit none
177  public
178  double precision :: arad=0.d0,betas=0.d0,beti=0.d0,gammas=0.d0,dtbyds=0.d0,&
179       deltas=0.d0,bet0=0.d0,bet0i=0.d0
180  logical :: dodamp=.false.,dorad=.false.,dorand=.false.,fsecarb=.false.
181end module trackfi
182module plotfi
183  implicit none
184  public
185  !--- m_adble is the number of different types of elements
186  integer mtype, m_adble
187  parameter (mtype = 50, m_adble = 20)
188
189  !--- mcnam adjusted to NAME_L
190  integer mcnam, maxpnt
191  parameter (mcnam = 48, maxpnt = 500)
192
193  !--- szcompar is the size of the arrays returned by the routine comm_para
194  integer szcompar
195  parameter (szcompar = 100)
196
197  !--- szchara is the size of the character strings char_a & version in
198  !--- the routine pesopt
199  integer szchara
200  parameter (szchara = 400)
201
202  !--- character sizes:
203  !   MLSIZE    label character height
204  !   MTSIZE    text   - " -
205  !   MASIZE    annotation - " -
206  integer mlsize, mtsize, masize
207  parameter  (mlsize = 13,mtsize = 13,masize = 20)
208
209  !--- parameters used in the routine pegacn in file plot.F
210
211  integer mposx, mposy, mpost
212  parameter (mposx = 8, mposy = 3, mpost = mposx * mposy)
213
214  !--- parameters used in the routine peschm in file plot.F
215
216  integer mobj, msize
217  parameter (mobj = 14, msize = 88)
218
219  integer mtitl, mxlabl, mnvar, mxdep, mqadd, mntmax, mksmax, mplred, mplout
220
221  parameter (mtitl  = 128, mxlabl = 160)
222  parameter (mnvar = 74, mxdep = 2)
223  parameter (mqadd = 100000)
224  parameter (mntmax = 20, mksmax = 10)
225  parameter (mplred = 46, mplout = 47)
226
227  integer maxseql, mtwcol, mpparm, mxcurv, mopt, mfile, marg, maxarg, mxdp, mxplot
228
229  parameter (maxseql = 20000, mtwcol = 46, mpparm = 10,             &
230       mxcurv = 10, mopt = 60, mfile = 120, marg = 60, maxarg = 1000,    &
231       mxdp = 25, mxplot = 100)
232
233  integer mintpl
234  parameter (mintpl = 18)
235
236  !--- Definition of common / peaddi /
237  !--- itbv set in routine pesopt, used in routines pemima, peplot
238  !---      and pesopt.
239  !--- nivvar set in routine pesopt, used in routines pefill, peintp,
240  !---        pemima, peplot and pesopt.
241  !--- nelmach set in routine pefill, used in routines pefill, peplot.
242  !--- numax set in routine pemima, used in routines pemima, peplot.
243  !--- interf set in routine pesopt, used in routines pecurv, pefill.
244  !--- noline set in routine pesopt, used in routines pefill, pesopt.
245  !--- nqval set in routine pefill and peintp, used in routines pefill,
246  !---        peintp, pemima and peplot.
247  !--- nvvar set in routine pesopt and pemima, used in routine pemima.
248  !--- nrrang set in routine pefill, used in routine pesopt and pefill.
249  !--- proc_flag set in routine pesopt, used in routine pefill, peintp
250  !---           and pesopt.
251  !--- ipparm set in routine peintp and pesopt, used in routine peplot
252  !---           and pesopt.
253  !--- naxref set in routine pemima and pesopt, used in routine pemima
254  !---           and pesopt.
255  !--- ieltyp set in routine pefill, used in routine psplot.
256
257  integer :: itbv=0,nivvar=0,nelmach=0,numax=0,interf=0,noline=0, &
258       nqval(mxcurv)=0,nvvar(4)=0,nrrang(2)=0,                         &
259       proc_flag(2,mxcurv)=0,ipparm(mpparm,mxcurv)=0,                  &
260       naxref(mxcurv)=0,ieltyp(maxseql)=0
261
262  !--- Definition of common / peaddr /
263  !--- qascl set in routine pesopt, used in routine peplot.
264  !--- qlscl set in routine pesopt, used in routine peplot.
265  !--- qsscl set in routine pesopt, used in routine peplot.
266  !--- qtscl set in routine pesopt, used in routine peplot.
267  !--- hrange set in routine pesopt, used in routines pefill, peplot.
268  !--- vrange set in routine pesopt, used in routine peplot.
269  !--- hmima set in routines pesopt and pemima,
270  !---       used in routines peplot, pesopt and pemima.
271  !--- vmima set in routines pesopt and pemima,
272  !---       used in routines peplot and pemima.
273  !--- qhval set in routines pefill and peintp,
274  !---       used in routines peplot, pefill and pemima.
275  !--- qvval set in routines pefill and peintp,
276  !---       used in routines peplot, pefill and pemima.
277  !--- estart set in routine pefill, and peintp,
278  !---        used in routines peplot and pefill.
279  !--- eend set in routine pefill, used in routine peplot.
280
281  real :: qascl=0.,qlscl=0.,qsscl=0.,qtscl=0.,                    &
282       hrange(2)=0.,vrange(2,4)=0.,hmima(2)=0.,vmima(2,4)=0.,          &
283       qhval(maxseql,mxcurv)=0.,qvval(maxseql,mxcurv)=0.,              &
284       estart(maxseql)=0.,eend(maxseql)=0.
285
286  !--- Definition of common / peaddc /
287  !--- horname set in routine pesopt,
288  !---         used in routines pefill, peplot and pesopt.
289  !--- tabname set in routine pesopt,
290  !---         used in routines pefill, peintp, pelfill and pesopt.
291  !--- toptitle set in routine pesopt, used in routine peplot.
292  !--- plfnam set in routine plginit, used in routines plotit and plginit.
293  !--- axlabel set in routine pemima, used in routine peplot.
294  !--- sname set in routine pesopt,
295  !---         used in routines pefill and pesopt.
296  !--- slabl set in routine pesplit,
297  !---         used in routines peplot, pemima and pesopt.
298
299  character(mfile) :: plfnam=' '
300  character(mcnam) :: horname=' ',tabname=' ',sname(mxcurv)=' ',slabl(mxcurv)=' '
301  character(mxlabl) :: axlabel(4)=' '
302  character(mtitl) :: toptitle=' '
303
304  save itbv,nivvar,nelmach,numax,interf,noline,nqval,nvvar,nrrang,proc_flag,&
305       ipparm,naxref,ieltyp,qascl,qlscl,qsscl,qtscl,hrange,vrange,hmima,&
306       vmima,qhval,qvval,estart,eend,horname,tabname,toptitle,plfnam,axlabel,&
307       sname,slabl
308end module plotfi
309module plot_bfi
310  implicit none
311  public
312  !--- Definition of common / peotcl /
313  !--- fpmach set in routines pesopt and pefill, used in routine peplot
314  !--- ddp_flag set in routine pefill, used in routine peplot
315  !--- ptc_flag set in routines pesopt, used in routine pefill
316
317  logical :: fpmach=.false.,dpp_flag=.false.,ptc_flag=.false.
318  save fpmach,dpp_flag,ptc_flag
319end module plot_bfi
320module plot_cfi
321  implicit none
322  public
323  !--- Definition of common / e2save /
324  !--- e2s initialised in routine pefill, used in routine peelma
325
326  double precision :: e2s=0.d0
327end module plot_cfi
328module plot_mathfi
329  implicit none
330  public
331  !--- Definitions of mathematical constants
332
333  double precision pi, zero, eps, one, two, twopi, half
334  parameter         (pi = 3.1415926535898d0)
335  parameter         (zero = 0.d0, half = 0.5d0, eps = 1.d-5)
336  parameter         (one = 1.d0, two = 2.d0, twopi = two * pi)
337end module plot_mathfi
338module resindexfi
339  implicit none
340  public
341  integer mnres,mymorder
342  parameter (mnres=1000,mymorder=20)
343end module resindexfi
344module gxx11_common
345  implicit none
346  public
347  !
348  integer madim1,madim2,maxset,mconid,merrun,metaun,miunit,mmetat,&
349       normt,mounit,mpaxs,mpcurv,mtermt,mtick,mtmeta,mtterm,mxaxs,mxpix,&
350       mxsize,myaxs,mypix,mysize,mnormt,mx11pr,mx11tf,mxxpix,mxypix,&
351       mcolor,mpspro,meppro,mdict,mlpro,&
352       mpsep,mepep,mhead,mline,msfact,mlbb1,mlbb2,mubb1,mubb2,mtfont,&
353       mwid1,mwid2
354  real toleps,versio
355  parameter (mxaxs = 4, myaxs = 4, mpaxs = 23, mpcurv = 10,&
356       maxset = 20, mtterm = 1, mmetat = 4,&
357       mtermt = 101, mtmeta = 2, mconid = 7, mtick = 10, metaun = 11,&
358       mxpix = 1000, mypix = 1000, mxsize = 27, mysize = 19,&
359       madim1 = 500, toleps = 1.e-5,&
360       merrun = 10, miunit = 5, mounit = 6, versio = 1.50,&
361       mx11pr = 10, mx11tf = 10, mxxpix = 1200, mxypix = 1000,&
362       mcolor = 6, mpspro = 8, meppro = 8, mdict = 24, mlpro = 68,&
363       mpsep = 3, mepep = 2, mhead = 4, mline = 72, msfact = 4,&
364       mlbb1 = 17, mlbb2 = 30, mubb1 = 573, mubb2 = 790, mtfont = 12,&
365       mwid1 = mubb1 - mlbb1, mwid2 = mubb2 - mlbb2 )
366  parameter (mnormt = 2, madim2 = 100)
367  !
368  integer :: &
369       itermt=0, interm=0, inmeta=0, ierrun=0, imetun=0, inunit=0, iounit=0, ipage=0,&
370       isfflg=0, isqflg=0, iwtflg=0, iclflg=0, inormt=0, ipseps=0, iepsop=0, itseop=0,&
371       iepscf=0, imetps=0, ipctct=0, iczebr=0, idinit=0, ipstyp=0, iclear=0, istotx=0,&
372       lmpict=0, ltermt=0, lnterm=0, lnmeta=0, lerrun=0, lmetun=0, lnunit=0, lounit=0,&
373       lsfflg=0, lsqflg=0, lwtflg=0, lclflg=0, lnormt=0, lmetax=0, lmetay=0, lmetnm=0,&
374       lerrnm=0, ldefnl=0, lerrop=0, lmetop=0, ltotin=0, lacttm=0, lpseps=0, lundef=0,&
375       lttime=0, ldinit=0, ltseop=0,&
376       ixapar(mpaxs,mxaxs)=0, iyapar(mpaxs,myaxs)=0, icvpar(mpcurv ,maxset)=0
377  !
378  integer :: &
379       nxpix=0, nypix=0, lxpix=0, lypix=0, icucol=0, iorips=0,      &
380       iutlps=0, ibbox(4)=0, ix11pr(mx11pr)=0, ix11tf(mx11tf)=0, ix11op(mx11tf)=0  !
381  !
382  real :: &
383       fxpix=0., fypix=0., rx11pr(mx11pr)=0., rgbcol(3,mcolor)=0.
384  !
385  real :: &
386       xmetaf=0., ymetaf=0., xsterm=0., ysterm=0., wfact=0., wttime=0., wxfact=0., wyfact=0.,&
387       vpfacx=0., vpfacy=0.,&
388       vptdef(4)=0., vploc(4)=0., actwnd(4)=0., rangex(2,mxaxs)=0., rangey(2 ,myaxs)=0.,&
389       cvwnwd(4,maxset)=0., axwndx(2,maxset), axwndy(2,maxset)=0.
390  !
391  character :: &
392       smetnm*256=" ", serrnm*256=" ", sxtext(mxaxs)*300=" ", sytext(myaxs)*300=" ",&
393       sxform(mxaxs)*20=" ", syform(myaxs)*20=" ", splotc*(maxset)=" ", stortx * 20=" ",&
394       sdefnl*1=" ",spsnam * 256=" ", colour(mcolor) * 16=" ", pshead(mhead) * 60=" "
395  !
396  real :: xp(madim2+1)=0.,xvp(madim2+1)=0,yp(madim2+1)=0.,yvp(madim2+1)=0
397  !
398  real :: p(madim1,2)=0.,s(madim1)=0.,yy1d(madim1,2)=0.,yy2d(madim1,2)=0.
399end module gxx11_common
400module gxx11_aux
401  implicit none
402  public
403  !
404  character(100) strloc
405  integer, dimension(14) :: ivals=(/ 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 1, 1, 0 /)
406  !   ivals(1)       marker type
407  !     (2)       fill area interior style
408  !     (3)       horizontal text alignment
409  !     (4)       vertical text alignment
410  !     (5)       text font
411  !     (6)       text precision
412  !     (7)       marker colour index
413  !     (8)       metafile status (0 closed, 1 open)
414  !     (9)       text colour index
415  !    (10)       free
416  !    (11)       polyline colour index
417  !    (12)       polyline style
418  !    (13)       current normalisation transformation number
419  !    (14)       last call type: 0 undef., 1 line, 2 text, 3 marker
420
421  real, dimension(14) :: rvals=(/ 0., 1., 0.01, 0., 1., 0., 1., 0., 1., 0., 1., 1., 1., 1. /)
422  !   rvals(1-2)  chup vector
423  !     3         character height
424  !     4-7       window
425  !     8-11      viewport
426  !     12        character expansion factor
427  !     13        line width scale factor
428  !     14        marker scale factor
429  save    ivals, rvals
430end module gxx11_aux
431subroutine fort_info(t1, t2)
432  implicit none
433
434
435  character(*) t1, t2
436  integer get_option
437  if (get_option('info ') .ne. 0 .and. get_option('warn ') .ne. 0)  &
438       print '(a,1x,a,1x,a)', '++++++ info:', t1, t2
439end subroutine fort_info
440subroutine fort_warn(t1, t2)
441  implicit none
442
443
444
445  character(*) t1, t2
446  integer get_option
447  if (get_option('warn ') .ne. 0) then
448     print '(a,1x,a,1x,a)', '++++++ warning:', t1, t2
449     call augmentfwarn()
450  endif
451end subroutine fort_warn
452subroutine getclor(orbit0, rt, tt, error)
453
454  !----------------------------------------------------------------------*
455  ! Purpose:
456  !   Get periodic closed orbit (e.g. at start of Twiss),
457  !   first + second order one-turn map
458  ! Input:
459  !   orbit0(6)   (real)  initial guess
460  ! Output:
461  !   rt(6,6)     (real)  one-turn matrix
462  !   tt(6,6,6)   (real)  one-turn second-order map
463  !   error       (int)   error flag (0: OK, else != 0)
464  !----------------------------------------------------------------------*
465  use twiss0fi
466  implicit none
467
468  double precision orbit0(6), rt(6,6), tt(6,6,6)
469  double precision opt(fundim)
470  integer error
471  call m66one(rt)
472  call dzero(opt,fundim)
473  call tmclor(orbit0, .true., .true., opt, rt, tt, error)
474end subroutine getclor
475subroutine m66add(term1,term2,target)
476  implicit none
477  !----------------------------------------------------------------------*
478  ! Purpose:
479  !   Add two matrices.
480  ! Input:
481  !   TERM1(6,6)  (real)  First term.
482  !   TERM2(6,6)  (real)  Second term.
483  ! Output:
484  !   TARGET(6,6) (real)  Sum: TARGET = TERM1 + TERM2.
485  !----------------------------------------------------------------------*
486  integer i,j
487  double precision target(6,6),term1(6,6),term2(6,6)
488
489  do i = 1, 6
490     do j = 1, 6
491        target(i,j) = term1(i,j) + term2(i,j)
492     enddo
493  enddo
494
495end subroutine m66add
496subroutine m66byv(amat,avec,target)
497  implicit none
498  !----------------------------------------------------------------------*
499  ! Purpose:
500  !   Multiply matrix times vector.
501  ! Input:
502  !   AMAT(6,6)   (real)  Input matrix.
503  !   AVEC(6)     (real)  Input vector.
504  ! Output:
505  !   TARGET(6)   (real)  Output vector: TARGET = AMAT * AVEC.
506  !----------------------------------------------------------------------*
507  integer i,j
508  double precision amat(6,6),avec(6),target(6),temp(6)
509
510  call dzero(temp,6)
511  do i = 1, 6
512     do j = 1, 6
513        temp(i) = temp(i) + amat(i,j) * avec(j)
514     enddo
515  enddo
516
517  call dcopy(temp,target,6)
518
519end subroutine m66byv
520subroutine m66cpy(source,target)
521  implicit none
522  !----------------------------------------------------------------------*
523  ! Purpose:
524  !   Copy matrix.
525  ! Input:
526  !   SOURCE(6,6) (real)  Input matrix.
527  ! Output:
528  !   TARGET(6,6) (real)  Output matrix: TARGET = SOURCE.
529  !----------------------------------------------------------------------*
530  integer i,j
531  double precision source(6,6),target(6,6)
532
533  do i = 1, 6
534     do j = 1, 6
535        target(i,j) = source(i,j)
536     enddo
537  enddo
538
539end subroutine m66cpy
540subroutine m66div(anum,aden,target,eflag)
541  implicit none
542  !----------------------------------------------------------------------*
543  ! Purpose:
544  !   "Divide" matrices, i. e. postmultiply with inverse of denominator.
545  ! Input:
546  !   ANUM(6,6)   (real)  "Numerator" matrix.
547  !   ADEN(6,6)   (real)  "Denominator" matrix.
548  ! Output:
549  !   TARGET(6,6) (real)  "Quotient" matrix: TARGET = ANUM * ADEN**(-1).
550  !   EFLAG       (logical) Error flag.
551  !----------------------------------------------------------------------*
552  logical eflag
553  integer i,irank,j
554  double precision aden(6,6),anum(6,6),augmat(6,12),target(6,6)
555
556  !---- Copy input to local array.
557  do i = 1, 6
558     do j = 1, 6
559        augmat(i,j)   = aden(i,j)
560        augmat(i,j+6) = anum(i,j)
561     enddo
562  enddo
563
564  !---- Solve resulting system.
565  call solver(augmat,6,6,irank)
566  if (irank .lt. 6) then
567     eflag = .true.
568
569     !---- Copy result.
570  else
571     eflag = .false.
572     do i = 1, 6
573        do j = 1, 6
574           target(i,j) = augmat(i,j+6)
575        enddo
576     enddo
577  endif
578
579end subroutine m66div
580subroutine m66exp(source,target,eflag)
581  implicit none
582  !----------------------------------------------------------------------*
583  ! Purpose:
584  !   "Exponentiate" matrix.
585  !   Original author:    Liam Healy.
586  ! Input:
587  !   SOURCE(6,6) (real)  Input matrix.
588  ! Output:
589  !   TARGET(6,6) (real)  Output matrix: TARGET = exp(SOURCE).
590  !   EFLAG     (logical) Error flag.
591  !----------------------------------------------------------------------*
592  logical eflag
593  integer i,j
594  double precision b(6,6),c(6,6),source(6,6),target(6,6),one,two,twelve
595  parameter(one=1d0,two=2d0,twelve=12d0)
596
597  call m66mpy(source,source,b)
598  call m66mpy(source,b,c)
599  do j = 1, 6
600     do i = 1, 6
601        b(i,j) = (source(i,j) - c(i,j) / twelve) / two
602        c(i,j) = - b(i,j)
603     enddo
604     b(j,j) = b(j,j) + one
605     c(j,j) = c(j,j) + one
606  enddo
607  call m66div(b,c,target,eflag)
608
609end subroutine m66exp
610subroutine m66inv(source,target)
611  implicit none
612  !----------------------------------------------------------------------*
613  ! Purpose:
614  !   Invert symplectic matrix.
615  ! Input:
616  !   SOURCE(6,6) (real)  Input matrix.
617  ! Output:
618  !   TARGET(6,6) (real)  Output matrix: TARGET = tr(J) * tr(SOURCE) * J.
619  !----------------------------------------------------------------------*
620  integer i
621  double precision source(6,6),target(6,6),temp(6,6)
622
623  !---- TEMP = transpose(SOURCE) * J.
624  do i = 1, 6
625     temp(i,1) = - source(2,i)
626     temp(i,2) = + source(1,i)
627     temp(i,3) = - source(4,i)
628     temp(i,4) = + source(3,i)
629     temp(i,5) = - source(6,i)
630     temp(i,6) = + source(5,i)
631  enddo
632
633  !---- TARGET = transpose(J) * TEMP.
634  do i = 1, 6
635     target(1,i) = - temp(2,i)
636     target(2,i) = + temp(1,i)
637     target(3,i) = - temp(4,i)
638     target(4,i) = + temp(3,i)
639     target(5,i) = - temp(6,i)
640     target(6,i) = + temp(5,i)
641  enddo
642
643end subroutine m66inv
644subroutine m66symp(r,nrm)
645  implicit none
646  !----------------------------------------------------------------------*
647  ! Purpose:
648  !   Check if a 6 by 6 matrix R is symplectic.
649  ! Input:
650  !   r(6,6)    (double)  Matrix R to check
651  ! Output:
652  !   nrm       (double)  The column norm of R'*J*R-J
653  !----------------------------------------------------------------------*
654  double precision R(6,6),J(6,6),T(6,6),nrm,z,o,n
655  parameter(z=0d0,o=1d0,n=-1d0)
656  J = reshape((/ z, o, z, z, z, z, &
657               & n, z, z, z, z, z, &
658               & z, z, z, o, z, z, &
659               & z, z, n, z, z, z, &
660               & z, z, z, z, z, o, &
661               & z, z, z, z, n, z /), shape(J))
662  call m66trm(R,J,T)
663  call m66mpy(T,R,T)
664  call m66sub(T,J,T)
665  call m66nrm(T,nrm)
666end subroutine m66symp
667subroutine m66mak(f2,target)
668  implicit none
669  !----------------------------------------------------------------------*
670  ! Purpose:
671  !   Compute matrix TARGET corresponding to Lie polynomial F2.
672  !   Original author:    Liam Healy.
673  ! Input:
674  !   F2          (poly)  Polynomial of order 2.
675  ! Output:
676  !   TARGET(6,6) (real)  Output matrix: TARGET * v = - [J,v].
677  !----------------------------------------------------------------------*
678  double precision f2(*),target(6,6),two
679  parameter(two=2d0)
680
681  target(1,1) = - f2(8)
682  target(1,2) = - two * f2(13)
683  target(1,3) = - f2(14)
684  target(1,4) = - f2(15)
685  target(1,5) = - f2(16)
686  target(1,6) = - f2(17)
687  target(2,1) = two * f2(7)
688  target(2,2) = f2(8)
689  target(2,3) = f2(9)
690  target(2,4) = f2(10)
691  target(2,5) = f2(11)
692  target(2,6) = f2(12)
693  target(3,1) = - f2(10)
694  target(3,2) = - f2(15)
695  target(3,3) = - f2(19)
696  target(3,4) = - two * f2(22)
697  target(3,5) = - f2(23)
698  target(3,6) = - f2(24)
699  target(4,1) = f2(9)
700  target(4,2) = f2(14)
701  target(4,3) = two * f2(18)
702  target(4,4) = f2(19)
703  target(4,5) = f2(20)
704  target(4,6) = f2(21)
705  target(5,1) = - f2(12)
706  target(5,2) = - f2(17)
707  target(5,3) = - f2(21)
708  target(5,4) = - f2(24)
709  target(5,5) = - f2(26)
710  target(5,6) = - two * f2(27)
711  target(6,1) = f2(11)
712  target(6,2) = f2(16)
713  target(6,3) = f2(20)
714  target(6,4) = f2(23)
715  target(6,5) = two * f2(25)
716  target(6,6) = f2(26)
717
718end subroutine m66mak
719subroutine m66mpy(fact1,fact2,target)
720  implicit none
721  !----------------------------------------------------------------------*
722  ! Purpose:
723  !   Multiply two matrices.
724  !   TARGET may coincide with one of the factors.
725  ! Input:
726  !   FACT1(6,6)  (real)  First factor.
727  !   FACT2(6,6)  (real)  Second factor.
728  ! Output:
729  !   TARGET(6,6) (real)  Product matrix: TARGET = FACT1 * FACT2.
730  !----------------------------------------------------------------------*
731  integer i,j,k
732  double precision fact1(6,6),fact2(6,6),target(6,6),temp(6,6)
733
734  call dzero(temp,36)
735  do k = 1, 6
736     do j = 1, 6
737        do i = 1, 6
738           temp(i,k) = temp(i,k) + fact1(i,j) * fact2(j,k)
739        enddo
740     enddo
741  enddo
742  call dcopy(temp,target,36)
743
744end subroutine m66mpy
745subroutine m66mtr(fact1,fact2,target)
746  implicit none
747  !----------------------------------------------------------------------*
748  ! Purpose:
749  !   Multiply a matrix with the transpose of another matrix.
750  !   TARGET must not coincide with either factor.
751  ! Input:
752  !   FACT1(6,6)  (real)  First factor.
753  !   FACT2(6,6)  (real)  Second factor (will be transposed).
754  ! Output:
755  !   TARGET(6,6) (real)  Product: TARGET = FACT1 * tr(FACT2).
756  !----------------------------------------------------------------------*
757  integer i,j,k
758  double precision fact1(6,6),fact2(6,6),target(6,6)
759
760  call dzero(target,36)
761  do j = 1, 6
762     do k = 1, 6
763        do i = 1, 6
764           target(i,j) = target(i,j) + fact1(i,k) * fact2(j,k)
765        enddo
766     enddo
767  enddo
768
769end subroutine m66mtr
770subroutine m66nrm(fm,res)
771  implicit none
772  !----------------------------------------------------------------------*
773  ! Purpose:
774  !   Computes the norm of a matrix.
775  !   Reference:          L. Collatz,
776  !                       Functional Analysis & Numerical Mathematics.
777  !   Source:             MARYLIE, Version 3.0.
778  ! Input:
779  !   FM(6,6)     (real)  Input matrix.
780  ! Output:
781  !   RES         (real)  Norm of FM: RES = max abs column sum.
782  !----------------------------------------------------------------------*
783  integer i,j
784  double precision fm(6,6),res,sum,zero
785  parameter(zero=0d0)
786
787  res = zero
788  do j = 1, 6
789     sum = zero
790     do i = 1, 6
791        sum = sum + abs(fm(i,j))
792     enddo
793     res = max(res,sum)
794  enddo
795
796end subroutine m66nrm
797subroutine m66one(target)
798  implicit none
799  !----------------------------------------------------------------------*
800  ! Purpose:
801  !   Set matrix to unity.
802  ! Output:
803  !   TARGET(6,6) (real)  Unit matrix: TARGET = I.
804  !----------------------------------------------------------------------*
805  integer j
806  double precision target(6,6),one
807  parameter(one=1d0)
808
809  call dzero(target,36)
810  do j = 1, 6
811     target(j,j) = one
812  enddo
813
814end subroutine m66one
815subroutine m66ref(source,target)
816  implicit none
817  !----------------------------------------------------------------------*
818  ! Purpose:
819  !   Reflect symplectic first order transform.
820  ! Input:
821  !   SOURCE(6,6) (real)  Input matrix.
822  ! Output:
823  !   TARGET(6,6) (real)  Reflected matrix.
824  !----------------------------------------------------------------------*
825  integer i
826  double precision source(6,6),target(6,6),temp(6,6)
827
828  !---- TEMP = transpose(SOURCE) * J * signs.
829  do i = 1, 6
830     temp(i,1) =   source(2,i)
831     temp(i,2) =   source(1,i)
832     temp(i,3) =   source(4,i)
833     temp(i,4) =   source(3,i)
834     temp(i,5) = - source(6,i)
835     temp(i,6) = - source(5,i)
836  enddo
837
838  !---- TARGET = signs * transpose(J) * TEMP.
839  do i = 1, 6
840     target(1,i) =   temp(2,i)
841     target(2,i) =   temp(1,i)
842     target(3,i) =   temp(4,i)
843     target(4,i) =   temp(3,i)
844     target(5,i) = - temp(6,i)
845     target(6,i) = - temp(5,i)
846  enddo
847
848end subroutine m66ref
849subroutine m66scl(scalar,source,target)
850  implicit none
851  !----------------------------------------------------------------------*
852  ! Purpose:
853  !   Multiply matrix by scalar.
854  ! Input:
855  !   SCALAR      (real)  Scale factor.
856  !   SOURCE(6,6) (real)  Input matrix.
857  ! Output:
858  !   TARGET(6,6) (real)  Scaled matrix: TARGET = SCALAR * SOURCE.
859  !----------------------------------------------------------------------*
860  integer i,j
861  double precision scalar,source(6,6),target(6,6)
862
863  do i = 1, 6
864     do j = 1, 6
865        target(i,j) = scalar * source(i,j)
866     enddo
867  enddo
868
869end subroutine m66scl
870logical function m66sta(amat)
871  implicit none
872  !----------------------------------------------------------------------*
873  ! Purpose:
874  !   Check effect of a matrix on momentum.
875  ! Input:
876  !   AMAT(6,6)   (real)  Input matrix.
877  ! Result:
878  !   .TRUE.              For static case     (constant p).
879  !   .FALSE.             For dynamic case    (variable p).
880  !----------------------------------------------------------------------*
881  integer j
882  double precision amat(6,6),tol,one
883  parameter(one=1d0,tol=1d-12)
884
885  m66sta = abs(amat(6,6) - one) .le. tol
886  do j = 1, 5
887     m66sta = m66sta .and. abs(amat(6,j)) .le. tol
888  enddo
889
890end function m66sta
891subroutine m66sub(term1,term2,target)
892  implicit none
893  !----------------------------------------------------------------------*
894  ! Purpose:
895  !   Subtract two matrices.
896  ! Input:
897  !   TERM1(6,6)  (real)  Minuend matrix.
898  !   TERM2(6,6)  (real)  Subtrahend matrix.
899  ! Output:
900  !   TARGET(6,6) (real)  Difference matrix: TARGET = TERM1 - TERM2.
901  !----------------------------------------------------------------------*
902  integer i,j
903  double precision target(6,6),term1(6,6),term2(6,6)
904
905  do j = 1, 6
906     do i = 1, 6
907        target(i,j) = term1(i,j) - term2(i,j)
908     enddo
909  enddo
910
911end subroutine m66sub
912subroutine m66trm(fact1,fact2,target)
913  implicit none
914  !----------------------------------------------------------------------*
915  ! Purpose:
916  !   Multiply the transpose of a matrix with another matrix.
917  !   TARGET must not coincide with either factor.
918  ! Input:
919  !   FACT1(6,6)  (real)  First factor (will be transposed).
920  !   FACT2(6,6)  (real)  Second factor.
921  ! Output:
922  !   TARGET(6,6) (real)  Product: TARGET = tr(FACT1) * FACT2.
923  !----------------------------------------------------------------------*
924  integer i,j,k
925  double precision fact1(6,6),fact2(6,6),target(6,6)
926
927  call dzero(target,36)
928  do j = 1, 6
929     do k = 1, 6
930        do i = 1, 6
931           target(i,j) = target(i,j) + fact1(k,i) * fact2(k,j)
932        enddo
933     enddo
934  enddo
935
936end subroutine m66trm
937subroutine m66tp(source,target)
938  implicit none
939  !----------------------------------------------------------------------*
940  ! Purpose:
941  !   Transpose a matrix.
942  !   TARGET and SOURCE may overlap.
943  ! Input:
944  !   SOURCE(6,6) (real)  Input matrix.
945  ! Output:
946  !   TARGET(6,6) (real)  Transposed matrix: TARGET = tr(SOURCE).
947  !----------------------------------------------------------------------*
948  integer i,j
949  double precision source(6,6),target(6,6),temp(6,6)
950
951  do i = 1, 6
952     do j = 1, 6
953        temp(j,i) = source(i,j)
954     enddo
955  enddo
956  call m66cpy(temp,target)
957
958end subroutine m66tp
959subroutine m66zro(target)
960  implicit none
961  !----------------------------------------------------------------------*
962  ! Purpose:
963  !   Clear a matrix to zero.
964  ! Output:
965  !   TARGET(6,6) (real)  Zero matrix: TARGET = 0.
966  !----------------------------------------------------------------------*
967  double precision target(6,6)
968
969  call dzero(target,36)
970
971end subroutine m66zro
972subroutine solver(augmat,ndim,mdim,irank)
973  implicit none
974  !----------------------------------------------------------------------*
975  ! Purpose:
976  !   Solve the linear equation  A * X = B.
977  ! Input:
978  !   AUGMAT(n,n+m)       A(n,n), augmented by B(n,m).
979  !   NDIM, MDIM          n, m.
980  ! Output:
981  !   AUGMAT(n,n+m)       Identity(n,n), augmented by X(n,m).
982  !   IRANK               Rank of A.
983  !----------------------------------------------------------------------*
984  integer ic,ip,ir,irank,it,mdim,nc,ndim,nr
985  double precision augmat(ndim,ndim+mdim),h,pivot,zero
986  parameter(zero=0d0)
987
988  nr = ndim
989  nc = ndim + mdim
990  irank = 0
991  do it = 1, nr
992     pivot = zero
993     ip = 0
994     do ir = it, nr
995        if (abs(augmat(ir,it)) .ge. abs(pivot)) then
996           pivot = augmat(ir,it)
997           ip = ir
998        endif
999     enddo
1000
1001     if (pivot .eq. zero) go to 9999
1002     irank = it
1003
1004     do ic = 1, nc
1005        augmat(ip,ic) = augmat(ip,ic) / pivot
1006     enddo
1007
1008     if (ip .ne. it) then
1009        do ic = 1, nc
1010           h = augmat(ip,ic)
1011           augmat(ip,ic) = augmat(it,ic)
1012           augmat(it,ic) = h
1013        enddo
1014     endif
1015
1016     do ir = 1, nr
1017        if (ir .ne. it) then
1018           h = augmat(ir,it)
1019           do ic = 1, nc
1020              augmat(ir,ic) = augmat(ir,ic) - h * augmat(it,ic)
1021           enddo
1022        endif
1023     enddo
1024  enddo
1025
1026  irank = ndim
1027
10289999 end subroutine solver
1029subroutine symsol(a,n,eflag,work_1,work_2,work_3)
1030  implicit none
1031  !----------------------------------------------------------------------*
1032  ! Purpose:
1033  !   Invert symmetric matrix.
1034  ! Input:
1035  !   A(*,*)    (real)    Matrix to be inverted.
1036  !   N         (integer) Actual size of A.
1037  ! Output:
1038  !   A(*,*)    (real)    Inverted matrix.
1039  !   EFLAG     (logical) Error flag.
1040  !----------------------------------------------------------------------*
1041  logical eflag
1042  integer i,j,k,n
1043  double precision a(n,n),si,work_1(n),work_2(n),work_3(n),zero,one
1044  parameter(zero=0d0,one=1d0)
1045
1046  !---- Scale upper triangle.
1047  eflag = .true.
1048  do i = 1, n
1049     si = a(i,i)
1050     if (si .le. zero) go to 100
1051     work_1(i) = one / sqrt(si)
1052  enddo
1053  do i = 1, n
1054     do j = i, n
1055        a(i,j) = a(i,j) * work_1(i) * work_1(j)
1056     enddo
1057  enddo
1058
1059  !---- Invert upper triangle.
1060  do i = 1, n
1061     if (a(i,i) .eq. zero) go to 100
1062     work_2(i) = one
1063     work_3(i) = one / a(i,i)
1064     a(i,i) = zero
1065     do j = 1, n
1066        if (j .lt. i) then
1067           work_2(j) = a(j,i)
1068           work_3(j) = work_2(j) * work_3(i)
1069           a(j,i) = zero
1070        else if (j .gt. i) then
1071           work_2(j) = a(i,j)
1072           work_3(j) = - work_2(j) * work_3(i)
1073           a(i,j) = zero
1074        endif
1075     enddo
1076     do j = 1, n
1077        do k = j, n
1078           a(j,k) = a(j,k) + work_2(j) * work_3(k)
1079        enddo
1080     enddo
1081  enddo
1082
1083  !---- Rescale upper triangle and symmetrize.
1084  do i = 1, n
1085     do j = i, n
1086        a(i,j) = a(i,j) * work_1(i) * work_1(j)
1087        a(j,i) = a(i,j)
1088     enddo
1089  enddo
1090  eflag = .false.
1091
1092100 continue
1093
1094end subroutine symsol
1095subroutine symeig(a,nd,n,eigen,nval,work)
1096  implicit none
1097  !----------------------------------------------------------------------*
1098  ! Purpose:                                                             *
1099  !   Eigenvalues of a real symmetric matrix in ascending order.         *
1100  ! Input:                                                               *
1101  !   A(ND,ND)  (real)    Symmetric input matrix; destroyed by call.     *
1102  !   N         (integer) Rank of matrix.                                *
1103  ! Output:                                                              *
1104  !   EIGEN(*)  (real)    Eigenvalues of A in descending order.          *
1105  !   NVAL      (integer) Number of eigenvalues found.                   *
1106  !----------------------------------------------------------------------*
1107  integer i,it,j,k,l,m,n,nd,nval,itmax
1108  parameter(itmax=15)
1109  double precision b,c,f,g,h,p,r,s,work(nd),a(nd,nd),eigen(nd),zero,&
1110       one,two,four,big,eps
1111  parameter(zero=0d0,one=1d0,two=2d0,four=4d0,big=1d10,eps=1d-20)
1112
1113  !---- Matrix is 1 * 1.
1114  nval = n
1115  if (n .le. 0) go to 300
1116  if (n .eq. 1) then
1117     eigen(1) = a(1,1)
1118     go to 300
1119  endif
1120
1121  !---- Matrix is 2 * 2.
1122  if (n .eq. 2) then
1123     f = a(1,1) + a(2,2)
1124     g = sqrt((a(1,1) - a(2,2))**2 + four * a(2,1)**2)
1125     eigen(1) = (f - g) / two
1126     eigen(2) = (f + g) / two
1127     go to 300
1128  endif
1129
1130  !---- N is at least 3, reduce to tridiagonal form.
1131  do i = n, 3, -1
1132     g = zero
1133     do k = 1, i-2
1134        g = g + a(i,k)**2
1135     enddo
1136     eigen(i) = a(i,i)
1137     if (g .eq. zero) then
1138        work(i) = a(i,i-1)
1139     else
1140        h = g + a(i,i-1)**2
1141        work(i) = sign(sqrt(h),a(i,i-1))
1142        h = h + a(i,i-1) * work(i)
1143        a(i,i-1) = a(i,i-1) + work(i)
1144        f = zero
1145        do j = 1, i-1
1146           g = zero
1147           do k = 1, i-1
1148              if (k .le. j) then
1149                 g = g + a(j,k) * a(i,k)
1150              else
1151                 g = g + a(k,j) * a(i,k)
1152              endif
1153           enddo
1154           work(j) = g / h
1155           f = f + work(j) * a(i,j)
1156        enddo
1157        do j = 1, i-1
1158           work(j) = work(j) - (f / (h + h)) * a(i,j)
1159           do k = 1, j
1160              a(j,k) = a(j,k) - a(i,j) * work(k) - work(j) * a(i,k)
1161           enddo
1162        enddo
1163     endif
1164  enddo
1165  work(2) = a(2,1)
1166  work(1) = zero
1167  eigen(2) = a(2,2)
1168  eigen(1) = a(1,1)
1169
1170  !---- Iterate on tridiagonal matrix.
1171  do i = 2, n
1172     work(i-1) = work(i)
1173  enddo
1174
1175  work(n) = zero
1176  f = zero
1177  b = zero
1178  do l = 1, n
1179     b = max(eps*(abs(eigen(l))+abs(work(l))),b)
1180     do m = l, n
1181        if (abs(work(m)) .le. b) go to 130
1182     enddo
1183     m = n
1184130  if (m .ne. l) then
1185        do it = 1, itmax
1186           p = (eigen(l+1) - eigen(l)) / (two * work(l))
1187           if (abs(p) .gt. big) then
1188              r = abs(p)
1189           else
1190              r = sqrt(p*p+one)
1191           endif
1192           h = eigen(l) - work(l) / (p + sign(r,p))
1193           do i = l, n
1194              eigen(i) = eigen(i) - h
1195           enddo
1196           f = f + h
1197           p = eigen(m)
1198           c = one
1199           s = zero
1200           do i = m-1, l, -1
1201              g = c * work(i)
1202              h = c * p
1203              r = sqrt(work(i)**2+p**2)
1204              work(i+1) = s * r
1205              s = work(i) / r
1206              c = p / r
1207              p = c * eigen(i) - s * g
1208              eigen(i+1) = h + s * (c * g + s * eigen(i))
1209           enddo
1210           work(l) = s * p
1211           eigen(l) = c * p
1212           if (abs(work(l)) .le. b) go to 170
1213        enddo
1214        nval = l - 1
1215        go to 300
1216     endif
1217170  p = eigen(l) + f
1218     do i = l, 2, -1
1219        if (p .ge. eigen(i-1)) go to 190
1220        eigen(i) = eigen(i-1)
1221     enddo
1222     i = 1
1223190  eigen(i) = p
1224  enddo
1225300 continue
1226
1227end subroutine symeig
1228subroutine dcopy(in,out,n)
1229  !----------------------------------------------------------------------*
1230  ! Purpose:                                                             *
1231  !   Copy arrays.                                                       *
1232  ! Input:                                                               *
1233  !   in  (double)    array to be copied.                                *
1234  !   n   (integer)   array length.                                      *
1235  ! Output:                                                              *
1236  !   out (double)    target array.                                      *
1237  !----------------------------------------------------------------------*
1238  implicit none
1239  integer n, i
1240  double precision in(*), out(*)
1241
1242  do i = 1, n
1243     out(i) = in(i)
1244  enddo
1245
1246end subroutine dcopy
1247subroutine dzero(vector,n)
1248  !----------------------------------------------------------------------*
1249  ! Purpose:                                                             *
1250  !   Zero an array.                                                     *
1251  ! Input:                                                               *
1252  !   n      (integer) array length.                                     *
1253  ! Input/output:                                                        *
1254  !   vector (double)  array to be zeroed.                               *
1255  !----------------------------------------------------------------------*
1256  implicit none
1257  integer n, i
1258  double precision vector(*),zero
1259  parameter(zero=0d0)
1260
1261  do i = 1, n
1262     vector(i) = zero
1263  enddo
1264
1265end subroutine dzero
1266subroutine aawarn(rout,text)
1267  implicit none
1268  !----------------------------------------------------------------------*
1269  ! Purpose:                                                             *
1270  !   Print warning message.                                             *
1271  ! Input:                                                               *
1272  !   ROUT      (char)    Calling routine name.                          *
1273  !   TEXT      (char)    Message.                                       *
1274  !----------------------------------------------------------------------*
1275  character(*) rout,text
1276
1277  print *,  '++++++ warning: ',rout,text
1278  call augmentfwarn()
1279
1280end subroutine aawarn
1281subroutine aafail(rout,text)
1282  implicit none
1283  !----------------------------------------------------------------------*
1284  ! Purpose:                                                             *
1285  !   Print fatal error message.                                         *
1286  ! Input:                                                               *
1287  !   ROUT      (char)    Calling routine name.                          *
1288  !   TEXT      (char)    Message.                                       *
1289  !----------------------------------------------------------------------*
1290  character(*) rout,text
1291  print *,' '
1292  print *,  '+-+-+- fatal: ',rout,text
1293  print *,' '
1294  print *,' '
1295  stop  1
1296
1297end subroutine aafail
1298double precision function proxim(x,y)
1299
1300
1301  !----------------------------------------------------------------------*
1302  !   Proximity function of x and y.                                     *
1303  !   If angle is larger than pi between vector x and y, 2pi is added to *
1304  !   to this angle                                                      *
1305  !----------------------------------------------------------------------*
1306  implicit none
1307  double precision x,y,twopi,get_variable
1308
1309  twopi=get_variable('twopi ')
1310  proxim = x+twopi*anint((y-x)/twopi)
1311
1312end function proxim
1313character(48) function charconv(tint)
1314  !----------------------------------------------------------------------*
1315  ! purpose:                                                             *
1316  !   converts integer array to string (based on ascii)                  *
1317  ! input:                                                               *
1318  !   tint  (int array)  1 = length, rest = string                       *
1319  !----------------------------------------------------------------------*
1320  implicit none
1321  integer tint(*)
1322  integer i, j, m, n
1323  parameter (m = 128)
1324  character(m) letter
1325  data letter / &
1326       '                                !"#$%&''()*+,-./0123456789:;<=>?@&
1327       &ABCDEFGHIJKLMNOPQRSTUVWXYZ[ ]^_`abcdefghijklmnopqrstuvwxyz{|}~'/
1328  charconv = ' '
1329  n = tint(1)
1330  do i = 1, n
1331     j = tint(i+1)
1332     if (j .lt. m)  charconv(i:i) = letter(j:j)
1333  enddo
1334end function charconv
1335subroutine laseig(fm,reeig,aieig,am)
1336  implicit none
1337  !----------------------------------------------------------------------*
1338  ! Purpose:                                                             *
1339  !   Return eigenvalues and eigenvectors of a 4x4 matrix.               *
1340  ! Input:                                                               *
1341  !   FM(6,6)   (real)    Matrix to be transformed.                      *
1342  ! Output:                                                              *
1343  !   REEIG(6)  (real)    Real parts of eigenvalues.                     *
1344  !   AIEIG(6)  (real)    Imaginary parts of eigenvalues.                *
1345  !   AM(6,6)   (real)    Transforming matrix, contains eigenvectors.    *
1346  !----------------------------------------------------------------------*
1347  integer i,ihi,ilo,info,ipind,iqind,j,k,mdim,nn,kpnt(6)
1348  double precision fm(6,6),reeig(6),aieig(6),am(6,6),aival(6),big,c,&
1349       d(6),dx,dy,pb,reval(6),s,tm(6,6),zero,one
1350  parameter(zero=0d0,one=1d0,ilo=1,ihi=4,mdim=6,nn=4)
1351
1352  !---- Compute eigenvalues and vectors.
1353  call m66cpy(fm,tm)
1354  call m66one(am)
1355  call orthes(mdim,nn,ilo,ihi,tm,d)
1356  call ortran(mdim,nn,ilo,ihi,tm,d,am)
1357  call hqr2(mdim,nn,ilo,ihi,tm,reval,aival,am,info)
1358  if (info .ne. 0) then
1359     write (6, 910) ((fm(i,k), k = 1, 6), i = 1, 6)
1360910  format('Unable to find eigenvalues for matrix:'/(6f12.6))
1361     call aafail('LASEIG',' Unable to find eigenvalues for matrix')
1362     go to 999
1363  endif
1364  !---- Normalize the eigenvectors.
1365  do k = 1, 5, 2
1366     pb = zero
1367     do ipind = 2, 6, 2
1368        iqind = ipind - 1
1369        pb = pb + am(iqind,k) * am(ipind,k+1) - am(ipind,k) * am(iqind,k+1)
1370     enddo
1371     s = sqrt(abs(pb))
1372     if (pb .lt. zero) then
1373        aival(k) = - aival(k)
1374        aival(k+1) = - aival(k+1)
1375     endif
1376     do i = 1, 6
1377        am(i,k)   = am(i,k) / s
1378        am(i,k+1) = am(i,k+1) * (s / pb)
1379     enddo
1380  enddo
1381  !---- Sort these eigenvectors.
1382  call m66cpy(am,tm)
1383  !---- Find the eigenvectors with the largest vertical component.
1384  big = zero
1385  kpnt(3) = 1
1386  do i = 1, 3, 2
1387     c = tm(3,i)**2 + tm(3,i+1)**2 + tm(4,i)**2 + tm(4,i+1)**2
1388     if (c .gt. big) then
1389        big = c
1390        kpnt(3) = i
1391     endif
1392  enddo
1393  !---- Find the remaining vector.
1394  do i = 1, 3, 2
1395     if (i .ne. kpnt(3)) kpnt(1) = i
1396  enddo
1397  !---- Reorder vectors.
1398  do i = 1, 3, 2
1399     k = kpnt(i)
1400     reeig(i) = reval(k)
1401     aieig(i) = aival(k)
1402     reeig(i+1) = reval(k+1)
1403     aieig(i+1) = aival(k+1)
1404     do j = 1, 6
1405        am(j,i) = tm(j,k)
1406        am(j,i+1) = tm(j,k+1)
1407     enddo
1408  enddo
1409  reeig(5) = one
1410  aieig(5) = zero
1411  reeig(6) = one
1412  aieig(6) = zero
1413  !---- Rephase the result.
1414  call m66one(tm)
1415  dx = sqrt(am(1,1)**2 + am(1,2)**2)
1416  tm(1,1) = am(1,1) / dx
1417  tm(2,1) = am(1,2) / dx
1418  tm(1,2) = - tm(2,1)
1419  tm(2,2) = tm(1,1)
1420  dy = sqrt(am(3,3)**2 + am(3,4)**2)
1421  tm(3,3) = am(3,3) / dy
1422  tm(4,3) = am(3,4) / dy
1423  tm(3,4) = - tm(4,3)
1424  tm(4,4) = tm(3,3)
1425  call m66mpy(am,tm,am)
1426999 end subroutine laseig
1427subroutine ladeig(fm,reeig,aieig,am)
1428  implicit none
1429  !----------------------------------------------------------------------*
1430  ! Purpose:                                                             *
1431  !   Return eigenvalues and eigenvectors of a 6x6 matrix.               *
1432  ! Input:                                                               *
1433  !   FM(6,6)   (real)    Matrix to be transformed.                      *
1434  ! Output:                                                              *
1435  !   REEIG(6)  (real)    Real parts of eigenvalues.                     *
1436  !   AIEIG(6)  (real)    Imaginary parts of eigenvalues.                *
1437  !   AM(6,6)   (real)    Transforming matrix, contains eigenvectors.    *
1438  !----------------------------------------------------------------------*
1439  integer i,ihi,ilo,info,j,k,mdim,nn,kpnt(6)
1440  double precision fm(6,6),reeig(6),aieig(6),am(6,6),aival(6),big,c,&
1441       d(6),dt,dx,dy,pb,reval(6),s,tm(6,6),zero
1442  parameter(zero=0d0,ilo=1,ihi=6,mdim=6,nn=6)
1443
1444  !---- Compute eigenvalues and eigenvectors.
1445  call m66cpy(fm,tm)
1446  call orthes(mdim,nn,ilo,ihi,tm,d)
1447  call ortran(mdim,nn,ilo,ihi,tm,d,am)
1448  call hqr2(mdim,nn,ilo,ihi,tm,reval,aival,am,info)
1449  if (info .ne. 0) then
1450     write (6, 910) ((fm(i,k), k = 1, 6), i = 1, 6)
1451910  format('Unable to find eigenvalues for matrix:'/(6f12.6))
1452     call aafail('LADEIG',' Unable to find eigenvalues for matrix')
1453     go to 9999
1454  endif
1455  !---- Normalize the eigenvectors.
1456  do k = 1, 5, 2
1457     pb = zero
1458     do i = 1, 5, 2
1459        pb = pb + am(i,k) * am(i+1,k+1) - am(i+1,k) * am(i,k+1)
1460     enddo
1461     s = sqrt(abs(pb))
1462     if (pb .lt. zero) then
1463        aival(k) = - aival(k)
1464        aival(k+1) = - aival(k+1)
1465     endif
1466     do i = 1, 6
1467        am(i,k)   = am(i,k) / s
1468        am(i,k+1) = am(i,k+1) * (s / pb)
1469     enddo
1470  enddo
1471  !---- Copy vectors to temporary array.
1472  call m66cpy(am,tm)
1473  !---- Find the vector with the largest vertical component.
1474  big = zero
1475  kpnt(3) = 1
1476  do i = 1, 5, 2
1477     c = tm(3,i)**2 + tm(3,i+1)**2 + tm(4,i)**2 + tm(4,i+1)**2
1478     if (c .gt. big) then
1479        big = c
1480        kpnt(3) = i
1481     endif
1482  enddo
1483  !---- Find  the vector with the largest horizontal component.
1484  kpnt(1) = 1
1485  big = zero
1486  do i = 1, 5, 2
1487     if (i .ne. kpnt(3)) then
1488        c = tm(1,i)**2 + tm(1,i+1)**2 + tm(2,i)**2 + tm(2,i+1)**2
1489        if (c .gt. big) then
1490           big = c
1491           kpnt(1) = i
1492        endif
1493     endif
1494  enddo
1495  !---- Find the remaining vector.
1496  do i = 1, 5, 2
1497     if (i .ne. kpnt(3)  .and.  i .ne. kpnt(1)) kpnt(5) = i
1498  enddo
1499  !---- Reorder vectors.
1500  do i = 1, 5, 2
1501     k = kpnt(i)
1502     reeig(i) = reval(k)
1503     aieig(i) = aival(k)
1504     reeig(i+1) = reval(k+1)
1505     aieig(i+1) = aival(k+1)
1506     do j = 1, 6
1507        am(j,i) = tm(j,k)
1508        am(j,i+1) = tm(j,k+1)
1509     enddo
1510  enddo
1511  !---- Rephase the result.
1512  call m66one(tm)
1513  dx = sqrt(am(1,1)**2 + am(1,2)**2)
1514  tm(1,1) = am(1,1) / dx
1515  tm(2,1) = am(1,2) / dx
1516  tm(1,2) = - tm(2,1)
1517  tm(2,2) = tm(1,1)
1518  dy = sqrt(am(3,3)**2 + am(3,4)**2)
1519  tm(3,3) = am(3,3) / dy
1520  tm(4,3) = am(3,4) / dy
1521  tm(3,4) = - tm(4,3)
1522  tm(4,4) = tm(3,3)
1523  dt = sqrt(am(5,5)**2 + am(5,6)**2)
1524  tm(5,5) = am(5,5) / dt
1525  tm(6,5) = am(5,6) / dt
1526  tm(5,6) = - tm(6,5)
1527  tm(6,6) = tm(5,5)
1528  call m66mpy(am,tm,am)
15299999 end subroutine ladeig
1530subroutine orthes(ndim,n,ilow,iupp,a,d)
1531  implicit none
1532  !----------------------------------------------------------------------*
1533  ! Purpose:                                                             *
1534  !   Converts an unsymmetric real matrix, A, to upper Hessenberg form   *
1535  !   applying successive orthogonal transformations.                    *
1536  !                                                                      *
1537  !   Translation of the ALGOL procedure ORTHES in:                      *
1538  !   Handbook Series Linear Algebra,                                    *
1539  !   Num. Math. 12, 349-368 (1968) by R. S. Martin and J. H. Wilkinson. *
1540  ! Input:                                                               *
1541  !   N         (integer) Order of the matrix A.                         *
1542  !   ILOW,IUPP (integer) Determine a submatrix, set by BALANC.          *
1543  !                       May be set to 1 and N respectively.            *
1544  !   A(NDIM,N) (real)    Input matrix.                                  *
1545  ! Output:                                                              *
1546  !   A(NDIM,N) (real)    The matrix A, converted to upper Hessenberg.   *
1547  !                       The lower triangle contains information        *
1548  !                       about the orthogonal transformations.          *
1549  !   D(N)      (real)    Further information.                           *
1550  !----------------------------------------------------------------------*
1551  integer i,ilow,iupp,j,m,n,ndim
1552  double precision a(ndim,n),d(n),f,g,h,scale,zero
1553  parameter(zero=0d0)
1554
1555  do m = ilow + 1, iupp - 1
1556     h = zero
1557     d(m) = zero
1558     !---- Find scale factor.
1559     scale = zero
1560     do i = m, iupp
1561        scale = scale + abs(a(i,m-1))
1562     enddo
1563     if (scale .ne. zero) then
1564        do i = iupp, m, - 1
1565           d(i) = a(i,m-1) / scale
1566           h = h + d(i) * d(i)
1567        enddo
1568        g = sign(sqrt(h),d(m))
1569        h = h + d(m) * g
1570        d(m) = d(m) + g
1571        !---- Form (I - (u*uT) / h) * A.
1572        do j = m, n
1573           f = zero
1574           do i = iupp, m, - 1
1575              f = f + d(i) * a(i,j)
1576           enddo
1577           f = f / h
1578           do i = m, iupp
1579              a(i,j) = a(i,j) - f * d(i)
1580           enddo
1581        enddo
1582        !---- Form (I - (u*uT) / h) * A * (I - (u*uT) / h).
1583        do i = 1, iupp
1584           f = zero
1585           do j = iupp, m, - 1
1586              f = f + d(j) * a(i,j)
1587           enddo
1588           f = f / h
1589           do j = m, iupp
1590              a(i,j) = a(i,j) - f * d(j)
1591           enddo
1592        enddo
1593        d(m) = scale * d(m)
1594        a(m,m-1) = - scale * g
1595     endif
1596  enddo
1597end subroutine orthes
1598subroutine ortran(ndim,n,ilow,iupp,h,d,v)
1599  implicit none
1600  !----------------------------------------------------------------------*
1601  ! Purpose:                                                             *
1602  !   Accumulate the orthogonal similarity transformation used by        *
1603  !   ORTHES to reduce a general real matrix A to upper Hessenberg form. *
1604  !                                                                      *
1605  !   Translation of the ALGOL procedure ORTRANS in:                     *
1606  !   Handbook Series Linear Algebra,                                    *
1607  !   Num. Math. 16, 181-204 (1970) by G. Peters and J. H. Wilkinson.    *
1608  ! Input:                                                               *
1609  !   N         (integer) Order of the matrices A and V.                 *
1610  !   ILOW,IUPP (integer) Determine a sub-matrix set by BALANC.          *
1611  !                       May be set to 1 and N respectively.            *
1612  !   H(NDIM,N) (real)    The matrix resulting from running ORTHES.      *
1613  !   D(N)      (real)    Further information about the transformation.  *
1614  ! Output:                                                              *
1615  !   V(NDIM,N) (real)    The accumulated transformation.                *
1616  !   D(N)      (real)    Destroyed.                                     *
1617  !----------------------------------------------------------------------*
1618  integer i,ilow,iupp,j,k,m,n,ndim
1619  double precision d(n),h(ndim,n),v(ndim,n),x,y,zero,one
1620  parameter(zero=0d0,one=1d0)
1621
1622  !---- Initialize V to identity matrix.
1623  do i = 1, n
1624     do j = 1, n
1625        v(i,j) = zero
1626     enddo
1627     v(i,i) = one
1628  enddo
1629  !---- Accumulate transformations.
1630  do k = iupp - 2, ilow, - 1
1631     m = k + 1
1632     y = h(m,k)
1633     if (y .ne. zero) then
1634        y = y * d(m)
1635        do i = k + 2, iupp
1636           d(i) = h(i,k)
1637        enddo
1638        do j = m, iupp
1639           x = zero
1640           do i = m, iupp
1641              x = x + d(i) * v(i,j)
1642           enddo
1643           x = x / y
1644           do i = m, iupp
1645              v(i,j) = v(i,j) + x * d(i)
1646           enddo
1647        enddo
1648     endif
1649  enddo
1650end subroutine ortran
1651subroutine hqr2(ndim,n,ilow,iupp,h,wr,wi,vecs,ierr)
1652  use max_iterate
1653  implicit none
1654  !----------------------------------------------------------------------*
1655  ! Purpose:                                                             *
1656  !   Finds eigenvalues and eigenvectors of an unsymmetric real matrix,  *
1657  !   A which has been reduced to upper Hessenberg form, H, by the       *
1658  !   subroutine ORTHES. The orthogonal transformations must be placed   *
1659  !   in the array VECS by subroutine ORTRAN.                            *
1660  !                                                                      *
1661  !   Translation of the ALGOL procedure HQR2 in:                        *
1662  !   Handbook Series Linear Algebra,                                    *
1663  !   Num. Math. 16, 181 - 204 (1970) by G. Peters and J. H. Wilkinson.  *
1664  ! Input:                                                               *
1665  !   N         (integer) Order of the Hessenberg matrix H.              *
1666  !   ILOW,IUPP (integer)                                                *
1667  !   H(NDIM,N) (real)    The Hessenberg matrix produced by ORTHES.      *
1668  !   VECS(NDIM,N) (real) A square matrix of order N containing the      *
1669  !                       similarity transformation from A to H          *
1670  ! Output:                                                              *
1671  !   H(NDIM,N) (real)    Modified.                                      *
1672  !   WR(N)     (real)    Real parts of eigenvalues of H (or A).         *
1673  !   WI(N)     (real)    Imaginary parts of eigenvalues of H (or A).    *
1674  !   VECS(NDIM,N) (real) The unnormalized eigenvectors of A.            *
1675  !                       Complex vectors are stored as pairs of reals.  *
1676  !----------------------------------------------------------------------*
1677  integer i,ien,ierr,ilow,its,iupp,j,k,l,m,n,na,ndim
1678  double precision den,h(ndim,n),hnorm,p,q,r,ra,s,sa,t,temp,tempi,  &
1679       tempr,vecs(ndim,n),vi,vr,w,wi(n),wr(n),x,y,z,epsmch,zero,one,two, &
1680       triqua,fac1
1681  parameter(epsmch=1d-16,zero=0d0,one=1d0,two=2d0,triqua=.75d0,fac1=.4375d0)
1682
1683  !Initialize
1684  z=zero
1685  s=zero
1686  p=zero
1687  q=zero
1688  r=zero
1689
1690  ierr = 0
1691  !---- Store isolated roots.
1692  do i = 1, n
1693     if (i .lt. ilow  .or.  i .gt. iupp) then
1694        wr(i) = h(i,i)
1695        wi(i) = zero
1696     endif
1697  enddo
1698  ien = iupp
1699  t = zero
1700  !---- Next eigenvalue.
170160 if (ien .ge. ilow) then
1702     its = 0
1703     na = ien - 1
1704     !---- Next iteration; look for single small sub-diagonal element.
170570   continue
1706     do l = ien, ilow + 1, -1
1707        if (abs(h(l,l-1)) .le. epsmch * (abs(h(l-1,l-1)) + abs(h(l,l)))) go to 100
1708     enddo
1709     l = ilow
1710100  continue
1711     x = h(ien,ien)
1712     if (l .eq. ien) go to 270
1713     y = h(na,na)
1714     w = h(ien,na) * h(na,ien)
1715     if (l .eq. na) go to 280
1716     if (its .eq. MAXITER) then
1717        write(6,*) "Maximum Iteration exceeded in HQR2, increase MAXITER: ",MAXITER
1718        ierr = ien
1719        go to 9999
1720     endif
1721     !---- Form exceptional shift.
1722     if (its .eq. 10  .or.  its .eq. 20) then
1723        t = t + x
1724        do i = ilow, ien
1725           h(i,i) = h(i,i) - x
1726        enddo
1727        s = abs(h(ien,na)) + abs(h(na,ien-2))
1728        x = triqua * s
1729        y = x
1730        w = - fac1 * s * s
1731     endif
1732     its = its + 1
1733     !---- Look for two consecutive small sub-diagonal elements.
1734     do m = ien - 2, l, - 1
1735        z = h(m,m)
1736        r = x - z
1737        s = y - z
1738        p = (r * s - w) / h(m+1,m) + h(m,m+1)
1739        q = h(m+1,m+1) - z - r - s
1740        r = h(m+2,m+1)
1741        s = abs(p) + abs(q) + abs(r)
1742        p = p / s
1743        q = q / s
1744        r = r / s
1745        if (m .eq. l) go to 150
1746        if (abs(h(m,m-1)) * (abs(q) + abs(r)) .le. epsmch * abs(p)    &
1747             * (abs(h(m-1,m-1)) + abs(z) + abs(h(m+1,m+1)))) go to 150
1748     enddo
1749150  continue
1750     h(m+2,m) = zero
1751     do i = m + 3, ien
1752        h(i,i-2) = zero
1753        h(i,i-3) = zero
1754     enddo
1755     !---- Double QR step involving rows L to IEN and columns M to IEN.
1756     do k = m, na
1757        if (k .ne. m) then
1758           p = h(k,k-1)
1759           q = h(k+1,k-1)
1760           if (k .ne. na) then
1761              r = h(k+2,k-1)
1762           else
1763              r = zero
1764           endif
1765           x = abs(p) + abs(q) + abs(r)
1766           if (x .eq. zero) go to 260
1767           p = p / x
1768           q = q / x
1769           r = r / x
1770        endif
1771        s = sign(sqrt(p**2+q**2+r**2),p)
1772        if (k .ne. m) then
1773           h(k,k-1) = - s * x
1774        else if (l .ne. m) then
1775           h(k,k-1) = - h(k,k-1)
1776        endif
1777        p = p + s
1778        x = p / s
1779        y = q / s
1780        z = r / s
1781        q = q / p
1782        r = r / p
1783        !---- Row modification.
1784        do j = k, n
1785           p = h(k,j) + q * h(k+1,j)
1786           if (k .ne. na) then
1787              p = p + r * h(k+2,j)
1788              h(k+2,j) = h(k+2,j) - p * z
1789           endif
1790           h(k+1,j) = h(k+1,j) - p * y
1791           h(k,j) = h(k,j) - p * x
1792        enddo
1793        !---- Column modification.
1794        j = min(ien,k+3)
1795        do i = 1, j
1796           p = x * h(i,k) + y * h(i,k+1)
1797           if (k .ne. na) then
1798              p = p + z * h(i,k+2)
1799              h(i,k+2) = h(i,k+2) - p * r
1800           endif
1801           h(i,k+1) = h(i,k+1) - p * q
1802           h(i,k) = h(i,k) - p
1803        enddo
1804        !---- Accumulate transformations.
1805        do i = ilow, iupp
1806           p = x * vecs(i,k) + y * vecs(i,k+1)
1807           if (k .ne. na) then
1808              p = p + z * vecs(i,k+2)
1809              vecs(i,k+2) = vecs(i,k+2) - p * r
1810           endif
1811           vecs(i,k+1) = vecs(i,k+1) - p * q
1812           vecs(i,k) = vecs(i,k) - p
1813        enddo
1814260     continue
1815     enddo
1816     !---- Go to next iteration.
1817     go to 70
1818     !==== One real root found.
1819270  h(ien,ien) = x + t
1820     wr(ien) = h(ien,ien)
1821     wi(ien) = zero
1822     ien = na
1823     go to 60
1824     !==== Two roots (real pair or complex conjugate) found.
1825280  p = (y - x) / two
1826     q = p**2 + w
1827     z = sqrt(abs(q))
1828     x = x + t
1829     h(ien,ien) = x
1830     h(na,na) = y + t
1831     !---- Real pair.
1832     if (q .gt. zero) then
1833        z = p + sign(z,p)
1834        wr(na) = x + z
1835        wr(ien) = x - w / z
1836        wi(na) = zero
1837        wi(ien) = zero
1838        x = h(ien,na)
1839        r = sqrt(x**2+z**2)
1840        p = x / r
1841        q = z / r
1842        !---- Row modification.
1843        do j = na, n
1844           z = h(na,j)
1845           h(na,j) = q * z + p * h(ien,j)
1846           h(ien,j) = q * h(ien,j) - p * z
1847        enddo
1848        !---- Column modification.
1849        do i = 1, ien
1850           z = h(i,na)
1851           h(i,na) = q * z + p * h(i,ien)
1852           h(i,ien) = q * h(i,ien) - p * z
1853        enddo
1854        !---- Accumulate transformations.
1855        do i = ilow, iupp
1856           z = vecs(i,na)
1857           vecs(i,na) = q * z + p * vecs(i,ien)
1858           vecs(i,ien) = q * vecs(i,ien) - p * z
1859        enddo
1860        !---- Complex pair.
1861     else
1862        wr(na) = x + p
1863        wr(ien) = x + p
1864        wi(na) = z
1865        wi(ien) = -z
1866     endif
1867     !----- Go to next root.
1868     ien = ien - 2
1869     go to 60
1870  endif
1871  !==== Compute matrix norm.
1872  hnorm = zero
1873  k = 1
1874  do i = 1, n
1875     do j = k, n
1876        hnorm = hnorm + abs(h(i,j))
1877     enddo
1878     k = i
1879  enddo
1880  !==== Back substitution.
1881  do ien = n, 1, -1
1882     p = wr(ien)
1883     q = wi(ien)
1884     na = ien - 1
1885     !---- Real vector.
1886     if (q .eq. zero) then
1887        m = ien
1888        h(ien,ien) = one
1889        do i = na, 1, -1
1890           w = h(i,i) - p
1891           r = h(i,ien)
1892           do j = m, na
1893              r = r + h(i,j) * h(j,ien)
1894           enddo
1895           if (wi(i) .lt. zero) then
1896              z = w
1897              s = r
1898           else
1899              m = i
1900              if (wi(i) .eq. zero) then
1901                 temp = w
1902                 if (w .eq. zero) temp = epsmch * hnorm
1903                 h(i,ien) = - r / temp
1904              else
1905                 x = h(i,i+1)
1906                 y = h(i+1,i)
1907                 q = (wr(i) - p)**2 + wi(i)**2
1908                 t = (x * s - z * r) / q
1909                 h(i,ien) = t
1910                 if (abs(x) .gt. abs(z)) then
1911                    h(i+1,ien) = - (r + w * t) / x
1912                 else
1913                    h(i+1,ien) = - (s + y * t) / z
1914                 endif
1915              endif
1916           endif
1917        enddo
1918        !---- Complex vector associated with lamda = P - i * Q.
1919     else if (q .lt. zero) then
1920        m = na
1921        if (abs(h(ien,na)) .gt. abs(h(na,ien))) then
1922           h(na,na) = - (h(ien,ien) - p) / h(ien,na)
1923           h(na,ien) = - q / h(ien,na)
1924        else
1925           den = (h(na,na) - p)**2 + q**2
1926           h(na,na) = - h(na,ien) * (h(na,na) - p) / den
1927           h(na,ien) = h(na,ien) * q / den
1928        endif
1929        h(ien,na) = one
1930        h(ien,ien) = zero
1931        do i = ien - 2, 1, - 1
1932           w = h(i,i) - p
1933           ra = h(i,ien)
1934           sa = zero
1935           do j = m, na
1936              ra = ra + h(i,j) * h(j,na)
1937              sa = sa + h(i,j) * h(j,ien)
1938           enddo
1939           if (wi(i) .lt. zero) then
1940              z = w
1941              r = ra
1942              s = sa
1943           else
1944              m = i
1945              if (wi(i) .eq. zero) then
1946                 den = w**2 + q**2
1947                 h(i,na) = - (ra * w + sa * q) / den
1948                 h(i,ien) = (ra * q - sa * w) / den
1949              else
1950                 x = h(i,i+1)
1951                 y = h(i+1,i)
1952                 vr = (wr(i) - p)**2 + wi(i)**2 - q**2
1953                 vi = two * (wr(i) - p) * q
1954                 if (vr .eq. zero  .and.  vi .eq. zero) then
1955                    vr = epsmch * hnorm                                   &
1956                         * (abs(w) + abs(q) + abs(x) + abs(y) + abs(z))
1957                 endif
1958                 tempr = x * r - z * ra + q * sa
1959                 tempi = x * s - z * sa - q * ra
1960                 den = vr**2 + vi**2
1961                 h(i,na) = (tempr * vr + tempi * vi) / den
1962                 h(i,ien) = (tempi * vr - tempr * vi) / den
1963                 if (abs(x) .gt. abs(z) + abs(q)) then
1964                    h(i+1,na) = (- ra - w * h(i,na) + q * h(i,ien)) / x
1965                    h(i+1,ien) = (- sa - w * h(i,ien) - q * h(i,na)) / x
1966                 else
1967                    tempr = - r - y * h(i,na)
1968                    tempi = - s - y * h(i,ien)
1969                    den = z**2 + q**2
1970                    h(i+1,na) = (tempr * z + tempi * q) / den
1971                    h(i+1,ien) = (tempi * z - tempr * q) / den
1972                 endif
1973              endif
1974           endif
1975        enddo
1976     endif
1977  enddo
1978  !==== Vectors of isolated roots.
1979  do i = 1, n
1980     if (i .lt. ilow  .or.  i .gt. iupp) then
1981        do j = i, n
1982           vecs(i,j) = h(i,j)
1983        enddo
1984     endif
1985  enddo
1986  !==== Multiply by transformation matrix to give eigenvectors of the
1987  !     original full matrix.
1988  do j = n, ilow, - 1
1989     m = min(j,iupp)
1990     if (wi(j) .lt. zero) then
1991        l = j - 1
1992        do i = ilow, iupp
1993           y = zero
1994           z = zero
1995           do k = ilow, m
1996              y = y + vecs(i,k) * h(k,l)
1997              z = z + vecs(i,k) * h(k,j)
1998           enddo
1999           vecs(i,l) = y
2000           vecs(i,j) = z
2001        enddo
2002     else if (wi(j) .eq. zero) then
2003        do i = ilow, iupp
2004           z = zero
2005           do k = ilow, m
2006              z = z + vecs(i,k) * h(k,j)
2007           enddo
2008           vecs(i,j) = z
2009        enddo
2010     endif
2011  enddo
20129999 end subroutine hqr2
2013
2014subroutine suelem(el, ve, we,tilt)
2015
2016  use twtrrfi
2017  implicit none
2018
2019  !----------------------------------------------------------------------*
2020  ! Purpose:                                                             *
2021  !   Compute Displacement and rotation for one element.                 *
2022  ! Output:                                                              *
2023  !   EL        (real)    Element length along design orbit.             *
2024  !   VE(3)     (real)    Displacement of exit w.r.t. entry.             *
2025  !   WE(3,3)   (real)    Rotation of exit w.r.t. entry.                 *
2026  ! Reference pointer used:                                              *
2027  !   LCELM     /REFER/   Current element bank.                          *
2028  ! Local links:                                                         *
2029  !   LSEQ                Beam lines sequence for a lump.                *
2030  !----------------------------------------------------------------------*
2031  ! Modified: 17-FEB-2005, F. Tecker                                     *
2032  !   Return tilt also for elements other than BENDs and MULTIPOLEs      *
2033  ! Modified: 28-DEC-1998, T. Raubenheimer (SLAC)                        *
2034  !   Added LCAVITY element at ISP 27                                    *
2035  !----------------------------------------------------------------------*
2036  integer code,nn
2037  double precision angle,cospsi,costhe,ds,dx,sinpsi,sinthe,tilt,    &
2038       ve(3),we(3,3),node_value,el,normal(0:maxmul),skew(0:maxmul)       &
2039       ,zero,one
2040  parameter(zero=0d0,one=1d0)
2041  !---- Branch on subprocess code.
2042  tilt = zero
2043  angle = zero
2044  code = node_value('mad8_type ')
2045  if(code.eq.39) code=15
2046  if(code.eq.38) code=24
2047  go to ( 10,  20,  20,  40,  50,  60,  70,  80,  90, 100,          &
2048       110, 120, 130, 140, 150, 160, 170, 180, 190, 200,                 &
2049       210, 220, 230, 240, 250,  20, 270, 280, 290, 300,                 &
2050       310, 310, 310, 310, 310, 310, 310, 310, 310, 310), code
2051
2052  !---- elements without tilt attribute
2053  !---- Drift space.
205410 continue
2055
2056  !---- Arbitrary matrix.
205740 continue
2058
2059  !---- Solenoid.
206090 continue
2061
2062  !---- RF cavity.
2063100 continue
2064
2065  !---- Monitors.
2066170 continue
2067180 continue
2068190 continue
2069
2070  !---- Marker.
2071250 continue
2072
2073  !---- Beam-beam.
2074220 continue
2075
2076  !---- lcavity
2077270 continue
2078
2079  !---- Reserved.
2080280 continue
2081290 continue
2082300 continue
2083
2084  !---- Lump.
2085230 continue
2086
2087  !---- User-defined elements.
2088310 continue
2089
2090  !     not necessary to distinguish between elements with or w/o tilt attribute
2091  !     see below (FT 17.2.05)
2092  !     goto 400
2093
2094
2095  !---- elements with tilt attribute
2096  !---- Quadrupole.
209750 continue
2098
2099  !---- Sextupole.
210060 continue
2101
2102  !---- Octupole.
210370 continue
2104
2105  !---- Electrostatic separator.
2106110 continue
2107
2108  !---- Kickers.
2109140 continue
2110150 continue
2111160 continue
2112
2113  !---- Apertures.
2114200 continue
2115210 continue
2116
2117  !---- Beam instrument.
2118240 continue
2119
2120  !---- get tilt attribute
2121  !---- checked to work for elements without this attribute (FT 17.2.05)
2122  !---- OK with F.Sschmidt
2123  tilt =  node_value('tilt ')
2124
2125  !---- calculate matrix for straight elements
2126  continue
2127  ve(1) = zero
2128  ve(2) = zero
2129  ve(3) = el
2130  we(1,1) = one
2131  we(2,1) = zero
2132  we(3,1) = zero
2133  we(1,2) = zero
2134  we(2,2) = one
2135  we(3,2) = zero
2136  we(1,3) = zero
2137  we(2,3) = zero
2138  we(3,3) = one
2139  go to 500
2140  !****** end of straight elements ***************
2141
2142  !---- multipoles , introduced  17.09.02 / AV
214380 continue
2144  call dzero(normal,maxmul+1)
2145  call dzero(skew,maxmul+1)
2146  call get_node_vector('knl ',nn,normal)
2147  !-----  dipole_bv introduced to suppress SU in MADX input (AV  7.10.02)
2148  !      angle = normal(0)*node_value('dipole_bv ')
2149  angle = normal(0)*node_value('other_bv ')
2150  if (abs(angle) .lt. 1d-13) then
2151     tilt = zero
2152  else
2153     tilt =  node_value('tilt ')
2154  endif
2155  ! As el=0, there is no dx and no ds
2156  dx = zero
2157  ds = zero
2158  go to 490
2159
2160  !---- Any kind of  bend.
216120 continue
2162  !--------------  dipole_bv introduced to suppress SU (AV  7.10.02)
2163  !      angle = node_value('angle ')*node_value('dipole_bv ')
2164  angle = node_value('angle ')*node_value('other_bv ')
2165  !      print *,"SUELEM dipole : angle =",angle
2166  if (abs(angle) .lt. 1d-13) then
2167     dx = zero
2168     ds = el
2169     tilt = zero
2170  else
2171     tilt =  node_value('tilt ')
2172     !      print *,"SUELEM dipole : tilt =",tilt," length= ",el
2173     dx = el * (cos(angle)-one)/angle
2174     ds = el * sin(angle)/angle
2175  endif
2176  !      print *,"SUELEM dipole : tilt =",tilt," length= ",&
2177  !     el," angv = ",angv," bv =",node_value('dipole_bv ')
2178  !     el," angv = ",angv," bv =",node_value('other_bv ')
2179  go to 490
2180
2181  !---- Rotation around S-axis. SPECIAL CASE
2182120 continue
2183  print *,"SROT"
2184  tilt = node_value('angle ')
2185  ve(1) = zero
2186  ve(2) = zero
2187  ve(3) = zero
2188  we(1,1) =  cos(tilt)
2189  we(2,1) =  sin(tilt)
2190  we(3,1) = zero
2191  we(1,2) = -sin(tilt)
2192  we(2,2) = cos(tilt)
2193  we(3,2) = zero
2194  we(1,3) = zero
2195  we(2,3) = zero
2196  we(3,3) = one
2197  go to 500
2198
2199  !---- Rotation around Y-axis.  QUESTIONABLE USEFULNESS  !!!!!!!!!!!!!
2200130 continue
2201  print *,"YROT"
2202  dx = node_value('angle ')
2203  print *,"angle",dx
2204  tilt = zero
2205  ve(1) = zero
2206  ve(2) = zero
2207  ve(3) = zero
2208  we(1,1) =  cos(dx)
2209  we(2,1) = zero
2210  we(3,1) = sin(dx)
2211  we(1,2) = zero
2212  we(2,2) = zero
2213  we(3,2) = one
2214  we(1,3) = -sin(dx)
2215  we(2,3) = zero
2216  we(3,3) = cos(dx)
2217  go to 500
2218
2219  !---- Common for bends and multipoles: Displacement and rotation matrix.
2220490 continue
2221  cospsi = cos(tilt)
2222  sinpsi = sin(tilt)
2223  costhe = cos(angle)
2224  sinthe = sin(angle)
2225  ve(1) = dx * cospsi
2226  ve(2) = dx * sinpsi
2227  ve(3) = ds
2228  we(1,1) = costhe * cospsi*cospsi + sinpsi*sinpsi
2229  we(2,1) = (costhe - one) * cospsi * sinpsi
2230  we(3,1) = sinthe * cospsi
2231  we(1,2) = we(2,1)
2232  we(2,2) = costhe * sinpsi*sinpsi + cospsi*cospsi
2233  we(3,2) =  sinthe * sinpsi
2234  we(1,3) = - we(3,1)
2235  we(2,3) = - we(3,2)
2236  we(3,3) = costhe
2237500 continue
2238end subroutine suelem
2239!-----------------  end of suelem subroutine --------------------------
2240
2241!**********************************************************************
2242subroutine sumtrx(the, phi, psi, w)
2243  implicit none
2244  !----------------------------------------------------------------------*
2245  ! Purpose:                                                             *
2246  !   Given three survey angles, compute rotation matrix.                *
2247  ! Input:                                                               *
2248  !   THE       (real)    Azimuthal angle.                               *
2249  !   PHI       (real)    Elevation angle.                               *
2250  !   PSI       (real)    Roll angle.                                    *
2251  ! Output:                                                              *
2252  !   W(3,3)    (real)    Rotation matrix.                               *
2253  !----------------------------------------------------------------------*
2254  double precision cosphi,cospsi,costhe,phi,psi,sinphi,sinpsi,sinthe,the,w(3,3)
2255
2256  costhe = cos(the)
2257  sinthe = sin(the)
2258  cosphi = cos(phi)
2259  sinphi = sin(phi)
2260  cospsi = cos(psi)
2261  sinpsi = sin(psi)
2262  w(1,1) = + costhe * cospsi - sinthe * sinphi * sinpsi
2263  w(1,2) = - costhe * sinpsi - sinthe * sinphi * cospsi
2264  w(1,3) =                     sinthe * cosphi
2265  w(2,1) =                              cosphi * sinpsi
2266  w(2,2) =                              cosphi * cospsi
2267  w(2,3) =                              sinphi
2268  w(3,1) = - sinthe * cospsi - costhe * sinphi * sinpsi
2269  w(3,2) = + sinthe * sinpsi - costhe * sinphi * cospsi
2270  w(3,3) =                     costhe * cosphi
2271
2272end subroutine sumtrx
2273!-----------------  end of sumtrx subroutine --------------------------
2274
2275!**********************************************************************
2276subroutine sutran(w, v, we)
2277  implicit none
2278  !----------------------------------------------------------------------*
2279  ! Purpose:                                                             *
2280  !   Transform rotation W and displacement V from entrance to exit.     *
2281  ! Input:                                                               *
2282  !   W(3,3)    (real)    Rotation matrix w.r.t. input system.           *
2283  !   V(3)      (real)    Displacement w.r.t. input system.              *
2284  !   WE(3,3)   (real)    Rotation matrix due to element.                *
2285  ! Output:                                                              *
2286  !   W(3,3)    (real)    Rotation matrix w.r.t. output system.          *
2287  !   V(3)      (real)    Displacement w.r.t. output system.             *
2288  !----------------------------------------------------------------------*
2289  integer i,k
2290  double precision v(3),vt(3),w(3,3),we(3,3),wt(3,3)
2291
2292  !---- VT := transpose(WE) * V;
2293  !     WT := transpose(WE) * W;
2294  do i = 1, 3
2295     vt(i) = we(1,i)*v(1) + we(2,i)*v(2) + we(3,i)*v(3)
2296     do k = 1, 3
2297        wt(i,k) = we(1,i)*w(1,k) + we(2,i)*w(2,k) + we(3,i)*w(3,k)
2298     enddo
2299  enddo
2300
2301  !---- V := VT       [= transpose(WE) * V];
2302  !     W := WT * WE  [= transpose(WE) * W * WE];
2303  do i = 1, 3
2304     v(i) = vt(i)
2305     do k = 1, 3
2306        w(i,k) = wt(i,1)*we(1,k) + wt(i,2)*we(2,k) + wt(i,3)*we(3,k)
2307     enddo
2308  enddo
2309end subroutine sutran
2310!-----------------  end of sutran subroutine --------------------------
2311integer function lastnb(t)
2312  !----------------------------------------------------------------------*
2313  ! Purpose:
2314  !   Find last non-blank in string
2315  !
2316  !----------------------------------------------------------------------*
2317  implicit none
2318
2319  character(*) t
2320  integer i
2321  do i = len(t), 1, -1
2322     if (t(i:i) .ne. ' ') goto 20
2323  enddo
2324  i = 1
232520 lastnb = i
2326end function lastnb
2327subroutine tmfoc(el,sk1,c,s,d,f)
2328  implicit none
2329  !----------------------------------------------------------------------*
2330  ! Purpose:                                                             *
2331  !   Compute linear focussing functions.                                *
2332  ! Input:                                                               *
2333  !   el        (double)  element length.                                *
2334  !   sk1       (double)  quadrupole strength.                           *
2335  ! Output:                                                              *
2336  !   c         (double)  cosine-like function.             c(k,l)       *
2337  !   s         (double)  sine-like function.               s(k,l)       *
2338  !   d         (double)  dispersion function.              d(k,l)       *
2339  !   f         (double)  integral of dispersion function.  f(k,l)       *
2340  !----------------------------------------------------------------------*
2341  double precision c,d,el,f,qk,qkl,qkl2,s,sk1,zero,one,two,six,     &
2342       twelve,twty,thty,foty2
2343  parameter(zero=0d0,one=1d0,two=2d0,six=6d0,twelve=12d0,twty=20d0, &
2344       thty=30d0,foty2=42d0)
2345
2346  !---- Initialize.
2347  qk = sqrt(abs(sk1))
2348  qkl = qk * el
2349  qkl2 = sk1 * el**2
2350  if (abs(qkl2) .le. 1e-2) then
2351     c = (one - qkl2 * (one - qkl2 / twelve) /  two)
2352     s = (one - qkl2 * (one - qkl2 / twty) /  six) * el
2353     d = (one - qkl2 * (one - qkl2 / thty) / twelve) * el**2 / two
2354     f = (one - qkl2 * (one - qkl2 / foty2) / twty) * el**3 / six
2355  else
2356     if (qkl2 .gt. zero) then
2357        c = cos(qkl)
2358        s = sin(qkl) / qk
2359     else
2360        c = cosh(qkl)
2361        s = sinh(qkl) / qk
2362     endif
2363     d = (one - c) / sk1
2364     f = (el  - s) / sk1
2365  endif
2366
2367end subroutine tmfoc
2368
2369subroutine f77flush(i,option)
2370  implicit none
2371  integer i,ios
2372  real a
2373  logical ostat, fexist,option
2374  character(20) faccess,fform
2375  character(255) fname
2376  character(1) c
2377  inquire(err=5,iostat=ios,unit=i,opened=ostat,exist=fexist)
2378  if (.not.ostat.or..not.fexist) return
2379  inquire(err=6,iostat=ios,unit=i,access=faccess,form=fform,name=fname)
2380  close (unit=i,err=7,iostat=ios)
2381  !     write (*,*) 'Re-opening ',i,' ',faccess,fform,fname
2382  open(err=8,iostat=ios,unit=i,access=faccess,form=fform,file=fname,status='old')
2383  if (option) then
2384     if (fform.eq.'FORMATTED') then
23853       read (i,100,err=9,iostat=ios,end=4) c
2386        go to 3
2387     else
23882       read (i,err=10,iostat=ios,end=1) a
2389        go to 2
2390     endif
23914    backspace i
23921    continue
2393  endif
2394  return
2395100 format (a1)
23965 write (*,*) ' F77FLUSH 1st INQUIRE FAILED with IOSTAT ',ios,' on UNIT ',i
2397  stop
23986 write (*,*) ' F77FLUSH 2nd INQUIRE FAILED with IOSTAT ', ios,' on UNIT ',i
2399  stop
24007 write (*,*) ' F77FLUSH CLOSE FAILED with IOSTAT ',ios,' on UNIT ',i
2401  stop
24028 write (*,*) ' F77FLUSH RE-OPEN FAILED with IOSTAT ',ios,' on UNIT ',i
2403  stop
24049 write (*,*) ' F77FLUSH FORMATTED READ FAILED with IOSTAT ',ios,' on UNIT ',i
2405  stop
240610 write (*,*) ' F77FLUSH UNFORMATTED READ FAILED with IOSTAT ',ios,' on UNIT ',i
2407  stop
2408end subroutine f77flush
2409
2410
2411subroutine seterrorflag(errorcode,from,descr)
2412  !----------------------------------------------------------------------*
2413  !     Purpose:                                                         *
2414  !     Puts global error flag in c code.                                *
2415  !     Input:                                                           *
2416  !     errorcode                                                        *
2417  !     from - name of a routine where the error occured                 *
2418  !     descr - description of the error that has occured                *
2419  !     Input/output:                                                    *
2420  !----------------------------------------------------------------------*
2421  implicit none
2422  integer :: errorcode
2423  character(*) :: from
2424  character(*) :: descr
2425  integer  n,m
2426
2427  n = LEN(from)
2428  m = LEN(descr)
2429  call seterrorflagfort(errorcode,from,n,descr,m)
2430
2431end subroutine seterrorflag
Note: See TracBrowser for help on using the repository browser.