Merge pull request #522 from bmad-sim/ptc/10

Ptc/10

forest/code/Se_status.f90

@@ -144,7 +144,7 @@ module S_status
  private track_TREE_G_complexr,track_TREE_G_complexp,track_TREE_probe_complexr,track_TREE_probe_complexp_new
  real(dp),TARGET ::INITIAL_CHARGE=1
  logical :: mcmillan=.false.
- real(dp) :: radfac=1 ! to fudge radiation (lower it)
+ real(dp) :: radfac=1,flucfac=1 ! to fudge radiation (lower it)
  TYPE B_CYL
  integer firsttime
  integer, POINTER :: nmul,n_mono !,nmul_e,n_mono_e
@@ -264,7 +264,7 @@ real(dp) function cflucf(p)
  ! if(junk_e) then 
  ! cflucf=cfluc*p%p0c**5
  ! else
- cflucf=cfluc0*twopii/p%GAMMA0I**5/p%MASS**2
+ cflucf=flucfac*cfluc0*twopii/p%GAMMA0I**5/p%MASS**2
  ! endif
  end function cflucf
 
@@ -1529,6 +1529,8 @@ subroutine S_init(STATE,NO1,NP1,pack,ND2,NPARA,number_of_clocks)
  INTEGER,optional :: ND2,NPARA,number_of_clocks
  INTEGER ND2l,NPARAl,n_acc,no1c,nv,i
  LOGICAL(lp) package
+
+
  do_damping=.false.
  do_spin=.false.
  if(state%radiation) do_damping=.true.
@@ -1589,7 +1591,7 @@ subroutine S_init(STATE,NO1,NP1,pack,ND2,NPARA,number_of_clocks)
  ! endif
  ! write(6,*) NO1,ND1,NP1,NDEL,NDPT1
  !pause 678
-
+ 
  CALL INIT(NO1,ND1,NP1+NDEL+2*n_acc,NDPT1,PACKAGE)
  nv=2*nd1+NP1+NDEL+2*n_acc
 
@@ -1611,12 +1613,14 @@ subroutine S_init(STATE,NO1,NP1,pack,ND2,NPARA,number_of_clocks)
 ND1=ND1+n_acc
  if(use_complex_in_ptc) call c_init(NO1c,nd1,np1+ndel,ndpt1,n_acc,ptc=my_false) ! PTC false because we will not use the real FPP for acc modulation
  n_rf=n_acc
-
+  previous_newtpsa=newtpsa
  END subroutine S_init
 
  subroutine kill_map_cp()
  implicit none
-
+ logical present_newtpsa
+ present_newtpsa=newtpsa
+ newtpsa=previous_newtpsa 
  if(associated(dz_8)) then
  call kill(dz_8)
  deallocate(dz_8)
@@ -1635,7 +1639,7 @@ subroutine kill_map_cp()
  deallocate(dz_c)
  nullify(dz_c)
  endif 
-
+ newtpsa=present_newtpsa
  end subroutine kill_map_cp
 
 

forest/code/Si_def_element.f90

@@ -2031,68 +2031,68 @@ SUBROUTINE transfer_ANBN(EL,ELP,VR,DVR,VP,DVP,T)
 
  if(EL%KIND==kind1) return
 
