source: PSPA/madxPSPA/libs/ptc/src/h_definition.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: 16.5 KB
Line 
1!The Full Polymorphic Package
2!Copyright (C) Etienne Forest
3
4module definition
5  !  use define_newda
6  !  use precision_constants   ! added to replace use define_newda
7  !  use scratch_size
8  !  use DABNEW
9  !  use my_own_1D_TPSA
10  use lielib_yang_berz, junk_no=>no,junk_nd=>nd,junk_nd2=>nd2,junk_ndpt=>ndpt,junk_nv=>nv
11  !  use newda
12  !  USE LIELIB_ETIENNE
13  implicit none
14  public
15  logical(lp) :: newread=.false. ,newprint =  .false. , first_time = .true.
16  logical(lp) :: print77=.true. ,read77 =  .true.
17  logical(lp) :: no_ndum_check = .false.
18  logical(lp),TARGET :: insane_PTC = .false.
19  logical(lp),TARGET :: setknob = .false. !@1 Real part of knobs cannot set
20  logical(lp),TARGET :: knob=.true. !@1 Knobs are effective
21  integer, target :: npara_fpp !@1 position of last non-parameter tpsa variable
22  !  complex(dp), parameter :: i_ = = cmplx(zero,one,kind=dp)
23  !  complex(dp)  :: i_ = (zero,one)    ! cmplx(zero,one,kind=dp)
24  complex(dp), parameter :: i_ = ( 0.0_dp,1.0_dp )    ! cmplx(zero,one,kind=dp)
25  integer master
26  !  integer,parameter::lnv=100
27  !  scratch variables
28  INTEGER iassdoluser(ndumt)
29  integer DUMMY,temp
30  integer iass0user(ndumt)
31  integer,parameter::ndim2=2*ndim
32  integer,parameter::mmmmmm1=1,mmmmmm2=2,mmmmmm3=3,mmmmmm4=4
33  !  type (taylorlow) DUMMYl,templ             !,DUMl(ndum)
34  !  private NDC,NDC2,NDT,IREF,itu,iflow,jtune,nres,ifilt  ,idpr
35  !  private nplane,idsta,ista
36  !  private xintex,dsta,sta,angle,rad,ps,rads,mx
37  ! numerical differentiation by knobs
38  logical(lp) :: knob_numerical=.false.
39  real(dp) ::  knob_eps(lnv)=1e-6_dp
40  integer ::  knob_i =0
41  INTEGER,PARAMETER::NMAX=20
42  integer,private,parameter::n_max=10   ! sagan stuff
43  INTEGER, PARAMETER :: CASE1=1,CASE2=2, CASE0=0, CASEP1=-1,CASEP2=-2
44  INTEGER, PARAMETER :: CASET=3,CASETF1=4,CASETF2=5
45  INTEGER,PARAMETER  :: ISPIN0R=1,ISPIN1R=3
46  logical(lp) :: doing_ac_modulation_in_ptc=.false.
47  integer, target :: nb_ =0   ! global group index
48  !
49  TYPE sub_taylor
50     INTEGER j(lnv)
51     INTEGER min,max
52  END TYPE sub_taylor
53
54  !!&1
55  TYPE taylor
56     INTEGER I !@1  integer I is a pointer in old da-package of Berz
57     !     type (taylorlow) j !@1   Taylorlow is an experimental type not supported
58  END TYPE taylor
59  !@2  UNIVERSAL_TAYLOR is used by Sagan in BMAD Code at Cornell
60  !@2  Also used by MAD-XP
61  TYPE UNIVERSAL_TAYLOR
62     INTEGER, POINTER:: N,NV    !  Number of coeeficients and number of variables
63     REAL(DP), POINTER,dimension(:)::C  ! Coefficients C(N)
64     INTEGER, POINTER,dimension(:,:)::J ! Exponents of each coefficients J(N,NV)
65  END TYPE UNIVERSAL_TAYLOR
66  !@3 ---------------------------------------------</br>
67  TYPE complextaylor
68     type (taylor) r  !@1 Real part
69     type (taylor) i  !@1 Imaginary part
70  END TYPE complextaylor
71
72  !&2
73  ! this is a real polymorphic type
74
75  TYPE REAL_8
76     TYPE (TAYLOR) T      !@1  USED IF TAYLOR
77     REAL(DP) R           !@1    USED IF REAL
78     !&2
79     INTEGER KIND  !@1  0,1,2,3 (1=REAL,2=TAYLOR,3=TAYLOR KNOB, 0=SPECIAL)
80     INTEGER I   !@1   USED FOR KNOBS AND SPECIAL KIND=0
81     REAL(DP) S   !@1   SCALING FOR KNOBS AND SPECIAL KIND=0
82     LOGICAL(lp) :: ALLOC  !@1 IF TAYLOR IS ALLOCATED IN DA-PACKAGE
83     integer g,nb  !  group index, number in group
84     !&2
85  END TYPE REAL_8
86  !@3 ---------------------------------------------</br>
87  TYPE double_complex
88     type (complextaylor) t
89     complex(dp) r
90     logical(lp) alloc
91     integer kind
92     integer i,j
93     complex(dp) s
94     integer g,nb  !  group index
95  END TYPE double_complex
96
97  type(taylor)        varf1,varf2
98  type(complextaylor) varc1,varc2
99
100  !Radiation
101  TYPE ENV_8
102     type (REAL_8) v
103     type (REAL_8) e(ndim2)
104     type (REAL_8) sigma0(ndim2)
105     type (REAL_8) sigmaf(ndim2)
106  END TYPE ENV_8
107
108  type spinor
109     real(dp) x(3)
110     !  real(dp) G
111  end type spinor
112
113  type spinor_8
114     type(real_8) x(3)
115  end type spinor_8
116
117  type res_spinor_8
118     type(double_complex) x(3)
119  end type res_spinor_8
120
121  !    scratch levels of DA using linked list
122
123  type dascratch
124     type(taylor), pointer :: t
125     TYPE (dascratch),POINTER :: PREVIOUS
126     TYPE (dascratch),POINTER :: NEXT
127  end type dascratch
128
129  TYPE dalevel
130     INTEGER,  POINTER :: N     ! TOTAL ELEMENT IN THE CHAIN
131     !
132     logical(lp),POINTER ::CLOSED
133     TYPE (dascratch), POINTER :: PRESENT
134     TYPE (dascratch), POINTER :: END
135     TYPE (dascratch), POINTER :: START
136     TYPE (dascratch), POINTER :: START_GROUND ! STORE THE GROUNDED VALUE OF START DURING CIRCULAR SCANNING
137     TYPE (dascratch), POINTER :: END_GROUND ! STORE THE GROUNDED VALUE OF END DURING CIRCULAR SCANNING
138  END TYPE dalevel
139
140  !&1
141  TYPE DAMAP
142     TYPE (TAYLOR) V(ndim2)    ! Ndim2=6 but allocated to nd2=2,4,6 ! etienne_oct_2004
143  END TYPE DAMAP
144  !&1
145  !@3 ---------------------------------------------</br>
146
147  TYPE GMAP
148     TYPE (TAYLOR) V(lnv)    !
149     integer N
150  END TYPE GMAP
151
152  !&4
153  TYPE vecfield
154     type (taylor) v(ndim2)          !@1 <font face="Times New Roman">V<sub>i</sub>&#8706;<sub>i</sub></font> Operator
155     integer ifac                    !@1 Type of Factorization 0,1,-1 (One exponent, Dragt-Finn, Reversed Dragt-Finn)
156  END TYPE vecfield
157  !@3 ---------------------------------------------</br>
158  TYPE pbfield
159     type (taylor) h
160     integer ifac
161     integer nd_used
162  END TYPE pbfield
163  !&4
164
165
166  TYPE tree
167     type (taylor) branch(ndim2)
168  END TYPE tree
169
170  !Radiation
171  TYPE radtaylor
172     type (taylor) v
173     type (taylor) e(ndim2)
174  END TYPE radtaylor
175
176  !&5
177  TYPE DRAGTFINN
178     real(dp)  constant(ndim2)
179     type (damap) Linear
180     type (vecfield) nonlinear
181     type (pbfield)  pb
182  END TYPE DRAGTFINN
183  !@3 ---------------------------------------------</br>
184  TYPE reversedragtfinn
185     real(dp)  CONSTANT(NDIM2)
186     type (damap) Linear
187     type (vecfield) nonlinear
188     type (pbfield)  pb
189  END TYPE reversedragtfinn
190  !@3 ---------------------------------------------</br>
191  TYPE ONELIEEXPONENT
192     real(dp) EPS
193     type (vecfield) VECTOR
194     type (pbfield)  pb
195  END TYPE ONELIEEXPONENT
196  !@3 ---------------------------------------------</br>
197
198  TYPE normalform
199     type (damap) A_t   ! Total A  :  A_t= A1 o A_rest
200     type (damap) A1    !@1 Dispersion
201     type (reversedragtfinn) A  !@1 Linear A_t and nonlinear A_t
202     type (dragtfinn) NORMAL    !@1 Normal is the Normal Form R
203     type (damap) DHDJ  !@1 Contains the tunes in convenient form: extracted from NORMAL (=R)
204     real(dp) TUNE(NDIM),DAMPING(NDIM)  !@1 linear tune and linear damping
205     integer nord,jtune                 !@1 nord=1 A1 first order in parameters
206     integer NRES,M(NDIM,NRESO),PLANE(NDIM) !@1 NRES,M(NDIM,NRESO) -> resonances left in the map
207     logical(lp) AUTO
208  END TYPE normalform
209  !@3 ---------------------------------------------</br>
210  TYPE genfield
211     type (taylor) h
212     type (damap) m
213     type (taylor) d(ndim,ndim)
214     type (damap) linear
215     type (damap) lineart
216     type (damap) mt
217     real(dp) constant(ndim2),eps
218     integer imax     !@1 imax=Maximum Number of Iteration (default=1000)
219     integer ifac     !@1 ifac = the map is raised to the power 1/ifac and iterated ifac times (default=1)
220     logical(lp) linear_in     !@1 Linear part is left in the map  (default=.false.)
221     integer no_cut   !@1 Original map is not symplectic on and above no_cut
222  END TYPE genfield
223
224  !@3 ---------------------------------------------</br>
225
226
227  TYPE pbresonance
228     type (pbfield)  cos,sin
229     integer ifac
230  END TYPE pbresonance
231  !@3 ---------------------------------------------</br>
232  TYPE vecresonance
233     type (vecfield)  cos,sin
234     integer ifac
235  END TYPE vecresonance
236  !@3 ---------------------------------------------</br>
237  TYPE taylorresonance
238     type (taylor)  cos,sin
239  END TYPE taylorresonance
240
241  TYPE beamenvelope   !@2 A kind of Normal Form for Radiative Envelope
242     ! radiation normalization
243     type (damap) transpose    ! Transpose of map which acts on polynomials
244     type (taylor) bij         !  Represents the stochastic kick at the end of the turn  Env_f=M Env_f M^t + B
245     TYPE (pbresonance) bijnr   !  Equilibrium beam sizes in resonance basis
246     real(dp) s_ij0(6,6)  !  equilibrium beam sizes
247     type (taylor) sij0  !  equilibrium beam sizes
248     real(dp) emittance(3),tune(3),damping(3)
249     logical(lp) AUTO,STOCHASTIC
250     real(dp)  KICK(3)   ! fake kicks for tracking stochastically
251     type (damap) STOCH  ! Diagonalized of stochastic part of map for tracking stochastically
252  END TYPE beamenvelope
253
254
255  type  tree_element   !@1  USED FOR FAST TRACKING IN O_TREE_ELEMENT.F90
256     real(dp) ,  DIMENSION(:), POINTER :: CC
257     real(dp) ,  DIMENSION(:), POINTER :: fix
258     integer,  DIMENSION(:), POINTER :: JL,JV
259     INTEGER,POINTER :: N,ND2,no
260  end  type tree_element
261
262  type spinmatrix
263     type(real_8) s(3,3)
264  end type spinmatrix
265
266  type damapspin
267     type(damap) M
268     type(spinmatrix) s
269     !     type(real_8) s(3,3)
270     real(dp) e_ij(6,6)
271  end type damapspin
272
273  type normal_spin
274     type(normalform) N   ! regular orbital normal form
275     type(damapspin) a1   ! brings to fixed point
276     type(damapspin) ar   ! normalises around the fixed point
277     type(damapspin) as   ! pure spin map
278     type(damapspin) a_t  ! !! (a_t%m,a_t%s) = (a1%m, I ) o (I ,as%s) o (ar%m,I)
279!!!  extra spin info
280     integer M(NDIM,NRESO),MS(NRESO),NRES  ! orbital and spin resonances to be left in the map
281     type(real_8) n0(3)     ! n0 vector
282     type(real_8) theta0    !  angle for the matrix around the orbit (analogous to linear tunes)
283     real(dp) nu    !  spin tune
284!!!Envelope radiation stuff
285     real(dp) s_ij0(6,6)  !  equilibrium beam sizes
286     complex(dp) s_ijr(6,6)  !  equilibrium beam sizes in resonance basis
287     real(dp) emittance(3),tune(3),damping(3)   ! equilibrium emittances (partially well defined only for infinitesimal damping)
288     logical(lp) AUTO,STOCHASTIC
289     real(dp)  KICK(3)   ! fake kicks for tracking stochastically
290     real(dp)  STOCH(6,6)  ! Diagonalized of stochastic part of map for tracking stochastically
291     real(dp)  STOCH_inv(6,6)  ! Diagonalized of stochastic part of map for tracking stochastically
292  end type normal_spin
293
294  include "a_def_frame_patch_chart.inc"
295  include "a_def_all_kind.inc"
296  include "a_def_sagan.inc"
297  include "a_def_element_fibre_layout.inc"
298
299  type(fibre), pointer :: lost_fibre
300  type(integration_node), pointer :: lost_node
301
302  type rf_phasor
303     real(dp) x(2)
304     real(dp) om
305     real(dp) t
306  end type rf_phasor
307
308  type rf_phasor_8
309     type(real_8)  x(2)
310     type(real_8) om
311     reaL(DP) t
312  end type rf_phasor_8
313
314  type probe
315     real(dp) x(6)
316     type(spinor) s(ISPIN0R:ISPIN1R)
317     type(rf_phasor) AC
318     logical u
319     type(integration_node),pointer :: lost_node
320  end type probe
321
322  type probe_8
323     type(real_8) x(6)     ! Polymorphic orbital ray
324     type(spinor_8) s(ISPIN0R:ISPIN1R)   ! Polymorphic spin
325     type(rf_phasor_8) AC  ! Modulation
326     real(dp) E_ij(6,6)   !  Envelope
327     !   stuff for exception
328     logical u
329     type(integration_node),pointer :: lost_node
330  end type probe_8
331
332  type TEMPORAL_PROBE
333     TYPE(probe)  XS
334     TYPE(INTEGRATION_NODE), POINTER :: NODE
335     real(DP)  DS,POS(6)
336  END type TEMPORAL_PROBE
337
338  type TEMPORAL_BEAM
339     TYPE(TEMPORAL_PROBE), pointer :: TP(:)
340     real(DP) a(3),ent(3,3),p0c,total_time
341     integer n
342     type(integration_node),pointer :: c   ! pointer close to a(3)
343     type(internal_state)  state
344  END type TEMPORAL_BEAM
345
346contains
347
348  SUBROUTINE RESET_APERTURE_FLAG(complete)
349    IMPLICIT NONE
350    logical(lp), optional :: complete
351    logical(lp)  comp
352    !    IF(c_%WATCH_USER) THEN
353    !       IF(.NOT.ALLOW_TRACKING) THEN
354    !          WRITE(6,*) "  EXECUTION OF THE CODE MUST BE INTERRUPTED AT YOUR REQUEST"
355    !          WRITE(6,*) "  YOU DID NOT CHECK THE APERTURE STATUS"
356    !          WRITE(6,*) "  USING A CALL TO PRODUCE_APERTURE_FLAG"
357    !          WRITE(6,*) "  BEFORE CALLING A TRACKING FUNCTION"
358    !          STOP 666
359    !       ENDIF
360    !       IF(.NOT.c_%check_stable) THEN
361    !          WRITE(6,*) "  EXECUTION OF THE CODE MUST BE INTERRUPTED AT YOUR REQUEST"
362    !          WRITE(6,*) " CODE MOST LIKELY DIED IN PURE DA/TPSA/LIE OPERATIONS  "
363    !          STOP 667
364    !       ENDIF
365    !    ENDIF
366    comp=.true.
367    if(present(complete)) comp=complete
368    c_%STABLE_DA =.TRUE.
369    c_%CHECK_STABLE =.TRUE.
370    c_%CHECK_MADX_APERTURE =.TRUE.
371    c_%stable_da =.true.
372    if(comp) then
373       xlost=0.0_dp
374       messagelost=" Aperture has been reset "
375       nullify(lost_fibre)
376       nullify(lost_node)
377
378    endif
379  END   SUBROUTINE RESET_APERTURE_FLAG
380
381  SUBROUTINE PRODUCE_APERTURE_FLAG(I)
382    IMPLICIT NONE
383    INTEGER I
384    I=0
385    IF(.NOT.c_%CHECK_STABLE) THEN
386       I=1
387    ENDIF
388
389
390    !    ALLOW_TRACKING=.TRUE.
391
392  END   SUBROUTINE PRODUCE_APERTURE_FLAG
393
394  ! moved here from sa_extend_poly.f90
395
396  REAL(DP) FUNCTION  ROOT(X)  ! REPLACES SQRT(X)
397    IMPLICIT NONE
398    REAL(DP),INTENT(IN)::X
399    IF(.NOT.c_%CHECK_STABLE) then
400       ROOT=1.0_dp
401       return
402    endif
403
404    IF((X<0.0_dp).AND.c_%ROOT_CHECK) THEN
405       ROOT=1.0_dp
406       c_%CHECK_STABLE=.FALSE.
407       messagelost="Root undefined "
408    ELSEIF(X>=0.0_dp) THEN
409       ROOT=SQRT(X)
410    ELSE      !  IF X IS NOT A NUMBER
411       ROOT=1.0_dp
412       c_%CHECK_STABLE=.FALSE.
413    ENDIF
414
415  END FUNCTION ROOT
416
417  REAL(DP) FUNCTION  ARCSIN(X)  ! REPLACES ASIN(X)
418    IMPLICIT NONE
419    REAL(DP),INTENT(IN)::X
420    IF(.NOT.c_%CHECK_STABLE) then
421       ARCSIN=0.0_dp
422       return
423    endif
424    IF((ABS(X)>1.0_dp).AND.c_%ROOT_CHECK) THEN
425       ARCSIN=0.0_dp
426       c_%CHECK_STABLE=.FALSE.
427       messagelost="Arcsin undefined "
428    ELSEIF(ABS(X)<=1.0_dp) THEN
429       ARCSIN=ASIN(X)
430    ELSE      !  IF X IS NOT A NUMBER
431       ARCSIN=0.0_dp
432       c_%CHECK_STABLE=.FALSE.
433    ENDIF
434
435  END FUNCTION ARCSIN
436
437  REAL(DP) FUNCTION  ARCCOS(X)  ! REPLACES ACOS(X)
438    IMPLICIT NONE
439    REAL(DP),INTENT(IN)::X
440    IF(.NOT.c_%CHECK_STABLE) then
441       ARCCOS=0.0_dp
442       return
443    endif
444    IF((ABS(X)>1.0_dp).AND.c_%ROOT_CHECK) THEN
445       ARCCOS=0.0_dp
446       c_%CHECK_STABLE=.FALSE.
447       messagelost="Arccos undefined "
448    ELSEIF(ABS(X)<=1.0_dp) THEN
449       ARCCOS=ACOS(X)
450    ELSE      !  IF X IS NOT A NUMBER
451       ARCCOS=0.0_dp
452       c_%CHECK_STABLE=.FALSE.
453    ENDIF
454
455  END FUNCTION ARCCOS
456
457  REAL(DP) FUNCTION  LOGE(X)  ! REPLACES ACOS(X)
458    IMPLICIT NONE
459    REAL(DP),INTENT(IN)::X
460    IF(.NOT.c_%CHECK_STABLE) then
461       LOGE=0.0_dp
462       return
463    endif
464
465    IF(X<=0.0_dp.AND.c_%ROOT_CHECK) THEN
466       LOGE=0.0_dp
467       c_%CHECK_STABLE=.FALSE.
468       messagelost="Log undefined "
469    ELSE
470       LOGE=LOG(X)
471    ENDIF
472
473  END FUNCTION LOGE
474
475
476
477  REAL(DP) FUNCTION  COSEH(X) ! REPLACES COSH(X)
478    IMPLICIT NONE
479    REAL(DP),INTENT(IN)::X
480    IF(.NOT.c_%CHECK_STABLE) then
481       COSEH=1.0_dp
482       return
483    endif
484
485    IF((ABS(X)>c_%hyperbolic_aperture).AND.c_%ROOT_CHECK) THEN
486       COSEH=1.0_dp
487       c_%CHECK_STABLE=.FALSE.
488       messagelost="Coseh undefined "
489    ELSEIF(ABS(X)<=c_%hyperbolic_aperture) THEN
490       COSEH=COSH(X)
491    ELSE      !  IF X IS NOT A NUMBER
492       COSEH=1.0_dp
493       c_%CHECK_STABLE=.FALSE.
494    ENDIF
495
496  END FUNCTION COSEH
497
498  REAL(DP) FUNCTION  SINEH(X) ! REPLACES SINH(X)
499    IMPLICIT NONE
500    REAL(DP),INTENT(IN)::X
501    IF(.NOT.c_%CHECK_STABLE) then
502       SINEH=0.0_dp
503       return
504    endif
505
506    IF((ABS(X)>c_%hyperbolic_aperture).AND.c_%ROOT_CHECK) THEN
507       SINEH=0.0_dp
508       c_%CHECK_STABLE=.FALSE.
509       messagelost="Sineh undefined "
510    ELSEIF(ABS(X)<=c_%hyperbolic_aperture) THEN
511       SINEH=SINH(X)
512    ELSE      !  IF X IS NOT A NUMBER
513       SINEH=0.0_dp
514       c_%CHECK_STABLE=.FALSE.
515    ENDIF
516
517  END FUNCTION SINEH
518
519  REAL(DP) FUNCTION  arctan(X) ! REPLACES SINH(X)
520    IMPLICIT NONE
521    REAL(DP),INTENT(IN)::X
522    IF(.NOT.c_%CHECK_STABLE) then
523       arctan=0.0_dp
524       return
525    endif
526
527    IF((ABS(X)>c_%hyperbolic_aperture).AND.c_%ROOT_CHECK) THEN
528       arctan=0.0_dp
529       c_%CHECK_STABLE=.FALSE.
530       messagelost="Arctan undefined "
531    ELSEIF(ABS(X)<=c_%hyperbolic_aperture) THEN
532       arctan=atan(X)
533    ELSE      !  IF X IS NOT A NUMBER
534       arctan=0.0_dp
535       c_%CHECK_STABLE=.FALSE.
536    ENDIF
537
538  END FUNCTION arctan
539
540
541
542end module definition
Note: See TracBrowser for help on using the repository browser.