Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

o_tree_element.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: 142.7 KB

Line
1	!The Full Polymorphic Package
2	!Copyright (C) Etienne Forest
3
4	module tree_element_MODULE
5	USE polymorphic_complextaylor, mon1=>mon
6	IMPLICIT NONE
7	public
8	integer,private,parameter::ndd=6
9
10	PRIVATE track_TREE,track_TREEP,KILL_TREE,KILL_TREE_N,SET_TREE
11	PRIVATE track_TREE_G,track_TREEP_g
12	PRIVATE ALLOC_SPINOR_8,ALLOC_probe_8
13	PRIVATE KILL_SPINOR_8,KILL_probe_8,KILL_DASPIN
14	PRIVATE EQUAL_SPINOR8_SPINOR8,EQUAL_IDENTITY_SPINOR_8 !,EQUAL_SPINOR8_RAY8,EQUAL_RAY8_SPINOR8,
15	PRIVATE EQUAL_IDENTITY_probe_8,ALLOC_daspin,EQUAL_DAMAPSPIN_RAY8 ,EQUAL_RAY8_DAMAPSPIN
16	private alloc_normal_spin,kill_normal_spin,EQUAL_IDENTITY_probe
17	private READ_DASPIN,PRINT_DASPIN,EQUAL_IDENTITY_SPINOR,EQUAL_PROBE_REAL6
18	PRIVATE EQUAL_SPINOR8_SPINOR,EQUAL_PROBE8_PROBE,EQUAL_PROBE8_REAL6
19	private EQUAL_IDENTITY_SPINOR_8_r3 ,EQUAL_SPINOR_SPINOR8,exp_spinor_8
20	private INTO_RES_SPIN8_eq,INTO_SPIN8_from_RES_eq,ALLOC_rf_phasor_8,KILL_rf_phasor_8
21
22	private find_ar
23	private find_ap,PRINT_spinor_8,dot_spinor_8,dot_spinor,print_res_spinor_8
24	private find_n_thetar,find_n_thetap,read_spinor_8
25	! private smatp,smatmulp
26	private inv_asr,inv_asp !,inv_as
27	PRIVATE EQUAL_DAmapSPIN_int,daddsc,scdadd,EQUAL_PROBE8_PROBE8,PRINT_probe8
28	PRIVATE concat,assmap,EQUAL_damapspin,CUTORDER,assprobe_8,POWMAP,cmul,addm,mmul,spin8_mul_map
29	private read_probe8,ALLOC_33t,ALLOC_33p,KILL_33t,KILL_33p,zero_33t
30	private purge_transverse
31	private get_spin_nx_r,get_spin_nx_t,get_spin_nx_rd,get_spin_nx_probe,get_spin_nx_spinor_8
32	private scdaddo,daddsco,damapspin_spinor8_mul,damapspin_spinor_mul,eval_spinor_8
33	private real_8REAL6,REAL6real_8,real_8REAL_8,PRINT6
34	private check_fix,test_jc,A_OPT_damap,K_OPT_damap,factor_am,factor_as,concatxp
35	private find_axisp,spin8_scal8_map,add_spin8_spin8,sub_spin8_spin8,mul_spin8_spin8
36	private find_perp_basisp,find_exponentp
37	private alloc_res_SPINOR_8,KILL_res_SPINOR_8,full_abst,EQUAL_RF8_RF8,extract_envelope_probe8
38	PRIVATE EQUAL_RF8_RF,EQUAL_RF_RF8,print_rf_phasor_8,concat_envelope,extract_envelope_damap
39	private EQUAL_DAMAP_RAY8,spimat_spinmat,EQUAL_damapspin_smat,EQUAL_smat_damapspin
40	private flip,dlie,dphase,phase_shift ! flip in lielib
41	integer, target :: spin_extra_tpsa = 0 ,n0_normal= 2
42	integer :: clockwise=1
43	logical(lp) :: force_positive=.false.
44	logical(lp) :: use_ptc_ac_position=.false.
45	integer, private :: nd_used,i_phase,i_plane, decal
46	logical :: onelie = my_false
47	logical :: firstfac=.true.
48	integer, private, parameter :: nfac=20
49	real(dp), private :: fac(0:nfac)
50
51	INTERFACE assignment (=)
52	!
53	MODULE PROCEDURE REAL_8REAL6
54	MODULE PROCEDURE REAL6REAL_8
55	MODULE PROCEDURE real_8REAL_8
56	module PROCEDURE spimat_spinmat
57	!
58	MODULE PROCEDURE EQUAL_IDENTITY_probe ! probe_8=0 p%s=Id p%x=x_i (tpsa) i=1,npara_fpp
59	MODULE PROCEDURE EQUAL_IDENTITY_probe_8 ! probe=0 p%s=Id p%x=0.d0
60	MODULE PROCEDURE EQUAL_IDENTITY_SPINOR_8 ! SPINOR_8%x(r)=1 if SPINOR_8=r otherwise zero
61	MODULE PROCEDURE EQUAL_IDENTITY_SPINOR ! spinor%x(r)=1 if spinor=r otherwise zero
62	MODULE PROCEDURE EQUAL_PROBE_REAL6
63	MODULE PROCEDURE EQUAL_PROBE8_REAL6
64	MODULE PROCEDURE EQUAL_PROBE8_PROBE
65
66	MODULE PROCEDURE EQUAL_SPINOR8_SPINOR8
67	MODULE PROCEDURE EQUAL_SPINOR8_SPINOR
68	MODULE PROCEDURE EQUAL_DAMAPSPIN_RAY8
69	MODULE PROCEDURE EQUAL_RAY8_DAMAPSPIN
70	MODULE PROCEDURE EQUAL_DAMAP_RAY8
71
72	MODULE PROCEDURE EQUAL_DAmapSPIN_int ! damapspin = Identity (or zero)
73	MODULE PROCEDURE EQUAL_PROBE8_PROBE8
74	MODULE PROCEDURE EQUAL_damapspin ! damapspin <= real spin matrix
75	MODULE PROCEDURE EQUAL_damapspin_smat ! damapspin => real spin matrix
76	MODULE PROCEDURE EQUAL_smat_damapspin
77	MODULE PROCEDURE EQUAL_SPINOR_SPINOR8
78	MODULE PROCEDURE EQUAL_PROBE_PROBE8
79	MODULE PROCEDURE normalise_spin
80	MODULE PROCEDURE INTO_RES_SPIN8_eq !
81	MODULE PROCEDURE INTO_SPIN8_from_RES_eq
82	MODULE PROCEDURE EQUAL_RF8_RF8
83	MODULE PROCEDURE EQUAL_RF8_RF
84	MODULE PROCEDURE EQUAL_RF_RF8
85	! MODULE PROCEDURE EQUAL_nn_damap
86	end INTERFACE
87
88	INTERFACE OPERATOR (*)
89	MODULE PROCEDURE concat ! damapspin= damapspin o damapspin
90	MODULE PROCEDURE cmul ! damapspin%s= real(dp) * damapspin%s orbital unchanged
91	MODULE PROCEDURE mmul ! damapspin%s=damapspin%s o damap
92	MODULE PROCEDURE damapspin_spinor_mul ! spinor_8%s(i)= damapspin%s(i,j) * spinor%s(j)
93	MODULE PROCEDURE damapspin_spinor8_mul ! spinor_8%s(i)= damapspin%s(i,j) * spinor_8%s(j)
94	MODULE PROCEDURE spin8_mul_map ! spinor_8=spinor_8 o damap
95	MODULE PROCEDURE eval_spinor_8 ! spinor=spinor_8 * xp(lnv)
96	MODULE PROCEDURE concatxp ! damapspin=damapspin*xp(lnv)
97	MODULE PROCEDURE spin8_scal8_map ! real_8*spinor_8
98	MODULE PROCEDURE mul_spin8_spin8 ! mul_spin8_spin8= spin8 X spin8 cross product
99	END INTERFACE
100
101	INTERFACE OPERATOR (**)
102	MODULE PROCEDURE POWMAP
103	END INTERFACE
104
105	INTERFACE OPERATOR (.cut.)
106	MODULE PROCEDURE CUTORDER
107	END INTERFACE
108
109	INTERFACE OPERATOR (.dot.)
110	MODULE PROCEDURE dot_spinor
111	MODULE PROCEDURE dot_spinor_8
112	END INTERFACE
113
114	INTERFACE operator (+)
115	MODULE PROCEDURE scdadd
116	MODULE PROCEDURE daddsc
117	MODULE PROCEDURE scdaddo
118	MODULE PROCEDURE daddsco
119	MODULE PROCEDURE addm ! (damapspin1%m, damapspin1%s + damapspin2%s)
120	MODULE PROCEDURE add_spin8_spin8 ! spinor_8 = spinor_8 + spinor_8
121	END INTERFACE
122
123	INTERFACE operator (-)
124	MODULE PROCEDURE sub_spin8_spin8
125	END INTERFACE
126
127	INTERFACE extract_beam_sizes !
128	MODULE PROCEDURE extract_envelope_damap !
129	MODULE PROCEDURE extract_envelope_probe8 !
130	END INTERFACE
131
132	INTERFACE exp ! damapspin = exp(spinor_8)
133	MODULE PROCEDURE exp_spinor_8 !
134	END INTERFACE
135
136	INTERFACE texp ! damapspin = exp(spinor_8)
137	MODULE PROCEDURE exp_spinor_8 !
138	END INTERFACE
139
140	INTERFACE find_n0 ! (s0(3,3),n0(3))
141	MODULE PROCEDURE find_n_thetar ! finds n0 the naive way
142	MODULE PROCEDURE find_n_thetap !(s0(3,3),n0(3), spinor_8 optional)
143	END INTERFACE
144
145	INTERFACE find_axis ! (ds,spinor_8)
146	MODULE PROCEDURE find_axisp ! s= exp(spinor_8.dot.L)
147	END INTERFACE
148
149	INTERFACE find_perp_basis !
150	MODULE PROCEDURE find_perp_basisp !
151	END INTERFACE
152
153	INTERFACE find_exponent !
154	MODULE PROCEDURE find_exponentp !
155	END INTERFACE
156
157	INTERFACE find_a ! (n(3),a(3,3)) if you have n you get a(3,3)
158	MODULE PROCEDURE find_ar ! such that ai.s.a = exp(theta L_y)
159	MODULE PROCEDURE find_ap
160	END INTERFACE
161
162	INTERFACE inv_as ! gets inverse of spin matrix by transposing
163	MODULE PROCEDURE inv_asr
164	MODULE PROCEDURE inv_asp
165	END INTERFACE
166
167	INTERFACE get_spin_n0 ! (S,theta0,n0)
168	MODULE PROCEDURE get_spin_nx_r ! real s(3,3) = exp(theta0 n0.L)
169	MODULE PROCEDURE get_spin_nx_rd ! takes DS (damapspin but returns real)
170	MODULE PROCEDURE get_spin_nx_t ! takes DS (damapspin but returns real_8)
171	MODULE PROCEDURE get_spin_nx_probe ! takes probe (damapspin but returns real)
172	MODULE PROCEDURE get_spin_nx_spinor_8
173	END INTERFACE
174
175	INTERFACE find_exponent_jet !
176	MODULE PROCEDURE find_exponent_jet_p
177	end interface
178
179	INTERFACE full_abs !
180	MODULE PROCEDURE full_abst
181	END INTERFACE !
182
183
184
185	INTERFACE PRINT
186	MODULE PROCEDURE PRINT6
187	!!!!
188	MODULE PROCEDURE PRINT_DASPIN
189	MODULE PROCEDURE PRINT_probe8
190	MODULE PROCEDURE PRINT_spinor_8
191	MODULE PROCEDURE print_res_spinor_8
192	MODULE PROCEDURE print_rf_phasor_8
193	END INTERFACE
194
195	INTERFACE daPRINT
196	MODULE PROCEDURE PRINT6
197	!!!!
198	END INTERFACE
199
200	INTERFACE READ
201	MODULE PROCEDURE READ_DASPIN
202	MODULE PROCEDURE read_probe8 ! a bit illegal : reading polymorphs as taylor...
203	MODULE PROCEDURE read_spinor_8
204	END INTERFACE
205
206	INTERFACE ALLOC
207	MODULE PROCEDURE ALLOC_SPINOR_8
208	MODULE PROCEDURE ALLOC_probe_8
209	MODULE PROCEDURE ALLOC_daspin
210	MODULE PROCEDURE A_OPT_damap
211	MODULE PROCEDURE alloc_normal_spin
212	MODULE PROCEDURE alloc_res_SPINOR_8
213	MODULE PROCEDURE ALLOC_rf_phasor_8
214	MODULE PROCEDURE SET_TREE
215	END INTERFACE
216
217	INTERFACE KILL
218	MODULE PROCEDURE KILL_SPINOR_8
219	MODULE PROCEDURE KILL_probe_8
220	MODULE PROCEDURE KILL_DASPIN
221	MODULE PROCEDURE K_OPT_damap
222	MODULE PROCEDURE KILL_normal_spin
223	MODULE PROCEDURE KILL_res_SPINOR_8
224	MODULE PROCEDURE KILL_rf_phasor_8
225	END INTERFACE
226
227	INTERFACE ALLOC_33
228	MODULE PROCEDURE ALLOC_33t
229	MODULE PROCEDURE ALLOC_33p
230	END INTERFACE
231
232	INTERFACE ALLOC_nn
233	MODULE PROCEDURE ALLOC_33t
234	MODULE PROCEDURE ALLOC_33p
235	END INTERFACE
236
237	INTERFACE matmul_nn
238	MODULE PROCEDURE matmul_33
239	END INTERFACE
240
241	INTERFACE zero_33
242	MODULE PROCEDURE zero_33t
243	MODULE PROCEDURE zero_33p
244	END INTERFACE
245
246	INTERFACE zero_nn
247	MODULE PROCEDURE zero_33t
248	MODULE PROCEDURE zero_33p
249	END INTERFACE
250
251	INTERFACE KILL_nn
252	MODULE PROCEDURE KILL_33t
253	MODULE PROCEDURE KILL_33p
254	END INTERFACE
255
256	INTERFACE KILL_33
257	MODULE PROCEDURE KILL_33t
258	MODULE PROCEDURE KILL_33p
259	END INTERFACE
260
261	INTERFACE track
262	MODULE PROCEDURE track_TREE
263	MODULE PROCEDURE track_TREEP
264	END INTERFACE
265
266
267	INTERFACE trackg
268	MODULE PROCEDURE track_TREE_G
269	MODULE PROCEDURE track_TREEP_g
270	END INTERFACE
271
272
273	INTERFACE KILL
274	MODULE PROCEDURE KILL_TREE
275	MODULE PROCEDURE KILL_TREE_N
276	MODULE PROCEDURE KILL_33t
277	MODULE PROCEDURE KILL_33p
278	END INTERFACE
279
280	INTERFACE factor
281	MODULE PROCEDURE factor_am
282	MODULE PROCEDURE factor_as
283	END INTERFACE
284
285	INTERFACE ass
286	MODULE PROCEDURE assmap
287	MODULE PROCEDURE assprobe_8
288	END INTERFACE
289
290
291	CONTAINS
292	SUBROUTINE A_OPT_damap(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10)
293	implicit none
294	type (damapspin),INTENT(INout)::S1,S2
295	type (damapspin),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
296	call alloc(s1)
297	call alloc(s2)
298	if(present(s3)) call alloc(s3)
299	if(present(s4)) call alloc(s4)
300	if(present(s5)) call alloc(s5)
301	if(present(s6)) call alloc(s6)
302	if(present(s7)) call alloc(s7)
303	if(present(s8)) call alloc(s8)
304	if(present(s9)) call alloc(s9)
305	if(present(s10))call alloc(s10)
306	END SUBROUTINE A_opt_damap
307
308	SUBROUTINE K_OPT_damap(S1,S2,s3,s4,s5,s6,s7,s8,s9,s10)
309	implicit none
310	type (damapspin),INTENT(INout)::S1,S2
311	type (damapspin),optional, INTENT(INout):: s3,s4,s5,s6,s7,s8,s9,s10
312	call KILL(s1)
313	call KILL(s2)
314	if(present(s3)) call KILL(s3)
315	if(present(s4)) call KILL(s4)
316	if(present(s5)) call KILL(s5)
317	if(present(s6)) call KILL(s6)
318	if(present(s7)) call KILL(s7)
319	if(present(s8)) call KILL(s8)
320	if(present(s9)) call KILL(s9)
321	if(present(s10))call KILL(s10)
322	END SUBROUTINE K_OPT_damap
323
324	!!! use to be in extend_poly
325
326	FUNCTION daddsco( S1, S2 )
327	implicit none
328	TYPE (real_8) daddsco(ndd)
329	TYPE (damap), INTENT (IN) :: S1
330	real(dp) , INTENT (IN) :: S2(ndd)
331	integer localmaster,iia(4),ico(4),nd2,i
332	call liepeek(iia,ico)
333	nd2=iia(4)
334
335
336	do i=1,nd2
337	localmaster=master
338	call ass(daddsco(i))
339	daddsco(i)=s1%v(i)+s2(i)-(s1%v(i).sub.'0')
340	master=localmaster
341	enddo
342	do i=nd2+1,ndd
343	localmaster=master
344	call ass(daddsco(i))
345	if(nd2==4.and.(c_%npara==5.or.c_%npara==8).and.i==5) then
346	daddsco(i)=s2(i)+(1.0_dp.mono.'00001')
347	else
348	daddsco(i)=s2(i)
349	endif
350	master=localmaster
351	enddo
352
353	END FUNCTION daddsco
354
355	FUNCTION scdaddo( S2,S1 )
356	implicit none
357	TYPE (real_8) scdaddo(ndd)
358	TYPE (damap), INTENT (IN) :: S1
359	real(dp) , INTENT (IN) :: S2(ndd)
360	integer localmaster,iia(4),ico(4),nd2,i
361	call liepeek(iia,ico)
362	nd2=iia(4)
363
364	do i=1,nd2
365	localmaster=master
366	call ass(scdaddo(i))
367	scdaddo(i)=s1%v(i)+s2(i)-(s1%v(i).sub.'0')
368	master=localmaster
369	enddo
370	do i=nd2+1,ndd
371	localmaster=master
372	call ass(scdaddo(i))
373	if(nd2==4.and.(c_%npara==5.or.c_%npara==8).and.i==5) then
374	scdaddo(i)=s2(i)+(1.0_dp.mono.'00001')
375	else
376	scdaddo(i)=s2(i)
377	endif
378	master=localmaster
379	enddo
380
381
382	END FUNCTION scdaddo
383
384	SUBROUTINE REAL6real_8(S2,S1)
385	implicit none
386	real(dp),INTENT(inOUT)::S2(ndd)
387	type (real_8),INTENT(IN)::S1(ndd)
388	integer i
389
390
391	do i=1,ndd
392	s2(i)=s1(i) !%t
393	enddo
394	END SUBROUTINE REAL6real_8
395
396	SUBROUTINE real_8REAL_8(S1,S2)
397	implicit none
398	type (real_8),INTENT(in)::S2(ndd)
399	type (real_8),INTENT(inOUT)::S1(ndd)
400	integer i
401
402
403	do i=1,ndd
404	s1(i)=s2(i)
405	enddo
406	END SUBROUTINE real_8REAL_8
407
408	SUBROUTINE real_8REAL6(S1,S2)
409	implicit none
410	real(dp),INTENT(in)::S2(ndd)
411	type (real_8),INTENT(inOUT)::S1(ndd)
412	integer i
413
414
415	do i=1,ndd
416	s1(i)=s2(i)
417	enddo
418	END SUBROUTINE real_8REAL6
419
420
421
422	SUBROUTINE spimat_spinmat(S1,S2)
423	implicit none
424	type (spinmatrix),INTENT(in)::S2
425	type (spinmatrix),INTENT(inOUT)::S1
426	integer i,j
427
428	do i=1,3
429	do j=1,3
430	s1%s(i,j)=s2%s(i,j)
431	enddo
432	enddo
433	END SUBROUTINE spimat_spinmat
434
435	SUBROUTINE print6(S1,mf)
436	implicit none
437	type (real_8),INTENT(INout)::S1(ndd)
438	integer mf,i
439
440	do i=1,ndd
441	call print(s1(i),mf)
442	enddo
443
444	END SUBROUTINE print6
445
446	!!! end of "use to be in extend_poly"
447
448	SUBROUTINE COPY_TREE(T,U)
449	IMPLICIT NONE
450	TYPE(TREE_ELEMENT), INTENT(IN) :: T
451	TYPE(TREE_ELEMENT), INTENT(INOUT) :: U
452
453	IF(.NOT.ASSOCIATED(T%CC).OR..NOT.ASSOCIATED(U%CC) ) RETURN
454	IF(SIZE(T%CC)/=SIZE(U%CC) ) STOP 888
455
456	U%CC=T%CC
457	U%JL=T%JL
458	U%JV=T%JV
459	U%N=T%N
460	U%ND2=T%ND2
461	U%no=T%no
462
463	END SUBROUTINE COPY_TREE
464
465	SUBROUTINE COPY_TREE_N(T,U)
466	IMPLICIT NONE
467	TYPE(TREE_ELEMENT), INTENT(IN) :: T(:)
468	TYPE(TREE_ELEMENT), INTENT(INOUT) :: U(:)
469	INTEGER I
470
471	DO I=1,SIZE(T)
472	CALL COPY_TREE(T(I),U(I))
473	ENDDO
474
475	END SUBROUTINE COPY_TREE_N
476
477
478
479	SUBROUTINE NULL_TREE(T)
480	IMPLICIT NONE
481	TYPE(TREE_ELEMENT), INTENT(INOUT) :: T
482
483	NULLIFY(T%CC,T%JL,T%JV,T%N,T%ND2,T%no)
484
485	END SUBROUTINE NULL_TREE
486
487
488	SUBROUTINE ALLOC_TREE(T,N,ND2)
489	IMPLICIT NONE
490	TYPE(TREE_ELEMENT), INTENT(INOUT) :: T
491	INTEGER , INTENT(IN) :: N,ND2
492
493	!IF(N==0) RETURN
494
495	ALLOCATE(T%CC(N),T%fix(nd2),T%JL(N),T%JV(N),T%N,T%ND2,T%no)
496	T%N=N
497	T%ND2=ND2
498	T%no=0
499	T%fix=0.0_dp
500	END SUBROUTINE ALLOC_TREE
501
502	SUBROUTINE SET_TREE(T,MA)
503	IMPLICIT NONE
504	TYPE(TREE_ELEMENT), INTENT(INOUT) :: T
505	TYPE(DAMAP), INTENT(INOUT) :: MA
506	INTEGER N
507	TYPE(DAMAP) M
508
509	CALL ALLOC(M)
510
511	CALL mtree(MA%v%I,C_%ND2,M%v%I,C_%ND2)
512	CALL ppushGETN(M%v%I,C_%ND2,N)
513
514
515	CALL ALLOC_TREE(T,N,C_%ND2)
516	T%no=c_%no
517
518	CALL ppushstore(M%v%I,C_%ND2,T%CC,T%JL,T%JV)
519
520	CALL KILL(M)
521
522	END SUBROUTINE SET_TREE
523
524	! FOR FAST B FIELD IN PACKAGE OF PTC
525	SUBROUTINE SET_TREE_G(T,MA)
526	IMPLICIT NONE
527	TYPE(TREE_ELEMENT), INTENT(INOUT) :: T
528	TYPE(taylor), INTENT(INOUT) :: MA(:)
529	INTEGER N,NP
530	TYPE(taylor), ALLOCATABLE :: M(:)
531
532	NP=SIZE(MA)
533
534	ALLOCATE(M(NP))
535	CALL ALLOC(M,NP)
536
537	CALL mtree(MA%I,NP,M%I,NP)
538	CALL ppushGETN(M%I,NP,N)
539
540
541	CALL ALLOC_TREE(T,N,NP)
542	T%no=c_%no
543
544	CALL ppushstore(M%I,NP,T%CC,T%JL,T%JV)
545
546	CALL KILL(M,NP)
547	deallocate(M)
548	END SUBROUTINE SET_TREE_G
549
550	SUBROUTINE track_TREE_G(T,XI)
551	use da_arrays
552	IMPLICIT NONE
553	TYPE(TREE_ELEMENT), INTENT(IN) :: T
554	REAL(DP), INTENT(INOUT) :: XI(:)
555	REAL(DP) XT(lno),XF(lnv),XM(lno+1),XX
556	INTEGER JC,I,IV
557
558	XT=0.0_dp
559	XF=0.0_dp
560	XM=0.0_dp
561
562	do i=1,T%ND2
563	xt(i)=xi(i)
564	enddo
565	do i=1,T%ND2
566	xf(i) = T%cc(i)
567	enddo
568
569	XM(1) = 1.0_dp
570	JC=T%ND2
571	do i=1,(T%N-T%ND2)/T%ND2
572	!
573	xx = xm(T%jl(JC+1))*xt(T%jV(JC+1))
574	xm(T%jl(JC+1)+1) = xx
575	!
576	do iv=1,T%ND2
577	jc=jc+1
578	xf(iv) = xf(iv) + t%cc(jc) * xx
579	enddo
580	enddo
581	xi=xf
582
583
584	END SUBROUTINE track_TREE_G
585
586
587
588	SUBROUTINE track_TREEP_g(T,XI)
589	use da_arrays
590	IMPLICIT NONE
591	TYPE(TREE_ELEMENT), INTENT(IN) :: T
592	TYPE(REAL_8), INTENT(INOUT) :: XI(:)
593	TYPE(REAL_8) XT(lno),XF(lnv),XM(lno+1),XX
594	INTEGER JC,I,IV
595
596	CALL ALLOC(XT,lno)
597	CALL ALLOC(XF,lnv)
598	CALL ALLOC(XM,lno+1)
599	CALL ALLOC(XX)
600
601
602
603
604	do i=1,T%ND2
605	xt(i)=xi(i)
606	enddo
607	do i=1,T%ND2
608	xf(i) = T%cc(i)
609	enddo
610
611	XM(1) = 1.0_dp
612	JC=T%ND2
613
614	do i=1,(T%N-T%ND2)/T%ND2
615	!
616	xx = xm(T%jl(JC+1))*xt(T%jV(JC+1))
617	xm(T%jl(JC+1)+1) = xx
618	!
619	do iv=1,T%ND2
620	jc=jc+1
621	xf(iv) = xf(iv) + t%cc(jc) * xx
622	enddo
623	enddo
624
625	do i=1,T%ND2
626	xI(i)=xF(i)
627	enddo
628
629	CALL KILL(XT,lno)
630	CALL KILL(XF,lnv)
631	CALL KILL(XM,lno+1)
632	CALL KILL(XX)
633
634	END SUBROUTINE track_TREEP_g
635
636
637
638
639
640	! END OF FAST B FIELD IN PACKAGE OF PTC
641
642	! type tree_element
643	! real(dp) , DIMENSION(:), POINTER :: CC
644	! integer, DIMENSION(:), POINTER :: JL,JV
645	! INTEGER,POINTER :: N,ND2
646	! end type tree_element
647
648
649	SUBROUTINE KILL_TREE(T)
650	IMPLICIT NONE
651	TYPE(TREE_ELEMENT), INTENT(INOUT) :: T
652
653
654	IF(ASSOCIATED(T%CC)) DEALLOCATE(T%CC,T%fix,T%JL,T%JV,T%N,T%ND2,T%No)
655
656
657	END SUBROUTINE KILL_TREE
658
659	SUBROUTINE KILL_TREE_N(T)
660	IMPLICIT NONE
661	TYPE(TREE_ELEMENT), INTENT(INOUT) :: T(:)
662	INTEGER I
663
664	DO I=1,SIZE(T)
665	CALL KILL(T(I))
666	ENDDO
667
668	END SUBROUTINE KILL_TREE_N
669
670
671
672
673	SUBROUTINE track_TREE(T,XI,n)
674	use da_arrays
675	IMPLICIT NONE
676	TYPE(TREE_ELEMENT), INTENT(IN) :: T
677	REAL(DP), INTENT(INOUT) :: XI(6)
678	integer, optional, INTENT(IN) :: n
679	integer n1,k
680	REAL(DP) XT(lno),XF(6),XM(lno+1),XX
681	INTEGER JC,I,IV
682	xi(1:t%nd2)=xi(1:t%nd2)-t%fix
683	n1=1
684	if(present(n)) n1=n
685	do k=1,n1
686	if(.not.c_%CHECK_STABLE) return
687	XT=0.0_dp
688	XF=0.0_dp
689	XM=0.0_dp
690
691	do i=1,T%ND2
692	xt(i)=xi(i)
693	enddo
694	do i=1,T%ND2
695	xf(i) = T%cc(i)
696	enddo
697
698	XM(1) = 1.0_dp
699	JC=T%ND2
700	do i=1,(T%N-T%ND2)/T%ND2
701	!
702	xx = xm(T%jl(JC+1))*xt(T%jV(JC+1))
703	xm(T%jl(JC+1)+1) = xx
704	!
705	do iv=1,T%ND2
706	jc=jc+1
707	xf(iv) = xf(iv) + t%cc(jc) * xx
708	enddo
709	enddo
710	xi=xf
711
712	if(abs(xi(1))>c_%absolute_aperture.or.abs(xi(3))>c_%absolute_aperture) then
713	c_%CHECK_STABLE=.FALSE.
714	xlost=xi
715	endif
716	enddo
717
718	END SUBROUTINE track_TREE
719
720	SUBROUTINE track_TREEP(T,XI,n)
721	use da_arrays
722	IMPLICIT NONE
723	TYPE(TREE_ELEMENT), INTENT(IN) :: T
724	TYPE(REAL_8), INTENT(INOUT) :: XI(6)
725	integer, optional, INTENT(IN) :: n
726	integer n1,k
727	TYPE(REAL_8) XT(lno),XF(6),XM(lno+1),XX
728	INTEGER JC,I,IV
729
730	n1=1
731	if(present(n)) n1=n
732	do k=1,n1
733
734	CALL ALLOC(XT,lno)
735	CALL ALLOC(XF,6)
736	CALL ALLOC(XM,lno+1)
737	CALL ALLOC(XX)
738
739
740
741
742	do i=1,T%ND2
743	xt(i)=xi(i)
744	enddo
745	do i=1,T%ND2
746	xf(i) = T%cc(i)
747	enddo
748
749	XM(1) = 1.0_dp
750	JC=T%ND2
751
752	do i=1,(T%N-T%ND2)/T%ND2
753	!
754	xx = xm(T%jl(JC+1))*xt(T%jV(JC+1))
755	xm(T%jl(JC+1)+1) = xx
756	!
757	do iv=1,T%ND2
758	jc=jc+1
759	xf(iv) = xf(iv) + t%cc(jc) * xx
760	enddo
761	enddo
762
763	do i=1,T%ND2
764	xI(i)=xF(i)
765	enddo
766
767	CALL KILL(XT,lno)
768	CALL KILL(XF,6)
769	CALL KILL(XM,lno+1)
770	CALL KILL(XX)
771
772	enddo
773	END SUBROUTINE track_TREEP
774
775	! READING COSY MAPS AND UNIVERSAL TAYLORS
776	SUBROUTINE dainput_SPECIAL6(S1,MFILE,COSY)
777	implicit none
778	INTEGER,INTENT(in)::MFILE,COSY
779	type (damap),INTENT(INOUT)::S1
780	INTEGER I,js(6),k1,k2,l,ncoef
781	real(dp) x
782	character*200 line
783	S1=0
784	SELECT CASE(COSY)
785	CASE(0:1)
786	call dainput(s1,mfile)
787	CASE(2) ! COSY INFINITY
788	do i=1,c_%nd2
789	read(mfile,'(a200)') line
790	read(mfile,'(a200)') line
791	do while(line(14:16)/='---')
792	read(line,*) k1,x,k2,js
793	s1%v(i)=s1%v(i)+(x.mono.js)
794	read(mfile,'(a200)') line
795	enddo
796	enddo
797
798	CASE(-1) ! sagan
799	do i=1,c_%nd2
800	read(mfile,*) ncoef
801	do l=1,ncoef
802	read(mfile,*) x,js
803	s1%v(i)=s1%v(i)+(x.mono.js)
804	enddo
805	enddo
806
807	CASE DEFAULT
808
809	WRITE(6,*) " NOT SUPPORTED IN DAREADMAP_SPECIAL"
810	STOP 111
811
812	END SELECT
813
814
815	END SUBROUTINE dainput_SPECIAL6
816
817	! SYMPLECTIFY A MAP NEAR THE IDENTITY
818	SUBROUTINE symplectic(m,eps,nst)
819	IMPLICIT NONE
820	type(damap), INTENT(INOUT) :: m
821	type(onelieexponent) uno
822	type(damap) id
823	real(dp), optional :: eps
824	integer nst
825
826	call alloc(uno); call alloc(id);
827
828	if(present(eps))then
829	if(eps>0.0_dp) uno%eps=eps
830	endif
831	uno=m
832	id=1
833	uno%pb%h=uno%pb%h/nst
834	m=texp(uno%pb,id)
835
836	call kill(uno); call kill(id);
837
838
839	end SUBROUTINE symplectic
840
841	! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
842	! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
843	! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
844	! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
845	! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
846	! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
847	! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
848	! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
849	! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
850	! SPIN STUFF IS HERE
851
852
853
854
855	subroutine assmap(s1)
856	implicit none
857	TYPE (damapspin) s1
858	integer i,j
859
860	select case(master)
861	case(0:ndumt-1)
862	master=master+1
863	case(ndumt)
864	w_p=0
865	w_p%nc=1
866	w_p=(/" cannot indent anymore "/)
867	w_p%fc='(1((1X,A72),/))'
868	! call !write_e(100)
869	end select
870
871	do i=1,3
872	do j=1,3
873	call assp_no_master(s1%s%s(i,j))
874	! call ass0(s1%s(i,j)%t) ! taylor
875	! s1%s(i,j)%alloc=my_true
876	! s1%s(i,j)%kind=1
877	! s1%s(i,j)%i=0
878	enddo
879	enddo
880	do i=1,c_%nd2
881	call ass0(s1%m%v(i))
882	enddo
883
884	end subroutine assmap
885
886	subroutine assprobe_8(s1)
887	implicit none
888	TYPE (probe_8) s1
889	integer i,j
890
891	select case(master)
892	case(0:ndumt-1)
893	master=master+1
894	case(ndumt)
895	w_p=0
896	w_p%nc=1
897	w_p=(/" cannot indent anymore "/)
898	w_p%fc='(1((1X,A72),/))'
899	! call !write_e(100)
900	end select
901
902	! if(C_%SPIN_POS/=0) then
903	! do i=1,3
904	! call ass0(s1%s%x(i)%t)
905	! s1%s%x(i)%alloc=my_true
906	! s1%s%x(i)%kind=2
907	! s1%s%x(i)%i=0
908	! enddo
909	do j=1,3
910	do i=1,3
911	call assp_no_master(s1%s(j)%x(i))
912	! call ass0(s1%s(j)%x(i)%t)
913	! s1%s(j)%x(i)%alloc=my_true
914	! s1%s(j)%x(i)%kind=1
915	! s1%s(j)%x(i)%i=0
916	enddo
917	enddo
918
919	! endif
920	do i=1,6 !
921	call assp_no_master(s1%x(i))
922	! call ass0(s1%x(i)%t)
923	! s1%x(i)%alloc=my_true
924	! s1%x(i)%kind=1
925	! s1%x(i)%i=0
926	enddo
927
928	do i=1,2 !
929	call assp_no_master(s1%AC%x(i))
930	enddo
931
932
933	end subroutine assprobe_8
934
935
936
937
938	FUNCTION FULL_ABST( S1 )
939	implicit none
940	REAL(DP) FULL_ABST
941	TYPE (damapspin), INTENT (IN) :: S1
942	integer i,j
943
944	FULL_ABST=0.0_dp
945
946	do i=1,3
947	do j=1,3
948	FULL_ABST=FULL_ABS(S1%s%s(i,j))+FULL_ABST
949	enddo
950	enddo
951
952
953	END FUNCTION FULL_ABST
954
955	FUNCTION CUTORDER( S1, S2 )
956	implicit none
957	TYPE (damapspin) CUTORDER
958	TYPE (damapspin), INTENT (IN) :: S1
959	integer , INTENT (IN) :: S2
960	integer localmaster
961	integer i,j
962
963	localmaster=master
964	call ass(CUTORDER)
965
966	CUTORDER=0
967
968	do i=1,3
969	do j=1,3
970	CUTORDER%s%s(i,j)=s1%s%s(i,j).cut.(s2-1)
971	enddo
972	enddo
973
974	CUTORDER%m=s1%m.cut.s2
975
976	master=localmaster
977
978	END FUNCTION CUTORDER
979
980	subroutine inv_damapspin(r,s)
981	implicit none
982	TYPE(damapspin), INTENT(INOUT) :: S
983	TYPE(damapspin), INTENT(IN) :: r
984	TYPE(damapspin) t
985	INTEGER I,j
986
987	call alloc(t)
988
989	t=r
990	s%m=t%m**(-1)
991
992	do i=1,3
993	do j=1,3
994	s%s%s(i,j)=t%s%s(j,i)
995	! if(s%s(i,j)%kind==2) then ! taylor
996	s%s%s(i,j)=s%s%s(i,j)*s%m
997	! s%s(i,j)=s%s(i,j)%t*s%m ! taylor
998	! endif ! taylor
999	enddo
1000	enddo
1001
1002	call kill(t)
1003
1004	END subroutine inv_damapspin
1005
1006
1007	! type damapspin
1008	! REAL(DP) X(6)
1009	! type(damap) M
1010	! type(real_8) s(3,3)
1011	! end type damapspin
1012
1013
1014
1015
1016	subroutine EQUAL_damapspin(S,R)
1017	implicit none
1018	TYPE(damapspin), INTENT(INOUT) :: S
1019	TYPE(damapspin), INTENT(IN) :: r
1020	INTEGER I,j
1021
1022	do i=1,3
1023	do j=1,3
1024	s%s%s(i,j)=r%s%s(i,j)
1025	enddo
1026	enddo
1027
1028	s%m=r%m
1029	s%e_ij=r%e_ij
1030	! s%s0=r%s0
1031
1032	! do i=1,6
1033	! s%x(i)=r%x(i)
1034	! enddo
1035
1036	END subroutine EQUAL_damapspin
1037
1038
1039	subroutine EQUAL_damapspin_smat(S,R)
1040	implicit none
1041	TYPE(damapspin), INTENT(INOUT) :: S
1042	real(dp), INTENT(IN) :: r(3,3)
1043	INTEGER I,j
1044
1045	do i=1,3
1046	do j=1,3
1047	s%s%s(i,j)=r(i,j)
1048	enddo
1049	enddo
1050
1051
1052	END subroutine EQUAL_damapspin_smat
1053
1054	subroutine EQUAL_smat_damapspin(R,S)
1055	implicit none
1056	real(dp), INTENT(INout) :: r(3,3)
1057	TYPE(damapspin), INTENT(IN) :: S
1058	INTEGER I,j
1059
1060	do i=1,3
1061	do j=1,3
1062	r(i,j)=s%s%s(i,j)
1063	enddo
1064	enddo
1065
1066
1067	END subroutine EQUAL_smat_damapspin
1068
1069
1070	subroutine EQUAL_IDENTITY_SPINOR_8(S,R)
1071	implicit none
1072	TYPE(SPINOR_8), INTENT(INOUT) :: S
1073	INTEGER, INTENT(IN) :: R
1074	INTEGER I
1075
1076	! S%G=A_PARTICLE
1077	IF(R==1.or.r==2.or.r==3) THEN
1078	DO I=1,3 !C_%NSPIN
1079	S%X(I)=0.0_dp
1080	enddo
1081	S%X(r)=1.0_dp
1082	! write(6,*) " in EQUAL_IDENTITY_SPINOR_8"
1083	! stop 888
1084	! DO I=1,C_%NSPIN
1085	! S%X(I)=ONE.MONO.(c_%npara_fpp-C_%NSPIN+I)
1086	! ENDDO
1087	ELSEIF(R==0) THEN
1088	DO I=1,3 !C_%NSPIN
1089	S%X(I)=0.0_dp
1090	enddo
1091	ELSE
1092	write(6,*) " stopped in EQUAL_IDENTITY_SPINOR_8"
1093	STOP 100
1094	ENDIF
1095
1096	END subroutine EQUAL_IDENTITY_SPINOR_8
1097
1098
1099	subroutine EQUAL_IDENTITY_SPINOR_8_r3(S,R)
1100	implicit none
1101	TYPE(SPINOR_8), INTENT(INOUT) :: S
1102	real(dp), INTENT(IN) :: R(3)
1103	INTEGER I
1104
1105	! S%G=A_PARTICLE
1106
1107	DO I=1,3 !C_%NSPIN
1108	S%X(I)=r(i) !+(ONE.MONO.(c_%npara_fpp-C_%NSPIN+I))
1109	ENDDO
1110
1111	END subroutine EQUAL_IDENTITY_SPINOR_8_r3
1112
1113	subroutine EQUAL_IDENTITY_SPINOR (S,R)
1114	implicit none
1115	TYPE(SPINOR), INTENT(INOUT) :: S
1116	INTEGER, INTENT(IN) :: R
1117	INTEGER I
1118
1119	! S%G=A_PARTICLE
1120	IF(R==1.or.r==2.or.r==3) THEN
1121
1122	DO I=1,3 !C_%NSPIN
1123	S%X(I)=0.0_dp
1124	enddo
1125	S%X(r)=1.0_dp
1126	ELSEIF(R==0) THEN
1127	DO I=1,3 !C_%NSPIN
1128	S%X(I)=0.0_dp
1129	enddo
1130	ELSE
1131	STOP 100
1132	ENDIF
1133
1134	END subroutine EQUAL_IDENTITY_SPINOR
1135
1136	subroutine EQUAL_PROBE_REAL6 (P,X)
1137	implicit none
1138	TYPE(PROBE), INTENT(INOUT) :: P
1139	REAL(DP), INTENT(IN) :: X(6)
1140	INTEGER I
1141	P%u=my_false
1142	! P%s(0)%x=zero
1143	! P%s(0)%x(n0_normal)=one
1144
1145	DO I=1,ISPIN1R
1146	P%s(i)%x=0.0_dp
1147	P%s(i)%x(i)=1.0_dp
1148	enddo
1149	P%X=X
1150	P%ac%t=0.0_dp
1151
1152	END subroutine EQUAL_PROBE_REAL6
1153
1154	subroutine EQUAL_PROBE8_REAL6 (P,X)
1155	implicit none
1156	TYPE(PROBE_8), INTENT(INOUT) :: P
1157	REAL(DP), INTENT(IN) :: X(6)
1158	INTEGER I
1159
1160	P%u=my_false
1161	! P%S=0
1162
1163	! P%s(0)=0
1164	! P%s(0)%x(n0_normal)=one
1165
1166	DO I=1,3
1167	P%s(i)=0
1168	P%s(i)%x(i)=1.0_dp
1169	enddo
1170
1171	DO I=1,6
1172	P%X(i)=X(i)
1173	enddo
1174
1175	P%AC%X(1)=0.0_dp
1176	P%AC%X(2)=0.0_dp
1177	P%AC%t=0.0_dp
1178	p%e_ij=0.0_dp
1179	END subroutine EQUAL_PROBE8_REAL6
1180
1181	subroutine EQUAL_PROBE8_PROBE8(P8,P)
1182	implicit none
1183	TYPE(PROBE_8), INTENT(IN) :: P
1184	TYPE(PROBE_8), INTENT(INOUT) :: P8
1185	INTEGER I,J
1186
1187	DO I=1,6
1188	P8%X(I)=P%X(I)
1189	ENDDO
1190
1191	DO I=1,6
1192	DO j=1,6
1193	P8%E_ij(I,j)=P%E_ij(I,j)
1194	ENDDO
1195	ENDDO
1196
1197	! P8%S=P%S
1198	DO I=1,3
1199	P8%S(i)=P%S(i)
1200	enddo
1201	P8%AC=P%AC
1202
1203	P8%u=P%u
1204
1205	END subroutine EQUAL_PROBE8_PROBE8
1206
1207
1208
1209	subroutine EQUAL_RF8_RF8(P8,P)
1210	implicit none
1211	TYPE(rf_phasor_8), INTENT(IN) :: P
1212	TYPE(rf_phasor_8), INTENT(INOUT) :: P8
1213	INTEGER I
1214
1215	DO I=1,2
1216	P8%X(I)=P%X(I)
1217	ENDDO
1218	P8%om=P%om
1219	P8%t=P%t
1220
1221	END subroutine EQUAL_RF8_RF8
1222
1223	subroutine EQUAL_RF8_RF(P8,P)
1224	implicit none
1225	TYPE(rf_phasor_8), INTENT(INOUT) :: P8
1226	TYPE(rf_phasor), INTENT(IN) :: P
1227	INTEGER I
1228
1229	DO I=1,2
1230	P8%X(I)=P%X(I)
1231	ENDDO
1232	P8%om=P%om
1233	P8%t=P%t
1234
1235	END subroutine EQUAL_RF8_RF
1236
1237	subroutine EQUAL_RF_RF8(P,P8)
1238	implicit none
1239	TYPE(rf_phasor), INTENT(INOUT) :: P
1240	TYPE(rf_phasor_8), INTENT(IN) :: P8
1241	INTEGER I
1242
1243	DO I=1,2
1244	P%X(I)=P8%X(I)
1245	ENDDO
1246	P%om=P8%om
1247	P%t=P8%t
1248
1249	END subroutine EQUAL_RF_RF8
1250
1251
1252
1253	subroutine EQUAL_PROBE8_PROBE (P8,P)
1254	implicit none
1255	TYPE(PROBE), INTENT(IN) :: P
1256	TYPE(PROBE_8), INTENT(INOUT) :: P8
1257	INTEGER I
1258	DO I=1,6
1259	P8%X(I)=P%X(I)
1260	ENDDO
1261	! P8%S=P%S
1262	! do i=1,3
1263	! P8%S(i)=0
1264	! P8%S(i)%x(i)=one
1265	! enddo
1266	! P8%S(0)=P%s(0)
1267	do I=1,ISPIN1R
1268	P8%S(I)=P%s(I)
1269	enddo
1270
1271
1272	P8%u=P%u
1273	P8%e_ij=0.0_dp
1274	END subroutine EQUAL_PROBE8_PROBE
1275
1276	subroutine EQUAL_PROBE_PROBE8 (P,P8)
1277	implicit none
1278	TYPE(PROBE), INTENT(INOUT) :: P
1279	TYPE(PROBE_8), INTENT(IN) :: P8
1280	INTEGER I
1281	DO I=1,6
1282	P%X(I)=P8%X(I)
1283	ENDDO
1284	! P%S(0)=P8%S(0)
1285	DO I=1,ISPIN1R
1286	P%S(I)=P8%S(I)
1287	ENDDO
1288	P%u=P8%u
1289
1290	END subroutine EQUAL_PROBE_PROBE8
1291
1292	subroutine EQUAL_IDENTITY_probe(R,S)
1293	implicit none
1294	TYPE(probe), INTENT(INOUT) :: R
1295	INTEGER, INTENT(IN) :: S
1296	INTEGER I
1297
1298	R%S(1)=0
1299	R%S(2)=0
1300	R%S(3)=0
1301	DO I=1,6
1302	R%X(I)=0.0_dp
1303	ENDDO
1304	IF(S==1) THEN
1305	! R%S(0)%x(n0_normal)=one
1306	R%S(1)=1
1307	R%S(2)=2
1308	R%S(3)=3
1309	ELSEif(s==0) then
1310	else
1311	STOP 100
1312	ENDIF
1313
1314	END subroutine EQUAL_IDENTITY_probe
1315
1316	subroutine EQUAL_IDENTITY_probe_8(R,S)
1317	implicit none
1318	TYPE(probe_8), INTENT(INOUT) :: R
1319	INTEGER, INTENT(IN) :: S
1320	INTEGER I
1321
1322	! R%S=S
1323	! R%S(0)=0
1324	R%S(1)=0
1325	R%S(2)=0
1326	R%S(3)=0
1327	DO I=1,6 !-C_%NSPIN
1328	R%X(I)=0.0_dp
1329	enddo
1330	IF(S==1) THEN
1331	DO I=1,c_%npara_fpp !-C_%NSPIN
1332	R%X(I)=1.0_dp.MONO.I
1333	ENDDO
1334	! R%S(0)%x(n0_normal)=one
1335	R%S(1)=1
1336	R%S(2)=2
1337	R%S(3)=3
1338	ELSEIF(S==0) THEN
1339
1340	! DO I=1,6 !-C_%NSPIN
1341	! R%X(I)=ZERO!
1342	! enddo
1343	ELSE
1344	STOP 100
1345	ENDIF
1346	R%e_ij=0.0_dp
1347	END subroutine EQUAL_IDENTITY_probe_8
1348
1349
1350	subroutine EQUAL_SPINOR8_SPINOR8(R,S)
1351	implicit none
1352	TYPE(SPINOR_8), INTENT(IN) :: S
1353	TYPE(SPINOR_8), INTENT(INOUT) :: R
1354	INTEGER I
1355
1356	DO I=1,3
1357	R%X(I)=S%X(I)
1358	ENDDO
1359	! R%G=S%G
1360	END subroutine EQUAL_SPINOR8_SPINOR8
1361
1362	subroutine EQUAL_SPINOR_SPINOR8(R,S)
1363	implicit none
1364	TYPE(SPINOR_8), INTENT(IN) :: S
1365	TYPE(SPINOR), INTENT(INOUT) :: R
1366	INTEGER I
1367
1368	DO I=1,3
1369	R%X(I)=S%X(I)
1370	ENDDO
1371	! R%G=S%G
1372	END subroutine EQUAL_SPINOR_SPINOR8
1373
1374	subroutine EQUAL_SPINOR8_SPINOR(R,S)
1375	implicit none
1376	TYPE(SPINOR), INTENT(IN) :: S
1377	TYPE(SPINOR_8), INTENT(INOUT) :: R
1378	INTEGER I
1379
1380	DO I=1,3
1381	R%X(I)=S%X(I)
1382	ENDDO
1383	! R%G=S%G
1384	END subroutine EQUAL_SPINOR8_SPINOR
1385
1386	subroutine EQUAL_DAMAP_RAY8(DS,R)
1387	implicit none
1388	TYPE(probe_8), INTENT(IN) :: R
1389	TYPE(damap), INTENT(INOUT) :: DS
1390	INTEGER I,nd2t
1391
1392	nd2t=C_%ND2
1393	if(doing_ac_modulation_in_ptc) then
1394	nd2t=C_%ND2-2
1395	endif
1396
1397
1398
1399	DO I=1,nd2t
1400	DS%V(I)=R%X(I)
1401	ENDDO
1402	DO I=nd2t+1,C_%ND2
1403	DS%V(I)=R%ac%x(i-nd2t)
1404	ENDDO
1405
1406	END subroutine EQUAL_DAMAP_RAY8
1407
1408
1409	subroutine EQUAL_DAMAPSPIN_RAY8(DS,R)
1410	implicit none
1411	TYPE(probe_8), INTENT(IN) :: R
1412	TYPE(damapspin), INTENT(INOUT) :: DS
1413	real(dp) s(3,3)
1414	INTEGER I,J,nd2t
1415
1416	nd2t=C_%ND2
1417	if(doing_ac_modulation_in_ptc) then
1418	nd2t=C_%ND2-2
1419	endif
1420
1421
1422
1423	DO I=1,nd2t
1424	DS%M%V(I)=R%X(I)
1425	ENDDO
1426	DO I=nd2t+1,C_%ND2
1427	DS%M%V(I)=R%ac%x(i-nd2t)
1428	ENDDO
1429
1430	DO I=1,3
1431	DO J=1,3
1432	DS%S%s(I,J)=R%S(J)%X(I)
1433	ENDDO
1434	ENDDO
1435	! DS%S(I,1)=R%SX%X(I)
1436	! DS%S(I,2)=R%SY%X(I)
1437	! DS%S(I,3)=R%SZ%X(I)
1438	ds%e_ij=r%e_ij
1439
1440	! DO I=1,3
1441	! ! n(i)=R%S(0)%X(I)
1442	! DO J=1,3
1443	! S(I,J)=DS%S(I,J)
1444	! ENDDO
1445	! ENDDO
1446
1447	! call inv_as(s,s)
1448
1449	! n= matmul(s,n)
1450
1451	! ds%s0=n
1452
1453
1454
1455	! DS%G=R%S%G
1456	END subroutine EQUAL_DAMAPSPIN_RAY8
1457
1458	subroutine EQUAL_RAY8_DAMAPSPIN(R,DS)
1459	implicit none
1460	TYPE(probe_8), INTENT(INOUT) :: R
1461	TYPE(damapspin), INTENT(IN) :: DS
1462	INTEGER I,J
1463
1464	! if(C_%ND2/=6) then
1465	! write(6,*) " not implemented in EQUAL_DASPIN_RAY8"
1466	! stop
1467	! endif
1468
1469
1470
1471
1472	DO I=1,C_%ND2
1473	R%X(I)= MORPH(DS%M%V(I)) !+R%X(I)
1474	ENDDO
1475
1476	DO J=1,3
1477	DO I=1,3
1478	R%S(J)%X(I)=DS%S%s(I,J) ! taylor MORPH(DS%S(I,J))
1479	ENDDO
1480	ENDDO
1481
1482	r%e_ij=ds%e_ij
1483
1484
1485
1486	! R%S%G=DS%G
1487
1488	END subroutine EQUAL_RAY8_DAMAPSPIN
1489
1490	subroutine EQUAL_DAmapSPIN_int(DS,I1)
1491	implicit none
1492	TYPE(damapspin), INTENT(INOUT) :: DS
1493	INTEGER, INTENT(IN) :: I1
1494	INTEGER I
1495
1496	! type daspin
1497	! REAL(DP) X(6)
1498	! type(damap) M
1499	! type(real_8) s(3,3)
1500	! end type daspin
1501
1502	! DS%X=ZERO
1503	DS%M=I1
1504	CALL zero_33(DS%S%s)
1505
1506	DO I=1,3
1507	if(i1==1) then
1508	DS%S%s(I,I)=1.0_dp
1509	endif
1510	ENDDO
1511	DS%e_ij=0.0_dp
1512
1513	END SUBROUTINE EQUAL_DAmapSPIN_int
1514
1515	FUNCTION scdadd( S2,S1 )
1516	implicit none
1517	TYPE (probe_8) scdadd
1518	TYPE (damapspin), INTENT (IN) :: S1
1519	type(probe) , INTENT (IN) :: S2
1520	integer localmaster,i,j,nd2t,dc
1521	type(taylor) d
1522
1523
1524	scdadd%u=my_false
1525	scdadd%E_ij=0.0_dp
1526
1527	if(doing_ac_modulation_in_ptc) then
1528	dc=2
1529	else
1530	dc=0
1531	endif
1532	if(c_%ndpt/=0) then ! delta-l without cavity
1533	nd2t=6
1534	elseif(c_%ndpt==0.and.c_%nd2==4+dc) then ! no delta-l maybe delta as a parameter
1535	nd2t=4
1536	elseif(c_%ndpt==0.and.c_%nd2==6+dc) then ! delta-l with cavity
1537	nd2t=6
1538	else
1539	write(6,*) " internal PTC error in scdadd o_tree_element.f90 "
1540	stop 666
1541	endif
1542
1543	call alloc(d)
1544	do i=1,nd2t ! from 1-4 or 1-6 (if ndpt=0)
1545	localmaster=master
1546	call ass(scdadd%x(i))
1547	! scdadd%x(i)=s1%m%v(i)+s2%x(i)
1548	scdadd%x(i)=s1%M%V(i)+s2%x(i)-(s1%M%V(i).sub.'0')
1549	master=localmaster
1550	enddo
1551
1552	do i=nd2t+1,6
1553	localmaster=master
1554	call ass(scdadd%x(i))
1555	if((c_%npara==5+dc).AND.I==5) then ! npr
1556	scdadd%x(i)=s2%x(i)+(1.0_dp.mono.c_%npara)
1557	else
1558	scdadd%x(i)=s2%x(i)
1559	endif
1560	master=localmaster
1561	enddo
1562
1563	! if(doing_ac_modulation_in_ptc) then ! useless now
1564	localmaster=master
1565	call ass(scdadd%AC%x(1))
1566	! scdadd%x(i)=s1%m%v(i)+s2%x(i)
1567	scdadd%ac%x(1)=s1%M%V(C_%ND2-1) +s2%AC%x(1)
1568	master=localmaster
1569	localmaster=master
1570	call ass(scdadd%AC%x(2))
1571	! scdadd%x(i)=s1%m%v(i)+s2%x(i)
1572	scdadd%ac%x(2)=s1%M%V(C_%ND2) +s2%AC%x(2)
1573	master=localmaster
1574	localmaster=master
1575	call ass(scdadd%AC%om)
1576	! call ass(scdadd%AC%t)
1577	! scdadd%x(i)=s1%m%v(i)+s2%x(i)
1578	scdadd%AC%om=s2%AC%om
1579	scdadd%AC%t=s2%AC%t
1580	master=localmaster
1581	! endif
1582
1583	DO I=1,3
1584	! call ass(scdadd%s%x(i))
1585	DO J=1,3
1586	localmaster=master
1587	call ass(scdadd%s(J)%x(i))
1588	scdadd%s(J)%x(i)=S1%S%s(I,J)
1589	master=localmaster
1590	ENDDO
1591
1592	ENDDO
1593
1594	scdadd%e_ij=s1%e_ij
1595	END FUNCTION scdadd
1596
1597
1598	FUNCTION daddsc( S1,S2 )
1599	implicit none
1600	TYPE (probe_8) daddsc
1601	TYPE (damapspin), INTENT (IN) :: S1
1602	type(probe) , INTENT (IN) :: S2
1603	integer localmaster,i,j,nd2t,dc
1604	type(taylor) d
1605
1606
1607
1608
1609	! call ass(daddsc)
1610	daddsc%u=my_false
1611	daddsc%E_ij=0.0_dp
1612
1613	if(doing_ac_modulation_in_ptc) then
1614	dc=2
1615	else
1616	dc=0
1617	endif
1618	if(c_%ndpt/=0) then ! delta-l without cavity
1619	nd2t=6
1620	elseif(c_%ndpt==0.and.c_%nd2==4+dc) then ! no delta-l maybe delta as a parameter
1621	nd2t=4
1622	elseif(c_%ndpt==0.and.c_%nd2==6+dc) then ! delta-l with cavity
1623	nd2t=6
1624	else
1625	write(6,*) " internal PTC error in daddsc o_tree_element.f90 "
1626	stop 666
1627	endif
1628
1629	call alloc(d)
1630	do i=1,nd2t ! from 1-4 or 1-6 (if ndpt=0)
1631	localmaster=master
1632	call ass(daddsc%x(i))
1633	! scdadd%x(i)=s1%m%v(i)+s2%x(i)
1634	daddsc%x(i)=s1%M%V(i)+s2%x(i)-(s1%M%V(i).sub.'0')
1635	master=localmaster
1636	enddo
1637
1638	do i=nd2t+1,6
1639	localmaster=master
1640	call ass(daddsc%x(i))
1641	if((c_%npara==5+dc).AND.I==5) then ! npr
1642	daddsc%x(i)=s2%x(i)+(1.0_dp.mono.c_%npara)
1643	else
1644	daddsc%x(i)=s2%x(i)
1645	endif
1646	master=localmaster
1647	enddo
1648
1649	! if(doing_ac_modulation_in_ptc) then ! useless now
1650	localmaster=master
1651	call ass(daddsc%AC%x(1))
1652	! scdadd%x(i)=s1%m%v(i)+s2%x(i)
1653	daddsc%ac%x(1)=s1%M%V(C_%ND2-1) +s2%AC%x(1)
1654	master=localmaster
1655	localmaster=master
1656	call ass(daddsc%AC%x(2))
1657	! scdadd%x(i)=s1%m%v(i)+s2%x(i)
1658	daddsc%ac%x(2)=s1%M%V(C_%ND2) +s2%AC%x(2)
1659	master=localmaster
1660	localmaster=master
1661	call ass(daddsc%AC%om)
1662	! call ass(daddsc%AC%t)
1663	! scdadd%x(i)=s1%m%v(i)+s2%x(i)
1664	daddsc%AC%om=s2%AC%om
1665	daddsc%AC%t=s2%AC%t
1666	master=localmaster
1667	! endif
1668
1669	DO I=1,3
1670	! call ass(scdadd%s%x(i))
1671	DO J=1,3
1672	localmaster=master
1673	call ass(daddsc%s(J)%x(i))
1674	daddsc%s(J)%x(i)=S1%S%s(I,J)
1675	master=localmaster
1676	ENDDO
1677
1678	ENDDO
1679
1680	daddsc%e_ij=s1%e_ij
1681
1682
1683	END FUNCTION daddsc
1684
1685	subroutine print_DASPIN(DS,MF,prec)
1686	implicit none
1687	TYPE(damapspin), INTENT(INOUT) :: DS
1688	real(dp), optional :: prec
1689	INTEGER MF,I,J
1690	logical(lp) spin_in,rad_in
1691
1692	! WRITE(MF,*) " ORBIT "
1693	! WRITE(MF,*) DS%X(1:3)
1694	! WRITE(MF,*) DS%X(4:6)
1695	WRITE(MF,*) c_%nd2 ," DIMENSIONAL ORBITAL MAP "
1696	CALL PRINT(DS%M,MF,prec)
1697	call check_spin(DS,spin_in)
1698	if(spin_in) then
1699	WRITE(MF,*) " SPIN MATRIX "
1700	DO I=1,3
1701	DO J=1,3
1702	WRITE(MF,*) " SPIN MATRIX COMPONENT ",I,J
1703	CALL PRINT(DS%S%s(I,J),MF,prec)
1704	ENDDO
1705	ENDDO
1706	else
1707
1708	WRITE(MF,*) "NO SPIN MATRIX OR IDENTITY "
1709
1710	endif
1711	call check_rad(DS%e_ij,rad_in)
1712	if(rad_in) then
1713	WRITE(MF,*) " STOCHASTIC KICK "
1714	DO I=1,6
1715	DO J=1,6
1716	WRITE(MF,*) " STOCHASTIC KICK COMPONENT ",I,J
1717	write(mf,*) ds%e_ij(i,j)
1718	ENDDO
1719	ENDDO
1720	else
1721	WRITE(MF,*) "NO STOCHASTIC KICK "
1722	endif
1723
1724	END subroutine print_DASPIN
1725
1726	subroutine print_probe8(DS,MF)
1727	implicit none
1728	TYPE(probe_8), INTENT(INOUT) :: DS
1729	INTEGER MF,I,j
1730	logical(lp) rad_in
1731
1732	WRITE(MF,*) " ORBIT "
1733	do i=1,6
1734	write(mf,*) ' Variable ',i
1735	call print(ds%x(i),mf)
1736	enddo
1737	! WRITE(MF,*) " SPIN 0 "
1738	! do i=1,3
1739	! write(mf,*) ' Spin Variable ',i
1740	! call print(ds%s(0)%x(i),mf)
1741	! enddo
1742	WRITE(MF,*) " SPIN X "
1743	call print(ds%s(1),mf)
1744	! do i=1,3
1745	! write(mf,*) ' Spin Variable ',i
1746	! call print(ds%s(1)%x(i),mf)
1747	! enddo
1748	WRITE(MF,*) " SPIN Y "
1749	call print(ds%s(2),mf)
1750	! do i=1,3
1751	! write(mf,*) ' Spin Variable ',i
1752	! call print(ds%s(2)%x(i),mf)
1753	! enddo
1754	WRITE(MF,*) " SPIN Z "
1755	call print(ds%s(3),mf)
1756	! do i=1,3
1757	! write(mf,*) ' Spin Variable ',i
1758	! call print(ds%s(3)%x(i),mf)
1759	! enddo
1760
1761
1762	call check_rad(DS%e_ij,rad_in)
1763
1764
1765
1766	if(rad_in) then
1767
1768	WRITE(MF,*) " STOCHASTIC KICK "
1769
1770	do i=1,6
1771	do j=1,6
1772	write(6,*) i,j,ds%e_ij(i,j)
1773	enddo
1774	enddo
1775	else
1776
1777	WRITE(MF,*) "NO STOCHASTIC KICK "
1778
1779	endif
1780	if(doing_ac_modulation_in_ptc) then
1781	call print(ds%ac,mf)
1782	else
1783	WRITE(MF,*) "NO MODULATION "
1784
1785	endif
1786
1787	END subroutine print_probe8
1788
1789	subroutine print_rf_phasor_8(S,MF)
1790	implicit none
1791	TYPE(rf_phasor_8), INTENT(INOUT) :: s
1792	INTEGER MF,I
1793
1794	write(mf,*) ' AC INFORMATION '
1795	call print(s%om,mf)
1796	call print(s%t,mf)
1797	do i=1,2
1798	call print(s%x(i),mf)
1799	enddo
1800
1801	END subroutine print_rf_phasor_8
1802
1803
1804	subroutine print_spinor_8(S,MF)
1805	implicit none
1806	TYPE(spinor_8), INTENT(INOUT) :: s
1807	INTEGER MF,I
1808
1809	do i=1,3
1810	write(mf,*) ' Spin Variable ',i
1811	call print(s%x(i),mf)
1812	enddo
1813
1814	END subroutine print_spinor_8
1815
1816	subroutine print_res_spinor_8(S,MF)
1817	implicit none
1818	TYPE(res_spinor_8), INTENT(INOUT) :: s
1819	INTEGER MF,I
1820
1821	do i=1,3
1822	write(mf,*) ' Spin Variable ',i
1823	call print(s%x(i),mf)
1824	enddo
1825
1826	END subroutine print_res_spinor_8
1827
1828	subroutine read_spinor_8(S,MF)
1829	implicit none
1830	TYPE(spinor_8), INTENT(INOUT) :: s
1831	INTEGER MF,I
1832	character*255 line
1833	type(taylor) t
1834	call alloc(t)
1835
1836	do i=1,3
1837	read(mf,'(a255)') line
1838	call read(t,mf)
1839	s%x(i)=morph(t)
1840	enddo
1841
1842	call kill(t)
1843
1844	END subroutine read_spinor_8
1845
1846	subroutine READ_DASPIN(DS,MF,file)
1847	implicit none
1848	TYPE(damapspin), INTENT(INOUT) :: DS
1849	INTEGER MF1,I,J,nd2
1850	CHARACTER*20 LINE
1851	CHARACTER(*), optional :: file
1852	INTEGER, optional :: MF
1853	TYPE(TAYLOR) T
1854
1855	if(present(mf)) mf1=mf
1856	if(present(file)) then
1857	call kanalnummer(mf1)
1858	open(unit=mf1,file=file)
1859	endif
1860	CALL ALLOC(T)
1861
1862	! READ(MF,*) LINE
1863	! READ(MF,*) DS%X(1:3)
1864	! READ(MF,*) DS%X(4:6)
1865	READ(MF1,*) nd2, LINE(1:7)
1866	if(nd2>=c_%nd2) then
1867	CALL READ(DS%M,MF1)
1868	do i=c_%nd2+1,nd2
1869	call rea(t,mf1)
1870	enddo
1871	else
1872	do i=1,nd2
1873	call read(ds%m%v(i),mf1)
1874	enddo
1875	do i=nd2+1,c_%nd2
1876	ds%m%v(i)=1.0_dp.mono.i
1877	enddo
1878	endif
1879	READ(MF1,*) LINE
1880	DO I=1,3
1881	DO J=1,3
1882	READ(MF1,*) LINE
1883	CALL READ(T,MF1)
1884	DS%S%s(I,J)=T !MORPH(T)
1885	ENDDO
1886	ENDDO
1887	READ(MF1,*) LINE
1888	DO I=1,3
1889	DO J=1,3
1890	READ(MF1,*) LINE
1891	read(mf,*) ds%e_ij(i,j)
1892	ENDDO
1893	ENDDO
1894	if(present(file)) close(mf1)
1895
1896	CALL KILL(T)
1897	END subroutine READ_DASPIN
1898
1899	subroutine read_probe8(DS,MF)
1900	implicit none
1901	TYPE(probe_8), INTENT(INOUT) :: DS
1902	INTEGER MF,I
1903	character*120 line
1904	type(taylor) t
1905	call alloc(t)
1906
1907	read(mf,*) line
1908	do i=1,6
1909	read(mf,*) line
1910	call read(t,mf)
1911	ds%x(i)=morph(t)
1912	enddo
1913	do i=ISPIN0R,ISPIN1R
1914	read(mf,'(a120)') line
1915	call read(ds%s(i),mf)
1916	enddo
1917
1918	call kill(t)
1919
1920	END subroutine read_probe8
1921
1922
1923
1924	!!! NOT USED NOW
1925	subroutine CHECK_RES_ORBIT(J,NRES,M,SKIP)
1926	implicit none
1927	INTEGER M(:,:),NRES
1928	LOGICAL(lp) SKIP,SKIP1,SKIP2
1929	INTEGER I,K,J(:)
1930
1931	SKIP=.FALSE.
1932	IF(NRES==0) RETURN
1933
1934	DO I=1,NRES
1935	SKIP1=.TRUE.
1936	SKIP2=.TRUE.
1937	DO K=1,C_%ND
1938	SKIP1=((J(2K)-J(2K-1))==M(k,i)).AND.SKIP1
1939	SKIP2=((J(2K)-J(2K-1))==-M(k,i)).AND.SKIP2
1940	ENDDO
1941	IF(SKIP1.OR.SKIP2) THEN
1942	SKIP=.TRUE.
1943	RETURN
1944	ENDIF
1945	ENDDO
1946
1947	END subroutine CHECK_RES_ORBIT
1948
1949
1950
1951	subroutine find_n_thetar(s0,n0)
1952	implicit none
1953	real(dp),intent(in) :: s0(3,3)
1954	real(dp) theta0,n0(3)
1955	real(dp) det,ss(3,3),detm
1956	integer i,is,j
1957
1958
1959
1960
1961
1962	do i=1,3
1963	do j=1,3
1964	ss(i,j)=s0(i,j)
1965	enddo
1966	enddo
1967
1968
1969
1970	do i=1,3
1971	ss(i,i)=ss(i,i)-1.0_dp
1972	enddo
1973	det=(ss(2,2)ss(3,3)-ss(2,3)ss(3,2))
1974	is=1
1975
1976	detm=(ss(1,1)ss(3,3)-ss(1,3)ss(3,1))
1977	if(abs(detm)>=abs(det)) then
1978	det=detm
1979	is=2
1980	endif
1981	detm=(ss(1,1)ss(2,2)-ss(1,2)ss(2,1))
1982	if(abs(detm)>=abs(det)) then
1983	det=detm
1984	is=3
1985	endif
1986
1987	n0(is)=1.0_dp
1988	if(is==1) then
1989	n0(2)=(-ss(3,3)ss(2,1)+ss(2,3)ss(3,1))/det
1990	n0(3)=(-ss(2,2)ss(3,1)+ss(2,1)ss(3,2))/det
1991	elseif(is==2) then
1992	n0(1)=(-ss(3,3)ss(1,2)+ss(3,2)ss(1,3))/det
1993	n0(3)=(-ss(1,1)ss(3,2)+ss(1,2)ss(3,1))/det
1994	else
1995	n0(1)=(-ss(2,2)ss(1,3)+ss(2,3)ss(1,2))/det
1996	n0(2)=(-ss(1,1)ss(2,3)+ss(1,3)ss(2,1))/det
1997	endif
1998
1999
2000	theta0=sqrt(n0(1)2+n0(2)2+n0(3)**2)
2001
2002	do i=1,3
2003	n0(i)=n0(i)/theta0
2004	enddo
2005
2006
2007
2008
2009	end subroutine find_n_thetar
2010
2011	subroutine find_n_thetap(s0,n0)
2012	implicit none
2013	type(real_8),intent(in) :: s0(3,3)
2014
2015	type(real_8) theta0,n0(3)
2016	type(real_8) det,ss(3,3),detm
2017	integer i,is,j
2018
2019
2020	call alloc(det)
2021	call alloc(detm)
2022	call alloc(theta0)
2023	call alloc_33(ss)
2024
2025
2026	do i=1,3
2027	do j=1,3
2028	ss(i,j)=s0(i,j)
2029	enddo
2030	enddo
2031
2032
2033
2034	do i=1,3
2035	ss(i,i)=ss(i,i)-1.0_dp
2036	enddo
2037
2038	det=(ss(2,2)ss(3,3)-ss(2,3)ss(3,2))
2039	is=1
2040	detm=(ss(1,1)ss(3,3)-ss(1,3)ss(3,1))
2041	if(abs(detm)>=abs(det)) then
2042	det=detm
2043	is=2
2044	endif
2045	detm=ss(1,1)ss(2,2)-ss(1,2)ss(2,1)
2046	if(abs(detm)>=abs(det)) then
2047	det=detm
2048	is=3
2049	endif
2050
2051
2052
2053	n0(is)=1.0_dp
2054	if(is==1) then
2055	n0(2)=(-ss(3,3)ss(2,1)+ss(2,3)ss(3,1))/det
2056	n0(3)=(-ss(2,2)ss(3,1)+ss(2,1)ss(3,2))/det
2057	elseif(is==2) then
2058	n0(1)=(-ss(3,3)ss(1,2)+ss(3,2)ss(1,3))/det
2059	n0(3)=(-ss(1,1)ss(3,2)+ss(1,2)ss(3,1))/det
2060	else
2061	n0(1)=(-ss(2,2)ss(1,3)+ss(2,3)ss(1,2))/det
2062	n0(2)=(-ss(1,1)ss(2,3)+ss(1,3)ss(2,1))/det
2063	endif
2064
2065
2066	theta0=sqrt(n0(1)2+n0(2)2+n0(3)**2)
2067
2068
2069	do i=1,3
2070	n0(i)=n0(i)/theta0
2071	enddo
2072
2073
2074	call kill(det)
2075	call kill(detm)
2076	call kill(theta0)
2077	call kill_33(ss)
2078
2079	end subroutine find_n_thetap
2080
2081	! find exponent of rotation routines
2082	subroutine find_axisp(ds,H_axis) !
2083	implicit none
2084	type(damapspin),intent(in) :: ds
2085	type(spinor_8),intent(inout) :: H_axis
2086
2087	type(real_8) theta0,n0(3)
2088	type(real_8) det,ss(3,3),detm
2089	integer i,is,j
2090
2091
2092	call alloc(det)
2093	call alloc(detm)
2094	call alloc(theta0)
2095	call alloc_33(ss)
2096	call alloc(n0)
2097
2098
2099	do i=1,3
2100	do j=1,3
2101	ss(i,j)=ds%s%s(i,j)
2102	enddo
2103	enddo
2104
2105
2106
2107	do i=1,3
2108	ss(i,i)=ss(i,i)-1.0_dp
2109	enddo
2110
2111	det=(ss(2,2)ss(3,3)-ss(2,3)ss(3,2))
2112	is=1
2113	detm=(ss(1,1)ss(3,3)-ss(1,3)ss(3,1))
2114	if(abs(detm)>=abs(det)) then
2115	det=detm
2116	is=2
2117	endif
2118	detm=ss(1,1)ss(2,2)-ss(1,2)ss(2,1)
2119	if(abs(detm)>=abs(det)) then
2120	det=detm
2121	is=3
2122	endif
2123
2124
2125
2126	n0(is)=1.0_dp
2127	if(is==1) then
2128	n0(2)=(-ss(3,3)ss(2,1)+ss(2,3)ss(3,1))/det
2129	n0(3)=(-ss(2,2)ss(3,1)+ss(2,1)ss(3,2))/det
2130	elseif(is==2) then
2131	n0(1)=(-ss(3,3)ss(1,2)+ss(3,2)ss(1,3))/det
2132	n0(3)=(-ss(1,1)ss(3,2)+ss(1,2)ss(3,1))/det
2133	else
2134	n0(1)=(-ss(2,2)ss(1,3)+ss(2,3)ss(1,2))/det
2135	n0(2)=(-ss(1,1)ss(2,3)+ss(1,3)ss(2,1))/det
2136	endif
2137
2138	theta0=sqrt(n0(1)2+n0(2)2+n0(3)**2)
2139
2140	do i=1,3
2141	n0(i)=n0(i)/theta0
2142	enddo
2143
2144	h_axis%x(1)=n0(1); h_axis%x(2)=n0(2) ; h_axis%x(3)=n0(3);
2145
2146
2147
2148
2149
2150	call kill(det)
2151	call kill(detm)
2152	call kill(theta0)
2153	call kill_33(ss)
2154	call kill(n0)
2155
2156	end subroutine find_axisp
2157
2158	subroutine find_perp_basisp(y_axis,x_axis,z_axis) !
2159	implicit none
2160	type(spinor_8),intent(inout) :: y_axis,x_axis,z_axis
2161	integer i,is
2162	type(real_8) norm
2163
2164	call alloc(norm)
2165
2166	is=1
2167	if(abs(y_axis%x(is))>abs(y_axis%x(2))) then
2168	is=2
2169	endif
2170	if(abs(y_axis%x(is))>abs(y_axis%x(3))) then
2171	is=3
2172	endif
2173
2174	! now is = smallest
2175
2176	x_axis=is
2177
2178	x_axis=x_axis-(x_axis.dot.y_axis)*y_axis
2179
2180	norm=sqrt(x_axis.dot.x_axis)
2181	x_axis=(1.d0/norm)*x_axis
2182
2183
2184
2185	z_axis=x_axis*y_axis
2186	norm=sqrt(z_axis.dot.z_axis)
2187	z_axis=(1.d0/norm)*z_axis
2188
2189
2190
2191	call kill(norm)
2192
2193	end subroutine find_perp_basisp
2194
2195	subroutine find_exponent_jet_p(ds,h_axis) !
2196	implicit none
2197	type(damapspin),intent(inout) :: ds
2198	type(spinor_8),intent(inout) :: h_axis
2199	type(damapspin) s,sf,st,SA
2200	integer i,nmax
2201	real(dp) eps1,EPS0
2202	LOGICAL DONOT
2203	nmax=1000
2204	call alloc(s,sf,st,SA)
2205
2206	DONOT=.FALSE.
2207	EPS0=MYBIG
2208
2209	st=1
2210	s=ds
2211	s%m=1
2212	sf=0
2213	do i=1,3
2214	s%s%s(i,i)=s%s%s(i,i)-1.0_dp
2215	enddo
2216
2217	st=s
2218	do i=1,nmax
2219	SA=SF
2220	sf=sf + (1.0_dp/i)*st
2221	st= s*st
2222	st=(-1.0_dp)*st
2223	SA=SA+((-1.0_dp)*SF)
2224	eps1=FULL_ABS(SA)
2225	IF(EPS1>deps_tracking) THEN
2226	EPS0=EPS1
2227	ELSE
2228	IF(EPS1>=EPS0.or.EPS1<=PUNY) THEN
2229	DONOT=.TRUE.
2230	ELSE
2231	EPS0=EPS1
2232	ENDIF
2233	ENDIF
2234
2235	IF(DONOT) EXIT
2236
2237	enddo
2238
2239	if(i>nmax-2) then
2240	write(6,*) "Did not converge in find_exponent_jet_p"
2241	stop 666
2242	endif
2243
2244	h_axis%x(1)=sf%s%s(3,2)
2245	h_axis%x(2)=sf%s%s(1,3)
2246	h_axis%x(3)=sf%s%s(2,1)
2247
2248	!etienne
2249
2250
2251	call kill(s,sf,st,SA)
2252	end subroutine find_exponent_jet_p
2253
2254	subroutine find_exponentp(ds,y_axis,x_axis,z_axis,h_axis,spin_tune) !
2255	implicit none
2256	type(damapspin),intent(inout) :: ds
2257	type(spinor_8),intent(inout) :: y_axis,x_axis,z_axis,h_axis
2258	type(spinor_8) s
2259	type(real_8) cos,sin,theta
2260	type(real_8) , optional :: spin_tune
2261
2262	call alloc(s)
2263	call alloc(cos,sin,theta)
2264
2265	call find_axis(ds,y_axis)
2266	! etienne
2267
2268	call find_perp_basis(y_axis,x_axis,z_axis)
2269
2270
2271
2272	s=ds*x_axis
2273
2274
2275	cos=s.dot.x_axis
2276	sin=-(s.dot.z_axis)
2277
2278	theta=clockwise*atan2(sin,cos)
2279	!if(force_positive.and.theta<zero) theta = theta + twopi !!!! allow negative theta
2280
2281	h_axis=(clockwisetheta)y_axis
2282
2283	if(present(spin_tune)) then
2284	spin_tune=(1.0_dp/twopi)*theta
2285	endif
2286
2287	call kill(cos,sin,theta)
2288	call kill(s)
2289	end subroutine find_exponentp
2290
2291	subroutine find_exponent_only(ds,h_axis,spin_tune) !
2292	implicit none
2293	type(damapspin),intent(inout) :: ds
2294	type(spinor_8),intent(inout) :: h_axis
2295	type(spinor_8) y_axis,x_axis,z_axis
2296	type(spinor_8) s
2297	type(real_8) cos,sin,theta
2298	type(real_8) , optional :: spin_tune
2299
2300	call alloc(s)
2301	call alloc(cos,sin,theta)
2302	call alloc(y_axis)
2303	call alloc(x_axis)
2304	call alloc(z_axis)
2305	call find_axis(ds,y_axis)
2306	! etienne
2307
2308	call find_perp_basis(y_axis,x_axis,z_axis)
2309
2310
2311
2312	s=ds*x_axis
2313
2314
2315	cos=s.dot.x_axis
2316	sin=-(s.dot.z_axis)
2317
2318	theta=clockwise*atan2(sin,cos)
2319	!if(force_positive.and.theta<zero) theta = theta + twopi !!!! allow negative theta
2320
2321	h_axis=(clockwisetheta)y_axis
2322
2323	if(present(spin_tune)) then
2324	spin_tune=(1.0_dp/twopi)*theta
2325	endif
2326
2327	call kill(y_axis)
2328	call kill(x_axis)
2329	call kill(z_axis)
2330	call kill(cos,sin,theta)
2331	call kill(s)
2332	end subroutine find_exponent_only
2333
2334	subroutine remove_y_rot(as_xyz,r_y)
2335	implicit none
2336	type(damapspin), intent(inout) :: as_xyz
2337	type(damapspin), optional , intent(inout) ::r_y
2338	type(damapspin) temp,as_y,as_nl,rot_y
2339	type(spinor_8) n_expo,n_tune
2340	type(taylorresonance) tr
2341	type(taylor) t
2342	integer i,j
2343	integer nmax
2344	real(dp) c,eps,norm1,norm2
2345	logical check
2346	!!! original as_xyz = as_xyzr_y = a_ya_nl*r_y on exit
2347	check=.true.
2348	eps=1.d-1
2349	nmax=1000
2350
2351	call alloc(n_expo)
2352	call alloc(n_tune)
2353	call alloc(temp,as_y,as_nl,rot_y)
2354	call alloc(tr)
2355	call alloc(t)
2356
2357	as_y=as_xyz
2358
2359	norm1=mybig
2360	rot_y=1
2361
2362	check=.true.
2363	norm1=mybig
2364	do i=1,nmax
2365	call find_exponent_only(as_y,n_expo)
2366	! call dalog_spinor_8(as_y,n_expo)
2367	norm2=0.0_dp
2368	n_tune%x(1)=0.0_dp
2369	n_tune%x(3)=0.0_dp
2370	! do j=1,3
2371	j=2 ! do loop was a theory error
2372	t=n_expo%x(j)
2373	tr=t
2374	call cfu(tr%cos,phase_shift,tr%cos)
2375	tr%sin=0.0_dp
2376	t=tr
2377	n_tune%x(j)=t
2378	! enddo
2379	temp=exp(n_tune)
2380
2381	rot_y=temp*rot_y
2382	call inv_as(temp%s%s,temp%s%s)
2383	as_y=as_y*temp
2384	if(check) then
2385	if(norm2<eps) then
2386	check=.false.
2387	endif
2388	else
2389	if(norm2>=norm1) exit
2390	endif
2391	norm1=norm2
2392
2393	enddo
2394	if(i>nmax-10) then
2395	write(6,*) "no convergence in remove_y_rot "
2396	stop 1067
2397	endif
2398	as_xyz=as_y
2399	if(present(r_y)) r_y=rot_y
2400
2401	call kill(n_expo)
2402	call kill(temp,as_y,as_nl,rot_y)
2403	call kill(tr)
2404	call kill(n_tune)
2405	call kill(t)
2406
2407	end subroutine remove_y_rot
2408
2409	subroutine remove_y_rot0(as_xyz,a_y,a_nl,r_y)
2410	implicit none
2411	type(damapspin), intent(inout) :: as_xyz
2412	type(damapspin), optional , intent(inout) :: a_y,a_nl,r_y
2413	type(damapspin) temp,as_y,as_nl,rot_y
2414	type(spinor_8) n_expo,n_tune
2415	type(taylorresonance) tr
2416	type(taylor) t
2417	integer i,j
2418	integer nmax
2419	real(dp) c,eps,norm1,norm2
2420	logical check
2421	!!! original as_xyz = as_xyzr_y = a_ya_nl*r_y on exit
2422	check=.true.
2423	eps=1.d-5
2424	nmax=1000
2425
2426	call alloc(n_expo)
2427	call alloc(n_tune)
2428	call alloc(temp,as_y,as_nl,rot_y)
2429	call alloc(tr)
2430	call alloc(t)
2431
2432	call factor_parameter_dependent_s0(as_xyz,as_y,as_nl,n_expo,1)
2433
2434	norm1=mybig
2435	rot_y=1
2436	do i=1,nmax
2437	call find_exponent_only(as_y,n_expo)
2438	n_expo%x(1)=0.0_dp
2439	n_expo%x(3)=0.0_dp
2440	temp=exp(n_expo)
2441	rot_y=temp*rot_y
2442	call inv_as(temp%s%s,temp%s%s)
2443	as_y=as_y*temp
2444	call norm_spinor_8(n_expo,norm2)
2445	if(check) then
2446	if(norm2<eps) then
2447	check=.false.
2448	endif
2449	else
2450	if(norm2>=norm1) exit
2451	endif
2452	norm1=norm2
2453	enddo
2454	if(i>nmax-10) then
2455	write(6,*) "no convergence in daexplogp "
2456	stop 1066
2457	endif
2458	as_nl=rot_yas_nlrot_y**(-1)
2459
2460	check=.true.
2461	norm1=mybig
2462	do i=1,nmax
2463	call dalog_spinor_8(as_nl,n_expo)
2464	norm2=0.0_dp
2465	n_tune%x(1)=0.0_dp
2466	n_tune%x(3)=0.0_dp
2467	! do j=1,3
2468	j=2 ! do loop was a theory error
2469	t=n_expo%x(j)
2470	tr=t
2471	call cfu(tr%cos,phase_shift,tr%cos)
2472	tr%sin=0.0_dp
2473	t=tr
2474	n_tune%x(j)=t
2475	norm2=norm2+full_abs(t)
2476	! enddo
2477	temp=exp(n_tune)
2478
2479	rot_y=temp*rot_y
2480	call inv_as(temp%s%s,temp%s%s)
2481	as_nl=as_nl*temp
2482	if(check) then
2483	if(norm2<eps) then
2484	check=.false.
2485	endif
2486	else
2487	if(norm2>=norm1) exit
2488	endif
2489	norm1=norm2
2490
2491	enddo
2492
2493	if(i>nmax-10) then
2494	write(6,*) "no convergence in daexplogp "
2495	stop 1067
2496	endif
2497	as_xyz=as_y*as_nl
2498	if(present(r_y)) r_y=rot_y
2499	if(present(a_y)) a_y=as_y
2500	if(present(a_nl)) a_nl=as_nl
2501
2502	call kill(n_expo)
2503	call kill(temp,as_y,as_nl,rot_y)
2504	call kill(tr)
2505	call kill(n_tune)
2506	call kill(t)
2507
2508	end subroutine remove_y_rot0
2509
2510	subroutine daexplogp(h_axis,ds)
2511	implicit none
2512	TYPE(damapspin), INTENT(INout) :: DS
2513	TYPE(spinor_8), INTENT(IN) :: h_axis
2514	integer nmax
2515	integer i
2516	TYPE(damapspin) dh,dhn,dr
2517	real(dp) c,eps,norm1,norm2
2518	logical check
2519
2520	check=.true.
2521	eps=1.d-5
2522	nmax=1000
2523
2524	call alloc(dh)
2525	call alloc(dhn)
2526	call alloc(dr)
2527
2528	! this only works with a da-map
2529	ds=1
2530	dh%m=1
2531	dh%s%s(2,1)=h_axis%x(3)
2532	dh%s%s(1,3)=h_axis%x(2)
2533	dh%s%s(3,2)=h_axis%x(1)
2534	dh%s%s(1,2)=-h_axis%x(3)
2535	dh%s%s(3,1)=-h_axis%x(2)
2536	dh%s%s(2,3)=-h_axis%x(1)
2537	dhn=1
2538	c=1.0_dp
2539	norm1=mybig
2540	do i=1,nmax
2541	dhn=dhn*dh
2542	c=c/i
2543	dr=ds
2544	ds=ds+c*dhn
2545	dr=ds+(-1.0_dp)*dr
2546	call norm_damapspin(dr,norm2)
2547	if(check) then
2548	if(norm2<eps) then
2549	check=.false.
2550	endif
2551	else
2552	if(norm2>=norm1) exit
2553	endif
2554	norm1=norm2
2555	enddo
2556
2557	if(i>nmax-10) then
2558	write(6,*) "no convergence in daexplogp "
2559	stop 1066
2560	endif
2561
2562	call kill(dh)
2563	call kill(dhn)
2564	call kill(dr)
2565
2566	end subroutine daexplogp
2567
2568	subroutine norm_damapspin(ds,norm)
2569	implicit none
2570	TYPE(damapspin), INTENT(INout) :: DS
2571	real(dp) norm
2572	integer i,j
2573
2574	norm=0.0_dp
2575
2576	do i=1,3
2577	do j=1,3
2578	norm=norm+full_abs( ds%s%s(i,j) )
2579	enddo
2580	enddo
2581
2582	end subroutine norm_damapspin
2583
2584	subroutine norm_spinor_8(s,norm)
2585	implicit none
2586	TYPE(spinor_8), INTENT(INout) :: s
2587	real(dp) norm
2588	integer i
2589
2590	norm=0.0_dp
2591
2592	do i=1,3
2593	norm=norm+full_abs( s%x(i) )
2594	enddo
2595
2596	end subroutine norm_spinor_8
2597
2598	! end of find exponent of rotation routines
2599
2600
2601	subroutine find_ar(n2,a)
2602	implicit none
2603	real(dp) n2(3),n1(3),n3(3)
2604	real(dp) a(3,3),s,n
2605	integer i,is
2606
2607	! here we find smallest value of n2
2608	is=2
2609	if(abs(n2(1))< abs(n2(2))) is=1
2610
2611	if(is==1) then
2612	if(abs(n2(3))<abs(n2(1))) is=3
2613	else
2614	if(abs(n2(3))<abs(n2(2))) is=3
2615	endif
2616
2617	! put n1 in along that value
2618	n1=0.0_dp
2619	n1(is)=1.0_dp
2620
2621	s=n2(is)*n1(is)
2622
2623	n=0.0_dp
2624	do i=1,3
2625	n1(i)=n1(i)-s*n2(i)
2626	n=n1(i)**2+n
2627	enddo
2628	n1=n1/sqrt(n)
2629
2630	n3(1)=n1(2)n2(3)-n1(3)n2(2)
2631	n3(2)=n1(3)n2(1)-n1(1)n2(3)
2632	n3(3)=n1(1)n2(2)-n1(2)n2(1)
2633
2634	n=0.0_dp
2635	do i=1,3
2636	n=n3(i)**2+n
2637	enddo
2638	n3=n3/sqrt(n)
2639	! spin_normal_position
2640
2641	! a=zero
2642	if(spin_normal_position==2) then
2643	a(:,1)=n1
2644	a(:,2)=n2
2645	a(:,3)=n3
2646	elseif(spin_normal_position==3) then
2647	a(:,2)=n1
2648	a(:,3)=n2
2649	a(:,1)=n3
2650	else
2651	a(:,3)=n1
2652	a(:,1)=n2
2653	a(:,2)=n3
2654	endif
2655
2656
2657	end subroutine find_ar
2658
2659	subroutine find_ap(n2,a)
2660	implicit none
2661	type(real_8) n2(3),n1(3),n3(3)
2662	type(real_8) a(3,3),s,n
2663	real(dp) aa(3,3)
2664	integer j
2665	integer i,is
2666	call alloc(n1,3)
2667	call alloc(n3,3)
2668	call alloc(s,n)
2669	! here we find smallest value of n2
2670	is=2
2671	if(abs(n2(1))< abs(n2(2))) is=1
2672
2673	if(is==1) then
2674	if(abs(n2(3))<abs(n2(1))) is=3
2675	else
2676	if(abs(n2(3))<abs(n2(2))) is=3
2677	endif
2678
2679	! put n1 in along that value
2680	do i=1,3
2681	n1(i)=0.0_dp
2682	enddo
2683	n1(is)=1.0_dp
2684
2685	s=n2(is)*n1(is)
2686
2687	n=0.0_dp
2688	do i=1,3
2689	n1(i)=n1(i)-s*n2(i)
2690	n=n1(i)**2+n
2691	enddo
2692	do i=1,3
2693	n1(i)=n1(i)/sqrt(n)
2694	enddo
2695
2696	n3(1)=n1(2)n2(3)-n1(3)n2(2)
2697	n3(2)=n1(3)n2(1)-n1(1)n2(3)
2698	n3(3)=n1(1)n2(2)-n1(2)n2(1)
2699
2700	n=0.0_dp
2701	do i=1,3
2702	n=n3(i)**2+n
2703	enddo
2704	do i=1,3
2705	n3(i)=n3(i)/sqrt(n)
2706	enddo
2707
2708	if(spin_normal_position==2) then
2709	do i=1,3
2710	a(i,1)=n1(i)
2711	a(i,2)=n2(i)
2712	a(i,3)=n3(i)
2713	enddo
2714	elseif(spin_normal_position==3) then
2715	do i=1,3
2716	a(i,2)=n1(i)
2717	a(i,3)=n2(i)
2718	a(i,1)=n3(i)
2719	enddo
2720	else
2721	do i=1,3
2722	a(i,3)=n1(i)
2723	a(i,1)=n2(i)
2724	a(i,2)=n3(i)
2725	enddo
2726	endif
2727
2728
2729	call kill(n1,3)
2730	call kill(n3,3)
2731	call kill(s,n)
2732	end subroutine find_ap
2733
2734
2735
2736	subroutine copy_damap_matrix(mi,a)
2737	implicit none
2738	type(taylor), intent(inout) :: a(:,:)
2739	type(damap), intent(in) :: mi
2740	type(damap) m
2741
2742	integer i,j,nt
2743	logical doflip
2744	integer, allocatable :: jl(:)
2745
2746	call alloc(m)
2747
2748	m=mi
2749
2750	if(perform_flip.and.new_ndpt.and.c_%ndpt/=0) then
2751	perform_flip=.false.
2752	call flip_damap(m,m)
2753	doflip=.true.
2754	else
2755	doflip=.false.
2756	endif
2757	nt=c_%nd2
2758	if(c_%ndpt/=0)nt=nt-2
2759
2760	allocate(jl(nt))
2761	jl=0
2762	do i=1,min(nt,size(a,dim=1))
2763	do j=1,min(nt,size(a,dim=2))
2764	jl(j)=1
2765	a(i,j)=m%v(i).par.jl
2766	jl(j)=0
2767	enddo
2768	enddo
2769
2770	if(doflip) then
2771	call flip_damap(m,m)
2772	perform_flip=.true.
2773	endif
2774
2775	call kill(m)
2776	deallocate(jl)
2777
2778	end subroutine copy_damap_matrix
2779
2780	subroutine copy_matrix_matrix(ma,a)
2781	implicit none
2782	type(taylor), intent(inout) :: a(:,:)
2783	type(taylor), intent(in) :: ma(:,:)
2784	integer i,j
2785
2786	do i=1,size(a,dim=1)
2787	do j=1,size(a,dim=2)
2788	a(i,j)=ma(i,j)
2789	enddo
2790	enddo
2791
2792	end subroutine copy_matrix_matrix
2793
2794
2795
2796	subroutine alloc_33t(a)
2797	implicit none
2798	type(taylor) a(:,:)
2799	integer i,j
2800
2801	do i=1,size(a,dim=1)
2802	do j=1,size(a,dim=2)
2803	call alloc(a(i,j))
2804	enddo
2805	enddo
2806
2807	end subroutine alloc_33t
2808
2809	subroutine kill_33t(a)
2810	implicit none
2811	type(taylor) a(:,:)
2812	integer i,j
2813
2814	do i=1,size(a,dim=1)
2815	do j=1,size(a,dim=2)
2816	call kill(a(i,j))
2817	enddo
2818	enddo
2819
2820	end subroutine kill_33t
2821
2822	subroutine zero_33t(a,r)
2823	implicit none
2824	type(taylor) a(:,:) ! taylor
2825	! type(taylor) a(3,3)
2826	integer i,j
2827	real(dp), optional :: r
2828
2829	do i=1,size(a,dim=1)
2830	do j=1,size(a,dim=2)
2831	a(i,j)=0.0_dp
2832	enddo
2833	if(present(r)) a(i,i)=1.0_dp
2834	enddo
2835	end subroutine zero_33t
2836
2837	subroutine zero_33p(a,r)
2838	implicit none
2839	type(real_8) a(:,:) ! taylor
2840	! type(taylor) a(3,3)
2841	integer i,j
2842	real(dp), optional :: r
2843
2844	do i=1,size(a,dim=1)
2845	do j=1,size(a,dim=2)
2846	a(i,j)=0.0_dp
2847	enddo
2848	if(present(r)) a(i,i)=1.0_dp
2849	enddo
2850
2851	end subroutine zero_33p
2852
2853	subroutine alloc_33p(a)
2854	implicit none
2855	type(real_8) a(:,:) ! taylor
2856	integer i,j
2857
2858	do i=1,size(a,dim=1)
2859	do j=1,size(a,dim=2)
2860	call alloc(a(i,j))
2861	enddo
2862	enddo
2863
2864	end subroutine alloc_33p
2865
2866	subroutine kill_33p(a)
2867	implicit none
2868	type(real_8) a(:,:)
2869	integer i,j
2870
2871	do i=1,size(a,dim=1)
2872	do j=1,size(a,dim=2)
2873	call kill(a(i,j))
2874	enddo
2875	enddo
2876
2877	end subroutine kill_33p
2878
2879	subroutine matmul_33(m,n,mo,sc)
2880	implicit none
2881	type(taylor) m(:,:),n(:,:),mo(:,:)
2882	type(taylor), allocatable :: a(:,:)
2883	real(dp), optional :: sc
2884	real(dp) sc0
2885	integer i,j,k
2886	sc0=1.0_dp
2887	allocate(a(size(m,dim=1),size(n,dim=2)))
2888
2889	call alloc_33(a)
2890
2891	do i=1,size(m,dim=1)
2892	do j=1,size(m,dim=2)
2893	do k=1,size(n,dim=2)
2894	a(i,k)=m(i,j)*n(j,k)+a(i,k)
2895	enddo
2896	enddo
2897	enddo
2898
2899	if(present(sc)) sc0=sc
2900	do i=1,size(mo,dim=1)
2901	do j=1,size(mo,dim=2)
2902	mo(i,j)=sc0*a(i,j)
2903	enddo
2904	enddo
2905	call kill_33(a)
2906	deallocate(a)
2907	end subroutine matmul_33
2908
2909	subroutine matmul_3344(m,n,mo)
2910	implicit none
2911	type(taylor) m(:,:),n(:,:),mo(:,:)
2912	type(taylor), allocatable :: a(:,:)
2913	integer i,j,k
2914
2915	allocate(a(size(m,dim=1),size(n,dim=2)))
2916
2917	call alloc_33(a)
2918
2919	do i=1,size(m,dim=1)
2920	do j=1,size(m,dim=2)
2921	do k=1,size(n,dim=2)
2922	a(i,k)=m(i,j)*n(j,k)+a(i,k)
2923	enddo
2924	enddo
2925	enddo
2926
2927
2928	do i=1,size(mo,dim=1)
2929	do j=1,size(mo,dim=2)
2930	mo(i,j)=a(i,j)
2931	enddo
2932	enddo
2933	call kill_33(a)
2934	deallocate(a)
2935	end subroutine matmul_3344
2936
2937	subroutine matmulp(m,n,mo)
2938	implicit none
2939	type(real_8) m(3,3),n(3,3),mo(3,3),a(3,3)
2940	integer i,j,k
2941
2942
2943	call alloc_33(a)
2944
2945	do i=1,3
2946	do j=1,3
2947	do k=1,3
2948	a(i,k)=m(i,j)*n(j,k)+a(i,k)
2949	enddo
2950	enddo
2951	enddo
2952
2953
2954	do i=1,3
2955	do j=1,3
2956	mo(i,j)=a(i,j)
2957	enddo
2958	enddo
2959	call kill_33(a)
2960
2961	end subroutine matmulp
2962
2963	subroutine smatmulp(sc,m,n,mo)
2964	implicit none
2965	type(real_8) m(3,3),n(3,3),mo(3,3),a(3,3)
2966	integer i,j
2967	real(dp) sc
2968
2969	call alloc_33(a)
2970
2971	do i=1,3
2972	do j=1,3
2973	a(i,j)=m(i,j)*sc+n(i,j)
2974	enddo
2975	enddo
2976
2977
2978	do i=1,3
2979	do j=1,3
2980	mo(i,j)=a(i,j)
2981	enddo
2982	enddo
2983	call kill_33(a)
2984
2985	end subroutine smatmulp
2986
2987	subroutine smatp(sc,m,mo)
2988	implicit none
2989	type(real_8) m(3,3),mo(3,3),a(3,3)
2990	integer i,j
2991	real(dp) sc
2992
2993	call alloc_33(a)
2994
2995	do i=1,3
2996	do j=1,3
2997	a(i,j)=m(i,j)*sc
2998	enddo
2999	enddo
3000
3001
3002	do i=1,3
3003	do j=1,3
3004	mo(i,j)=a(i,j)
3005	enddo
3006	enddo
3007	call kill_33(a)
3008
3009	end subroutine smatp
3010
3011
3012	subroutine make_unitary_p(m)
3013	! inverse 3x3 rrotation
3014	implicit none
3015	integer i, nmax
3016	type(real_8) m(3,3)
3017	type(real_8) mt(3,3),id(3,3)
3018	real(dp) r,eps,rb
3019	logical doit
3020
3021	doit=.true.
3022	nmax=100
3023	eps=1.e-4_dp
3024	call alloc_33(mt)
3025	call alloc_33(id)
3026
3027	do i=1,3
3028	id(i,i)=1.5e0_dp
3029	enddo
3030
3031	rb=mybig
3032
3033	do i=1,nmax
3034	call transpose_p(m,mt)
3035
3036	call matmulp(m,mt,mt)
3037
3038	call smatp(-0.5_dp,mt,mt)
3039
3040	call smatmulp(1.0_dp,mt,id,mt)
3041
3042	call matmulp(mt,m,m)
3043
3044	call transpose_p(m,mt)
3045	call matmulp(mt,m,mt)
3046	call smatmulp(-1.0_dp/1.5e0_dp,id,mt,mt)
3047	call absolute_p(mt,r)
3048	if(r<eps.and.doit) then
3049	doit=.false.
3050	rb=r
3051	else
3052	if(rb<=r) exit
3053	rb=r
3054	endif
3055	enddo
3056
3057	if(i>=nmax) then
3058	write(6,*) " Too many iterations "
3059	write(6,*) " Perhaps did not converged in make_unitary_p "
3060
3061	endif
3062
3063	call kill_33(id)
3064	call kill_33(mt)
3065
3066	end subroutine make_unitary_p
3067
3068	subroutine check_unitary_p(m,r,rt,cc)
3069	! inverse 3x3 rrotation
3070	implicit none
3071	integer i, j,c
3072	type(real_8) m(3,3)
3073	type(real_8) mt(3,3)
3074	real(dp) r,rt
3075	integer, optional :: cc
3076
3077	c=c_%no+1
3078	if(present(cc)) c=cc
3079
3080	call alloc_33(mt)
3081
3082	r=0.0_dp
3083	rt=0.0_dp
3084	call transpose_p(m,mt)
3085
3086	call matmulp(m,mt,mt)
3087
3088	do i=1,3
3089	do j=1,3
3090
3091	r=r+abs(mt(i,j))
3092	rt=rt+full_abs(mt(i,j).cut.c)
3093	enddo
3094	enddo
3095
3096
3097
3098	call kill_33(mt)
3099
3100	end subroutine check_unitary_p
3101
3102
3103	subroutine absolute_p(m,r)
3104	implicit none
3105	integer i,j
3106	type(real_8) m(3,3)
3107	real(dp) r
3108
3109	r=0.0_dp
3110
3111	do i=1,3
3112	do j=1,3
3113	r=r+full_abs(m(i,j))
3114	enddo
3115	enddo
3116
3117	end subroutine absolute_p
3118
3119	subroutine transpose_p(m,mi)
3120	! inverse 3x3 rrotation
3121	implicit none
3122	integer i,j
3123	type(real_8), INTENT(INOUT) :: m(3,3),mi(3,3)
3124	type(real_8) mt(3,3)
3125
3126	call alloc_33(mt)
3127
3128	call smatp(1.0_dp,m,mt)
3129
3130	do i=1,3
3131	do j=1,3
3132	mi(j,i)=mt(i,j)
3133	enddo
3134	enddo
3135
3136	call kill_33(mt)
3137
3138	end subroutine transpose_p
3139
3140	subroutine inv_asr(m,mi)
3141	! inverse 3x3 rrotation
3142	implicit none
3143	integer i,j
3144	real(dp) m(3,3),mi(3,3)
3145	real(dp) n(3,3)
3146
3147
3148	do i=1,3
3149	do j=1,3
3150	n(j,i)=m(i,j)
3151	enddo
3152	enddo
3153
3154	do i=1,3
3155	do j=1,3
3156	mi(i,j)=n(i,j)
3157	enddo
3158	enddo
3159
3160	end subroutine inv_asr
3161
3162
3163	subroutine inv_asp(m,mi)
3164	! inverse 3x3 rrotation
3165	implicit none
3166	integer i,j
3167	type(real_8) m(3,3),mi(3,3)
3168	type(real_8) n(3,3)
3169
3170	call alloc_33(n)
3171
3172	do i=1,3
3173	do j=1,3
3174	n(j,i)=m(i,j)
3175	enddo
3176	enddo
3177
3178	do i=1,3
3179	do j=1,3
3180	mi(i,j)=n(i,j)
3181	enddo
3182	enddo
3183
3184	call kill_33(n)
3185	end subroutine inv_asp
3186
3187	subroutine trans_mat(m,ri,mi)
3188	implicit none
3189	integer i,j
3190	type(real_8) m(3,3),mi(3,3)
3191	type(real_8) n(3,3)
3192	type(damap) ri
3193	type(real_8) r
3194
3195	call alloc_33(n)
3196	call alloc(r)
3197
3198	do i=1,3
3199	do j=1,3
3200	if(c_%nd2/=0) then
3201	r=(m(i,j)*ri)
3202	else
3203	r=m(i,j)
3204	endif
3205	n(i,j)=r !morph(r)
3206	! advances by r**(-1)
3207	enddo
3208	enddo
3209
3210
3211	do i=1,3
3212	do j=1,3
3213	mi(i,j)=n(i,j)
3214	enddo
3215	enddo
3216
3217	call kill(r)
3218	call kill_33(n)
3219	end subroutine trans_mat
3220
3221
3222
3223
3224
3225
3226	subroutine ALLOC_DASPIN(D)
3227	implicit none
3228	TYPE(damapspin), INTENT(INOUT) :: D
3229	INTEGER I,J
3230	! D%s0=ZERO
3231	CALL ALLOC(D%M)
3232	DO I=1,3
3233	DO J=1,3
3234	CALL ALLOC(D%S%s(I,J))
3235	d%e_ij(i,j)=0.0_dp
3236	ENDDO
3237	ENDDO
3238
3239	END subroutine ALLOC_DASPIN
3240
3241	subroutine KILL_DASPIN(D)
3242	implicit none
3243	TYPE(damapspin), INTENT(INOUT) :: D
3244	INTEGER I,J
3245
3246	! D%s0=ZERO
3247	CALL KILL(D%M)
3248	DO I=1,3
3249	DO J=1,3
3250	CALL KILL(D%S%s(I,J))
3251	d%e_ij(i,j)=0.0_dp
3252	ENDDO
3253	ENDDO
3254
3255	END subroutine KILL_DASPIN
3256
3257
3258
3259
3260	subroutine ALLOC_SPINOR_8(S)
3261	implicit none
3262	TYPE(SPINOR_8), INTENT(INOUT) :: S
3263
3264	CALL ALLOC(S%X,3)
3265	! S%G=A_PARTICLE
3266	END subroutine ALLOC_SPINOR_8
3267
3268	subroutine ALLOC_res_SPINOR_8(S)
3269	implicit none
3270	TYPE(RES_SPINOR_8), INTENT(INOUT) :: S
3271
3272	CALL ALLOC(S%X,3)
3273	! S%G=A_PARTICLE
3274	END subroutine ALLOC_res_SPINOR_8
3275
3276	subroutine ALLOC_probe_8(R)
3277	implicit none
3278	TYPE(probe_8), INTENT(INOUT) :: R
3279	INTEGER I !,J
3280
3281	! CALL ALLOC(R%S)
3282	DO I=1,3
3283	CALL ALLOC(R%S(I))
3284	ENDDO
3285	CALL ALLOC(R%X,6)
3286	! R%S(0)%X(N0_NORMAL)=ONE
3287	DO I=1,3
3288	R%S(I)=0
3289	ENDDO
3290	CALL ALLOC(R%ac)
3291	r%e_ij=0.0_dp
3292
3293	END subroutine ALLOC_probe_8
3294
3295	subroutine ALLOC_rf_phasor_8(R)
3296	implicit none
3297	TYPE(rf_phasor_8), INTENT(INOUT) :: R
3298	INTEGER I !,J
3299
3300	DO I=1,2
3301	CALL alloc(R%X(I))
3302	ENDDO
3303	CALL alloc(R%om)
3304	! CALL alloc(R%t)
3305
3306	END subroutine ALLOC_rf_phasor_8
3307
3308	subroutine KILL_res_SPINOR_8(S)
3309	implicit none
3310	TYPE(RES_SPINOR_8), INTENT(INOUT) :: S
3311
3312	CALL KILL(S%X,3)
3313	! S%G=A_PARTICLE
3314	END subroutine KILL_res_SPINOR_8
3315
3316	subroutine KILL_SPINOR_8(S)
3317	implicit none
3318	TYPE(SPINOR_8), INTENT(INOUT) :: S
3319
3320	CALL KILL(S%X,3)
3321
3322
3323	END subroutine KILL_SPINOR_8
3324
3325	subroutine KILL_probe_8(R)
3326	implicit none
3327	TYPE(probe_8), INTENT(INOUT) :: R
3328	INTEGER I !,J
3329
3330	! CALL KILL(R%S)
3331	DO I=1,3
3332	CALL KILL(R%S(I))
3333	ENDDO
3334	CALL KILL(R%X,6)
3335
3336	CALL KILL(R%ac)
3337	r%e_ij=0.0_dp
3338	END subroutine KILL_probe_8
3339
3340	subroutine kill_rf_phasor_8(R)
3341	implicit none
3342	TYPE(rf_phasor_8), INTENT(INOUT) :: R
3343	INTEGER I !,J
3344
3345	DO I=1,2
3346	CALL KILL(R%X(I))
3347	ENDDO
3348	CALL KILL(R%om)
3349	! CALL KILL(R%t)
3350
3351	END subroutine kill_rf_phasor_8
3352
3353
3354
3355
3356	!!!!!!!!!!!!!!!! new stuff
3357	subroutine get_spin_nx_spinor_8(DS,theta0,n0)
3358	implicit none
3359	TYPE(damapspin), INTENT(INout) :: DS
3360	type(real_8), intent(inout) :: theta0
3361	type(spinor_8), intent(inout) :: n0
3362
3363	call get_spin_n0(DS,theta0,n0%x) !get_spin_nx_t
3364
3365	end subroutine get_spin_nx_spinor_8
3366
3367
3368	subroutine get_spin_nx_t(DS,theta0,n0)
3369	implicit none
3370	TYPE(damapspin), INTENT(INout) :: DS
3371	type(damapspin) a1i,ds0
3372	type(real_8), intent(inout) :: theta0,n0(3)
3373	type(real_8) s(3,3),a(3,3),ai(3,3)
3374	type(real_8) a11,a13
3375
3376	call alloc(a1i)
3377	call alloc(ds0)
3378	call alloc_33(s)
3379	call alloc_33(a)
3380	call alloc_33(ai)
3381	call alloc(a11,a13)
3382
3383	call find_n0(ds%s%s,n0)
3384
3385	! call print(theta0,6)
3386
3387
3388	call find_a(n0,a)
3389
3390	call inv_as(a,ai)
3391
3392
3393
3394	call matmulp(ai,ds%s%s,ai) !
3395	call matmulp(ai,a,s) !
3396
3397	a11=(s(1,1))
3398	a13=(s(1,3))
3399
3400	theta0=clockwise*atan2(a13,a11)
3401	if(force_positive.and.theta0<0.0_dp) theta0 = theta0 + twopi !!!! allow negative theta0
3402
3403	! at this stage the spin map is:
3404
3405	! exp( theta0 n0 L ) where n0 is not the real n but it depends on (x,p)
3406
3407
3408	! call print(theta0,6)
3409
3410	! theta0=theta0*nsc
3411
3412	! call print(theta0,6)
3413
3414	! write(6,*) mod((theta0.sub.'0'),twopi)+twopi
3415
3416	call kill(a1i)
3417	call kill(ds0)
3418	call kill_33(s)
3419	call kill_33(a)
3420	call kill_33(ai)
3421
3422	call kill(a11,a13)
3423
3424	end subroutine get_spin_nx_t
3425
3426	subroutine get_spin_nx_r(S,theta0,n0)
3427	implicit none
3428	real(dp), intent(inout) :: theta0,n0(3)
3429	real(dp) a(3,3),ai(3,3)
3430	real(dp), INTENT(INout) :: S(3,3)
3431
3432	call find_n0(s,n0)
3433
3434
3435
3436
3437
3438	call find_a(n0,a)
3439	call inv_as(a,ai)
3440	ai=matmul(ai,s) !
3441	s=matmul(ai,a) !
3442
3443	theta0=clockwise*atan2(s(1,3),s(1,1))
3444	if(force_positive.and.theta0<0.0_dp) theta0 = theta0 + twopi !!!! allow negative theta0
3445
3446
3447	end subroutine get_spin_nx_r
3448
3449	subroutine get_spin_nx_probe(xs0,theta0,n0)
3450	implicit none
3451	type(probe), intent(in) :: xs0
3452	real(dp), intent(inout) :: theta0,n0(3)
3453	integer i,j
3454	real(dp) S(3,3)
3455
3456	do i=1,3
3457	do j=1,3
3458	s(i,j)=xs0%s(j)%x(i)
3459	enddo
3460	enddo
3461
3462	call get_spin_n0(S,theta0,n0)
3463
3464
3465	end subroutine get_spin_nx_probe
3466
3467
3468	subroutine get_spin_nx_rd(dS,theta0,n0)
3469	implicit none
3470	TYPE(damapspin), INTENT(INout) :: DS
3471	real(dp), intent(inout) :: theta0,n0(3)
3472	real(dp) a(3,3),ai(3,3),S(3,3)
3473	integer i,j
3474
3475	do i=1,3
3476	do j=1,3
3477	s(i,j)=ds%s%s(i,j)
3478	enddo
3479	enddo
3480
3481	call find_n0(s,n0)
3482
3483	call find_a(n0,a)
3484	call inv_as(a,ai)
3485	ai=matmul(ai,s) !
3486	s=matmul(ai,a) !
3487
3488	theta0=clockwise*atan2(s(1,3),s(1,1))
3489	if(force_positive.and.theta0<0.0_dp) theta0 = theta0 + twopi !!!! allow negative theta0
3490
3491
3492
3493	end subroutine get_spin_nx_rd
3494
3495
3496
3497	!!!!!!!!!!!!!!
3498	FUNCTION POWMAP( S1, R2 )
3499	implicit none
3500	TYPE (damapspin) POWMAP
3501	TYPE (damapspin), INTENT (IN) :: S1
3502	INTEGER, INTENT (IN) :: R2
3503	TYPE (damapspin) S11
3504	INTEGER I,R22
3505	integer localmaster
3506	IF(.NOT.C_%STABLE_DA) RETURN
3507	localmaster=master
3508
3509
3510	! call checkdamap(s1)
3511
3512	call ass(POWMAP)
3513
3514	call alloc(s11)
3515
3516	s11=1
3517
3518
3519	R22=IABS(R2)
3520	DO I=1,R22
3521	s11=s1*s11
3522	ENDDO
3523
3524	IF(R2.LT.0) THEN
3525	! if(old) then
3526
3527	s11%m=s11%m**(-1)
3528	call trans_mat(s11%s%s,s11%m,s11%s%s)
3529	call inv_as(s11%s%s,s11%s%s)
3530
3531	! CALL etinv(S11%v%i,S11%v%i)
3532	! else
3533	! CALL newetinv(S11%v%j,S11%v%j)
3534	! endif
3535	ENDIF
3536
3537	powmap=s11
3538
3539
3540	! powmap=junk
3541	call kill(s11)
3542
3543	master=localmaster
3544
3545	END FUNCTION POWMAP
3546
3547	FUNCTION dot_spinor( S1, S2 )
3548	implicit none
3549	real(dp) dot_spinor
3550	TYPE (SPINOR), INTENT (IN) :: S1,S2
3551
3552	INTEGER I
3553
3554	dot_spinor=0.0_dp
3555
3556	DO I=1,3
3557	dot_spinor=dot_spinor+s1%x(i)*s2%x(i)
3558	ENDDO
3559
3560
3561
3562	END FUNCTION dot_spinor
3563
3564	FUNCTION dot_spinor_8( S1, S2 )
3565	implicit none
3566	TYPE (real_8) dot_spinor_8
3567	TYPE (SPINOR_8), INTENT (IN) :: S1,S2
3568
3569	INTEGER I
3570	integer localmaster
3571	IF(.NOT.C_%STABLE_DA) RETURN
3572	localmaster=master
3573
3574
3575	! call checkdamap(s1)
3576
3577	call ass(dot_spinor_8)
3578
3579	dot_spinor_8=0.0_dp
3580
3581	DO I=1,3
3582	dot_spinor_8=dot_spinor_8+s1%x(i)*s2%x(i)
3583	ENDDO
3584
3585
3586	master=localmaster
3587
3588	END FUNCTION dot_spinor_8
3589
3590	FUNCTION concat(S2,S1)
3591	implicit none
3592	TYPE (damapspin) concat,t2
3593	TYPE (damapspin), INTENT (IN) :: S1, S2
3594	integer i,j,k
3595	integer localmaster
3596	type(real_8) s(3,3)
3597
3598
3599	IF(.NOT.C_%STABLE_DA) RETURN
3600	localmaster=master
3601
3602	call ass(concat)
3603	call alloc(t2);
3604	call alloc_33(s);
3605
3606	t2=s2
3607	concat=s2
3608
3609
3610	do i=1,3
3611	do j=1,3
3612	! if(s2%s(i,j)%kind==2) then
3613	! if(c_%nd2/=0) then
3614	t2%s%s(i,j)=s2%s%s(i,j)*s1%m
3615	! else
3616	! t2%s(i,j)=s2%s(i,j)
3617	! endif
3618	! else
3619	! t2%s(i,j)=s2%s(i,j)
3620	! endif
3621	enddo
3622	enddo
3623
3624	if(c_%nd2/=0) concat%m=s2%m*s1%m
3625
3626	! do i=1,6
3627	! concat%x(i)=s2%x(i)
3628	! enddo
3629	do i=1,3
3630	do j=1,3
3631	do k=1,3
3632	s(i,j)= t2%s%s(i,k)*s1%s%s(k,j)+ s(i,j)
3633	enddo
3634	enddo
3635	enddo
3636	call smatp(1.0_dp,s,concat%s%s)
3637	! concat%s0=s1%s0
3638
3639	call concat_envelope(S2,S1,concat)
3640
3641
3642	call kill_33(s);
3643	call kill(t2);
3644	master=localmaster
3645
3646	END FUNCTION concat
3647
3648	subroutine concat_envelope(S2,S1,S3)
3649	implicit none
3650	TYPE (damapspin), INTENT (IN) :: S1
3651	TYPE (damapspin), INTENT (IN) :: S2
3652	TYPE (damapspin), INTENT (INout) :: s3
3653	real(dp) s1mi(ndim2,ndim2)
3654	real(dp) e(ndim2,ndim2)
3655
3656	s1mi=0.0_dp
3657	e=0.0_dp
3658
3659	s1mi=(s1%m.sub.1)**(-1)
3660
3661	e(1:c_%nd2,1:c_%nd2)=s2%e_ij
3662	e=matmul(matmul(s1mi,e),transpose(s1mi))
3663
3664	s3%e_ij=e(1:c_%nd2,1:c_%nd2)+s1%e_ij
3665
3666	end subroutine concat_envelope
3667
3668	subroutine extract_envelope_damap(S1,E0_ij,E_ij)
3669	implicit none
3670	TYPE (damapspin), INTENT (IN) :: S1
3671	real(dp), INTENT (in) :: E0_ij(6,6)
3672	real(dp), INTENT (out) :: E_ij(6,6)
3673	real(dp) s1m(ndim2,ndim2)
3674	real(dp) e(ndim2,ndim2)
3675
3676	s1m=0.0_dp
3677	e=0.0_dp
3678
3679
3680	s1m=s1%m.sub.1
3681
3682
3683	e(1:c_%nd2,1:c_%nd2)=s1%e_ij+E0_ij
3684	e=matmul(matmul(s1m,e),transpose(s1m))
3685
3686	E_ij=e(1:c_%nd2,1:c_%nd2)
3687
3688	end subroutine extract_envelope_damap
3689
3690	subroutine extract_envelope_probe8(p,E0_ij,E_ij)
3691	implicit none
3692	TYPE (probe_8), INTENT (IN) :: p
3693	TYPE (damapspin) S1
3694	real(dp), INTENT (in) :: E0_ij(6,6)
3695	real(dp), INTENT (out) :: E_ij(6,6)
3696	real(dp) s1m(ndim2,ndim2)
3697	real(dp) e(ndim2,ndim2)
3698
3699	call alloc(s1)
3700	s1=p
3701
3702	call extract_envelope_damap(s1,E0_ij,E_ij)
3703
3704	call kill(s1)
3705	end subroutine extract_envelope_probe8
3706
3707
3708	FUNCTION cmul(S2,S1) ! multiply spin part with real(dp) s1
3709	implicit none
3710	TYPE (damapspin) cmul
3711	TYPE (damapspin), INTENT (IN) :: S1
3712	real(dp), INTENT (IN) :: S2
3713	integer i,j
3714	integer localmaster
3715
3716
3717	IF(.NOT.C_%STABLE_DA) RETURN
3718	localmaster=master
3719
3720	call ass(cmul)
3721
3722	if(c_%nd2/=0) cmul%m=s1%m
3723
3724
3725
3726	do i=1,3
3727	do j=1,3
3728	cmul%s%s(i,j)=s2*s1%s%s(i,j)
3729	enddo
3730	enddo
3731	master=localmaster
3732
3733	END FUNCTION cmul
3734
3735	FUNCTION mmul(S1,S2) ! transform spin part with damap s2
3736	implicit none
3737	TYPE (damapspin) mmul
3738	TYPE (damapspin), INTENT (IN) :: S1
3739	type(damap), INTENT (IN) :: s2
3740	! integer i,j,k
3741	integer localmaster
3742
3743
3744	IF(.NOT.C_%STABLE_DA) RETURN
3745	localmaster=master
3746
3747	call ass(mmul)
3748
3749	mmul=s1
3750	call trans_mat( mmul%s%s,s2, mmul%s%s)
3751
3752	! do i=1,3
3753	! do j=1,3
3754	! if(s2%s(i,j)%kind==2) then
3755	! mmul%s(i,j)=s1%s(i,j)*s2
3756	! enddo
3757	! enddo
3758
3759	master=localmaster
3760
3761	END FUNCTION mmul
3762
3763	FUNCTION damapspin_spinor8_mul(S1,S2) ! transform spin part with spinor_8 s2 and returns spinor_8
3764	implicit none
3765	type(spinor_8) damapspin_spinor8_mul
3766	TYPE (damapspin), INTENT (IN) :: S1
3767	type(spinor_8), INTENT (IN) :: S2
3768	integer i,j
3769	integer localmaster
3770
3771
3772	IF(.NOT.C_%STABLE_DA) RETURN
3773	localmaster=master
3774
3775	call assp_master(damapspin_spinor8_mul%x(1))
3776	call assp_no_master(damapspin_spinor8_mul%x(2))
3777	call assp_no_master(damapspin_spinor8_mul%x(3))
3778
3779	damapspin_spinor8_mul%x(1)=0.0_dp
3780	damapspin_spinor8_mul%x(2)=0.0_dp
3781	damapspin_spinor8_mul%x(3)=0.0_dp
3782
3783	do i=1,3
3784	do j=1,3
3785	damapspin_spinor8_mul%x(i)=(s1%s%s(i,j))*S2%x(j)+damapspin_spinor8_mul%x(i)
3786	enddo
3787	enddo
3788
3789	master=localmaster
3790
3791	END FUNCTION damapspin_spinor8_mul
3792
3793
3794	FUNCTION damapspin_spinor_mul(S1,S2) ! transform spin part with spinor_8 s2 and returns spinor_8
3795	implicit none
3796	type(spinor_8) damapspin_spinor_mul
3797	TYPE (damapspin), INTENT (IN) :: S1
3798	type(spinor), INTENT (IN) :: S2
3799	integer i,j
3800	integer localmaster
3801
3802
3803
3804	IF(.NOT.C_%STABLE_DA) RETURN
3805	localmaster=master
3806
3807	call assp_master(damapspin_spinor_mul%x(1))
3808	call assp_no_master(damapspin_spinor_mul%x(2))
3809	call assp_no_master(damapspin_spinor_mul%x(3))
3810
3811	damapspin_spinor_mul%x(1)=0.0_dp
3812	damapspin_spinor_mul%x(2)=0.0_dp
3813	damapspin_spinor_mul%x(3)=0.0_dp
3814
3815	do i=1,3
3816	do j=1,3
3817	damapspin_spinor_mul%x(i)=(s1%s%s(i,j))*S2%x(j)+damapspin_spinor_mul%x(i)
3818	enddo
3819	enddo
3820
3821	master=localmaster
3822
3823
3824	END FUNCTION damapspin_spinor_mul
3825
3826	FUNCTION exp_spinor_8(S1) ! transform spin part with damap s2
3827	implicit none
3828	TYPE (damapspin) exp_spinor_8
3829	TYPE (spinor_8), INTENT (IN) :: S1
3830	! integer i,j,k
3831	integer localmaster
3832
3833
3834	IF(.NOT.C_%STABLE_DA) RETURN
3835	localmaster=master
3836
3837	call ass(exp_spinor_8)
3838
3839	call daexplogp(S1,exp_spinor_8)
3840
3841	master=localmaster
3842
3843	END FUNCTION exp_spinor_8
3844
3845
3846	FUNCTION spin8_mul_map(S2,S1) ! transforms spinor_8 s2 with damap s1
3847	implicit none
3848	type(spinor_8) spin8_mul_map
3849	TYPE (damap), INTENT (IN) :: S1
3850	type(spinor_8), INTENT (IN) :: S2
3851	integer i
3852	integer localmaster
3853
3854
3855	IF(.NOT.C_%STABLE_DA) RETURN
3856	localmaster=master
3857
3858	call assp_master(spin8_mul_map%x(1))
3859	call assp_no_master(spin8_mul_map%x(2))
3860	call assp_no_master(spin8_mul_map%x(3))
3861
3862	spin8_mul_map%x(1)=0.0_dp
3863	spin8_mul_map%x(2)=0.0_dp
3864	spin8_mul_map%x(3)=0.0_dp
3865
3866	do i=1,3
3867	if(S2%x(i)%kind==2) then
3868	spin8_mul_map%x(i)=S2%x(i)%t*s1
3869	else
3870	spin8_mul_map%x(i)=S2%x(i)
3871	endif
3872	enddo
3873
3874	master=localmaster
3875
3876	END FUNCTION spin8_mul_map
3877
3878	FUNCTION spin8_scal8_map(S1,S2) ! transforms spinor_8 s2 with damap s1
3879	implicit none
3880	type(spinor_8) spin8_scal8_map
3881	TYPE (real_8), INTENT (IN) :: S1
3882	type(spinor_8), INTENT (IN) :: S2
3883	integer i
3884	integer localmaster
3885
3886
3887	IF(.NOT.C_%STABLE_DA) RETURN
3888	localmaster=master
3889
3890	call assp_master(spin8_scal8_map%x(1))
3891	call assp_no_master(spin8_scal8_map%x(2))
3892	call assp_no_master(spin8_scal8_map%x(3))
3893
3894	spin8_scal8_map%x(1)=0.0_dp
3895	spin8_scal8_map%x(2)=0.0_dp
3896	spin8_scal8_map%x(3)=0.0_dp
3897
3898	do i=1,3
3899	spin8_scal8_map%x(i)=s1*S2%x(i)
3900	enddo
3901
3902	master=localmaster
3903
3904	END FUNCTION spin8_scal8_map
3905
3906
3907	FUNCTION add_spin8_spin8(S1,S2) ! transforms spinor_8 s2 with damap s1
3908	implicit none
3909	type(spinor_8) add_spin8_spin8
3910	type(spinor_8), INTENT (IN) :: S1,S2
3911	integer i
3912	integer localmaster
3913
3914
3915	IF(.NOT.C_%STABLE_DA) RETURN
3916	localmaster=master
3917
3918	call assp_master(add_spin8_spin8%x(1))
3919	call assp_no_master(add_spin8_spin8%x(2))
3920	call assp_no_master(add_spin8_spin8%x(3))
3921
3922	add_spin8_spin8%x(1)=0.0_dp
3923	add_spin8_spin8%x(2)=0.0_dp
3924	add_spin8_spin8%x(3)=0.0_dp
3925
3926	do i=1,3
3927	add_spin8_spin8%x(i)=s1%x(i)+S2%x(i)
3928	enddo
3929
3930	master=localmaster
3931
3932	END FUNCTION add_spin8_spin8
3933
3934	FUNCTION mul_spin8_spin8(S1,S2) ! transforms spinor_8 s2 with damap s1
3935	implicit none
3936	type(spinor_8) mul_spin8_spin8
3937	type(spinor_8), INTENT (IN) :: S1,S2
3938	integer i
3939	integer localmaster
3940
3941	IF(.NOT.C_%STABLE_DA) RETURN
3942	localmaster=master
3943
3944	call assp_master(mul_spin8_spin8%x(1))
3945	call assp_no_master(mul_spin8_spin8%x(2))
3946	call assp_no_master(mul_spin8_spin8%x(3))
3947
3948	mul_spin8_spin8%x(1)=0.0_dp
3949	mul_spin8_spin8%x(2)=0.0_dp
3950	mul_spin8_spin8%x(3)=0.0_dp
3951
3952	mul_spin8_spin8%x(1)=s1%x(2)S2%x(3)-s1%x(3)S2%x(2)
3953	mul_spin8_spin8%x(2)=s1%x(3)S2%x(1)-s1%x(1)S2%x(3)
3954	mul_spin8_spin8%x(3)=s1%x(1)S2%x(2)-s1%x(2)S2%x(1)
3955
3956	master=localmaster
3957
3958	END FUNCTION mul_spin8_spin8
3959
3960	FUNCTION sub_spin8_spin8(S1,S2) ! transforms spinor_8 s2 with damap s1
3961	implicit none
3962	type(spinor_8) sub_spin8_spin8
3963	type(spinor_8), INTENT (IN) :: S1,S2
3964	integer i
3965	integer localmaster
3966
3967
3968	IF(.NOT.C_%STABLE_DA) RETURN
3969	localmaster=master
3970
3971	call assp_master(sub_spin8_spin8%x(1))
3972	call assp_no_master(sub_spin8_spin8%x(2))
3973	call assp_no_master(sub_spin8_spin8%x(3))
3974
3975	sub_spin8_spin8%x(1)=0.0_dp
3976	sub_spin8_spin8%x(2)=0.0_dp
3977	sub_spin8_spin8%x(3)=0.0_dp
3978
3979	do i=1,3
3980	sub_spin8_spin8%x(i)=s1%x(i)-S2%x(i)
3981	enddo
3982
3983	master=localmaster
3984
3985	END FUNCTION sub_spin8_spin8
3986
3987	function eval_spinor_8(s,x)
3988	implicit none
3989	TYPE(spinor) eval_spinor_8
3990	TYPE(spinor_8), INTENT(IN) :: s
3991	real(dp), INTENT(IN) :: x(lnv)
3992	integer i
3993
3994
3995	do i=1,3
3996	if(s%x(i)%kind==2) then
3997	eval_spinor_8%x(i)=s%x(i)%t*x
3998	endif
3999	enddo
4000
4001	end function eval_spinor_8
4002
4003	FUNCTION concatxp(S2,x)
4004	implicit none
4005	TYPE (damapspin) concatxp
4006	TYPE (damapspin), INTENT (IN) :: S2
4007	real(dp), INTENT (IN) :: x(lnv)
4008	integer localmaster
4009
4010
4011	IF(.NOT.C_%STABLE_DA) RETURN
4012	localmaster=master
4013
4014	call ass(concatxp)
4015
4016	call eval_spin_matrix(S2,x,concatxp)
4017
4018	master=localmaster
4019
4020	END FUNCTION concatxp
4021
4022
4023	FUNCTION addm(S2,S1) ! adds spin part of s1 and s2
4024	implicit none
4025	TYPE (damapspin) addm
4026	TYPE (damapspin), INTENT (IN) :: S1, S2
4027	integer i,j
4028	integer localmaster
4029
4030
4031	IF(.NOT.C_%STABLE_DA) RETURN
4032	localmaster=master
4033
4034	call ass(addm)
4035
4036	addm=s2
4037
4038	do i=1,3
4039	do j=1,3
4040	addm%s%s(i,j)=s2%s%s(i,j)+s1%s%s(i,j)
4041	enddo
4042	enddo
4043
4044
4045	master=localmaster
4046
4047	END FUNCTION addm
4048
4049
4050	real(dp) function purge_transverse(j)
4051	implicit none
4052	integer i
4053	! INTEGER J(NTT)
4054	integer,dimension(:)::j
4055	if(.not.c_%stable_da) return
4056
4057	purge_transverse=1.0_dp
4058
4059
4060	do i=1,c_%nd2
4061
4062	if(i/=c_%ndpt) then
4063	if(j(i)/=0) then
4064	purge_transverse=0.0_dp
4065	exit
4066	endif
4067	endif
4068	enddo
4069
4070	return
4071	end function purge_transverse
4072
4073	subroutine clean_orbital_33(s,sf)
4074	implicit none
4075	TYPE (real_8), intent(inout) :: s(3,3),sf(3,3)
4076	integer i,j
4077	type(taylor) t
4078	logical doflip
4079	call alloc(t)
4080
4081	if(perform_flip.and.new_ndpt.and.c_%ndpt/=0) then
4082	perform_flip=.false.
4083	do i=1,3
4084	do j=1,3
4085	if(s(i,j)%kind==2) call fliptaylor(s(i,j)%t,s(i,j)%t,1)
4086	enddo
4087	enddo
4088	doflip=.true.
4089	else
4090	doflip=.false.
4091	endif
4092
4093	do i=1,3
4094	do j=1,3
4095	if(s(i,j)%kind==2) then
4096	call cfu(s(i,j)%t,purge_transverse,t)
4097	sf(i,j)=t
4098	else
4099	sf(i,j)=s(i,j)
4100	endif
4101	enddo
4102	enddo
4103	call kill(t)
4104
4105	if(doflip) then
4106	do i=1,3
4107	do j=1,3
4108	if(s(i,j)%kind==2) call fliptaylor(s(i,j)%t,s(i,j)%t,-1)
4109	if(s(i,j)%kind==2.and.s(i,j)%t%i/=sf(i,j)%t%i) then
4110	call fliptaylor(sf(i,j)%t,sf(i,j)%t,-1)
4111	endif
4112	enddo
4113	enddo
4114	perform_flip=.true.
4115	endif
4116
4117
4118	end subroutine clean_orbital_33
4119
4120	!!!!!!!!!! Normal form !!!!!!!!!!!!!!!
4121
4122	subroutine alloc_normal_spin(D)
4123	implicit none
4124	TYPE(normal_spin), INTENT(INOUT) :: D
4125
4126	CALL alloc(D%n0,3)
4127	CALL alloc(D%theta0)
4128	CALL alloc(D%n)
4129	CALL alloc(D%a_t)
4130	CALL alloc(D%as)
4131	CALL alloc(D%a1)
4132	CALL alloc(D%ar)
4133	D%NRES=0
4134	D%M=0
4135	D%Ms=0
4136	D%s_ij0=0.0_dp
4137	D%s_ijr=0.0_dp
4138	D%emittance=0.0_dp
4139	D%tune=0.0_dp
4140	D%damping=0.0_dp
4141	D%AUTO=my_true
4142	D%STOCHASTIC=my_false
4143	D%STOCH=0.0_dp
4144	D%STOCH_inv=0.0_dp
4145	D%nu=0.0_dp
4146
4147	END subroutine alloc_normal_spin
4148
4149	subroutine KILL_normal_spin(D)
4150	implicit none
4151	TYPE(normal_spin), INTENT(INOUT) :: D
4152
4153	CALL KILL(D%n0,3)
4154	CALL KILL(D%theta0)
4155	CALL KILL(D%n)
4156	CALL KILL(D%a_t)
4157	CALL KILL(D%as)
4158	CALL KILL(D%a1)
4159	CALL KILL(D%ar)
4160
4161	END subroutine KILL_normal_spin
4162
4163	!!!!!!!!!!!!!!!! new stuff
4164	subroutine Go_to_closed(ns,DS,spin_in)
4165	implicit none
4166	TYPE(normal_spin), INTENT(INOUT) :: ns
4167	TYPE(damapspin), INTENT(INout) :: DS
4168	logical(lp), INTENT(IN) :: spin_in
4169	logical(lp) rad_in
4170	type(damapspin) a1i,ds0
4171	type(damapspin)s ,a ,ai
4172	type(taylor) nn
4173	type(real_8) a11,a13
4174	integer i,j,jj(lnv)
4175	real(dp) ss(3,3)
4176	!!!!!!!!!! wrapping for radiation !!!!!!!!!!!!!!!!!!!
4177	type(radtaylor) ys(ndim2)
4178	type(beamenvelope) env
4179
4180
4181	call alloc(a1i)
4182	call alloc(ds0)
4183	call alloc(s)
4184	call alloc(a)
4185	call alloc(ai)
4186	call alloc(nn)
4187	call alloc(a11,a13)
4188
4189	ns%as=1 ! damapspin to fix point
4190
4191	ns%n=ds%m ! ORBITAL NORMALIZATION
4192
4193	ns%as%m=ns%n%a1 ! fix point map a1 of orbital normal form
4194
4195	a1i=ns%as**(-1) ! fix point map a1 of orbital normal form
4196
4197
4198	DS0=a1idsns%as
4199	! DS0=a1idsns%as ! at this stage ds0 is the map around the fixed point
4200	! but spin matrix unchanged except for its dependence on transverse
4201
4202	!!!!!!!!!! wrapping for radiation !!!!!!!!!!!!!!!!!!!
4203	!!!!!!!!!! later to be change changed !!!!!!!!!!!!!!!
4204	call check_rad(ds%e_ij,rad_in)
4205	if(rad_in) then
4206	if(c_%no==1) then
4207	call normalize_envelope(ns,DS)
4208	Write(6,*) "New envelope calculation attempted with NO=1 "
4209	else
4210	call alloc(ys,6)
4211	call alloc(env)
4212	jj=0
4213	do i=1,6
4214	ys(i)%v=ds%m%v(i)
4215	do j=1,6
4216	ys(i)%e(j)=ds%e_ij(i,j)
4217	enddo
4218	enddo
4219	env%stochastic=ns%stochastic
4220	env%auto=ns%auto
4221	env=ys
4222
4223	ns%s_ij0 = env%s_ij0
4224	ns%emittance = env%emittance
4225	ns%KICK = env%KICK
4226	do i=1,6
4227	do j=1,6
4228	jj(j)=1
4229	ns%STOCH(i,j) = env%STOCH%v(i).sub.jj
4230	jj(j)=0
4231	enddo
4232	enddo
4233
4234	call kill(env)
4235	call kill(ys,6)
4236	endif
4237	endif
4238	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4239	s=1
4240	a=1
4241	ai=1
4242
4243	call clean_orbital_33(ds0%s%s,s%s%s)
4244	! at this stage s is the spin map without dependence on orbital
4245
4246	if(spin_in) then !!!! spin in
4247	call find_n0(s%s%s,ns%n0)
4248
4249	call find_a(ns%n0,a%s%s)
4250
4251	ai=a**(-1)
4252
4253	s=(ais)a
4254
4255	!write(6,*) " diagonal "
4256	! do i=1,3
4257	! do j=1,3
4258	! ss(i,j)=s%s(i,j)
4259	! enddo
4260	! write(6,*) ss(i,1:3)
4261	! enddo
4262
4263	!pause 777
4264
4265	a11=s%s%s(1,1)
4266	a13=s%s%s(1,3)
4267
4268	ns%theta0=clockwise*atan2(a13,a11)
4269	ns%nu=ns%theta0/twopi
4270	if(force_positive.and.ns%theta0<0.0_dp) ns%theta0 = ns%theta0 + twopi !!!! allow negative theta0
4271
4272	ns%as=ns%as*a
4273
4274	ns%ar=1
4275	ns%ar%m=ns%as%m*(-1)ns%n%a_t
4276	else
4277	ns%ar=1
4278	ns%ar%m=ns%as%m*(-1)ns%n%a_t
4279
4280	ns%theta0=0.0_dp ! arbitrary junk
4281	ns%n0(1)=0.0_dp ! arbitrary junk
4282	ns%n0(2)=1.0_dp ! arbitrary junk
4283	ns%n0(3)=0.0_dp ! arbitrary junk
4284	endif
4285
4286	call kill(a1i)
4287	call kill(ds0)
4288	call kill(s)
4289	call kill(a)
4290	call kill(ai)
4291
4292	call kill(a11,a13)
4293	call kill(nn)
4294
4295	end subroutine Go_to_closed
4296
4297	subroutine normalize_envelope(norm_spin,m_spin)
4298	implicit none
4299	type(normal_spin) norm_spin
4300	type(damapspin) m_spin
4301	type(normalform) norm
4302	type(damap) m
4303	integer i,j,i1,i2
4304	real(dp) a(6,6),ai(6,6),ait(6,6),mat(6,6), sigma_inf(6,6),at(6,6),br(6,6)
4305	complex(dp) c(6,6),ci(6,6),cit(6,6), b(6,6),ct(6,6)
4306	complex(dp) coef,ba(6,6),b_phasor(6,6)
4307	complex(dp) R(6,6),r_phasor(6,6), sigma_inf_phasor(6,6)
4308	real(dp) xj(6,6),mj(6,6),xn,jb(6,6),bs(6,6)
4309
4310	b=m_spin%e_ij
4311	bs=m_spin%e_ij
4312
4313	c=0.0_dp
4314	ci=0.0_dp
4315	do i=1,3
4316	do j=1,3
4317	xj(2i,2i-1)=-1.0_dp
4318	xj(2i-1,2i)=1.0_dp
4319	c(2i-1,2i-1)=0.5_dp
4320	c(2i-1,2i)=0.5_dp
4321	c(2i,2i-1)=0.5_dp/i_
4322	c(2i,2i)=-0.5_dp/i_
4323	ci(2i-1,2i-1)=1.0_dp
4324	ci(2i-1,2i)=i_
4325	ci(2i,2i-1)=1.0_dp
4326	ci(2i,2i)=-i_
4327	enddo
4328	enddo
4329
4330	mat=m_spin%m
4331	a=norm_spin%n%a_t
4332	ai=norm_spin%n%a_t**(-1)
4333	ait=transpose(ai)
4334	at=transpose(a)
4335	cit=transpose(ci)
4336	ct=transpose(c)
4337
4338	R=matmul(matmul(ai,b),ait)
4339
4340	ba=matmul(matmul(ai,b),ait)
4341
4342	b_phasor=matmul(matmul(ci,ba),cit)
4343
4344	r=matmul(matmul(ai,mat),a)
4345	r_phasor=matmul(matmul(ci,r),c)
4346
4347	do i=1,6
4348	do j=1,6
4349	sigma_inf_phasor(i,j)= r_phasor(i,i)r_phasor(j,j)/(1.0_dp-r_phasor(i,i)r_phasor(j,j))*b_phasor(i,j)
4350	enddo
4351	enddo
4352	do i=1,3
4353	norm_spin%emittance(i)=sigma_inf_phasor(2i-1,2i)/2.0_dp
4354	enddo
4355
4356	sigma_inf=matmul(matmul(c,sigma_inf_phasor),ct)
4357	sigma_inf=matmul(matmul(a,sigma_inf),at)
4358
4359	norm_spin%s_ij0=sigma_inf
4360	norm_spin%s_ijr=sigma_inf_phasor
4361
4362	norm_spin%tune=norm_spin%n%tune(1:3)
4363	norm_spin%damping=norm_spin%n%damping(1:3)
4364
4365
4366
4367	if(norm_spin%STOCHASTIC) then
4368	call diagonalise_envelope_a(bs,br,a,ai,norm_spin%kick) !diagonalise_envelope_a(b,br,a,ai,kick)
4369	norm_spin%STOCH=a
4370	norm_spin%STOCH_inv=ai
4371	endif
4372
4373	end subroutine normalize_envelope
4374
4375	subroutine fetch_s0(DS,s)
4376	implicit none
4377	TYPE(damapspin), INTENT(INout) :: DS
4378	TYPE(damapspin), INTENT(INout) :: s
4379	TYPE(damapspin) s0
4380	integer i,j
4381
4382	call alloc(s0)
4383
4384	s0=1
4385
4386	do i=1,3
4387	do j=1,3
4388	s0%s%s(i,j)=ds%s%s(i,j).sub.'0'
4389	enddo
4390	enddo
4391	s=s0
4392	call kill(s0)
4393
4394	end subroutine fetch_s0
4395
4396
4397	subroutine dalog_spinor_8(DS,n)
4398	implicit none
4399	TYPE(damapspin), INTENT(INout) :: DS
4400	TYPE(spinor_8), INTENT(INout) :: n
4401	TYPE(taylor) om(3)
4402	integer i
4403	call alloc(om)
4404
4405	call dalog(DS,om)
4406
4407	do i=1,3
4408	n%x(i)=morph(om(i))
4409	enddo
4410
4411	call kill(om)
4412	end subroutine dalog_spinor_8
4413
4414	subroutine factor_parameter_dependent_s0(DS,s0,NS,N_AXIS,DIR)
4415	implicit none
4416	TYPE(damapspin), INTENT(INout) :: DS,s0,NS
4417	TYPE(damapspin)A_F,A_S
4418	TYPE(spinor_8), INTENT(INout) :: N_AXIS
4419	INTEGER DIR
4420
4421	CALL ALLOC(A_F,A_S)
4422
4423	A_f=DS
4424	a_f%m=1
4425	CALL clean_orbital_33(A_f%S%s,A_f%S%s)
4426	A_S=DS
4427	a_s%m=1
4428	IF(DIR==1) THEN
4429	a_S=a_f*(-1)a_s ! =s0*ns
4430	ELSE
4431	a_S=a_sa_f(-1) ! =nss0 of course s0 is different
4432	ENDIF
4433	CALL dalog_spinor_8(a_S,N_AXIS) ! S=exp(N_AXIS)
4434	S0=A_F
4435	NS=A_S
4436
4437
4438	CALL KILL(A_F,A_S)
4439
4440	END subroutine factor_parameter_dependent_s0
4441
4442
4443	subroutine factor_s0(DS,s0,NS,N_AXIS,DIR)
4444	implicit none
4445	TYPE(damapspin), INTENT(INout) :: DS,s0,NS
4446	TYPE(damapspin)A_F,A_S
4447	TYPE(spinor_8), INTENT(INout) :: N_AXIS
4448	INTEGER DIR
4449
4450	CALL ALLOC(A_F,A_S)
4451
4452	A_f=DS
4453	a_f%m=1
4454	CALL fetch_s0(A_f,A_f)
4455	A_S=DS
4456	a_s%m=1
4457	IF(DIR==1) THEN
4458	a_S=a_fa_s*(-1) ! angle has no constant part (DA) BUT PARAMETERS
4459	ELSE
4460	a_S=a_f*(-1)a_s ! angle has no constant part (DA)
4461	ENDIF
4462	CALL dalog_spinor_8(a_S,N_AXIS) ! S=exp(N_AXIS)
4463	S0=A_F
4464	NS=A_S
4465
4466
4467	CALL KILL(A_F,A_S)
4468
4469	END subroutine factor_s0
4470
4471
4472	subroutine dalog(DS,om)
4473	implicit none
4474	TYPE(damapspin), INTENT(INout) :: DS
4475	TYPE(taylor), INTENT(INout) :: om(3)
4476	TYPE(damapspin) h
4477	integer i
4478	TYPE(damapspin) dh,dhn
4479	real(dp) c,cl
4480	call alloc(h)
4481	call alloc(dh)
4482	call alloc(dhn)
4483	! this only works with a da-map
4484
4485	dh=ds
4486	dh%m=1
4487	h=0
4488	dhn=1
4489	dh%s%s(1,1)=dh%s%s(1,1)-1.0_dp
4490	dh%s%s(2,2)=dh%s%s(2,2)-1.0_dp
4491	dh%s%s(3,3)=dh%s%s(3,3)-1.0_dp
4492
4493	c=1.0_dp
4494	do i=1,c_%no
4495	dhn=dhn*dh
4496	cl=c/i
4497	h=h+cl*dhn
4498	c=-c
4499	enddo
4500
4501	om(1)=h%s%s(3,2)
4502	om(2)=h%s%s(1,3)
4503	om(3)=h%s%s(2,1)
4504
4505	call kill(h)
4506	call kill(dh)
4507	call kill(dhn)
4508
4509	end subroutine dalog
4510
4511	subroutine daexplog(om,ds)
4512	implicit none
4513	TYPE(damapspin), INTENT(INout) :: DS
4514	TYPE(taylor), INTENT(INout) :: om(3)
4515	integer i
4516	TYPE(damapspin) dh,dhn
4517	real(dp) c
4518	call alloc(dh)
4519	call alloc(dhn)
4520	! this only works with a da-map
4521	ds=1
4522	dh%m=1
4523	dh%s%s(2,1)=om(3)
4524	dh%s%s(1,3)=om(2)
4525	dh%s%s(3,2)=om(1)
4526	dh%s%s(1,2)=-om(3)
4527	dh%s%s(3,1)=-om(2)
4528	dh%s%s(2,3)=-om(1)
4529	dhn=1
4530	c=1.0_dp
4531	do i=1,c_%no
4532	dhn=dhn*dh
4533	c=c/i
4534	ds=ds+c*dhn
4535	enddo
4536
4537
4538	call kill(dh)
4539	call kill(dhn)
4540
4541	end subroutine daexplog
4542
4543	subroutine get_kernel(ns,om,oma)
4544	implicit none
4545	TYPE(normal_spin), INTENT(INout) :: ns
4546	TYPE(taylor), INTENT(INout) :: om(3),oma(3)
4547	type(taylorresonance) t
4548	TYPE(complextaylor) omr(3),omc(3)
4549	integer N,i,nd,j
4550	integer, allocatable :: jc(:)
4551	real(dp) value,ang,tune(4)
4552	complex(dp) denom
4553	logical doit,doflip
4554
4555	if(perform_flip.and.new_ndpt.and.c_%ndpt/=0) then
4556	perform_flip=.false.
4557	call fliptaylor(om(1),om(1),1)
4558	call fliptaylor(om(2),om(2),1)
4559	call fliptaylor(om(3),om(3),1)
4560	call flip_real_array(ns%n%tune,ns%n%tune,1)
4561	if(use_ptc_ac_position) then
4562	call flip_resonance(ns%m,ns%m,1)
4563	endif
4564	doflip=.true.
4565	else
4566	doflip=.false.
4567	endif
4568
4569
4570	nd=c_%nd2/2
4571	if(c_%ndpt/=0) nd=nd-1
4572
4573	call alloc(omc,3)
4574	call alloc(omr,3)
4575	call alloc(t)
4576	! write(6,*) c_%nd2,c_%ndpt,nd
4577	! write(6,*) ns%n%tune
4578
4579
4580	omc(1)= om(1)-i_*om(3) ! coeff of lambdda
4581	omc(3)= om(1)+i_om(3) ! coeff of lambdda
4582	omc(2)= om(2)
4583
4584	allocate(jc(c_%nv))
4585
4586
4587	! om(2) is always real!!!
4588	t=omc(2)%r
4589	! cos part of omc(2) !!!
4590
4591	call taylor_cycle(t%cos,size=n)
4592
4593	do i=1,n
4594	call taylor_cycle(t%cos,ii=i,value=value,j=jc)
4595
4596	call test_jc(ns,jc,nd,doit)
4597
4598	if(doit) then
4599	ang=0.0_dp
4600	do j=1,nd
4601	ang=(jc(j2-1)-jc(j2))twopins%n%tune(j)+ang
4602	enddo
4603	denom=value/(exp(-i_*ang)-1.0_dp)
4604	omr(2)=omr(2)+ (denom.mono.jc)
4605	endif
4606	enddo
4607
4608	! sin part of om(2)
4609	call taylor_cycle(t%sin,size=N)
4610
4611	do i=1,n
4612	call taylor_cycle(t%sin,ii=i,value=value,j=jc)
4613
4614	call test_jc(ns,jc,nd,doit)
4615
4616	if(doit) then
4617	ang=0.0_dp
4618	do j=1,nd
4619	ang=(jc(j2-1)-jc(j2))twopins%n%tune(j)+ang
4620	enddo
4621	denom=value/(exp(-i_*ang)-1.0_dp)
4622	omr(2)=omr(2)+ i_*(denom.mono.jc)
4623	endif
4624	enddo
4625	doit=.true.
4626	! real part of om(1)
4627	t=omc(1)%r
4628	call taylor_cycle(t%cos,N)
4629
4630	do i=1,n
4631	call taylor_cycle(t%cos,ii=i,value=value,j=jc)
4632
4633	call test_jc_spin(ns,jc,clockwise,nd,doit)
4634
4635	if(doit) then
4636	ang=0.0_dp
4637	do j=1,nd
4638	ang=(jc(j2-1)-jc(j2))twopins%n%tune(j)+ang
4639	enddo
4640	denom=value/(exp(-i_ang)-exp(clockwisei_*ns%theta0))
4641	omr(1)=omr(1)+ (denom.mono.jc)
4642	endif
4643	enddo
4644
4645	call taylor_cycle(t%SIN,N)
4646
4647	do i=1,n
4648	call taylor_cycle(t%SIN,ii=i,value=value,j=jc)
4649
4650	call test_jc_spin(ns,jc,clockwise,nd,doit)
4651
4652	if(doit) then
4653	ang=0.0_dp
4654	do j=1,nd
4655	ang=(jc(j2-1)-jc(j2))twopins%n%tune(j)+ang
4656	enddo
4657	denom=value/(exp(-i_ang)-exp(clockwisei_*ns%theta0))
4658	omr(1)=omr(1)+ I_*(denom.mono.jc)
4659	endif
4660	enddo
4661
4662	! imaginary part of om(1)
4663	t=omc(1)%i
4664	call taylor_cycle(t%cos,N)
4665
4666	do i=1,n
4667	call taylor_cycle(t%cos,ii=i,value=value,j=jc)
4668
4669	call test_jc_spin(ns,jc,clockwise,nd,doit)
4670
4671	if(doit) then
4672	ang=0.0_dp
4673	do j=1,nd
4674	ang=(jc(j2-1)-jc(j2))twopins%n%tune(j)+ang
4675	enddo
4676	denom=value/(exp(-i_ang)-exp(clockwisei_*ns%theta0))
4677	omr(1)=omr(1)+ i_*(denom.mono.jc)
4678	endif
4679	enddo
4680
4681	call taylor_cycle(t%SIN,N)
4682
4683	do i=1,n
4684	call taylor_cycle(t%SIN,ii=i,value=value,j=jc)
4685
4686	call test_jc_spin(ns,jc,clockwise,nd,doit)
4687
4688	if(doit) then
4689	ang=0.0_dp
4690	do j=1,nd
4691	ang=(jc(j2-1)-jc(j2))twopins%n%tune(j)+ang
4692	enddo
4693	denom=value/(exp(-i_ang)-exp(clockwisei_*ns%theta0))
4694	omr(1)=omr(1)-(denom.mono.jc)
4695	endif
4696	enddo
4697
4698	! real part of om(3)
4699	t=omc(3)%r
4700	call taylor_cycle(t%cos,N)
4701
4702	do i=1,n
4703	call taylor_cycle(t%cos,ii=i,value=value,j=jc)
4704
4705	call test_jc_spin(ns,jc,-clockwise,nd,doit)
4706
4707	if(doit) then
4708	ang=0.0_dp
4709	do j=1,nd
4710	ang=(jc(j2-1)-jc(j2))twopins%n%tune(j)+ang
4711	enddo
4712	denom=value/(exp(-i_ang)-exp(-clockwisei_*ns%theta0))
4713	omr(3)=omr(3)+ (denom.mono.jc)
4714	endif
4715	enddo
4716
4717	call taylor_cycle(t%SIN,N)
4718
4719	do i=1,n
4720	call taylor_cycle(t%SIN,ii=i,value=value,j=jc)
4721
4722	call test_jc_spin(ns,jc,-clockwise,nd,doit)
4723
4724	if(doit) then
4725	ang=0.0_dp
4726	do j=1,nd
4727	ang=(jc(j2-1)-jc(j2))twopins%n%tune(j)+ang
4728	enddo
4729	denom=value/(exp(-i_ang)-exp(-clockwisei_*ns%theta0))
4730	omr(3)=omr(3)+ I_*(denom.mono.jc)
4731	endif
4732	enddo
4733
4734	! imaginary part of om(3)
4735	t=omc(3)%i
4736	call taylor_cycle(t%cos,N)
4737
4738	do i=1,n
4739	call taylor_cycle(t%cos,ii=i,value=value,j=jc)
4740
4741	call test_jc_spin(ns,jc,-clockwise,nd,doit)
4742
4743	if(doit) then
4744	ang=0.0_dp
4745	do j=1,nd
4746	ang=(jc(j2-1)-jc(j2))twopins%n%tune(j)+ang
4747	enddo
4748	denom=value/(exp(-i_ang)-exp(-clockwisei_*ns%theta0))
4749	omr(3)=omr(3)+ i_*(denom.mono.jc)
4750	endif
4751	enddo
4752
4753	call taylor_cycle(t%SIN,N)
4754
4755	do i=1,n
4756	call taylor_cycle(t%SIN,ii=i,value=value,j=jc)
4757
4758	call test_jc_spin(ns,jc,-clockwise,nd,doit)
4759
4760	if(doit) then
4761	ang=0.0_dp
4762	do j=1,nd
4763	ang=(jc(j2-1)-jc(j2))twopins%n%tune(j)+ang
4764	enddo
4765	denom=value/(exp(-i_ang)-exp(-clockwisei_*ns%theta0))
4766	omr(3)=omr(3)-(denom.mono.jc)
4767	endif
4768	enddo
4769	1111 continue
4770	call kill(omc,3)
4771	call alloc(omc,3)
4772
4773	omc(1)= (omr(3)+omr(1))/2.0_dp
4774	omc(3)= (omr(3)-omr(1))/2.0_dp/i_
4775	omc(2)=omr(2)
4776
4777	do i=1,3
4778	t%cos=omc(i)%r
4779	t%sin=omc(i)%i
4780	oma(i)=t
4781	enddo
4782
4783	if(doflip) then
4784	call flip_real_array(ns%n%tune,ns%n%tune,-1)
4785	call fliptaylor(om(1),om(1),-1)
4786	call fliptaylor(om(2),om(2),-1)
4787	call fliptaylor(om(3),om(3),-1)
4788	call fliptaylor(oma(1),oma(1),-1)
4789	call fliptaylor(oma(2),oma(2),-1)
4790	call fliptaylor(oma(3),oma(3),-1)
4791	if(use_ptc_ac_position) then
4792	call flip_resonance(ns%m,ns%m,-1)
4793	endif
4794	perform_flip=.true.
4795	endif
4796
4797	deallocate(jc)
4798	call kill(t)
4799	call kill(omr,3)
4800	call kill(omc,3)
4801
4802	end subroutine get_kernel
4803
4804	subroutine test_jc(ns,jc,nd,doit)
4805	implicit none
4806	logical doit
4807	integer i,nd,k,l,j
4808	integer jc(:)
4809	TYPE(normal_spin), INTENT(INout) :: ns
4810
4811	doit=.true.
4812
4813	k=0
4814	do i=1,nd
4815	k=iabs(jc(i2-1)-jc(i2))+k
4816	enddo
4817	if(k==0) doit=.false.
4818
4819	if(.not.doit) return
4820
4821	do j=1,ns%n%nres
4822	k=0
4823	l=0
4824	do i=1,nd
4825	k=iabs(jc(i2-1)-jc(i2)-ns%n%m(i,j))+k
4826	l=iabs(jc(i2-1)-jc(i2)+ns%n%m(i,j))+l
4827	enddo
4828	if(k==0.or.l==0) then
4829	doit=.false.
4830	exit
4831	endif
4832	enddo
4833
4834
4835	end subroutine test_jc
4836
4837	subroutine test_jc_spin(ns,jc,is,nd,doit)
4838	implicit none
4839	logical doit
4840	integer i,nd,k,l,j,is
4841	integer jc(:)
4842	TYPE(normal_spin), INTENT(INout) :: ns
4843
4844	doit=.true.
4845
4846	! write(6,*) jc(1:4),is
4847	! write(6,*) ns%m(1,1),ns%m(2,1),ns%ms(1)
4848
4849	do j=1,ns%nres
4850	k=0
4851	l=0
4852	do i=1,nd
4853	k=iabs(jc(i2-1)-jc(i2)-ns%m(i,j))+k
4854	l=iabs(jc(i2-1)-jc(i2)+ns%m(i,j))+l
4855	enddo
4856	k=k+iabs(is-ns%ms(j))
4857	l=l+iabs(is+ns%ms(j))
4858	if(k==0.or.l==0) then
4859	doit=.false.
4860	exit
4861	endif
4862	enddo
4863
4864
4865	end subroutine test_jc_spin
4866
4867	subroutine check_spin(DS,spin_in)
4868	implicit none
4869	logical(lp), INTENT(INout) :: spin_in
4870	TYPE(damapspin), INTENT(IN) :: DS
4871	integer i,j
4872	real(dp) norm
4873
4874	spin_in=.true.
4875	norm=0.0_dp
4876	do i=1,3
4877	do j=1,3
4878	norm=norm+full_abs(DS%s%s(i,j))
4879	enddo
4880	enddo
4881	norm=abs(norm-3.0_dp)
4882
4883	if(norm<=eps_tpsalie) then
4884	write(6,*) " Spin Map is identity : not normalized "
4885	spin_in=.false.
4886	endif
4887
4888	end subroutine check_spin
4889
4890	subroutine check_rad(e_ij,rad_in)
4891	implicit none
4892	logical(lp), INTENT(INout) :: rad_in
4893	real(dp) e_ij(6,6)
4894	integer i,j
4895	real(dp) norm
4896
4897	rad_in=.true.
4898	norm=0.0_dp
4899	do i=1,6
4900	do j=1,6
4901	norm=norm+abs(e_ij(i,j))
4902	enddo
4903	enddo
4904
4905	if(norm==0.0_dp) then
4906	if(global_verbose) write(6,*) " Radiation Envelope is 0.0_dp : not printed "
4907	rad_in=.false.
4908	endif
4909
4910	end subroutine check_rad
4911
4912	subroutine normalise_spin(ns,DS_in)
4913	implicit none
4914	TYPE(normal_spin), INTENT(INout) :: ns
4915	TYPE(damapspin), INTENT(IN) :: DS_in
4916	TYPE(damapspin) ds,ds0,dst,dsi
4917	TYPE(damap) r0
4918	type(taylor) om(3),oma(3)
4919	integer i
4920	logical(lp) spin_in
4921
4922	call alloc(ds)
4923	call alloc(dst)
4924	call alloc(dsi)
4925	call alloc(ds0)
4926	call alloc(r0)
4927	call alloc(om)
4928	call alloc(oma)
4929	spin_in=.false.
4930
4931	call check_spin(DS_in,spin_in)
4932	ds=ds_in
4933	! step 1
4934	! Normalize to the parameter dependent n0 of the theory
4935	call Go_to_closed(ns,DS,spin_in)
4936
4937	! step 2
4938	! Apply the transformation found in step 1 to ds
4939	ds=ns%as*(-1)ds*ns%as
4940	! Apply the transformation found in step 1 to ds
4941	ds=ns%ar*(-1)ds*ns%ar ! around fixed point and transversely normalised
4942
4943	call fetch_s0(DS,ds0) ! ds0=(I,exp(theta0 L_y) )
4944
4945
4946	ns%a_t=1
4947	r0=ns%n%normal
4948	if(spin_in) then
4949
4950	do i=1,c_%no+1
4951
4952	dst=dsds0(-1) ! ds=exp(small)ds0
4953	! dst=exp(small)
4954
4955	call dalog(DSt,om) ! small=log(dst) finite # of steps
4956
4957	call get_kernel(ns,om,oma)
4958
4959	call daexplog(oma,dsi) ! dsi=exp(oma(i)*L_i)
4960	ns%a_t=ns%a_tdsi ! updated a_t= a_tdsi
4961
4962	ds=dsdsi ! ds=dsi(-1) ds dsi
4963	dsi=dsi**(-1) ! eventually ds is normalized
4964	ds=dsi*ds
4965	enddo ! keep going
4966
4967	endif
4968	ns%a_t=ns%asns%arns%a_t
4969
4970
4971	dsi=ns%a_t
4972
4973	ns%a1=1
4974	ns%a1%m=ns%n%a1
4975
4976	dsi=ns%a1*(-1)ns%a_t
4977	ns%as=dsi
4978	ns%as%m=1
4979	ns%as=ns%asdsi%m*(-1)
4980	ns%ar=1
4981	ns%ar%m=dsi%m ! a_t= a1 o A_s o A_r
4982	call kill(ds)
4983	call kill(dst)
4984	call kill(dsi)
4985	call kill(ds0)
4986	call kill(r0)
4987	call kill(om)
4988	call kill(oma)
4989	end subroutine normalise_spin
4990
4991	subroutine eval_spin_matrix(a_in,x,a_out)
4992	implicit none
4993	TYPE(damapspin), INTENT(IN) :: a_in
4994	TYPE(damapspin), INTENT(out) :: a_out
4995	TYPE(damapspin) a_temp
4996	real(dp) x(lnv)
4997	integer i,j
4998
4999	call alloc(a_temp)
5000
5001	do i=c_%nv+1,lnv
5002	x(i)=0.0_dp
5003	enddo
5004
5005	a_temp=a_in
5006
5007	do i=1,3
5008	do j=1,3
5009	if(a_in%s%s(i,j)%kind==2) then
5010	a_temp%s%s(i,j)=a_in%s%s(i,j)%t*x
5011	endif
5012	enddo
5013	enddo
5014	a_out=a_temp
5015
5016	call kill(a_temp)
5017
5018	end subroutine eval_spin_matrix
5019
5020
5021	subroutine factor_as(a_t,a_f,a_s,a_l,a_nl,DR,R_TE,CS_TE,COSLIKE,s0,s_nl)
5022	implicit none
5023	TYPE(damapspin), INTENT(INout) :: a_t,a_f,a_s,a_l,a_nl
5024	logical(lp) lagrange0,factor_spin0
5025	TYPE(damapspin), optional, intent(inout) ::R_TE,CS_TE,DR,s0,s_nl
5026	logical(lp) , optional, intent(inout) :: COSLIKE
5027	TYPE(damapspin) rot_y,temp,tempi
5028
5029	factor_spin0=my_false
5030	if(present(s0).and.present(s_nl)) then
5031	factor_spin0=my_true
5032	elseif(.not.(present(s0)).and.(.not.present(s_nl))) then
5033	! do nothing
5034	else
5035	write(6,*) " error in factor_as "
5036	stop 333
5037	endif
5038	a_f=1
5039	a_l=1
5040	a_nl=1
5041	lagrange0=my_false
5042	if(present(dr)) then
5043	lagrange0=my_true
5044	dr=1
5045	endif
5046	if(present(R_TE).and.present(CS_TE)) then
5047	R_TE=1
5048	CS_TE=1
5049	elseif((.not.present(R_TE)).and.(.not.present(CS_TE))) then
5050	! nothing to be done
5051	else
5052	write(6,*) " factor_as has an error: R_TE and CS_TE must be both present or absent "
5053	stop 1068
5054
5055	endif
5056
5057	call factor(a_t%m,a_f%m,a_l%m,a_nl%m,DR%m,R_TE%m,CS_TE%m,COSLIKE)
5058
5059	if(lagrange0) then
5060	a_s=1
5061	a_s%s=a_t%s
5062	a_s=a_sdr*(-1)
5063
5064	call alloc(rot_y,temp,tempi)
5065	if(factor_spin0) then
5066	call remove_y_rot0(a_s,s0,s_nl,r_y=rot_y)
5067	else
5068	call remove_y_rot(a_s,r_y=rot_y)
5069	endif
5070	if(present(dr)) dr=dr*rot_y
5071	call kill(rot_y,temp,tempi)
5072
5073	a_t%s=a_s%s
5074	endif
5075
5076
5077
5078	a_s=a_f*(-1)a_t
5079	a_s%m=1
5080	a_s=a_sa_t%m*(-1)
5081
5082
5083
5084
5085	!!! this creates
5086
5087	!! (a_t%m,a_t%s) = (a_f%m, I ) o (I ,a_s%s) o (a_l%m,I) o (a_nl%m,I)
5088
5089	!
5090	!if(present(a_s0)) then
5091	!a_s0=1
5092	!call clean_orbital_33(a_s%s,a_s0%s)
5093
5094	!a_s=a_s0*(-1)a_s
5095	!endif
5096
5097
5098
5099	end subroutine factor_as
5100
5101	subroutine CANONIZE( a_t,A_cs,PHASE_ADVANCE,R_TE,CS_TE,COSLIKE )
5102	implicit none
5103	TYPE(damap), INTENT(INout) :: a_t,A_cs
5104	type(taylor),optional, INTENT(INout) :: PHASE_ADVANCE(:)
5105	TYPE(damap), optional, intent(inout) ::R_TE,CS_TE
5106	logical(lp) , optional, intent(inout) :: COSLIKE
5107	TYPE(damap) a_f,a_l,a_nl,dr1,a_tt
5108	type(onelieexponent) uno
5109	logical(lp) doflip
5110	integer i
5111
5112
5113	if(perform_flip.and.new_ndpt.and.c_%ndpt/=0) then
5114	! write(6,*) " flipping ",c_%ndpt,c_%nd2-1
5115	! pause 7123
5116
5117	perform_flip=.false.
5118	call flip_damap(a_t,a_t)
5119	doflip=.true.
5120	else
5121	doflip=.false.
5122	endif
5123
5124	call alloc(a_f,a_l,a_nl,dr1,a_tt)
5125	call alloc(uno)
5126	dr1=1
5127	a_tt=a_t
5128	call factor(a_tt,a_f,a_l,a_nl,DR1,R_TE,CS_TE,COSLIKE)
5129
5130	A_cs=a_fa_la_nl
5131	if(present(PHASE_ADVANCE)) then
5132	uno=dr1
5133
5134	if(c_%ndpt==0) then
5135	do i=1,c_%nd
5136	PHASE_ADVANCE(i)=PHASE_ADVANCE(i)+((uno%VECTOR%v(2i-1)).k.(2i))/twopi
5137	enddo
5138	else
5139	if(c_%ndpt>c_%nd2-2) then
5140	do i=1,c_%nd-1
5141	PHASE_ADVANCE(i)=PHASE_ADVANCE(i)+((uno%VECTOR%v(2i-1)).k.(2i))/twopi
5142	enddo
5143	else
5144	do i=1,c_%nd-2
5145	PHASE_ADVANCE(i)=PHASE_ADVANCE(i)+((uno%VECTOR%v(2i-1)).k.(2i))/twopi
5146	enddo
5147	i=c_%nd-1
5148	PHASE_ADVANCE(i+1)=PHASE_ADVANCE(i)+((uno%VECTOR%v(2i-1)).k.(2i))/twopi
5149	endif
5150	endif
5151
5152	endif
5153
5154
5155	if(doflip) then
5156	call flip_damap(a_t,a_t)
5157	call flip_damap(a_cs,a_cs)
5158	call flip_damap(dr1,dr1)
5159	if(present(PHASE_ADVANCE)) then
5160	do i=1,c_%nd
5161	call flip_taylor(PHASE_ADVANCE(i),PHASE_ADVANCE(i),-1)
5162	enddo
5163	endif
5164	perform_flip=.true.
5165	endif
5166	call kill(uno)
5167	call kill(a_f,a_l,a_nl,dr1,a_tt)
5168
5169	end subroutine CANONIZE
5170
5171	subroutine factor_am(a_t,a_f,a_l,a_nl,DR,R_TE,CS_TE,COSLIKE)
5172	implicit none
5173	TYPE(damap), INTENT(INout) :: a_t,a_nl,a_l,a_f
5174	integer i,n,k,nt,j
5175	integer, allocatable :: jc(:)
5176	real(dp) value,alpha0,sip
5177	logical doit
5178	TYPE(damap) atemp,s1,s1i,m1
5179	TYPE(taylor)a12,a11,p(ndim)
5180	logical(lp) doflip,dote,lagrange0
5181	TYPE(damap), optional, intent(inout) ::R_TE,CS_TE,DR
5182	logical(lp) , optional, intent(inout) :: COSLIKE
5183	type(taylor) m(ndim2,ndim2)
5184	type(taylor) at(2,2),bt(2,2),ct(2,2),dt(2,2),ati(2,2),bti(2,2),alpha,det
5185	type(gmap) g
5186	type(onelieexponent) un
5187	type(reversedragtfinn) rdf
5188	type(vecresonance) vr
5189	logical(lp) t_e
5190
5191	t_e=my_true
5192	lagrange0=my_false
5193	if(present(dr)) lagrange0=my_true
5194
5195	call alloc(atemp,s1,s1i,m1)
5196	call alloc(a12,a11)
5197	call alloc(p,ndim)
5198	call alloc_nn(m)
5199	call alloc_nn(at)
5200	call alloc_nn(bt)
5201	call alloc_nn(ct)
5202	call alloc_nn(dt)
5203	call alloc_nn(ati)
5204	call alloc_nn(bti)
5205	call alloc(alpha,det)
5206	allocate(jc(c_%nv))
5207
5208
5209
5210	if(perform_flip.and.new_ndpt.and.c_%ndpt/=0) then
5211	! write(6,*) " flipping ",c_%ndpt,c_%nd2-1
5212	! pause 7123
5213
5214	perform_flip=.false.
5215	call flip_damap(a_t,a_t)
5216	doflip=.true.
5217	else
5218	doflip=.false.
5219	endif
5220
5221	atemp=1
5222	do k=1,c_%nd2
5223	call taylor_cycle(a_t%v(k),N)
5224
5225	do i=1,n
5226	call taylor_cycle(a_t%v(k),ii=i,value=value,j=jc)
5227	call check_fix(jc,0,doit)
5228	if(doit) atemp%v(k)= atemp%v(k) + (value.mono.jc)
5229	enddo
5230	enddo
5231
5232	! in case we have ndpt/=0
5233
5234	if(c_%ndpt==c_%nd2-1) then ! ptc convention
5235	atemp%v(c_%nd2-1)=atemp%v(c_%nd2-1)-(1.0_dp.mono.(c_%nd2-1))
5236	atemp%v(c_%nd2)=atemp%v(c_%nd2)-(1.0_dp.mono.(c_%nd2))
5237	do i=1,c_%nd-1
5238	atemp%v(c_%ndpt+1)=atemp%v(c_%ndpt+1) &
5239	+(atemp%v(2i).d.c_%ndpt)(1.0_dp.mono.2i-1)-(atemp%v(2i-1).d.c_%ndpt)(1.0_dp.mono.2i)
5240	enddo
5241	elseif(c_%ndpt==c_%nd2) then ! Marylie convention
5242	atemp%v(c_%nd2-1)=atemp%v(c_%nd2-1)-(1.0_dp.mono.(c_%nd2-1))
5243	atemp%v(c_%nd2)=atemp%v(c_%nd2)-(1.0_dp.mono.(c_%nd2))
5244	do i=1,c_%nd-1
5245	atemp%v(c_%ndpt-1)=atemp%v(c_%ndpt-1) &
5246	-(atemp%v(2i).d.c_%ndpt)(1.0_dp.mono.2i-1)+(atemp%v(2i-1).d.c_%ndpt)(1.0_dp.mono.2i)
5247	enddo
5248	endif
5249
5250
5251	a_f=atemp !!! not exact ! temporal part could be wrong!!!
5252
5253	a_l=a_f*(-1)a_t
5254
5255	atemp=0
5256	do k=1,c_%nd2
5257	call taylor_cycle(a_l%v(k),N)
5258
5259	do i=1,n
5260	call taylor_cycle(a_l%v(k),ii=i,value=value,j=jc)
5261	call check_fix(jc,1,doit)
5262	if(doit) atemp%v(k)= atemp%v(k) + (value.mono.jc)
5263	enddo
5264	enddo
5265
5266
5267	! in case we have ndpt/=0
5268
5269	if(c_%ndpt/=0) then ! ptc convention
5270	atemp%v(c_%nd2-1)=(1.0_dp.mono.(c_%nd2-1))
5271	atemp%v(c_%nd2)= (1.0_dp.mono.(c_%nd2))
5272
5273	if(c_%ndpt==c_%nd2) then
5274	k=c_%nd2-1
5275	else
5276	k=c_%nd2
5277	endif
5278
5279	call taylor_cycle(a_l%v(k),N)
5280
5281	do i=1,n
5282	call taylor_cycle(a_l%v(k),ii=i,value=value,j=jc)
5283	call check_fix(jc,2,doit)
5284	if(doit) atemp%v(k)= atemp%v(k) + (value.mono.jc)
5285	enddo
5286
5287	endif
5288
5289	a_nl=atemp*(-1)a_l
5290	a_l=atemp
5291
5292	!!! a_t=a_fa_la_nl at this stage
5293
5294	! if(present(lagrange)) then
5295	if(lagrange0) then
5296	! enforce Teng-Edwards to all orders in parameters i.e. A_12=0 and A_3_4=0 etc....
5297	nt=c_%nd
5298	if(C_%ndpt/=0) nt=nt-1
5299
5300	s1=1
5301	s1i=1
5302
5303	if(C_%ndpt/=0) then
5304	s1%v(C_%ndpt)=1.0_dp.mono.C_%ndpt
5305	s1i%v(C_%ndpt)=1.0_dp.mono.C_%ndpt
5306	endif
5307
5308	do i=1,nt
5309	a11=a_l%v(2i-1).d.(2i-1)
5310	a12=a_l%v(2i-1).d.(2i)
5311	p(i)=-a12/a11
5312	p(i)=atan(p(i))
5313	sip=a11cos(p(i))-a12sin(p(i))
5314	if(sip<0) p(i)=p(i)+pi
5315	if(C_%ndpt/=0) then
5316	if(mod(C_%ndpt,2)==1) then
5317	s1%v(C_%ndpt+1)=s1%v(C_%ndpt+1)-(p(i).d.C_%ndpt)((1.0_dp.mono.(2i))*2+(1.0_dp.mono.(2i-1))*2)0.5_dp
5318	s1i%v(C_%ndpt+1)=s1i%v(C_%ndpt+1)+(p(i).d.C_%ndpt)((1.0_dp.mono.(2i))*2+(1.0_dp.mono.(2i-1))*2)0.5_dp
5319	else
5320	s1%v(C_%ndpt-1)=s1%v(C_%ndpt-1)+(p(i).d.C_%ndpt)((1.0_dp.mono.(2i))*2+(1.0_dp.mono.(2i-1))*2)0.5_dp
5321	s1i%v(C_%ndpt-1)=s1i%v(C_%ndpt-1)-(p(i).d.C_%ndpt)((1.0_dp.mono.(2i))*2+(1.0_dp.mono.(2i-1))*2)0.5_dp
5322	endif
5323	endif
5324	s1%v(2i-1)=COS(p(i))(1.0_dp.mono.(2i-1))+SIN(p(i))(1.0_dp.mono.(2*i))
5325	s1%v(2i) =COS(p(i))(1.0_dp.mono.(2i))-SIN(p(i))(1.0_dp.mono.(2*i-1))
5326	s1i%v(2i-1)=COS(p(i))(1.0_dp.mono.(2i-1))-SIN(p(i))(1.0_dp.mono.(2*i))
5327	s1i%v(2i) =COS(p(i))(1.0_dp.mono.(2i))+SIN(p(i))(1.0_dp.mono.(2*i-1))
5328	enddo
5329	if(.not.courant_snyder) then
5330	s1i=1
5331	s1=1
5332	endif
5333	a_nl=s1ia_nls1
5334	a_l=a_l*s1
5335	dr=s1i
5336	! a_t_cs= a_fa_l_cs a_nl ! a_nl is not yet " standard "
5337	!!!!! nonlinear part if s1i !!!!!!
5338	if(c_%no>1.and.courant_snyder) then
5339	if(onelie) then
5340	nd_used=c_%nd2/2
5341	if(c_%ndpt/=0) nd_used=nd_used-1
5342	call alloc(vr)
5343	call alloc(un)
5344
5345	s1i=a_nl ! preserve it
5346	s1=1
5347	do i=1,c_%no ! i=1,c_%no
5348	un%eps=-c_%no
5349	un=s1i
5350	vr=un%vector
5351	decal=1 ! not coasting plane
5352	do k=1,nd_used*2
5353	i_phase=k
5354	if(mod(k,2)==0) then
5355	i_plane=k/2
5356	else
5357	i_plane=(k+1)/2
5358	endif
5359	call cfu(vr%cos%v(k),dphase,vr%cos%v(k))
5360	call cfu(vr%sin%v(k),dphase,vr%sin%v(k))
5361	enddo
5362
5363	do k=nd_used*2+1,c_%nd2 ! fix coasting later
5364	vr%cos%v(k)=0.0_dp
5365	vr%sin%v(k)=0.0_dp
5366	enddo
5367
5368	un%vector=vr
5369	s1=texp(un%vector,s1)
5370	s1i=texp(un%vector,s1i)
5371	enddo ! i=1,c_%no
5372
5373	!!!! now fix coasting !!!!!
5374	if( nd_used == c_%nd) then
5375	a_nl=s1i
5376	else ! coasting
5377	un%eps=-c_%no
5378	un=s1
5379	alpha0=1.0_dp
5380	if(mod(c_%ndpt,2)==0) alpha0=-1.0_dp
5381	do i=1,c_%nd2-2
5382	un%vector%v(i)=alpha0*(un%vector%v(i).d.c_%ndpt)
5383	enddo
5384	call int_partial(un%vector,un%pb,2)
5385
5386	if(mod(c_%ndpt,2)==0) then
5387	s1%v(c_%ndpt-1)=s1%v(c_%ndpt-1)+un%pb%h
5388	else
5389	s1%v(c_%ndpt+1)=s1%v(c_%ndpt+1)+un%pb%h
5390	endif
5391
5392	a_nl=a_nl*s1
5393	endif ! end of nd_used == c_%nd
5394
5395	!!!! end of now fix coasting !!!!!
5396
5397	dr=drs1*(-1)
5398	call kill(vr)
5399	call kill(un)
5400	else ! reverse Dragt-Finn
5401	nd_used=c_%nd2/2
5402	if(c_%ndpt/=0) nd_used=nd_used-1
5403	call alloc(vr)
5404	call alloc(rdf)
5405
5406	s1i=a_nl ! preserve it
5407	s1=1
5408	do i=1,c_%no ! i=1,c_%no
5409	rdf=s1i
5410	vr=rdf%nonlinear
5411	decal=1 ! not coasting plane
5412	do k=1,nd_used*2
5413	i_phase=k
5414	if(mod(k,2)==0) then
5415	i_plane=k/2
5416	else
5417	i_plane=(k+1)/2
5418	endif
5419	call cfu(vr%cos%v(k),dphase,vr%cos%v(k))
5420	call cfu(vr%sin%v(k),dphase,vr%sin%v(k))
5421	enddo
5422
5423	do k=nd_used*2+1,c_%nd2 ! fix coasting later
5424	vr%cos%v(k)=0.0_dp
5425	vr%sin%v(k)=0.0_dp
5426	enddo
5427
5428	rdf%nonlinear=vr
5429	s1=texp(rdf%nonlinear,s1)
5430	s1i=texp(rdf%nonlinear,s1i)
5431	enddo ! i=1,c_%no
5432
5433	!!!! now fix coasting !!!!!
5434	if( nd_used == c_%nd) then
5435	a_nl=s1i
5436	else ! coasting
5437	rdf=s1
5438	alpha0=1.0_dp
5439	if(mod(c_%ndpt,2)==0) alpha0=-1.0_dp
5440	do i=1,c_%nd2-2
5441	rdf%nonlinear%v(i)=alpha0*(rdf%nonlinear%v(i).d.c_%ndpt)
5442	enddo
5443	call int_partial(rdf%nonlinear,rdf%pb,2)
5444
5445	if(mod(c_%ndpt,2)==0) then
5446	s1%v(c_%ndpt-1)=s1%v(c_%ndpt-1)+rdf%pb%h
5447	else
5448	s1%v(c_%ndpt+1)=s1%v(c_%ndpt+1)+rdf%pb%h
5449	endif
5450
5451	a_nl=a_nl*s1
5452	endif ! end of nd_used == c_%nd
5453
5454	!!!! end of now fix coasting !!!!!
5455
5456	dr=drs1*(-1)
5457	call kill(vr)
5458	call kill(rdf)
5459	endif
5460	endif ! no>1
5461
5462
5463	!!!!
5464	endif ! it (te)
5465	! endif ! present(te)
5466
5467
5468	dote=present(R_TE).and.present(CS_TE)
5469	if(doing_ac_modulation_in_ptc) then
5470	dote=dote.and.(c_%nd2==6.or.c_%ndpt>=5)
5471	else
5472	dote=dote.and.(c_%nd2==4.or.c_%ndpt>=5)
5473	endif
5474	if(dote) then
5475
5476	call copy_damap_matrix(a_l,m)
5477	call copy_matrix_matrix(m(1:2,1:2),at)
5478	call copy_matrix_matrix(m(1:2,3:4),ct)
5479	call copy_matrix_matrix(m(3:4,1:2),dt)
5480	call copy_matrix_matrix(m(3:4,3:4),bt)
5481
5482	call invert_22(at,ati)
5483	call invert_22(bt,bti)
5484	if(.not.c_%STABLE_DA) then
5485	t_e=my_false
5486	endif
5487
5488
5489	call matmul_nn(dt,ati,ati,sc=-1.0_dp)
5490	call matmul_nn(ati,ct,ct)
5491	call matmul_nn(ct,bti,ct)
5492	if(.not.c_%STABLE_DA) then
5493	t_e=my_false
5494	goto 888
5495	endif
5496
5497
5498	alpha=ct(1,1)
5499	alpha0=alpha
5500
5501	if(alpha0<=-1.0_dp) then
5502	t_e=my_false
5503	goto 888
5504	endif
5505
5506	det=sqrt(1.0_dp/(1.0_dp+alpha))
5507
5508
5509
5510
5511	if(alpha0>=0.0_dp) then
5512	COSLIKE=my_true
5513	else
5514	! det=sqrt(one/(one-alpha))
5515	COSLIKE=my_false
5516	endif
5517
5518	CS_TE=0
5519	do i=1,2
5520	do j=1,2
5521	CS_TE%v(i)=at(i,j)*(1.0_dp.mono.j)/det+CS_TE%v(i)
5522	CS_TE%v(i+2)=bt(i,j)*(1.0_dp.mono.(j+2))/det+CS_TE%v(i+2)
5523	enddo
5524	enddo
5525
5526
5527
5528	! The rotation matrix is created but it may not have the correct path length
5529	!dependence
5530	if(c_%ndpt/=0.and.t_e) then
5531	call alloc(g,c_%nv)
5532	call alloc(un)
5533	! write(6,*) " epseone "
5534	! read(5,*) un%eps
5535	un%eps=-c_%no
5536	do i=1,c_%nv
5537	g%v(i)=1.0_dp.mono.i
5538	enddo
5539	g%v(c_%ndpt)=0.0_dp
5540
5541	do i=1,c_%nd2
5542	s1%v(i) = CS_TE%v(i)*g
5543	enddo
5544	! s1%v(c_%ndpt)=one.mono.c_%ndpt
5545	!!!!!new
5546	s1%v(c_%nd2)=1.0_dp.mono.c_%nd2
5547	s1%v(c_%nd2-1)=1.0_dp.mono.c_%nd2-1
5548
5549
5550	CS_TE%v(c_%nd2)=1.0_dp.mono.c_%nd2
5551	CS_TE%v(c_%nd2-1)=1.0_dp.mono.c_%nd2-1
5552	!!!!!
5553	s1i=s1**(-1)
5554	s1i=s1i*CS_TE ! s1i is completely nonlinear.
5555	un=s1i
5556
5557	alpha0=1.0_dp
5558	if(mod(c_%ndpt,2)==0) alpha0=-1.0_dp
5559	do i=1,c_%nd2-2
5560	un%vector%v(i)=alpha0*(un%vector%v(i).d.c_%ndpt)
5561	enddo
5562	call int_partial(un%vector,un%pb,2)
5563	! un%pb=un%vector ! this is the longitudinal part
5564	if(mod(c_%ndpt,2)==0) then
5565	s1i%v(c_%ndpt-1)=s1i%v(c_%ndpt-1)+un%pb%h
5566	else
5567	s1i%v(c_%ndpt+1)=s1i%v(c_%ndpt+1)+un%pb%h
5568	endif
5569
5570
5571	CS_TE=s1*s1i
5572
5573
5574
5575	call kill(un)
5576	call kill(g)
5577	!endif !eps_te
5578	endif
5579
5580	888 continue
5581	if(.not.t_e) then
5582	c_%STABLE_DA=my_true
5583	cs_te=0
5584	R_TE=0
5585	write(6,*) " Teng-Edwards is crap !"
5586	else
5587	R_TE=a_lcs_TE*(-1)
5588	endif
5589
5590	endif ! end of T-E done
5591
5592	if(lagrange0) then
5593	a_t=a_fa_la_nl
5594	endif
5595
5596	if(doing_ac_modulation_in_ptc.and.present(CS_TE)) then ! removing useless tiny numbers
5597	CS_TE%v(c_%nd2-1)=1.0_dp.mono.(c_%nd2-1)
5598	CS_TE%v(c_%nd2)=1.0_dp.mono.(c_%nd2)
5599	endif
5600
5601	if(doflip) then
5602	call flip_damap(a_t,a_t)
5603	call flip_damap(a_nl,a_nl)
5604	call flip_damap(a_l,a_l)
5605	call flip_damap(a_f,a_f)
5606	if(present(dr)) call flip_damap(dr,dr)
5607	if(present(CS_TE)) call flip_damap(CS_TE,CS_TE)
5608	if(present(R_TE)) call flip_damap(R_TE,R_TE)
5609	perform_flip=.true.
5610	endif
5611
5612	deallocate(jc)
5613	call kill(atemp,s1,s1i,m1)
5614	call kill(a12,a11)
5615	call kill(p,ndim)
5616	call kill_nn(m)
5617	call kill_nn(at)
5618	call kill_nn(bt)
5619	call kill_nn(ct)
5620	call kill_nn(dt)
5621	call kill_nn(ati)
5622	call kill_nn(bti)
5623	call kill(alpha,det)
5624
5625
5626	end subroutine factor_am
5627
5628	subroutine factor_am_special(a_t,a_f,a_l,a_nl,DR)
5629	implicit none
5630	TYPE(damap), INTENT(INout) :: a_t,a_nl,a_l,a_f
5631	integer i,n,k,nt,j
5632	integer, allocatable :: jc(:)
5633	real(dp) value,alpha0
5634	logical doit
5635	TYPE(damap) atemp,s1,s1i,m1
5636	TYPE(taylor)a12,a11,p(ndim)
5637	logical(lp) doflip,dote,lagrange0
5638	TYPE(damap), optional, intent(inout) ::DR
5639	! type(taylor) m(ndim2,ndim2)
5640	! type(taylor) at(2,2),bt(2,2),ct(2,2),dt(2,2),ati(2,2),bti(2,2),alpha,det
5641	! type(gmap) g
5642	! type(onelieexponent) un
5643	! type(reversedragtfinn) rdf
5644	! type(vecresonance) vr
5645
5646
5647	lagrange0=my_false
5648	if(present(dr)) lagrange0=my_true
5649
5650	call alloc(atemp,s1,s1i,m1)
5651	call alloc(a12,a11)
5652	call alloc(p,ndim)
5653	! call alloc_nn(m)
5654	! call alloc_nn(at)
5655	! call alloc_nn(bt)
5656	! call alloc_nn(ct)
5657	! call alloc_nn(dt)
5658	! call alloc_nn(ati)
5659	! call alloc_nn(bti)
5660	! call alloc(alpha,det)
5661	allocate(jc(c_%nv))
5662
5663
5664
5665	if(perform_flip.and.new_ndpt.and.c_%ndpt/=0) then
5666	! write(6,*) " flipping ",c_%ndpt,c_%nd2-1
5667	! pause 7123
5668
5669	perform_flip=.false.
5670	call flip_damap(a_t,a_t)
5671	doflip=.true.
5672	else
5673	doflip=.false.
5674	endif
5675
5676	atemp=1
5677	do k=1,c_%nd2
5678	call taylor_cycle(a_t%v(k),N)
5679
5680	do i=1,n
5681	call taylor_cycle(a_t%v(k),ii=i,value=value,j=jc)
5682	call check_fix(jc,0,doit)
5683	if(doit) atemp%v(k)= atemp%v(k) + (value.mono.jc)
5684	enddo
5685	enddo
5686
5687	! in case we have ndpt/=0
5688
5689	if(c_%ndpt==c_%nd2-1) then ! ptc convention
5690	atemp%v(c_%nd2-1)=atemp%v(c_%nd2-1)-(1.0_dp.mono.(c_%nd2-1))
5691	atemp%v(c_%nd2)=atemp%v(c_%nd2)-(1.0_dp.mono.(c_%nd2))
5692	do i=1,c_%nd-1
5693	atemp%v(c_%ndpt+1)=atemp%v(c_%ndpt+1) &
5694	+(atemp%v(2i).d.c_%ndpt)(1.0_dp.mono.2i-1)-(atemp%v(2i-1).d.c_%ndpt)(1.0_dp.mono.2i)
5695	enddo
5696	elseif(c_%ndpt==c_%nd2) then ! Marylie convention
5697	atemp%v(c_%nd2-1)=atemp%v(c_%nd2-1)-(1.0_dp.mono.(c_%nd2-1))
5698	atemp%v(c_%nd2)=atemp%v(c_%nd2)-(1.0_dp.mono.(c_%nd2))
5699	do i=1,c_%nd-1
5700	atemp%v(c_%ndpt-1)=atemp%v(c_%ndpt-1) &
5701	-(atemp%v(2i).d.c_%ndpt)(1.0_dp.mono.2i-1)+(atemp%v(2i-1).d.c_%ndpt)(1.0_dp.mono.2i)
5702	enddo
5703	endif
5704
5705
5706	a_f=atemp !!! not exact ! temporal part could be wrong!!!
5707
5708	a_l=a_f*(-1)a_t
5709
5710	atemp=0
5711	do k=1,c_%nd2
5712	call taylor_cycle(a_l%v(k),N)
5713
5714	do i=1,n
5715	call taylor_cycle(a_l%v(k),ii=i,value=value,j=jc)
5716	call check_fix(jc,1,doit)
5717	if(doit) atemp%v(k)= atemp%v(k) + (value.mono.jc)
5718	enddo
5719	enddo
5720
5721
5722	! in case we have ndpt/=0
5723
5724	if(c_%ndpt/=0) then ! ptc convention
5725	atemp%v(c_%nd2-1)=(1.0_dp.mono.(c_%nd2-1))
5726	atemp%v(c_%nd2)= (1.0_dp.mono.(c_%nd2))
5727
5728	if(c_%ndpt==c_%nd2) then
5729	k=c_%nd2-1
5730	else
5731	k=c_%nd2
5732	endif
5733
5734	call taylor_cycle(a_l%v(k),N)
5735
5736	do i=1,n
5737	call taylor_cycle(a_l%v(k),ii=i,value=value,j=jc)
5738	call check_fix(jc,2,doit)
5739	if(doit) atemp%v(k)= atemp%v(k) + (value.mono.jc)
5740	enddo
5741
5742	endif
5743
5744	a_nl=atemp*(-1)a_l
5745	a_l=atemp
5746
5747	!!! a_t=a_fa_la_nl at this stage
5748
5749	! if(present(lagrange)) then
5750	if(lagrange0) then
5751	! enforce Teng-Edwards to all orders in parameters i.e. A_12=0 and A_3_4=0 etc....
5752	nt=c_%nd
5753	if(C_%ndpt/=0) nt=nt-1
5754
5755	s1=1
5756	s1i=1
5757
5758	if(C_%ndpt/=0) then
5759	s1%v(C_%ndpt)=1.0_dp.mono.C_%ndpt
5760	s1i%v(C_%ndpt)=1.0_dp.mono.C_%ndpt
5761	endif
5762
5763	do i=1,nt
5764	a11=a_l%v(2i-1).d.(2i-1)
5765	a12=a_l%v(2i-1).d.(2i)
5766	p(i)=-a12/a11
5767	p(i)=atan(p(i))
5768	if(C_%ndpt/=0) then
5769	if(mod(C_%ndpt,2)==1) then
5770	s1%v(C_%ndpt+1)=s1%v(C_%ndpt+1)-(p(i).d.C_%ndpt)((1.0_dp.mono.(2i))*2+(1.0_dp.mono.(2i-1))*2)0.5_dp
5771	s1i%v(C_%ndpt+1)=s1i%v(C_%ndpt+1)+(p(i).d.C_%ndpt)((1.0_dp.mono.(2i))*2+(1.0_dp.mono.(2i-1))*2)0.5_dp
5772	else
5773	s1%v(C_%ndpt-1)=s1%v(C_%ndpt-1)+(p(i).d.C_%ndpt)((1.0_dp.mono.(2i))*2+(1.0_dp.mono.(2i-1))*2)0.5_dp
5774	s1i%v(C_%ndpt-1)=s1i%v(C_%ndpt-1)-(p(i).d.C_%ndpt)((1.0_dp.mono.(2i))*2+(1.0_dp.mono.(2i-1))*2)0.5_dp
5775	endif
5776	endif
5777	s1%v(2i-1)=COS(p(i))(1.0_dp.mono.(2i-1))+SIN(p(i))(1.0_dp.mono.(2*i))
5778	s1%v(2i) =COS(p(i))(1.0_dp.mono.(2i))-SIN(p(i))(1.0_dp.mono.(2*i-1))
5779	s1i%v(2i-1)=COS(p(i))(1.0_dp.mono.(2i-1))-SIN(p(i))(1.0_dp.mono.(2*i))
5780	s1i%v(2i) =COS(p(i))(1.0_dp.mono.(2i))+SIN(p(i))(1.0_dp.mono.(2*i-1))
5781	enddo
5782	if(.not.courant_snyder) then
5783	s1i=1
5784	s1=1
5785	endif
5786	a_nl=s1ia_nls1
5787	a_l=a_l*s1
5788	dr=s1i
5789	! a_t_cs= a_fa_l_cs a_nl ! a_nl is not yet " standard "
5790	!!!!! nonlinear part if s1i !!!!!!
5791
5792
5793
5794	!!!!
5795	endif ! it (te)
5796
5797
5798
5799
5800	if(lagrange0) then
5801	a_t=a_fa_la_nl
5802	endif
5803
5804
5805	if(doflip) then
5806	call flip_damap(a_t,a_t)
5807	call flip_damap(a_nl,a_nl)
5808	call flip_damap(a_l,a_l)
5809	call flip_damap(a_f,a_f)
5810	if(present(dr)) call flip_damap(dr,dr)
5811
5812	perform_flip=.true.
5813	endif
5814
5815	deallocate(jc)
5816	call kill(atemp,s1,s1i,m1)
5817	call kill(a12,a11)
5818	call kill(p,ndim)
5819	! call kill_nn(m)
5820	! call kill_nn(at)
5821	! call kill_nn(bt)
5822	! call kill_nn(ct)
5823	! call kill_nn(dt)
5824	! call kill_nn(ati)
5825	! call kill_nn(bti)
5826	! call kill(alpha,det)
5827
5828
5829	end subroutine factor_am_special
5830
5831	subroutine int_partial(v,h,nd0)
5832	implicit none
5833	! IF SCA=-one
5834	! \VEC{V}.GRAD = J GRAD H . GRAD = :H:
5835	!
5836	! IF SCA=one
5837	! \VEC{V}.GRAD = GRAD H . GRAD
5838	integer i,nd0
5839	type(vecfield) v
5840	type(pbfield) h
5841	type(taylor) b4,b3,b2,b1
5842	type(damap) x
5843	logical doflip
5844	if(.not.c_%stable_da) return
5845
5846	if(perform_flip.and.new_ndpt.and.c_%ndpt/=0) then
5847	! write(6,*) " flipping ",c_%ndpt,c_%nd2-1
5848	! pause 7123
5849
5850	perform_flip=.false.
5851	call flip_vecfield(v,v,1)
5852	call flip_taylor(h%h,h%h,1)
5853	doflip=.true.
5854	else
5855	doflip=.false.
5856	endif
5857
5858	nd_used=nd0
5859	call alloc(x)
5860	call alloc(b4,b3,b2,b1)
5861
5862	x=1
5863
5864	do i=1,nd_used
5865	call cfu(v%v(2*i-1),dlie,b3)
5866	call cfu(v%v(2*i),dlie,b1)
5867	b2=b1x%v(2i-1)
5868	b1=b3x%v(2i)
5869	b3=b2-b1
5870	b2=b3+b4
5871	b4=b2
5872	enddo
5873	h%h=b4
5874
5875
5876	call kill(b4,b3,b2,b1)
5877	call kill(x)
5878
5879	if(doflip) then
5880	call flip_vecfield(v,v,-1)
5881	call flip_taylor(h%h,h%h,-1)
5882	perform_flip=.true.
5883	endif
5884	end subroutine int_partial
5885
5886	real(dp) function dlie(j)
5887	implicit none
5888	integer i
5889	! INTEGER J(NTT)
5890	integer,dimension(:)::j
5891	if(.not.c_%stable_da) return
5892
5893	dlie=0.0_dp
5894	do i=1,nd_used
5895	dlie=REAL(j(2i-1)+j(2i),kind=DP)+dlie
5896	enddo
5897	dlie=dlie+1.0_dp
5898	dlie=1.0_dp/dlie
5899	return
5900	end function dlie
5901
5902
5903	real(dp) function dphase(j)
5904	implicit none
5905	integer i
5906	! INTEGER J(NTT)
5907	integer,dimension(:)::j
5908	integer t,tu
5909	if(.not.c_%stable_da) return
5910
5911	t=-decal
5912	tu=0
5913	dphase=0.0_dp
5914	do i=1,nd_used
5915	t=abs(j(2i-1)-j(2i))+t
5916	enddo
5917	if(t==0.and.decal/=0) then
5918	tu=(j(2i_plane-1)-j(2i_plane))
5919	if(mod(i_phase,2)==1) then
5920	tu=tu-1
5921	else
5922	tu=tu+1
5923	endif
5924	if(tu==0) dphase=-1.0_dp
5925	else
5926	if(t==0) dphase=-1.0_dp
5927	endif
5928
5929	return
5930	end function dphase
5931
5932	real(dp) function phase_shift(j)
5933	implicit none
5934	integer i
5935	! INTEGER J(NTT)
5936	integer,dimension(:)::j
5937	integer nd,t
5938	if(.not.c_%stable_da) return
5939	nd=c_%nd2/2
5940	if(c_%ndpt/=0) nd=nd-1
5941
5942
5943	phase_shift=0.0_dp
5944	t=0
5945	do i=1,nd
5946	t=abs(j(2i-1)-j(2i))+t
5947	enddo
5948	if(t==0) phase_shift=1.0_dp
5949
5950	return
5951	end function phase_shift
5952
5953
5954	subroutine invert_22(a,ai)
5955	implicit none
5956	type(taylor) a(2,2),ai(2,2),t(2,2)
5957	type(taylor) det
5958	call alloc_nn(t)
5959	call alloc(det)
5960
5961	det=a(1,1)a(2,2)-a(1,2)a(2,1)
5962	t(1,1)=a(2,2)/det
5963	t(2,2)=a(1,1)/det
5964	t(1,2)=-a(1,2)/det
5965	t(2,1)=-a(2,1)/det
5966
5967	ai(1,1)=t(1,1)
5968	ai(1,1)=t(1,1)
5969	ai(1,1)=t(1,1)
5970	ai(1,1)=t(1,1)
5971
5972	call copy_matrix_matrix(t,ai)
5973
5974	call kill_nn(t)
5975	call kill(det)
5976
5977	end subroutine invert_22
5978
5979
5980	subroutine check_fix(jc,ord,doit)
5981	implicit none
5982	logical doit
5983	integer i,nd2,l,ord
5984	integer jc(:)
5985
5986	nd2=c_%nd2
5987	if(c_%ndpt/=0) nd2=c_%nd2-2
5988
5989	doit=.false.
5990
5991
5992	l=0
5993	do i=1,nd2
5994	l=jc(i)+l
5995	enddo
5996	if(l==ord) doit=.true.
5997
5998	end subroutine check_fix
5999
6000
6001	subroutine INTO_RES_SPIN8_eq(H_RES,H)
6002	implicit none
6003	TYPE(spinor_8), INTENT(in) :: h
6004	TYPE(res_spinor_8), INTENT(inout) :: H_RES
6005	call INTO_RES_SPIN8(H,H_RES)
6006	end subroutine INTO_RES_SPIN8_eq
6007
6008	subroutine INTO_RES_SPIN8(H,H_RES)
6009	implicit none
6010	TYPE(spinor_8), INTENT(in) :: h
6011	TYPE(res_spinor_8), INTENT(inout) :: H_RES
6012	type(taylorresonance) tr,ti
6013
6014	call alloc( tr)
6015	call alloc( ti)
6016
6017
6018	H_RES%X(1)= H%X(1)-i_*H%X(3)
6019	H_RES%X(2)= H%X(1)+i_*H%X(3)
6020	H_RES%X(3)= H%X(2)
6021
6022	tr=H_RES%X(3)%t%r
6023
6024	H_RES%X(3)= tr%cos + i_*tr%sin
6025
6026	tr=H_RES%X(1)%t%r
6027	ti=H_RES%X(1)%t%i
6028
6029	H_RES%X(1)=tr%cos + i_tr%sin + i_ (ti%cos + i_*ti%sin)
6030
6031	tr=H_RES%X(2)%t%r
6032	ti=H_RES%X(2)%t%i
6033
6034	H_RES%X(2)=tr%cos + i_tr%sin + i_ (ti%cos + i_*ti%sin)
6035
6036
6037	call kill( tr)
6038	call kill( ti)
6039
6040
6041	END SUBROUTINE INTO_RES_SPIN8
6042
6043	subroutine INTO_SPIN8_from_RES_eq(H,H_RES)
6044	implicit none
6045	TYPE(spinor_8), INTENT(inout) :: h
6046	TYPE(res_spinor_8), INTENT(in) :: H_RES
6047	type(taylorresonance) tr
6048	type(taylor) t
6049	type(complextaylor) c
6050
6051	call alloc( tr)
6052	call alloc( t)
6053	call alloc( c)
6054
6055	tr%cos=H_RES%X(3)%t%r
6056	tr%sin=H_RES%X(3)%t%i
6057	t=tr
6058	H%X(2)=morph(t)
6059	c=(H_RES%X(1)+H_RES%X(2))/2.0_dp
6060
6061	tr%cos=c%r
6062	tr%sin=c%i
6063	t=tr
6064	H%X(1)=morph(t)
6065
6066	c=i_*(H_RES%X(1)-H_RES%X(2))/2.0_dp
6067
6068	tr%cos=c%r
6069	tr%sin=c%i
6070	t=tr
6071	H%X(3)=morph(t)
6072
6073	call kill( c)
6074	call kill( t)
6075	call kill( tr)
6076
6077
6078	END SUBROUTINE INTO_SPIN8_from_RES_eq
6079
6080	!!!!!!!!!! normal form on theta into theta(H) !!!!!!!!!!!!!!!!
6081
6082	subroutine normal_thetaH(ds,ns)
6083	implicit none
6084	TYPE(damapspin), INTENT(INout) :: ds
6085	TYPE(normal_spin), INTENT(INout) :: ns
6086	TYPE(onelieexponent) uno
6087	TYPE(pbfield) h
6088	TYPE(damapspin) nc
6089	TYPE(normalform) nf
6090	TYPE(taylorresonance) tr,theta0r
6091	TYPE(taylor) theta0
6092	TYPE(complextaylor) gam
6093	real(dp) muc,epsy,ath,vvb
6094	complex(dp) dely,bth,beta,alpha
6095
6096
6097	muc=0.25e0_dp*twopi
6098	call alloc(h)
6099	call alloc(uno)
6100	call alloc(nc)
6101	call alloc(tr)
6102	call alloc(theta0r)
6103	call alloc(theta0)
6104	call alloc(gam)
6105	call alloc(nf)
6106
6107	nc=ns%a_t*(-1)ds*ns%a_t
6108
6109	h%h=muc*((1.0_dp.mono.'002')+(1.0_dp.mono.'0002'))/2.0_dp
6110	h%h=h%h+muc*((1.0_dp.mono.'2')+(1.0_dp.mono.'02'))/2.0_dp
6111
6112	nc%m=texp(h,nc%m)
6113
6114
6115	! call print(nc%m,6)
6116	! nf=nc%m
6117
6118	uno=nc%m
6119	uno%pb%h=uno%pb%h-muc*((1.0_dp.mono.'2')+(1.0_dp.mono.'02'))/2.0_dp
6120	tr=uno%pb%h
6121
6122	call print(ns%theta0,6)
6123
6124	! a11=(s(1,1))
6125	! a13=(s(1,3))
6126
6127	theta0=atan2(nc%s%s(1,3),nc%s%s(1,1))
6128	if((theta0.sub.'0')<0.0_dp) theta0 = theta0 + twopi
6129
6130	theta0r=theta0
6131
6132	!! call taylor_clean(theta0r%cos,1.d-1)
6133	! call taylor_clean(theta0r%sin,1.d-1)
6134	! call taylor_clean(tr%cos,1.d-5)
6135	! call taylor_clean(tr%sin,1.d-5)
6136
6137	call print(theta0r%cos,6)
6138	call print(tr%cos,6)
6139
6140	call print(theta0r%sin,6)
6141	call print(tr%sin,6)
6142
6143	epsy=tr%cos.sub.'0011'
6144	dely=(tr%cos.sub.'0040')+i_*(tr%sin.sub.'0040')
6145	ath=theta0r%cos.sub.'0011'
6146	bth=(theta0r%cos.sub.'0040')+i_*(theta0r%sin.sub.'0040')
6147	write(6,*) epsy
6148	write(6,*) dely
6149	write(6,*) ath
6150	write(6,*) bth
6151	vvb=(0.2d-3)**2
6152
6153	alpha=(athdely/epsy-bth)i_/8.d0/dely
6154	write(6,*) alpha
6155	alpha=abs(alpha)
6156	write(6,*) alpha
6157	alpha=0
6158	beta=((athdely/epsy-bth)i_/8.d0-dely*alpha)/epsy
6159	beta=abs(beta)
6160	write(6,*) beta
6161	alpha=0.1d0
6162	beta=((athdely/epsy-bth)i_/8.d0-dely*alpha)/epsy
6163	beta=abs(beta)
6164	write(6,*) beta
6165	! alpha=beta(dely/epsy)vvb
6166
6167
6168	call kill(uno)
6169	call kill(h)
6170	call kill(nc)
6171	call kill(tr)
6172	call kill(theta0r)
6173	call kill(theta0)
6174	call kill(nf)
6175	call kill(gam)
6176
6177	end subroutine normal_thetaH
6178
6179
6180	!!! Some useful routines
6181
6182	subroutine AVERAGE(F,A,F_floquet,F_xp,use_J)
6183	implicit none
6184	type(damap) A
6185	TYPE(TAYLOR), intent(inout):: F
6186	TYPE(TAYLOR), intent(inout):: F_floquet
6187	TYPE(TAYLOR), optional, intent(inout):: F_xp
6188	type(taylor) tt,tt_xp
6189	TYPE(taylorresonance) fq
6190	real(dp) value,valuexp
6191	integer, allocatable :: jc(:)
6192	logical, optional :: use_J
6193	logical usej
6194	integer i,n,j,it,nd,iu
6195	logical doflip,uj
6196
6197	if(perform_flip.and.new_ndpt.and.c_%ndpt/=0) then
6198	perform_flip=.false.
6199	call fliptaylor(F,F,1)
6200	call flip_damap(A,A)
6201	doflip=.true.
6202	else
6203	doflip=.false.
6204	endif
6205
6206	uj=.false.
6207	if(present(use_j)) uj=use_j
6208
6209	nd=c_%nd2/2
6210	if(c_%ndpt/=0) nd=nd-1
6211
6212	allocate(jc(c_%nv))
6213	call alloc(tt,tt_xp)
6214	call alloc(fq)
6215	!! USE_J is false by default
6216
6217
6218	fq=F*A
6219	fq%sin=0.0_dp
6220
6221	call taylor_cycle(fq%cos,size=n)
6222
6223	do i=1,n
6224	call taylor_cycle(fq%cos,ii=i,value=value,j=jc)
6225
6226	it=0
6227	iu=0
6228	do j=1,nd
6229	it=iabs(jc(j2-1)-jc(j2))+it
6230	iu=iabs(jc(j2-1)+jc(j2))+iu
6231	enddo
6232	if(it==0) then
6233	iu=iu/2
6234	valuexp=value
6235	if(uj) then
6236	value=valuexp2.0_dp*iu
6237	endif
6238	tt=(value.mono.jc)+tt
6239	tt_xp=(valuexp.mono.jc)+tt_xp
6240	endif
6241
6242	enddo
6243
6244	fq%cos=tt
6245	F_floquet=tt
6246
6247	if(present(F_xp)) then
6248	fq%cos=tt_xp
6249	F_xp=fq
6250	F_xp=F_xpA*(-1)
6251	endif
6252
6253	deallocate(jc)
6254	call kill(tt,tt_xp)
6255	call kill(fq)
6256
6257	if(doflip) then
6258	call fliptaylor(F,F,-1)
6259	call flip_damap(A,A)
6260	call fliptaylor(F_floquet,F_floquet,-1)
6261	if(present(F_xp)) call fliptaylor(F_xp,F_xp,-1)
6262	perform_flip=.true.
6263	endif
6264
6265	end subroutine AVERAGE
6266
6267	! remove small numbers
6268
6269	SUBROUTINE clean_damapspin(S1,S2,prec)
6270	implicit none
6271	type (damapspin),INTENT(INOUT)::S2
6272	type (damapspin), intent(INOUT):: s1
6273	real(dp) prec
6274	integer i,j
6275
6276	call clean_damap(s1%m,s2%m,prec)
6277	do i=1,3
6278	do j=1,3
6279	call clean_real_8(s1%s%s(i,j),s2%s%s(i,j),prec)
6280	enddo
6281	enddo
6282
6283
6284	END SUBROUTINE clean_damapspin
6285
6286	SUBROUTINE clean_spinor_8(S1,S2,prec)
6287	implicit none
6288	type (spinor_8),INTENT(INOUT)::S2
6289	type (spinor_8), intent(INOUT):: s1
6290	real(dp) prec
6291	integer i
6292
6293	do i=1,3
6294	call clean_real_8(s1%x(i),s2%x(i),prec)
6295	enddo
6296
6297	END SUBROUTINE clean_spinor_8
6298
6299	SUBROUTINE clean_res_spinor_8(S1,S2,prec)
6300	implicit none
6301	type (res_spinor_8),INTENT(INOUT)::S2
6302	type (res_spinor_8), intent(INOUT):: s1
6303	real(dp) prec
6304	integer i
6305
6306	do i=1,3
6307	call clean_double_complex(s1%x(i),s2%x(i),prec)
6308	enddo
6309
6310	END SUBROUTINE clean_res_spinor_8
6311
6312
6313	! TABLE STUFF FOR FUTURE DA
6314
6315	integer function number_mon(n,m)
6316	implicit none
6317	integer i,n,m
6318
6319	number_mon=1
6320
6321
6322	do i=n+m,max(n,m)+1,-1
6323
6324	number_mon=number_mon*i
6325	enddo
6326
6327	do i=2,min(n,m)
6328	number_mon=number_mon/i
6329	enddo
6330
6331	end function number_mon
6332
6333	integer function mul_fac(ju)
6334	implicit none
6335	integer ju(:),nv
6336	integer i,k,no
6337
6338	mul_fac=1.0_dp
6339	if(firstfac) then
6340	call make_fac
6341	firstfac=.false.
6342	endif
6343	nv=size(ju)
6344
6345	k=0
6346	do i=1,nv
6347	k=k+ju(i)
6348	enddo
6349
6350
6351	do i=1,nv
6352	if(ju(i)==0) cycle
6353
6354	mul_fac=(fac(k)/fac(k-ju(i))/fac(ju(i)))*mul_fac
6355	k=k-ju(i)
6356
6357	enddo
6358
6359
6360	end function mul_fac
6361
6362	subroutine make_fac()
6363	implicit none
6364	integer i
6365
6366	fac(0)=1.0_dp
6367	do i=1,nfac
6368	fac(i)=i*fac(i-1)
6369	enddo
6370
6371	end subroutine make_fac
6372
6373
6374	integer function pos_mon(ju,nomax,nv)
6375	implicit none
6376	integer ju(:),no,nv,nomax
6377	integer i,k,nk
6378
6379	pos_mon=0
6380	no=0
6381	do i=1,nv
6382	no=no+ju(i)
6383	enddo
6384
6385	nk=no
6386
6387	if(nk>nomax) then
6388	pos_mon=0
6389	return
6390	endif
6391
6392	do k=1,nv-1
6393	if(ju(k)/=0) then
6394	pos_mon=pos_mon+number_mon(nk,nv-k)-number_mon(nk-ju(k),nv-k)
6395	nk=nk-ju(k)
6396	endif
6397	enddo
6398	pos_mon=pos_mon+1
6399
6400	if(no>0) pos_mon=pos_mon+number_mon(no-1,nv)
6401
6402	end function pos_mon
6403
6404	integer function pos_no(no,nomax,nv)
6405	implicit none
6406	integer ju(lnv),no,nv,nomax
6407
6408	if(no==0) then
6409	pos_no=0
6410	endif
6411	if(no>nomax) then
6412	pos_no=-1
6413	endif
6414	ju=0
6415	ju(nv)=no
6416	pos_no=pos_mon(ju,nomax,nv)
6417	end function pos_no
6418
6419	subroutine find_exp(p,ju,no,nv)
6420	implicit none
6421	integer ju(:),no,nv
6422	integer i,nk,nvk,p,p0,p1,pg
6423
6424	ju=0
6425
6426	if(p==1) then
6427	return
6428	endif
6429
6430	if(p<=nv+1) then
6431	ju(nv-p+2)=1
6432	return
6433	endif
6434
6435
6436	do i=1,no
6437	p1=number_mon(i,nv)
6438	if(p1>=p) then
6439	nk=i
6440	exit
6441	endif
6442	p0=p1
6443	enddo
6444
6445
6446
6447	nvk=nv-1
6448	pg=p0
6449	do while(nvk>0)
6450
6451	p1=pg
6452	do i=0,nk
6453	p1=number_mon(nk-i,nvk-1)+p1
6454	if(p1>=p) then
6455	nk=nk-i
6456	ju(nv-nvk)=i
6457	nvk=nvk-1
6458	exit
6459	endif
6460	pg=p1
6461	enddo
6462	enddo
6463	ju(nv)=nk
6464	end subroutine find_exp
6465
6466
6467
6468	end module tree_element_MODULE

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: