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

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

import madx-5.01.00
File size: 66.6 KB
Line 
1!The Polymorphic Tracking Code
2!Copyright (C) Etienne Forest and CERN
3
4module S_fitting_new
5  USE ptc_spin
6  IMPLICIT NONE
7  public
8  integer:: m_turn,m_skip=0
9  integer :: with_c=1
10  TYPE fibre_monitor_data
11     type(fibre), pointer :: p    ! fibre location
12     integer, pointer ::  turn,kind  ! kind=1 x, kind = 2 y
13     real(dp), pointer :: bpm(:,:)  ! store fake experiment from alex_track_monitors
14     real(dp), pointer :: r(:,:)  ! store fake experiment from alex_track_monitors
15     real(dp), pointer :: xf(:,:)  ! real data put here
16     real(dp), pointer :: xn(:,:)  ! real data put here
17     real(dp), pointer :: mom(:,:)
18     real(dp), pointer :: A(:,:)
19     real(dp), pointer :: At(:,:)
20     logical full
21  END TYPE fibre_monitor_data
22
23  TYPE(fibre_monitor_data), allocatable :: monitors(:)
24
25contains
26
27  subroutine special_alex_main_ring(r,n_name,targ,sc)
28    implicit none
29    TYPE(layout), target, intent(inout):: R
30    integer  i1,i2,it1,it2,it3,it4
31    INTEGER I,N,NU,N2,NP2,mf,nt,NP,J,n_name
32    type(fibre), pointer :: p1,p2
33    TYPE(POL_BLOCK) QC(11)
34    TYPE(REAL_8) Y(6)
35    TYPE(DAMAP) ID
36    REAL(DP) X(6),targ(2),xx
37    TYPE(INTERNAL_STATE) STATE
38    LOGICAL(LP) U
39    type(normalform) nf
40    type(gmap) g
41    type(TAYLOR) T,eq(5)
42
43    REAL(DP) ALX,ALY,TA(2),sc
44
45    !targ(1)=22.43d0
46    !targ(2)=20.82d0
47
48    N=10
49    NU=11
50    if(.not.associated(r%t)) then
51       write(6,*) " thin lens lattice not made "
52       CALL MAKE_node_LAYOUT(r)
53    endif
54    p1=>r%start
55    do i=1,r%n
56       if(p1%mag%name(1:3)=='QDX') then
57          i1=i
58          exit
59       endif
60       p1=>p1%next
61    enddo
62    p2=>p1%next
63    do i=i1+1,r%n
64       if(p2%mag%name(1:3)=='QDX') then
65          i2=i
66          exit
67       endif
68       p2=>p2%next
69    enddo
70    !    call move_to(r,p1,i1)
71    !    call move_to(r,p2,i2)
72
73    write(6,*) p1%mag%name,p1%mag%p%nst
74    write(6,*) p2%mag%name,p2%mag%p%nst
75
76    IT1=p1%T1%POS+2 + (p1%mag%p%nst/2 )
77    IT2=p2%T1%POS+2 + (p2%mag%p%nst/2 )
78
79    write(6,*) i1,IT1,i2,IT2
80
81    if(mod(p1%mag%p%nst,2)/=0.or.mod(p2%mag%p%nst,2)/=0) then
82       write(6,*) " Even number of split needed for fitting in Alex_special "
83       stop 100
84    endif
85    !    call move_to(r,p1,i3)
86    !    call move_to(r,p2,i4)
87    p1=>p2%next
88    do i=i2+1,r%n
89       if(p1%mag%name(1:3)=='QFS') then
90          i1=i
91          exit
92       endif
93       p1=>p1%next
94    enddo
95
96    do i=i1,1,-1
97       if(p1%mag%name(1:3)=='QDX') then
98          i1=i
99          exit
100       endif
101       p1=>p1%previous
102    enddo
103
104
105    p2=>p1%next
106    do i=i1+1,r%n
107       if(p2%mag%name(1:3)=='QDX') then
108          i2=i
109          exit
110       endif
111       p2=>p2%next
112    enddo
113
114    write(6,*) p1%mag%name,p1%mag%p%nst
115    write(6,*) p2%mag%name,p2%mag%p%nst
116
117    IT3=p1%T1%POS+2 + (p1%mag%p%nst/2 )
118    IT4=p2%T1%POS+2 + (p2%mag%p%nst/2 )
119
120    write(6,*) i1,IT3,i2,IT4
121    if(mod(p1%mag%p%nst,2)/=0.or.mod(p2%mag%p%nst,2)/=0) then
122       write(6,*) " Even number of split needed for fitting in Alex_special "
123       stop 101
124    endif
125
126    DO I=1,NU
127       QC(I)=0
128       QC(I)%n_name=n_name
129    ENDDO
130    QC(1)%NAME='QFX'
131    QC(2)%NAME='QDX'
132    QC(3)%NAME='QFN'
133    QC(4)%NAME='QDN'
134    QC(5)%NAME='QFS'
135    QC(6)%NAME='QDS'
136    QC(7)%NAME='QFT'
137    QC(8)%NAME='QFP'
138    QC(9)%NAME='QDT'
139    QC(10)%NAME='QFR'
140    QC(11)%NAME='QDR'
141
142    DO I=1,NU
143       QC(I)%IBN(2)=I
144    ENDDO
145
146    DO I=1,NU
147       R=QC(I)
148    ENDDO
149
150111 STATE=DEFAULT0+ONLY_4D0
151
152    X=0.0_dp
153    CALL INIT(STATE,2,NU,BERZ,N2,NP2)
154    CALL ALLOC(ID); call alloc(nf);call alloc(Y);call alloc(eq,5);
155    ID=1
156    Y=X+ID
157    !( R,X,U,K,POS1,POS2,T1,T2,P1,P2,IN_P1,IN_P2,POS1_FIBRE,POS2_FIBRE)
158    !   CALL TRACK_BEAM_x(R,Y,U,+STATE,POS1=IT1,POS2=IT2)
159    CALL TRACK_probe_x(R,Y,+STATE,node1=IT1,node2=IT2)
160    U=.NOT.CHECK_STABLE
161    nf=y
162    TA(1)=0.75_dp
163    TA(2)=0.68_dp
164    write(6,*) " arc tunes ",nf%tune(1:2)
165    eq(1)=nf%dhdj%v(1) !-TA(1)
166    eq(2)=nf%dhdj%v(2) !-TA(2)
167    eq(3)=nf%dhdj%v(1)
168
169    !   call print(nf%dhdj%v(1),6)
170    !   call print(nf%dhdj%v(2),6)
171
172    ! nf%dhdj%v(1)=(nf%dhdj%v(1)<=4)
173    ! nf%dhdj%v(2)=(nf%dhdj%v(2)<=4)
174
175    call kanalnummer(mf)
176    open(unit=mf,file='eq.txt')
177    eq(4)=(nf%A_T%V(1).par.'1000')**2+(nf%A_T%V(1).par.'0100')**2
178    eq(5)=(nf%A_T%V(3).par.'0010')**2+(nf%A_T%V(3).par.'0001')**2
179    !call print(eq(4),6)
180    !call print(eq(5),6)
181    ALX=-(nf%A_T%V(1).SUB.'1')*(nf%A_T%V(2).SUB.'1')-(nf%A_T%V(1).SUB.'01')*(nf%A_T%V(2).SUB.'01')
182    ALY=-(nf%A_T%V(3).SUB.'001')*(nf%A_T%V(4).SUB.'001')-(nf%A_T%V(3).SUB.'0001')*(nf%A_T%V(4).SUB.'0001')
183
184    !WRITE(6,*) BEX,BEY
185    !WRITE(6,*) ALX,ALY
186    x=0.0_dp
187    ID=1
188    Y=X+ID
189    !( R,X,U,K,POS1,POS2,T1,T2,P1,P2,IN_P1,IN_P2,POS1_FIBRE,POS2_FIBRE)
190    CALL TRACK_probe_x(R,Y,+STATE,node1=IT3,node2=IT4)
191    U=.NOT.CHECK_STABLE
192
193    !    CALL TRACK_BEAM_x(R,Y,U,+STATE,POS1=IT3,POS2=IT4)
194    nf=y
195    eq(1)=(eq(1)*8 + (1.0_dp+nf%dhdj%v(1)))*3-targ(1)
196    eq(2)=(eq(2)*8 + (1.0_dp+nf%dhdj%v(2)))*3-targ(2)
197    eq(3)=eq(3)-ta(1)
198
199    eq(4)=eq(4)-(nf%A_T%V(1).par.'1000')**2-(nf%A_T%V(1).par.'0100')**2
200    eq(5)=eq(5)-(nf%A_T%V(3).par.'0010')**2-(nf%A_T%V(3).par.'0001')**2
201    !call print(eq(4),6)
202    !call print(eq(5),6)
203
204    do i=1,5
205       eq(i)=eq(i)<=4
206       xx=eq(i)
207       write(6,*) i,xx
208    enddo
209    do i=1,5
210       call print(eq(i),mf)
211    enddo
212
213    close(mf)
214    CALL kill(ID); call kill(nf);call kill(Y);call kill(eq,5);
215
216
217    NP=nu
218    nt=NP+5
219
220    CALL INIT(1,nt)
221    call alloc(g,nt)
222    call alloc(T)
223
224    call kanalnummer(mf)
225    open(unit=mf,file='eq.txt')
226    do i=np+1,nt
227       call read(g%v(i),mf)
228       g%v(i)=g%v(i)-(1.0_dp-sc)*(g%v(i).sub.'0')
229    enddo
230
231    call alloc(t)
232    do i=1,np
233       g%v(i)=1.0_dp.mono.i
234       do j=np+1,nt
235          t=g%v(j).d.i
236          g%v(i)=g%v(i)+(1.0_dp.mono.j)*t
237       enddo
238    enddo
239    CALL KILL(t)
240    write(6,*) "stable ", check_stable
241    g=g.oo.(-1)
242    write(6,*) "stable ", check_stable,nu
243    tpsafit(1:nt)=g
244    write(6,*) tpsafit(1:nu)
245    CALL KILL(G)
246
247    SET_TPSAFIT=.true.
248    DO I=1,NU
249       R=QC(I)
250    ENDDO
251    SET_TPSAFIT=.false.
252    close(mf)
253
254    WRITE(6,*) " MORE "
255    READ(5,*) MF
256    IF(MF==1) GOTO 111
257
258
259    CALL KILL_PARA(R)
260
261    CALL ELP_TO_EL(R)
262
263  end   subroutine special_alex_main_ring
264
265  subroutine special_alex_main_ring_auto(r,n_name,targ,sc,epsa)
266    implicit none
267    TYPE(layout), target, intent(inout):: R
268    integer  i1,i2,it1,it2,it3,it4
269    INTEGER I,N,NU,N2,NP2,mf,nt,NP,J,n_name
270    type(fibre), pointer :: p1,p2
271    TYPE(POL_BLOCK) QC(11)
272    TYPE(REAL_8) Y(6)
273    TYPE(DAMAP) ID
274    REAL(DP) X(6),targ(2),xx,epsa,nxx
275    TYPE(INTERNAL_STATE) STATE
276    LOGICAL(LP) U
277    type(normalform) nf
278    type(gmap) g
279    type(TAYLOR) T,eq(5)
280
281    REAL(DP) ALX,ALY,TA(2),sc
282
283    !targ(1)=22.43d0
284    !targ(2)=20.82d0
285
286    N=10
287    NU=11
288    if(.not.associated(r%t)) then
289       write(6,*) " thin lens lattice not made "
290       CALL MAKE_node_LAYOUT(r)
291    endif
292    p1=>r%start
293    do i=1,r%n
294       if(p1%mag%name(1:3)=='QDX') then
295          i1=i
296          exit
297       endif
298       p1=>p1%next
299    enddo
300    p2=>p1%next
301    do i=i1+1,r%n
302       if(p2%mag%name(1:3)=='QDX') then
303          i2=i
304          exit
305       endif
306       p2=>p2%next
307    enddo
308    !    call move_to(r,p1,i1)
309    !    call move_to(r,p2,i2)
310
311    write(6,*) p1%mag%name,p1%mag%p%nst
312    write(6,*) p2%mag%name,p2%mag%p%nst
313
314    IT1=p1%T1%POS+2 + (p1%mag%p%nst/2 )
315    IT2=p2%T1%POS+2 + (p2%mag%p%nst/2 )
316
317    write(6,*) i1,IT1,i2,IT2
318
319    if(mod(p1%mag%p%nst,2)/=0.or.mod(p2%mag%p%nst,2)/=0) then
320       write(6,*) " Even number of split needed for fitting in Alex_special "
321       stop 100
322    endif
323    !    call move_to(r,p1,i3)
324    !    call move_to(r,p2,i4)
325    p1=>p2%next
326    do i=i2+1,r%n
327       if(p1%mag%name(1:3)=='QFS') then
328          i1=i
329          exit
330       endif
331       p1=>p1%next
332    enddo
333
334    do i=i1,1,-1
335       if(p1%mag%name(1:3)=='QDX') then
336          i1=i
337          exit
338       endif
339       p1=>p1%previous
340    enddo
341
342
343    p2=>p1%next
344    do i=i1+1,r%n
345       if(p2%mag%name(1:3)=='QDX') then
346          i2=i
347          exit
348       endif
349       p2=>p2%next
350    enddo
351
352    write(6,*) p1%mag%name,p1%mag%p%nst
353    write(6,*) p2%mag%name,p2%mag%p%nst
354
355    IT3=p1%T1%POS+2 + (p1%mag%p%nst/2 )
356    IT4=p2%T1%POS+2 + (p2%mag%p%nst/2 )
357
358    write(6,*) i1,IT3,i2,IT4
359    if(mod(p1%mag%p%nst,2)/=0.or.mod(p2%mag%p%nst,2)/=0) then
360       write(6,*) " Even number of split needed for fitting in Alex_special "
361       stop 101
362    endif
363
364    DO I=1,NU
365       QC(I)=0
366       QC(I)%n_name=n_name
367    ENDDO
368    QC(1)%NAME='QFX'
369    QC(2)%NAME='QDX'
370    QC(3)%NAME='QFN'
371    QC(4)%NAME='QDN'
372    QC(5)%NAME='QFS'
373    QC(6)%NAME='QDS'
374    QC(7)%NAME='QFT'
375    QC(8)%NAME='QFP'
376    QC(9)%NAME='QDT'
377    QC(10)%NAME='QFR'
378    QC(11)%NAME='QDR'
379
380    DO I=1,NU
381       QC(I)%IBN(2)=I
382    ENDDO
383
384    DO I=1,NU
385       R=QC(I)
386    ENDDO
387
388111 STATE=DEFAULT0+ONLY_4D0
389
390    X=0.0_dp
391    CALL INIT(STATE,2,NU,BERZ,N2,NP2)
392    CALL ALLOC(ID); call alloc(nf);call alloc(Y);call alloc(eq,5);
393    ID=1
394    Y=X+ID
395    !( R,X,U,K,POS1,POS2,T1,T2,P1,P2,IN_P1,IN_P2,POS1_FIBRE,POS2_FIBRE)
396    !   CALL TRACK_BEAM_x(R,Y,U,+STATE,POS1=IT1,POS2=IT2)
397    CALL TRACK_probe_x(R,Y,+STATE,node1=IT1,node2=IT2)
398    U=.NOT.CHECK_STABLE
399
400    nf=y
401    TA(1)=0.75_dp
402    TA(2)=0.68_dp
403    write(6,*) " arc tunes ",nf%tune(1:2)
404    eq(1)=nf%dhdj%v(1) !-TA(1)
405    eq(2)=nf%dhdj%v(2) !-TA(2)
406    eq(3)=nf%dhdj%v(1)
407
408    !   call print(nf%dhdj%v(1),6)
409    !   call print(nf%dhdj%v(2),6)
410
411    ! nf%dhdj%v(1)=(nf%dhdj%v(1)<=4)
412    ! nf%dhdj%v(2)=(nf%dhdj%v(2)<=4)
413
414    call kanalnummer(mf)
415    open(unit=mf,file='eq.txt')
416    eq(4)=(nf%A_T%V(1).par.'1000')**2+(nf%A_T%V(1).par.'0100')**2
417    eq(5)=(nf%A_T%V(3).par.'0010')**2+(nf%A_T%V(3).par.'0001')**2
418    !call print(eq(4),6)
419    !call print(eq(5),6)
420    ALX=-(nf%A_T%V(1).SUB.'1')*(nf%A_T%V(2).SUB.'1')-(nf%A_T%V(1).SUB.'01')*(nf%A_T%V(2).SUB.'01')
421    ALY=-(nf%A_T%V(3).SUB.'001')*(nf%A_T%V(4).SUB.'001')-(nf%A_T%V(3).SUB.'0001')*(nf%A_T%V(4).SUB.'0001')
422
423    !WRITE(6,*) BEX,BEY
424    !WRITE(6,*) ALX,ALY
425    x=0.0_dp
426    ID=1
427    Y=X+ID
428    !( R,X,U,K,POS1,POS2,T1,T2,P1,P2,IN_P1,IN_P2,POS1_FIBRE,POS2_FIBRE)
429    CALL TRACK_probe_x(R,Y,+STATE,node1=IT3,node2=IT4)
430    U=.NOT.CHECK_STABLE
431
432    !    CALL TRACK_BEAM_x(R,Y,U,+STATE,POS1=IT3,POS2=IT4)
433    nf=y
434    eq(1)=(eq(1)*8 + (1.0_dp+nf%dhdj%v(1)))*3-targ(1)
435    eq(2)=(eq(2)*8 + (1.0_dp+nf%dhdj%v(2)))*3-targ(2)
436    eq(3)=eq(3)-ta(1)
437
438    eq(4)=eq(4)-(nf%A_T%V(1).par.'1000')**2-(nf%A_T%V(1).par.'0100')**2
439    eq(5)=eq(5)-(nf%A_T%V(3).par.'0010')**2-(nf%A_T%V(3).par.'0001')**2
440    !call print(eq(4),6)
441    !call print(eq(5),6)
442    nxx=0.0_dp
443    do i=1,5
444       eq(i)=eq(i)<=4
445       xx=eq(i)
446       write(6,*) i,xx
447       nxx=nxx+abs(xx)
448    enddo
449    do i=1,5
450       call print(eq(i),mf)
451    enddo
452
453    close(mf)
454    CALL kill(ID); call kill(nf);call kill(Y);call kill(eq,5);
455
456
457    NP=nu
458    nt=NP+5
459
460    CALL INIT(1,nt)
461    call alloc(g,nt)
462    call alloc(T)
463
464    call kanalnummer(mf)
465    open(unit=mf,file='eq.txt')
466    do i=np+1,nt
467       call read(g%v(i),mf)
468       g%v(i)=g%v(i)-(1.0_dp-sc)*(g%v(i).sub.'0')
469    enddo
470
471    call alloc(t)
472    do i=1,np
473       g%v(i)=1.0_dp.mono.i
474       do j=np+1,nt
475          t=g%v(j).d.i
476          g%v(i)=g%v(i)+(1.0_dp.mono.j)*t
477       enddo
478    enddo
479    CALL KILL(t)
480    write(6,*) "stable ", check_stable
481    g=g.oo.(-1)
482    write(6,*) "stable ", check_stable,nu
483    tpsafit(1:nt)=g
484    write(6,*) tpsafit(1:nu)
485    CALL KILL(G)
486
487    SET_TPSAFIT=.true.
488    DO I=1,NU
489       R=QC(I)
490    ENDDO
491    SET_TPSAFIT=.false.
492    close(mf)
493    mf=1
494    if(  nxx<=epsa) mf=0
495    !    WRITE(6,*) " MORE "
496    !    READ(5,*) MF
497    IF(MF==1) GOTO 111
498
499
500    CALL KILL_PARA(R)
501
502    CALL ELP_TO_EL(R)
503
504  end   subroutine special_alex_main_ring_auto
505
506
507  subroutine special_alex_main_ring1(r,i1,i2,I3,I4,targ,sc)
508    implicit none
509    TYPE(layout), target, intent(inout):: R
510    integer  i1,i2,I3,I4,it1,it2,it3,it4
511    INTEGER I,N,NU,N2,NP2,mf,nt,NP,J
512    type(fibre), pointer :: p1,p2
513    TYPE(POL_BLOCK) QC(10)
514    TYPE(REAL_8) Y(6)
515    TYPE(DAMAP) ID
516    REAL(DP) X(6),targ(2),xx  ,tas(6)
517    TYPE(INTERNAL_STATE) STATE
518    LOGICAL(LP) U
519    type(normalform) nf
520    type(gmap) g
521    type(TAYLOR) T,eq(6)
522
523    REAL(DP) ALX,ALY,BEX,BEY,NUX,NUY,TA(2),sc
524
525    !targ(1)=22.43d0
526    !targ(2)=20.82d0
527
528    N=10
529    NU=4
530    if(.not.associated(r%t)) then
531       write(6,*) " thin lens lattice not made "
532       stop 300
533    endif
534
535    call move_to(r,p1,i1)
536    call move_to(r,p2,i2)
537
538    write(6,*) p1%mag%name,p1%mag%p%nst
539    write(6,*) p2%mag%name,p2%mag%p%nst
540
541    IT1=p1%T1%POS+2 + (p1%mag%p%nst/2 )
542    IT2=p2%T1%POS+2 + (p2%mag%p%nst/2 )
543
544    write(6,*) IT1,IT2
545
546    call move_to(r,p1,i3)
547    call move_to(r,p2,i4)
548
549    write(6,*) p1%mag%name,p1%mag%p%nst
550    write(6,*) p2%mag%name,p2%mag%p%nst
551
552    IT3=p1%T1%POS+2 + (p1%mag%p%nst/2 )
553    IT4=p2%T1%POS+2 + (p2%mag%p%nst/2 )
554
555    write(6,*) IT3,IT4
556
557    DO I=1,NU
558       QC(I)=0
559    ENDDO
560    QC(1)%NAME='QFX'
561    QC(2)%NAME='QDX'
562    QC(3)%NAME='QFN'
563    QC(4)%NAME='QDN'
564
565    DO I=1,NU
566       QC(I)%IBN(2)=I
567    ENDDO
568
569    DO I=1,NU
570       R=QC(I)
571    ENDDO
572
573111 STATE=DEFAULT0+ONLY_4D0
574
575    X=0.0_dp
576    CALL INIT(STATE,2,NU,BERZ,N2,NP2)
577    CALL ALLOC(ID); call alloc(nf);call alloc(Y);
578    ID=1
579    Y=X+ID
580    !( R,X,U,K,POS1,POS2,T1,T2,P1,P2,IN_P1,IN_P2,POS1_FIBRE,POS2_FIBRE)
581    !    CALL TRACK_BEAM_x(R,Y,U,+STATE,POS1=IT1,POS2=IT2)
582    CALL TRACK_probe_x(R,Y,+STATE,node1=IT1,node2=IT2)
583    U=.NOT.CHECK_STABLE
584
585    nf=y
586    TA(1)=0.75_dp
587    TA(2)=0.68_dp
588    write(6,*) " arc tunes ",nf%tune(1:2)
589    nf%dhdj%v(1)=nf%dhdj%v(1)-TA(1)
590    nf%dhdj%v(2)=nf%dhdj%v(2)-TA(2)
591    call print(nf%dhdj%v(1),6)
592    call print(nf%dhdj%v(2),6)
593
594    nf%dhdj%v(1)=(nf%dhdj%v(1)<=4)
595    nf%dhdj%v(2)=(nf%dhdj%v(2)<=4)
596
597    call kanalnummer(mf)
598    open(unit=mf,file='eq.txt')
599    call print(nf%dhdj%v(1),mf)
600    call print(nf%dhdj%v(2),mf)
601    BEX=(nf%A_T%V(1).SUB.'1')**2+(nf%A_T%V(1).SUB.'01')**2
602    BEY=(nf%A_T%V(3).SUB.'001')**2+(nf%A_T%V(3).SUB.'0001')**2
603    ALX=-(nf%A_T%V(1).SUB.'1')*(nf%A_T%V(2).SUB.'1')-(nf%A_T%V(1).SUB.'01')*(nf%A_T%V(2).SUB.'01')
604    ALY=-(nf%A_T%V(3).SUB.'001')*(nf%A_T%V(4).SUB.'001')-(nf%A_T%V(3).SUB.'0001')*(nf%A_T%V(4).SUB.'0001')
605
606    WRITE(6,*) BEX,BEY
607    WRITE(6,*) ALX,ALY
608    tas(1:2)=TARG(1:2)
609    tas(3)=BEX
610    tas(4)=BEY
611    tas(5)=ALX
612    tas(6)=ALY
613    close(mf)
614    CALL kill(ID); call kill(nf);call kill(Y);
615
616    NP=nu
617    nt=NP+2
618
619    CALL INIT(1,nt)
620    call alloc(g,nt)
621    call alloc(T)
622
623    call kanalnummer(mf)
624    open(unit=mf,file='eq.txt')
625    do i=np+1,nt
626       call read(g%v(i),mf)
627       g%v(i)=g%v(i)-(1.0_dp-sc)*(g%v(i).sub.'0')
628    enddo
629
630
631    call alloc(t)
632    do i=1,np
633       g%v(i)=1.0_dp.mono.i
634       do j=np+1,nt
635          t=g%v(j).d.i
636          g%v(i)=g%v(i)+(1.0_dp.mono.j)*t
637       enddo
638    enddo
639    CALL KILL(t)
640
641    write(6,*) "stable ", check_stable
642    g=g.oo.(-1)
643    write(6,*) "stable ", check_stable,nu
644    tpsafit(1:nt)=g
645    CALL KILL(G)
646    SET_TPSAFIT=.true.
647    DO I=1,NU
648       R=QC(I)
649    ENDDO
650    SET_TPSAFIT=.false.
651    close(mf)
652
653    WRITE(6,*) " MORE "
654    READ(5,*) MF
655    IF(MF==1) GOTO 111
656    CALL KILL_PARA(R)
657
658    NU=7
659    DO I=1,NU
660       QC(I)=0
661    ENDDO
662    QC(1)%NAME='QFS'
663    QC(2)%NAME='QDS'
664    QC(3)%NAME='QFT'
665    QC(4)%NAME='QFP'
666    QC(5)%NAME='QDT'
667    QC(6)%NAME='QFR'
668    QC(7)%NAME='QDR'
669
670    DO I=1,NU
671       QC(I)%IBN(2)=I
672    ENDDO
673
674    DO I=1,NU
675       R=QC(I)
676    ENDDO
677
678112 continue
679    X=0.0_dp
680    CALL INIT(STATE,2,NU,BERZ,N2,NP2)
681    CALL ALLOC(ID); call alloc(nf);call alloc(Y);
682    ID=1
683    Y=X+ID
684    !( R,X,U,K,POS1,POS2,T1,T2,P1,P2,IN_P1,IN_P2,POS1_FIBRE,POS2_FIBRE)
685    !    CALL TRACK_BEAM_x(R,Y,U,+STATE,POS1=IT3,POS2=IT4)
686    CALL TRACK_probe_x(R,Y,+STATE,node1=IT3,node2=IT4)
687    U=.NOT.CHECK_STABLE
688
689    write(6,*) "before stable nf=y", check_stable
690    !    global_verbose=.true.
691    nf=y
692    !   nf%dhdj%v(1)=nf%dhdj%v(1)-0.75_DP
693    !   nf%dhdj%v(2)=nf%dhdj%v(2)-0.68_DP
694    write(6,*) "stable nf=y", check_stable
695
696
697
698    call alloc(eq,4)
699
700    eq(1)=(TA(1)*8 + (1.0_dp+nf%dhdj%v(1)))*3
701    eq(2)=(TA(2)*8 + (1.0_dp+nf%dhdj%v(2)))*3
702
703    nux=eq(1)
704    nuy=eq(2)
705
706    eq(1)=eq(1)-targ(1)
707    eq(2)=eq(2)-targ(2)
708
709    WRITE(6,*) NUX,NUY
710    WRITE(6,*) NUX-targ(1),NUY-targ(2)
711
712    eq(3)=(nf%A_T%V(1).par.'1000')**2+(nf%A_T%V(1).par.'0100')**2-bex
713    eq(4)=(nf%A_T%V(3).par.'0010')**2+(nf%A_T%V(3).par.'0001')**2-bey
714    !eq(5)=-alx-(nf%A_T%V(1).par.'1000')*(nf%A_T%V(2).par.'1')-(nf%A_T%V(1).par.'0100')*(nf%A_T%V(2).par.'0100')
715    !eq(6)=-aly-(nf%A_T%V(3).par.'0010')*(nf%A_T%V(4).par.'0010')-(nf%A_T%V(3).par.'0001')*(nf%A_T%V(4).par.'0001')
716
717    do i=1,4
718       eq(i)=(eq(i)<=4)
719       xx=eq(i)
720       write(6,*) i,xx,TAS(I)
721    enddo
722    call kanalnummer(mf)
723    open(unit=mf,file='eq.txt')
724    do i=1,4
725       call print(eq(i),mf)
726    enddo
727    close(mf)
728    CALL kill(ID); call kill(nf);call kill(Y);
729    call kill(eq,4)
730
731    NP=nu
732    nt=NP+4
733
734    CALL INIT(1,nt)
735    call alloc(g,nt)
736    call alloc(T)
737
738    call kanalnummer(mf)
739    open(unit=mf,file='eq.txt')
740    do i=np+1,nt
741       call read(g%v(i),mf)
742       g%v(i)=g%v(i)-(1.0_dp-sc)*(g%v(i).sub.'0')
743    enddo
744
745    call alloc(t)
746    do i=1,np
747       g%v(i)=1.0_dp.mono.i
748       do j=np+1,nt
749          t=g%v(j).d.i
750          g%v(i)=g%v(i)+(1.0_dp.mono.j)*t
751       enddo
752    enddo
753    CALL KILL(t)
754    write(6,*) "stable ", check_stable
755    g=g.oo.(-1)
756    write(6,*) "stable ", check_stable,nu
757    tpsafit(1:nt)=g
758    write(6,*) tpsafit(1:nu)
759    CALL KILL(G)
760
761    SET_TPSAFIT=.true.
762    DO I=1,NU
763       R=QC(I)
764    ENDDO
765    SET_TPSAFIT=.false.
766    close(mf)
767
768    WRITE(6,*) " MORE "
769    READ(5,*) MF
770    IF(MF==1) GOTO 112
771
772
773    CALL ELP_TO_EL(R)
774
775    CALL KILL_PARA(R)
776
777  end   subroutine special_alex_main_ring1
778
779  subroutine lattice_linear_res_gmap(R,my_state,EPSF,POLY,NPOLY,TARG,NP)
780    IMPLICIT NONE
781    TYPE(layout), intent(inout):: R
782    TYPE(POL_BLOCK), intent(inout),dimension(:)::POLY
783    INTEGER, intent(in):: NPOLY,NP
784    real(dp) , intent(IN),dimension(:)::TARG
785    real(dp) CLOSED(6)
786    TYPE(INTERNAL_STATE), intent(IN):: my_STATE
787    TYPE(INTERNAL_STATE) STATE
788    INTEGER I,SCRATCHFILE
789    TYPE(TAYLOR), allocatable:: EQ(:)
790    TYPE(REAL_8) Y(6)
791    TYPE(NORMALFORM) NORM
792    integer :: neq=6, no=2,nt,j,it
793    type(damap) id
794    type(gmap) g
795    TYPE(TAYLOR)t
796    real(dp) epsf,epsr,epsnow,tune(2),co(4)
797    !    EPSF=.0001
798    epsr=abs(epsf)
799
800    allocate(eq(neq))
801
802    nt=neq+np
803    STATE=((((my_state+nocavity0))+only_4d0)-RADIATION0)
804
805    CALL INIT(STATE,no,NP,BERZ)
806
807    SET_TPSAFIT=.FALSE.
808
809    DO I=1,NPOLY
810       R=POLY(i)
811    ENDDO
812    CLOSED(:)=0.0_dp
813    it=0
814100 continue
815    it=it+1
816
817    CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
818    write(6,*) "closed orbit ", CHECK_STABLE
819    write(6,*) CLOSED
820
821
822    CALL INIT(STATE,no,NP,BERZ)
823    CALL ALLOC(NORM)
824    CALL ALLOC(Y)
825    CALL ALLOC(EQ)
826    call alloc(id)
827
828    id=1
829    Y=CLOSED+id
830
831    CALL TRACK(R,Y,1,+STATE)
832    NORM=Y
833    write(6,*) " tunes ",NORM%TUNE(1), NORM%TUNE(2), CHECK_STABLE
834    tune(1)=(NORM%dhdj%v(1)).SUB.'0000'
835    tune(2)=(NORM%dhdj%v(2)).SUB.'0000'
836    co(1)=(y(1).SUB.'0010')
837    co(2)=(y(1).SUB.'0001')
838    co(3)=(y(2).SUB.'0010')
839    co(4)=(y(2).SUB.'0001')
840    id=y
841    id=NORM%a1**(-1)*id*NORM%a1
842
843    eq(1)=       ((NORM%dhdj%v(1)).par.'00000')-targ(1)
844    eq(2)=       ((NORM%dhdj%v(2)).par.'00000')-targ(2)
845    eq(3)=       (id%v(1).par.'0010') !-targ(3)
846    eq(4)=       (id%v(1).par.'0001') !-targ(4)
847    eq(5)=       (id%v(2).par.'0010') !-targ(5)
848    eq(6)=       (id%v(2).par.'0001') !-targ(6)
849    epsnow=abs(eq(1))+abs(eq(2))+abs(eq(3))+abs(eq(4))+abs(eq(5))+abs(eq(6))
850    write(6,*) " closeness ",epsnow
851    call kanalnummer(SCRATCHFILE)
852    OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
853    rewind scratchfile
854
855    do i=1,neq
856       eq(i)=eq(i)<=c_%npara
857    enddo
858    do i=1,neq
859       call daprint(eq(i),scratchfile)
860    enddo
861    close(SCRATCHFILE)
862    CALL KILL(NORM)
863    CALL KILL(Y)
864    CALL KILL(id)
865    CALL KILL(EQ)
866
867
868
869    CALL INIT(1,nt)
870    call alloc(g,nt)
871    call kanalnummer(SCRATCHFILE)
872    OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
873    rewind scratchfile
874    do i=np+1,nt
875       call read(g%v(i),scratchfile)
876    enddo
877    close(SCRATCHFILE)
878
879    call alloc(t)
880    do i=1,np
881       g%v(i)=1.0_dp.mono.i
882       do j=np+1,nt
883          t=g%v(j).d.i
884          g%v(i)=g%v(i)+(1.0_dp.mono.j)*t
885       enddo
886    enddo
887    CALL KILL(t)
888
889    g=g.oo.(-1)
890    tpsafit(1:nt)=g
891
892    SET_TPSAFIT=.true.
893
894    DO I=1,NPOLY
895       R=POLY(i)
896    ENDDO
897    SET_TPSAFIT=.false.
898
899    CALL ELP_TO_EL(R)
900
901    !    write(6,*) " more "
902    !    read(5,*) more
903    if(it>=max_fit_iter) goto 101
904    if(epsnow<=epsr) goto 102
905    GOTO 100
906
907101 continue
908    write(6,*) " warning did not converge "
909
910102 continue
911    CALL KILL_PARA(R)
912    deallocate(eq)
913
914  end subroutine lattice_linear_res_gmap
915
916
917  subroutine special_alex_main_ring_removal(main)
918    IMPLICIT NONE
919
920!!!!!!!!!
921    type(layout),pointer :: main
922    type(fibre),pointer :: p,p1,p2,pp
923    real(dp) an,bn
924    logical doit
925    integer ij,mf,i,k,jk
926
927    call kanalnummer(mf,"diag.txt")
928
929
930    if(associated(main%t)) call kill(main%t)
931
932    ij=0
933    jk=0
934    p=>main%start
935    do i=1,main%n
936       p1=>p%next
937       p2=>p1%next
938       doit=p%mag%name(1:len_trim(p%mag%name)-1)==p1%mag%name(1:len_trim(p1%mag%name)-1)
939       doit=doit.and.(p%mag%name(1:len_trim(p%mag%name)-1)==p2%mag%name(1:len_trim(p1%mag%name)-1))
940       if(doit) then
941          ij=ij+1
942          do k=p%mag%p%nmul,1,-1
943             an=p%mag%an(k)/p1%mag%l*2.0_dp
944             bn=p%mag%bn(k)/p1%mag%l*2.0_dp
945             call add(p1,-k,1,an)
946             call add(p1,k,1,bn)
947             !   call add(p,-k,0,zero)
948             !   call add(p,k,0,zero)
949             !   call add(p2,-k,0,zero)
950             !   call add(p2,k,0,zero)
951          enddo
952          ! el%k3%thin_h_foc,el%k3%thin_v_foc,el%k3%thin_h_angle,el%k3%thin_v_angle
953          if(p%mag%kind==kind3) then
954             bn=p%mag%k3%thin_h_angle/p1%mag%l*2.0_dp
955             call add(p1,1,1,bn)
956             an=p%mag%k3%thin_v_angle/p1%mag%l*2.0_dp
957             call add(p1,-1,1,an)
958             doit=(p%mag%k3%thin_h_foc/=0.0_dp.or.p%mag%k3%thin_v_foc/=0.0_dp)  !.or.p%mag%k3%thin_v_angle/=zero)
959             if(doit) then
960                write(6,*) " cannot handle the stuff in kind3 related to fake thinlens tracking of MAD8 "
961                stop 8
962             endif
963          endif
964
965          write(mf,*) p%mag%name,p1%mag%name,p2%mag%name
966          !  goto 200
967          pp=>p%previous
968          pp%next=>p1
969          p1%previous=>pp
970          pp=>p2%next
971          pp%previous=>p1
972          p1%next=>pp
973
974          call super_kill(p)
975          call super_kill(p2)
976          p=>p1
977          !   200 continue
978          !   p=>p2
979       else
980          if(p%mag%kind==kind3) then
981             jk=jk+1
982             bn=p%mag%k3%thin_h_angle
983             call add(p,1,1,bn)
984             an=p%mag%k3%thin_v_angle
985             call add(p,-1,1,an)
986             p%mag%k3%thin_v_angle=0.0_dp
987             p%mag%k3%thin_h_angle=0.0_dp
988             p%magp%k3%thin_v_angle=0.0_dp
989             p%magp%k3%thin_h_angle=0.0_dp
990             doit=(p%mag%k3%thin_h_foc/=0.0_dp.or.p%mag%k3%thin_v_foc/=0.0_dp)
991             if(doit) then
992                write(6,*) " cannot handle the stuff in kind3 related to fake thinlens tracking of MAD8 "
993                stop 9
994             endif
995          endif
996       endif
997       p=>p%next
998       if(associated(p,main%start)) exit
999
1000
1001    enddo
1002
1003    write(mf,*) ij, " magnet sandwiches done"
1004    write(mf,*) jk, " single magnets done"
1005
1006    ij=0
1007    p=>main%start
1008    do i=1,main%n
1009       ij=ij+1
1010       p%pos=i
1011       if(p%mag%l==0.0_dp.and.p%mag%kind>kind1) then
1012          write(mf,*) i,p%mag%name,p%mag%kind,p%next%mag%name
1013       endif
1014       p=>p%next
1015       if(associated(p,main%start)) exit
1016    enddo
1017    write(6,*) ij,main%n
1018    main%n=ij
1019
1020    close(mf)
1021  end subroutine special_alex_main_ring_removal
1022
1023!!!!!!!!!!!!!!!!!!!!!!!!!!!!  real stuff   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1024  subroutine  alex_track_monitors(r,x,state) ! one particle tracking m_turn
1025    implicit none
1026    type(layout), target :: r
1027    integer nturn,i,j,nm
1028    real(dp) x(6),y(6)
1029    type(fibre),pointer::p1
1030    type(internal_state) state
1031
1032    p1=>r%start
1033
1034    nturn=SIZE (monitors(1)%r,dim=2)
1035    nm=size(monitors)
1036    y=x
1037    CALL FIND_ORBIT(R,y,1,STATE,1e-5_dp)
1038    write(6,*)" closed orbit "
1039    write(6,*) y(1:3)
1040    write(6,*) y(4:6)
1041    x(1:4)=x(1:4)+y(1:4)
1042    do i=1,nturn
1043       call track(x,state,p1,monitors(1)%p)
1044
1045       call track(y,state,p1,monitors(1)%p)
1046       monitors(1)%r(1:4,i)=x(1:4)-y(1:4)
1047       monitors(1)%bpm(1,i)=monitors(1)%r(1,i)
1048       monitors(1)%bpm(2,i)=monitors(1)%r(3,i)
1049
1050       do j=1,nm-1
1051          call track(x,state,monitors(j)%p,monitors(j+1)%p)
1052          call track(y,state,monitors(j)%p,monitors(j+1)%p)
1053          monitors(j+1)%r(1:4,i)=x(1:4)-y(1:4)
1054          monitors(j+1)%bpm(1,i)=monitors(j)%r(1,i)
1055          monitors(j+1)%bpm(2,i)=monitors(j)%r(3,i)
1056
1057       enddo
1058       call track(x,state,monitors(nm)%p,p1)
1059       call track(y,state,monitors(nm)%p,p1)
1060
1061    enddo
1062
1063
1064
1065
1066
1067  end subroutine  alex_track_monitors
1068
1069
1070
1071  subroutine invert_monitors(kind,zf0,g)
1072    implicit none
1073    type(damap) ex
1074    type(gmap) gi,g
1075    real(dp) zf(4),zx(6),zf0(4)
1076    integer kind,i
1077
1078    zf=zf0
1079
1080    call alloc(ex)
1081    call alloc(gi)
1082    ! zf contains measuments xf(x_1,x_11,y_2,y_22)
1083    gi=g
1084
1085    if(kind==1) then
1086       !write(6,*) zf(1:4)
1087       zx=0.0_dp
1088       zx(1)=zf(1)
1089       zx(2)=zf(2)   !!!
1090       ex=0
1091       ex=zx
1092       ex%v(3)=ex%v(3)+(1.0_dp.mono.3)
1093       ex%v(4)=ex%v(4)+(1.0_dp.mono.4)
1094       gi%v(2)=(gi%v(2).o.ex)-(1.0_dp.mono.2)
1095       zx=0.0_dp
1096       zx(1)=zf(1)
1097       zx(3)=zf(3)   !!!
1098       ex=0
1099       ex=zx
1100       ex%v(2)=ex%v(2)+(1.0_dp.mono.2)
1101       ex%v(4)=ex%v(4)+(1.0_dp.mono.4)
1102
1103       gi%v(3)=(gi%v(3).o.ex)-(1.0_dp.mono.3)
1104       zx=0.0_dp
1105       zx(1)=zf(1)
1106       zx(4)=zf(4)   !!!
1107       ex=0
1108       ex=zx
1109       ex%v(2)=ex%v(2)+(1.0_dp.mono.2)
1110       ex%v(3)=ex%v(3)+(1.0_dp.mono.3)
1111       gi%v(4)=(gi%v(4).o.ex)-(1.0_dp.mono.4)
1112       !call print(gi,6)
1113       gi%v(1)=gi%v(1)-zf(1)
1114
1115    else
1116       !write(6,*) zf(1:4)
1117
1118       zx=0.0_dp
1119       zx(3)=zf(3)
1120       zx(4)=zf(4)   !!!
1121       ex=0
1122       ex=zx
1123       ex%v(1)=ex%v(1)+(1.0_dp.mono.1)
1124       ex%v(2)=ex%v(2)+(1.0_dp.mono.2)
1125       gi%v(4)=(gi%v(4).o.ex)-(1.0_dp.mono.4)
1126       zx=0.0_dp
1127       zx(1)=zf(1)
1128       zx(3)=zf(3)   !!!
1129       ex=0
1130       ex=zx
1131       ex%v(2)=ex%v(2)+(1.0_dp.mono.2)
1132       ex%v(4)=ex%v(4)+(1.0_dp.mono.4)
1133
1134       gi%v(1)=(gi%v(1).o.ex)-(1.0_dp.mono.1)
1135       zx=0.0_dp
1136       zx(3)=zf(3)
1137       zx(2)=zf(2)   !!!
1138       ex=0
1139       ex=zx
1140       ex%v(4)=ex%v(4)+(1.0_dp.mono.4)
1141       ex%v(1)=ex%v(1)+(1.0_dp.mono.1)
1142       gi%v(2)=(gi%v(2).o.ex)-(1.0_dp.mono.2)
1143       !call print(gi,6)
1144
1145       gi%v(3)=gi%v(3)-zf(3)
1146    endif
1147
1148    gi=gi.oo.(-1)
1149
1150
1151    zf0=gi
1152
1153    call kill(ex)
1154    call kill(gi)
1155
1156  end subroutine invert_monitors
1157
1158
1159  subroutine  alex_mom_monitors(jm,nt)
1160    implicit none
1161    integer nturn,i,j,nm,k,jm,nt
1162    real(dp), pointer :: mom(:,:),r(:,:),a(:,:),at(:,:)
1163    real(dp)ar, momt(6,6),kicke(3),a6(6,6),ai6(6,6),br(6,6)
1164    type(damap) id,m12,ex
1165    type(pbfield) h
1166    type(normalform) norm
1167
1168    CALL INIT(2,2,0,0)    ! fpp
1169    call alloc(id,m12,ex)
1170    call alloc(norm)
1171    call alloc(h)
1172
1173
1174    nturn=SIZE (monitors(1)%r,dim=2)
1175    nm=size(monitors)
1176
1177    !do jm=1,nm
1178
1179    mom=>monitors(jm)%mom
1180    r=>monitors(jm)%xf
1181    a=>monitors(jm)%a
1182    at=>monitors(jm)%at
1183    mom=0.0_dp
1184    do i=1,4
1185       do j=i,4
1186          do k=1,nt
1187             mom(i,j)=r(i,k)*r(j,k)+mom(i,j)
1188          enddo
1189       enddo
1190    enddo
1191
1192    mom=mom/nt
1193    do i=1,4
1194       do j=i,4
1195          mom(j,i)=mom(i,j)
1196       enddo
1197    enddo
1198
1199    i=0
1200
1201    if(i==0) then
1202
1203       momt=0.0_dp
1204       do i=1,2
1205          do j=1,2
1206             if(with_c==0.and.i/=j) cycle
1207             ! momt(2*i-1,2*j-1)=mom(2*i-1,2*j-1)
1208             ! momt(2*i-1,2*j)=mom(2*i-1,2*j)
1209             ! momt(2*i,2*j-1)=mom(2*i,2*j-1)
1210             ! momt(2*i,2*j)=mom(2*i,2*j)
1211             momt(2*i,2*j)=mom(2*i-1,2*j-1)
1212             momt(2*i,2*j-1)=-mom(2*i-1,2*j)
1213             momt(2*i-1,2*j)=-mom(2*i,2*j-1)
1214             momt(2*i-1,2*j-1)=mom(2*i,2*j)
1215          enddo
1216       enddo
1217
1218       momt(5,5)=mom(1,1)
1219       momt(6,6)=mom(1,1)
1220
1221       call diagonalise_envelope_a(momt,br,a6,ai6,kicke)
1222       a=a6(1:4,1:4)
1223       at=matmul(at,a)
1224    else
1225       h%h=0.0_dp
1226
1227       do i=1,2
1228          do j=1,2
1229             if(with_c==0.and.i/=j) cycle
1230             h%h=h%h + mom(2*i-1,2*j-1)*(1.0_dp.mono.(2*i))*(1.0_dp.mono.(2*j))
1231             h%h=h%h + mom(2*i,2*j)*(1.0_dp.mono.(2*i-1))*(1.0_dp.mono.(2*j-1))
1232             h%h=h%h - mom(2*i,2*j-1)*(1.0_dp.mono.(2*i-1))*(1.0_dp.mono.(2*j))
1233             h%h=h%h - mom(2*i-1,2*j)*(1.0_dp.mono.(2*i))*(1.0_dp.mono.(2*j-1))
1234          enddo
1235       enddo
1236
1237       h%h=-h%h
1238       ar=full_abs(h%h)
1239       ar=ar*20.0_dp
1240       h%h=h%h/ar
1241
1242       id=1
1243       ex=texp(h,id)
1244       norm=ex
1245
1246       a=norm%a_t
1247       at=matmul(at,a)
1248       !enddo   ! jm=1,nm
1249
1250    endif
1251
1252    !write(6,*) "betax", a(1,1)**2+a(1,2)**2, sqrt(a(1,1)**2+a(1,2)**2)
1253    !write(6,*) "alphax", -a(1,1)*a(2,1)-a(1,2)*a(2,2)
1254    !write(6,*) "betay", a(3,3)**2+a(3,4)**2, sqrt(a(3,3)**2+a(3,4)**2)
1255    !write(6,*) "alphay", -a(3,3)*a(4,3)-a(3,4)*a(4,4)
1256    !write(6,*) "betaxy", a(1,3)**2+a(1,4)**2
1257    !write(6,*) "betayx", a(3,1)**2+a(3,2)**2
1258
1259
1260    call kill(id,m12,ex)
1261    call kill(norm)
1262    call kill(h)
1263
1264  end subroutine  alex_mom_monitors
1265
1266  subroutine  alex_apply_A_on_data(jm)
1267    implicit none
1268    integer jm,ier,nm
1269    real(dp), pointer :: mom(:,:),r(:,:),a(:,:),xn(:,:)
1270    real(dp) ai(4,4)
1271
1272
1273
1274    nm=size(monitors)
1275
1276    !do jm=1,nm
1277
1278    r=>monitors(jm)%xf
1279    xn=>monitors(jm)%xn
1280    a=>monitors(jm)%a
1281    call matinv(a,ai,4,4,ier)
1282    if(ier/=0) then
1283       write(6,*) " error in matinv ",jm
1284       stop 999
1285    endif
1286
1287    xn=matmul(ai,r)
1288
1289
1290
1291    !enddo   ! jm=1,nm
1292
1293
1294
1295  end subroutine  alex_apply_A_on_data
1296
1297  subroutine  alex_count_monitors(r,n,kind) !locate monitors and create array monitors(n)
1298    implicit none
1299    type(layout), target :: r
1300    type(fibre), pointer :: p
1301    integer i,n,nturn
1302    integer, optional :: kind
1303!!! find qf and qds
1304    !qf are horizontal monitors
1305    !qd are vertical monitors
1306    !!
1307
1308    p=>r%start
1309    n=0
1310    do i=1,r%n
1311
1312       if(p%mag%name(1:2)=="QF") then
1313          n=n+1
1314          !write(6,*) n,p%mag%name, "     x"
1315       endif
1316       if(p%mag%name(1:2)=="QD") then
1317          n=n+1
1318          !write(6,*) n,p%mag%name, "     y"
1319       endif
1320
1321       p=>p%next
1322    enddo
1323
1324    allocate( monitors(n))
1325
1326    p=>r%start
1327    n=0
1328    do i=1,r%n
1329
1330       if(p%mag%name(1:2)=="QF") then
1331          n=n+1
1332          call alloc_fibre_monitor_data(monitors(n),m_turn,p)  ! m_)turn is number of turns
1333          if(present(kind)) then
1334             monitors(n)%kind=kind
1335          else
1336             monitors(n)%kind=1
1337          endif
1338          ! write(6,*) n,p%mag%name
1339       endif
1340       if(p%mag%name(1:2)=="QD") then
1341          n=n+1
1342          call alloc_fibre_monitor_data(monitors(n),m_turn,p)
1343          if(present(kind)) then
1344             monitors(n)%kind=kind
1345          else
1346             monitors(n)%kind=2
1347          endif
1348          ! write(6,*) n,p%mag%name
1349       endif
1350
1351       p=>p%next
1352    enddo
1353
1354  end subroutine alex_count_monitors
1355
1356  subroutine  alex_count_jparc_monitors(r,n,kind) !locate monitors and create array monitors(n)
1357    implicit none
1358    type(layout), target :: r
1359    type(fibre), pointer :: p
1360    integer i,n,nturn,mf,i1,i2,j,i3,kind
1361    integer, allocatable :: ind(:)
1362    character(255) line
1363!!! find qf and qds
1364    !qf are horizontal monitors
1365    !qd are vertical monitors
1366    !!
1367
1368    call kanalnummer(mf,"J_parc_monitor_id.txt")
1369
1370    n=0
1371    do i=1,r%n
1372       read(mf,*,end=100) i1,i2
1373
1374       n=n+1
1375    enddo
1376100 write(6,*) n
1377    allocate(ind(n))
1378    ind=0
1379
1380    rewind mf
1381
1382    do i=1,n
1383       read(mf,*) i1,i2
1384       ind(i1)=i2
1385    enddo
1386
1387
1388    close(mf)
1389
1390
1391    allocate( monitors(n))
1392
1393    p=>r%start
1394    do i=1,r%n
1395
1396       if(p%mag%name(1:2)=="QF") then
1397
1398
1399          i1=0
1400          line=p%mag%name(4:4)
1401          read(line,*) i2
1402          i1=i2*100+i1
1403          line=p%mag%name(5:5)
1404          read(line,*) i2
1405          i1=i2*10+i1
1406          line=p%mag%name(6:6)
1407          read(line,*) i2
1408          i1=i2+i1
1409          i3=i1
1410          if(i1>n) then
1411             i3=n
1412          endif
1413          i2=0
1414          do j=i3,1,-1
1415             if(ind(j)==i1) then
1416                i2=ind(j)
1417                exit
1418             endif
1419          enddo
1420          if(i2/=0) then
1421             call alloc_fibre_monitor_data(monitors(j),m_turn,p)
1422             monitors(j)%kind=kind
1423          endif
1424       endif
1425
1426       if(p%mag%name(1:2)=="QD") then
1427
1428          i1=0
1429          line=p%mag%name(4:4)
1430          read(line,*) i2
1431          i1=i2*100+i1
1432          line=p%mag%name(5:5)
1433          read(line,*) i2
1434          i1=i2*10+i1
1435          line=p%mag%name(6:6)
1436          read(line,*) i2
1437          i1=i2+i1
1438
1439          i3=i1
1440          if(i1>n) then
1441             i3=n
1442          endif
1443          i2=0
1444          do j=i3,1,-1
1445             if(ind(j)==i1) then
1446                i2=ind(j)
1447                exit
1448             endif
1449          enddo
1450
1451          if(i2/=0) then
1452             call alloc_fibre_monitor_data(monitors(j),m_turn,p)
1453             monitors(j)%kind=kind
1454          endif
1455
1456
1457       endif
1458
1459
1460       p=>p%next
1461    enddo
1462
1463
1464
1465    deallocate(ind)
1466
1467  end subroutine alex_count_jparc_monitors
1468
1469  subroutine  alex_read_r_jparc(filename1,filename2,n_max,nav)
1470    implicit none
1471    integer jm,mf1,mf2,n,i,na,i1,i2
1472    CHARACTER(*) filename1,filename2
1473    integer, optional :: n_max,nav
1474    real(dp) x,y,jj(4),xa(2),ya(2),norm(2)
1475    real(dp), allocatable :: bpmx(:), bpmy(:)
1476
1477    allocate(bpmx(size(monitors)),bpmy(size(monitors)))
1478
1479    call kanalnummer(mf1,filename1)
1480    call kanalnummer(mf2,filename2)
1481
1482    n=m_turn
1483    na=0
1484    if(present(n_max)) n=n_max
1485    if(present(nav)) na=nav
1486
1487
1488    x=0.0_dp
1489    y=0.0_dp
1490    do i=1,m_turn+m_skip
1491       read(mf1,*)i1,i2,bpmx
1492       read(mf2,*)i1,i2,bpmy
1493       if(i>m_skip) then
1494          do jm=1,size(monitors)
1495             monitors(jm)%r(1,i-m_skip)=bpmx(jm)/1.d3
1496             monitors(jm)%r(3,i-m_skip)=bpmy(jm)/1.d3
1497             monitors(jm)%bpm(1,i-m_skip)=bpmx(jm)/1.d3
1498             monitors(jm)%bpm(2,i-m_skip)=bpmy(jm)/1.d3
1499          enddo
1500       endif
1501    enddo
1502
1503write(6,*) size(monitors), "monitors "
1504    do jm=1,size(monitors)
1505       norm=0.0_dp
1506       do i=1,m_turn
1507          norm(1)=norm(1)+abs(monitors(jm)%r(1,i))
1508          norm(2)=norm(2)+abs(monitors(jm)%r(3,i))
1509       enddo
1510       if(norm(1)<1.e-10.or. norm(2)<1.e-10) then
1511          write(6,*) " monitor ",jm, " has no data "
1512          monitors(jm)%full=.false.
1513       endif
1514    enddo
1515
1516    close(mf1)
1517    close(mf2)
1518    deallocate(bpmx,bpmy)
1519    ! 171683 59 68 1 2.6794 -1.1248
1520
1521  end   subroutine  alex_read_r_jparc
1522
1523
1524  subroutine scale_bpm
1525    implicit none
1526    integer i,n,ib
1527    real x,y,x0,y0,sx,sy,ax,ay
1528
1529    if(.not.allocated(monitors)) return
1530
1531    do ib=1,size(monitors)
1532       n=m_turn/10
1533
1534
1535       do i=1,n
1536          x0=abs(monitors(ib)%r(1,i))+x0
1537          y0=abs(monitors(ib)%r(3,i))+y0
1538       enddo
1539
1540       do i=m_turn-n+1,m_turn
1541          x=abs(monitors(ib)%r(1,i))+x
1542          y=abs(monitors(ib)%r(3,i))+y
1543       enddo
1544
1545       sx=(x0/x-1)/(m_turn-1)
1546       ax=1-sx
1547       sy=(y0/y-1)/(m_turn-1)
1548       ay=1-sy
1549
1550       do i=1,m_turn
1551          monitors(ib)%r(1,i)=monitors(ib)%r(1,i)*(i*sx+ax)
1552          monitors(ib)%r(3,i)=monitors(ib)%r(3,i)*(i*sy+ay)
1553       enddo
1554
1555    enddo
1556
1557  end subroutine scale_bpm
1558
1559
1560
1561  subroutine  alex_average_r_jparc
1562    implicit none
1563    integer jm,mf1,mf2,n,i,na,i1,i2
1564    real(dp) x,y,jj(4),xa(2),ya(2)
1565
1566
1567    n=m_turn
1568
1569    do jm=1,size(monitors)
1570       monitors(jm)%r(1,:)=monitors(jm)%bpm(1,:)
1571       monitors(jm)%r(3,:)=monitors(jm)%bpm(2,:)
1572    enddo
1573    do jm=1,size(monitors)
1574
1575       x=0.0_dp
1576       y=0.0_dp
1577       do i=1,n
1578          x=monitors(jm)%r(1,i)+x
1579          y=monitors(jm)%r(3,i)+y
1580       enddo
1581
1582       x=x/n
1583       y=y/n
1584       write(6,*) "Averages ",x,y
1585       do i=1,n
1586          monitors(jm)%r(1,i)=(monitors(jm)%r(1,i)-x)
1587          monitors(jm)%r(3,i)=(monitors(jm)%r(3,i)-y)
1588       enddo
1589
1590    enddo
1591
1592
1593    ! 171683 59 68 1 2.6794 -1.1248
1594
1595  end   subroutine  alex_average_r_jparc
1596
1597
1598
1599  subroutine alloc_fibre_monitor_data(c,turn,p)
1600    implicit none
1601    TYPE(fibre_monitor_data), target :: c
1602    TYPE(fibre), target :: p
1603    integer turn,i
1604
1605    nullify(c%p)
1606    c%p=>p
1607    allocate(c%turn)
1608    allocate(c%kind)
1609    allocate(c%r(4,turn))
1610    allocate(c%bpm(2,turn))
1611    allocate(c%xf(4,turn))
1612    allocate(c%xn(4,turn))
1613    allocate(c%mom(4,4))
1614    allocate(c%A(4,4))
1615    allocate(c%At(4,4))
1616    c%turn=0
1617    c%kind=0
1618    c%r=0.0_dp
1619    c%bpm=0.0_dp
1620    c%xf=0.0_dp
1621    c%xn=0.0_dp
1622    c%mom=0.0_dp
1623    c%A=0.0_dp
1624    c%full=.true.
1625    do i=1,4
1626       c%a(i,i)=1.0_dp
1627    enddo
1628    c%At=0.0_dp
1629    do i=1,4
1630       c%at(i,i)=1.0_dp
1631    enddo
1632  end subroutine alloc_fibre_monitor_data
1633
1634  subroutine kill_fibre_monitor_data(c)
1635    implicit none
1636    TYPE(fibre_monitor_data), pointer :: c
1637    integer turn
1638
1639    deallocate(c%turn)
1640    deallocate(c%kind)
1641    deallocate(c%r)
1642    deallocate(c%bpm)
1643    deallocate(c%xf)
1644    deallocate(c%xn)
1645    deallocate(c%mom)
1646    deallocate(c%a)
1647    deallocate(c%at)
1648    nullify(c%turn)
1649    nullify(c%kind)
1650    nullify(c%r)
1651    nullify(c%xf)
1652    nullify(c%xn)
1653    nullify(c%mom)
1654    nullify(c%a)
1655    nullify(c%at)
1656
1657  end subroutine kill_fibre_monitor_data
1658
1659
1660!!!1 this is the correct one
1661  subroutine  alex_mom_real_monitors(ring,jm,jn,x,state_in,no,nt)
1662    implicit none
1663    integer jm,nm,i,jinv(4),i1,i11,i2,i22,jn(4),mf1,mf2,no,nt
1664    type(layout), target :: ring
1665    real(dp), pointer :: mom(:,:),r(:,:),a(:,:),at(:,:)
1666    real(dp)ar,x(6),z(4),zz(lnv),zx(4),z0(4),zf(4)
1667    type(damap) id,m11,m2,m22,m1,ex,map1,fin
1668    type(pbfield) h
1669    type(normalform) norm
1670    type(fibre),pointer::p1,p11,p2,p22
1671    type(real_8) y(6)
1672    type(internal_state) state,state_in
1673    type(gmap) g
1674!!! will produce the real data at jm
1675
1676    z0=0.0_dp
1677
1678    nm=size(monitors)
1679
1680    state=state_in+only_4d0
1681
1682    CALL INIT(STATE,no,0)
1683
1684    call alloc(y)
1685    call alloc(id,m11,m2,m22,m1,ex)
1686    call alloc(map1)
1687    call alloc(g)
1688
1689    CALL FIND_ORBIT(ring,x,1,STATE,1e-5_dp)
1690    write(6,*) monitors(jm)%kind
1691    if(monitors(jm)%kind==1.or.monitors(jm)%kind==3) then
1692       jn(1)=jm
1693
1694       p1=>monitors(jm)%p
1695       if(monitors(jm)%kind==3) then
1696          p2=>monitors(jm)%p
1697          i2=jm
1698          jn(3)=jm
1699       else
1700          do i=jm+1,jm+nm
1701             if(monitors(mod_n(i,nm))%kind==2) then
1702                p2=>monitors(mod_n(i,nm))%p
1703                jn(3)=mod_n(i,nm)
1704                i2=i
1705                exit
1706             endif
1707          enddo
1708
1709       endif
1710
1711       do i=jm+1,jm+nm
1712          if(monitors(mod_n(i,nm))%kind==1.or.monitors(mod_n(i,nm))%kind==3) then
1713             p11=>monitors(mod_n(i,nm))%p
1714             jn(2)=mod_n(i,nm)
1715             i11=i
1716             exit
1717          endif
1718       enddo
1719
1720       if(monitors(mod_n(jn(2),nm))%kind==3) then
1721          p22=>monitors(mod_n(jn(2),nm))%p
1722          jn(4)=jn(2)
1723          i22=i2+1
1724       else
1725          do i=i2+1,jm+nm
1726             if(monitors(mod_n(i,nm))%kind==2) then
1727                p22=>monitors(mod_n(i,nm))%p
1728                jn(4)=mod_n(i,nm)
1729                i22=i
1730                exit
1731             endif
1732          enddo
1733       endif
1734
1735       write(6,*)
1736       write(6,*) p1%mag%name,p11%mag%name
1737       write(6,*) p2%mag%name,p22%mag%name
1738       write(6,*) p1%pos,p11%pos
1739       write(6,*) p2%pos,p22%pos
1740
1741       ! track from ip to p1 : jm x monitor
1742
1743       call track(x,state,ring%start,p1)
1744       !
1745       ! x is closed orbit at p1
1746
1747
1748
1749       id=1
1750
1751       y=x+id
1752       call track(y,state,p1,p11)
1753       m11=y
1754       ! m11 is map from p1 to p11  (second x monitor)
1755       m11=z0  ! removed zeroth order
1756
1757       y=x+id
1758       call track(y,state,p1,p2)
1759       m2=y
1760       m2=z0
1761       ! m2 is map from p1 to p2  (first y monitor)
1762
1763
1764       y=x+id
1765       call track(y,state,p1,p22)
1766       m22=y
1767       m22=z0
1768       ! m22 is map from p1 to p22  (second y monitor)
1769
1770!!!!!!!!!!!!   testing and computing m11 !!!!!!!!!!!!!!!!!
1771       g=1
1772
1773
1774       jinv=0
1775       jinv(1)=1
1776       ex%v(1)=1.0_dp.mono.2
1777       ex%v(2)=1.0_dp.mono.1
1778       ex%v(3)=1.0_dp.mono.3
1779       ex%v(4)=1.0_dp.mono.4
1780       m11=m11*ex
1781       m11=m11**jinv
1782       ex%v(1)=1.0_dp.mono.2
1783       ex%v(2)=1.0_dp.mono.1
1784       ex%v(3)=1.0_dp.mono.3
1785       ex%v(4)=1.0_dp.mono.4
1786       m11=ex*m11*ex
1787
1788       ! m11 is now the map from x,x11,y,py
1789
1790       g%v(2)=m11%v(2)     !   px_0(x_0,x_p11,y_0,py_0)
1791
1792       ! testing tracking from jm to p2
1793
1794       jinv=0
1795       jinv(3)=1
1796       ex%v(1)=1.0_dp.mono.1
1797       ex%v(2)=1.0_dp.mono.2
1798       ex%v(3)=1.0_dp.mono.4
1799       ex%v(4)=1.0_dp.mono.3
1800       m22=m22*ex
1801       m22=m22**jinv
1802       ex%v(1)=1.0_dp.mono.1
1803       ex%v(2)=1.0_dp.mono.2
1804       ex%v(3)=1.0_dp.mono.4
1805       ex%v(4)=1.0_dp.mono.3
1806       m22=ex*m22*ex
1807       ! m22 is now the map from x,px,y,y22
1808
1809
1810       g%v(4)=m22%v(4)
1811
1812!!!!!!!!!!!!   testing and computing m2  !!!!!!!!!!!!!!!!!
1813
1814       jinv=0
1815       jinv(3)=1
1816       ex%v(1)=1.0_dp.mono.1
1817       ex%v(2)=1.0_dp.mono.2
1818       ex%v(3)=1.0_dp.mono.3
1819       ex%v(4)=1.0_dp.mono.4
1820       m2=m2*ex
1821       m2=m2**jinv
1822       !m11=m11**(-1)
1823       ex%v(1)=1.0_dp.mono.1
1824       ex%v(2)=1.0_dp.mono.2
1825       ex%v(3)=1.0_dp.mono.3
1826       ex%v(4)=1.0_dp.mono.4
1827       m2=ex*m2*ex
1828
1829       ! m2 is now the map from x,px,y2,py
1830
1831       g%v(3)=m2%v(3)
1832
1833
1834    else  !!!! if(monitors(jm)%kind==2)
1835       jn(3)=jm
1836
1837       !write(6,*) " crotte "
1838       !pause
1839       p2=>monitors(jm)%p
1840
1841       do i=jm+1,jm+nm
1842          if(monitors(mod_n(i,nm))%kind==2.or.monitors(mod_n(i,nm))%kind==3) then
1843             p22=>monitors(mod_n(i,nm))%p
1844             jn(4)=mod_n(i,nm)
1845             i22=i
1846             exit
1847          endif
1848       enddo
1849
1850       if(monitors(mod_n(jn(4),nm))%kind==3) then
1851
1852          p1=>monitors(mod_n(jn(4),nm))%p
1853          i1=i22
1854          jn(2)=jn(4)
1855       else
1856          do i=jm+1,jm+nm
1857             if(monitors(mod_n(i,nm))%kind==1) then
1858                p1=>monitors(mod_n(i,nm))%p
1859                jn(2)=mod_n(i,nm)
1860                i1=i
1861                exit
1862             endif
1863          enddo
1864
1865       endif
1866
1867
1868       do i=i22+1,jm+nm
1869          if(monitors(mod_n(i,nm))%kind==1.or.monitors(mod_n(i,nm))%kind==3) then
1870             p11=>monitors(mod_n(i,nm))%p
1871             jn(1)=mod_n(i,nm)
1872             i11=i
1873             exit
1874          endif
1875       enddo
1876       write(6,*)
1877       write(6,*) p1%mag%name,p11%mag%name
1878       write(6,*) p2%mag%name,p22%mag%name
1879       write(6,*) p1%pos,p11%pos
1880       write(6,*) p2%pos,p22%pos
1881
1882       ! track from ip to p1 : jm x monitor
1883
1884       call track(x,state,ring%start,p2)
1885       !
1886       ! x is closed orbit at p2
1887
1888
1889
1890       id=1
1891
1892       y=x+id
1893       call track(y,state,p2,p22)
1894       m22=y
1895       ! m22 is map from p2 to p22  (second y monitor)
1896       m22=z0  ! removed zeroth order
1897
1898       y=x+id
1899       call track(y,state,p2,p1)
1900       m1=y
1901       m1=z0
1902       ! m1 is map from p2 to p1  (first x monitor)
1903
1904
1905       y=x+id
1906       call track(y,state,p2,p11)
1907       m11=y
1908       m11=z0
1909       ! m11 is map from p2 to p11  (second x monitor)
1910
1911!!!!!!!!!!!!   testing and computing m22 !!!!!!!!!!!!!!!!!  asdad
1912       g=1
1913
1914
1915       jinv=0
1916       jinv(3)=1
1917       ex%v(1)=1.0_dp.mono.1
1918       ex%v(2)=1.0_dp.mono.2
1919       ex%v(3)=1.0_dp.mono.4
1920       ex%v(4)=1.0_dp.mono.3
1921       m22=m22*ex
1922       m22=m22**jinv
1923       ex%v(1)=1.0_dp.mono.1
1924       ex%v(2)=1.0_dp.mono.2
1925       ex%v(3)=1.0_dp.mono.4
1926       ex%v(4)=1.0_dp.mono.3
1927       m22=ex*m22*ex
1928
1929       g%v(4)=m22%v(4)     !   py_0(x_0,p_x0,y_0,y_p22)
1930
1931       ! m22 is now the map from x,px,y,y22
1932
1933!!!!!!!!!!!!   testing and computing m1  !!!!!!!!!!!!!!!!!
1934
1935       jinv=0
1936       jinv(1)=1
1937       ex%v(1)=1.0_dp.mono.2
1938       ex%v(2)=1.0_dp.mono.1
1939       ex%v(3)=1.0_dp.mono.3
1940       ex%v(4)=1.0_dp.mono.4
1941       m1=m1*ex
1942       m1=m1**jinv
1943       ex%v(1)=1.0_dp.mono.2
1944       ex%v(2)=1.0_dp.mono.1
1945       ex%v(3)=1.0_dp.mono.3
1946       ex%v(4)=1.0_dp.mono.4
1947       m1=ex*m1*ex
1948       ! m1 is now the map from x,x1,y,py
1949
1950       g%v(2)=m1%v(2)
1951
1952
1953!!!!!!!!!!!!   testing and computing m11  !!!!!!!!!!!!!!!!!
1954       jinv=0
1955       jinv(1)=1
1956       ex%v(1)=1.0_dp.mono.1
1957       ex%v(2)=1.0_dp.mono.2
1958       ex%v(3)=1.0_dp.mono.3
1959       ex%v(4)=1.0_dp.mono.4
1960       m11=m11*ex
1961       m11=m11**jinv
1962
1963       ex%v(1)=1.0_dp.mono.1
1964       ex%v(2)=1.0_dp.mono.2
1965       ex%v(3)=1.0_dp.mono.3
1966       ex%v(4)=1.0_dp.mono.4
1967       m11=ex*m11*ex
1968
1969       ! m11 is now the map from x2,px,y,py
1970
1971       g%v(1)=m11%v(1)
1972
1973
1974    endif !!!! if(monitors(jm)%kind==2)
1975
1976    !call kanalnummer(mf1,"original.dat")
1977    !call kanalnummer(mf2,"new.dat")
1978    zf=0.0_dp
1979    zf(1)=monitors(jn(1))%r(1,1)
1980    zf(2)=monitors(jn(2))%r(1,1)
1981    zf(3)=monitors(jn(3))%r(3,1)
1982    zf(4)=monitors(jn(4))%r(3,1)
1983    write(6,*) zf(1:4)
1984    write(6,*) monitors(jm)%r(1:4,1)
1985    !pause 200
1986    call invert_monitors(monitors(jm)%kind,zf,g)
1987    write(6,*) zf(1:4)
1988    !pause 201
1989
1990    write(6,*) " M_turn used in Eikonal",m_turn
1991    do i=1,nt
1992       zf=0.0_dp
1993       zf(1)=monitors(jn(1))%r(1,i)
1994       zf(2)=monitors(jn(2))%r(1,i)
1995       zf(3)=monitors(jn(3))%r(3,i)
1996       zf(4)=monitors(jn(4))%r(3,i)
1997       call invert_monitors(monitors(jm)%kind,zf,g)
1998       ! WRITE(mf1,'(4(1x,E15.8))')monitors(jm)%r(1:4,i)
1999       ! WRITE(mf2,'(4(1x,E15.8))')zf(1:4)
2000       monitors(jm)%xf(1:4,i)=zf(1:4)
2001    enddo
2002
2003    !close(mf1)
2004    !close(mf2)
2005
2006    call kill(y)
2007    call kill(id,m11,m2,m22,m1,ex)
2008    call kill(map1)
2009    call kill(g)
2010
2011
2012
2013  end subroutine  alex_mom_real_monitors
2014
2015  subroutine  alex_print_xf(jm,filename,n_max)
2016    implicit none
2017    integer jm,mf1,n,i
2018    CHARACTER(*) FILENAME
2019    integer, optional :: n_max
2020    call kanalnummer(mf1,filename)
2021
2022    n=m_turn
2023    if(present(n_max)) n=n_max
2024
2025    do i=1,n
2026       WRITE(mf1,'(4(1x,E15.8))')monitors(jm)%xf(1:4,i)
2027    enddo
2028
2029    close(mf1)
2030
2031  end   subroutine  alex_print_xf
2032
2033  subroutine  alex_read_r(sc,jm,filename,n_max,nav)
2034    implicit none
2035    integer jm,mf1,n,i,na
2036    CHARACTER(*), optional ::  FILENAME
2037    integer, optional :: n_max,nav
2038    real(dp) x,y,jj(4),xa(2),ya(2),sc
2039    call kanalnummer(mf1,filename)
2040
2041    n=m_turn
2042    na=0
2043    if(present(n_max)) n=n_max
2044    if(present(nav)) na=nav
2045
2046    if(present(filename)) then
2047       x=0.0_dp
2048       y=0.0_dp
2049       do i=1,n
2050          read(mf1,*)jj,monitors(jm)%r(1,i),monitors(jm)%r(3,i)
2051          monitors(jm)%r(3,i)=sc*monitors(jm)%r(3,i)
2052          x=monitors(jm)%r(1,i)+x
2053          y=monitors(jm)%r(3,i)+y
2054       enddo
2055
2056       x=x/n
2057       y=y/n
2058       write(6,*) "Averages ",x,y
2059       do i=1,n
2060          monitors(jm)%r(1,i)=(monitors(jm)%r(1,i)-x)/1000.0_dp
2061          monitors(jm)%r(3,i)=(monitors(jm)%r(3,i)-y)/1000.0_dp
2062       enddo
2063
2064    endif  ! filename
2065
2066
2067    if(na>0) then
2068       xa=0.0_dp
2069       ya=0.0_dp
2070
2071       do i=1,na
2072          xa(1)=monitors(jm)%r(1,i)**2+xa(1)
2073          ya(1)=monitors(jm)%r(3,i)**2+ya(1)
2074          xa(2)=monitors(jm)%r(1,n+1-i)**2+xa(2)
2075          ya(2)=monitors(jm)%r(3,n+1-i)**2+ya(2)
2076       enddo
2077       xa(2)=sqrt(xa(1)/xa(2))
2078       ya(2)=sqrt(ya(1)/ya(2))
2079
2080       write(6,*) " scales ",xa(2),ya(2)
2081       xa(1)=(-1.0_dp+xa(2) )/(n-1)
2082       ya(1)=(-1.0_dp+ya(2) )/(n-1)
2083       xa(2)=(n-xa(2))/(n-1)
2084       ya(2)=(n-ya(2))/(n-1)
2085
2086       do i=1,n
2087          monitors(jm)%r(1,i)=monitors(jm)%r(1,i) * (i*xa(1)+xa(2))
2088          monitors(jm)%r(3,i)=monitors(jm)%r(3,i) * (i*ya(1)+ya(2))
2089       enddo
2090
2091       x=0.0_dp
2092       y=0.0_dp
2093       do i=1,n
2094          x=monitors(jm)%r(1,i)+x
2095          y=monitors(jm)%r(3,i)+y
2096       enddo
2097
2098       x=x/n
2099       y=y/n
2100       write(6,*) "Averages ",x,y
2101       do i=1,n
2102          monitors(jm)%r(1,i)=(monitors(jm)%r(1,i)-x)
2103          monitors(jm)%r(3,i)=(monitors(jm)%r(3,i)-y)
2104       enddo
2105
2106    endif
2107
2108    close(mf1)
2109
2110    ! 171683 59 68 1 2.6794 -1.1248
2111
2112
2113  end   subroutine  alex_read_r
2114
2115  SUBROUTINE compute_twiss(r,my_state,filename,pos,del,no,thinlens,name,teng,notusingteng)
2116    IMPLICIT NONE
2117    TYPE(layout),target,INTENT(INOUT):: r
2118    TYPE(internal_state), intent(in):: my_state
2119    integer pos,no,I,mf,ind,j
2120    TYPE(internal_state) state
2121    real(dp) closed(6),del,bx,by,ax,ay,s,px,py,ex,exp,dpx,dpy,phi,E1(6,6),E2(6,6)
2122    type(DAMAP) ID,a_f,a_l,a_nl,DR,R_TE,cs_te
2123    logical(lp) COSLIKE,thinlens, doname, printres
2124    TYPE(REAL_8) Y(6)
2125    CHARACTER(*) FILENAME
2126    TYPE(FIBRE), POINTER :: P
2127    TYPE(integration_node), POINTER :: t
2128    TYPE(NORMALFORM) NORM
2129    character(nlp) name
2130    character(9) typec
2131    logical(lp) printmap,teng,notusingteng
2132    printmap=my_false
2133    if(.not.associated(r%t)) then
2134       call make_node_layout(r)
2135    endif
2136    printres=my_false
2137    doname=my_true
2138    call context(name)
2139    if(index(name,'###')/=0) then
2140       doname=my_false
2141       printres=my_true
2142    endif
2143    typec=' '
2144    call kanalnummer(mf,filename)
2145
2146    STATE=my_state+nocavity0
2147    call print(state,6)
2148   ! pause 777
2149    closed=0.0_dp
2150    closed(5)=del
2151    CALL FIND_ORBIT(R,CLOSED,pos,STATE,1e-5_dp)
2152    write(6,*) "closed orbit "
2153    write(6,*) CLOSED
2154
2155    CALL INIT(STATE,no,0)
2156    CALL ALLOC(Y);   CALL ALLOC(NORM)
2157    call alloc(id,a_f,a_l,a_nl,DR,R_TE,cs_te)
2158    id=1
2159    Y=CLOSED+id
2160
2161    CALL TRACK(R,Y,pos,STATE)
2162    id=y
2163    if(index(filename,'MAP')/=0.or.index(filename,'map')/=0) then
2164      printmap=my_true
2165    endif
2166
2167    NORM=ID
2168    write(mf,*) " fractional tunes ",norm%tune(1:2)
2169    p=>r%start
2170    s=0.0_dp
2171    do i=1,pos-1
2172       s=s+p%mag%p%ld
2173       p=>p%next
2174    enddo
2175    px=0.0_dp
2176    py=0.0_dp
2177    if(teng) then
2178     call factor(norm%a_t,a_f,a_l,a_nl,DR=dr,r_te=r_te,cs_te=cs_te,COSLIKE=COSLIKE)
2179    else
2180     if(notusingteng) then
2181      call factor_am_special(norm%a_t,a_f,a_l,a_nl,DR=dr)
2182       else
2183      call factor(norm%a_t,a_f,a_l,a_nl,DR=dr)
2184     endif
2185    endif
2186   
2187    Y=CLOSED+norm%a_t
2188    !    bx,by,ay,ay,s,px,py
2189    bx=(norm%a_t%v(1).sub.'1')**2+(norm%a_t%v(1).sub.'01')**2
2190    by=(norm%a_t%v(3).sub.'001')**2+(norm%a_t%v(3).sub.'0001')**2
2191    ax=-(norm%a_t%v(1).sub.'1')*(norm%a_t%v(2).sub.'1')-(norm%a_t%v(1).sub.'01')*(norm%a_t%v(2).sub.'01')
2192    ay=-(norm%a_t%v(3).sub.'001')*(norm%a_t%v(4).sub.'001')-(norm%a_t%v(3).sub.'0001')*(norm%a_t%v(4).sub.'0001')
2193    ex=a_f%v(1).sub.'00001'
2194    exp=a_f%v(2).sub.'00001'
2195    if(teng) then
2196     phi=(r_te%v(1).sub.'1')-1.0_dp   !acos
2197    else
2198     phi=(norm%a_t%v(1).sub.'001')**2+(norm%a_t%v(1).sub.'0001')**2
2199    endif
2200    if(abs(phi)<epsflo) phi=0.0_dp
2201    if(teng) then
2202     if(COSLIKE) then
2203       typec=' = COS-1 '
2204     else
2205       typec=' = COSH-1'
2206     endif
2207    endif
2208    if(doname) then
2209       if(index(p%mag%name,name(1:len_trim(name)))/=0) then
2210          printres=my_true
2211       endif
2212    endif
2213    if(teng) then
2214    write(mf,'(6x,a4,6x,(6x,a1,6x),2(4x,a5,4x),2(4x,a6,3x),2(4x,a6,3x),2(3x,a7,3x),a25)') &
2215         "name", "s", "betax","betay","alphax","alphay","eta_x ","etap_x","Phase x","Phase y","Teng-Edwards Cos or Cosh "
2216    else
2217    write(mf,'(6x,a4,6x,(6x,a1,6x),2(4x,a5,4x),2(4x,a6,3x),2(4x,a6,3x),2(3x,a7,3x),a25)') &
2218         "name", "s", "betax","betay","alphax","alphay","eta_x ","etap_x","Phase x","Phase y","  d<x**2>/de_y (Ripken)  "
2219    endif
2220    if(printres) write(mf,'(a16,10(1x,g12.5),1x,a9)') p%mag%name(1:16),s,bx,by,ax,ay,ex,exp,px,py,phi,typec
2221    if(doname) printres=my_false
2222
2223    do i=pos,pos+r%n-1
2224       t=>p%t1
2225       do while(.not.associated(t,p%t2%next))
2226          if(thinlens) then
2227             call track_probe_x(y,STATE, node1=t,node2=t%next)
2228          else
2229             call track_probe_x(y,STATE, node1=t,node2=p%t2%next)
2230          endif
2231          closed=y
2232          norm%a_t=y
2233        if(teng) then
2234         call factor(norm%a_t,a_f,a_l,a_nl,DR=dr,r_te=r_te,cs_te=cs_te,COSLIKE=COSLIKE)
2235        else
2236          if(notusingteng) then
2237           call factor_am_special(norm%a_t,a_f,a_l,a_nl,DR=dr)
2238            else
2239           call factor(norm%a_t,a_f,a_l,a_nl,DR=dr)
2240          endif       
2241        endif
2242
2243!          call factor(norm%a_t,a_f,a_l,a_nl,DR=dr,r_te=r_te,cs_te=cs_te,COSLIKE=COSLIKE)
2244          px=px+asin(dr%v(1).sub.'01')/twopi
2245          py=py+asin(dr%v(3).sub.'0001')/twopi
2246          Y=CLOSED+norm%a_t
2247
2248          if(thinlens) then
2249             t=>t%next
2250          else
2251             t=>p%t2%next
2252          endif
2253       enddo
2254       s=s+p%mag%p%ld
2255       p=>p%next
2256
2257       bx=(norm%a_t%v(1).sub.'1')**2+(norm%a_t%v(1).sub.'01')**2
2258       by=(norm%a_t%v(3).sub.'001')**2+(norm%a_t%v(3).sub.'0001')**2
2259       ax=-(norm%a_t%v(1).sub.'1')*(norm%a_t%v(2).sub.'1')-(norm%a_t%v(1).sub.'01')*(norm%a_t%v(2).sub.'01')
2260       ay=-(norm%a_t%v(3).sub.'001')*(norm%a_t%v(4).sub.'001')-(norm%a_t%v(3).sub.'0001')*(norm%a_t%v(4).sub.'0001')
2261       ex=a_f%v(1).sub.'00001'
2262       exp=a_f%v(2).sub.'00001'
2263       if(teng) then
2264        phi=(r_te%v(1).sub.'1')-1.0_dp   !acos
2265       else
2266        phi=(norm%a_t%v(1).sub.'001')**2+(norm%a_t%v(1).sub.'0001')**2
2267       endif
2268       if(abs(phi)<epsflo) phi=0.0_dp
2269        if(teng) then
2270         if(COSLIKE) then
2271           typec=' = COS-1 '
2272         else
2273           typec=' = COSH-1'
2274         endif
2275        endif
2276       if(doname) then
2277          if(index(p%mag%name,name(1:len_trim(name)))/=0) then
2278             printres=my_true
2279          endif
2280       endif
2281       if(printres) write(mf,'(a16,10(1x,g12.5),1x,a9)') p%mag%name(1:16),s,bx,by,ax,ay,ex,exp,px,py,phi,typec
2282       if(doname) printres=my_false
2283    enddo
2284    if(printmap) then
2285        write(6,*) "stable_da, check_stable",stable_da, check_stable
2286        stable_da=my_true
2287        check_stable=my_true
2288        call print(id,mf)
2289    endif
2290
2291    a_l=norm%a_t
2292    a_f=a_l**(-1)
2293    a_nl=0
2294    a_nl%v(1)=1.0_dp.mono.2
2295    a_nl%v(2)=-(1.0_dp.mono.1)
2296    a_nl%v(3)=1.0_dp.mono.4
2297    a_nl%v(4)=-(1.0_dp.mono.3)
2298   
2299    id=1
2300    id%v(3)=0
2301    id%v(4)=0   
2302    DR=A_l*a_nl*id*a_f*a_nl
2303    e1=dr
2304    e1=-e1
2305    id=1
2306    id%v(1)=0
2307    id%v(2)=0   
2308    DR=A_l*a_nl*id*a_f*a_nl
2309    e2=dr
2310    e2=-e2
2311   
2312  !! write(mf,'(6x,a4,6x,(6x,a1,6x),2(4x,a5,4x),2(4x,a6,3x),2(4x,a6,3x),2(3x,a7,3x),a25)')
2313          write(mf,*)  "   "
2314          write(mf,*)  "  All 4x4 Ripken Lattice Functions  "
2315          write(mf,*)  "   "
2316          write(mf,*)  " i , j, d<x_i x_j>/dJ_1  , d<x_i x_j>/dJ_2  "
2317    do i=1,4
2318    do j=i,4
2319        write(mf,'(1x,i2,2x,i2,2x,(1x,g12.5),(9x,g12.5))')  i,j,e1(i,j),e2(i,j)
2320    enddo
2321    enddo
2322          write(mf,*)  "   "
2323          write(mf,*)  "   4x4 Dispersions  "
2324          write(mf,*)  "   "
2325     do i=1,4
2326         e1(1,1)=a_l%v(i).sub.'00001'
2327        write(mf,'(1x,i2,(1x,g12.5))')  i,e1(1,1)
2328    enddo
2329   
2330    CLOSE(mf)
2331 
2332    CALL kill(Y)
2333    call kill(NORM)
2334    call kill(id,a_f,a_l,a_nl,DR,R_TE,cs_te)
2335
2336  end SUBROUTINE compute_twiss
2337 
2338!!!! special rcs
2339
2340  subroutine lattice_fit_bump_rcs(R,EPSF)
2341    IMPLICIT NONE
2342    TYPE(layout), target,intent(inout):: R
2343    integer, parameter :: neq=6,np=7
2344    real(dp)  TARG(neq)
2345    real(dp) CLOSED(6)
2346    TYPE(INTERNAL_STATE) STATE
2347    INTEGER I,SCRATCHFILE, MF
2348    TYPE(TAYLOR) EQ(neq)
2349    TYPE(REAL_8) Y(6)
2350    integer ::  no=2,nt,j,it,pos1,pos2
2351    type(damap) id
2352    type(gmap) g
2353    TYPE(TAYLOR)t
2354    real(dp) epsf,epsr,epsnow,gam(2)
2355    type(fibre), pointer:: p
2356    TYPE(POL_BLOCK) poly(np)
2357   
2358do i=1,np
2359 poly(i)=0
2360enddo
2361
2362poly(1)%name='QDX2701'
2363poly(2)%name='QFL0101'
2364poly(3)%name='QDL0101'
2365poly(4)%name='QFM0201'
2366poly(5)%name='QDL0201'
2367poly(6)%name='QFL0301'
2368poly(7)%name='QDX0301'
2369
2370do i=1,np
2371 poly(i)%ibn(2)=i
2372 r=poly(i)
2373enddo
2374!r%lastpos=1
2375!r%last=>R%start
2376call move_to( R,p,poly(1)%name,pos1,reset=my_true)
2377call move_to( R,p,poly(7)%name,pos2)
2378pos2=pos2+1
2379
2380SET_TPSAFIT=.FALSE.
2381
2382 
2383targ(1)=-1.216091585893243e0_dp     !  1 1
2384targ(2)=-6.061686871135153e0_dp     !  1 2
2385targ(3)=-0.7900090444471683E-01_dp  ! 2 1
2386targ(4)=2.120414509984705e0_dp      ! 3 3
2387targ(5)=-19.02200463334060e0_dp     ! 3 4
2388targ(6)=-0.1837954391003473e0_dp    ! 4 3
2389
2390
2391
2392
2393
2394    epsr=abs(epsf)
2395
2396    nt=neq+np
2397    STATE=only_4d0
2398
2399    CALL INIT(STATE,no,NP)
2400
2401
2402    it=0
2403100 continue
2404    it=it+1
2405
2406    CLOSED=0.0_dp
2407
2408
2409    CALL INIT(STATE,no,NP)
2410    CALL ALLOC(Y)
2411    CALL ALLOC(EQ)
2412    call alloc(id)
2413
2414    id=1
2415    Y=CLOSED+id
2416
2417    CALL TRACK(R,Y,pos1,pos2,+STATE)
2418
2419    !stop
2420
2421 !   write(6,*) "c_%no,c_%nv,c_%nd,c_%nd2"
2422 !   write(6,*) c_%no,c_%nv,c_%nd,c_%nd2
2423 !   write(6,*) "c_%ndpt,c_%npara,c_%npara,c_%np_pol"
2424 !   write(6,*)  c_%ndpt,c_%npara,c_%npara,c_%np_pol
2425
2426    eq(1)=(y(1)%t.par.'1000')-targ(1)
2427    eq(2)=(y(1)%t.par.'0100')-targ(2)
2428    eq(3)=(y(2)%t.par.'1000')-targ(3)
2429    eq(4)=(y(3)%t.par.'0010')-targ(4)
2430    eq(5)=(y(3)%t.par.'0001')-targ(5)
2431    eq(6)=(y(4)%t.par.'0010')-targ(6)
2432
2433    epsnow=0.0_dp
2434    do i=1,neq
2435     epsnow=epsnow+abs(eq(i))
2436    enddo
2437     write(6,*) " deviation ",epsnow
2438     
2439    call kanalnummer(SCRATCHFILE)
2440    OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
2441    rewind scratchfile
2442
2443    do i=1,neq
2444       eq(i)=eq(i)<=c_%npara
2445    enddo
2446    do i=1,neq
2447       call daprint(eq(i),scratchfile)
2448    enddo
2449    close(SCRATCHFILE)
2450
2451    CALL KILL(Y)
2452    CALL KILL(id)
2453    CALL KILL(EQ)
2454
2455
2456
2457    CALL INIT(1,nt)
2458    call alloc(g,nt)
2459    call kanalnummer(SCRATCHFILE)
2460    OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
2461    rewind scratchfile
2462    do i=np+1,nt
2463       call read(g%v(i),scratchfile)
2464    enddo
2465    close(SCRATCHFILE)
2466
2467    call alloc(t)
2468    do i=1,np
2469       g%v(i)=1.0_dp.mono.i
2470       do j=np+1,nt
2471          t=g%v(j).d.i
2472          g%v(i)=g%v(i)+(1.0_dp.mono.j)*t
2473       enddo
2474    enddo
2475    CALL KILL(t)
2476
2477    g=g.oo.(-1)
2478    tpsafit=0.0_dp
2479    tpsafit(1:nt)=g
2480
2481    SET_TPSAFIT=.true.
2482
2483    do i=1,7
2484     r=poly(i)
2485    enddo
2486
2487
2488    SET_TPSAFIT=.false.
2489
2490    CALL ELP_TO_EL(R)
2491
2492    !    write(6,*) " more "
2493    !    read(5,*) more
2494    if(it>=max_fit_iter) goto 101
2495    if(epsnow<=epsr) goto 102
2496    GOTO 100
2497
2498101 continue
2499    write(6,*) " warning did not converge "
2500
2501102 continue
2502    CALL KILL_PARA(R)
2503
2504
2505  end subroutine lattice_fit_bump_rcs
2506 
2507    subroutine lattice_fit_bump_min_rcs(R0,R1,EPSF,poly,np,sca)
2508    IMPLICIT NONE
2509    TYPE(layout), target,intent(inout):: R0,R1
2510    integer, parameter :: neq=2
2511    real(dp)  TARG(neq),bety0m,bety1m,bety1,bety0,dbx,dby,betx1ma,bety1ma
2512    real(dp) X0(6),X1(6),CLOSED(6),betx0m,betx1m,betx1,betx0,dbyr,dbyrm,dbxr,dbxrm
2513    TYPE(INTERNAL_STATE) STATE
2514    INTEGER I,SCRATCHFILE, MF
2515    TYPE(TAYLOR) EQ(neq)
2516    TYPE(REAL_8) Y(6),Y1(6),Y0(6)
2517    integer ::  no=2,nt,j,it,pos1,pos2,np
2518    type(damap) id1,ID0
2519    type(gmap) g
2520    TYPE(TAYLOR)t,bp(2),v
2521    type(taylor), allocatable :: dr(:)
2522    real(dp) epsf,epsr,epsnow,gam(2),sca
2523    type(fibre), pointer:: p,px,py
2524    TYPE(POL_BLOCK) poly(:)
2525    type(normalform) norm
2526    real(dp) s0,ds
2527 !   write(6,*)" correct ? "
2528 !   read(5,*) j
2529 !   if(j==1)    call lattice_fit_bump_rcs(R1,EPSF)
2530!do i=1,np
2531! poly(i)=0
2532!enddo
2533
2534!poly(1)%name='QDX2701'
2535!poly(2)%name='QFL0101'
2536!poly(3)%name='QDL0101'
2537!poly(4)%name='QFM0201'
2538!poly(5)%name='QDL0201'
2539!poly(6)%name='QFL0301'
2540!poly(7)%name='QDX0301'
2541
2542!do i=1,np
2543! poly(i)%ibn(2)=i
2544! r1=poly(i)
2545!enddo
2546
2547    do i=1,np
2548     r1=poly(i)
2549    enddo
2550
2551SET_TPSAFIT=.FALSE.
2552    STATE=only_4d0
2553
2554
2555    CALL FIND_ORBIT(R0,X0,1,STATE,1e-5_dp)
2556    write(6,*) "closed orbit ", CHECK_STABLE
2557    write(6,*) X0
2558    CALL FIND_ORBIT(R1,X1,1,STATE,1e-5_dp)
2559    write(6,*) "closed orbit ", CHECK_STABLE
2560    write(6,*) X1
2561
2562   CALL INIT(STATE,1,0)
2563   
2564   CALL ALLOC(Y0)
2565   CALL ALLOC(Y1)
2566   CALL ALLOC(norm)
2567   CALL ALLOC(ID1,ID0)
2568   id1=1
2569   id0=1
2570   
2571   y0=x0+id0
2572   y1=x1+id1
2573
2574    CALL TRACK(R0,Y0,1,STATE)
2575    CALL TRACK(R1,Y1,1,STATE)
2576   
2577    norm=y0
2578    id0=norm%a_t
2579    targ(1:2) = norm%tune(1:2)
2580    norm=y1
2581    id1=norm%a_t
2582    y0=x0+id0
2583    y1=x1+id1
2584 
2585    betx0m=0
2586    betx1m=0
2587    bety0m=0
2588    bety1m=0
2589    dbxrm=0
2590    dbyrm=0
2591   
2592    dbxr=0
2593    dbyr=0
2594   
2595    dbx=0
2596    dby=0
2597    s0=0
2598    ds=0
2599   
2600    p=>r0%start
2601   
2602    do i=1,r0%n
2603   
2604    ds=p%mag%p%ld
2605   
2606      call TRACK(R0,Y0,i,i+1,STATE)     
2607      betx0=(y0(1)%t.sub.'1')**2+(y0(1)%t.sub.'01')**2
2608      bety0=(y0(3)%t.sub.'001')**2+(y0(3)%t.sub.'0001')**2
2609        call TRACK(R1,Y1,i,i+1,STATE)         
2610      betx1=(y1(1)%t.sub.'1')**2+(y1(1)%t.sub.'01')**2
2611      bety1=(y1(3)%t.sub.'001')**2+(y1(3)%t.sub.'0001')**2
2612     if(betx0>betx0m) betx0m=betx0
2613     if(betx1>betx1m) betx1m=betx1   
2614     if(bety0>bety0m) bety0m=bety0
2615     if(bety1>bety1m) bety1m=bety1 
2616     dbxr= abs(betx1-betx0)  !/bety0
2617     dbyr= abs(bety1-bety0)  !/bety0
2618     if(dbxr>dbxrm) then
2619      dbxrm=dbxr
2620     ! write(6,*)"x", betx0,betx1
2621      px=>p
2622     endif
2623     if(dbyr>dbyrm) then
2624      dbyrm=dbyr
2625      !write(6,*)"y", bety0,bety1
2626      py=>p
2627     endif
2628     dbx=ds*(betx1-betx0)**2+dbx
2629     dby=ds*(bety1-bety0)**2+dby
2630     p=>p%next
2631     s0=s0+ds
2632    enddo
2633   
2634     dbx=dbx/s0
2635     dby=dby/s0
2636    write(6,*) " maximum "
2637    write(6,*) sqrt(betx0m),sqrt(betx1m)
2638    write(6,*) sqrt(bety0m),sqrt(bety1m)
2639 !   write(6,*) sqrt(dbx),sqrt(dby),s0
2640 !   write(6,*) px%mag%name,py%mag%name
2641 !   write(6,*) dbxr,dbyr
2642       
2643   
2644   CALL kill(Y0)
2645   CALL kill(Y1)
2646   CALL kill(norm)
2647   CALL kill(ID1,ID0)
2648
2649
2650
2651    epsr=abs(epsf)
2652
2653    nt=neq+np
2654
2655    CALL INIT(STATE,no,NP)
2656
2657
2658    it=0
2659100 continue
2660    it=it+1
2661
2662    x1=0.0_dp
2663
2664
2665    CALL INIT(STATE,no,NP)
2666   
2667    CALL ALLOC(Y0)
2668    CALL ALLOC(Y1)
2669    CALL ALLOC(EQ)
2670    call alloc(id0)
2671    call alloc(id1)
2672    call alloc(bp)
2673    call alloc(norm)
2674    call alloc(t,v)
2675
2676    CALL FIND_ORBIT(R1,X1,1,STATE,1e-5_dp)
2677
2678    id0=1
2679    Y0=x0+id0
2680    id1=1
2681    Y1=x1+id1
2682
2683    CALL TRACK(R0,Y0,1,STATE)
2684    CALL TRACK(R1,Y1,1,+STATE)
2685
2686    norm=y0
2687    id0=norm%a_t
2688    targ(1:2) = norm%tune(1:2)
2689    targ(1)=targ(1)+0.02d0
2690    targ(2)=targ(2)+0.01d0
2691    norm=y1
2692    id1=norm%a_t
2693     
2694    y0=x0+id0
2695    y1=x1+id1
2696
2697
2698    betx1ma=0
2699    bety1ma=0
2700
2701    s0=0
2702    p=>r0%start
2703   
2704    do i=1,r0%n
2705   
2706    ds=p%mag%p%ld
2707   
2708      call TRACK(R0,Y0,i,i+1,STATE)     
2709      betx0=(y0(1)%t.sub.'1')**2+(y0(1)%t.sub.'01')**2
2710      bety0=(y0(3)%t.sub.'001')**2+(y0(3)%t.sub.'0001')**2
2711
2712
2713      call TRACK(R1,Y1,i,i+1,+STATE)         
2714      bp(1)=ds*((y1(1)%t.par.'1000')**2+(y1(1)%t.par.'0100')**2-betx0)**2+bp(1)
2715      bp(2)=ds*((y1(3)%t.par.'0010')**2+(y1(3)%t.par.'0001')**2-bety0)**2+bp(2)
2716     
2717      betx1=(y1(1)%t.sub.'1000')**2+(y1(1)%t.sub.'0100')**2
2718      bety1=(y1(3)%t.sub.'0010')**2+(y1(3)%t.sub.'0001')**2
2719      if(betx1>betx1ma) betx1ma=betx1
2720      if(bety1>bety1ma) bety1ma=bety1
2721
2722
2723
2724     p=>p%next
2725     s0=s0+ds
2726    enddo
2727
2728    write(6,*) " maximum "
2729    write(6,*) sqrt(betx0m),sqrt(betx1m)
2730    write(6,*) sqrt(bety0m),sqrt(bety1m)
2731    write(6,*) " maximum now"
2732    write(6,*) sqrt(betx0m),sqrt(betx1ma)
2733    write(6,*) sqrt(bety0m),sqrt(bety1ma)
2734
2735    !stop
2736
2737 !   write(6,*) "c_%no,c_%nv,c_%nd,c_%nd2"
2738 !   write(6,*) c_%no,c_%nv,c_%nd,c_%nd2
2739 !   write(6,*) "c_%ndpt,c_%npara,c_%npara,c_%np_pol"
2740 !   write(6,*)  c_%ndpt,c_%npara,c_%npara,c_%np_pol
2741
2742    eq(1)=       ((NORM%dhdj%v(1)).par.'0000')-targ(1)
2743    eq(2)=       ((NORM%dhdj%v(2)).par.'0000')-targ(2)
2744    bp(1)=bp(1)/s0
2745    bp(2)=bp(2)/s0
2746    epsnow=0.0_dp
2747    do i=1,neq
2748     epsnow=epsnow+abs(eq(i))
2749    enddo
2750   
2751     write(6,*) " deviation ",epsnow
2752!write(6,*) " scale "
2753!read(5,*) sca
2754    do i=1,neq
2755     eq(i)=eq(i)-(1.0_dp-sca)*(eq(i).sub.'0')
2756    enddo
2757    epsnow=0.0_dp
2758    do i=1,neq
2759     epsnow=epsnow+abs(eq(i))
2760    enddo
2761      write(6,*) " deviation ",epsnow
2762   
2763    call kanalnummer(SCRATCHFILE)
2764    OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
2765    rewind scratchfile
2766   
2767     t=(bp(1)+bp(2))
2768         ds=t
2769     write(6,*) " Merit function  ",ds
2770
2771    ! t=0
2772    ! do i=1,np
2773    !  t=t+half*(one.mono.(i+c_%nd2))**2
2774    ! enddo
2775     
2776    do i=1,np
2777     v=t.d.(i+c_%nd2)
2778     v=v<=c_%npara
2779       call daprint(v,scratchfile) 
2780    enddo
2781 
2782    do i=1,neq
2783        eq(i)=eq(i)<=c_%npara
2784      call daprint(eq(i),scratchfile)
2785    enddo
2786    close(SCRATCHFILE)
2787
2788    CALL KILL(Y0)
2789    CALL KILL(Y1)
2790    call KILL(id0)
2791    CALL KILL(id1)
2792    CALL KILL(EQ)
2793    call KILL(bp)
2794    call KILL(norm)
2795    call KILL(t,v)
2796
2797
2798     
2799    CALL INIT(1,nt)
2800   
2801    call alloc(g,nt)
2802    call alloc(t,v)
2803    allocate(dr(nt))
2804    call alloc(dr,nt)
2805   
2806    call kanalnummer(SCRATCHFILE)
2807    OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
2808    rewind scratchfile
2809    do i=1,nt
2810       call read(dr(i),scratchfile)
2811    enddo
2812
2813    do i=1,np
2814      t=0.0_dp
2815      do j=1,neq
2816       v=(dr(np+j).d.i)
2817       t=(1.0_dp.mono.(np+j))*v+t
2818      enddo
2819       g%v(i)=dr(i)+t
2820    enddo
2821        do j=1,neq
2822         g%v(np+j)=dr(np+j)
2823       enddo 
2824    close(SCRATCHFILE)
2825
2826
2827    g=g.oo.(-1)
2828    tpsafit=0.0_dp
2829    tpsafit(1:nt)=g
2830
2831    SET_TPSAFIT=.true.
2832
2833    do i=1,np
2834     r1=poly(i)
2835    enddo
2836
2837
2838    SET_TPSAFIT=.false.
2839
2840    CALL ELP_TO_EL(R1)
2841    CALL KILL(t,v)
2842    CALL KILL(g)
2843    call KILL(dr,nt)
2844    deallocate(dr)
2845
2846     !   write(6,*) " more "
2847     !   read(5,*) i
2848     !   if(i==0) goto 102
2849    if(it>=max_fit_iter/sca**2) goto 101
2850    if(epsnow<=epsr) goto 102
2851    GOTO 100
2852
2853101 continue
2854    write(6,*) " warning did not converge "
2855
2856102 continue
2857    CALL KILL_PARA(R1)
2858
2859
2860
2861
2862  end subroutine lattice_fit_bump_min_rcs
2863 
2864
2865 
2866end module S_fitting_new
Note: See TracBrowser for help on using the repository browser.