- if(associated(EL%ramp)) then
-
- if(EL%KIND/=kind15) then
- do n=1,EL%P%NMUL
- EL%BN(N)= EL%ramp%table(0)%bn(n)
- EL%AN(N)= EL%ramp%table(0)%an(n)
- ELP%BN(N)= ELP%ramp%table(0)%bn(n)
- ELP%AN(N)= ELP%ramp%table(0)%an(n)
- enddo 
- else
- EL%VOLT=EL%ramp%table(0)%bn(1)*COS(twopi*EL%ramp%table(0)%an(1)*T/clight+EL%ramp%table(0)%bn(2))+EL%ramp%table(0)%an(2)
- ELP%VOLT=EL%ramp%table(0)%bn(1)*COS(twopi*EL%ramp%table(0)%an(1)*T/clight+EL%ramp%table(0)%bn(2))+EL%ramp%table(0)%an(2)
- write(6,*) " volt ",el%volt,EL%ramp%table(0)%bn(1)
- endif
-
- if(EL%ramp%table(0)%b_t/=0.0_dp) then
- if(EL%parent_fibre%PATCH%TIME==0) EL%parent_fibre%PATCH%TIME=2
- if(EL%parent_fibre%PATCH%TIME==1) EL%parent_fibre%PATCH%TIME=3
- EL%parent_fibre%PATCH%b_T=EL%ramp%table(0)%b_t
- else
- if(EL%parent_fibre%PATCH%TIME==2) EL%parent_fibre%PATCH%TIME=0
- if(EL%parent_fibre%PATCH%TIME==3) EL%parent_fibre%PATCH%TIME=1
- EL%parent_fibre%PATCH%b_T=0.0_dp
- endif
+! if(associated(EL%ramp)) then
+! 
+! if(EL%KIND/=kind15) then
+! do n=1,EL%P%NMUL
+! EL%BN(N)= EL%ramp%table(0)%bn(n)
+! EL%AN(N)= EL%ramp%table(0)%an(n)
+! ELP%BN(N)= ELP%ramp%table(0)%bn(n)
+! ELP%AN(N)= ELP%ramp%table(0)%an(n)
+! enddo 
+! else
+! EL%VOLT=EL%ramp%table(0)%bn(1)*COS(twopi*EL%ramp%table(0)%an(1)*T/clight+EL%ramp%table(0)%bn(2))+EL%ramp%table(0)%an(2)
+! ELP%VOLT=EL%ramp%table(0)%bn(1)*COS(twopi*EL%ramp%table(0)%an(1)*T/clight+EL%ramp%table(0)%bn(2))+EL%ramp%table(0)%an(2)
+! write(6,*) " volt ",el%volt,EL%ramp%table(0)%bn(1)
+! endif
+! 
+! if(EL%ramp%table(0)%b_t/=0.0_dp) then
+! if(EL%parent_fibre%PATCH%TIME==0) EL%parent_fibre%PATCH%TIME=2
+! if(EL%parent_fibre%PATCH%TIME==1) EL%parent_fibre%PATCH%TIME=3
+! EL%parent_fibre%PATCH%b_T=EL%ramp%table(0)%b_t
+! else
+! if(EL%parent_fibre%PATCH%TIME==2) EL%parent_fibre%PATCH%TIME=0
+! if(EL%parent_fibre%PATCH%TIME==3) EL%parent_fibre%PATCH%TIME=1
+! EL%parent_fibre%PATCH%b_T=0.0_dp
+ ! endif
 
- else
+ ! else
 
  IF(EL%P%NMUL>=1.and.associated(EL%D_BN)) THEN
  if(present(VR))then
  do n=1,EL%P%NMUL
  EL%BN(N)= vR*EL%D0_BN(N)+DVR*EL%D_BN(N) 
  EL%AN(N)= vR*EL%D0_AN(N)+DVR*EL%D_AN(N)
- ELP%BN(N)= vR*EL%D0_BN(N)+DVR*EL%D_BN(N)
- ELP%AN(N)= vR*EL%D0_AN(N)+DVR*EL%D_AN(N)
+ ELP%BN(N)= vR*ELp%D0_BN(N)+DVR*ELp%D_BN(N)
+ ELP%AN(N)= vR*ELp%D0_AN(N)+DVR*ELp%D_AN(N)
  enddo
  else
  do n=1,EL%P%NMUL
  EL%BN(N)= vp*EL%D0_BN(N)+DVp*EL%D_BN(N)
  EL%AN(N)= vp*EL%D0_AN(N)+DVp*EL%D_AN(N)
- ELP%BN(N)= vp*EL%D0_BN(N)+DVp*EL%D_BN(N)
- ELP%AN(N)= vp*EL%D0_AN(N)+DVp*EL%D_AN(N)
+ ELP%BN(N)= vp*ELp%D0_BN(N)+DVp*ELp%D_BN(N)
+ ELP%AN(N)= vp*ELp%D0_AN(N)+DVp*ELp%D_AN(N)
  enddo
  endif
 
 
- endif 
+ ! endif 
  if(associated(el%volt)) then
  if(present(VR))then
  EL%volt= vR*EL%D0_Volt+DVR*EL%D_Volt
- ELP%volt= vR*EL%D0_Volt+DVR*EL%D_Volt
+ ELP%volt= vR*ELp%D0_Volt+DVR*ELp%D_Volt
  else
  EL%volt= vp*EL%D0_Volt+DVp*EL%D_Volt
- ELP%volt= vp*EL%D0_Volt+DVp*EL%D_Volt
+ ELP%volt= vp*ELp%D0_Volt+DVp*ELp%D_Volt
  endif
  endif
  if(associated(el%phas)) then
  if(present(VR))then
  EL%phas= vR*EL%D0_phas+DVR*EL%D_phas
- ELP%phas= vR*EL%D0_phas+DVR*EL%D_phas
+ ELP%phas= vR*ELp%D0_phas+DVR*ELp%D_phas
  else
  EL%phas= vp*EL%D0_phas+DVp*EL%D_phas
- ELP%phas= vp*EL%D0_phas+DVp*EL%D_phas
+ ELP%phas= vp*ELp%D0_phas+DVp*ELp%D_phas
  endif
  endif
 

forest/code/Sr_spin.f90

@@ -1572,8 +1572,10 @@ SUBROUTINE TRACK_NODE_FLAG_probe_QUAP(C,XS,K)
  ! endif
 
  IF(K%MODULATION.and.xs%nac/=0) then
