module field_emitter_mod use beam_mod use wall3d_mod use em_field_mod !use field_emitter_struct implicit none type field_emitter_struct real(rp) :: x, y, s, s_rel ! Position on wall real(rp) :: normal_x, normal_y, normal_s ! Wall normal vector real(rp) :: I_av ! average current real(rp) :: freq ! frequency real(rp) :: min_ratio real(rp) :: A_fn, beta_fn ! Fowler-Nordheim A and beta character(40) :: ele_name integer :: ix_ele type(coord_struct), allocatable :: particles(:) endtype contains !+ ! subroutine wall_point_finder(phi, s_rel, ele, wall_point, normals) ! ! Subroutine gets x, y, s_rel coordinates of location of the wall at a particular ! s_rel and angle phi and the normal vector to the wall at that point. ! ! Input: ! phi -- real(rp): specifies the angle at which the field emitter is located ! ! s_rel -- real(rp): specifies the s-coordinate of the field emitter relative ! to the upstream end of the element (ele). ! If ele%key == e_gun, then s_rel is the radius on the cathode ! ! ele -- type(ele_struct): specifies the element in which the field emitter ! is located ! ! Output: ! wall_point(3) -- real(rp): (/x, y, s_rel/) ! ! normals(3) -- real(rp): the normal vector at the point specified by ! wall_point !- subroutine wall_point_finder(phi, s_rel, ele, FE, wall_point, normals) type(field_emitter_struct) :: FE real(rp) :: d_radius,r_wall, r_particle, wall_point0(3), wall_point(3), normals(3), r0, phi real(rp) :: x, y, s_rel, r0_gun type(ele_struct), pointer :: ele character(40) :: spec_file if (ele%key == lcavity$ .or. ele%key == rfcavity$) then r0 = 0.001 x = r0*cos(phi) y = r0*sin(phi) d_radius = wall3d_d_radius([x, 0.0_rp, y, 0.0_rp, s_rel, 0.0_rp],& ele, 1, normals) r_wall = r0 - d_radius x = r_wall*cos(phi) y = r_wall*sin(phi) wall_point = [x, y, s_rel] elseif (ele%key == e_gun$) then ! s_rel is used as radius r0_gun = s_rel x = r0_gun*cos(phi) y = r0_gun*sin(phi) normals = [0, 0, 1] wall_point = [x, y, 1e-6_rp] else print *, 'Element must be either lcavity or e_gun.' stop endif end subroutine !+ ! function field_emitter_current(ele, t, FE) nonzero_q(I_fn) ! ! Routine to calculate the instantaneous current at time t of a field-emitter (I_fn) on the wall ! of an element (position of wall specified by x, y, and s_rel in FE). ! Uses the Fowler-Nordheim equation given on pg. 94 of RF Superconductivity, ! by HasanPadamsee to determine I_fn. ! ! Input: ! lat -- type (lat_struct): lat_struct of lattice which contains the element ! for which we are sampling field emitters. ! ! t -- real(rp): time at which I_fn is to be calculated. ! ! FE -- type(field_emitter_struct): contains the position x, y, s_rel of ! particle, its normal to the wall, and normal_x, normal_y, normal_s. ! ! Output: ! I_fn -- real(rp): instantanteous emission current at point x, y, s_rel, at time t. !- function field_emitter_current(ele, t, FE) result(I_fn) type(ele_struct), pointer :: ele type(field_emitter_struct) :: FE real(rp) ::lar, A_fn, beta_fn, t, q, I_fn, E_dot_n, total_q, I_temp type(coord_struct) :: particle type(em_field_struct) :: field particle%species = electron$ !default particle particle%vec(1) = FE%x !Need to define xyz in terms of vec(1:6) particle%vec(3) = FE%y !for valid argument type in em_field_calc particle%vec(5) = FE%s_rel call em_field_calc(ele, ele%branch%lat%param,& !FindE-field at given zxy, t, FE%s_rel, particle, .false., field, .false., rf_time = t) lar = sqrt( field%E(1)**2 + field%E(2)**2 + field%E(3)**2 ) E_dot_n = FE%normal_x*field%E(1) + FE%normal_y*field%E(2) + FE%normal_s*field%E(3) I_fn = 3.85e-7*FE%A_fn*( (FE%beta_fn*E_dot_n)**2)*exp(-5.464e10/(FE%beta_fn*E_dot_n) ) !print *, 'I_fn=', I_fn if (E_dot_n <= 0) then I_fn = 0 endif end function !+ ! subroutine field_emitter_init(ele, FE, t_a, t_b, n_t, min_ratio) ! ! Subroutine takes FE and populates the array (of type (coord_struct)) FE%particles; each item of which ! having the proper weighting of charge at a time t. This will be done n_t times ! at times t_a through t_b. ! ! Note: Particles that have charge weight zero are thrown out unless the ! optional argument min_ratio = 0. ! ! Input: ! ele -- type(ele_struct): the ele_struct of the element in which the field ! emitter will be weighted. ! ! FE -- type(field_emitter_struct): FE%x, FE%y, FE%s_rel specify the position ! of the field emitter. ! ! t_a -- real(rp): specifies the lower bound on the times at which the field ! emitter will be weighted. ! ! t_b -- real(rp): specifies the upper bound on the times at which the ! field emitter will be weighted ! ! n_t -- integer: specifies the number of time samples, or number of ! particles that will be weigted. ! ! min_ratio -- real(rp), optional: specifies the minimun ratio of the maximum ! charge to the charge of a given particle. If min_ratio = 0 then ! all particles are used to populate FE%particles. If min_ratio ! is absent then min_ratio is set to 1e9. ! ! Output: ! FE -- type(field_emitter_struct): the populated array FE%particles, each ! item of which being a particle with an associated charge at a given ! time. !- subroutine field_emitter_init(ele, FE, t_a, t_b, n_t, min_ratio) type(ele_struct), pointer :: ele real(rp) :: dt, max_q, min_r, q_scale_ratio, total_q, f, I_temp, t_a, t_b real(rp), optional :: min_ratio integer :: n_t, i, j, k, nonzero_q, part_size type(field_emitter_struct) :: FE type(coord_struct) :: particle type(coord_struct), allocatable :: temp_particles(:) !Particle at x, y, s_rel at dt = (t_b - t_a) / n_t !defining particle with coord struct for convenience particle%species = electron$ !particle is assumed to be electron particle%vec(1) = FE%x particle%vec(3) = FE%y particle%vec(5) = FE%s_rel allocate(temp_particles(n_t)) do i=1, n_t particle%charge = field_emitter_current(ele, dt*(i-1) + t_a, FE)*dt particle%t = dt*(i-1) + t_a temp_particles(i) = particle enddo total_q = sum(temp_particles(:)%charge) I_temp = total_q*FE%freq if (total_q <= 1e-30_rp) return q_scale_ratio = abs(( FE%I_av / I_temp )) ! Set A_FN FE%A_FN = FE%A_FN * q_scale_ratio do i=1, n_t temp_particles(i)%charge = (temp_particles(i)%charge)*q_scale_ratio enddo nonzero_q = count(temp_particles(:)%charge /= 0) allocate(FE%particles(nonzero_q)) k=0 do i=1, n_t !do loop first populates FE%particles with particles if (temp_particles(i)%charge /= 0) then !with nonzero charge k=k+1 FE%particles(k) = temp_particles(i) endif enddo if (present(min_ratio) .and. min_ratio == 0) then return else if (present(min_ratio) .and. min_ratio /= 0) then max_q = maxval(temp_particles(:)%charge) min_ratio = FE%min_ratio deallocate(FE%particles) allocate(FE%particles(count( max_q / min_ratio <=& temp_particles(:)%charge ))) k=0 do i=1,size(temp_particles) if ( max_q / temp_particles(i)%charge <= min_ratio) then k=k+1 FE%particles(k) = temp_particles(i) endif enddo return else max_q = maxval(temp_particles(:)%charge) min_r = 1e9 !min_ratio value if none is provided deallocate(FE%particles) allocate(FE%particles(count( max_q / min_r <=& temp_particles(:)%charge ))) k=0 do i=1,size(temp_particles) if ( max_q / temp_particles(i)%charge <= min_r) then k=k+1 FE%particles(k) = temp_particles(i) endif enddo return endif end subroutine function particle_counter(particles) result(n_particles) integer :: n_particles type(coord_struct) :: particles(:) n_particles = size(particles) end function !+ ! subroutine particle_writer(particles, unit_num) ! ! Subroutine writes the coord_struct array, particles, into a text file with ! unit number unit_num. ! ! Input: ! particles(:) -- type(coord_struct), allocatable: list of particles to be ! written. Format of a single line is (x, cpx, y, cpy, z, cpz, t, q) ! ! unit_num -- integer: specifies the unit number of the file to which the ! particles will be written. !- subroutine particle_writer(particles, unit_num) !particle is an element of particles(:) type(coord_struct) :: particles(:) integer :: unit_num, i if (ubound(particles, 1) /= 0) then do i=1, size(particles) write(unit_num, '(5ES25.10,$)') particles(i)%vec(1), particles(i)%vec(2), & particles(i)%vec(3), particles(i)%vec(4), particles(i)%vec(5), & particles(i)%vec(6), particles(i)%t, particles(i)%charge enddo else return endif end subroutine !+ ! subroutine write_fe_info(fe, iu) ! ! Write info about a field emitter to a file subroutine write_fe_info(fe, iu) type(field_emitter_struct) :: fe integer :: iu character(40) :: fmt = '(a,x,a,es15.8)' character(3) :: delim = ' = ' ! write(iu, '(3a)') 'ele_name', delim, trim(fe%ele_name) write(iu, '(a,x,a,x,i0)') 'ix_ele', delim, fe%ix_ele write(iu, fmt) 'x', delim, fe%x write(iu, fmt) 'y', delim, fe%y write(iu, fmt) 's', delim, fe%s write(iu, fmt) 's_rel', delim, fe%s_rel write(iu, fmt) 'I_av', delim, fe%I_av write(iu, fmt) 'freq', delim, fe%freq write(iu, fmt) 'A_FN', delim, fe%A_FN write(iu, fmt) 'beta_FN', delim, fe%beta_FN write(iu, '(a,x,a,x,i0)') 'n_particles',delim, size(fe%particles) end subroutine !+ ! subroutine read_FE(unit_spec, lat, FE, ix_FE) ! ! Subroutine reads in the .dat text file (with unit number unit_spec), and ! populates the coord_struct array FE%particles such that each of its elements ! is a particle of properly weighted charge at a corresponding time t. ! ! Note: Text file format is the following. ! ! line 1 = n (integer numer of field emitters) ! line 2 = s_rel angle I_av beta_fn n_t f ele_name ! ... ! ... ! line n+1 = s_rel angle I_av beta_fn n_t f ele_name ! ! Input: ! unit_spec -- integer: specifies the unit number of the .dat text file that contains ! the specifications of the field emitters. ! ! lat -- type(lat_struct): speicifies the lattice that contains the ! elements whose names are in the .dat input text file. ! ! FE -- type(field_emitter_struct): the struct to be populated ! ! ! Output: ! FE -- type(field_emitter_struct): contains the fully populated ! coord_struct array FE%particles. !- subroutine read_FE(unit_spec, lat, FE) type(lat_struct) :: lat type(field_emitter_struct) :: FE type(ele_struct), pointer :: ele type(ele_pointer_struct), allocatable :: eles(:) real(rp) :: s, min_ratio, t_a, t_b, phi, s_rel, wall_point(3), normals(3), angle, offset integer ::j,ix, i, unit_spec, n_t, num_loc, io logical :: err character(40) :: spec_file FE%A_fn = 1e-9 !dummy A_fn value offset = 1e-6 ! offset of field emitter from the wall io = 0 read(unit_spec, *, iostat = io) s_rel, phi, FE%I_av,& FE%beta_fn, n_t, FE%freq, FE%ele_name if (io /= 0 ) then print *, 'Input text file error!' print *, 'Text file must be of the form' print *, 'line 1 = n (integer number of field emitters)' print *, 'line 2 = s_rel angle I_av beta_fn n_t f ele_name' print *, '...' print *, '...' print *, 'line n+1 = s_rel angle I_av beta_fn n_t f ele_name' endif print *, 'Locating ', trim(adjustl(FE%ele_name)), '..' call lat_ele_locator(FE%ele_name, lat, eles, num_loc, err) If(num_loc == 0) then print *, 'Number of located elements is zero' print *, 'Stopping..' stop else print *, 'Found: ', trim(adjustl(FE%ele_name)) endif ele => eles(1)%ele FE%ix_ele = ele%ix_ele if(FE%freq == 0) then print *, 'Frequency input is 0' print *, 'Finding frequency of ', trim(adjustl(FE%ele_name)) if (ele%key == lcavity$ .or. ele%key == e_gun$) then FE%freq = ele%value(rf_frequency$) if (FE%freq == 0) then print *, 'Element must have nonzero rf_frequency' stop endif else print *, 'Element must be an lcavity if no frequency is defined.' stop endif print *,'Frequency of ', trim(adjustl(FE%ele_name)), ' found: ', FE%freq endif print *, 'Finding Cartesian postion of field emitter..' call wall_point_finder(phi, s_rel, ele, FE, wall_point, normals) FE%x = wall_point(1) - offset*normals(1) FE%y = wall_point(2) - offset*normals(2) FE%s_rel = wall_point(3)- offset*normals(3) FE%s = FE%s_rel + ele%s - ele%value(L$) FE%normal_x = normals(1) FE%normal_y = normals(2) FE%normal_s = normals(3) print *, 'Found position of field field emitter.' t_a = 0 t_b = 1 / FE%freq FE%min_ratio = 1e9 print *, 'Sampling field emitter over one period..' call field_emitter_init(ele, FE, t_a, t_b, n_t) !changing to global s for dark current tracker print *, 'Obtained data for ', size(FE%particles), ' particles.' FE%particles(:)%vec(5) = ele%s - ele%value(l$) + FE%s_rel end subroutine !+ ! subroutine wall_contour(ele, angle, file_num) ! ! Subroutine calculates the maximum E-field along the contour of the wall traced ! out along the angle given in the argument of the subroutine. ! ! Input: ! ele -- type(ele_struct): element for which the max E-field is to be ! calculated ! angle -- real(rp): angle that specifies the contour traced out along the ! wall ! file_num -- integer: the unit file number to which the results will be ! written ! Output: ! 'contour.ex' -- file containing the calculated data, written such that a ! line is [x, y, s_rel, E_max] !- subroutine wall_contour(ele, angle, file_num, n_contour_points, s_init) real(rp), optional :: s_init integer :: j, n, file_num, n_contour_points type(ele_struct) :: ele real(rp) :: length, L, length_upto, ds, angle, dr_dtheta, rad(2), E_max real(rp) :: point(3), d_radius, perp(3), delta, total_len, r0, r1, s0, s1 type(coord_struct) :: c_point type(em_field_struct) :: field0, field1 !get the total contour length of the element ! !Note: contour length is not ele%value(ele_num$) call calc_wall_radius(ele%wall3d(1)%section(1)%v, cos(angle), sin(angle), r0, dr_dtheta) s0 = ele%wall3d(1)%section(1)%s length = 0 print *, 'Finding contour length of ', trim(adjustl(ele%name)), '..' do j=2, ubound(ele%wall3d(1)%section, 1) call calc_wall_radius(ele%wall3d(1)%section(j)%v, cos(angle), sin(angle), r1, dr_dtheta) s1 = ele%wall3d(1)%section(j)%s length = length + hypot( (r1 - r0), (s1 - s0) ) s0 = s1 r0 = r1 enddo total_len = length print *, "Total contour length of ", trim(adjustl(ele%name)), ":", length,"m" !n_contour_points is the number of sample points at which E_max will be calculated ds = total_len / n_contour_points call calc_wall_radius(ele%wall3d(1)%section(1)%v, cos(angle), sin(angle),r0, dr_dtheta) s0 = ele%wall3d(1)%section(1)%s length = 0 c_point%vec(2) = 0 c_point%vec(4) = 0 c_point%vec(6) = 0 n=1 do j=2, ubound(ele%wall3d(1)%section, 1) call calc_wall_radius(ele%wall3d(1)%section(j)%v, cos(angle), sin(angle), r1, dr_dtheta) s1 = ele%wall3d(1)%section(j)%s length_upto = length length = length + hypot( ( r1 - r0), (s1 - s0) ) do if (n*ds < length ) then L = n*ds - length_upto point = [r0*cos(angle), r0*sin(angle), s0]& +[(r1-r0)*cos(angle), (r1-r0)*sin(angle), & (s1 - s0) ]*( L / (length - length_upto)) c_point%species = electron$ c_point%vec(1) = point(1) c_point%vec(3) = point(2) c_point%vec(5) = point(3) d_radius = wall3d_d_radius(c_point%vec, ele, 1, perp) delta = 1e-6 point = point - delta*perp call em_field_calc(ele, ele%branch%lat%param, point(3), & c_point, .false., field0, .false., rf_time = 0.0_rp) call em_field_calc(ele, ele%branch%lat%param, point(3), & c_point, .false., field1, .false., rf_time = 1 / (4*ele%value(rf_frequency$))) E_max = sqrt(field0%E(1)**2 + field0%E(2)**2 + field0%E(3)**2& + field1%E(1)**2 + field1%E(2)**2 + field1%E(3)**2) point(3) = point(3) + s_init write(file_num,'(4ES25.10, A25)') point(1), point(2), point(3), E_max,& trim(adjustl(ele%name)) n=n+1 else exit endif enddo r0 = r1 s0 = s1 enddo endsubroutine end module