module ecloudmod use bmad_struct use spline_mod implicit none type :: DiscreteEcloud ! 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 :: 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] real(rp), allocatable :: qdenstimes(:) ! 1D array of all [times] used by qdens type(spline_struct), allocatable :: qdens_spline(:) ! spline_struct for interpolating qdens end type interface allocatecopy module procedure allocatecopy1, allocatecopy2, allocatecopy3, allocatecopy4 end interface logical, save :: on = .true. real(rp), save :: qslope = 3.e-6 real(rp), save :: qcslope = 1.e-5 real(rp), save :: qslopeslope = 0 real(rp), save :: qcslopeslope = 0 integer, save :: start_bunch contains subroutine setstate(state) implicit none logical :: state print *, on, ":=", state on = state end subroutine subroutine setSlopes(s, cs, ss,css) implicit none real(rp) :: s, cs, ss, css qslope = s qcslope = cs qslopeslope = ss qcslopeslope = css end subroutine subroutine setStartBunch(bunchnum) implicit none integer :: bunchnum start_bunch = bunchnum end subroutine subroutine calculateEM(ele, time, orb, field) implicit none type(ele_struct) :: ele type(coord_struct), intent(in) :: orb type(DiscreteEcloud) :: ecloud_data type(em_field_struct) :: field real(rp), intent(in) :: time logical :: err_flag ! === Parameters === integer, parameter :: reset_value$ = custom_attribute2$ integer, parameter :: multi_val$ = custom_attribute3$ integer, parameter :: first_particle$ = custom_attribute4$ integer, parameter :: latest_particle$ = custom_attribute5$ ! === ========== === real(rp) :: qs,qcs, t logical :: time_elapsed = .false. logical :: first_zero = .false. logical :: time_dc = .false. ! Find time since head of train time_elapsed = (time - ele%value(latest_particle$)) .GE. ele%value(reset_value$) ! if enough time has elapsed since the last particle first_zero = ele%value(first_particle$) .LE. 0 ! or if the first particle was at time zero time_dc = ele%value(first_particle$) .GT. time ! or if the first particle is in the future if ( time_elapsed .or. first_zero .or. time_dc ) then ele%value(first_particle$) = time ! then consider this to be the first particle !print *, "[*] New start of train at time: ", time end if ele%value(latest_particle$) = time t = time - ele%value(first_particle$) qs = qslope + (t*qslopeslope) qcs = qcslope + (t*qcslopeslope) field%E(1) = orb%vec(1)*qs + (orb%vec(1)*orb%vec(1)*orb%vec(1))*qcs field%E(2) = orb%vec(3)*qs + (orb%vec(3)*orb%vec(3)*orb%vec(3))*qcs end subroutine !===============================================================! ! Calculates the E fields on the xy plane from the Qdens struct ! !===============================================================! subroutine calculateEMFromDiscrete_QDens(ele, time, orb, ecloud_data, field, err_flag) implicit none ! === SUBROUTINE I/O === type(ele_struct) :: ele type(coord_struct), intent(in) :: orb type(DiscreteEcloud), intent(in) :: 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$ ! === ===================================== === ! === UTILITY VARIABLES === integer :: lun, ios ! === ================= === ! Calculated special values real(rp) :: train_time real(rp) :: extra_qdens ! interpolated qdensity real(rp) :: qdens_length_multiplier ! From decreasing the size of the element ! === INPUT ALIASES === real(rp) :: theta,r,x,y,Ex,Ey ! Input and output temps ! === ============= === logical :: time_elapsed = .false. logical :: first_zero = .false. logical :: time_dc = .false. ! === TEMP SLOPE === ! real(rp) :: slope = 3.e-6 ! This should be filled in by user input or a file later ! real(rp) :: cslope = 1.e-5 ! This is a similar constant for the cubic terms ! === ========== === ! If the lattice element is not on, return 0 immediately if (.not. on) then field%E(1) = 0 field%E(2) = 0 return end if ! 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. ! Find time since head of train time_elapsed = (time - ele%value(latest_particle$)) .GE. ele%value(reset_value$) ! if enough time has elapsed since the last particle first_zero = ele%value(first_particle$) .LE. 0 ! or if the first particle was at time zero time_dc = ele%value(first_particle$) .GT. time ! or if the first particle is in the future if ( time_elapsed .or. first_zero .or. time_dc ) then ele%value(first_particle$) = time ! then consider this to be the first particle !print *, "[*] New start of train at time: ", time end if ! This is now the most recent particle ele%value(latest_particle$) = time ! Now, time since head of train is current time minus time of first particle train_time = time - ele%value(first_particle$) ! Throw an error if the train is longer than allowed if (train_time > ele%value(reset_value$)) then print *, "Train time exceeds reset time" print *, train_time, ">", ele%value(reset_value$) print *, "Train head: ", ele%value(first_particle$) print *, "Latest: ", ele%value(latest_particle$) print *, "Time: ", time !err_flag = .True. end if ! === FIND QDENS === extra_qdens = interpolateQdens(ecloud_data, train_time) * qdens_length_multiplier ! === FIRST ADD LINEAR PARTS === Ex = qslope * x * extra_qdens Ey = qslope * y * extra_qdens !print *, "Fields:", Ex, Ey !print *, "At: ", train_time, extra_qdens ! === CUBIC TERMS === ! Simplified formula: ! Eq = (...) + c*r^3*. = (...) + c * r^2 * q^2 Ex = Ex + (qcslope * x * (x*x + y*y) * extra_qdens* extra_qdens) Ey = Ey + (qcslope * y * (x*x + y*y) * extra_qdens* extra_qdens) if (ISNAN(Ex)) then Ex = 0 end if if (ISNAN(Ey)) then Ey = 0 end if ! Return the fields for this lattice element field%E(1) = Ex field%E(2) = Ey ! Write out debug information !write(1708,'(11es12.4)') time, train_time, 0, & ! 1, 2, 3 ! orb%vec(1), orb%vec(3), orb%vec(5), & ! 4, 5, 6 ! extra_qdens, & ! 7 ! Ex, Ey, field%E(1), field%E(2) ! 8, 9, 10, 11 end subroutine !========================================================================! ! Calculates the E fields on the xy plane from the DiscreteEcloud struct ! !========================================================================! subroutine calculateEMFromDiscrete_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(DiscreteEcloud) :: 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_multiplier ! From decreasing the size of the element ! === INPUT ALIASES === real(rp) :: x,y,Ex,Ey ! Input and output temps ! === ============= === ! === TEMP SLOPE === !real(rp) :: slope = 3170.9 ! This should be filled in by user input or a file later ! === ========== === ! Helper variables integer :: bunch_index logical :: time_elapsed = .false. logical :: first_zero = .false. logical :: time_dc = .false. ! Setting up ! [x xp y yp z zp] x = orb%vec(1) y = orb%vec(3) ! L != 1 qdens_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. ! Find time since head of train time_elapsed = (time - ele%value(latest_particle$)) .GE. ele%value(reset_value$) ! if enough time has elapsed since the last particle first_zero = ele%value(first_particle$) .LE. 0 ! or if the first particle was at time zero time_dc = ele%value(first_particle$) .GT. time ! or if the first particle is in the future if ( time_elapsed .or. first_zero .or. time_dc ) then ele%value(first_particle$) = time ! then consider this to be the first particle !print *, "[*] New start of train at time: ", time end if ! This is now the most recent particle ele%value(latest_particle$) = time ! Now, time since head of train is current time minus time of first particle train_time = time - ele%value(first_particle$) if (train_time > ele%value(reset_value$)) then print *, "train_time > reset_value" print *, "This is very bad" end if ! Add the apropriate bunch offset (if you want to start at, say bunch 10) train_time = train_time + ecloud_data%bunchtimes(start_bunch) ! === FIND QDENS === extra_qdens = interpolateQdens(ecloud_data, train_time) * qdens_multiplier ! This should average stuff call interpolateField(ecloud_data, x, y, train_time, Ex, Ey) field%E(1) = (Ex * extra_qdens) field%E(2) = (Ey * extra_qdens) end subroutine !========================================================================! ! Calculates the E fields on the xy plane from the DiscreteEcloud struct ! !========================================================================! subroutine calculateEMFromDiscrete_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(DiscreteEcloud) :: 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 logical :: time_elapsed = .false. logical :: first_zero = .false. logical :: time_dc = .false. if (.not. on) then field%E(1) = 0 field%E(2) = 0 return end if ! Setting up ! vec=[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. ! Find time since head of train time_elapsed = (time - ele%value(latest_particle$)) .GE. ele%value(reset_value$) ! if enough time has elapsed since the last particle first_zero = ele%value(first_particle$) .LE. 0 ! or if the first particle was at time zero time_dc = ele%value(first_particle$) .GT. time ! or if the first particle is in the future if ( time_elapsed .or. first_zero .or. time_dc ) then ele%value(first_particle$) = time ! then consider this to be the first particle !print *, "[*] New start of train at time: ", time end if ! This is now the most recent particle ele%value(latest_particle$) = time ! Now, time since head of train is current time minus time of first particle train_time = time-ele%value(first_particle$) if (train_time > ele%value(reset_value$)) then print *, "train_time > reset_value" print *, "This is very bad" end if ! Add the apropriate bunch offset (if you want to start at, say bunch 10) train_time = train_time + ecloud_data%bunchtimes(start_bunch) ! === 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) ! 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! DEBUG ONLY! !write(1708,'(11es12.4)') time, train_time, bunch_index, & ! 1, 2, 3 ! orb%vec(1), orb%vec(3), orb%vec(5), & ! 4, 5, 6 ! interpolateQdens(ecloud_data, train_time), & ! 7 ! Ex, Ey, field%E(1), field%E(2) ! 8, 9, 10, 11 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 !====================================! ! 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) 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 Discrete struct ! !========================================================! ! ecloud_data - data struct initialized with calculateDiscrete*** call ! train_time - real time since beginning of train of particles ! [spline] - optional: logical whether to use splinging or linear interpolation, defaults true real(rp) function interpolateQdens(ecloud_data, train_time, usespline) type(DiscreteEcloud) :: ecloud_data real(rp) :: train_time, extra_qdens, d_ex_q integer :: i=0 logical :: spline_method logical, optional :: usespline logical :: ok if (present(usespline)) then spline_method = usespline else spline_method = .true. endif ! Extrapolate qdens extra_qdens = 0 if (spline_method) then !if (train_time > ecloud_data%qdens_spline(size(ecloud_data%qdens_spline))%x ) then ! extra_qdens = 0 !else call spline_evaluate(ecloud_data%qdens_spline, train_time, ok, extra_qdens, d_ex_q) !end if else if ( train_time < ecloud_data%qdenstimes(1) ) then ! Extend slope from first and second points extra_qdens = cy( ecloud_data%qdenstimes(1),ecloud_data%qdens(1),ecloud_data%qdenstimes(2),ecloud_data%qdens(2),train_time) else do i=1, size(ecloud_data%qdenstimes) if (i==size(ecloud_data%qdenstimes)) then ! Extend slope between last two points extra_qdens = cy( ecloud_data%qdenstimes(i-1),ecloud_data%qdens(i-1),ecloud_data%qdenstimes(i),ecloud_data%qdens(i),train_time) exit else if (train_time < ecloud_data%qdenstimes(i+1)) then ! Implicitely (AND 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 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 Discrete struct ! !========================================================! subroutine interpolateField(ecloud_data, x, y, bunch_time, Ex, Ey) use bmad use bmad_struct implicit none type(DiscreteEcloud), 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(DiscreteEcloud) :: ecloud_data integer :: lowb = 0 real(rp) :: exX, exVal real(rp) :: test_val 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 *, " ", "[*] SPECIFIED START BUNCH: ", start_bunch print *, "[+] FIRST TESTS: QDENS, LOWBOUND, CY" exX = 0.0 ! Add the apropriate bunch offset (if you want to start at, say bunch 10) test_val = exX + ecloud_data%bunchtimes(start_bunch) lowb = lowbound(ecloud_data%qdenstimes, test_val) !exVal = cy(ecloud_data%qdenstimes(lowb),ecloud_data%qdens(lowb),ecloud_data%qdenstimes(lowb+1),ecloud_data%qdens(lowb+1),exX) exVal = interpolateQdens(ecloud_data, test_val) 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),test_val exX = 1.54e-7 ! Right between 2 and 3 ! Add the apropriate bunch offset (if you want to start at, say bunch 10) test_val = exX + ecloud_data%bunchtimes(start_bunch) lowb = lowbound(ecloud_data%qdenstimes, test_val) !exVal = cy(ecloud_data%qdenstimes(lowb),ecloud_data%qdens(lowb),ecloud_data%qdenstimes(lowb+1),ecloud_data%qdens(lowb+1),exX) exVal = interpolateQdens(ecloud_data, test_val) 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),test_val 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(DiscreteEcloud) :: ecloud_data character(len=4) :: mode real(rp) :: valmin, valmax, valstep real(rp) :: x, y, t, Ex, Ey, qdens, temp_t 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 ! Add the apropriate bunch offset (if you want to start at, say bunch 10) t = t + ecloud_data%bunchtimes(start_bunch) x = valmin print *, "[*] x ", "Ex ", "Ey ", "Qdens" do while (x<=valmax) call interpolateField(ecloud_data, x, y, t, Ex, Ey) qdens = interpolateQdens(ecloud_data, t) print '(4es12.4)', x,Ex,Ey,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 ! Add the apropriate bunch offset (if you want to start at, say bunch 10) t = t + ecloud_data%bunchtimes(start_bunch) y = valmin print *, "[*] y ", "Ex ", "Ey ", "Qdens" do while (y<=valmax) call interpolateField(ecloud_data, x, y, t, Ex, Ey) qdens = interpolateQdens(ecloud_data, t) print '(4es12.4)', y,Ex,Ey,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 print *, "[*] t ", "Ex ", "Ey ", "Qdens" do while (t<=valmax) ! Add the apropriate bunch offset (if you want to start at, say bunch 10) temp_t = t + ecloud_data%bunchtimes(start_bunch) call interpolateField(ecloud_data, x, y, temp_t, Ex, Ey) qdens = interpolateQdens(ecloud_data, temp_t) print '(4es12.4)', t,Ex,Ey,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 calculateDiscreteFromBeamfieldFile(fileunit, ecloud_data) implicit none integer, intent(in) :: fileunit type (DiscreteEcloud) :: ecloud_data end subroutine !==========================================================! ! Creates an ecloud_data struct from a snapshots directory ! !==========================================================! subroutine calculateDiscreteFromSnapshotsFiles(snapfileunit, densfileunit, bunchunit, ecloud_data) implicit none integer, intent(in) :: snapfileunit, densfileunit, bunchunit integer :: RetCode character(len=500) :: line type (DiscreteEcloud) :: 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(:) type(spline_struct), allocatable :: splines(:) logical :: ok ! 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)) allocate(splines(qdens_num+1)) 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) splines(t)%x = times(t) splines(t)%y = qdensity(t) t = t+1 end do qdens_dens splines(t)%x = splines(t-1)%x * 2 splines(t)%y = 0 close(densfileunit) ! Make splines call spline_akima(splines, ok) !========================! ! 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) allocate(ecloud_data%qdens_spline(size(splines))) ecloud_data%qdens_spline = splines end subroutine subroutine calculateDiscreteFromQdensFile(densfileunit, ecloud_data) implicit none integer, intent(in) :: densfileunit type (DiscreteEcloud) :: ecloud_data integer :: qdens_num, RetCode, t character(len=500) :: line real(rp), allocatable :: times(:) real(rp), allocatable :: qdensity(:) type(spline_struct), allocatable :: splines(:) logical :: ok !============================! ! 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)) allocate(splines(qdens_num+1)) 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) splines(t)%x = times(t) splines(t)%y = qdensity(t) t = t+1 end do qdens_dens splines(t)%x = splines(t-1)%x * 2 splines(t)%y = 0 close(densfileunit) ! Make splines call spline_akima(splines, ok) call allocateCopy(ecloud_data%qdens, qdensity) call allocateCopy(ecloud_data%qdenstimes, times) allocate(ecloud_data%qdens_spline(size(splines))) ecloud_data%qdens_spline = splines 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