-
Notifications
You must be signed in to change notification settings - Fork 39
/
Ci_tpsa.f90
19478 lines (15343 loc) · 416 KB
/
Ci_tpsa.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!The Full Polymorphic Package
!Copyright (C) Etienne Forest
! ndct=iabs(ndpt-ndptb) ! 1 if coasting, otherwise 0
! ndc2t=2*ndct ! 2 if coasting, otherwise 0
! nd2t=nd2-2*rf-ndc2t ! size of harmonic oscillators minus modulated clocks
! ndt=nd2t/2 ! ndt number of harmonic oscillators minus modulated clocks
! nd2harm=nd2t+2*rf !!!! total dimension of harmonic phase space
! ndharm=ndt+rf !!!! total number of harmonic planes
!
MODULE c_TPSA
!use newda
use definition
use file_handler
use tree_element_MODULE
use duan_zhe_map, probe_zhe=>probe,tree_element_zhe=>tree_element,dp_zhe=>dp, &
DEFAULT0_zhe=>DEFAULT0,TOTALPATH0_zhe=>TOTALPATH0,TIME0_zhe=>TIME0,ONLY_4d0_zhe=>ONLY_4d0,RADIATION0_zhe=>RADIATION0, &
NOCAVITY0_zhe=>NOCAVITY0,FRINGE0_zhe=>FRINGE0,STOCHASTIC0_zhe=>STOCHASTIC0,ENVELOPE0_zhe=>ENVELOPE0, &
DELTA0_zhe=>DELTA0,SPIN0_zhe=>SPIN0,MODULATION0_zhe=>MODULATION0,only_2d0_zhe=>only_2d0 , &
INTERNAL_STATE_zhe=>INTERNAL_STATE,ALLOC_TREE_zhe=>ALLOC_TREE
IMPLICIT NONE
public
private ety2,etdiv,ety,etyt,checksympn
integer,private::nd2par,nd2part,nd2partt
integer,private,target ::pos_of_delta
integer,private,dimension(lnv)::jfil,jfilt
private equal,DAABSEQUAL,Dequaldacon ,equaldacon ,Iequaldacon,derive,DEQUALDACONS !,AABSEQUAL 2002.10.17
private pow, GETORDER,CUTORDER,getchar,GETint,GETORDERMAP !, c_bra_v_spinmatrix
private getdiff,getdATRA ,mul,dmulsc,dscmul,GETintmat !,c_spinor_spinmatrix
private mulsc,scmul,imulsc,iscmul,DAREADTAYLORS,c_pri_c_ray,EQUALql_c_spin
private div,ddivsc,dscdiv,divsc,scdiv,idivsc,iscdiv,equalc_ray_r6r
private unaryADD,add,daddsca,dscadd,addsc,scadd,iaddsc,iscadd,print_ql
private unarySUB,subs,dsubsc,dscsub,subsc,scsub,isubsc,iscsub,c_clean_taylors
private c_allocda,c_killda,c_a_opt,K_opt,c_,c_allocdas,filter_part
private dexpt,dcost,dsint,dtant,DAPRINTTAYLORS,c_clean_yu_w,mul_ql_m,mul_ql_cm
PRIVATE GETCHARnd2,GETintnd2,dputchar,dputint, filter,check_j,c_dputint0,c_dputint0r
private GETintnd2t,equalc_cspinor_cspinor,c_AIMAG,c_real,equalc_ray_ray,EQUALql_q,EQUALq_ql,EQUALql_i,EQUALql_ql
PRIVATE DEQUAL,REQUAL,varf,varf001,equalc_spinor_cspinor,EQUALql_r !,CHARINT,pbbrav,cpbbrav
! PUBLIC VAR,ASS
private pbbra,liebra,full_absT,c_asstaylor,getcharnd2s,GETintnd2s,GETintk
private shiftda,shift000,cDEQUAL,pri,rea,cfu000,alloc_DA,alloc_c_spinmatrix,cpbbra,alloc_c_damaps
private alloc_c_damap,c_DPEKMAP,c_DPOKMAP,kill_c_damap,kill_c_spinmatrix,c_etcct,c_spinmatrix_mul_cray
private EQUALspinmatrix,c_trxtaylor,powmap,POWMAPs,alloc_c_vector_field,kill_c_vector_field,kill_c_damaps
private alloc_c_normal_form,kill_c_normal_form,c_EQUALVEC,c_spinmatrix_spinmatrix,c_IdentityEQUALVEC,qua_ql
private liebraquaternion,pow_tpsaMAP,c_concat_quaternion_ray,matrix_to_quaternion_in_c_damap,iexp_ad
private EQUALql_cmap,EQUALcmap_ql,EQUAL_complex_quaternion_c_quaternion,EQUAL_c_quaternion_complex_quaternion
private NO,ND,ND2,NP,NDPT,NV,ndptb,rf
integer, target :: NP,NO,ND,ND2,NV,NDPT,ndptb,rf
private nd_used
integer nd_used
logical(lp):: do_linear_ac_longitudinal=.true.
private old
!private map_mul_vec
logical(lp) :: old
logical(lp),target :: c_real_warning =.true.
logical(lp) :: c_mess_up_vector=.false.
real(dp) :: a_mess=0.d0 , b_mess=1.d0
integer :: i_piotr(3)= (/0,0,0/)
logical :: c_skip_gofix=.false.
PRIVATE null_it,Set_Up,de_Set_Up,LINE_L,RING_L,kill_DALEVEL,dealloc_DASCRATCH,set_up_level
private insert_da,append_da,GETINTegrate,c_pek000,c_pok000,cDEQUALDACON
private cdaddsc,cdscadd,cdsubsc,cdscsub,cdmulsc,cdscmul,cddivsc,cdscdiv
private equalc_t_ct,equalc_ct_c,matrixMAPr,MAPmatrixr,r_MAPmatrixr,c_EQUALMAP,c_IdentityEQUALMAP
private c_IdentityEQUALSPIN,c_pri_spinmatrix,c_pri_map,r_matrixMAPr
private equalc_cmap_map,c_bra_v_ct,c_bra_v_dm,equalc_cvec_vec,c_expflo,c_expflo_map
private alloc_c_factored_lie,kill_c_factored_lie,c_expflo_fac
private c_trxspinmatrix,c_inv_as,sqrtt,alloc_c_spinor,kill_c_spinor,c_complex_spinmatrix,c_trxspinmatrixda
private c_spinmatrix_add_spinmatrix,c_exp_spinmatrix,unarySUB_spinor,c_spinor_add_spinor,c_taylor_spinor
private c_IdentityEQUALSPINOR,c_spinmatrix_spinor,c_logt,c_pri_factored_lie,equalc_map_cmap
private c_expflo_fac_inv,c_logc,c_complex_spinor,c_real_spinor,GETORDERSPINMATRIX,c_pri_spinor
private c_spinor_cmap,c_adjoint_vec,c_adjoint,c_trxtaylor_da,c_spinmatrix_sub_spinmatrix,c_spinor_cmap_tpsa
PRIVATE CUTORDERMAP,CUTORDERspin,CUTORDERspinor,c_concat_tpsa,c_concat_spinor_ray,GETORDERquaternion
type(C_dalevel) c_scratchda(ndumt) !scratch levels of DA using linked list
integer, private,target :: nd2t=6,ndt=3,ndc2t=2,ndct=1,nd2harm,ndharm
!integer, private, parameter :: ndim2t=10
logical(lp), private :: c_similarity=my_false
logical(lp) :: symp =my_false
logical(lp) :: c_normal_auto=my_true,c_verbose=my_true
integer :: spin_def_tune=1 !, private
integer :: order_gofix=1
logical(lp) :: time_lie_choice=my_false,courant_snyder_teng_edwards=my_true,dosymp=my_false
private copy_damap_matrix,copy_matrix_matrix,invert_22,ALLOC_33t,kill_33t,matmul_33,print_33t
private A_OPT_C_damap,K_OPT_c_damap,equalc_t,equalt_c,daddsco,scdaddo,daddsc,scdadd,matmulr_33
private equal_real8_cmap,equal_cmap_real8,EQUAL_c_map_RAY8,EQUAL_RAY8_c_map,c_add_vf,real_mul_vec
private c_sub_vf,c_spinor_sub_spinor,matmult_33,EQUALq_i
private c_IdentityEQUALfactored,c_log_spinmatrix,c_concat_c_ray,equalc_ray_r6,equalc_r6_ray
private dotc_spinor,c_spinor_spinor,c_read_spinmatrix,c_read_map,c_concat_spinmatrix_ray
integer,private,parameter::ndd=6
private c_concat_vector_field_ray,CUTORDERVEC,kill_c_vector_field_fourier,alloc_c_vector_field_fourier
private complex_mul_vec,equal_c_vector_field_fourier,c_IdentityEQUALVECfourier
private c_add_map,c_sub_map,c_read_spinor,flatten_c_factored_lie_r,c_EQUALcray,c_read_quaternion
integer :: n_fourier=12,n_extra=0
logical :: remove_tune_shift=.false.
complex(dp) :: n_cai=-2*i_
integer :: complex_extra_order=0
logical :: special_extra_order_1=.true.
real(dp) :: epso_factor =1000.d0 ! for log
logical(lp):: extra_terms_log=.false.
logical :: add_constant_part_concat=.true.,assume_c_quaternion_normalised=.true.
private EQUAL_c_spinmatrix_probe,EQUAL_c_spinmatrix_3_by_3,EQUAL_3_by_3_probe,EQUAL_probe_c_spinmatrix
private EQUAL_probe_3_by_3,equalc_cspinor_spinor,EQUAL_3_by_3_c_spinmatrix
private EQUALq_r,EQUALq_8_c,EQUALq_c_8,EQUALq,POWq,c_invq,subq,mulq,addq,alloc_c_quaternion,kill_c_quaternion
private c_pri_quaternion,CUTORDERquaternion,c_trxquaternion,EQUALq_c_r,EQUALq_r_c,mulcq,c_exp_quaternion
private equalc_quaternion_c_spinor,equalc_spinor_c_quaternion,unarySUB_q,c_trxquaternion_tpsa
private c_exp_vectorfield_on_quaternion,c_vector_field_quaternion,addql,subql,mulqdiv,powql,quaternion_to_matrix
private copy_tree_into_tree_zhe,absq2,absq
!private equal_map_real8,equal_map_complex8,equal_real8_map,equal_complex8_map
real(dp) dts
real(dp), private :: sj(6,6)
logical :: use_new_stochastic_normal_form=.true.
logical :: qphase=.false.,qphasedef=.false.
integer :: size_tree=15
integer :: ind_spin0(3,3),ind_spin(3,3),k1_spin(9),k2_spin(9)
logical :: hypercube_integration = .true.
!private ip_mat,jp_mat,jt_mat
real(dp) ip_mat(3,6,6),jp_mat(3,6,6),jt_mat(6,6)
character(24) formatlf(6)
!-----------------------------------
logical :: inside_normal=.false.,bmad_automatic=.false.,spin_automatic=.false.
integer i_alloc
logical :: sagan_gen =.false.
type c_linear_map
complex(dp) mat(6,6)
complex(dp) q(0:3,0:6)
end type c_linear_map
type c_lattice_function
real(dp) :: E(3,6,6) =0 ,K(3,6,6) =0,B(3,6,6) =0
real(dp) :: H(3,6,6) = 0
real(dp) :: S(1:3,1:3,0:6) =0
real(dp) ::phase(3) =0 ,damping(3) =0 , spin(2) =0,fix(6) =0
type(fibre), pointer :: f => null()
type(integration_node), pointer :: t => null()
logical :: symplectic = .true.
!!!! radiation quantity
real(dp) :: sigmas(6,6)
end type c_lattice_function
type(c_linear_map) q_phasor,qi_phasor
INTERFACE abs_square
MODULE PROCEDURE absq2
END INTERFACE
INTERFACE abs
MODULE PROCEDURE absq
END INTERFACE
INTERFACE assignment (=)
MODULE PROCEDURE EQUAL
MODULE PROCEDURE EQUALq_r
MODULE PROCEDURE EQUALq_8_c
MODULE PROCEDURE EQUALq_c_8
MODULE PROCEDURE EQUALql_q
MODULE PROCEDURE EQUALql_c_spin
MODULE PROCEDURE EQUALq_ql
MODULE PROCEDURE EQUALql_i
MODULE PROCEDURE EQUALq_i
MODULE PROCEDURE EQUALql_r
MODULE PROCEDURE EQUALq
MODULE PROCEDURE EQUALq_c_r
MODULE PROCEDURE EQUALq_r_c
MODULE PROCEDURE EQUALql_ql
MODULE PROCEDURE EQUAL_c_l_f ! c_lattice_function = true or false (symplectic)
MODULE PROCEDURE EQUAL_complex_quaternion_c_quaternion
MODULE PROCEDURE EQUAL_c_quaternion_complex_quaternion
MODULE PROCEDURE EQUALql_cmap
MODULE PROCEDURE EQUALcmap_ql
MODULE PROCEDURE cDEQUAL
MODULE PROCEDURE DEQUAL ! added 2002.10.17 ! check2002.10.17
MODULE PROCEDURE REQUAL ! added 2002.10.17 ! check2002.10.17
MODULE PROCEDURE cDEQUALDACON
MODULE PROCEDURE Dequaldacon
MODULE PROCEDURE DEQUALDACONS ! C_taylor(:) = real(dp))
MODULE PROCEDURE Iequaldacon
MODULE PROCEDURE equalc_t_ct !# complex Taylor
MODULE PROCEDURE equalc_ct_c
MODULE PROCEDURE equalc_t !# Taylor
MODULE PROCEDURE equalt_c
MODULE PROCEDURE equal_real8_cmap !# put cmap in real_8(6)
MODULE PROCEDURE equal_cmap_real8 ! put real_8(6) in cmap
MODULE PROCEDURE c_EQUALMAP !# c_damap=c_damap
MODULE PROCEDURE c_IdentityEQUALMAP
MODULE PROCEDURE c_IdentityEQUALSPIN
MODULE PROCEDURE c_IdentityEQUALSPINOR !# c_spinor = 0,1,2, or 3
MODULE PROCEDURE c_IdentityEQUALVEC
MODULE PROCEDURE c_IdentityEQUALfactored
MODULE PROCEDURE c_IdentityEQUALVECfourier
MODULE PROCEDURE EQUAL_c_spinmatrix_probe
MODULE PROCEDURE EQUAL_c_spinmatrix_3_by_3
module procedure EQUAL_3_by_3_c_spinmatrix
MODULE PROCEDURE EQUAL_3_by_3_probe
MODULE PROCEDURE EQUAL_probe_3_by_3
MODULE PROCEDURE EQUAL_probe_c_spinmatrix
MODULE PROCEDURE matrixMAPr
MODULE PROCEDURE r_matrixMAPr
MODULE PROCEDURE MAPmatrixr
MODULE PROCEDURE r_MAPmatrixr
! MODULE PROCEDURE c_DPEKMAP
! MODULE PROCEDURE c_DPOKMAP
MODULE PROCEDURE EQUALspinmatrix
MODULE PROCEDURE EQUAL_c_map_RAY8 !# c_damap=probe_8
MODULE PROCEDURE EQUAL_RAY8_c_map !# probe_8=c_damap
MODULE PROCEDURE c_EQUALcray !# c_ray=integer
MODULE PROCEDURE equalc_cmap_map
MODULE PROCEDURE equalc_map_cmap
MODULE PROCEDURE equalc_cvec_vec
MODULE PROCEDURE c_EQUALVEC !# c_vector_field=c_vector_field
MODULE PROCEDURE equalc_cspinor_cspinor
MODULE PROCEDURE equalc_ray_r6
MODULE PROCEDURE equalc_ray_r6r
MODULE PROCEDURE equalc_r6_ray
MODULE PROCEDURE equalc_ray_ray
MODULE PROCEDURE equal_c_vector_field_fourier
MODULE PROCEDURE equalc_spinor_cspinor
MODULE PROCEDURE equalc_cspinor_spinor
MODULE PROCEDURE flatten_c_factored_lie_r !# same as flatten_c_factored_lie
MODULE PROCEDURE equalc_quaternion_c_spinor
MODULE PROCEDURE equalc_spinor_c_quaternion
MODULE PROCEDURE equal_map_real8
MODULE PROCEDURE equal_map_complex8 ! replaces c_dpokmap
MODULE PROCEDURE equal_real8_map
MODULE PROCEDURE equal_complex8_map ! replaces c_dpekmap
MODULE PROCEDURE copy_tree_into_tree_zhe
end INTERFACE
INTERFACE OPERATOR (+)
MODULE PROCEDURE unaryADD !@2 This is a unary operation
MODULE PROCEDURE add
MODULE PROCEDURE addq
MODULE PROCEDURE addql
MODULE PROCEDURE daddsco !# c_damap + real(6)
MODULE PROCEDURE scdaddo !# real(6) + c_damap
MODULE PROCEDURE daddsc !# c_damap + probe_8
MODULE PROCEDURE scdadd !# probe_8 + c_damap
MODULE PROCEDURE cdaddsc
MODULE PROCEDURE daddsca
MODULE PROCEDURE cdscadd
MODULE PROCEDURE iscadd
MODULE PROCEDURE dscadd
MODULE PROCEDURE addsc
MODULE PROCEDURE scadd
MODULE PROCEDURE iaddsc
MODULE PROCEDURE c_spinmatrix_add_spinmatrix
MODULE PROCEDURE c_spinor_add_spinor ! adding c_spinor
MODULE PROCEDURE c_add_vf !# adding vector field for exp(ad VF)
MODULE PROCEDURE c_add_map
END INTERFACE
INTERFACE OPERATOR (-)
MODULE PROCEDURE unarySUB
MODULE PROCEDURE subq
MODULE PROCEDURE subql
MODULE PROCEDURE subs
MODULE PROCEDURE dsubsc
MODULE PROCEDURE dscsub
MODULE PROCEDURE cdsubsc
MODULE PROCEDURE cdscsub
MODULE PROCEDURE subsc
MODULE PROCEDURE scsub
MODULE PROCEDURE isubsc
MODULE PROCEDURE iscsub
MODULE PROCEDURE unarySUB_vec
MODULE PROCEDURE unarySUB_q
MODULE PROCEDURE unarySUB_spinor
MODULE PROCEDURE c_spinor_sub_spinor
MODULE PROCEDURE c_spinmatrix_sub_spinmatrix
MODULE PROCEDURE c_sub_vf
MODULE PROCEDURE c_sub_map
END INTERFACE
INTERFACE OPERATOR (*)
! MODULE PROCEDURE transform_vector_field_by_map
MODULE PROCEDURE mul
MODULE PROCEDURE mulq
MODULE PROCEDURE mulql
MODULE PROCEDURE mulcq
MODULE PROCEDURE mul_ql_m
MODULE PROCEDURE mul_ql_cm
MODULE PROCEDURE dmulsc
MODULE PROCEDURE cdmulsc
MODULE PROCEDURE dscmul
MODULE PROCEDURE cdscmul
MODULE PROCEDURE mulsc
MODULE PROCEDURE scmul
MODULE PROCEDURE imulsc
MODULE PROCEDURE iscmul
MODULE PROCEDURE c_concat !# c_damap o c_damap
MODULE PROCEDURE c_trxtaylor ! c_taylor o c_damap
MODULE PROCEDURE c_bra_v_ct !# F.grad taylor !v1
MODULE PROCEDURE c_bra_v_dm !# c_damap=(F.grad) c_damap !v2
! MODULE PROCEDURE c_bra_v_v ! (exp(F.grad) H) . grad
MODULE PROCEDURE c_spinmatrix_spinmatrix !# Spinmatrix*Spinmatrix
MODULE PROCEDURE c_complex_spinmatrix !# c*Spinmatrix
MODULE PROCEDURE c_taylor_spinor !# taylor * spinor
MODULE PROCEDURE c_complex_spinor !# complex * spinor
MODULE PROCEDURE c_real_spinor !# real(dp) * spinor
MODULE PROCEDURE c_spinor_spinor !# spinor x spinor
MODULE PROCEDURE c_trxspinmatrix !# c_spinmatrix= c_spinmatrix * c_damap
MODULE PROCEDURE c_trxquaternion !# c_quaternion= c_quaternion * c_damap
MODULE PROCEDURE c_spinmatrix_spinor !# matrix * spinor
MODULE PROCEDURE c_spinor_cmap !# spinor * c_damap
MODULE PROCEDURE real_mul_vec ! real(dp)*vf
MODULE PROCEDURE complex_mul_vec !# complex(dp)*vf
MODULE PROCEDURE real_mul_map !# real(dp)*c_damap
! MODULE PROCEDURE c_spinor_spinmatrix !# spinor.L spinmatrix
MODULE PROCEDURE c_vector_field_quaternion !# (f.grad + q)quaternion
! bra_v_q(f,quaternion) = f.grad quaternion
! uses map_mul_vec
MODULE PROCEDURE map_mul_vec_q ! c_damap * c_vector_field means "transform field by map"
END INTERFACE
INTERFACE OPERATOR (.o.)
!
module procedure c_concat_c_ray !# c_taylor .o. c_ray
module procedure c_concat_map_ray !# c_ray= c_damap .o. c_ray
module procedure c_trxtaylor_da !# c_taylor= c_taylor .o. c_damap
module procedure c_concat_tpsa !# c_damap .o. c_damap
module procedure c_concat_spinor_ray !# c_spinor= c_spinor .o. c_ray
! not overloaded
! module procedure c_concat_spinmatrix_ray !# c_spinmatrix= c_spinmatrix .o. c_ray
! module procedure c_concat_quaternion_ray !# c_quaternion= c_quaternion .o. c_ray
! Instead these overloaded
MODULE PROCEDURE c_spinmatrix_mul_cray !# c_ray%s = spin_matrix.o.c_ray%x c_ray%S(1:3)
MODULE PROCEDURE c_quaternion_mul_cray !# c_ray%q = c_quaternion.o.c_ray%x c_ray%q c_quaternion.o.c_ray%x**(-1)
! module procedure c_concat_vector_field_ray !# c_vector_field .o. cray
module procedure c_trxspinmatrixda !# c_spinmatrix=c_spinmatrix .o. c_damap
module procedure c_trxquaternion_tpsa !# c_quaternion= c_quaternion .o. c_damap
MODULE PROCEDURE c_spinor_cmap_tpsa !# spinor * c_damap
END INTERFACE
INTERFACE OPERATOR (.oo.)
module procedure pow_tpsaMAP
END INTERFACE
INTERFACE OPERATOR (/)
MODULE PROCEDURE div
MODULE PROCEDURE cddivsc
MODULE PROCEDURE ddivsc
MODULE PROCEDURE cdscdiv
MODULE PROCEDURE dscdiv
MODULE PROCEDURE divsc
MODULE PROCEDURE scdiv
MODULE PROCEDURE idivsc
MODULE PROCEDURE iscdiv
MODULE PROCEDURE mulqdiv
END INTERFACE
INTERFACE OPERATOR (**)
MODULE PROCEDURE POW
MODULE PROCEDURE POWQ
MODULE PROCEDURE POWql
MODULE PROCEDURE powmap
MODULE PROCEDURE powmaps
END INTERFACE
! New Operators
INTERFACE OPERATOR (.cmono.)
MODULE PROCEDURE c_dputint0 !@1 single integer
MODULE PROCEDURE c_dputint0r !@1 single integer
MODULE PROCEDURE dputintr !@1 Accepts J(nv)
MODULE PROCEDURE dputcharr !@1 Accepts String such as '12'
MODULE PROCEDURE dputint !@1 Accepts J(nv)
MODULE PROCEDURE dputchar !@1 Accepts String such as '12'
END INTERFACE
INTERFACE OPERATOR (.var.)
MODULE PROCEDURE varf !@1 replaces var (overloads DAVAR)
MODULE PROCEDURE varf001 !@1 replaces var001
END INTERFACE
INTERFACE OPERATOR (.d.)
MODULE PROCEDURE getdiff !@1 takes derivatives
END INTERFACE
INTERFACE OPERATOR (.i.)
MODULE PROCEDURE GETINTegrate !@1 takes anti-derivatives
END INTERFACE
INTERFACE OPERATOR (.SUB.)
MODULE PROCEDURE GETORDER
MODULE PROCEDURE getchar
MODULE PROCEDURE GETint
MODULE PROCEDURE GETORDERMAP ! with negative integer map.sub.i the spin is handled with iabs(i)-1
MODULE PROCEDURE GETORDERSPINMATRIX ! FOR SPIN MATRICES
MODULE PROCEDURE GETORDERquaternion
END INTERFACE
INTERFACE OPERATOR (.index.)
MODULE PROCEDURE GETintmat
END INTERFACE
INTERFACE OPERATOR (.PAR.)
MODULE PROCEDURE getcharnd2
MODULE PROCEDURE GETintnd2
END INTERFACE
INTERFACE OPERATOR (.part.)
MODULE PROCEDURE GETintnd2t
END INTERFACE
INTERFACE OPERATOR (<=)
MODULE PROCEDURE getcharnd2s
MODULE PROCEDURE GETintnd2s
MODULE PROCEDURE GETintk
END INTERFACE
INTERFACE OPERATOR (.CUT.)
MODULE PROCEDURE CUTORDER
MODULE PROCEDURE CUTORDERMAP
MODULE PROCEDURE CUTORDERvec
MODULE PROCEDURE CUTORDERspin
MODULE PROCEDURE CUTORDERspinor
MODULE PROCEDURE CUTORDERquaternion
END INTERFACE
INTERFACE OPERATOR (.dot.)
MODULE PROCEDURE dotc_spinor ! Used internally primarily
END INTERFACE
INTERFACE OPERATOR (.K.)
MODULE PROCEDURE getdATRA ! Used internally primarily
END INTERFACE
INTERFACE OPERATOR (.pb.)
MODULE PROCEDURE pbbra
END INTERFACE
INTERFACE OPERATOR (.cpb.)
MODULE PROCEDURE cpbbra
END INTERFACE
INTERFACE OPERATOR (.lb.)
!
MODULE PROCEDURE liebra !# F=<G,H> includes spin now
MODULE PROCEDURE liebraquaternion !# used by liebra
END INTERFACE
! INTERFACE getvf
! MODULE PROCEDURE pbbrav
!end INTERFACE getvf
!INTERFACE cgetvf
! MODULE PROCEDURE cpbbrav
!end INTERFACE cgetvf
! intrisic functions overloaded
INTERFACE checksymp
MODULE PROCEDURE checksympn
end INTERFACE checksymp
INTERFACE q_part
MODULE PROCEDURE qua_ql
end INTERFACE q_part
INTERFACE c_phasor
MODULE PROCEDURE from_phasor
end INTERFACE c_phasor
INTERFACE ci_phasor
MODULE PROCEDURE to_phasor
end INTERFACE ci_phasor
! Exponential of Lie Operators
INTERFACE clean
! MODULE PROCEDURE c_clean
MODULE PROCEDURE c_clean_spinor
MODULE PROCEDURE c_clean_taylor
MODULE PROCEDURE c_clean_spinmatrix
MODULE PROCEDURE c_clean_damap
MODULE PROCEDURE c_clean_vector_field
MODULE PROCEDURE c_clean_yu_w
MODULE PROCEDURE c_clean_quaternion
MODULE PROCEDURE c_clean_taylors
MODULE PROCEDURE c_clean_vector_field_fourier
end INTERFACE clean
! Exponential of Lie Operators
INTERFACE c_simil
MODULE PROCEDURE c_adjoint
MODULE PROCEDURE c_adjoint_vec
end INTERFACE c_simil
INTERFACE texp_inv
MODULE PROCEDURE c_expflo_fac_inv
END INTERFACE
INTERFACE exp_inv
MODULE PROCEDURE c_expflo_fac_inv ! v9
END INTERFACE
INTERFACE real
MODULE PROCEDURE c_real
END INTERFACE
INTERFACE aimag
MODULE PROCEDURE c_AIMAG
END INTERFACE
!exp_mat(f,m) is a subroutine to exponentiate matrices
INTERFACE exp
MODULE PROCEDURE c_expflo ! flow on c_taylor !v3
MODULE PROCEDURE c_expflo_map !v4
MODULE PROCEDURE c_expflo_fac !v7
MODULE PROCEDURE c_exp_spinmatrix
MODULE PROCEDURE c_exp_vectorfield_on_quaternion !v6
MODULE PROCEDURE c_exp_quaternion
MODULE PROCEDURE exp_ad ! exp(<F,>)F F is a vector field !v5
END INTERFACE
INTERFACE iexp
MODULE PROCEDURE iexp_ad
END INTERFACE iexp
INTERFACE texp
MODULE PROCEDURE c_expflo
MODULE PROCEDURE c_expflo_map
MODULE PROCEDURE c_expflo_fac
MODULE PROCEDURE c_exp_spinmatrix
MODULE PROCEDURE c_exp_vectorfield_on_quaternion
MODULE PROCEDURE c_exp_quaternion
MODULE PROCEDURE exp_ad ! exp(<F,>)F F is a vector field
END INTERFACE
INTERFACE abs
MODULE PROCEDURE DAABSEQUAL ! remove 2002.10.17
END INTERFACE
INTERFACE dabs
MODULE PROCEDURE DAABSEQUAL ! remove 2002.10.17
END INTERFACE
INTERFACE exp
MODULE PROCEDURE dexpt
END INTERFACE
INTERFACE dexp
MODULE PROCEDURE dexpt
END INTERFACE
INTERFACE cexp
MODULE PROCEDURE dexpt
END INTERFACE
INTERFACE cdexp
MODULE PROCEDURE dexpt
END INTERFACE
INTERFACE log
MODULE PROCEDURE c_logt
MODULE PROCEDURE c_logc
MODULE PROCEDURE c_logf !# log of a map see subroutine c_flofacg
MODULE PROCEDURE c_log_spinmatrix !# spinor=log(s)
! c_logf_spin is not overloaded
END INTERFACE
INTERFACE cos
MODULE PROCEDURE dcost
END INTERFACE
INTERFACE cdcos
MODULE PROCEDURE dcost
END INTERFACE
INTERFACE dcos
MODULE PROCEDURE dcost
END INTERFACE
INTERFACE ccos
MODULE PROCEDURE dcost
END INTERFACE
INTERFACE sin
MODULE PROCEDURE dsint
END INTERFACE
INTERFACE cdsin
MODULE PROCEDURE dsint
END INTERFACE
INTERFACE ccsin
MODULE PROCEDURE dsint
END INTERFACE
INTERFACE dsin
MODULE PROCEDURE dsint
END INTERFACE
INTERFACE sqrt
MODULE PROCEDURE sqrtt
END INTERFACE
INTERFACE tan
MODULE PROCEDURE dtant
END INTERFACE
INTERFACE dtan
MODULE PROCEDURE dtant
END INTERFACE
! Non-intrisic Functions
INTERFACE c_pek
MODULE PROCEDURE c_pek000 !
END INTERFACE
INTERFACE c_pok
MODULE PROCEDURE c_pok000 !
END INTERFACE
INTERFACE shiftda
MODULE PROCEDURE shift000 ! not private
END INTERFACE
! INTERFACE var
! MODULE PROCEDURE var000 ! not private
! MODULE PROCEDURE var001 ! not private
! END INTERFACE
INTERFACE cfu
MODULE PROCEDURE c_cfu000 ! not private
END INTERFACE
INTERFACE full_abs
MODULE PROCEDURE full_absT
END INTERFACE
INTERFACE daread
MODULE PROCEDURE c_rea
module procedure c_read_spinmatrix
module procedure c_read_map
MODULE PROCEDURE c_read_spinor
MODULE PROCEDURE DAREADTAYLORS
MODULE PROCEDURE c_read_quaternion
END INTERFACE
INTERFACE read
MODULE PROCEDURE c_rea
module procedure c_read_spinmatrix
module procedure c_read_map
MODULE PROCEDURE c_read_spinor
MODULE PROCEDURE DAREADTAYLORS
MODULE PROCEDURE c_read_quaternion
END INTERFACE
INTERFACE daprint
MODULE PROCEDURE c_pri
MODULE PROCEDURE c_pri_spinmatrix
MODULE PROCEDURE c_pri_map
MODULE PROCEDURE c_pri_vec
MODULE PROCEDURE c_pri_factored_lie
MODULE PROCEDURE c_pri_spinor
MODULE PROCEDURE print_33t
MODULE PROCEDURE c_pri_quaternion
MODULE PROCEDURE DAPRINTTAYLORS
MODULE PROCEDURE c_pri_c_ray
END INTERFACE
INTERFACE print
MODULE PROCEDURE c_pri
MODULE PROCEDURE c_pri_spinmatrix
MODULE PROCEDURE c_pri_map
MODULE PROCEDURE c_pri_vec
MODULE PROCEDURE c_pri_factored_lie
MODULE PROCEDURE c_pri_spinor
MODULE PROCEDURE print_33t
MODULE PROCEDURE DAPRINTTAYLORS
MODULE PROCEDURE c_pri_quaternion
MODULE PROCEDURE print_ql
MODULE PROCEDURE c_pri_c_ray
END INTERFACE
! Constructors and Destructors
INTERFACE alloc
MODULE PROCEDURE c_allocda
MODULE PROCEDURE c_a_opt
MODULE PROCEDURE A_OPT_c_damap
MODULE PROCEDURE c_allocdas
MODULE PROCEDURE alloc_c_spinmatrix
MODULE PROCEDURE alloc_c_vector_field_fourier
MODULE PROCEDURE alloc_c_damap
MODULE PROCEDURE alloc_c_damaps
MODULE PROCEDURE alloc_c_vector_field
MODULE PROCEDURE alloc_c_factored_lie
MODULE PROCEDURE alloc_c_normal_form
MODULE PROCEDURE alloc_c_spinor
MODULE PROCEDURE alloc_c_yu_w
MODULE PROCEDURE alloc_c_quaternion
END INTERFACE
INTERFACE ALLOC_nn
MODULE PROCEDURE ALLOC_33t
END INTERFACE
INTERFACE kill_nn
MODULE PROCEDURE kill_33t
END INTERFACE
INTERFACE matmul_nn
MODULE PROCEDURE matmul_33
END INTERFACE
INTERFACE matmulr_nn
MODULE PROCEDURE matmulr_33
MODULE PROCEDURE matmult_33
END INTERFACE
INTERFACE KILL
MODULE PROCEDURE c_killda
MODULE PROCEDURE c_killdas
MODULE PROCEDURE K_opt
MODULE PROCEDURE K_OPT_c_damap
MODULE PROCEDURE kill_c_damap
MODULE PROCEDURE kill_c_damaps
MODULE PROCEDURE kill_c_spinmatrix
MODULE PROCEDURE kill_c_vector_field
MODULE PROCEDURE kill_c_factored_lie
MODULE PROCEDURE kill_c_vector_field_fourier
MODULE PROCEDURE kill_c_normal_form
MODULE PROCEDURE kill_c_spinor
MODULE PROCEDURE kill_c_yu_w
MODULE PROCEDURE kill_c_quaternion
END INTERFACE
INTERFACE alloctpsa
MODULE PROCEDURE c_allocda
END INTERFACE
INTERFACE KILLtpsa
MODULE PROCEDURE c_killda
END INTERFACE
INTERFACE MAKEquaternion
MODULE PROCEDURE matrix_to_quaternion_in_c_damap
END INTERFACE
INTERFACE MAKESO3
MODULE PROCEDURE quaternion_to_matrix_in_c_damap
MODULE PROCEDURE c_linear_map_to_matrix
MODULE PROCEDURE c_linear_map_to_3_by_3_by_6
MODULE PROCEDURE quaternion_to_matrix
END INTERFACE
! management routines
INTERFACE ass
MODULE PROCEDURE c_asstaylor !2000.12.25
END INTERFACE
INTERFACE AVERAGE
MODULE PROCEDURE c_AVERAGE !2000.12.25
END INTERFACE
CONTAINS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine c_get_indices(n,mf)
!#general: informationc_%
!# In the arrary n(>11), important parameters of the normal
!# form can be retrieved.
!# If mf/=0, they are printed on file mf.
!# n(1)=NO
!# n(2)=ND
!# n(3)=ND2
!# n(4)=NV
!# n(5)=Ndpt
!# n(6)=ndptb
!# n(7)=np
!# n(8)=rf*2
!# n(9)=ndc2t
!# n(10)=nd2t
!# n(11)=nd2harm
implicit none
integer n(:),mf
if(size(n)<11) then
write(6,*) " index array to small in c_get_indices "
stop
endif
n(1)=NO
n(2)=ND
n(3)=ND2
n(4)=NV
n(5)=Ndpt
n(6)=ndptb
n(7)=np
n(8)=rf*2
n(9)=ndc2t
n(10)=nd2t
n(11)=nd2harm
if(mf/=0) then
write(mf,"(11(a7))") " NO "," ND "," ND2"," NV "," NDPT "," NDPTB "," NP "," RF "," NDC2T ", &
" ND2T "," HARM "
write(mf,"(11(3x,i2,2x))") n(1:11)
endif
end subroutine c_get_indices
subroutine locally_set_da_pointers()
use da_arrays,only : c_
implicit none
c_%no => no
c_%nv => nv
c_%nd => nd
c_%nd2 => nd2
c_%np => np
c_%ndpt => ndpt
end subroutine locally_set_da_pointers
subroutine c_count_taylor(n,ns,ne)
!#restricted : information
!# Counts number of c_taylor allocated;
!# it is used for debugging purposes.
implicit none
integer n,ns,ne,i
call c_count_da(n)
ns=0
do i=1,ndumt
ns=c_scratchda(i)%n+ns
enddo
ne=n-ns
end subroutine c_count_taylor
FUNCTION scdadd( S2,S1 )
implicit none
TYPE (probe_8) scdadd
TYPE (c_damap), INTENT (IN) :: S1
type(probe) , INTENT (IN) :: S2
integer localmaster,i,j,nd2h,dc
type(taylor) d
call alloc(d)
! call ass(scdadd)
scdadd%u=my_false
scdadd%E_ij=0.0_dp
scdadd%nac=s2%nac
scdadd%use_q=s2%use_q
if(doing_ac_modulation_in_ptc) then
dc=2*rf
else
dc=0
endif
if(c_%ndpt==0) then ! 1
nd2h=nd2t
do i=1,nd2h ! from 1-4 or 1-6 (if ndpt=0)
localmaster=master
call ass(scdadd%x(i))
! scdadd%x(i)=s1%m%v(i)+s2%x(i)
d=s1%V(i)+s2%x(i)-(s1%V(i).sub.'0')
scdadd%x(i)=d
master=localmaster
enddo
do i=nd2h+1,6
localmaster=master
call ass(scdadd%x(i))
! ! if((c_%npara==5+dc).AND.I==5) then ! npr
if(((c_%npara==5+dc).AND.I==5+ndpt_bmad).or.((c_%npara==3+dc).AND.I==5+ndpt_bmad)) then ! npr
scdadd%x(i)=s2%x(i)+(1.0_dp.mono.c_%npara)
else
scdadd%x(i)=s2%x(i)
endif
master=localmaster
enddo
else ! 1
nd2h=nd2t+2
do i=1,nd2h ! from 1-4 or 1-6 (if ndpt=0)
localmaster=master
call ass(scdadd%x(i))
! scdadd%x(i)=s1%m%v(i)+s2%x(i)
d=s1%V(i)+s2%x(i)-(s1%V(i).sub.'0')
scdadd%x(i)=d
master=localmaster
enddo
endif ! 1
do i=1,scdadd%nac
localmaster=master
call ass(scdadd%AC(i)%x(1))
! scdadd%x(i)=s1%m%v(i)+s2%x(i)
j=2*scdadd%nac-(2*i-1)
d=s1%V(C_%ND2-j) +addclock*s2%AC(i)%x(1)
scdadd%ac(i)%x(1)=d
master=localmaster
localmaster=master
j=2*scdadd%nac-(2*i)
call ass(scdadd%AC(i)%x(2))
! scdadd%x(i)=s1%m%v(i)+s2%x(i)
d=s1%V(C_%ND2-j) +addclock*s2%AC(i)%x(2)
scdadd%ac(i)%x(2)=d
master=localmaster
localmaster=master
call ass(scdadd%AC(i)%om)
! call ass(scdadd%AC%t)
! scdadd%x(i)=s1%m%v(i)+s2%x(i)
scdadd%AC(i)%om=s2%AC(i)%om
scdadd%AC(i)%t=s2%AC(i)%t
master=localmaster
enddo
! endif
! if(use_quaternion) THEN
DO J=0,3
localmaster=master
call ass(scdadd%q%x(J))
d=S1%q%x(j)
scdadd%q%x(J)=d
master=localmaster
ENDDO
!else
DO I=1,3
! call ass(scdadd%s%x(i))
DO J=1,3
localmaster=master
call ass(scdadd%s(J)%x(i))
d=S1%S%s(I,J)
scdadd%s(J)%x(i)=d
master=localmaster
ENDDO
ENDDO
!endif
scdadd%damps=s1%damps
scdadd%d_spin=s1%d_spin
scdadd%b_kin=s1%b_kin
scdadd%e_ij=s1%e_ij
scdadd%x0(1:6)=s2%x
call kill(d)
END FUNCTION scdadd
FUNCTION daddsc( S1,S2 )
implicit none
TYPE (probe_8) daddsc
TYPE (c_damap), INTENT (IN) :: S1
type(probe) , INTENT (IN) :: S2
integer localmaster,i,j,nd2h,dc
type(taylor) d