+ knob=.true.
  if(c%parent_fibre%mag%slow_ac/=0) CALL MODULATE(C,XS,K) !modulate
  CALL TRACK_MODULATION(C,XS,K) !modulate
+ knob=.false.
  ENDIF !modulate
 
 
@@ -1939,8 +1941,10 @@ SUBROUTINE TRACK_NODE_FLAG_probe_P(C,XS,K)
  ! endif
 
  IF(K%MODULATION.and.xs%nac/=0) then
+ knob=.true.
  if(c%parent_fibre%mag%slow_ac/=0) CALL MODULATE(C,XS,K) !modulate
  CALL TRACK_MODULATION(C,XS,K) !modulate
+ knob=.false.
  ENDIF !modulate
 
 

forest/code/Su_duan_zhe_map.f90

@@ -1073,18 +1073,19 @@ end SUBROUTINE orthonormalise
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! new zhe tracking !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
- SUBROUTINE track_TREE_probe_complex_zhe(T,xs,spin,rad,stoch,linear,slim)
+ SUBROUTINE track_TREE_probe_complex_zhe(T,xs,spin,rad,stoch,linear,slim,flucfac)
 ! use da_arrays
  IMPLICIT NONE
  TYPE(TREE_ELEMENT),target, INTENT(INout) :: T(3)
 
  type(probe) xs
  real(dp) x(size_tree),x0(size_tree),s0(3,3),r(3,3),dx6,beta,q(3),p(3),qg(3),qf(3)
- real(dp) normb,norm,x0_begin(size_tree),xr(6),normbb 
+ real(dp) normb,norm,x0_begin(size_tree),xr(6),normbb,flucfac1
  integer i,j,k,ier,is
  logical, optional :: spin,stoch,rad,linear,slim
  logical spin0,stoch0,rad0,doit,as_is0,slim0
  integer no1
+ real(dp), optional :: flucfac
  type(quaternion) qu
  as_is0=t(1)%usenonsymp
  spin0=.true.
@@ -1099,7 +1100,8 @@ SUBROUTINE track_TREE_probe_complex_zhe(T,xs,spin,rad,stoch,linear,slim)
  if(present(spin)) spin0=spin
  if(present(stoch)) stoch0=stoch
  if(present(rad)) rad0=rad
-
+ flucfac1=1
+ if(present(flucfac)) flucfac1=flucfac
  if(as_is0) rad0=.true.
  doit=rad0.or.stoch0
  x=0
@@ -1117,11 +1119,11 @@ SUBROUTINE track_TREE_probe_complex_zhe(T,xs,spin,rad,stoch,linear,slim)
  xr=0.0_dp
  if(use_gaussian_zhe) then
  do i=1,6
- xr(i)=GRNF_zhe_gaussian()*t(2)%fix0(i) 
+ xr(i)=flucfac1*GRNF_zhe_gaussian()*t(2)%fix0(i) 
  enddo
  else
  do i=1,6
- xr(i)=GRNF_zhe()*t(2)%fix0(i) 
+ xr(i)=flucfac1*GRNF_zhe()*t(2)%fix0(i) 
  enddo
  endif
  xr =matmul(t(2)%rad,xr)

forest/code/a_scratch_size.f90

@@ -226,6 +226,8 @@ module precision_constants
  character(18) :: format8="(8(1x,g23.16,1x))"
  character(18) :: format9="(9(1x,g23.16,1x))"
  character(19) :: format10="(10(1x,g23.16,1x))"
+ character(19) :: format11="(11(1x,g23.16,1x))"
+ character(19) :: format12="(12(1x,g23.16,1x))"
  ! logical(lp) :: fixed_found
  ! lielib_print(1)=1 lieinit prints info
  ! lielib_print(2)=1 expflo warning if no convergence
@@ -246,7 +248,16 @@ module precision_constants
  ! lielib_print(17)=1 print magnets with excessive cutting
  integer , target :: old_integrator =1 ! before making spin high order
  character*255 :: preffile="pref.txt"
-
+ logical :: newtpsa=.false., assume_da_map = .true.,previous_newtpsa=.false.,use_np=.true.,newspin=.true.
+ integer :: with_para=2
+ integer :: nphere=0
+ integer(4), pointer :: inds(:,:),ind1(:),ind2(:) => null()
+ integer(4), pointer :: nind1(:),nind2(:) => null()
+ integer(4) combien, poscombien,ninds
+ real(dp), pointer:: reel(:),finds(:,:,:) => null(),finds1(:) => null()
+ complex(dp), pointer:: reelc(:) => null()
+ integer :: ipric=0,ipric2=0
+
  INTERFACE read
  MODULE PROCEDURE read_d
  MODULE PROCEDURE read_int