module ecloudmod !use bmad use bmad_struct implicit none type :: DiscreetEcloud ! In descriptions, [x] means x is one of the parameters adding a dimension real(rp), allocatable :: exfield(:,:,:) ! 3D snapshot array of 2D [[Ex fields]] for each [timeslice] (t:x:y) real(rp), allocatable :: eyfield(:,:,:) ! 3D snapshot array of 2D Ey fields ... (t:x:y) real(rp), allocatable :: qdenstimes(:) ! 1D array of all [times] used by qdens real(rp), allocatable :: qdens(:) ! 1D array of [electron cloud density] real(rp), allocatable :: timeslices(:) ! 1D array of [timeslices] in a bunch real(rp), allocatable :: bunchtimes(:) ! 1D array of [times] at which each bunch started real(rp), allocatable :: xgrid(:) ! 1D array of the grid values on [x] real(rp), allocatable :: ygrid(:) ! 1D array for [y] end type interface allocatecopy module procedure allocatecopy1, allocatecopy2, allocatecopy3, allocatecopy4 end interface contains subroutine calculateEM(ele, time, orb, field) implicit none type(ele_struct) :: ele type(coord_struct), intent(in) :: orb type(DiscreetEcloud) :: ecloud_data type(em_field_struct) :: field real(rp), intent(in) :: time logical :: err_flag end subroutine !========================================================================! ! Calculates the E fields on the xy plane from the DiscreetEcloud struct ! !========================================================================! subroutine calculateEMFromDiscreet_MultiBunch(ele, time, orb, ecloud_data, field, err_flag) implicit none ! === SUBROUTINE I/O === type(ele_struct) :: ele type(coord_struct), intent(in) :: orb type(DiscreetEcloud) :: ecloud_data type(em_field_struct) :: field real(rp), intent(in) :: time logical :: err_flag ! === ============== === ! === PARAMETERS FOR RETRIEVING INFORMATION === integer, parameter :: reset_value$ = custom_attribute2$ integer, parameter :: multi_val$ = custom_attribute3$ integer, parameter :: first_particle$ = custom_attribute4$ integer, parameter :: latest_particle$ = custom_attribute5$ ! === ===================================== === ! Calculated special values real(rp) :: train_time real(rp) :: bunch_time real(rp) :: qdens real(rp) :: extra_qdens ! interpolated qdensity real(rp) :: qdens_length_multiplier ! From decreasing the size of the element ! === INPUT ALIASES === real(rp) :: x,y,Ex,Ey ! Input and output temps ! === ============= === ! Helper variables integer :: bunch_index ! Setting up ! [x xp y yp z zp] x = orb%vec(1) y = orb%vec(3) ! L != 1 qdens_length_multiplier = ele%value(multi_val$) !==========================! ! Find closest t and qdens ! !==========================! ! Find the closest t and approx. qdens ! === FIND TIME SINCE HEAD OF TRAIN === ! The element will save the first particle time (all particles sorted) and ! compare everything to that. If no particles arrive after [RESET] time then ! if will reset and take a new value next time a particle enters. ! First: Has enough time elapsed to consider this a new train head? if (time-ele%value(latest_particle$) > ele%value(reset_value$)) then ele%value(first_particle$) = time end if ! Now, time since head of train is current time minus time of first particle train_time = time-ele%value(first_particle$) ! === FIND QDENS === extra_qdens = interpolateQdens(ecloud_data, train_time) * qdens_length_multiplier ! TODO end subroutine !========================================================================! ! Calculates the E fields on the xy plane from the DiscreetEcloud struct ! !========================================================================! subroutine calculateEMFromDiscreet_SingleBunch(ele, time, orb, ecloud_data, field, err_flag) implicit none ! === SUBROUTINE I/O === type(ele_struct) :: ele type(coord_struct), intent(in) :: orb type(DiscreetEcloud) :: ecloud_data type(em_field_struct) :: field real(rp), intent(in) :: time logical :: err_flag ! === ============== === ! === PARAMETERS FOR RETRIEVING INFORMATION === integer, parameter :: reset_value$ = custom_attribute2$ integer, parameter :: multi_val$ = custom_attribute3$ integer, parameter :: first_particle$ = custom_attribute4$ integer, parameter :: latest_particle$ = custom_attribute5$ ! === ===================================== === ! === INPUT ALIASES === real(rp) :: x,y,Ex,Ey ! Input and output temps ! === ============= === ! Calculated special values real(rp) :: train_time real(rp) :: bunch_time real(rp) :: qdens real(rp) :: extra_qdens ! interpolated qdensity real(rp) :: qdens_length_multiplier ! From decreasing the size of the element ! Helper variables integer :: bunch_index ! Setting up ! [x xp y yp z zp] x = orb%vec(1) y = orb%vec(3) ! L != 1 qdens_length_multiplier = ele%value(multi_val$) !==========================! ! Find closest t and qdens ! !==========================! ! Find the closest t and approx. qdens ! === FIND TIME SINCE HEAD OF TRAIN === ! The element will save the first particle time (all particles sorted) and ! compare everything to that. If no particles arrive after [RESET] time then ! if will reset and take a new value next time a particle enters. ! First: Has enough time elapsed to consider this a new train head? if (time-ele%value(latest_particle$) > ele%value(reset_value$)) then ele%value(first_particle$) = time end if ! Now, time since head of train is current time minus time of first particle train_time = time-ele%value(first_particle$) ! === FIND QDENS === extra_qdens = interpolateQdens(ecloud_data, train_time) * qdens_length_multiplier ! === FIND CLOSEST BUNCH HEAD TIME === ! Determine bunch index from time ! The particle really shouldn't be before bunch 1... bunch_index = lowbound( ecloud_data%bunchtimes, train_time) if (bunch_index==-1) then print *, " ", "[!] Bunch index: ", bunch_index, "Particle time: ", train_time print *, " ", "INFORMATION:" print *, " ", " ", "Particle time: ", time print *, " ", " ", "Particle vec: ", orb%vec err_flag = .true. return end if bunch_time = train_time-ecloud_data%bunchtimes(bunch_index) !print *, "Time, BunchIndex, BunchHead", train_time, bunch_index, ecloud_data%bunchtimes(bunch_index) ! This should average stuff call interpolateField(ecloud_data, x, y, bunch_time, Ex, Ey) !print *, "t, Ex, Ey, qdens", time, train_time, Ex, Ey, interpolateQdens(ecloud_data, train_time) field%E(1) = (Ex * extra_qdens) field%E(2) = (Ey * extra_qdens) ! WRITE OUTPUT FOR EACH PARTICLE ! write(15,'(9es12.4)') train_time, orb%vec(1), orb%vec(3), orb%vec(5), interpolateQdens(ecloud_data, train_time), Ex, Ey, field%E(1), field%E(2) end subroutine !==============================! ! Extends slope and returns cy ! !==============================! ! Assumes no stupid values (i.e ax=bx) real function cy(ax,ay,bx,by,cx) real(rp) :: ax, ay, bx, by, cx real(rp) :: m if (bx == ax) then cy = (ay+by)/2 else m = (by-ay)/(bx-ax) cy = ((cx-ax)*m) + ay end if return end function !=====================! ! Interpolates points ! !=====================! ! Uses spline method to interpolate points ! ax, ay, ka = x, y, gradient of point a ! bx, by, kb = x, y, gradient of point b ! cx = x position of desired point c real function splineCy(ax, ay, ka, bx, by, kb, cx) implicit None real(rp) :: ax, ay, bx, by, cx, ka, kb real(rp) :: a, b a = ka*( bx - ax ) - ( by - ay ) b = kb*( ax - bx ) - ( ay - by ) splineCy = (1-cx)*ay + cx*by + cx*(1-cx)*( a*(1-cx) + b*cx ) return end function !====================================! ! Finds lowbound of bounding indices ! !====================================! ! This assumes array is sorted integer function lowbound(array, val) real(rp) :: val real(rp) :: array(:) integer :: i = 0 if ( val < array(1) ) then lowbound = 1 return else if (array(size(array)) <= val) then lowbound = size(array) -1 return else do i=1,size(array)-1 if ( array(i) <= val .and. val < array(i+1)) then lowbound = i return end if end do end if ! Error state lowbound = -1 return end function !========================================================! ! Finds the Qdens expected at time given Discreet struct ! !========================================================! real(rp) function interpolateQdens(ecloud_data, train_time) type(DiscreetEcloud) :: ecloud_data real(rp) :: train_time, extra_qdens integer :: i=0 ! Extrapolare qdens extra_qdens = 0 if ( train_time[i]) ! Extend slope on either side of train_time extra_qdens = cy( ecloud_data%qdenstimes(i),ecloud_data%qdens(i),ecloud_data%qdenstimes(i+1),ecloud_data%qdens(i+1),train_time) exit end if end do end if if (extra_qdens > 0) then interpolateQdens = extra_qdens else interpolateQdens = 0 end if return end function !========================================================! ! Finds the Field expected at time given Discreet struct ! !========================================================! subroutine interpolateField(ecloud_data, x, y, bunch_time, Ex, Ey) use bmad use bmad_struct implicit none type(DiscreetEcloud), intent(in) :: ecloud_data type(em_field_struct) :: field integer :: time_low, x_low, y_low real(rp) bunch_time, x, y, qdens ! Cube Extrapolators real(rp) :: ta, tb real(rp) :: xa, xb real(rp) :: ya, yb real(rp) :: Exa,Exb,Exc,Exd, Eya,Eyb,Eyc,Eyd ! Extra between ta and tb real(rp) :: Exe,Exf, Eye,Eyf ! Extra between xa and xb real(rp) :: Ex, Ey ! === FIND CLOSEST BUNCH TIMESLICE === time_low = lowbound( ecloud_data%timeslices, bunch_time ) ! We should now have bunch_index and time_low set to 0 or more ! representing the bunch index and the lower index for the particle ! time in the bunch. !========================! ! Linearly interpolate E ! !========================! ! First find the bounds around which the values will be interpolated (low_x[+1] and low_y[+1]) ! We need to find the extrapolation of E between the two bounding time steps, we will call the ! 8 values E[x/y][a-d] x_low = lowbound(ecloud_data%xgrid, x) y_low = lowbound(ecloud_data%ygrid, y) ta = ecloud_data%timeslices(time_low) tb = ecloud_data%timeslices(time_low + 1) xa = ecloud_data%xgrid(x_low) xb = ecloud_data%xgrid(x_low + 1) ya = ecloud_data%ygrid(y_low) yb = ecloud_data%ygrid(y_low + 1) !==================================! ! * -- * * -- * a -- b ! ! | | -> | | => | | ! ! * -- * * -- * d -- c ! !@ ta tb b_t ! !==================================! Exa = cy( ta, ecloud_data%exfield(time_low,x_low,y_low), tb, ecloud_data%exfield(time_low+1,x_low,y_low), bunch_time ) Exb = cy( ta, ecloud_data%exfield(time_low,x_low+1,y_low), tb, ecloud_data%exfield(time_low+1,x_low+1,y_low), bunch_time ) Exc = cy( ta, ecloud_data%exfield(time_low,x_low+1,y_low+1), tb, ecloud_data%exfield(time_low+1,x_low+1,y_low+1), bunch_time ) Exd = cy( ta, ecloud_data%exfield(time_low,x_low,y_low+1), tb, ecloud_data%exfield(time_low+1,x_low,y_low+1), bunch_time ) Eya = cy( ta, ecloud_data%eyfield(time_low,x_low,y_low), tb, ecloud_data%eyfield(time_low+1,x_low,y_low), bunch_time ) Eyb = cy( ta, ecloud_data%eyfield(time_low,x_low+1,y_low), tb, ecloud_data%eyfield(time_low+1,x_low+1,y_low), bunch_time ) Eyc = cy( ta, ecloud_data%eyfield(time_low,x_low+1,y_low+1), tb, ecloud_data%eyfield(time_low+1,x_low+1,y_low+1), bunch_time ) Eyd = cy( ta, ecloud_data%eyfield(time_low,x_low,y_low+1), tb, ecloud_data%eyfield(time_low+1,x_low,y_low+1), bunch_time ) !==================! ! a b e ! ! | -> | => | ! ! d c f ! !@ xa xb x ! !==================! Exe = cy( xa, Exa, xb, Exb, x) Exf = cy( xa, Exd, xb, Exc, x) Eye = cy( xa, Eya, xb, Eyb, x) Eyf = cy( xa, Eyd, xb, Eyc, x) !==================! ! e -> f => * ! !@ ya yb y ! !==================! Ex = cy( ya, Exe, yb, Exf, y ) Ey = cy( ya, Eye, yb, Eyf, y ) end subroutine subroutine performTests(ele, ecloud_data) use bmad use bmad_struct implicit none integer, parameter :: calculation_index$ = custom_attribute1$ integer, parameter :: reset_time$ = custom_attribute2$ type(ele_struct) :: ele type(DiscreetEcloud) :: ecloud_data integer :: lowb = 0 real(rp) :: exX, exVal character(len=16) :: usr_in print *, "[*] === STARTING TESTS ===" ! Info print *, "[+] INFORMATION:" print *, " ", "[*] QDENS(TIMES) length: ", size(ecloud_data%qdens),",", size(ecloud_data%qdenstimes) print *, " ", "[*] (X/Y) gridsize: ", size(ecloud_data%xgrid),",", size(ecloud_data%ygrid) print *, " ", "[*] E(X/Y) field dim: ", shape(ecloud_data%exfield),",", shape(ecloud_data%eyfield) print *, " ", "[*] TIMESLICES: ", size(ecloud_data%timeslices) print *, " ", "[*] BUNCHTIMES: ", size(ecloud_data%bunchtimes) print *, " ", "[*] CALC_METHOD: ", int(ele%value(calculation_index$)) print *, " ", "[*] TRAIN RESET TIME: ", ele%value(reset_time$) print *, "[+] FIRST TESTS: QDENS, LOWBOUND, CY" exX = 0.0 lowb = lowbound(ecloud_data%qdenstimes, exX) exVal = cy(ecloud_data%qdenstimes(lowb),ecloud_data%qdens(lowb),ecloud_data%qdenstimes(lowb+1),ecloud_data%qdens(lowb+1),exX) print *, " ", "[*] Qdens at time: ", exX, " is ", exVal print *, " ", "[*] lowb: ", lowb print *, " ", "[*] From inputs: ",ecloud_data%qdenstimes(lowb),ecloud_data%qdens(lowb),& ecloud_data%qdenstimes(lowb+1),ecloud_data%qdens(lowb+1),exX exX = 3.32e-11 ! Right between 2 and 3 lowb = lowbound(ecloud_data%qdenstimes, exX) exVal = cy(ecloud_data%qdenstimes(lowb),ecloud_data%qdens(lowb),ecloud_data%qdenstimes(lowb+1),ecloud_data%qdens(lowb+1),exX) print *, " ", "[*] Qdens at time: ", exX, " is ", exVal print *, " ", "[*] lowb: ", lowb print *, " ", "[*] From inputs: ",ecloud_data%qdenstimes(lowb),ecloud_data%qdens(lowb),& ecloud_data%qdenstimes(lowb+1),ecloud_data%qdens(lowb+1),exX print *, "!=====================!" print *, "! INTERACTIVE TESTING !" print *, "!=====================!" print *, "! Type 'y' to enter interactive testing" print *, "! or anything else to quit" read(*,*) usr_in if (usr_in == "y") then call interactiveTesting(ecloud_data) end if end subroutine subroutine interactiveTesting(ecloud_data) use bmad use bmad_struct implicit none type(DiscreetEcloud) :: ecloud_data character(len=4) :: mode real(rp) :: valmin, valmax, valstep real(rp) :: x, y, t, Ex, Ey, qdens print *, "!=========================================!" print *, "! STARTING INTERACTIVE TESTING PROMPT !" print *, "! SET MODE TO 'HELP' FOR MORE INFORMATION !" print *, "!=========================================!" print *, " " print *, "[*] Enter: mode" interactive_loop: do read(*,*) mode if (mode=="exit") then return end if if (mode=="help") then print *, "Help:" print *, " Interactive testing allows you to view columns of data" print *, " from the extrapolation output. Enter a mode to begin." print *, " The mode defines what variable will vary while the" print *, " other two stay the same. The console will output the" print *, " E field at each value of [mode] between the minimum" print *, " and maximum value, incrementing by the step" print *, " Mode can be 'x', 'y', or 't'" else if (mode == "x") then print *, "[*] Enter: y, t, x_min, x_max, x_step" read(*,*) y, t, valmin, valmax, valstep x = valmin do while (x<=valmax) call interpolateField(ecloud_data, x, y, t, Ex, Ey) qdens = interpolateQdens(ecloud_data, t) print *, "[*] x:",x,"Ex:", Ex * qdens, "Ey:", Ey * qdens, "Qdens:", qdens x = x + valstep end do else if (mode == "y") then print *, "[*] Enter: x, t, y_min, y_max, y_step" read(*,*) x, t, valmin, valmax, valstep y = valmin do while (y<=valmax) call interpolateField(ecloud_data, x, y, t, Ex, Ey) qdens = interpolateQdens(ecloud_data, t) print *, "[*] y:",y,"Ex:", Ex * qdens, "Ey:", Ey * qdens, "Qdens:", qdens y = y + valstep end do else if (mode == "t") then print *, "[*] Enter: x, y, t_min, t_max, t_step" read(*,*) x, y, valmin, valmax, valstep t = valmin do while (t<=valmax) call interpolateField(ecloud_data, x, y, t, Ex, Ey) qdens = interpolateQdens(ecloud_data, t) print *, "[*] t:",t,"Ex:", Ex * qdens, "Ey:", Ey * qdens, "Qdens:", qdens t = t + valstep end do else print *, "Unknown mode, try 'help'" end if print *, "!==============!" print *, "[*] Enter mode:" end do interactive_loop end subroutine subroutine calculateDiscreetFromBeamfieldFile(fileunit, ecloud_data) implicit none integer, intent(in) :: fileunit type (DiscreetEcloud) :: ecloud_data end subroutine ! Closes files by itself subroutine calculateDiscreetFromSnapshotsFiles(snapfileunit, densfileunit, bunchunit, ecloud_data) implicit none integer, intent(in) :: snapfileunit, densfileunit, bunchunit integer :: RetCode character(len=500) :: line type (DiscreetEcloud) :: ecloud_data logical :: alt_header, alt_size, bunch_num_head integer :: xgridsize, ygridsize, nslices, qdens_num, num_bunches real(rp), allocatable :: xgrid(:) real(rp), allocatable :: ygrid(:) real(rp), allocatable :: timeslices(:) ! Slices in bunch real(rp), allocatable :: times(:) ! All times in qdens real(rp), allocatable :: bunchheads(:) ! times bunches start real(rp), allocatable :: ex(:,:,:) ! (t:x:y) real(rp), allocatable :: ey(:,:,:) ! (t:x:y) real(rp), allocatable :: qdensity(:) ! Helpers integer :: t = 1 integer :: y = 1 alt_header = .true. alt_size = .true. bunch_num_head = .true. !======================================! ! First read in the header information ! !======================================! header_loop: do read(snapfileunit,'(A)',iostat=RetCode) line line = trim(line) if (RetCode > 0) then print *,"[!] Major error in reading snapshot(header) file" exit header_loop else if (RetCode < 0) then print *,"[!] Reached end of snapshot file" print *,"[!] ERROR: End of file while reading header: ", RetCode exit header_loop end if if (index(line,'!')/=0) then ! Skip comment lines (with ! in them) cycle header_loop end if ! First read the gridsize if (alt_header) then if (alt_size) then read(line,*) xgridsize ygridsize = xgridsize print *, " ", "[+] expected xgridsize: ", xgridsize print *, " ", "[+] expected ygridsize: ", ygridsize alt_size = .false. ! Next time read the nslices else read(line,*) nslices print *, " ", "[+] nslices: ", nslices alt_header = .false. ! Next time read the grid information end if else ! Then read the x and y grids if (allocated(xgrid)) then! do y allocate( ygrid(ygridsize) ) read(line,*) ygrid print *, " ", "[+] ygridsize: ", size(ygrid) exit header_loop else ! do x allocate( xgrid(xgridsize) ) read(line,*) xgrid print *, " ", "[+] xgridsize: ", size(xgrid) ! Since allocated, next time do above end if end if end do header_loop ! Now I should have xgrid, ygrid, and nslices allocate( timeslices(nslices) ) allocate( ex(nslices,xgridsize,ygridsize) ) allocate( ey(nslices,xgridsize,ygridsize) ) !======================! ! Now read in the data ! !======================! read_loop: do ! First skip all the comments up to the time declaration read(snapfileunit,'(A)',iostat=RetCode) line line = trim(line) if (RetCode > 0) then print *,"[!] Major error in reading snapshot file" exit read_loop else if (RetCode < 0) then print *,"[+] Reached end of snapshot file with no problems" exit read_loop end if if (index(line,'!')/=0) then ! Skip comment lines (with ! in them) cycle read_loop end if ! Now read in the time read(line,*) timeslices(t) print *, " ", "[+] Reading E field for time: ", timeslices(t) ! And the ex do y=1,ygridsize read(snapfileunit,*) ex(t,:,y) end do ! And the ey do y=1,ygridsize read(snapfileunit,*) ey(t,:,y) end do t = t+1 end do read_loop ! Now I should have timeslices, ex and ey ! So we're done with the snapshot file close(snapfileunit) ! And now qdens :D !============================! ! Now read the cloud density ! !============================! ! First read header info qdens_loop: do read(densfileunit,'(A)',iostat=RetCode) line if (RetCode > 0) then print *,"[!] Major error in reading qdens file" exit qdens_loop else if (RetCode < 0) then print *,"[!] Reached end of qdens file" print *,"[!] ERROR: end of file in header" exit qdens_loop end if if (index(line,'!')/=0) then cycle qdens_loop end if ! Not a comment read(line,*) qdens_num exit qdens_loop end do qdens_loop ! Now allocate arrays allocate(qdensity(qdens_num)) allocate(times(qdens_num)) t = 1 ! Like do t=1, qdens_num but supports cycles qdens_dens: do while (t<=qdens_num) read(densfileunit,'(A)',iostat=RetCode) line if (RetCode > 0) then print *,"[!] Major error in reading qdens file" exit qdens_dens else if (RetCode < 0) then print *,"[!] Reached end of qdens file" print *,"[!] ERROR: Should not have reached this point without filling array" exit qdens_dens end if if (index(line,'!')/=0) then cycle qdens_dens end if read(line,*) times(t), qdensity(t) t = t+1 end do qdens_dens close(densfileunit) !========================! ! Read bunch information ! !========================!` bunchheads_do: do read(bunchunit,'(A)',iostat=RetCode) line if (RetCode > 0) then print *,"[!] Major error in reading qdens file" exit bunchheads_do else if (RetCode < 0) then print *,"[+] Reached end of bunchheads file with no problems" exit bunchheads_do end if if (index(line,'!')/=0) then cycle bunchheads_do else ! Else not necessary due to cycle if (bunch_num_head) then read(line, *) num_bunches print *, " ", "[+] number of bunches: ", num_bunches allocate(bunchheads(num_bunches)) bunch_num_head = .false. else read(line,*) bunchheads exit bunchheads_do end if end if end do bunchheads_do close(bunchunit) !============================! ! Now we have all the arrays ! !============================! ! Allocate with source is unavailable with gfortran ! allocate(ecloud_data%exfield, source=ex) call allocateCopy(ecloud_data%exfield, ex) call allocateCopy(ecloud_data%eyfield, ey) call allocateCopy(ecloud_data%qdens, qdensity) call allocateCopy(ecloud_data%xgrid, xgrid) call allocateCopy(ecloud_data%ygrid, ygrid) call allocateCopy(ecloud_data%timeslices, timeslices) call allocateCopy(ecloud_data%qdenstimes, times) call allocateCopy(ecloud_data%bunchtimes, bunchheads) end subroutine subroutine calculateDiscreetFromQdensFile(fileunit, ecloud_data) use bmad use bmad_struct implicit none integer, intent(in) :: fileunit type (DiscreetEcloud) :: ecloud_data end subroutine subroutine allocatecopy1(target, src) real(rp), allocatable :: target(:) real(rp), allocatable :: src(:) allocate(target(size(src))) target = src end subroutine subroutine allocatecopy2(target, src) real(rp), allocatable :: target(:,:) real(rp), allocatable :: src(:,:) allocate(target(size(src, 1), size(src, 2))) target = src end subroutine subroutine allocatecopy3(target, src) real(rp), allocatable :: target(:,:,:) real(rp), allocatable :: src(:,:,:) allocate(target(size(src, 1), size(src, 2), size(src, 3))) target = src end subroutine subroutine allocatecopy4(target, src) real(rp), allocatable :: target(:,:,:,:) real(rp), allocatable :: src(:,:,:,:) allocate(target(size(src, 1), size(src, 2), size(src, 3), size(src, 4))) target = src end subroutine end module