source: PSPA/madxPSPA/libs/ptc/src/Se_status.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: 68.9 KB
Line 
1!The Polymorphic Tracking Code
2!Copyright (C) Etienne Forest and CERN
3module S_status
4  use s_frame
5  USE S_extend_poly
6!  use anbn
7  ! use my_own_1D_TPSA
8  !  USE S_pol_user1
9  !  USE S_pol_user2
10  Use S_pol_sagan
11  implicit none
12  public
13
14  integer,private,parameter::ST=30
15  !
16  integer, parameter :: KIND0 = ST
17  integer, parameter :: KIND1 = ST+1
18  integer, parameter :: KIND2 = ST+2
19  integer, parameter :: KIND3 = ST+3
20  integer, parameter :: KIND4 = ST+4
21  integer, parameter :: KIND5 = ST+5
22  integer, parameter :: KIND6 = ST+6
23  integer, parameter :: KIND7 = ST+7
24  integer, parameter :: KIND8 = ST+8
25  integer, parameter :: KIND9 = ST+9
26  integer, parameter :: KIND10 =ST+10
27  integer, parameter :: KIND11 =ST+11
28  integer, parameter :: KIND12 =ST+12
29  integer, parameter :: KIND13 =ST+13
30  integer, parameter :: KIND14 =ST+14
31  integer, parameter :: KIND15 =ST+15
32  integer, parameter :: KIND16 =ST+16
33  integer, parameter :: KIND17 =ST+17
34  integer, parameter :: KIND18 =ST+18
35  integer, parameter :: KIND19 =ST+19
36  integer, parameter :: KIND20 =ST+20     !  MADLIKE wedges on RBEND
37  integer, parameter :: KIND21 =ST+21     !  travelling wave cavity
38  integer, parameter :: KIND22 =ST+22     !
39  integer, parameter :: KIND23 =ST+23     !
40  !  integer, parameter :: KINDFITTED = KIND23+1
41  !  integer, parameter :: KINDUSER1 = KIND23+2
42  !  integer, parameter :: KINDUSER2 = KIND23+3
43  integer, parameter :: KINDwiggler = KIND23+2
44  !  integer, parameter :: KINDmu      = KIND23+3
45  integer, parameter :: KINDpa     = KIND23+3
46  integer, parameter :: drift_kick_drift = kind2
47  integer, parameter :: matrix_kick_matrix = kind7
48  integer, parameter :: kick_sixtrack_kick = kind6
49  character(6) ind_stoc(ndim2)
50  !  Making PTC leaner is false
51  !
52
53  integer, target :: MADKIND2=KIND2
54  integer, target :: MADKIND3N=KIND3
55  integer, target :: MADKIND3S=KIND3
56  integer, pointer :: MADTHICK
57  integer, pointer :: MADTHIN_NORMAL
58  integer, pointer :: MADTHIN_SKEW
59  LOGICAL(lp), target :: MADLENGTH=.false.
60  LOGICAL(lp), target :: MAD=.false.
61  LOGICAL(lp), target :: EXACT_MODEL = .false.
62  INTEGER, target:: NSTD,METD
63  ! TYPE(B_CYL) SECTOR_B
64  !TYPE(B_CYL),ALLOCATABLE ::  S_B(:)
65  !  INTEGER, TARGET :: NDPT_OTHER = 0
66  real(dp) CRAD,CFLUC
67  !  real(dp) YOSK(0:4), YOSD(4)    ! FIRST 6TH ORDER OF YOSHIDA
68  !  real(dp),PARAMETER::AAA=-0.25992104989487316476721060727823e0_dp  ! fourth order integrator
69  !  real(dp),PARAMETER::FD1=half/(one+AAA),FD2=AAA*FD1,FK1=one/(one+AAA),FK2=(AAA-one)*FK1
70  INTEGER , target :: CAVITY_TOTALPATH=1   !  default is fake
71
72  LOGICAL(lp) :: firsttime_coef=.true.
73
74  PRIVATE EQUALt,ADD,PARA_REMA,EQUALtilt
75  !PRIVATE DTILTR,DTILTP,DTILTS
76  PRIVATE DTILTR_EXTERNAL,DTILTP_EXTERNAL
77  PRIVATE CHECK_APERTURE_R,CHECK_APERTURE_P !,CHECK_APERTURE_S
78  LOGICAL(lp), target:: electron
79  real(dp), target :: muon=1.0_dp
80  LOGICAL(lp),PRIVATE,PARAMETER::T=.TRUE.,F=.FALSE.
81  ! include "a_def_all_kind.inc"    ! sept 2007
82  ! include "a_def_sagan.inc"
83  ! include "a_def_element_fibre_layout.inc"
84  !  !  include "a_def_user1.inc"
85  !  !!  include "a_def_arbitrary.inc"
86  !  !  include "a_def_user2.inc"
87
88  TYPE(INTERNAL_STATE),PARAMETER::DEFAULT0=INTERNAL_STATE   (0,f,f,f,f,f,f,f,f,f,f,f)
89  TYPE(INTERNAL_STATE),PARAMETER::TOTALPATH0=INTERNAL_STATE (1,f,f,f,f,f,f,f,f,f,f,f)
90  TYPE(INTERNAL_STATE),PARAMETER::TIME0=INTERNAL_STATE      (0,t,f,f,f,f,f,f,f,f,f,f)
91  TYPE(INTERNAL_STATE),PARAMETER::RADIATION0=INTERNAL_STATE (0,f,t,f,f,f,f,f,f,f,f,f)
92  TYPE(INTERNAL_STATE),PARAMETER::NOCAVITY0=INTERNAL_STATE  (0,f,f,t,f,f,f,f,f,f,f,f)
93  TYPE(INTERNAL_STATE),PARAMETER::FRINGE0=INTERNAL_STATE    (0,f,f,f,t,f,f,f,f,f,f,f)
94  TYPE(INTERNAL_STATE),PARAMETER::STOCHASTIC0=INTERNAL_STATE(0,f,f,f,f,t,f,f,f,f,f,f)
95  TYPE(INTERNAL_STATE),PARAMETER::ENVELOPE0=INTERNAL_STATE  (0,f,f,f,f,f,t,f,f,f,f,f)
96  TYPE(INTERNAL_STATE),PARAMETER::ONLY_4d0=INTERNAL_STATE   (0,f,f,t,f,f,f,f,t,f,f,f)
97  TYPE(INTERNAL_STATE),PARAMETER::DELTA0=INTERNAL_STATE     (0,f,f,t,f,f,f,f,t,t,f,f)
98  TYPE(INTERNAL_STATE),PARAMETER::SPIN0=INTERNAL_STATE      (0,f,f,f,f,f,f,f,f,f,t,f)
99  TYPE(INTERNAL_STATE),PARAMETER::MODULATION0=INTERNAL_STATE(0,f,f,f,f,f,f,f,f,f,f,t)
100
101  TYPE(INTERNAL_STATE), target ::  DEFAULT=DEFAULT0
102  TYPE(INTERNAL_STATE), target ::  TOTALPATH=TOTALPATH0
103  TYPE(INTERNAL_STATE), target ::  RADIATION=RADIATION0
104  TYPE(INTERNAL_STATE), target ::  NOCAVITY=NOCAVITY0
105  TYPE(INTERNAL_STATE), target ::  FRINGE=FRINGE0
106  TYPE(INTERNAL_STATE), target ::  TIME=TIME0
107  TYPE(INTERNAL_STATE), target ::  STOCHASTIC=STOCHASTIC0
108  TYPE(INTERNAL_STATE), target ::  ENVELOPE=ENVELOPE0
109  TYPE(INTERNAL_STATE), target ::  ONLY_4D=ONLY_4D0
110  TYPE(INTERNAL_STATE), target ::  DELTA=DELTA0
111  TYPE(INTERNAL_STATE), target ::  SPIN=SPIN0
112  TYPE(INTERNAL_STATE), target ::  MODULATION=MODULATION0
113   type(acceleration), pointer :: acc
114   type(acceleration), pointer :: accFIRST
115   type(fibre), pointer :: paccfirst
116   type(fibre), pointer :: paccthen
117
118  !  private s_init,S_init_berz,MAKE_STATES_0,MAKE_STATES_m,print_s,CONV
119  private s_init,MAKE_STATES_0,MAKE_STATES_m,print_s,CONV
120  LOGICAL(lp), target :: compute_stoch_kick = .true.
121  private alloc_p,equal_p,dealloc_p,alloc_A,equal_A,dealloc_A
122  PRIVATE KILL_S_APERTURE,ALLOC_S_APERTURE
123  !,NULL_p
124  PRIVATE B2PERPR,B2PERPP !,S_init_berz0
125  type(tilting) tilt
126  private minu
127  real(dp) MADFAC(NMAX)
128  CHARACTER(24) MYTYPE(-100:100)
129  private check_S_APERTURE_r,check_S_APERTURE_out_r
130  private check_S_APERTURE_p,check_S_APERTURE_out_p
131  type(my_1D_taylor) val_del
132  logical(lp) :: ramp=my_false
133  logical(lp) :: accelerate=my_false, first_particle=my_false
134
135  TYPE B_CYL
136     integer firsttime
137     integer, POINTER ::  nmul,n_mono
138     integer, DIMENSION(:), POINTER   :: i,j
139     real(dp), DIMENSION(:,:), POINTER   :: a_x,a_y,b_x,b_y
140  END  TYPE B_CYL
141
142  TYPE(B_CYL),ALLOCATABLE ::  S_B(:)
143
144  INTERFACE OPERATOR (.min.)
145     MODULE PROCEDURE minu                       ! to define the minus of Schmidt
146  END INTERFACE
147
148  INTERFACE assignment (=)
149     MODULE PROCEDURE EQUALt
150     MODULE PROCEDURE EQUALtilt
151     MODULE PROCEDURE equal_p
152     MODULE PROCEDURE equal_A
153  end  INTERFACE
154
155  INTERFACE OPERATOR (+)
156     MODULE PROCEDURE add
157     MODULE PROCEDURE PARA_REMA
158  END INTERFACE
159
160  INTERFACE OPERATOR (-)
161     MODULE PROCEDURE sub
162  END INTERFACE
163
164  INTERFACE MAKE_STATES
165     MODULE PROCEDURE MAKE_STATES_m
166     MODULE PROCEDURE MAKE_STATES_0
167  END INTERFACE
168
169  INTERFACE CHECK_APERTURE
170     MODULE PROCEDURE CHECK_APERTURE_R
171     MODULE PROCEDURE CHECK_APERTURE_P
172     !     MODULE PROCEDURE CHECK_APERTURE_S
173  END INTERFACE
174
175  INTERFACE check_S_APERTURE
176     MODULE PROCEDURE check_S_APERTURE_r
177     MODULE PROCEDURE check_S_APERTURE_p
178  END INTERFACE
179
180  INTERFACE check_S_APERTURE_out
181     MODULE PROCEDURE check_S_APERTURE_out_r
182     MODULE PROCEDURE check_S_APERTURE_out_p
183  END INTERFACE
184
185  INTERFACE init
186     MODULE PROCEDURE s_init
187     !     MODULE PROCEDURE S_init_berz0
188     !     MODULE PROCEDURE S_init_berz
189  END INTERFACE
190
191  INTERFACE print
192     MODULE PROCEDURE print_s
193  END INTERFACE
194
195  INTERFACE alloc
196     MODULE PROCEDURE alloc_p
197     MODULE PROCEDURE alloc_A
198     MODULE PROCEDURE ALLOC_S_APERTURE
199  END INTERFACE
200
201  INTERFACE kill
202     MODULE PROCEDURE dealloc_p
203     MODULE PROCEDURE dealloc_A
204     MODULE PROCEDURE KILL_S_APERTURE
205  END INTERFACE
206
207  INTERFACE B2PERP
208     MODULE PROCEDURE B2PERPR
209     MODULE PROCEDURE B2PERPP
210  END INTERFACE
211
212  INTERFACE DTILTD
213     MODULE PROCEDURE DTILTR_EXTERNAL
214     MODULE PROCEDURE DTILTP_EXTERNAL       ! EXTERNAL
215  END INTERFACE
216
217
218
219CONTAINS
220
221  real(dp) function cradf(p)
222    implicit none
223    type (MAGNET_CHART), pointer:: P
224    cradf=crad*p%p0c**3
225  end function cradf
226
227  real(dp) function cflucf(p)
228    implicit none
229    type (MAGNET_CHART), pointer:: P
230    cflucf=cfluc*p%p0c**5
231  end function cflucf
232
233
234  SUBROUTINE  NULL_A(p)
235    implicit none
236    type (MADX_APERTURE), pointer:: P
237
238    nullify(P%KIND);nullify(P%R);nullify(P%X);nullify(P%Y);nullify(P%dX);nullify(P%dY);
239  end subroutine NULL_A
240
241  SUBROUTINE  alloc_A(p)
242    implicit none
243    type (MADX_APERTURE), pointer:: P
244
245    nullify(p)
246    allocate(p)
247    CALL NULL_A(p)
248    ALLOCATE(P%R(2));ALLOCATE(P%X);ALLOCATE(P%Y);ALLOCATE(P%KIND);
249    P%KIND=0; P%R=0.0_dp;P%X=0.0_dp;P%Y=0.0_dp;
250    ALLOCATE(P%DX);ALLOCATE(P%DY);
251    P%DX=0.0_dp;P%DY=0.0_dp;
252  end subroutine alloc_A
253
254  SUBROUTINE  dealloc_A(p)
255    implicit none
256    type (MADX_APERTURE), pointer:: P
257
258    if(associated(p%R)) then
259       DEALLOCATE(P%R);DEALLOCATE(P%X);DEALLOCATE(P%Y);DEALLOCATE(P%KIND);
260       DEALLOCATE(P%DX);DEALLOCATE(P%DY);
261    endif
262  end SUBROUTINE  dealloc_A
263
264
265
266  SUBROUTINE  NULL_p(p)
267    implicit none
268    type (MAGNET_CHART), pointer:: P
269
270    nullify(P%LD);nullify(P%B0);nullify(P%LC);
271    nullify(P%TILTD);  nullify(P%dir);
272    !    nullify(P%beta0);nullify(P%gamma0I);nullify(P%gambet);nullify(P%charge);
273    nullify(P%P0C);
274    nullify(P%EDGE)
275    !    nullify(P%TOTALPATH)
276    nullify(P%EXACT);  !nullify(P%RADIATION);nullify(P%NOCAVITY);
277    nullify(P%permFRINGE);
278    nullify(P%KILL_ENT_FRINGE);nullify(P%KILL_EXI_FRINGE);nullify(P%bend_fringe);  !nullify(P%TIME);
279    nullify(P%METHOD);nullify(P%NST);
280    nullify(P%NMUL);  !nullify(P%spin);
281    nullify(P%F);
282    nullify(P%APERTURE);
283    nullify(P%A);
284
285  end subroutine NULL_p
286
287
288
289
290  SUBROUTINE  alloc_p(p)
291    implicit none
292    type (MAGNET_CHART), pointer:: P
293
294    nullify(P);
295    ALLOCATE(P);
296    CALL NULL_P(P)
297
298    ALLOCATE(P%LD);ALLOCATE(P%B0);ALLOCATE(P%LC);
299    P%LD=0.0_dp;P%B0=0.0_dp;P%LC=0.0_dp;
300    ALLOCATE(P%TILTD);P%TILTD=0.0_dp;
301    ! ALLOCATE(P%beta0);
302    ! ALLOCATE(P%MASS0);
303    ! ALLOCATE(P%gamma0I);
304    ! ALLOCATE(P%gambet);
305    ALLOCATE(P%P0C);
306    ! P%beta0 =one;
307    ! P%MASS0 =one;
308    !P%gamma0I=zero;P%gambet =zero;
309    P%P0C =0.0_dp;
310    ALLOCATE(P%EDGE(2));P%EDGE(1)=0.0_dp;P%EDGE(2)=0.0_dp;
311    !    ALLOCATE(P%TOTALPATH); ! PART OF A STATE INITIALIZED BY EL=DEFAULT
312    ALLOCATE(P%EXACT);  !ALLOCATE(P%RADIATION);ALLOCATE(P%NOCAVITY);
313        ALLOCATE(P%permFRINGE);
314    ALLOCATE(P%KILL_ENT_FRINGE);ALLOCATE(P%KILL_EXI_FRINGE);ALLOCATE(P%bend_fringe); !ALLOCATE(P%TIME);
315    ALLOCATE(P%METHOD);ALLOCATE(P%NST);P%METHOD=2;P%NST=1;
316    ALLOCATE(P%NMUL);P%NMUL=0;
317    !    ALLOCATE(P%spin);
318    !   ALLOCATE(P%TRACK);P%TRACK=.TRUE.;
319    P%permFRINGE=.FALSE.
320    P%KILL_ENT_FRINGE=.FALSE.
321    P%KILL_EXI_FRINGE=.FALSE.
322    P%bend_fringe=.false.
323    call alloc(p%f)
324    ! if(junk) ccc=ccc+1
325
326  end subroutine alloc_p
327
328
329  SUBROUTINE  dealloc_p(p)
330    implicit none
331    type (MAGNET_CHART), pointer:: P
332    if(.not.associated(p)) return
333
334    !    if(associated(P%dir)) then
335    !    endif
336    if(associated(p%LD)) DEALLOCATE(P%LD);
337    if(associated(p%B0)) DEALLOCATE(P%B0);
338    if(associated(p%LC)) DEALLOCATE(P%LC);
339    if(associated(p%TILTD)) DEALLOCATE(P%TILTD);
340    !    if(associated(p%beta0)) DEALLOCATE(P%beta0);
341    !    if(associated(p%gamma0I)) DEALLOCATE(P%gamma0I);
342    !    if(associated(p%gambet)) DEALLOCATE(P%gambet);
343    if(associated(p%P0C)) DEALLOCATE(P%P0C);
344    if(associated(p%f)) then
345       call kill(p%f)
346       DEALLOCATE(P%f);
347    endif
348    if(associated(p%APERTURE)) then
349       CALL kill(p%APERTURE)
350       DEALLOCATE(p%APERTURE);
351    endif
352    if(associated(p%A)) then
353       CALL KILL(P%A)
354    ENDIF
355
356    if(associated(p%EDGE))DEALLOCATE(P%EDGE);
357    !    DEALLOCATE(P%TOTALPATH);
358    if(associated(p%EXACT))DEALLOCATE(P%EXACT);
359    !DEALLOCATE(P%RADIATION);DEALLOCATE(P%NOCAVITY);
360
361    if(associated(p%PERMFRINGE))DEALLOCATE(P%PERMFRINGE);
362    if(associated(p%KILL_ENT_FRINGE))DEALLOCATE(P%KILL_ENT_FRINGE);
363    if(associated(p%KILL_EXI_FRINGE))DEALLOCATE(P%KILL_EXI_FRINGE);
364    if(associated(p%bend_fringe))DEALLOCATE(P%bend_fringe); !DEALLOCATE(P%TIME);
365    if(associated(p%METHOD))DEALLOCATE(P%METHOD);
366    !DEALLOCATE(P%spin);
367    if(associated(p%NST))DEALLOCATE(P%NST);
368    if(associated(p%NMUL))DEALLOCATE(P%NMUL)
369    !    CALL NULL_P(P)
370    DEALLOCATE(P)
371    nullify(p);
372    ! if(junk) ccc=ccc-1
373  end subroutine dealloc_p
374
375
376  SUBROUTINE  KILL_S_APERTURE(A)
377    implicit none
378    TYPE(S_APERTURE), POINTER :: A(:)
379    INTEGER I
380    DO I=1,SIZE(A)
381       CALL kill(A(I)%APERTURE)
382       DEALLOCATE(A(I)%APERTURE)
383    ENDDO
384    DEALLOCATE(A)
385  END SUBROUTINE  KILL_S_APERTURE
386
387  SUBROUTINE  ALLOC_S_APERTURE(A,N,APERTURE)
388    implicit none
389    TYPE(S_APERTURE), POINTER :: A(:)
390    TYPE(MADX_APERTURE), OPTIONAL :: APERTURE
391    INTEGER I,N
392    ALLOCATE(A(N))
393    DO I=1,SIZE(A)
394       CALL ALLOC(A(I)%APERTURE)
395       IF(PRESENT(APERTURE))A(I)%APERTURE=APERTURE
396    ENDDO
397  END SUBROUTINE  ALLOC_S_APERTURE
398
399  SUBROUTINE  check_S_APERTURE_r(p,N,x)
400    implicit none
401    TYPE(magnet_chart), POINTER :: p
402    real(dp),intent(inout):: x(6)
403    INTEGER N
404    if(p%dir==1) then
405       call CHECK_APERTURE(p%a(n)%aperture,X)
406    else
407       call CHECK_APERTURE(p%a(p%nst+2-n)%aperture,X)
408    endif
409
410    !!   if(s_aperture_CHECK.and.associated(el%p%A)) call check_S_APERTURE(el%p,t%POS_IN_FIBRE-2)
411    !       if(associated(T%PARENT_FIBRE%MAG%p%aperture)) call CHECK_APERTURE(T%PARENT_FIBRE%MAG%p%aperture,X)
412  END SUBROUTINE  check_S_APERTURE_r
413
414  SUBROUTINE  check_S_APERTURE_p(p,N,x)
415    implicit none
416    TYPE(magnet_chart), POINTER :: p
417    type(real_8),intent(inout):: x(6)
418    INTEGER N
419
420    if(p%dir==1) then
421       call CHECK_APERTURE(p%a(n)%aperture,X)
422    else
423       call CHECK_APERTURE(p%a(p%nst+2-n)%aperture,X)
424    endif
425
426    !!   if(s_aperture_CHECK.and.associated(el%p%A)) call check_S_APERTURE(el%p,t%POS_IN_FIBRE-2)
427    !       if(associated(T%PARENT_FIBRE%MAG%p%aperture)) call CHECK_APERTURE(T%PARENT_FIBRE%MAG%p%aperture,X)
428  END SUBROUTINE  check_S_APERTURE_p
429
430  SUBROUTINE  check_S_APERTURE_out_r(p,N,x)
431    implicit none
432    TYPE(magnet_chart), POINTER :: p
433    real(dp),intent(inout):: x(6)
434    INTEGER N
435    if(p%dir==1.and.n==p%nst) then
436       call CHECK_APERTURE(p%a(n+1)%aperture,X)
437    elseif(p%dir==-1.and.n==p%nst) then
438       call CHECK_APERTURE(p%a(1)%aperture,X)
439    endif
440  END SUBROUTINE  check_S_APERTURE_out_r
441
442  SUBROUTINE  check_S_APERTURE_out_p(p,N,x)
443    implicit none
444    TYPE(magnet_chart), POINTER :: p
445    type(real_8),intent(inout):: x(6)
446    INTEGER N
447    if(p%dir==1.and.n==p%nst) then
448       call CHECK_APERTURE(p%a(n+1)%aperture,X)
449    elseif(p%dir==-1.and.n==p%nst) then
450       call CHECK_APERTURE(p%a(1)%aperture,X)
451    endif
452  END SUBROUTINE  check_S_APERTURE_out_p
453
454  SUBROUTINE  equal_A(elp,el)
455    implicit none
456    type (MADX_APERTURE),INTENT(inOUT)::elP
457    type (MADX_APERTURE),INTENT(IN)::el
458    ELP%KIND=EL%KIND
459    ELP%R=EL%R
460    ELP%X=EL%X
461    ELP%Y=EL%Y
462    ELP%DX=EL%DX
463    ELP%DY=EL%DY
464  END SUBROUTINE  equal_A
465
466
467
468  SUBROUTINE  equal_p(elp,el)
469    implicit none
470    type (MAGNET_CHART),INTENT(inOUT)::elP
471    type (MAGNET_CHART),INTENT(IN)::el
472
473    !    elp%beta0 =el%beta0
474    !    elp%gamma0I=el%gamma0I
475    !    elp%gambet =el%gambet
476    elp%P0C =el%P0C
477    elp%EXACT=el%EXACT
478    !   elp%RADIATION=el%RADIATION
479    !   elp%TIME=el%TIME
480    !   elp%NOCAVITY=el%NOCAVITY
481    !   elp%spin=el%spin
482       elp%permFRINGE=el%permFRINGE
483    elp%KILL_ENT_FRINGE=el%KILL_ENT_FRINGE
484    elp%KILL_EXI_FRINGE=el%KILL_EXI_FRINGE
485    elp%bend_fringe=el%bend_fringe
486
487    elp%LD=el%LD
488    elp%LC=el%LC
489    elp%TILTD=el%TILTD
490    elp%B0=el%B0
491    elp%EDGE(1)=el%EDGE(1)
492    elp%EDGE(2)=el%EDGE(2)
493    elp%METHOD=el%METHOD
494    elp%NST=el%NST
495    !   ELP%TOTALPATH=EL%TOTALPATH
496    ELP%nmul=EL%nmul
497    !    ELP%TRACK=EL%TRACK
498
499    if(associated(el%f)) then
500       elp%f=el%f
501    endif
502    if(associated(el%APERTURE)) then
503       if(.not.associated(elp%APERTURE)) call alloc(elp%APERTURE)
504       elp%APERTURE=el%APERTURE
505    endif
506
507  end SUBROUTINE  equal_p
508
509  SUBROUTINE  CHECK_APERTURE_R(E,X)
510    implicit none
511    type (MADX_APERTURE),INTENT(IN)::E
512    REAL(DP), INTENT(IN):: X(6)
513    !    real(dp) xx,yy,dx,dy,ex,ey
514
515    !  real(dp) :: xlost(6)=zero
516    !  character(120) :: messagelost
517
518    IF(CHECK_MADX_APERTURE.AND.APERTURE_FLAG) THEN
519       SELECT CASE(E%KIND)
520       CASE(1)  ! ellipse circles
521          IF((X(1)-E%DX)**2/E%R(1)**2+(X(3)-E%DY)**2/E%R(2)**2>1.0_dp) THEN
522             CHECK_STABLE=.FALSE.
523             STABLE_DA=.false.
524             xlost=0.0_dp
525             xlost=x
526             messagelost="Lost in real kind=1 elliptic Aperture"
527          ENDIF
528       CASE(2)  ! rectangle
529          IF(ABS(X(1)-E%DX)>E%X.OR.ABS(X(3)-E%DY)>E%Y) THEN
530             CHECK_STABLE=.FALSE.
531             STABLE_DA=.false.
532             xlost=0.0_dp
533             xlost=x
534             messagelost="Lost in real kind=2 rectangular Aperture"
535          ENDIF
536       CASE(3)  ! RECTANGLE + ELLIPSE (CIRCLE)
537          IF((ABS(X(1)-E%DX)>E%X).OR.(ABS(X(3)-E%DY)>E%Y).OR.  &
538               ((X(1)-E%DX)**2/E%R(1)**2+(X(3)-E%DY)**2**2/E%R(2)**2>1.0_dp)) THEN
539             CHECK_STABLE=.FALSE.
540             STABLE_DA=.false.
541             xlost=0.0_dp
542             xlost=x
543             messagelost="Lost in real kind=3 rect-ellipse Aperture"
544          ENDIF
545       CASE(4) ! MARGUERITE
546          IF(((X(1)-E%DX)**2/E%R(2)**2+(X(3)-E%DY)**2/E%R(1)**2>1.0_dp).OR.  &
547               ((X(1)-E%DX)**2/E%R(1)**2+(X(3)-E%DY)**2/E%R(2)**2>1.0_dp)) THEN
548             CHECK_STABLE=.FALSE.
549             STABLE_DA=.false.
550             xlost=0.0_dp
551             xlost=x
552             messagelost="Lost in real kind=4 marguerite Aperture"
553          ENDIF
554       CASE(5) ! RACETRACK
555          IF( (abs(x(1)-e%dx)) > (e%r(1)+e%x)                  &
556               .or. abs(x(3)-e%dy) .gt. (e%y+e%r(1)) .or.                &
557               ((((abs(x(1)-e%dx)-e%x)**2+                            &
558               (abs(x(3)-e%dy)-e%y)**2) .gt. e%r(1)**2)                  &
559               .and. (abs(x(1)-e%dx)) .gt. e%x                        &
560               .and. abs(x(3)-e%dy) .gt. e%y)) THEN
561             CHECK_STABLE=.FALSE.
562             STABLE_DA=.false.
563             xlost=0.0_dp
564             xlost=x
565             messagelost="Lost in real kind=5 racetrack Aperture"
566          ENDIF
567
568       CASE(6) ! PILES OF POINTS
569          STOP 222
570       CASE DEFAULT
571          !   STOP 223
572       END SELECT
573
574    ENDIF
575
576
577  END SUBROUTINE  CHECK_APERTURE_R
578
579  SUBROUTINE  CHECK_APERTURE_P(E,X)
580    implicit none
581    type (MADX_APERTURE),INTENT(IN)::E
582    TYPE(REAL_8), INTENT(IN):: X(6)
583    REAL(DP) Y(6)
584
585    Y=X
586    CALL CHECK_APERTURE(E,Y)
587
588  END SUBROUTINE  CHECK_APERTURE_P
589
590  FUNCTION minu( S1,S2  )
591    implicit none
592    logical(lp) minu
593    logical(lp), INTENT (IN) :: S1
594    logical(lp), INTENT (IN) :: S2
595
596    minu=.false.
597    if(s1.and.(.not.s2)) minu=my_true
598
599  END FUNCTION minu
600
601
602
603  SUBROUTINE  EQUALTILT(S2,S1)
604    implicit none
605    INTEGER I
606    type (TILTING),INTENT(OUT)::S2
607    type (TILTING),INTENT(IN)::S1
608
609    S2%NATURAL=S1%NATURAL
610    DO I=0,NMAX
611       S2%TILT(I)=S1%TILT(I)
612    ENDDO
613
614  END SUBROUTINE EQUALTILT
615
616
617  SUBROUTINE MAKE_STATES_0(particle)
618    USE   definition
619    USE   da_arrays
620    IMPLICIT NONE
621    LOGICAL(lp) particle
622    integer i,lda_old
623
624    W_P=>W_I
625    NULLIFY(ACC);       
626    NULLIFY(ACCfirst);       
627    NULLIFY(paccfirst);       
628    NULLIFY(paccthen);       
629
630
631    insane_PTC=.true.
632    NEWFILE%MF=.FALSE.
633    CLOSEFILE%MF=.FALSE.
634    KNOB=.FALSE.
635    SETKNOB=.TRUE.
636
637    CALL MAKE_YOSHIDA
638    MYTYPE(0)=" MARKER"
639    MYTYPE(kind1)=" DRIFT"
640    MYTYPE(kind2)=" DRIFT-KICK-DRIFT"
641    MYTYPE(kind3)=" THIN ELEMENT"
642    MYTYPE(kind4)=" RF CAVITY"
643    MYTYPE(kind5)=" SOLENOID "
644    MYTYPE(kind6)=" KICK-SixTrack-KICK "
645    MYTYPE(kind7)=" MATRIX-KICK-MATRIX "
646    MYTYPE(kind8)=" NORMAL SMI "
647    MYTYPE(kind9)=" SKEW   SMI "
648    MYTYPE(kind10)=" EXACT SECTOR "
649    MYTYPE(kind11)=" MONITOR "
650    MYTYPE(kind12)=" HORIZONTAL MONITOR "
651    MYTYPE(kind13)=" VERTICAL MONITOR "
652    MYTYPE(kind14)=" INSTRUMENT "
653    MYTYPE(kind15)=" ELECTRIC SEPTUM "
654    !               123456789012345678901234
655    MYTYPE(kind16)=" TRUE PARAELLEL  BEND  "
656    MYTYPE(kind20)=" STRAIGHT EXACT (BEND) "
657    MYTYPE(kind17)=" SOLENOID SIXTRACK"
658    !    MYTYPE(KINDFITTED)=" FITTED "
659    !    MYTYPE(KINDUSER1)=" USER_1 "
660    !    MYTYPE(KINDUSER2)=" USER_2 "
661    ind_stoc(1)='100000'
662    ind_stoc(2)='010000'
663    ind_stoc(3)='001000'
664    ind_stoc(4)='000100'
665    ind_stoc(5)='000010'
666    ind_stoc(6)='000001'
667    MADTHICK=>MADKIND2
668    MADTHIN_NORMAL=>MADKIND3N
669    MADTHIN_SKEW=>MADKIND3S
670
671    MADFAC=1.0_dp   ! to prevent overflow (David Sagan)
672    DO I=2,NMAX
673       MADFAC(I)=(I-1)*MADFAC(I-1)
674    ENDDO
675
676    !    w_p=0
677    !    w_p%nc=1
678    !    w_p%fc='(2((1X,A72,/)))'
679    !    w_p%c(1) = " RBEND ARE READ DIFFERENTLY DEPENDING ON VARIABLE MADLENGTH"
680    !    IF(MADLENGTH) THEN
681    !       w_p%c(2) = " INPUT OF CARTESIAN LENGTH FOR RBEND  "
682    !    ELSE
683    !       w_p%c(2) = " INPUT OF POLAR LENGTH FOR RBEND  "
684    !    ENDIF
685    !    ! call ! WRITE_I
686
687    NSTD=1
688    METD=2
689    electron=PARTICLE
690    if(electron) then
691       A_particle = A_ELECTRON
692    else
693       A_particle = A_PROTON
694    endif
695    !    w_p=0
696    !    w_p%nc=1
697    !    w_p%fc='(1((1X,A72)))'
698    !    IF(ELECTRON) THEN
699    !       w_p%c(1) = " THIS IS A ELECTRON (POSITRON ACTUALLY) "
700    !    ELSE
701    !       w_p%c(1) = " THIS IS A PROTON "
702    !    ENDIF
703    tilt%natural=.true.
704    tilt%tilt(0)=0.0_dp
705    do i=1,nmax
706       tilt%tilt(i)=pih/i
707    enddo
708    !  SECTOR_B AND SECTOR_NMUL FOR TYPE TEAPOT
709
710    change_sector=my_false
711
712    IF(SECTOR_NMUL>0.and.firsttime_coef) THEN
713     call alloc(e_muon_scale)
714     e_muon_scale=1.0_dp
715       !  verb=global_verbose
716       !  global_verbose=.false.
717       if(firsttime_coef.or.(.not.allocated(S_B))) then
718       
719         if(   SECTOR_NMUL==10.and. SECTOR_NMUL_MAX==10) then
720             call set_s_b
721         else
722         
723         
724         
725          !          SECTOR_B%firsttime=0   !slightly unsafe
726          ALLOCATE(S_B(SECTOR_NMUL_MAX))
727          lda_old=lda_used
728          lda_used=3000
729       !        ALLOCATE(S_B0)
730       !        S_B0%firsttime=0
731       !        call nul_coef(S_B0)
732       !        call make_coef(S_B0,SECTOR_NMUL_MAX)
733       !        call get_bend_coeff(S_B0,SECTOR_NMUL_MAX)
734          DO I=1,SECTOR_NMUL_MAX
735             !             if(i==SECTOR_NMUL_MAX)     global_verbose=.true.
736             S_B(I)%firsttime=0
737             call nul_coef(S_B(I))
738             call make_coef(S_B(I),I)
739             call get_bend_coeff(S_B(I),I)
740          ENDDO
741          lda_used=lda_old
742          call print_curv("Maxwellian_bend_for_ptc.txt")
743       endif
744
745       firsttime_coef=.FALSE.
746    endif  !  calling preset routine
747    else
748         call kill(e_muon_scale)
749         call alloc(e_muon_scale)
750     e_muon_scale=1.0_dp
751    ENDIF
752    call clear_states
753    !  global_verbose=verb
754
755  END  SUBROUTINE MAKE_STATES_0
756 
757 SUBROUTINE print_curv(filename)
758    IMPLICIT NONE
759    INTEGER I,J,nmul,mf
760    character(*) filename
761    character(255) line
762    call kanalnummer(mf,filename)
763   
764   
765    nmul=SECTOR_NMUL
766    write(mf,*) SECTOR_NMUL ,SECTOR_NMUL_MAX
767         DO J=1,S_B(NMUL)%N_MONO
768         write(line,*) "S_B(" ,NMUL, ")%I(",J,")=",S_B(NMUL)%I(J),";","S_B(" ,NMUL, ")%J(",J,")=",S_B(NMUL)%J(J),";"
769         call context(line)
770         write(mf,*) line(1:len_trim(line))
771         enddo
772    DO I=1,NMUL
773       DO J=1,S_B(NMUL)%N_MONO
774        if(S_B(NMUL)%A_X(I,J)/=0.0_dp) then
775         write(line,*) "S_B(" , NMUL , ")%A_X(" ,I,",",J, ")=",S_B(NMUL)%A_X(I,J),"e0_dp"
776         call context(line)
777         if(index(line,"E-")/=0) then
778          line(len_trim(line)-4:len_trim(line))="_DP  "
779          call context(line)
780         endif
781         write(mf,*) line(1:len_trim(line))
782        endif
783        if(S_B(NMUL)%B_X(I,J)/=0.0_dp) then
784         write(line,*) "S_B(" , NMUL , ")%B_X(" ,I,",",J, ")=",S_B(NMUL)%B_X(I,J),"e0_dp"
785         call context(line)
786         if(index(line,"E-")/=0) then
787          line(len_trim(line)-4:len_trim(line))="_DP  "
788          call context(line)
789         endif
790         write(mf,*) line(1:len_trim(line))
791        endif
792        if(S_B(NMUL)%A_y(I,J)/=0.0_dp) then
793         write(line,*) "S_B(" , NMUL , ")%A_y(" ,I,",",J, ")=",S_B(NMUL)%A_y(I,J),"e0_dp"
794         call context(line)
795         if(index(line,"E-")/=0) then
796          line(len_trim(line)-4:len_trim(line))="_DP  "
797          call context(line)
798         endif
799         write(mf,*) line(1:len_trim(line))
800        endif
801        if(S_B(NMUL)%B_y(I,J)/=0.0_dp) then
802         write(line,*) "S_B(" , NMUL , ")%B_y(" ,I,",",J, ")=",S_B(NMUL)%B_y(I,J),"e0_dp"
803         call context(line)
804         if(index(line,"E-")/=0) then
805          line(len_trim(line)-4:len_trim(line))="_DP  "
806          call context(line)
807         endif
808         write(mf,*) line(1:len_trim(line))
809        endif
810       ENDDO
811    ENDDO
812   
813    close(mf)
814
815  END SUBROUTINE print_curv
816
817  subroutine make_set_coef(b,no)
818    implicit none
819    integer no
820    integer i,j,k,a,m , ic
821
822    type(B_CYL) b
823
824    ic=1
825    b%firsttime=-100
826    allocate(b%nmul)
827    allocate(b%n_mono)
828    b%nmul=no
829    b%n_mono=((no+2-ic)*(no+1-ic))/2
830    allocate(b%i(b%n_mono),b%j(b%n_mono))
831    allocate(b%a_x(no,b%n_mono),b%a_y(no,b%n_mono))
832    allocate(b%b_x(no,b%n_mono),b%b_y(no,b%n_mono))
833    b%i=0
834    b%j=0
835    b%a_x=0.0_dp
836    b%b_x=0.0_dp
837    b%a_y=0.0_dp
838    b%b_y=0.0_dp
839  end subroutine make_set_coef
840 
841 
842  subroutine clear_states     !%nxyz
843    implicit none
844    DEFAULT=DEFAULT0
845    TOTALPATH=TOTALPATH0
846    RADIATION=RADIATION0
847    NOCAVITY=NOCAVITY0
848    ENVELOPE=ENVELOPE0
849    FRINGE=FRINGE0
850    TIME=TIME0
851    ! EXACTMIS=EXACTMIS0
852    !  exactmis%exactmis=exactmis%exactmis.and.(.not.sixtrack_compatible)
853    ONLY_4D=ONLY_4d0
854    DELTA=DELTA0
855    SPIN=SPIN0
856    MODULATION=MODULATION0
857  end subroutine clear_states  !%nxyz
858
859  subroutine print_s(S,MF)
860    implicit none
861    type (INTERNAL_STATE) S
862    INTEGER MF
863
864
865    write(mf,*) "************ State Summary ****************"
866    write(mf,'((1X,a16,1x,i4,1x,a24))' ) "MADTHICK=>KIND =", MADKIND2,Mytype(MADKIND2)
867    if(MADLENGTH) then
868       write(mf,*) ' Rectangular Bend: input cartesian length '
869    else
870       write(mf,*) ' Rectangular Bend: input arc length (rho alpha) '
871    endif
872    write(mf,'((1X,a28,1x,i4))' ) ' Default integration method ',METD
873    write(mf,'((1X,a28,1x,i4))' ) ' Default integration steps  ',NSTD
874 !   if(CAVITY_TOTALPATH==1) write(mf,'((1X,a24))' ) ' Real Pill Box Cavities '
875 !   if(CAVITY_TOTALPATH==0) write(mf,'((1X,a24))' ) ' Fake Pill Box Cavities '
876
877    If(electron) then
878       if(muon==1.0_dp)  then
879          write(mf,*)"This is an electron (positron actually if charge=1) "
880       else
881          write(mf,'((1X,a21,1x,G21.14,1x,A24))' ) "This a particle with ",muon, "times the electron mass "
882       endif
883    else
884       write(mf,*) "This is a proton "
885    endif
886    write(mf, '((1X,a20,1x,a5))' )  "      EXACT_MODEL = ", CONV(EXACT_MODEL    )
887    write(mf, '((1X,a20,1x,i4))' )  "      TOTALPATH   = ", S%TOTALPATH
888    !    write(mf, '((1X,a20,1x,a5))' )  "      EXACTMIS    = ", CONV(S%EXACTMIS    )
889    write(mf,'((1X,a20,1x,a5))' ) "      RADIATION   = ", CONV(S%RADIATION  )
890    write(mf,'((1X,a20,1x,a5))' ) "      STOCHASTIC  = ", CONV(S%STOCHASTIC  )
891    write(mf,'((1X,a20,1x,a5))' ) "      ENVELOPE    = ", CONV(S%ENVELOPE  )
892    write(mf,'((1X,a20,1x,a5))' ) "      NOCAVITY    = ", CONV(S%NOCAVITY )
893    write(mf,'((1X,a20,1x,a5))' ) "      TIME        = ", CONV(S%TIME )
894    write(mf,'((1X,a20,1x,a5))' ) "      FRINGE      = ", CONV(S%FRINGE   )
895    write(mf,'((1X,a20,1x,a5))' ) "      PARA_IN     = ", CONV(S%PARA_IN  )
896    write(mf,'((1X,a20,1x,a5))' ) "      ONLY_4D     = ", CONV(S%ONLY_4D   )
897    write(mf,'((1X,a20,1x,a5))' ) "      DELTA       = ", CONV(S%DELTA    )
898    write(mf,'((1X,a20,1x,a5))' ) "      SPIN        = ", CONV(S%SPIN    )
899    write(mf,'((1X,a20,1x,a5))' ) "      MODULATION  = ", CONV(S%MODULATION    )
900    write(mf, '((1X,a20,1x,a5))' )"      RAMPING     = "  , CONV(ramp    )
901    write(mf, '((1X,a20,1x,a5))' )"      ACCELERATE  = "  , CONV(ACCELERATE    )
902    !    write(mf,'((1X,a20,1x,I4))' ) " SPIN DIMENSION   = ", S%SPIN_DIM
903    !   ! call ! WRITE_I
904  end subroutine print_s
905
906  FUNCTION CONV(LOG)
907    IMPLICIT NONE
908    CHARACTER(5) CONV
909    logical(lp) LOG
910    CONV="FALSE"
911    IF(LOG) CONV="TRUE "
912  END FUNCTION CONV
913
914  SUBROUTINE MAKE_STATES_m(muonfactor)
915    USE   definition
916    IMPLICIT NONE
917    logical(lp) :: doneitt=.true.
918    real(dp) muonfactor,MASSF
919    CALL MAKE_YOSHIDA
920    muon=muonfactor
921    MASSF=muon*pmae
922    call MAKE_STATES_0(doneitt)
923
924    IF(ABS(MASSF-pmap)/PMAP<0.01E0_DP) THEN
925       A_PARTICLE=A_PROTON
926    ELSEIF(ABS(MASSF-pmae)/pmae<0.01E0_DP) THEN
927       A_PARTICLE=A_ELECTRON
928    ELSEIF(ABS(MASSF-pmaMUON)/pmaMUON<0.01E0_DP) THEN
929       A_PARTICLE=A_MUON
930    ENDIF
931
932  END  SUBROUTINE MAKE_STATES_m
933
934  SUBROUTINE update_STATES
935    USE   definition
936    IMPLICIT NONE
937
938    TOTALPATH=  TOTALPATH+DEFAULT
939
940    RADIATION=  RADIATION+DEFAULT
941
942    NOCAVITY=  NOCAVITY+DEFAULT
943
944    STOCHASTIC=  STOCHASTIC+DEFAULT
945
946    ENVELOPE=  ENVELOPE+DEFAULT
947
948    TIME=  TIME+DEFAULT
949
950    FRINGE=  FRINGE+DEFAULT
951
952    ONLY_4D= ONLY_4D+DEFAULT
953
954    DELTA= DELTA+DEFAULT
955
956    ! EXACTMIS= EXACTMIS+DEFAULT
957
958    SPIN= SPIN+DEFAULT
959
960    MODULATION= MODULATION+DEFAULT
961
962  END  SUBROUTINE update_STATES
963
964
965  SUBROUTINE  EQUALt(S2,S1)
966    implicit none
967    type (INTERNAL_STATE),INTENT(OUT)::S2
968    type (INTERNAL_STATE),INTENT(IN)::S1
969
970    S2%TOTALPATH=   S1%TOTALPATH
971    !   S2%EXACTMIS=       S1%EXACTMIS
972    S2%RADIATION=     S1%RADIATION
973    S2%NOCAVITY=    S1%NOCAVITY
974    S2%TIME=        S1%TIME
975    S2%FRINGE=           S1%FRINGE
976    S2%stochastic=           S1%stochastic
977    S2%ENVELOPE=           S1%ENVELOPE
978    S2%PARA_IN=     S1%PARA_IN
979    S2%ONLY_4D=      S1%ONLY_4D
980    S2%DELTA=       S1%DELTA
981    S2%SPIN=       S1%SPIN
982    S2%MODULATION=       S1%MODULATION
983    !    S2%spin_dim=       S1%spin_dim
984  END SUBROUTINE EQUALt
985
986
987
988  FUNCTION add( S1, S2 )
989    implicit none
990    TYPE (INTERNAL_STATE) add
991    TYPE (INTERNAL_STATE), INTENT (IN) :: S1, S2
992
993
994    add%TOTALPATH=0
995    if((S1%TOTALPATH==1).OR.(S2%TOTALPATH==1)) add%TOTALPATH=1
996
997    !    add%EXACT    =       S1%EXACT.OR.S2%EXACT
998    !   add%EXACTMIS    =       (S1%EXACTMIS.OR.S2%EXACTMIS).and.(.not.sixtrack_compatible)
999    add%RADIATION  =  S1%RADIATION.OR.S2%RADIATION
1000    add%NOCAVITY =  S1%NOCAVITY.OR.S2%NOCAVITY
1001    add%TIME     =  S1%TIME.OR.S2%TIME
1002    add%FRINGE   =       S1%FRINGE.OR.S2%FRINGE
1003    add%stochastic   =       S1%stochastic.OR.S2%stochastic
1004    add%ENVELOPE   =       S1%ENVELOPE.OR.S2%ENVELOPE
1005    add%ONLY_4D  =       S1%ONLY_4D.OR.S2%ONLY_4D
1006    add%DELTA  =       S1%DELTA.OR.S2%DELTA
1007    add%SPIN  =       S1%SPIN.OR.S2%SPIN
1008    add%MODULATION  =       S1%MODULATION.OR.S2%MODULATION
1009    add%PARA_IN  =       S1%PARA_IN.OR.S2%PARA_IN.or.ALWAYS_knobs
1010    !    add%SPIN_DIM  =       MAX(S1%SPIN_DIM,S2%SPIN_DIM)
1011    IF(add%stochastic) THEN
1012       add%RADIATION=T
1013    ENDIF
1014  !  IF(add%ENVELOPE) THEN
1015  !     add%radiation=T
1016  !  ENDIF
1017    IF(add%stochastic) THEN
1018       add%radiation=T
1019    ENDIF
1020    IF(add%DELTA) THEN
1021       add%ONLY_4D=T
1022       add%NOCAVITY =  T
1023    ENDIF
1024    IF(add%ONLY_4D) THEN
1025       add%TOTALPATH=  0
1026       add%RADIATION  =  F
1027       add%NOCAVITY =  T
1028       add%stochastic   =  F
1029       add%ENVELOPE   =  F
1030    ENDIF
1031  END FUNCTION add
1032
1033  FUNCTION sub( S1, S2 )
1034    implicit none
1035    TYPE (INTERNAL_STATE) sub
1036    TYPE (INTERNAL_STATE), INTENT (IN) :: S1, S2
1037    logical(lp) dum1,dum2
1038
1039    sub%TOTALPATH=0
1040    dum1=S1%TOTALPATH==1
1041    dum2=S2%TOTALPATH==1
1042    if(dum1.min.dum2) sub%TOTALPATH=1
1043
1044    !    sub%TOTALPATH=  S1%TOTALPATH.min.S2%TOTALPATH
1045
1046    !   sub%EXACTMIS    =       (S1%EXACTMIS.min.S2%EXACTMIS).and.(.not.sixtrack_compatible)
1047    sub%RADIATION  =  S1%RADIATION.min.S2%RADIATION
1048    sub%NOCAVITY =  S1%NOCAVITY.min.S2%NOCAVITY
1049    sub%TIME     =  S1%TIME.min.S2%TIME
1050    sub%FRINGE   =       S1%FRINGE.min.S2%FRINGE
1051    sub%stochastic   =       S1%stochastic.min.S2%stochastic
1052    sub%ENVELOPE   =       S1%ENVELOPE.min.S2%ENVELOPE
1053    sub%ONLY_4D  =       S1%ONLY_4D.min.S2%ONLY_4D
1054    sub%DELTA  =       S1%DELTA.min.S2%DELTA
1055    sub%SPIN  =       S1%SPIN.min.S2%SPIN
1056    sub%MODULATION  = S1%MODULATION.min.S2%MODULATION
1057    sub%PARA_IN  =       (S1%PARA_IN.MIN.S2%PARA_IN).or.ALWAYS_knobs
1058    !    sub%SPIN_DIM  =       MAX(S1%SPIN_DIM,S2%SPIN_DIM)
1059    IF(sub%stochastic) THEN
1060       sub%RADIATION=T
1061    ENDIF
1062 !   IF(sub%ENVELOPE) THEN
1063 !      sub%RADIATION=T
1064 !   ENDIF
1065    IF(sub%DELTA) THEN
1066       sub%ONLY_4D=T
1067       sub%NOCAVITY =  T
1068    ENDIF
1069    IF(sub%ONLY_4D) THEN
1070       sub%TOTALPATH=  0
1071       sub%RADIATION  =  F
1072       sub%ENVELOPE  =  F
1073       sub%stochastic  =  F
1074       sub%NOCAVITY =  T
1075       sub%stochastic   =  F
1076    ENDIF
1077
1078  END FUNCTION sub
1079
1080  FUNCTION PARA_REMA(S1)   ! UNARY +
1081    implicit none
1082    TYPE (INTERNAL_STATE) PARA_REMA
1083    TYPE (INTERNAL_STATE), INTENT (IN) :: S1
1084
1085    PARA_REMA         =    S1
1086    PARA_REMA%PARA_IN =     T
1087
1088  END FUNCTION PARA_REMA
1089
1090
1091
1092  subroutine S_init(STATE,NO1,NP1,pack,ND2,NPARA)
1093    !  subroutine S_init(STATE,NO1,NP1,PACKAGE,MAPINT,ND2,NPARA)
1094    implicit none
1095    TYPE (INTERNAL_STATE), INTENT(IN):: STATE
1096    LOGICAL(lp), optional, INTENT(IN):: pack
1097    INTEGER, INTENT(IN):: NO1,NP1
1098    INTEGER ND1,NDEL,NDPT1
1099    INTEGER,optional :: ND2,NPARA
1100    INTEGER  ND2l,NPARAl
1101    LOGICAL(lp) package
1102    call dd_p !valishev
1103    doing_ac_modulation_in_ptc=.false.
1104    package=my_true
1105    if(present(pack))     package=my_true
1106
1107
1108    NDEL=0
1109    NDPT1=0
1110    IF(STATE%NOCAVITY)  THEN
1111       IF(STATE%ONLY_4D) THEN
1112          IF(STATE%DELTA) THEN
1113             ND1=2
1114             NDEL=1
1115             !             MAPINT=5
1116          ELSE
1117             ND1=2
1118             NDEL=0
1119             !            MAPINT=4
1120          ENDIF
1121       ELSE
1122          ND1=3
1123          NDEL=0
1124          NDPT1=5  !+C_%NDPT_OTHER
1125          !         MAPINT=6
1126       ENDIF
1127    ELSE              ! CAVITY IN RING
1128       ND1=3
1129       NDPT1=0
1130       !       MAPINT=6
1131    ENDIF
1132
1133    IF(STATE%modulation)  then
1134       doing_ac_modulation_in_ptc=.true.
1135       ND1=ND1+1
1136    endif
1137
1138    !    write(6,*) NO1,ND1,NP1,NDEL,NDPT1
1139    !pause 678
1140    CALL INIT(NO1,ND1,NP1+NDEL,NDPT1,PACKAGE)
1141
1142    ND2l=ND1*2
1143    NPARAl=ND2l+NDEL
1144    C_%NPARA=NPARAl
1145    C_%ND2=ND2l
1146    C_%npara_fpp=NPARAl
1147    C_%SPIN_POS=0
1148    C_%NSPIN=0
1149
1150    if(present(nd2)) nd2=nd2l
1151    if(present(npara)) npara=nparal
1152
1153  END  subroutine S_init
1154
1155
1156  subroutine init_default(STATE,NO1,NP1)
1157    implicit none
1158    TYPE (INTERNAL_STATE),OPTIONAL, INTENT(IN):: STATE
1159    INTEGER, OPTIONAL, INTENT(IN):: NO1,NP1
1160    INTEGER   ND2,NPARA,NO2,NP2
1161    TYPE (INTERNAL_STATE)   STATE2
1162    STATE2=DEFAULT
1163    NO2=1;NP2=0;
1164    IF(PRESENT(STATE))      STATE2=STATE
1165    IF(PRESENT(NO1))      NO2=NO1
1166    IF(PRESENT(NP1))      NP2=NP1
1167    call init(STATE2,NO2,NP2,my_true,ND2,NPARA)
1168    C_%NPARA=NPARA
1169    C_%ND2=ND2
1170    C_%npara_fpp=NPARA
1171  END  subroutine init_default
1172
1173  SUBROUTINE B2PERPR(P,B,X,X5,B2)
1174    IMPLICIT NONE
1175    real(dp), INTENT(INOUT) :: B2
1176    real(dp), INTENT(IN) :: X(6),B(3),X5
1177    TYPE(MAGNET_CHART), INTENT(IN) :: P
1178    real(dp) E(3)
1179    !---  GOOD FOR TIME=FALSE
1180    E(1)=P%DIR*X(2)/(1.0_dp+X5)
1181    E(2)=P%DIR*X(4)/(1.0_dp+X5)
1182    E(3)=P%DIR*ROOT((1.0_dp+X5)**2-X(2)**2-X(4)**2)/(1.0_dp+X5)
1183
1184    B2=0.0_dp
1185    B2=(B(2)*E(3)-B(3)*E(2))**2+B2
1186    B2=(B(1)*E(2)-B(2)*E(1))**2+B2
1187    B2=(B(3)*E(1)-B(1)*E(3))**2+B2
1188
1189    RETURN
1190  END       SUBROUTINE B2PERPR
1191
1192
1193  SUBROUTINE B2PERPP(P,B,X,X5,B2)
1194    IMPLICIT NONE
1195    TYPE(REAL_8) , INTENT(INOUT) :: B2
1196    TYPE(REAL_8) , INTENT(IN) :: X(6),B(3),X5
1197    TYPE(MAGNET_CHART), INTENT(IN) :: P
1198    TYPE(REAL_8)  E(3)
1199    CALL ALLOC(E,3)
1200    !---  GOOD FOR TIME=FALSE
1201    E(1)=P%DIR*X(2)/(1.0_dp+X5)
1202    E(2)=P%DIR*X(4)/(1.0_dp+X5)
1203    E(3)=P%DIR*SQRT((1.0_dp+X5)**2-X(2)**2-X(4)**2)/(1.0_dp+X5)
1204
1205
1206    B2=0.0_dp
1207    B2=(B(2)*E(3)-B(3)*E(2))**2+B2
1208    B2=(B(1)*E(2)-B(2)*E(1))**2+B2
1209    B2=(B(3)*E(1)-B(1)*E(3))**2+B2
1210    CALL KILL(E,3)
1211    RETURN
1212  END       SUBROUTINE B2PERPP
1213
1214  SUBROUTINE DTILTR_EXTERNAL(DIR,TILTD,I,X)
1215    IMPLICIT NONE
1216    real(dp),INTENT(INOUT):: X(6)
1217    INTEGER,INTENT(IN):: I,DIR
1218    REAL(DP),INTENT(IN) :: TILTD
1219    real(dp) YS
1220
1221
1222    IF(TILTD==0.0_dp) RETURN
1223    IF(I==1) THEN
1224       ys=COS(DIR*TILTD)*x(1)+SIN(DIR*TILTD)*x(3)
1225       x(3)=COS(DIR*TILTD)*x(3)-SIN(DIR*TILTD)*x(1)
1226       x(1)=ys
1227       ys=COS(DIR*TILTD)*x(2)+SIN(DIR*TILTD)*x(4)
1228       x(4)=COS(DIR*TILTD)*x(4)-SIN(DIR*TILTD)*x(2)
1229       x(2)=ys
1230    ELSE
1231       ys=COS(DIR*TILTD)*x(1)-SIN(DIR*TILTD)*x(3)
1232       x(3)=COS(DIR*TILTD)*x(3)+SIN(DIR*TILTD)*x(1)
1233       x(1)=ys
1234       ys=COS(DIR*TILTD)*x(2)-SIN(DIR*TILTD)*x(4)
1235       x(4)=COS(DIR*TILTD)*x(4)+SIN(DIR*TILTD)*x(2)
1236       x(2)=ys
1237    ENDIF
1238
1239  END SUBROUTINE DTILTR_EXTERNAL
1240
1241  SUBROUTINE DTILTP_EXTERNAL(DIR,TILTD,I,X)
1242    IMPLICIT NONE
1243    TYPE(REAL_8),INTENT(INOUT):: X(6)
1244    INTEGER,INTENT(IN):: I,DIR
1245    REAL(DP),INTENT(IN) :: TILTD
1246    TYPE(REAL_8) YS
1247
1248    IF(TILTD==0.0_dp) RETURN
1249    CALL ALLOC(YS)
1250
1251    IF(I==1) THEN
1252       ys=COS(DIR*TILTD)*x(1)+SIN(DIR*TILTD)*x(3)
1253       x(3)=COS(DIR*TILTD)*x(3)-SIN(DIR*TILTD)*x(1)
1254       x(1)=ys
1255       ys=COS(DIR*TILTD)*x(2)+SIN(DIR*TILTD)*x(4)
1256       x(4)=COS(DIR*TILTD)*x(4)-SIN(DIR*TILTD)*x(2)
1257       x(2)=ys
1258    ELSE
1259       ys=COS(DIR*TILTD)*x(1)-SIN(DIR*TILTD)*x(3)
1260       x(3)=COS(DIR*TILTD)*x(3)+SIN(DIR*TILTD)*x(1)
1261       x(1)=ys
1262       ys=COS(DIR*TILTD)*x(2)-SIN(DIR*TILTD)*x(4)
1263       x(4)=COS(DIR*TILTD)*x(4)+SIN(DIR*TILTD)*x(2)
1264       x(2)=ys
1265    ENDIF
1266    CALL KILL(YS)
1267
1268  END SUBROUTINE DTILTP_EXTERNAL
1269
1270  SUBROUTINE dd_p   !(u,dd)    !valishev
1271    !use my_own_1D_TPSA
1272
1273    ! Tracking subroutine for elliptical lens
1274    ! A.Valishev (valishev@fnal.gov) October 19, 2010
1275    ! Modified by E. Forest for PTC
1276    implicit none
1277    type(real_8) u,dd
1278    type(my_1D_taylor) del,logdel
1279    integer i,n
1280
1281    call set_my_taylor_no(N_my_1D_taylor)
1282    del=0.0_dp
1283    del%a(1)=1.0_dp
1284    logdel=0.0_dp
1285    val_del=0.0_dp
1286    logdel=log(1.0_dp+del**2+del*sqrt(2.0_dp)*sqrt(1.0_dp+del**2/2.0_dp))
1287
1288    do i=0,N_my_1D_taylor
1289       n=2*i+1
1290       if(n>N_my_1D_taylor) cycle
1291       val_del%a(i)=logdel%a(n)
1292    enddo
1293  end SUBROUTINE dd_p  !valishev
1294
1295subroutine set_s_b
1296  implicit none
1297  integer i,s1,s2
1298 
1299  s1=SECTOR_NMUL_MAX
1300  s2=SECTOR_NMUL
1301  SECTOR_NMUL_MAX=10
1302  SECTOR_NMUL=10
1303
1304  ALLOCATE(S_B(SECTOR_NMUL_MAX))
1305           DO I=1,SECTOR_NMUL_MAX
1306             call nul_coef(S_B(I))
1307             call make_set_coef(S_B(I),I)
1308             S_B(I)%firsttime=0
1309          ENDDO
1310
1311 S_B(10)%I(1)=9;S_B(10)%J(1)=0;
1312 S_B(10)%I(2)=8;S_B(10)%J(2)=1;
1313 S_B(10)%I(3)=8;S_B(10)%J(3)=0;
1314 S_B(10)%I(4)=7;S_B(10)%J(4)=2;
1315 S_B(10)%I(5)=7;S_B(10)%J(5)=1;
1316 S_B(10)%I(6)=7;S_B(10)%J(6)=0;
1317 S_B(10)%I(7)=6;S_B(10)%J(7)=3;
1318 S_B(10)%I(8)=6;S_B(10)%J(8)=2;
1319 S_B(10)%I(9)=6;S_B(10)%J(9)=1;
1320 S_B(10)%I(10)=6;S_B(10)%J(10)=0;
1321 S_B(10)%I(11)=5;S_B(10)%J(11)=4;
1322 S_B(10)%I(12)=5;S_B(10)%J(12)=3;
1323 S_B(10)%I(13)=5;S_B(10)%J(13)=2;
1324 S_B(10)%I(14)=5;S_B(10)%J(14)=1;
1325 S_B(10)%I(15)=5;S_B(10)%J(15)=0;
1326 S_B(10)%I(16)=4;S_B(10)%J(16)=5;
1327 S_B(10)%I(17)=4;S_B(10)%J(17)=4;
1328 S_B(10)%I(18)=4;S_B(10)%J(18)=3;
1329 S_B(10)%I(19)=4;S_B(10)%J(19)=2;
1330 S_B(10)%I(20)=4;S_B(10)%J(20)=1;
1331 S_B(10)%I(21)=4;S_B(10)%J(21)=0;
1332 S_B(10)%I(22)=3;S_B(10)%J(22)=6;
1333 S_B(10)%I(23)=3;S_B(10)%J(23)=5;
1334 S_B(10)%I(24)=3;S_B(10)%J(24)=4;
1335 S_B(10)%I(25)=3;S_B(10)%J(25)=3;
1336 S_B(10)%I(26)=3;S_B(10)%J(26)=2;
1337 S_B(10)%I(27)=3;S_B(10)%J(27)=1;
1338 S_B(10)%I(28)=3;S_B(10)%J(28)=0;
1339 S_B(10)%I(29)=2;S_B(10)%J(29)=7;
1340 S_B(10)%I(30)=2;S_B(10)%J(30)=6;
1341 S_B(10)%I(31)=2;S_B(10)%J(31)=5;
1342 S_B(10)%I(32)=2;S_B(10)%J(32)=4;
1343 S_B(10)%I(33)=2;S_B(10)%J(33)=3;
1344 S_B(10)%I(34)=2;S_B(10)%J(34)=2;
1345 S_B(10)%I(35)=2;S_B(10)%J(35)=1;
1346 S_B(10)%I(36)=2;S_B(10)%J(36)=0;
1347 S_B(10)%I(37)=1;S_B(10)%J(37)=8;
1348 S_B(10)%I(38)=1;S_B(10)%J(38)=7;
1349 S_B(10)%I(39)=1;S_B(10)%J(39)=6;
1350 S_B(10)%I(40)=1;S_B(10)%J(40)=5;
1351 S_B(10)%I(41)=1;S_B(10)%J(41)=4;
1352 S_B(10)%I(42)=1;S_B(10)%J(42)=3;
1353 S_B(10)%I(43)=1;S_B(10)%J(43)=2;
1354 S_B(10)%I(44)=1;S_B(10)%J(44)=1;
1355 S_B(10)%I(45)=1;S_B(10)%J(45)=0;
1356 S_B(10)%I(46)=0;S_B(10)%J(46)=9;
1357 S_B(10)%I(47)=0;S_B(10)%J(47)=8;
1358 S_B(10)%I(48)=0;S_B(10)%J(48)=7;
1359 S_B(10)%I(49)=0;S_B(10)%J(49)=6;
1360 S_B(10)%I(50)=0;S_B(10)%J(50)=5;
1361 S_B(10)%I(51)=0;S_B(10)%J(51)=4;
1362 S_B(10)%I(52)=0;S_B(10)%J(52)=3;
1363 S_B(10)%I(53)=0;S_B(10)%J(53)=2;
1364 S_B(10)%I(54)=0;S_B(10)%J(54)=1;
1365 S_B(10)%I(55)=0;S_B(10)%J(55)=0;
1366 S_B(10)%A_Y(1,4)=-0.500000000000000E0_DP
1367 S_B(10)%A_X(1,7)=-1.16666666666667E0_DP
1368 S_B(10)%A_Y(1,8)=0.499999999999999E0_DP
1369 S_B(10)%A_Y(1,11)=2.62500000000000E0_DP
1370 S_B(10)%A_X(1,12)=0.999999999999999E0_DP
1371 S_B(10)%A_Y(1,13)=-0.500000000000000E0_DP
1372 S_B(10)%A_X(1,16)=2.62500000000000E0_DP
1373 S_B(10)%A_Y(1,17)=-1.87500000000000E0_DP
1374 S_B(10)%A_X(1,18)=-0.833333333333334E0_DP
1375 S_B(10)%A_Y(1,19)=0.500000000000000E0_DP
1376 S_B(10)%A_Y(1,22)=-2.18750000000000E0_DP
1377 S_B(10)%A_X(1,23)=-1.50000000000000E0_DP
1378 S_B(10)%A_Y(1,24)=1.25000000000000E0_DP
1379 S_B(10)%A_X(1,25)=0.666666666666667E0_DP
1380 S_B(10)%A_Y(1,26)=-0.500000000000000E0_DP
1381 S_B(10)%A_X(1,29)=-0.937500000000000E0_DP
1382 S_B(10)%A_Y(1,30)=0.937500000000000E0_DP
1383 S_B(10)%A_X(1,31)=0.750000000000000E0_DP
1384 S_B(10)%A_Y(1,32)=-0.750000000000000E0_DP
1385 S_B(10)%A_X(1,33)=-0.500000000000000E0_DP
1386 S_B(10)%A_Y(1,34)=0.500000000000000E0_DP
1387 S_B(10)%A_Y(1,37)=0.273437500000000E0_DP
1388 S_B(10)%A_X(1,38)=0.267857142857143E0_DP
1389 S_B(10)%A_Y(1,39)=-0.312500000000000E0_DP
1390 S_B(10)%A_X(1,40)=-0.300000000000000E0_DP
1391 S_B(10)%A_Y(1,41)=0.375000000000000E0_DP
1392 S_B(10)%A_X(1,42)=0.333333333333333E0_DP
1393 S_B(10)%A_Y(1,43)=-0.500000000000000E0_DP
1394 S_B(10)%B_X(1,45)=-1.00000000000000E0_DP
1395 S_B(10)%A_Y(1,45)=1.00000000000000E0_DP
1396 S_B(10)%A_X(1,46)=3.038194444444443E-002_DP
1397 S_B(10)%A_Y(1,47)=-3.906249999999999E-002_DP
1398 S_B(10)%A_X(1,48)=-4.464285714285716E-002_DP
1399 S_B(10)%A_Y(1,49)=6.250000000000000E-002_DP
1400 S_B(10)%A_X(1,50)=7.500000000000000E-002_DP
1401 S_B(10)%A_Y(1,51)=-0.125000000000000E0_DP
1402 S_B(10)%A_X(1,52)=-0.166666666666667E0_DP
1403 S_B(10)%A_Y(1,53)=0.500000000000000E0_DP
1404 S_B(10)%A_X(1,54)=1.00000000000000E0_DP
1405 S_B(10)%B_X(1,55)=-1.00000000000000E0_DP
1406 S_B(10)%A_Y(1,55)=1.00000000000000E0_DP
1407 S_B(10)%A_Y(2,4)=0.500000000000001E0_DP
1408 S_B(10)%A_X(2,7)=1.16666666666667E0_DP
1409 S_B(10)%B_Y(2,7)=0.166666666666667E0_DP
1410 S_B(10)%A_Y(2,8)=-0.499999999999998E0_DP
1411 S_B(10)%B_X(2,11)=0.250000000000001E0_DP
1412 S_B(10)%A_Y(2,11)=-2.62500000000000E0_DP
1413 S_B(10)%A_X(2,12)=-0.999999999999996E0_DP
1414 S_B(10)%B_Y(2,12)=-0.166666666666667E0_DP
1415 S_B(10)%A_Y(2,13)=0.500000000000001E0_DP
1416 S_B(10)%A_X(2,16)=-2.62500000000000E0_DP
1417 S_B(10)%B_Y(2,16)=-0.375000000000001E0_DP
1418 S_B(10)%B_X(2,17)=-0.208333333333333E0_DP
1419 S_B(10)%A_Y(2,17)=1.87500000000000E0_DP
1420 S_B(10)%A_X(2,18)=0.833333333333335E0_DP
1421 S_B(10)%B_Y(2,18)=0.166666666666666E0_DP
1422 S_B(10)%A_Y(2,19)=-0.500000000000001E0_DP
1423 S_B(10)%B_X(2,22)=-0.250000000000001E0_DP
1424 S_B(10)%A_Y(2,22)=2.18750000000000E0_DP
1425 S_B(10)%A_X(2,23)=1.50000000000000E0_DP
1426 S_B(10)%B_Y(2,23)=0.250000000000000E0_DP
1427 S_B(10)%B_X(2,24)=0.166666666666666E0_DP
1428 S_B(10)%A_Y(2,24)=-1.25000000000000E0_DP
1429 S_B(10)%A_X(2,25)=-0.666666666666667E0_DP
1430 S_B(10)%B_Y(2,25)=-0.166666666666667E0_DP
1431 S_B(10)%A_Y(2,26)=0.499999999999999E0_DP
1432 S_B(10)%A_X(2,29)=0.937500000000000E0_DP
1433 S_B(10)%B_Y(2,29)=0.133928571428572E0_DP
1434 S_B(10)%B_X(2,30)=0.125000000000000E0_DP
1435 S_B(10)%A_Y(2,30)=-0.937499999999998E0_DP
1436 S_B(10)%A_X(2,31)=-0.750000000000001E0_DP
1437 S_B(10)%B_Y(2,31)=-0.150000000000000E0_DP
1438 S_B(10)%B_X(2,32)=-0.125000000000001E0_DP
1439 S_B(10)%A_Y(2,32)=0.750000000000001E0_DP
1440 S_B(10)%A_X(2,33)=0.499999999999999E0_DP
1441 S_B(10)%B_Y(2,33)=0.166666666666667E0_DP
1442 S_B(10)%A_Y(2,34)=-0.500000000000000E0_DP
1443 S_B(10)%B_X(2,36)=-1.00000000000000E0_DP
1444 S_B(10)%A_Y(2,36)=1.00000000000000E0_DP
1445 S_B(10)%B_X(2,37)=3.348214285714293E-002_DP
1446 S_B(10)%A_Y(2,37)=-0.273437500000000E0_DP
1447 S_B(10)%A_X(2,38)=-0.267857142857142E0_DP
1448 S_B(10)%B_Y(2,38)=-4.464285714285710E-002_DP
1449 S_B(10)%B_X(2,39)=-4.999999999999985E-002_DP
1450 S_B(10)%A_Y(2,39)=0.312500000000000E0_DP
1451 S_B(10)%A_X(2,40)=0.300000000000000E0_DP
1452 S_B(10)%B_Y(2,40)=7.500000000000023E-002_DP
1453 S_B(10)%B_X(2,41)=8.333333333333333E-002_DP
1454 S_B(10)%A_Y(2,41)=-0.375000000000000E0_DP
1455 S_B(10)%A_X(2,42)=-0.333333333333333E0_DP
1456 S_B(10)%B_Y(2,42)=-0.166666666666667E0_DP
1457 S_B(10)%A_Y(2,43)=0.500000000000000E0_DP
1458 S_B(10)%A_X(2,44)=2.00000000000000E0_DP
1459 S_B(10)%B_Y(2,44)=1.00000000000000E0_DP
1460 S_B(10)%B_X(2,45)=-1.00000000000000E0_DP
1461 S_B(10)%A_Y(2,45)=1.00000000000000E0_DP
1462 S_B(10)%A_X(2,46)=-3.038194444444445E-002_DP
1463 S_B(10)%B_Y(2,46)=-4.340277777777785E-003_DP
1464 S_B(10)%B_X(2,47)=-5.580357142857138E-003_DP
1465 S_B(10)%A_Y(2,47)=3.906249999999994E-002_DP
1466 S_B(10)%A_X(2,48)=4.464285714285721E-002_DP
1467 S_B(10)%B_Y(2,48)=8.928571428571412E-003_DP
1468 S_B(10)%B_X(2,49)=1.250000000000004E-002_DP
1469 S_B(10)%A_Y(2,49)=-6.250000000000003E-002_DP
1470 S_B(10)%A_X(2,50)=-7.499999999999996E-002_DP
1471 S_B(10)%B_Y(2,50)=-2.500000000000000E-002_DP
1472 S_B(10)%B_X(2,51)=-4.166666666666665E-002_DP
1473 S_B(10)%A_Y(2,51)=0.125000000000000E0_DP
1474 S_B(10)%A_X(2,52)=0.166666666666667E0_DP
1475 S_B(10)%B_Y(2,52)=0.166666666666667E0_DP
1476 S_B(10)%B_X(2,53)=0.500000000000000E0_DP
1477 S_B(10)%A_Y(2,53)=-0.500000000000000E0_DP
1478 S_B(10)%A_X(2,54)=1.00000000000000E0_DP
1479 S_B(10)%B_Y(2,54)=1.00000000000000E0_DP
1480 S_B(10)%A_Y(3,4)=-0.499999999999996E0_DP
1481 S_B(10)%A_X(3,7)=-1.16666666666666E0_DP
1482 S_B(10)%B_Y(3,7)=-0.333333333333335E0_DP
1483 S_B(10)%A_Y(3,8)=0.500000000000001E0_DP
1484 S_B(10)%B_X(3,11)=-0.500000000000002E0_DP
1485 S_B(10)%A_Y(3,11)=2.74999999999999E0_DP
1486 S_B(10)%A_X(3,12)=1.00000000000000E0_DP
1487 S_B(10)%B_Y(3,12)=0.333333333333334E0_DP
1488 S_B(10)%A_Y(3,13)=-0.500000000000000E0_DP
1489 S_B(10)%A_X(3,16)=2.74999999999999E0_DP
1490 S_B(10)%B_Y(3,16)=0.750000000000003E0_DP
1491 S_B(10)%B_X(3,17)=0.416666666666667E0_DP
1492 S_B(10)%A_Y(3,17)=-2.00000000000000E0_DP
1493 S_B(10)%A_X(3,18)=-0.833333333333334E0_DP
1494 S_B(10)%B_Y(3,18)=-0.333333333333336E0_DP
1495 S_B(10)%A_Y(3,19)=0.499999999999999E0_DP
1496 S_B(10)%B_X(3,22)=0.500000000000002E0_DP
1497 S_B(10)%A_Y(3,22)=-2.31249999999999E0_DP
1498 S_B(10)%A_X(3,23)=-1.60000000000000E0_DP
1499 S_B(10)%B_Y(3,23)=-0.500000000000001E0_DP
1500 S_B(10)%B_X(3,24)=-0.333333333333336E0_DP
1501 S_B(10)%A_Y(3,24)=1.37500000000000E0_DP
1502 S_B(10)%A_X(3,25)=0.666666666666665E0_DP
1503 S_B(10)%B_Y(3,25)=0.333333333333333E0_DP
1504 S_B(10)%A_Y(3,26)=-0.500000000000000E0_DP
1505 S_B(10)%B_X(3,28)=-1.00000000000000E0_DP
1506 S_B(10)%A_Y(3,28)=1.00000000000000E0_DP
1507 S_B(10)%A_X(3,29)=-0.991071428571425E0_DP
1508 S_B(10)%B_Y(3,29)=-0.267857142857144E0_DP
1509 S_B(10)%B_X(3,30)=-0.250000000000000E0_DP
1510 S_B(10)%A_Y(3,30)=1.01250000000000E0_DP
1511 S_B(10)%A_X(3,31)=0.825000000000000E0_DP
1512 S_B(10)%B_Y(3,31)=0.300000000000001E0_DP
1513 S_B(10)%B_X(3,32)=0.249999999999999E0_DP
1514 S_B(10)%A_Y(3,32)=-0.874999999999999E0_DP
1515 S_B(10)%A_X(3,33)=-0.500000000000000E0_DP
1516 S_B(10)%B_Y(3,33)=-0.333333333333334E0_DP
1517 S_B(10)%A_Y(3,34)=0.500000000000000E0_DP
1518 S_B(10)%A_X(3,35)=3.00000000000000E0_DP
1519 S_B(10)%B_Y(3,35)=2.00000000000000E0_DP
1520 S_B(10)%B_X(3,36)=-1.00000000000000E0_DP
1521 S_B(10)%A_Y(3,36)=1.00000000000000E0_DP
1522 S_B(10)%B_X(3,37)=-6.696428571428595E-002_DP
1523 S_B(10)%A_Y(3,37)=0.290178571428571E0_DP
1524 S_B(10)%A_X(3,38)=0.289285714285715E0_DP
1525 S_B(10)%B_Y(3,38)=8.928571428571447E-002_DP
1526 S_B(10)%B_X(3,39)=0.100000000000000E0_DP
1527 S_B(10)%A_Y(3,39)=-0.350000000000000E0_DP
1528 S_B(10)%A_X(3,40)=-0.349999999999999E0_DP
1529 S_B(10)%B_Y(3,40)=-0.150000000000000E0_DP
1530 S_B(10)%B_X(3,41)=-0.166666666666667E0_DP
1531 S_B(10)%A_Y(3,41)=0.500000000000000E0_DP
1532 S_B(10)%A_X(3,42)=0.333333333333333E0_DP
1533 S_B(10)%B_Y(3,42)=0.333333333333333E0_DP
1534 S_B(10)%B_X(3,43)=2.00000000000000E0_DP
1535 S_B(10)%A_Y(3,43)=-2.00000000000000E0_DP
1536 S_B(10)%A_X(3,44)=2.00000000000000E0_DP
1537 S_B(10)%B_Y(3,44)=2.00000000000000E0_DP
1538 S_B(10)%A_X(3,46)=3.224206349206341E-002_DP
1539 S_B(10)%B_Y(3,46)=8.680555555555582E-003_DP
1540 S_B(10)%B_X(3,47)=1.116071428571431E-002_DP
1541 S_B(10)%A_Y(3,47)=-4.241071428571432E-002_DP
1542 S_B(10)%A_X(3,48)=-5.000000000000000E-002_DP
1543 S_B(10)%B_Y(3,48)=-1.785714285714292E-002_DP
1544 S_B(10)%B_X(3,49)=-2.499999999999997E-002_DP
1545 S_B(10)%A_Y(3,49)=7.499999999999990E-002_DP
1546 S_B(10)%A_X(3,50)=0.100000000000000E0_DP
1547 S_B(10)%B_Y(3,50)=5.000000000000010E-002_DP
1548 S_B(10)%B_X(3,51)=8.333333333333336E-002_DP
1549 S_B(10)%A_Y(3,51)=-0.250000000000000E0_DP
1550 S_B(10)%A_X(3,52)=-0.666666666666667E0_DP
1551 S_B(10)%B_Y(3,52)=-0.333333333333333E0_DP
1552 S_B(10)%B_X(3,53)=1.00000000000000E0_DP
1553 S_B(10)%A_Y(3,53)=-1.00000000000000E0_DP
1554 S_B(10)%A_Y(4,4)=0.499999999999999E0_DP
1555 S_B(10)%A_X(4,7)=1.16666666666667E0_DP
1556 S_B(10)%B_Y(4,7)=0.499999999999997E0_DP
1557 S_B(10)%A_Y(4,8)=-0.500000000000000E0_DP
1558 S_B(10)%B_X(4,11)=0.749999999999995E0_DP
1559 S_B(10)%A_Y(4,11)=-3.00000000000000E0_DP
1560 S_B(10)%A_X(4,12)=-1.00000000000000E0_DP
1561 S_B(10)%B_Y(4,12)=-0.499999999999999E0_DP
1562 S_B(10)%A_Y(4,13)=0.499999999999999E0_DP
1563 S_B(10)%A_X(4,16)=-3.00000000000000E0_DP
1564 S_B(10)%B_Y(4,16)=-1.19999999999999E0_DP
1565 S_B(10)%B_X(4,17)=-0.624999999999999E0_DP
1566 S_B(10)%A_Y(4,17)=2.25000000000000E0_DP
1567 S_B(10)%A_X(4,18)=0.833333333333331E0_DP
1568 S_B(10)%B_Y(4,18)=0.499999999999999E0_DP
1569 S_B(10)%A_Y(4,19)=-0.500000000000001E0_DP
1570 S_B(10)%B_X(4,21)=-1.00000000000000E0_DP
1571 S_B(10)%A_Y(4,21)=1.00000000000000E0_DP
1572 S_B(10)%B_X(4,22)=-0.799999999999996E0_DP
1573 S_B(10)%A_Y(4,22)=2.56250000000000E0_DP
1574 S_B(10)%A_X(4,23)=1.80000000000000E0_DP
1575 S_B(10)%B_Y(4,23)=0.824999999999999E0_DP
1576 S_B(10)%B_X(4,24)=0.499999999999999E0_DP
1577 S_B(10)%A_Y(4,24)=-1.62500000000000E0_DP
1578 S_B(10)%A_X(4,25)=-0.666666666666668E0_DP
1579 S_B(10)%B_Y(4,25)=-0.500000000000000E0_DP
1580 S_B(10)%A_Y(4,26)=0.500000000000000E0_DP
1581 S_B(10)%A_X(4,27)=4.00000000000000E0_DP
1582 S_B(10)%B_Y(4,27)=3.00000000000000E0_DP
1583 S_B(10)%B_X(4,28)=-1.00000000000000E0_DP
1584 S_B(10)%A_Y(4,28)=1.00000000000000E0_DP
1585 S_B(10)%A_X(4,29)=1.09821428571429E0_DP
1586 S_B(10)%B_Y(4,29)=0.433928571428570E0_DP
1587 S_B(10)%B_X(4,30)=0.412500000000000E0_DP
1588 S_B(10)%A_Y(4,30)=-1.16250000000000E0_DP
1589 S_B(10)%A_X(4,31)=-0.974999999999999E0_DP
1590 S_B(10)%B_Y(4,31)=-0.524999999999999E0_DP
1591 S_B(10)%B_X(4,32)=-0.375000000000000E0_DP
1592 S_B(10)%A_Y(4,32)=1.12500000000000E0_DP
1593 S_B(10)%A_X(4,33)=0.500000000000000E0_DP
1594 S_B(10)%B_Y(4,33)=0.499999999999998E0_DP
1595 S_B(10)%B_X(4,34)=4.50000000000000E0_DP
1596 S_B(10)%A_Y(4,34)=-4.50000000000000E0_DP
1597 S_B(10)%A_X(4,35)=3.00000000000000E0_DP
1598 S_B(10)%B_Y(4,35)=3.00000000000000E0_DP
1599 S_B(10)%B_X(4,37)=0.108482142857142E0_DP
1600 S_B(10)%A_Y(4,37)=-0.323660714285714E0_DP
1601 S_B(10)%A_X(4,38)=-0.332142857142857E0_DP
1602 S_B(10)%B_Y(4,38)=-0.150000000000000E0_DP
1603 S_B(10)%B_X(4,39)=-0.175000000000000E0_DP
1604 S_B(10)%A_Y(4,39)=0.425000000000000E0_DP
1605 S_B(10)%A_X(4,40)=0.450000000000000E0_DP
1606 S_B(10)%B_Y(4,40)=0.300000000000000E0_DP
1607 S_B(10)%B_X(4,41)=0.249999999999999E0_DP
1608 S_B(10)%A_Y(4,41)=-0.750000000000000E0_DP
1609 S_B(10)%A_X(4,42)=-3.00000000000000E0_DP
1610 S_B(10)%B_Y(4,42)=-2.00000000000000E0_DP
1611 S_B(10)%B_X(4,43)=3.00000000000000E0_DP
1612 S_B(10)%A_Y(4,43)=-3.00000000000000E0_DP
1613 S_B(10)%A_X(4,46)=-3.596230158730157E-002_DP
1614 S_B(10)%B_Y(4,46)=-1.413690476190471E-002_DP
1615 S_B(10)%B_X(4,47)=-1.874999999999998E-002_DP
1616 S_B(10)%A_Y(4,47)=4.910714285714286E-002_DP
1617 S_B(10)%A_X(4,48)=6.071428571428568E-002_DP
1618 S_B(10)%B_Y(4,48)=3.214285714285710E-002_DP
1619 S_B(10)%B_X(4,49)=4.999999999999998E-002_DP
1620 S_B(10)%A_Y(4,49)=-0.100000000000000E0_DP
1621 S_B(10)%A_X(4,50)=-0.150000000000000E0_DP
1622 S_B(10)%B_Y(4,50)=-0.150000000000000E0_DP
1623 S_B(10)%B_X(4,51)=-0.500000000000000E0_DP
1624 S_B(10)%A_Y(4,51)=0.500000000000000E0_DP
1625 S_B(10)%A_X(4,52)=-1.00000000000000E0_DP
1626 S_B(10)%B_Y(4,52)=-1.00000000000000E0_DP
1627 S_B(10)%A_Y(5,4)=-0.500000000000002E0_DP
1628 S_B(10)%A_X(5,7)=-1.16666666666667E0_DP
1629 S_B(10)%B_Y(5,7)=-0.666666666666668E0_DP
1630 S_B(10)%A_Y(5,8)=0.500000000000000E0_DP
1631 S_B(10)%B_X(5,11)=-1.00000000000000E0_DP
1632 S_B(10)%A_Y(5,11)=3.37500000000001E0_DP
1633 S_B(10)%A_X(5,12)=1.00000000000000E0_DP
1634 S_B(10)%B_Y(5,12)=0.666666666666664E0_DP
1635 S_B(10)%A_Y(5,13)=-0.500000000000000E0_DP
1636 S_B(10)%B_X(5,15)=-1.00000000000000E0_DP
1637 S_B(10)%A_Y(5,15)=1.00000000000000E0_DP
1638 S_B(10)%A_X(5,16)=3.37500000000001E0_DP
1639 S_B(10)%B_Y(5,16)=1.80000000000000E0_DP
1640 S_B(10)%B_X(5,17)=0.833333333333330E0_DP
1641 S_B(10)%A_Y(5,17)=-2.62500000000000E0_DP
1642 S_B(10)%A_X(5,18)=-0.833333333333334E0_DP
1643 S_B(10)%B_Y(5,18)=-0.666666666666667E0_DP
1644 S_B(10)%A_Y(5,19)=0.500000000000000E0_DP
1645 S_B(10)%A_X(5,20)=5.00000000000000E0_DP
1646 S_B(10)%B_Y(5,20)=4.00000000000000E0_DP
1647 S_B(10)%B_X(5,21)=-1.00000000000000E0_DP
1648 S_B(10)%A_Y(5,21)=1.00000000000000E0_DP
1649 S_B(10)%B_X(5,22)=1.20000000000000E0_DP
1650 S_B(10)%A_Y(5,22)=-3.00000000000001E0_DP
1651 S_B(10)%A_X(5,23)=-2.10000000000000E0_DP
1652 S_B(10)%B_Y(5,23)=-1.30000000000000E0_DP
1653 S_B(10)%B_X(5,24)=-0.666666666666667E0_DP
1654 S_B(10)%A_Y(5,24)=2.00000000000000E0_DP
1655 S_B(10)%A_X(5,25)=0.666666666666667E0_DP
1656 S_B(10)%B_Y(5,25)=0.666666666666667E0_DP
1657 S_B(10)%B_X(5,26)=8.00000000000000E0_DP
1658 S_B(10)%A_Y(5,26)=-8.00000000000000E0_DP
1659 S_B(10)%A_X(5,27)=4.00000000000000E0_DP
1660 S_B(10)%B_Y(5,27)=4.00000000000000E0_DP
1661 S_B(10)%A_X(5,29)=-1.28571428571429E0_DP
1662 S_B(10)%B_Y(5,29)=-0.664285714285714E0_DP
1663 S_B(10)%B_X(5,30)=-0.649999999999999E0_DP
1664 S_B(10)%A_Y(5,30)=1.45000000000000E0_DP
1665 S_B(10)%A_X(5,31)=1.20000000000000E0_DP
1666 S_B(10)%B_Y(5,31)=0.900000000000000E0_DP
1667 S_B(10)%B_X(5,32)=0.500000000000000E0_DP
1668 S_B(10)%A_Y(5,32)=-1.50000000000000E0_DP
1669 S_B(10)%A_X(5,33)=-8.00000000000000E0_DP
1670 S_B(10)%B_Y(5,33)=-6.00000000000000E0_DP
1671 S_B(10)%B_X(5,34)=6.00000000000000E0_DP
1672 S_B(10)%A_Y(5,34)=-6.00000000000000E0_DP
1673 S_B(10)%B_X(5,37)=-0.166071428571429E0_DP
1674 S_B(10)%A_Y(5,37)=0.383928571428572E0_DP
1675 S_B(10)%A_X(5,38)=0.414285714285714E0_DP
1676 S_B(10)%B_Y(5,38)=0.242857142857143E0_DP
1677 S_B(10)%B_X(5,39)=0.300000000000000E0_DP
1678 S_B(10)%A_Y(5,39)=-0.600000000000000E0_DP
1679 S_B(10)%A_X(5,40)=-0.600000000000000E0_DP
1680 S_B(10)%B_Y(5,40)=-0.600000000000000E0_DP
1681 S_B(10)%B_X(5,41)=-3.00000000000000E0_DP
1682 S_B(10)%A_Y(5,41)=3.00000000000000E0_DP
1683 S_B(10)%A_X(5,42)=-4.00000000000000E0_DP
1684 S_B(10)%B_Y(5,42)=-4.00000000000000E0_DP
1685 S_B(10)%A_X(5,46)=4.265873015873022E-002_DP
1686 S_B(10)%B_Y(5,46)=2.182539682539682E-002_DP
1687 S_B(10)%B_X(5,47)=3.035714285714282E-002_DP
1688 S_B(10)%A_Y(5,47)=-6.250000000000000E-002_DP
1689 S_B(10)%A_X(5,48)=-8.571428571428574E-002_DP
1690 S_B(10)%B_Y(5,48)=-5.714285714285715E-002_DP
1691 S_B(10)%B_X(5,49)=-0.100000000000000E0_DP
1692 S_B(10)%A_Y(5,49)=0.200000000000000E0_DP
1693 S_B(10)%A_X(5,50)=0.600000000000000E0_DP
1694 S_B(10)%B_Y(5,50)=0.400000000000000E0_DP
1695 S_B(10)%B_X(5,51)=-1.00000000000000E0_DP
1696 S_B(10)%A_Y(5,51)=1.00000000000000E0_DP
1697 S_B(10)%A_Y(6,4)=0.499999999999998E0_DP
1698 S_B(10)%A_X(6,7)=1.16666666666666E0_DP
1699 S_B(10)%B_Y(6,7)=0.833333333333330E0_DP
1700 S_B(10)%A_Y(6,8)=-0.500000000000004E0_DP
1701 S_B(10)%B_X(6,10)=-1.00000000000000E0_DP
1702 S_B(10)%A_Y(6,10)=1.00000000000000E0_DP
1703 S_B(10)%B_X(6,11)=1.24999999999999E0_DP
1704 S_B(10)%A_Y(6,11)=-3.87499999999999E0_DP
1705 S_B(10)%A_X(6,12)=-1.00000000000001E0_DP
1706 S_B(10)%B_Y(6,12)=-0.833333333333331E0_DP
1707 S_B(10)%A_Y(6,13)=0.499999999999999E0_DP
1708 S_B(10)%A_X(6,14)=6.00000000000000E0_DP
1709 S_B(10)%B_Y(6,14)=5.00000000000000E0_DP
1710 S_B(10)%B_X(6,15)=-1.00000000000000E0_DP
1711 S_B(10)%A_Y(6,15)=1.00000000000000E0_DP
1712 S_B(10)%A_X(6,16)=-3.87499999999999E0_DP
1713 S_B(10)%B_Y(6,16)=-2.62499999999999E0_DP
1714 S_B(10)%B_X(6,17)=-1.04166666666666E0_DP
1715 S_B(10)%A_Y(6,17)=3.12500000000001E0_DP
1716 S_B(10)%A_X(6,18)=0.833333333333331E0_DP
1717 S_B(10)%B_Y(6,18)=0.833333333333337E0_DP
1718 S_B(10)%B_X(6,19)=12.5000000000000E0_DP
1719 S_B(10)%A_Y(6,19)=-12.5000000000000E0_DP
1720 S_B(10)%A_X(6,20)=5.00000000000000E0_DP
1721 S_B(10)%B_Y(6,20)=5.00000000000000E0_DP
1722 S_B(10)%B_X(6,22)=-1.75000000000000E0_DP
1723 S_B(10)%A_Y(6,22)=3.75000000000000E0_DP
1724 S_B(10)%A_X(6,23)=2.50000000000001E0_DP
1725 S_B(10)%B_Y(6,23)=2.00000000000000E0_DP
1726 S_B(10)%B_X(6,24)=0.833333333333337E0_DP
1727 S_B(10)%A_Y(6,24)=-2.50000000000000E0_DP
1728 S_B(10)%A_X(6,25)=-16.6666666666667E0_DP
1729 S_B(10)%B_Y(6,25)=-13.3333333333333E0_DP
1730 S_B(10)%B_X(6,26)=10.0000000000000E0_DP
1731 S_B(10)%A_Y(6,26)=-10.0000000000000E0_DP
1732 S_B(10)%A_X(6,29)=1.60714285714286E0_DP
1733 S_B(10)%B_Y(6,29)=1.03571428571428E0_DP
1734 S_B(10)%B_X(6,30)=0.999999999999999E0_DP
1735 S_B(10)%A_Y(6,30)=-2.00000000000000E0_DP
1736 S_B(10)%A_X(6,31)=-1.50000000000000E0_DP
1737 S_B(10)%B_Y(6,31)=-1.50000000000000E0_DP
1738 S_B(10)%B_X(6,32)=-10.0000000000000E0_DP
1739 S_B(10)%A_Y(6,32)=10.0000000000000E0_DP
1740 S_B(10)%A_X(6,33)=-10.0000000000000E0_DP
1741 S_B(10)%B_Y(6,33)=-10.0000000000000E0_DP
1742 S_B(10)%B_X(6,37)=0.258928571428571E0_DP
1743 S_B(10)%A_Y(6,37)=-0.491071428571428E0_DP
1744 S_B(10)%A_X(6,38)=-0.571428571428572E0_DP
1745 S_B(10)%B_Y(6,38)=-0.428571428571428E0_DP
1746 S_B(10)%B_X(6,39)=-0.500000000000001E0_DP
1747 S_B(10)%A_Y(6,39)=1.00000000000000E0_DP
1748 S_B(10)%A_X(6,40)=4.00000000000000E0_DP
1749 S_B(10)%B_Y(6,40)=3.00000000000000E0_DP
1750 S_B(10)%B_X(6,41)=-5.00000000000000E0_DP
1751 S_B(10)%A_Y(6,41)=5.00000000000000E0_DP
1752 S_B(10)%A_X(6,46)=-5.456349206349204E-002_DP
1753 S_B(10)%B_Y(6,46)=-3.472222222222217E-002_DP
1754 S_B(10)%B_X(6,47)=-5.357142857142854E-002_DP
1755 S_B(10)%A_Y(6,47)=8.928571428571441E-002_DP
1756 S_B(10)%A_X(6,48)=0.142857142857143E0_DP
1757 S_B(10)%B_Y(6,48)=0.142857142857143E0_DP
1758 S_B(10)%B_X(6,49)=0.500000000000000E0_DP
1759 S_B(10)%A_Y(6,49)=-0.500000000000000E0_DP
1760 S_B(10)%A_X(6,50)=1.00000000000000E0_DP
1761 S_B(10)%B_Y(6,50)=1.00000000000000E0_DP
1762 S_B(10)%A_Y(7,4)=-0.499999999999998E0_DP
1763 S_B(10)%B_X(7,6)=-1.00000000000000E0_DP
1764 S_B(10)%A_Y(7,6)=1.00000000000000E0_DP
1765 S_B(10)%A_X(7,7)=-1.16666666666666E0_DP
1766 S_B(10)%B_Y(7,7)=-0.999999999999999E0_DP
1767 S_B(10)%A_Y(7,8)=0.499999999999999E0_DP
1768 S_B(10)%A_X(7,9)=7.00000000000000E0_DP
1769 S_B(10)%B_Y(7,9)=6.00000000000000E0_DP
1770 S_B(10)%B_X(7,10)=-1.00000000000000E0_DP
1771 S_B(10)%A_Y(7,10)=1.00000000000000E0_DP
1772 S_B(10)%B_X(7,11)=-1.50000000000000E0_DP
1773 S_B(10)%A_Y(7,11)=4.49999999999999E0_DP
1774 S_B(10)%A_X(7,12)=0.999999999999998E0_DP
1775 S_B(10)%B_Y(7,12)=1.00000000000000E0_DP
1776 S_B(10)%B_X(7,13)=18.0000000000000E0_DP
1777 S_B(10)%A_Y(7,13)=-18.0000000000000E0_DP
1778 S_B(10)%A_X(7,14)=6.00000000000000E0_DP
1779 S_B(10)%B_Y(7,14)=6.00000000000000E0_DP
1780 S_B(10)%A_X(7,16)=4.49999999999999E0_DP
1781 S_B(10)%B_Y(7,16)=3.75000000000000E0_DP
1782 S_B(10)%B_X(7,17)=1.25000000000000E0_DP
1783 S_B(10)%A_Y(7,17)=-3.75000000000000E0_DP
1784 S_B(10)%A_X(7,18)=-30.0000000000000E0_DP
1785 S_B(10)%B_Y(7,18)=-25.0000000000000E0_DP
1786 S_B(10)%B_X(7,19)=15.0000000000000E0_DP
1787 S_B(10)%A_Y(7,19)=-15.0000000000000E0_DP
1788 S_B(10)%B_X(7,22)=2.50000000000000E0_DP
1789 S_B(10)%A_Y(7,22)=-4.99999999999999E0_DP
1790 S_B(10)%A_X(7,23)=-3.00000000000000E0_DP
1791 S_B(10)%B_Y(7,23)=-3.00000000000000E0_DP
1792 S_B(10)%B_X(7,24)=-25.0000000000000E0_DP
1793 S_B(10)%A_Y(7,24)=25.0000000000000E0_DP
1794 S_B(10)%A_X(7,25)=-20.0000000000000E0_DP
1795 S_B(10)%B_Y(7,25)=-20.0000000000000E0_DP
1796 S_B(10)%A_X(7,29)=-2.14285714285714E0_DP
1797 S_B(10)%B_Y(7,29)=-1.71428571428571E0_DP
1798 S_B(10)%B_X(7,30)=-1.50000000000000E0_DP
1799 S_B(10)%A_Y(7,30)=3.00000000000000E0_DP
1800 S_B(10)%A_X(7,31)=15.0000000000000E0_DP
1801 S_B(10)%B_Y(7,31)=12.0000000000000E0_DP
1802 S_B(10)%B_X(7,32)=-15.0000000000000E0_DP
1803 S_B(10)%A_Y(7,32)=15.0000000000000E0_DP
1804 S_B(10)%B_X(7,37)=-0.428571428571428E0_DP
1805 S_B(10)%A_Y(7,37)=0.714285714285714E0_DP
1806 S_B(10)%A_X(7,38)=0.857142857142857E0_DP
1807 S_B(10)%B_Y(7,38)=0.857142857142857E0_DP
1808 S_B(10)%B_X(7,39)=4.00000000000000E0_DP
1809 S_B(10)%A_Y(7,39)=-4.00000000000000E0_DP
1810 S_B(10)%A_X(7,40)=6.00000000000000E0_DP
1811 S_B(10)%B_Y(7,40)=6.00000000000000E0_DP
1812 S_B(10)%A_X(7,46)=7.936507936507931E-002_DP
1813 S_B(10)%B_Y(7,46)=5.952380952380952E-002_DP
1814 S_B(10)%B_X(7,47)=0.107142857142857E0_DP
1815 S_B(10)%A_Y(7,47)=-0.178571428571429E0_DP
1816 S_B(10)%A_X(7,48)=-0.571428571428571E0_DP
1817 S_B(10)%B_Y(7,48)=-0.428571428571428E0_DP
1818 S_B(10)%B_X(7,49)=1.00000000000000E0_DP
1819 S_B(10)%A_Y(7,49)=-1.00000000000000E0_DP
1820 S_B(10)%B_X(8,3)=-1.00000000000000E0_DP
1821 S_B(10)%A_Y(8,3)=1.00000000000000E0_DP
1822 S_B(10)%A_Y(8,4)=0.500000000000004E0_DP
1823 S_B(10)%A_X(8,5)=8.00000000000000E0_DP
1824 S_B(10)%B_Y(8,5)=7.00000000000000E0_DP
1825 S_B(10)%B_X(8,6)=-1.00000000000000E0_DP
1826 S_B(10)%A_Y(8,6)=1.00000000000000E0_DP
1827 S_B(10)%A_X(8,7)=1.16666666666667E0_DP
1828 S_B(10)%B_Y(8,7)=1.16666666666667E0_DP
1829 S_B(10)%B_X(8,8)=24.5000000000000E0_DP
1830 S_B(10)%A_Y(8,8)=-24.5000000000000E0_DP
1831 S_B(10)%A_X(8,9)=7.00000000000000E0_DP
1832 S_B(10)%B_Y(8,9)=7.00000000000000E0_DP
1833 S_B(10)%B_X(8,11)=1.75000000000001E0_DP
1834 S_B(10)%A_Y(8,11)=-5.25000000000001E0_DP
1835 S_B(10)%A_X(8,12)=-49.0000000000000E0_DP
1836 S_B(10)%B_Y(8,12)=-42.0000000000000E0_DP
1837 S_B(10)%B_X(8,13)=21.0000000000000E0_DP
1838 S_B(10)%A_Y(8,13)=-21.0000000000000E0_DP
1839 S_B(10)%A_X(8,16)=-5.25000000000001E0_DP
1840 S_B(10)%B_Y(8,16)=-5.25000000000001E0_DP
1841 S_B(10)%B_X(8,17)=-52.5000000000000E0_DP
1842 S_B(10)%A_Y(8,17)=52.5000000000000E0_DP
1843 S_B(10)%A_X(8,18)=-35.0000000000000E0_DP
1844 S_B(10)%B_Y(8,18)=-35.0000000000000E0_DP
1845 S_B(10)%B_X(8,22)=-3.50000000000001E0_DP
1846 S_B(10)%A_Y(8,22)=7.00000000000001E0_DP
1847 S_B(10)%A_X(8,23)=42.0000000000000E0_DP
1848 S_B(10)%B_Y(8,23)=35.0000000000000E0_DP
1849 S_B(10)%B_X(8,24)=-35.0000000000000E0_DP
1850 S_B(10)%A_Y(8,24)=35.0000000000000E0_DP
1851 S_B(10)%A_X(8,29)=3.00000000000000E0_DP
1852 S_B(10)%B_Y(8,29)=3.00000000000000E0_DP
1853 S_B(10)%B_X(8,30)=17.5000000000000E0_DP
1854 S_B(10)%A_Y(8,30)=-17.5000000000000E0_DP
1855 S_B(10)%A_X(8,31)=21.0000000000000E0_DP
1856 S_B(10)%B_Y(8,31)=21.0000000000000E0_DP
1857 S_B(10)%B_X(8,37)=0.750000000000001E0_DP
1858 S_B(10)%A_Y(8,37)=-1.25000000000000E0_DP
1859 S_B(10)%A_X(8,38)=-5.00000000000000E0_DP
1860 S_B(10)%B_Y(8,38)=-4.00000000000000E0_DP
1861 S_B(10)%B_X(8,39)=7.00000000000000E0_DP
1862 S_B(10)%A_Y(8,39)=-7.00000000000000E0_DP
1863 S_B(10)%A_X(8,46)=-0.138888888888889E0_DP
1864 S_B(10)%B_Y(8,46)=-0.138888888888889E0_DP
1865 S_B(10)%B_X(8,47)=-0.500000000000000E0_DP
1866 S_B(10)%A_Y(8,47)=0.500000000000000E0_DP
1867 S_B(10)%A_X(8,48)=-1.00000000000000E0_DP
1868 S_B(10)%B_Y(8,48)=-1.00000000000000E0_DP
1869 S_B(10)%B_X(9,1)=-1.00000000000000E0_DP
1870 S_B(10)%A_Y(9,1)=1.00000000000000E0_DP
1871 S_B(10)%A_X(9,2)=9.00000000000000E0_DP
1872 S_B(10)%B_Y(9,2)=8.00000000000000E0_DP
1873 S_B(10)%B_X(9,3)=-1.00000000000000E0_DP
1874 S_B(10)%A_Y(9,3)=1.00000000000000E0_DP
1875 S_B(10)%B_X(9,4)=32.0000000000000E0_DP
1876 S_B(10)%A_Y(9,4)=-32.0000000000000E0_DP
1877 S_B(10)%A_X(9,5)=8.00000000000000E0_DP
1878 S_B(10)%B_Y(9,5)=8.00000000000000E0_DP
1879 S_B(10)%A_X(9,7)=-74.6666666666667E0_DP
1880 S_B(10)%B_Y(9,7)=-65.3333333333333E0_DP
1881 S_B(10)%B_X(9,8)=28.0000000000000E0_DP
1882 S_B(10)%A_Y(9,8)=-28.0000000000000E0_DP
1883 S_B(10)%B_X(9,11)=-98.0000000000000E0_DP
1884 S_B(10)%A_Y(9,11)=98.0000000000000E0_DP
1885 S_B(10)%A_X(9,12)=-56.0000000000000E0_DP
1886 S_B(10)%B_Y(9,12)=-56.0000000000000E0_DP
1887 S_B(10)%A_X(9,16)=98.0000000000000E0_DP
1888 S_B(10)%B_Y(9,16)=84.0000000000000E0_DP
1889 S_B(10)%B_X(9,17)=-70.0000000000000E0_DP
1890 S_B(10)%A_Y(9,17)=70.0000000000000E0_DP
1891 S_B(10)%B_X(9,22)=56.0000000000000E0_DP
1892 S_B(10)%A_Y(9,22)=-56.0000000000000E0_DP
1893 S_B(10)%A_X(9,23)=56.0000000000000E0_DP
1894 S_B(10)%B_Y(9,23)=56.0000000000000E0_DP
1895 S_B(10)%A_X(9,29)=-24.0000000000000E0_DP
1896 S_B(10)%B_Y(9,29)=-20.0000000000000E0_DP
1897 S_B(10)%B_X(9,30)=28.0000000000000E0_DP
1898 S_B(10)%A_Y(9,30)=-28.0000000000000E0_DP
1899 S_B(10)%B_X(9,37)=-5.00000000000000E0_DP
1900 S_B(10)%A_Y(9,37)=5.00000000000000E0_DP
1901 S_B(10)%A_X(9,38)=-8.00000000000000E0_DP
1902 S_B(10)%B_Y(9,38)=-8.00000000000000E0_DP
1903 S_B(10)%A_X(9,46)=0.555555555555555E0_DP
1904 S_B(10)%B_Y(9,46)=0.444444444444444E0_DP
1905 S_B(10)%B_X(9,47)=-1.00000000000000E0_DP
1906 S_B(10)%A_Y(9,47)=1.00000000000000E0_DP
1907 S_B(10)%B_X(10,1)=-1.00000000000000E0_DP
1908 S_B(10)%A_Y(10,1)=1.00000000000000E0_DP
1909 S_B(10)%A_X(10,2)=9.00000000000000E0_DP
1910 S_B(10)%B_Y(10,2)=9.00000000000000E0_DP
1911 S_B(10)%B_X(10,4)=36.0000000000000E0_DP
1912 S_B(10)%A_Y(10,4)=-36.0000000000000E0_DP
1913 S_B(10)%A_X(10,7)=-84.0000000000000E0_DP
1914 S_B(10)%B_Y(10,7)=-84.0000000000000E0_DP
1915 S_B(10)%B_X(10,11)=-126.000000000000E0_DP
1916 S_B(10)%A_Y(10,11)=126.000000000000E0_DP
1917 S_B(10)%A_X(10,16)=126.000000000000E0_DP
1918 S_B(10)%B_Y(10,16)=126.000000000000E0_DP
1919 S_B(10)%B_X(10,22)=84.0000000000000E0_DP
1920 S_B(10)%A_Y(10,22)=-84.0000000000000E0_DP
1921 S_B(10)%A_X(10,29)=-36.0000000000000E0_DP
1922 S_B(10)%B_Y(10,29)=-36.0000000000000E0_DP
1923 S_B(10)%B_X(10,37)=-9.00000000000000E0_DP
1924 S_B(10)%A_Y(10,37)=9.00000000000000E0_DP
1925 S_B(10)%A_X(10,46)=1.00000000000000E0_DP
1926 S_B(10)%B_Y(10,46)=1.00000000000000E0_DP
1927
1928end subroutine set_s_b 
1929 
1930
1931
1932!!!!!!!!!!!!!!!  New routines for solving Maxwell's Equations !!!!!!!!!!!!!!!
1933
1934 SUBROUTINE  get_bend_coeff(s_b0t,NO1)
1935    implicit none
1936    integer no,n,i,k,j(2),NO1,l
1937    type(taylor) x,y,h,df,ker,sol
1938    type(complextaylor) z
1939    type(taylor) f,kick_x,kick_y
1940    type(damap) y0
1941    real(dp) h0,cker,cker0
1942     TYPE(B_CYL), intent(inout) :: s_b0t
1943 
1944    no=NO1
1945
1946
1947    call init(no,1,0,0)
1948
1949
1950    call alloc(x,y,kick_x,kick_y)
1951    call alloc(z)
1952    call alloc(f,h,df,ker,sol)
1953    call alloc(y0)
1954
1955
1956      x=1.d0.mono.1
1957      y=1.d0.mono.2
1958      y0=1
1959      y0%v(2)=0
1960     z=x-i_*y   
1961   
1962!!! 
1963
1964
1965
1966   
1967    h0=1.d0
1968    h=(1.d0+h0*x)
1969    do k=1,no
1970    !  erect multipole
1971    f=dreal(-z**K/K)
1972   
1973    df=f
1974
1975    do i=k,no-1 !k+1
1976     df=h0*(df.d.1)/h
1977     call  invert_laplace(df)
1978      sol=f+df
1979      sol=-(sol.d.1)/h
1980      j=0
1981      j(1)=i
1982      cker=(sol.sub.j)
1983      df=df-cker*dreal(-z**(I+1)/(I+1))
1984      f=f+df
1985     enddo
1986
1987  !     write(16,*) " Erect ",k
1988  !     call clean_taylor(f,f,1.d-6)
1989  !     call print(f,16)
1990  !   y=((f.d.1).d.1)+((f.d.2).d.2)-h0*(f.d.1)/h
1991  !  y=y.cut.(no-1)
1992!       call print(y,6)
1993!       write(6,*) full_abs(y)
1994       kick_x=(f.d.1)
1995       kick_y=(f.d.2)
1996       call clean_taylor(kick_x,kick_x,1.d-6)
1997       call clean_taylor(kick_y,kick_y,1.d-6)
1998       do l=1,s_b0t%n_mono
1999        j(1)=s_b0t%i(l)
2000        j(2)=s_b0t%j(l)
2001        s_b0t%b_x(k,l)=kick_x.sub.j
2002        s_b0t%b_y(k,l)=kick_y.sub.j
2003       enddo
2004
2005   !    bx=bx*y0
2006   !    call clean_taylor(bx,bx,1.d-6)
2007   !    call print(bx,6)
2008enddo   
2009
2010
2011    do k=1,no
2012    !  Skew multipole
2013    f=aimag(-z**K/K)
2014    df=f
2015    do i=k,no-1 !k+1
2016     df=h0*(df.d.1)/h
2017     call  invert_laplace(df)
2018      sol=f+df
2019      sol=(sol.d.2)/h
2020      j=0
2021      j(1)=i
2022      cker=(sol.sub.j)
2023      df=df-cker*aimag(-z**(I+1)/(I+1))
2024      f=f+df
2025     enddo
2026  !     write(16,*) " Skew ",k
2027  !     call clean_taylor(f,f,1.d-6)
2028  !     call print(f,16)
2029  !   y=((f.d.1).d.1)+((f.d.2).d.2)-h0*(f.d.1)/h
2030  !  y=y.cut.(no-1)
2031 !      call print(y,6)
2032  !     write(6,*) full_abs(y)
2033       kick_x=(f.d.1)
2034       kick_y=(f.d.2)
2035       call clean_taylor(kick_x,kick_x,1.d-6)
2036       call clean_taylor(kick_y,kick_y,1.d-6)
2037       do l=1,s_b0t%n_mono
2038        j(1)=s_b0t%i(l)
2039        j(2)=s_b0t%j(l)
2040        s_b0t%a_x(k,l)=kick_x.sub.j
2041        s_b0t%a_y(k,l)=kick_y.sub.j
2042       enddo
2043
2044     !  write(6,*) k
2045     !  write(6,*) S_B0%A_X(k,7)
2046     !  write(6,*) S_B0%A_Y(k,8)    !   bx=bx*y0
2047    !   call clean_taylor(bx,bx,1.d-6)
2048    !   call print(bx,6)
2049enddo   
2050
2051
2052    call kill(x,y,kick_x,kick_y)
2053    call kill(z)
2054    call kill(f,h,df,ker,sol)
2055    call kill(y0)
2056
2057    end subroutine get_bend_coeff
2058
2059
2060    subroutine invert_laplace(df)
2061    implicit none
2062    type(taylor), intent(inout) :: df
2063    type(taylorresonance) dfr
2064     call alloc(dfr)
2065
2066     dfr=df
2067       dfr%cos=((dfr%cos.i.1).i.2)/4.0_dp
2068       dfr%sin=((dfr%sin.i.1).i.2)/4.0_dp
2069      df=dfr
2070
2071     call kill(dfr)
2072    end subroutine invert_laplace
2073
2074  subroutine make_coef(b,no)
2075    implicit none
2076    integer no
2077    integer i,j,k,a,m , ic
2078
2079    type(B_CYL) b
2080
2081    ic=1
2082    b%firsttime=-100
2083    allocate(b%nmul)
2084    allocate(b%n_mono)
2085    b%nmul=no
2086    b%n_mono=((no+2-ic)*(no+1-ic))/2
2087    allocate(b%i(b%n_mono),b%j(b%n_mono))
2088    allocate(b%a_x(no,b%n_mono),b%a_y(no,b%n_mono))
2089    allocate(b%b_x(no,b%n_mono),b%b_y(no,b%n_mono))
2090
2091    do i=1,no
2092       do j=1,b%n_mono
2093          b%a_x(i,j)=0.0_dp
2094          b%a_y(i,j)=0.0_dp
2095          b%b_x(i,j)=0.0_dp
2096          b%b_y(i,j)=0.0_dp
2097       enddo
2098    enddo
2099
2100
2101    k=0
2102    m=no-1
2103    do a=m,1,-1
2104       do j=m-a,1,-1
2105          k=k+1
2106          b%i(k)=a
2107          b%j(k)=j
2108       enddo
2109       k=k+1
2110       b%i(k)=a
2111       b%j(k)=0
2112    enddo
2113    do j=m,1,-1
2114       k=k+1
2115       b%i(k)=0
2116       b%j(k)=j
2117    enddo
2118    k=k+1
2119    b%i(k)=0
2120    b%j(k)=0
2121
2122  end subroutine make_coef
2123
2124  subroutine nul_coef(b)
2125    implicit none
2126    type(B_CYL) b
2127    if(b%firsttime/=-100) then
2128       nullify(b%nmul)
2129       nullify(b%n_mono)
2130       nullify(b%i)
2131       nullify(b%j)
2132       nullify(b%a_x)
2133       nullify(b%a_y)
2134       nullify(b%b_x)
2135       nullify(b%b_y)
2136       b%firsttime=-100
2137    else
2138       deallocate(b%nmul)
2139       deallocate(b%n_mono)
2140       deallocate(b%i)
2141       deallocate(b%j)
2142       deallocate(b%a_x)
2143       deallocate(b%a_y)
2144       deallocate(b%b_x)
2145       deallocate(b%b_y)
2146
2147    endif
2148
2149
2150  end subroutine nul_coef
2151
2152
2153end module S_status
Note: See TracBrowser for help on using the repository browser.