Long v2
MODULE long_v2
! Copyright Mats Lindroos and Steve Hancock, CERN, Switzerland, 1998, v1
! Copyright Mats Lindroos and Steve Hancock, CERN, Switzerland, 2000, v2
USE nrstuff
TYPE jlim !Type declaration for upper limits of j (buildmaps)
INTEGER u
INTEGER l
END TYPE jlim
!-----------------------------------------
! There are two reference systems: i) the reconstructed phase space
! coordinate system which has its origin at some minimum energy and minimum
! phase, with x-bins dtbin wide and y-bins dEbin wide; and ii) the physical
! coordinates expressed in energy and rf phase wrt the synchronous particle.
!
! xat0: synchronous phase in bins in beam_ref_frame
! yat0: synchronous energy (0 in relative terms) in reconstructed phase
! space coordinate system
! rebin: rebinning factor
! dturns: number of machine turns between each measurement
! preskiplength: subtract this number of bins from the beginning of
! the 'raw' input traces
! postskiplength: subtract this number of bins from the end of
! the 'raw' input traces
! skipcount: skip this number of traces from the beginning of the 'raw'
! input file
! framecount: number of traces in the 'raw' input file
! framelength: length of each trace in the 'raw' input file
! alldata: total number of datapoints in the 'raw' input file
! allbinmin,allbinmax: imin and imax as calculated from jmax and jmin
! Npt: Npt**2 is the number of test particles tracked from each pixel of
! reconstructed phase space
! niter: number of iterations in reconstruction process
! reconstruct_p: profile at which phase space is reconstructed
! machine_ref_frame: frame to which machine parameters are referenced (B0,VRF1,VRF2)
! beam_ref_frame: frame to which beam parameters are referenced (baseline, phi0)
! filmstep: step between consecutive reconstructions for the profiles from
! filmstart to filmstop
!-----------------------------------------
REAL(SP) xat0,yat0,xorigin
INTEGER preskiplength,postskiplength,skipcount
INTEGER framecount,framelength
INTEGER profilecount,profilelength
INTEGER alldata,iminin,imaxin
INTEGER, DIMENSION(:,:), ALLOCATABLE :: jmax,jmin
INTEGER, DIMENSION(:), ALLOCATABLE :: imin,imax,allbinmin,allbinmax
INTEGER Npt,niter,reconstruct_p,machine_ref_frame,beam_ref_frame
INTEGER filmstart,filmstop,filmstep,rebin,dturns,self_field_flag
INTERFACE trajectoryheight
MODULE PROCEDURE trajectoryheight_var,trajectoryheight_arr
END INTERFACE
INTERFACE rfvolt
MODULE PROCEDURE rfvolt_r,rfvolt_arr
END INTERFACE
INTERFACE realft
MODULE PROCEDURE realft_sp,realft_dp
END INTERFACE
INTERFACE four1
MODULE PROCEDURE four1_sp,four1_dp
END INTERFACE
INTERFACE fourrow
MODULE PROCEDURE fourrow_sp,fourrow_dp
END INTERFACE
!-----------------------------------------
! B0: B-field at machine_ref_frame
! Bdot: time derivative of B-field (considered constant)
! Rnom: mean orbit radius
! rhonom: bending radius
! gammatnom: transition gamma
! Erest: rest energy of accelerated particle
! q: charge state of accelerated particle
! c: speed of light in vacuum
! eunit: elementary electric charge
! h: principal harmonic number
! hratio: ratio of harmonics between the two RF systems
! phi12: phase difference between the two RF systems (considered constant)
! self_field_flag: flag to include self-fields in the tracking
! zpu: effective pick-up sensitivity
! gcoupling: space charge coupling coefficient
! zwallovern: magnitude of Zwall/n
! eperimage: total charge in beam_ref_profile
! vself: space charge voltage
! VRF1: peak voltage of first RF system at machine_ref_frame
! VRF2: peak voltage of second RF system at machine_ref_frame
! VRF1dot,VRF2dot: time derivatives of the RF voltages (considered constant)
! dtbin: pixel width
! dEbin: pixel height
! dEmax: maximum energy of reconstructed phase space
! E0: total energy of synchronous particle at each turn
! beta0: Lorenz beta factor (v/c) at each turn
! eta0: phase slip factor at each turn
! omegarev0: revolution frequency at each turn
! tatturn: time at each turn relative to machine_ref_frame
! phi0: synchronous phase angle at each turn
! fitxat0: synchronous phase in bins from foot tangent fit
! proffile: file with 'raw' input frames
! odir: directory into which the ouput is directed
! full_pp_flag: if set, all pixels in reconstructed phase space will be tracked
! wraplength: maximum number of bins to cover an integer number of rf periods
!-----------------------------------------
REAL(SP), PRIVATE :: h, Bdot, Rnom, rhonom, gammatnom, Erest, q, eperimage
REAL(SP), PRIVATE :: VRF1,VRF2,eunit,fitxat0,bunchphaselength
REAL(SP), PRIVATE :: VRF1dot,VRF2dot,zwallovern,gcoupling
REAL(SP), PRIVATE :: dtbin,dEbin,energy,B0,c,hratio,phi12,dEmax,zpu
REAL(SP), PRIVATE :: tangentfootl,tangentfootu,phiwrap
REAL(SP), DIMENSION(:), ALLOCATABLE, PRIVATE :: beta0,eta0,&
omegarev0,tatturn,phi0,E0,c1,c2,c3
REAL(SP), DIMENSION(:,:), ALLOCATABLE :: vself
INTEGER, PRIVATE :: full_pp_flag,wraplength
CHARACTER*60, PRIVATE :: proffile,odir
CONTAINS
SUBROUTINE find_xat0(profile)
! Subroutine to find the synchronous phase from the beam reference profile.
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(IN) :: profile
INTEGER i,i01,tangentbinl,tangentbinu,turnnow
REAL(SP) al,bl,au,bu,chil,chiu,ql,qu,sigal,sigbl,sigau,sigbu
REAL(SP) thold,bunchduration,phil,dradbin
REAL(SP), DIMENSION(profilelength) :: indarr
INTEGER, DIMENSION(1) :: maxbin
!
IF (xat0.LT.0) THEN
DO i=1,profilelength
indarr(i)=REAL(i,SP)-0.5_SP
END DO
turnnow=(beam_ref_frame-1)*dturns
thold=0.15_SP*MAXVAL(profile)
maxbin=MAXLOC(profile)
DO i=maxbin(1),1,-1
IF (profile(i).LT.thold) THEN
tangentbinl=i+1
EXIT
END IF
END DO
DO i=maxbin(1),profilelength,1
IF (profile(i).LT.thold) THEN
tangentbinu=i-1
EXIT
END IF
END DO
CALL fit(indarr(tangentbinl-2:tangentbinl+1),&
profile(tangentbinl-2:tangentbinl+1),al,bl,sigal,sigbl,chil,ql)
CALL fit(indarr(tangentbinu-1:tangentbinu+2),&
profile(tangentbinu-1:tangentbinu+2),au,bu,sigau,sigbu,chiu,qu)
tangentfootl=-1.0_SP*al/bl
tangentfootu=-1.0_SP*au/bu
bunchduration=(tangentfootu-tangentfootl)*dtbin
bunchphaselength=h*omegarev0(turnnow)*bunchduration
phil=rtnewt(phi0(turnnow)-bunchphaselength,phi0(turnnow),&
phi0(turnnow)-bunchphaselength/2_SP,0.0001_SP,phaselow,turnnow)
fitxat0=tangentfootl+(phi0(turnnow)-phil)/(h*omegarev0(turnnow)*dtbin)
xat0=fitxat0
ELSE
tangentfootl=0.0_SP
tangentfootu=0.0_SP
fitxat0=0.0_SP
END IF
! Calculate the absolute difference in bins between phase=0 and origin of
! the reconstructed phase space coordinate system.
i01=(beam_ref_frame-1)*dturns
xorigin=phi0(i01)/(h*omegarev0(i01)*dtbin)-xat0
! Calculate the number of bins in the first integer number of rf periods
! larger than the image width.
IF (Bdot .GT. 0.0_sp) THEN
dradbin=h*omegarev0((profilecount-1)*dturns)*dtbin
ELSE
dradbin=h*omegarev0(0)*dtbin
END IF
phiwrap=REAL(CEILING(profilelength*dradbin/twopi),sp)*twopi
wraplength=CEILING(phiwrap/dradbin)
END SUBROUTINE find_xat0
SUBROUTINE phaselow(phase,hamil,dhamil,rfturn)
!-----------------------------------------------
! This function is needed by the Newton-Raphson
! root finder to calculate fitxat0.
!-----------------------------------------------
USE nrstuff
IMPLICIT NONE
REAL(SP), INTENT(IN) :: phase
REAL(SP), INTENT(OUT) :: hamil,dhamil
INTEGER, INTENT(IN) :: rfturn
hamil=VRF2*(COS(hratio*(phase+bunchphaselength-phi12))-&
COS(hratio*(phase-phi12)))/hratio+&
VRF1*(COS(phase+bunchphaselength)-COS(phase))+&
bunchphaselength*rfvolt(phi0(rfturn),rfturn)
dhamil=-1.0_SP*VRF2*(SIN(hratio*(phase+bunchphaselength-phi12))-&
SIN(hratio*(phase-phi12)))-&
VRF1*(SIN(phase+bunchphaselength)-SIN(phase))
END SUBROUTINE phaselow
SUBROUTINE perimage(profile)
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(IN) :: profile
!
! Calculate the total charge in profile.
eperimage=SUM(profile)*dtbin/(rebin*eunit*zpu)
END SUBROUTINE perimage
SUBROUTINE calculate_self(dsprofiles)
IMPLICIT NONE
REAL(SP), DIMENSION(:,:), INTENT(IN) :: dsprofiles
INTEGER p
ALLOCATE(vself(profilecount-1,wraplength))
!
! Calculate self-field voltage.
vself=0.0_sp
DO p=1,profilecount-1
vself(p,1:profilelength)=0.5_sp*eperimage*(&
c3(p)*dsprofiles(p,1:profilelength)+&
c3(p+1)*dsprofiles(p+1,1:profilelength))
END DO
CALL out_profiles(vself(:,1:profilelength),'vself.data')
END SUBROUTINE calculate_self
SUBROUTINE subtract_baseline(indata)
! Find the baseline from the first 5% percent of the beam reference profile.
IMPLICIT NONE
INTEGER barray
REAL(SP), DIMENSION(:), INTENT(INOUT) :: indata
REAL(SP) baseline
!
! Skip some of the frame data to the start of the reference profile.
barray=(skipcount+beam_ref_frame-1)*framelength + 1 + preskiplength
!
baseline=SUM(indata(barray:&
FLOOR(barray+0.05*REAL(profilelength,SP))))&
/REAL(FLOOR(0.05*REAL(profilelength,SP)+1.0_SP),SP)
! Offset by baseline.
indata=indata-baseline
END SUBROUTINE subtract_baseline
SUBROUTINE get_indata(indata)
USE nrstuff
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(INOUT) :: indata
INTEGER i
OPEN(1,file=proffile,status='old',access='sequential',&
form='formatted')
DO i=1,alldata
read(1,*) indata(i)
END DO
CLOSE(1)
END SUBROUTINE get_indata
SUBROUTINE get_parameters()
!-----------------------------------------
! Read input data and initialise parameters.
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
CHARACTER*60 comment
INTEGER allturn,i,iminskip,imaxskip,turnhere
OPEN(9,file="input_v2.dat",status='old',&
access='sequential',form='formatted')
DO i=1,11
READ(9,'(a60)') comment
END DO
!
READ(9,'(a60)') comment
READ(9,'(a60)') proffile
!
READ(9,'(a60)') comment
READ(9,'(a60)') odir
!
READ(9,'(a60)') comment
READ(9,*) framecount
!
READ(9,'(a60)') comment
READ(9,*) skipcount
!
READ(9,'(a60)') comment
READ(9,*) framelength
!
READ(9,'(a60)') comment
READ(9,*) dtbin
!
READ(9,'(a60)') comment
READ(9,*) dturns
!
READ(9,'(a60)') comment
READ(9,*) preskiplength
!
READ(9,'(a60)') comment
READ(9,*) postskiplength
!
READ(9,'(a60)') comment
READ(9,'(a60)') comment
READ(9,*) iminskip
!
READ(9,'(a60)') comment
READ(9,'(a60)') comment
READ(9,*) imaxskip
!
READ(9,'(a60)') comment
READ(9,*) rebin
!
READ(9,'(a60)') comment
READ(9,'(a60)') comment
READ(9,*) xat0
!
READ(9,'(a60)') comment
READ(9,*) dEmax
!
READ(9,'(a60)') comment
READ(9,*) filmstart
!
READ(9,'(a60)') comment
READ(9,*) filmstop
!
READ(9,'(a60)') comment
READ(9,*) filmstep
!
READ(9,'(a60)') comment
READ(9,*) niter
!
READ(9,'(a60)') comment
READ(9,*) Npt
!
READ(9,'(a60)') comment
READ(9,*) full_pp_flag
!
READ(9,'(a60)') comment
READ(9,*) beam_ref_frame
!
READ(9,'(a60)') comment
READ(9,*) machine_ref_frame
!
READ(9,'(a60)') comment
READ(9,'(a60)') comment
!
READ(9,'(a60)') comment
READ(9,*) VRF1
!
READ(9,'(a60)') comment
READ(9,*) VRF1dot
!
READ(9,'(a60)') comment
READ(9,*) VRF2
!
READ(9,'(a60)') comment
READ(9,*) VRF2dot
!
READ(9,'(a60)') comment
READ(9,*) h
!
READ(9,'(a60)') comment
READ(9,*) hratio
!
READ(9,'(a60)') comment
READ(9,*) phi12
!
READ(9,'(a60)') comment
READ(9,*) B0
!
READ(9,'(a60)') comment
READ(9,*) Bdot
!
READ(9,'(a60)') comment
READ(9,*) Rnom
!
READ(9,'(a60)') comment
READ(9,*) rhonom
!
READ(9,'(a60)') comment
READ(9,*) gammatnom
!
READ(9,'(a60)') comment
READ(9,*) Erest
!
READ(9,'(a60)') comment
READ(9,*) q
!
READ(9,'(a60)') comment
READ(9,'(a60)') comment
READ(9,'(a60)') comment
READ(9,*) self_field_flag
!
READ(9,'(a60)') comment
READ(9,*) gcoupling
!
READ(9,'(a60)') comment
READ(9,*) zwallovern
!
READ(9,'(a60)') comment
READ(9,*) zpu
!
CLOSE(9)
!
allturn=(framecount-skipcount-1)*dturns
ALLOCATE(tatturn(0:allturn),omegarev0(0:allturn),phi0(0:allturn),&
c1(0:allturn),c2(0:allturn),&
beta0(0:allturn),eta0(0:allturn),E0(0:allturn))
! Speed of light in vacuum.
c=2.99792458e8_SP
! Elementary electric charge.
eunit=1.60217733e-19_sp
! Lorenz beta (relativistic factor).
CALL pbttrack(phi0,beta0,tatturn,E0,Erest,allturn)
! Phase slip factor.
eta0=(1.0_SP-beta0**2)-gammatnom**(-2)
c1=twopi*h*eta0/(E0*beta0**2)
! Revolution frequency.
omegarev0=beta0*c/Rnom
! Time binning.
dtbin=REAL(rebin,SP)*dtbin
xat0=xat0/REAL(rebin,SP)
profilecount=framecount-skipcount
! Here, profilelength is still in frame bins.
profilelength=framelength-preskiplength-postskiplength
iminin=iminskip/rebin + 1
IF (MOD(profilelength-imaxskip,rebin).EQ.0) THEN
imaxin=(profilelength-imaxskip)/rebin
ELSE
imaxin=(profilelength-imaxskip)/rebin + 1
END IF
alldata=framecount*framelength
! Self-field coefficient.
ALLOCATE(c3(profilecount))
DO i=1,profilecount
turnhere=(i-1)*dturns
c3(i)=(eunit/omegarev0(turnhere))*&
((1_sp/beta0(turnhere)-beta0(turnhere))*&
gcoupling*pi*2.e-7_sp*c-zwallovern)/dtbin**2
END DO
END SUBROUTINE get_parameters
SUBROUTINE pbttrack(phi0,beta0,tatturn,E0,Erest,allturn)
!-----------------------------------------
! Subroutine which calculates synchronous phase, Lorenz beta, energy and time
! for the synchronous particle as a function of turn. Values are calculated
! immediately after the 'single' cavity of the ring.
!
! Erest: rest energy of accelerated particle
! E0: total energy at the end of the turn
! beta0: Lorenz beta at the end of the turn
! phi0: synchronous phase at the end of the turn
! tatturn: time relative to machine_ref_frame at the end of the turn
! allturn: total number of turns
! i,i0: integer variables for loops
!-----------------------------------------
IMPLICIT NONE
REAL(SP), INTENT(IN) :: Erest
REAL(SP) phistart,phil,phiu
INTEGER allturn,i,i0
REAL(SP), DIMENSION(0:allturn), INTENT(OUT) :: phi0,beta0,tatturn,E0
i0=(machine_ref_frame-1)*dturns
tatturn(i0)=0
E0(i0)=BtoEnom(B0)
beta0(i0)=SQRT(1.0_SP-(Erest/E0(i0))**2)
IF (q*(gammatnom-E0(i0)/Erest)>0) THEN
phil=-1.0_SP*pi
phiu=pi
ELSE
phil=0.0_SP
phiu=twopi
END IF
phistart=rtnewt(phil,phiu,(phil+phiu)/2.0_SP,0.001_SP,rfvrf1only,i0)
! Synchronous phase of a particle on the nominal orbit.
phi0(i0)=rtnewt(phil,phiu,phistart,0.001_SP,rfvanddv,i0)
DO i=i0+1,allturn
tatturn(i)=tatturn(i-1)+twopi*Rnom/(beta0(i-1)*c)
phi0(i)=rtnewt(phil,phiu,phi0(i-1),0.001_SP,rfvanddv,i)
E0(i)=E0(i-1)+q*rfvolt(phi0(i),i)
beta0(i)=SQRT(1.0_SP-(Erest/E0(i))**2)
c2(i)=E0(i)-E0(i-1)
END DO
DO i=i0-1,0,-1
E0(i)=E0(i+1)-q*rfvolt(phi0(i+1),i+1)
beta0(i)=SQRT(1.0_SP-(Erest/E0(i))**2)
c2(i)=E0(i+1)-E0(i)
tatturn(i)=tatturn(i+1)-twopi*Rnom/(beta0(i)*c)
phi0(i)=rtnewt(phil,phiu,phi0(i+1),0.001_SP,rfvanddv,i)
END DO
END SUBROUTINE pbttrack
SUBROUTINE out_plotinfo
!-----------------------------------------
! Writes data needed for plots on a file.
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
INTEGER p
OPEN(8,file=TRIM(odir)//'plotinfo.data',status='new',access='sequential',&
form='formatted')
write(8,'(a)') 'Number of profiles used in each reconstruction,'
write(8,'(a,I3)') ' profilecount = ',profilecount
write(8,'(a)') 'Width (in pixels) of each image = length (in bins) of each profile,'
write(8,'(a,I3)') ' profilelength = ',profilelength
write(8,'(a)') 'Width (in s) of each pixel = width of each profile bin,'
write(8,'(a,1P,E11.4)') ' dtbin =',dtbin
write(8,'(a)') 'Height (in eV) of each pixel,'
write(8,'(a,1P,E11.4)') ' dEbin =',dEbin
write(8,'(a)') 'Number of elementary charges in each image,'
write(8,'(a,1P,E10.3)') ' eperimage =',eperimage
write(8,'(a)') 'Position (in pixels) of the reference synchronous point:'
write(8,'(a,F8.3)') ' xat0 =',xat0
write(8,'(a,F8.3)') ' yat0 =',yat0
write(8,'(a)') 'Foot tangent fit results (in bins):'
write(8,'(a,F8.3)') ' tangentfootl = ',tangentfootl
write(8,'(a,F8.3)') ' tangentfootu = ',tangentfootu
write(8,'(a,F8.3)') ' fit xat0 =',fitxat0
write(8,'(a)') 'Synchronous phase (in radians):'
DO p=filmstart,filmstop,filmstep
write(8,'(a,I3,a,F7.4)') ' phi0(',p,') =',phi0((p-1)*dturns)
END DO
write(8,'(a)') 'Horizontal range (in pixels) of the region in phase space of map elements:'
DO p=filmstart,filmstop,filmstep
write(8,'(a,i3,a,i3,a,i3,a,i3)') ' imin(',p,') = ',imin(p),&
' and imax(',p,') = ',imax(p)
END DO
CLOSE(8)
END SUBROUTINE out_plotinfo
SUBROUTINE out_profiles(profiles,fname)
!-----------------------------------------
! Writes profiles to a file.
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
INTEGER P,I
CHARACTER*(*) :: fname
REAL(SP), DIMENSION(:,:), INTENT(IN) :: profiles
OPEN(8,file=TRIM(odir)//fname,status='new',access='sequential',&
form='formatted')
DO P=1,SIZE(profiles)/profilelength
DO I=1,profilelength
write(8,'(E14.7)') profiles(P,I)
END DO
END DO
CLOSE(8)
END SUBROUTINE out_profiles
SUBROUTINE out_picture(phasespace,fname)
!-----------------------------------------
! Writes phasepace to a file.
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
REAL(SP), DIMENSION(:,:), INTENT(IN) :: phasespace
INTEGER i,j
CHARACTER*(*) :: fname
OPEN(8,file=TRIM(odir)//fname,status='new',access='sequential',&
form='formatted')
DO i=1,profilelength
DO j=1,profilelength
write(8,*) phasespace(i,j)
END DO
END DO
CLOSE(8)
END SUBROUTINE out_picture
SUBROUTINE out_array(onearray,fname)
!-----------------------------------------
! Writes any array to a file.
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(IN) :: onearray
INTEGER i
CHARACTER*(*) :: fname
CLOSE(8)
OPEN(8,file=TRIM(odir)//fname,status='new',access='sequential',&
form='formatted')
DO i=1,SIZE(onearray)
write(8,*) i-1,onearray(i)
END DO
CLOSE(8)
END SUBROUTINE out_array
SUBROUTINE ijlimits
!-----------------------------------------
! This procedure calculates and sets the limits in i (phase) and j (energy).
! The i-j coordinate system is the one used locally (reconstructed phase space).
!
! LOCAL VARIABLES:
! turnnow: machine turn (=0 at profile=1)
! jmaxul: array holding calculated max values (left and right of bin)
! phases: phase at the edge of each bin along the i-axis
! energiesl,energiesu: maximum energy starting from lowest phase and
! uppermost phase respectively
! indarray: local array holding real numbers from 1 to profilelength
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
INTEGER i,p,turnnow
TYPE(jlim), DIMENSION(profilelength+1) :: jmaxul
REAL(SP),DIMENSION(profilelength+1)::phases,energiesl,energiesu
REAL(SP), DIMENSION(profilelength) :: indarr
turnnow=(beam_ref_frame-1)*dturns
DO i=1,profilelength
indarr(i)=REAL(i,SP)
END DO
phases(1)=xorigin*dtbin*h*omegarev0(turnnow)
phases(2:profilelength+1)=(xorigin+indarr)*dtbin*h*omegarev0(turnnow)
! Energy binning.
IF (dEmax.LT.0.0_SP) THEN
IF (VRFt(VRF2,VRF2dot,turnnow).NE.0.0_SP) THEN
energiesl=trajectoryheight(phases,phases(1),0.0_SP,turnnow)
energiesu=trajectoryheight(phases,phases(profilelength+1),&
0.0_SP,turnnow)
dEbin=MIN(MAXVAL(energiesl),MAXVAL(energiesu))/(profilelength-yat0)
ELSE
dEbin=beta0(turnnow)*SQRT(E0(turnnow)*q*VRFt(VRF1,VRF1dot,turnnow)*&
COS(phi0(turnnow))/(twopi*h*eta0(turnnow)))*dtbin*h*omegarev0(turnnow)
END IF
ELSE
dEbin=dEmax/(profilelength-yat0)
ENDIF
! Extrema of active pixels in energy (j-index or y-index) direction.
energy=0.0_SP
IF (full_pp_flag.EQ.1) THEN
jmax=profilelength
jmin=CEILING(2.0_SP*yat0-jmax+0.5_SP)
allbinmin=1
allbinmax=profilelength
ELSE
DO p=filmstart,filmstop,filmstep
turnnow=(p-1)*dturns
phases(1)=xorigin*dtbin*h*omegarev0(turnnow)
phases(2:profilelength+1)=(xorigin+indarr)*dtbin*h*omegarev0(turnnow)
DO i=1,profilelength+1
jmaxul(i)%l=FLOOR(yat0+trajectoryheight(phases(i),phases(1),&
energy,turnnow)/dEbin)
jmaxul(i)%u=FLOOR(yat0+trajectoryheight(phases(i),&
phases(profilelength+1),energy,turnnow)/dEbin)
END DO
DO i=1,profilelength
jmax(p,i)=MIN(jmaxul(i)%u,jmaxul(i+1)%u,jmaxul(i)%l,jmaxul(i+1)%l,&
profilelength)
END DO
jmin(p,:)=MAX(1,CEILING(2.0_SP*yat0-jmax(p,:)+0.5_SP))
DO i=1,profilelength
IF ((jmax(p,i)-jmin(p,i)).GE.0) THEN
allbinmin(p)=i
EXIT
END IF
END DO
DO i=profilelength,1,-1
IF ((jmax(p,i)-jmin(p,i)).GE.0) THEN
allbinmax(p)=i
EXIT
END IF
END DO
END DO
END IF
OPEN(8,file=TRIM(odir)//'jmax.data',status='new',&
access='sequential',form='formatted')
DO p=filmstart,filmstop,filmstep
imin(p)=allbinmin(p)
IF (iminin.GT.allbinmin(p) .OR. full_pp_flag.EQ.1) THEN
imin(p)=iminin
jmax(p,1:iminin-1)=FLOOR(yat0)
jmin(p,:)=CEILING(2.0_SP*yat0-jmax(p,:)+0.5_SP)
END IF
imax(p)=allbinmax(p)
IF (imaxin.LT.allbinmax(p) .OR. full_pp_flag.EQ.1) THEN
imax(p)=imaxin
jmax(p,imaxin+1:profilelength)=FLOOR(yat0)
jmin(p,:)=CEILING(2.0_SP*yat0-jmax(p,:)+0.5_SP)
END IF
DO i=1,profilelength
write(8,*) i,jmax(p,i)
END DO
END DO
CLOSE(8)
DEALLOCATE(E0,eta0,beta0)
END SUBROUTINE ijlimits
FUNCTION trajectoryheight_var(phi,phiknown,deltaEknown,turnnow)
!-----------------------------------------
! This is a trajectory height calculator
! given a phase and energy.
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
COMPLEX complexheight
REAL(SP) trajectoryheight_var,phi,phiknown,deltaEknown,aa,bb,cc,dd
INTEGER turnnow
aa=deltaEknown**2
bb=2.0_SP*q/c1(turnnow)
cc=VRFt(VRF1,VRF1dot,turnnow)*(COS(phi)-COS(phiknown))+&
VRFt(VRF2,VRF2dot,turnnow)*(COS(hratio*(phi-phi12))-&
COS(hratio*(phiknown-phi12)))/hratio+(phi-phiknown)*&
rfvolt(phi0(turnnow),turnnow)
dd=aa+bb*cc
complexheight=SQRT(CMPLX(dd))
trajectoryheight_var=REAL(complexheight,SP)
END FUNCTION
FUNCTION trajectoryheight_arr(phi,phiknown,deltaEknown,turnnow)
!-----------------------------------------
! This is a trajectory height calculator
! given a phase and energy.
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
REAL(SP) phiknown,deltaEknown,aa,bb
INTEGER turnnow
REAL(SP), DIMENSION(:) :: phi
REAL(SP),DIMENSION(SIZE(phi))::trajectoryheight_arr,cc,dd
COMPLEX, DIMENSION(SIZE(phi))::complexheight
aa=deltaEknown**2
bb=2.0_SP*q/c1(turnnow)
cc=VRFt(VRF1,VRF1dot,turnnow)*(COS(phi)-COS(phiknown))+&
VRFt(VRF2,VRF2dot,turnnow)*(COS(hratio*(phi-phi12))-&
COS(hratio*(phiknown-phi12)))/hratio+(phi-phiknown)*&
rfvolt(phi0(turnnow),turnnow)
dd=aa+bb*cc
complexheight=SQRT(CMPLX(dd))
trajectoryheight_arr=REAL(complexheight,SP)
END FUNCTION
FUNCTION BtoEnom(B)
!-----------------------------------------
! Calculates the energy for a particle in
! a circular machine at dipole field B.
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
REAL(SP) BtoEnom,B
BtoEnom=SQRT((q*B*rhonom*c)**2 + Erest**2)
END FUNCTION
SUBROUTINE rfvanddv(phi,rfv,drfv,rfturn)
!-----------------------------------------
! The rfvolt function and its derivative.
! This function is needed by the Newton-Raphson
! root finder to calculate phi0.
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
REAL(SP), INTENT(IN) :: phi
REAL(SP), INTENT(OUT) :: rfv,drfv
INTEGER, INTENT(IN) :: rfturn
REAL(SP) V1,V2
V1=VRFt(VRF1,VRF1dot,rfturn)
V2=VRFt(VRF2,VRF2dot,rfturn)
rfv=V1*SIN(phi)+V2*SIN(hratio*(phi-phi12))-twopi*Rnom*rhonom*Bdot*SIGN(1.0_SP,q)
drfv=V1*COS(phi)+hratio*V2*COS(hratio*(phi-phi12))
END SUBROUTINE rfvanddv
SUBROUTINE rfvrf1only(phi,rfv,drfv,rfturn)
!-----------------------------------------
! The rfvolt function and its derivative.
! This function is needed by the Newton-Raphson
! root finder to calculate phi0.
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
REAL(SP), INTENT(IN) :: phi
REAL(SP), INTENT(OUT) :: rfv,drfv
INTEGER, INTENT(IN) :: rfturn
REAL(SP) V1
V1=VRFt(VRF1,VRF1dot,rfturn)
rfv=V1*SIN(phi)-twopi*Rnom*rhonom*Bdot*SIGN(1.0_SP,q)
drfv=V1*COS(phi)
END SUBROUTINE rfvrf1only
SUBROUTINE longtrack(direction,nrep,yp,xp,turnnow)
!---------------------------------------------------------------------------
! h: principal harmonic number
! eta0: phase slip factor
! E0: energy of synchronous particle
! beta0: relativistic beta of synchronous particle
! phi0: synchronous phase
! q: charge state of particles
! dphi: phase difference between considered particle and synchronous one
! denergy: energy difference between considered particle and synchronous one
! nrep: pass cavity nrep times before returning data
! direction: to inverse the time advance (rotation in the bucket), 1 or -1
! xp and yp: time and energy in pixels
! dtbin and dEbin: GLOBAL time and energy pixel size in s and MeV
! omegarev0: revolution frequency
! VRF1,VRF2,VRF1dot,VRF2dot: GLOBAL RF voltages and derivatives of volts
! turnnow: present turn
!---------------------------------------------------------------------------
USE nrstuff
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(INOUT) :: xp,yp
REAL(SP), DIMENSION(SIZE(xp)) :: dphi,denergy
!HPF$ distribute dphi(block)
!HPF$ align with dphi :: denergy
integer :: mm
INTEGER :: i,nrep,direction,turnnow
dphi=(xp+xorigin)*h*omegarev0(turnnow)*dtbin-phi0(turnnow)
denergy=(yp-yat0)*dEbin
IF (direction.GT.0) THEN
DO i=1,nrep
forall(mm=1:size(xp)) dphi(mm)=dphi(mm)-c1(turnnow)*denergy(mm)
turnnow=turnnow+1
forall(mm=1:size(xp)) denergy(mm)=denergy(mm)+q*(&
(VRF1+VRF1dot*tatturn(turnnow))*SIN(dphi(mm)+phi0(turnnow))+&
(VRF2+VRF2dot*tatturn(turnnow))*&
SIN(hratio*(dphi(mm)+phi0(turnnow)-phi12)))-c2(turnnow)
END DO
ELSE
DO i=1,nrep
forall(mm=1:size(xp)) denergy(mm)=denergy(mm)-q*(&
(VRF1+VRF1dot*tatturn(turnnow))*SIN(dphi(mm)+phi0(turnnow))+&
(VRF2+VRF2dot*tatturn(turnnow))*&
SIN(hratio*(dphi(mm)+phi0(turnnow)-phi12)))+c2(turnnow)
turnnow=turnnow-1
forall(mm=1:size(xp)) dphi(mm)=dphi(mm)+c1(turnnow)*denergy(mm)
END DO
END IF
xp=(dphi+phi0(turnnow))/(h*omegarev0(turnnow)*dtbin)-xorigin
yp=denergy/dEbin+yat0
END SUBROUTINE longtrack
SUBROUTINE longtrack_self(direction,nrep,yp,xp,turnnow)
!-------------------------------------------------------------------------
! h: principal harmonic number
! eta0: phase slip factor
! E0: energy of synchronous particle
! beta0: relativistic beta of synchronous particle
! phi0: synchronous phase
! q: charge state of particles
! dphi: phase difference between considered particle and synchronous one
! denergy: energy difference between considered particle and synchronous one
! nrep: pass cavity nrep times before returning data
! direction: to inverse the time advance (rotation in the bucket), 1 or -1
! xp and yp: time and energy in pixels
! dtbin and dEbin: GLOBAL time and energy pixel size in s and MeV
! omegarev0: revolution frequency
! VRF1,VRF2,VRF1dot,VRF2dot: GLOBAL RF voltages and derivatives of volts
! turnnow: present turn
!---------------------------------------------------------------------------
USE nrstuff
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(INOUT) :: xp,yp
REAL(SP), DIMENSION(SIZE(xp)) :: dphi,denergy,selfvolt
!HPF$ distribute dphi(block)
!HPF$ align with dphi :: denergy,selfvolt,xp
integer :: mm
INTEGER :: i,j,p,nrep,direction,turnnow,xx
dphi=(xp+xorigin)*h*omegarev0(turnnow)*dtbin-phi0(turnnow)
denergy=(yp-yat0)*dEbin
IF (direction.GT.0) THEN
p=turnnow/dturns+1
DO i=1,nrep
forall(mm=1:size(xp)) dphi(mm)=dphi(mm)-c1(turnnow)*denergy(mm)
turnnow=turnnow+1
forall(mm=1:size(xp)) xp(mm)=dphi(mm)+phi0(turnnow)-&
xorigin*h*omegarev0(turnnow)*dtbin
forall(mm=1:size(xp)) xp(mm)=(xp(mm)-&
phiwrap*FLOOR(xp(mm)/phiwrap))/(h*omegarev0(turnnow)*dtbin)
forall(mm=1:size(xp)) selfvolt(mm)=vself(p,FLOOR(xp(mm))+1)
forall(mm=1:size(xp)) denergy(mm)=denergy(mm)+q*((&
(VRF1+VRF1dot*tatturn(turnnow))*SIN(dphi(mm)+phi0(turnnow))+&
(VRF2+VRF2dot*tatturn(turnnow))*&
SIN(hratio*(dphi(mm)+phi0(turnnow)-phi12)))+selfvolt(mm))-c2(turnnow)
END DO
ELSE
p=turnnow/dturns
DO i=1,nrep
forall(mm=1:size(xp)) selfvolt(mm)=vself(p,FLOOR(xp(mm))+1)
forall(mm=1:size(xp)) denergy(mm)=denergy(mm)-q*((&
(VRF1+VRF1dot*tatturn(turnnow))*SIN(dphi(mm)+phi0(turnnow))+&
(VRF2+VRF2dot*tatturn(turnnow))*&
SIN(hratio*(dphi(mm)+phi0(turnnow)-phi12)))+selfvolt(mm))+c2(turnnow)
turnnow=turnnow-1
forall(mm=1:size(xp)) dphi(mm)=dphi(mm)+c1(turnnow)*denergy(mm)
forall(mm=1:size(xp)) xp(mm)=dphi(mm)+phi0(turnnow)-&
xorigin*h*omegarev0(turnnow)*dtbin
forall(mm=1:size(xp)) xp(mm)=(xp(mm)-&
phiwrap*FLOOR(xp(mm)/phiwrap))/(h*omegarev0(turnnow)*dtbin)
END DO
END IF
yp=denergy/dEbin+yat0
END SUBROUTINE longtrack_self
FUNCTION rfvolt_r(phi,rfturn)
!-----------------------------------------
! For a new RF function change rfvolt here,
! but pay attention for optmization purpose
! rfvolt is implemeted inline as well.
!-----------------------------------------
USE nrstuff
IMPLICIT NONE
REAL(SP) phi,rfvolt_r
INTEGER rfturn
rfvolt_r=VRFt(VRF1,VRF1dot,rfturn)*SIN(phi)+&
VRFt(VRF2,VRF2dot,rfturn)*SIN(hratio*(phi-phi12))
END FUNCTION
FUNCTION rfvolt_arr(phi,rfturn)
! For a new RF function change rfvolt here.
USE nrstuff
IMPLICIT NONE
INTEGER rfturn
REAL(SP), DIMENSION(:) :: phi
REAL(SP), DIMENSION(SIZE(phi)) :: rfvolt_arr
rfvolt_arr=VRFt(VRF1,VRF1dot,rfturn)*SIN(phi)+&
VRFt(VRF2,VRF2dot,rfturn)*SIN(hratio*(phi-phi12))
END FUNCTION
FUNCTION rtnewt(x1,x2,xstart,xacc,funcd,rfturn)
!--------------------------------------------
! This is a Newton-Raphson root finder.
! Modified from Numerical Recipies to start at
! xstart and calculate rfvolt at rfturn.
!--------------------------------------------
USE nrstuff
IMPLICIT NONE
REAL(SP) x1,x2,xstart,xacc
REAL(SP) :: rtnewt
INTEGER rfturn
INTERFACE
SUBROUTINE funcd(x,fval,fderiv,rfturn)
USE nrstuff
IMPLICIT NONE
REAL(SP), INTENT(IN) :: x
REAL(SP), INTENT(OUT) :: fval,fderiv
INTEGER, INTENT(IN) :: rfturn
END SUBROUTINE funcd
END INTERFACE
INTEGER(I4B), PARAMETER :: MAXIT=20
INTEGER(I4B) :: j
REAL(SP) :: df,dx,f
rtnewt=xstart
do j=1,MAXIT
call funcd(rtnewt,f,df,rfturn)
dx=f/df
rtnewt=rtnewt-dx
if ((x1-rtnewt)*(rtnewt-x2) < 0.0)&
call nrerror('rtnewt: values jumped out of brackets')
if (abs(dx) < xacc) RETURN
end do
call nrerror('rtnewt exceeded maximum iterations')
END FUNCTION rtnewt
FUNCTION VRFt(VRF,VRFdot,rfturn)
!--------------------------------------------
! This function calculates the RF peak voltage at turn rfturn
! assuming a linear voltage function VRFt=VRFdot*time+VRF.
! time=0 at machine_ref_frame.
!
! rfturn: turn for which the RF voltage should be calculated
! VRF: Volts
! VRFdot: Volts/s
!--------------------------------------------
REAL(SP) VRFt,VRF,VRFdot
INTEGER rfturn
VRFt=VRF+VRFdot*tatturn(rfturn)
END FUNCTION VRFt
SUBROUTINE fit(x,y,a,b,siga,sigb,chi2,q,sig)
USE nrstuff
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(IN) :: x,y
REAL(SP), INTENT(OUT) :: a,b,siga,sigb,chi2,q
REAL(SP), DIMENSION(:), OPTIONAL, INTENT(IN) :: sig
INTEGER(I4B) :: ndum
REAL(SP) :: sigdat,ss,sx,sxoss,sy,st2
REAL(SP), DIMENSION(size(x)), TARGET :: t
REAL(SP), DIMENSION(:), POINTER :: wt
if (present(sig)) then
ndum=assert_eq(size(x),size(y),size(sig),'fit')
wt=>t
wt(:)=1.0_sp/(sig(:)**2)
ss=sum(wt(:))
sx=dot_product(wt,x)
sy=dot_product(wt,y)
else
ndum=assert_eq(size(x),size(y),'fit')
ss=real(size(x),sp)
sx=sum(x)
sy=sum(y)
end if
sxoss=sx/ss
t(:)=x(:)-sxoss
if (present(sig)) then
t(:)=t(:)/sig(:)
b=dot_product(t/sig,y)
else
b=dot_product(t,y)
end if
st2=dot_product(t,t)
b=b/st2
a=(sy-sx*b)/ss
siga=sqrt((1.0_sp+sx*sx/(ss*st2))/ss)
sigb=sqrt(1.0_sp/st2)
t(:)=y(:)-a-b*x(:)
if (present(sig)) then
t(:)=t(:)/sig(:)
chi2=dot_product(t,t)
q=gammq(0.5_sp*(size(x)-2),0.5_sp*chi2)
else
chi2=dot_product(t,t)
q=1.0
sigdat=sqrt(chi2/(size(x)-2))
siga=siga*sigdat
sigb=sigb*sigdat
end if
END SUBROUTINE fit
FUNCTION savgol(nl,nr,ld,m)
IMPLICIT NONE
INTEGER(I4B), INTENT(IN) :: nl,nr,ld,m
REAL(SP), DIMENSION(nl+nr+1) :: savgol
INTEGER(I4B) :: imj,ipj,mm,np
INTEGER(I4B), DIMENSION(m+1) :: indx
REAL(SP) :: d,sm
REAL(SP), DIMENSION(m+1) :: b
REAL(SP), DIMENSION(m+1,m+1) :: a
INTEGER(I4B) :: irng(nl+nr+1)
call assert(nl >= 0, nr >= 0, ld <= m, nl+nr >= m, 'Bad savgol args!')
do ipj=0,2*m
sm=sum(arth(1.0_sp,1.0_sp,nr)**ipj)+&
sum(arth(-1.0_sp,-1.0_sp,nl)**ipj)
if (ipj == 0) sm=sm+1.0_sp
mm=min(ipj,2*m-ipj)
do imj=-mm,mm,2
a(1+(ipj+imj)/2,1+(ipj-imj)/2)=sm
end do
end do
call ludcmp(a(:,:),indx(:),d)
b(:)=0.0
b(ld+1)=1.0
call lubksb(a(:,:),indx(:),b(:))
savgol(:)=0.0
np=nl+nr+1
irng(:)=arth(-nl,1,np)
savgol(mod(np-irng(:),np)+1)=poly(real(irng(:),sp),b(:))
END FUNCTION savgol
SUBROUTINE lubksb(a,indx,b)
IMPLICIT NONE
REAL(SP), DIMENSION(:,:), INTENT(IN) :: a
INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indx
REAL(SP), DIMENSION(:), INTENT(INOUT) :: b
INTEGER(I4B) :: i,n,ii,ll
REAL(SP) :: summ
n=assert_eq(size(a,1),size(a,2),size(indx),'lubksb')
ii=0
do i=1,n
ll=indx(i)
summ=b(ll)
b(ll)=b(i)
if (ii /= 0) then
summ=summ-dot_product(a(i,ii:i-1),b(ii:i-1))
else if (summ /= 0.0) then
ii=i
end if
b(i)=summ
end do
do i=n,1,-1
b(i) = (b(i)-dot_product(a(i,i+1:n),b(i+1:n)))/a(i,i)
end do
END SUBROUTINE lubksb
SUBROUTINE ludcmp(a,indx,d)
IMPLICIT NONE
REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a
INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: indx
REAL(SP), INTENT(OUT) :: d
REAL(SP), DIMENSION(size(a,1)) :: vv
REAL(SP), PARAMETER :: TINY=1.0e-20_sp
INTEGER(I4B) :: j,n,imax
n=assert_eq(size(a,1),size(a,2),size(indx),'ludcmp')
d=1.0
vv=maxval(abs(a),dim=2)
if (any(vv == 0.0)) call nrerror('Singular matrix in ludcmp!')
vv=1.0_sp/vv
do j=1,n
imax=(j-1)+imaxloc(vv(j:n)*abs(a(j:n,j)))
if (j /= imax) then
call swap(a(imax,:),a(j,:))
d=-d
vv(imax)=vv(j)
end if
indx(j)=imax
if (a(j,j) == 0.0) a(j,j)=TINY
a(j+1:n,j)=a(j+1:n,j)/a(j,j)
a(j+1:n,j+1:n)=a(j+1:n,j+1:n)-outerprod(a(j+1:n,j),a(j,j+1:n))
end do
END SUBROUTINE ludcmp
FUNCTION convlv(data,respns,isign)
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(INOUT) :: data
REAL(SP), DIMENSION(:), INTENT(IN) :: respns
INTEGER(I4B), INTENT(IN) :: isign
REAL(SP), DIMENSION(size(data)) :: convlv
INTEGER(I4B) :: no2,n,m
COMPLEX(SPC), DIMENSION(size(data)/2) :: tmpd,tmpr
n=size(data)
m=size(respns)
call assert(iand(n,n-1)==0, 'n must be a power of 2 in convlv!')
call assert(mod(m,2)==1, 'm must be odd in convlv!')
convlv(1:m)=respns(:)
convlv(n-(m-3)/2:n)=convlv((m+3)/2:m)
convlv((m+3)/2:n-(m-1)/2)=0.0
no2=n/2
call realft(data,1,tmpd)
call realft(convlv,1,tmpr)
if (isign == 1) then
tmpr(1)=cmplx(real(tmpd(1))*real(tmpr(1))/no2, &
aimag(tmpd(1))*aimag(tmpr(1))/no2, kind=spc)
tmpr(2:)=tmpd(2:)*tmpr(2:)/no2
else if (isign == -1) then
if (any(abs(tmpr(2:)) == 0.0) .or. real(tmpr(1)) == 0.0 &
.or. aimag(tmpr(1)) == 0.0) call nrerror &
('deconvolving at response zero in convlv')
tmpr(1)=cmplx(real(tmpd(1))/real(tmpr(1))/no2, &
aimag(tmpd(1))/aimag(tmpr(1))/no2, kind=spc)
tmpr(2:)=tmpd(2:)/tmpr(2:)/no2
else
call nrerror('No meaning for isign in convlv!')
end if
call realft(convlv,-1,tmpr)
END FUNCTION convlv
SUBROUTINE realft_sp(data,isign,zdata)
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(INOUT) :: data
INTEGER(I4B), INTENT(IN) :: isign
COMPLEX(SPC), DIMENSION(:), OPTIONAL, TARGET :: zdata
INTEGER(I4B) :: n,ndum,nh,nq
COMPLEX(SPC), DIMENSION(size(data)/4) :: w
COMPLEX(SPC), DIMENSION(size(data)/4-1) :: h1,h2
COMPLEX(SPC), DIMENSION(:), POINTER :: cdata
COMPLEX(SPC) :: z
REAL(SP) :: c1=0.5_sp,c2
n=size(data)
call assert(iand(n,n-1)==0, 'n must be a power of 2 in realft_sp!')
nh=n/2
nq=n/4
if (present(zdata)) then
ndum=assert_eq(n/2,size(zdata),'realft_sp')
cdata=>zdata
if (isign == 1) cdata=cmplx(data(1:n-1:2),data(2:n:2),kind=spc)
else
allocate(cdata(n/2))
cdata=cmplx(data(1:n-1:2),data(2:n:2),kind=spc)
end if
if (isign == 1) then
c2=-0.5_sp
call four1(cdata,+1)
else
c2=0.5_sp
end if
w=zroots_unity(sign(n,isign),n/4)
w=cmplx(-aimag(w),real(w),kind=spc)
h1=c1*(cdata(2:nq)+conjg(cdata(nh:nq+2:-1)))
h2=c2*(cdata(2:nq)-conjg(cdata(nh:nq+2:-1)))
cdata(2:nq)=h1+w(2:nq)*h2
cdata(nh:nq+2:-1)=conjg(h1-w(2:nq)*h2)
z=cdata(1)
if (isign == 1) then
cdata(1)=cmplx(real(z)+aimag(z),real(z)-aimag(z),kind=spc)
else
cdata(1)=cmplx(c1*(real(z)+aimag(z)),c1*(real(z)-aimag(z)),kind=spc)
call four1(cdata,-1)
end if
if (present(zdata)) then
if (isign /= 1) then
data(1:n-1:2)=real(cdata)
data(2:n:2)=aimag(cdata)
end if
else
data(1:n-1:2)=real(cdata)
data(2:n:2)=aimag(cdata)
deallocate(cdata)
end if
END SUBROUTINE realft_sp
SUBROUTINE realft_dp(data,isign,zdata)
IMPLICIT NONE
REAL(DP), DIMENSION(:), INTENT(INOUT) :: data
INTEGER(I4B), INTENT(IN) :: isign
COMPLEX(DPC), DIMENSION(:), OPTIONAL, TARGET :: zdata
INTEGER(I4B) :: n,ndum,nh,nq
COMPLEX(DPC), DIMENSION(size(data)/4) :: w
COMPLEX(DPC), DIMENSION(size(data)/4-1) :: h1,h2
COMPLEX(DPC), DIMENSION(:), POINTER :: cdata
COMPLEX(DPC) :: z
REAL(DP) :: c1=0.5_dp,c2
n=size(data)
call assert(iand(n,n-1)==0, 'n must be a power of 2 in realft_dp!')
nh=n/2
nq=n/4
if (present(zdata)) then
ndum=assert_eq(n/2,size(zdata),'realft_dp')
cdata=>zdata
if (isign == 1) cdata=cmplx(data(1:n-1:2),data(2:n:2),kind=spc)
else
allocate(cdata(n/2))
cdata=cmplx(data(1:n-1:2),data(2:n:2),kind=spc)
end if
if (isign == 1) then
c2=-0.5_dp
call four1(cdata,+1)
else
c2=0.5_dp
end if
w=zroots_unity(sign(n,isign),n/4)
w=cmplx(-aimag(w),real(w),kind=dpc)
h1=c1*(cdata(2:nq)+conjg(cdata(nh:nq+2:-1)))
h2=c2*(cdata(2:nq)-conjg(cdata(nh:nq+2:-1)))
cdata(2:nq)=h1+w(2:nq)*h2
cdata(nh:nq+2:-1)=conjg(h1-w(2:nq)*h2)
z=cdata(1)
if (isign == 1) then
cdata(1)=cmplx(real(z)+aimag(z),real(z)-aimag(z),kind=dpc)
else
cdata(1)=cmplx(c1*(real(z)+aimag(z)),c1*(real(z)-aimag(z)),kind=dpc)
call four1(cdata,-1)
end if
if (present(zdata)) then
if (isign /= 1) then
data(1:n-1:2)=real(cdata)
data(2:n:2)=aimag(cdata)
end if
else
data(1:n-1:2)=real(cdata)
data(2:n:2)=aimag(cdata)
deallocate(cdata)
end if
END SUBROUTINE realft_dp
SUBROUTINE four1_sp(data,isign)
IMPLICIT NONE
COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: data
INTEGER(I4B), INTENT(IN) :: isign
COMPLEX(SPC), DIMENSION(:,:), ALLOCATABLE :: dat,temp
COMPLEX(DPC), DIMENSION(:), ALLOCATABLE :: w,wp
REAL(DP), DIMENSION(:), ALLOCATABLE :: theta
INTEGER(I4B) :: n,m1,m2,j
n=size(data)
call assert(iand(n,n-1)==0, 'n must be a power of 2 in four1_sp!')
m1=2**ceiling(0.5_sp*log(real(n,sp))/log(2.0_sp))
m2=n/m1
allocate(dat(m1,m2),theta(m1),w(m1),wp(m1),temp(m2,m1))
dat=reshape(data,shape(dat))
call fourrow(dat,isign)
theta=arth(0,isign,m1)*TWOPI_D/n
wp=cmplx(-2.0_dp*sin(0.5_dp*theta)**2,sin(theta),kind=dpc)
w=cmplx(1.0_dp,0.0_dp,kind=dpc)
do j=2,m2
w=w*wp+w
dat(:,j)=dat(:,j)*w
end do
temp=transpose(dat)
call fourrow(temp,isign)
data=reshape(temp,shape(data))
deallocate(dat,w,wp,theta,temp)
END SUBROUTINE four1_sp
SUBROUTINE four1_dp(data,isign)
IMPLICIT NONE
COMPLEX(DPC), DIMENSION(:), INTENT(INOUT) :: data
INTEGER(I4B), INTENT(IN) :: isign
COMPLEX(DPC), DIMENSION(:,:), ALLOCATABLE :: dat,temp
COMPLEX(DPC), DIMENSION(:), ALLOCATABLE :: w,wp
REAL(DP), DIMENSION(:), ALLOCATABLE :: theta
INTEGER(I4B) :: n,m1,m2,j
n=size(data)
call assert(iand(n,n-1)==0, 'n must be a power of 2 in four1_dp!')
m1=2**ceiling(0.5_sp*log(real(n,sp))/log(2.0_sp))
m2=n/m1
allocate(dat(m1,m2),theta(m1),w(m1),wp(m1),temp(m2,m1))
dat=reshape(data,shape(dat))
call fourrow(dat,isign)
theta=arth(0,isign,m1)*TWOPI_D/n
wp=cmplx(-2.0_dp*sin(0.5_dp*theta)**2,sin(theta),kind=dpc)
w=cmplx(1.0_dp,0.0_dp,kind=dpc)
do j=2,m2
w=w*wp+w
dat(:,j)=dat(:,j)*w
end do
temp=transpose(dat)
call fourrow(temp,isign)
data=reshape(temp,shape(data))
deallocate(dat,w,wp,theta,temp)
END SUBROUTINE four1_dp
SUBROUTINE fourrow_sp(data,isign)
IMPLICIT NONE
COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: data
INTEGER(I4B), INTENT(IN) :: isign
INTEGER(I4B) :: n,i,istep,j,m,mmax,n2
REAL(DP) :: theta
COMPLEX(SPC), DIMENSION(size(data,1)) :: temp
COMPLEX(DPC) :: w,wp
COMPLEX(SPC) :: ws
n=size(data,2)
call assert(iand(n,n-1)==0, 'n must be a power of 2 in fourrow_sp!')
n2=n/2
j=n2
do i=1,n-2
if (j > i) call swap(data(:,j+1),data(:,i+1))
m=n2
do
if (m < 2 .or. j < m) exit
j=j-m
m=m/2
end do
j=j+m
end do
mmax=1
do
if (n <= mmax) exit
istep=2*mmax
theta=PI_D/(isign*mmax)
wp=cmplx(-2.0_dp*sin(0.5_dp*theta)**2,sin(theta),kind=dpc)
w=cmplx(1.0_dp,0.0_dp,kind=dpc)
do m=1,mmax
ws=w
do i=m,n,istep
j=i+mmax
temp=ws*data(:,j)
data(:,j)=data(:,i)-temp
data(:,i)=data(:,i)+temp
end do
w=w*wp+w
end do
mmax=istep
end do
END SUBROUTINE fourrow_sp
SUBROUTINE fourrow_dp(data,isign)
IMPLICIT NONE
COMPLEX(DPC), DIMENSION(:,:), INTENT(INOUT) :: data
INTEGER(I4B), INTENT(IN) :: isign
INTEGER(I4B) :: n,i,istep,j,m,mmax,n2
REAL(DP) :: theta
COMPLEX(DPC), DIMENSION(size(data,1)) :: temp
COMPLEX(DPC) :: w,wp
COMPLEX(DPC) :: ws
n=size(data,2)
call assert(iand(n,n-1)==0, 'n must be a power of 2 in fourrow_dp!')
n2=n/2
j=n2
do i=1,n-2
if (j > i) call swap(data(:,j+1),data(:,i+1))
m=n2
do
if (m < 2 .or. j < m) exit
j=j-m
m=m/2
end do
j=j+m
end do
mmax=1
do
if (n <= mmax) exit
istep=2*mmax
theta=PI_D/(isign*mmax)
wp=cmplx(-2.0_dp*sin(0.5_dp*theta)**2,sin(theta),kind=dpc)
w=cmplx(1.0_dp,0.0_dp,kind=dpc)
do m=1,mmax
ws=w
do i=m,n,istep
j=i+mmax
temp=ws*data(:,j)
data(:,j)=data(:,i)-temp
data(:,i)=data(:,i)+temp
end do
w=w*wp+w
end do
mmax=istep
end do
END SUBROUTINE fourrow_dp
END MODULE long_v2