diff --git a/examples/gk-ions/2D-periodic-gk.toml b/examples/gk-ions/2D-periodic-gk.toml
new file mode 100644
index 000000000..df3f04e7d
--- /dev/null
+++ b/examples/gk-ions/2D-periodic-gk.toml
@@ -0,0 +1,87 @@
+n_ion_species = 1
+n_neutral_species = 0
+boltzmann_electron_response = true
+evolve_moments_density = false
+evolve_moments_parallel_flow = false
+evolve_moments_parallel_pressure = false
+evolve_moments_conservation = false
+gyrokinetic_ions = true
+T_e = 1.0
+T_wall = 1.0
+initial_density1 = 1.0
+initial_temperature1 = 1.0
+z_IC_option1 = "gaussian"
+z_IC_density_amplitude1 = 0.001
+z_IC_density_phase1 = 0.0
+z_IC_upar_amplitude1 = 1.0
+z_IC_upar_phase1 = 0.0
+z_IC_temperature_amplitude1 = 0.0
+z_IC_temperature_phase1 = 0.0
+vpa_IC_option1 = "gaussian"
+vpa_IC_density_amplitude1 = 1.0
+vpa_IC_density_phase1 = 0.0
+vpa_IC_upar_amplitude1 = 0.0
+vpa_IC_upar_phase1 = 0.0
+vpa_IC_temperature_amplitude1 = 0.0
+vpa_IC_temperature_phase1 = 0.0
+initial_density2 = 1.0
+initial_temperature2 = 1.0
+z_IC_option2 = "gaussian"
+z_IC_density_amplitude2 = 0.001
+z_IC_density_phase2 = 0.0
+z_IC_upar_amplitude2 = -1.0
+z_IC_upar_phase2 = 0.0
+z_IC_temperature_amplitude2 = 0.0
+z_IC_temperature_phase2 = 0.0
+vpa_IC_option2 = "gaussian"
+vpa_IC_density_amplitude2 = 1.0
+vpa_IC_density_phase2 = 0.0
+vpa_IC_upar_amplitude2 = 0.0
+vpa_IC_upar_phase2 = 0.0
+vpa_IC_temperature_amplitude2 = 0.0
+vpa_IC_temperature_phase2 = 0.0
+charge_exchange_frequency = 0.5
+ionization_frequency = 0.05
+constant_ionization_rate = true
+nstep = 50
+dt = 1.0e-3
+nwrite = 10
+nwrite_dfns = 10
+use_semi_lagrange = false
+n_rk_stages = 4
+split_operators = false
+r_ngrid = 5
+r_nelement = 2
+r_bc = "periodic"
+z_ngrid = 5
+z_nelement = 2
+z_bc = "periodic"
+z_discretization = "chebyshev_pseudospectral"
+vpa_ngrid = 5
+vpa_nelement = 2
+vpa_L = 6.0
+vpa_bc = "zero"
+vpa_discretization = "chebyshev_pseudospectral"
+vperp_ngrid = 5
+vperp_nelement = 1
+vperp_L = 3.0
+vperp_bc = "zero"
+vperp_discretization = "chebyshev_pseudospectral"
+vz_ngrid = 9
+vz_nelement = 64
+vz_L = 18.0
+vz_bc = "both_zero"
+vz_discretization = "chebyshev_pseudospectral"
+
+[numerical_dissipation]
+vpa_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+vperp_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+#r_disspipation_coefficient = 1.0e-3
+#force_minimum_pdf_value = 0.0
+
+[geometry]
+#option="1D-mirror"
+DeltaB=0.0
+option="constant-helical"
+pitch=0.1
+rhostar= 0.1
diff --git a/makie_post_processing/makie_post_processing/src/shared_utils.jl b/makie_post_processing/makie_post_processing/src/shared_utils.jl
index 4a7b32b7b..031bcda75 100644
--- a/makie_post_processing/makie_post_processing/src/shared_utils.jl
+++ b/makie_post_processing/makie_post_processing/src/shared_utils.jl
@@ -88,6 +88,7 @@ function get_composition(scan_input)
     phi_wall = get(scan_input, "phi_wall", 0.0)
     # if false use true Knudsen cosine for neutral wall bc
     use_test_neutral_wall_pdf = get(scan_input, "use_test_neutral_wall_pdf", false)
+    gyrokinetic_ions = get(scan_input, "gyrokinetic_ions", false)
     # constant to be used to test nonzero Er in wall boundary condition
     Er_constant = get(scan_input, "Er_constant", 0.0)
     recycling_fraction = get(scan_input, "recycling_fraction", 1.0)
@@ -106,7 +107,7 @@ function get_composition(scan_input)
     me_over_mi = 1.0/1836.0
     composition = species_composition(n_species, n_ion_species, n_neutral_species,
         electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant,
-        mn_over_mi, me_over_mi, recycling_fraction, allocate_float(n_species))
+        mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species))
     return composition
 
 end
diff --git a/moment_kinetics/debug_test/gyroaverage_inputs.jl b/moment_kinetics/debug_test/gyroaverage_inputs.jl
new file mode 100644
index 000000000..d9564238e
--- /dev/null
+++ b/moment_kinetics/debug_test/gyroaverage_inputs.jl
@@ -0,0 +1,87 @@
+test_type = "gyroaverage"
+
+# default inputs for tests
+test_input = Dict(
+    "run_name" => "gyroaverage",
+    "n_ion_species" => 1,
+    "n_neutral_species" => 0,
+    "evolve_moments_density" => false,
+    "evolve_moments_parallel_flow" => false,
+    "evolve_moments_parallel_pressure" => false,
+    "evolve_moments_conservation" => false,
+    "gyrokinetic_ions" => true,
+    "T_e" => 1.0,
+    "T_wall" => 1.0,
+    "initial_density1" => 1.0,
+    "initial_temperature1" => 1.0,
+    "z_IC_option1" => "gaussian",
+    "z_IC_density_amplitude1" => 0.001,
+    "z_IC_density_phase1" => 0.0,
+    "z_IC_upar_amplitude1" => 1.0,
+    "z_IC_upar_phase1" => 0.0,
+    "z_IC_temperature_amplitude1" => 0.0,
+    "z_IC_temperature_phase1" => 0.0,
+    "vpa_IC_option1" => "gaussian",
+    "vpa_IC_density_amplitude1" => 1.0,
+    "vpa_IC_density_phase1" => 0.0,
+    "vpa_IC_upar_amplitude1" => 0.0,
+    "vpa_IC_upar_phase1" => 0.0,
+    "vpa_IC_temperature_amplitude1" => 0.0,
+    "vpa_IC_temperature_phase1" => 0.0,
+    "initial_density2" => 1.0,
+    "initial_temperature2" => 1.0,
+    "z_IC_option2" => "gaussian",
+    "z_IC_density_amplitude2" => 0.001,
+    "z_IC_density_phase2" => 0.0,
+    "z_IC_upar_amplitude2" => -1.0,
+    "z_IC_upar_phase2" => 0.0,
+    "z_IC_temperature_amplitude2" => 0.0,
+    "z_IC_temperature_phase2" => 0.0,
+    "vpa_IC_option2" => "gaussian",
+    "vpa_IC_density_amplitude2" => 1.0,
+    "vpa_IC_density_phase2" => 0.0,
+    "vpa_IC_upar_amplitude2" => 0.0,
+    "vpa_IC_upar_phase2" => 0.0,
+    "vpa_IC_temperature_amplitude2" => 0.0,
+    "vpa_IC_temperature_phase2" => 0.0,
+    "charge_exchange_frequency" => 0.5,
+    "ionization_frequency" => 0.05,
+    "constant_ionization_rate" => true,
+    "nstep" => 3,
+    "dt" => 1.0e-12,
+    "nwrite" => 2,
+    "nwrite_dfns" => 2,
+    "n_rk_stages" => 4,
+    "r_ngrid" => 5,
+    "r_nelement" => 2,
+    "r_bc" => "periodic",
+    "z_ngrid" => 5,
+    "z_nelement" => 2,
+    "z_bc" => "periodic",
+    "z_discretization" => "chebyshev_pseudospectral",
+    "vpa_ngrid" => 5,
+    "vpa_nelement" => 2,
+    "vpa_L" => 6.0,
+    "vpa_bc" => "zero",
+    "vpa_discretization" => "chebyshev_pseudospectral",
+    "vperp_ngrid" => 5,
+    "vperp_nelement" => 1,
+    "vperp_L" => 3.0,
+    "vperp_bc" => "zero",
+    "vperp_discretization" => "chebyshev_pseudospectral",
+    "vz_ngrid" => 5,
+    "vz_nelement" => 2,
+    "vz_L" => 6.0,
+    "vz_bc" => "zero",
+    "vz_discretization" => "chebyshev_pseudospectral",
+    "numerical_dissipation" => Dict("vpa_dissipation_coefficient" => 1.0e-3,
+                                    "vperp_dissipation_coefficient" => 1.0e-3),
+    "geometry" => Dict("DeltaB"=>0.0,
+                       "option"=>"constant-helical",
+                       "pitch"=>0.1,
+                       "rhostar"=> 0.1),
+)
+
+test_input_list = [
+     test_input,
+    ]
diff --git a/moment_kinetics/debug_test/gyroaverage_tests.jl b/moment_kinetics/debug_test/gyroaverage_tests.jl
new file mode 100644
index 000000000..80d84ed17
--- /dev/null
+++ b/moment_kinetics/debug_test/gyroaverage_tests.jl
@@ -0,0 +1,23 @@
+module GyroaverageDebug
+
+# Debug test using wall boundary conditions.
+
+include("setup.jl")
+
+# Create a temporary directory for test output
+test_output_directory = get_MPI_tempdir()
+mkpath(test_output_directory)
+
+
+# Input parameters for the test
+include("gyroaverage_inputs.jl")
+
+# Defines the test functions, using variables defined in the *_inputs.jl file
+include("runtest_template.jl")
+
+end # GyroaverageDebug
+
+
+using .GyroaverageDebug
+
+GyroaverageDebug.runtests()
diff --git a/moment_kinetics/debug_test/runtests.jl b/moment_kinetics/debug_test/runtests.jl
index 55550595e..6c2d69d39 100644
--- a/moment_kinetics/debug_test/runtests.jl
+++ b/moment_kinetics/debug_test/runtests.jl
@@ -10,6 +10,7 @@ function runtests()
         include(joinpath(@__DIR__, "harrisonthompson.jl"))
         include(joinpath(@__DIR__, "mms_tests.jl"))
         include(joinpath(@__DIR__, "restart_interpolation_tests.jl"))
+        include(joinpath(@__DIR__, "gyroaverage_tests.jl"))
     end
 end
 
diff --git a/moment_kinetics/src/chebyshev.jl b/moment_kinetics/src/chebyshev.jl
index f63f22603..e3c789aec 100644
--- a/moment_kinetics/src/chebyshev.jl
+++ b/moment_kinetics/src/chebyshev.jl
@@ -667,7 +667,7 @@ function chebyshev_radau_weights(moments::Array{mk_float,1}, n)
     # create array for moments on extended [0,2π] domain in theta = ArcCos[z]
     fext = allocate_complex(nfft)
     # make fft plan
-    forward_transform = plan_fft!(fext, flags=FFTW.WISDOM_ONLY)
+    forward_transform = plan_fft!(fext, flags=FFTW.ESTIMATE)
     # assign values of fext from moments 
     @inbounds begin
         for j ∈ 1:n
diff --git a/moment_kinetics/src/coordinates.jl b/moment_kinetics/src/coordinates.jl
index 95aeef094..76d998ae6 100644
--- a/moment_kinetics/src/coordinates.jl
+++ b/moment_kinetics/src/coordinates.jl
@@ -98,6 +98,8 @@ struct coordinate
     element_shift::Array{mk_float,1}
     # option used to set up element spacing
     element_spacing_option::String
+    # list of element boundaries
+    element_boundaries::Array{mk_float,1}
 end
 
 """
@@ -167,7 +169,7 @@ function define_coordinate(input, parallel_io::Bool=false; run_directory=nothing
         cell_width, igrid, ielement, imin, imax, igrid_full, input.discretization, input.fd_option, input.cheb_option,
         input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch),
         scratch_2d, copy(scratch_2d), advection, send_buffer, receive_buffer, input.comm,
-        local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option)
+        local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option, element_boundaries)
 
     if coord.n == 1 && occursin("v", coord.name)
         spectral = null_velocity_dimension_info()
diff --git a/moment_kinetics/src/em_fields.jl b/moment_kinetics/src/em_fields.jl
index ebbdba5eb..6591e9062 100644
--- a/moment_kinetics/src/em_fields.jl
+++ b/moment_kinetics/src/em_fields.jl
@@ -14,21 +14,24 @@ using ..moment_kinetics_structs: em_fields_struct
 using ..velocity_moments: update_density!
 #using ..calculus: derivative!
 using ..derivatives: derivative_r!, derivative_z!
-
+using ..gyroaverages: gyro_operators, gyroaverage_field!
 """
 """
-function setup_em_fields(nz, nr, force_phi, drive_amplitude, drive_frequency, force_Er_zero)
+function setup_em_fields(nvperp, nz, nr, n_ion_species, force_phi, drive_amplitude, drive_frequency, force_Er_zero)
     phi = allocate_shared_float(nz,nr)
     phi0 = allocate_shared_float(nz,nr)
     Er = allocate_shared_float(nz,nr)
     Ez = allocate_shared_float(nz,nr)
-    return em_fields_struct(phi, phi0, Er, Ez, force_phi, drive_amplitude, drive_frequency, force_Er_zero)
+    gphi = allocate_shared_float(nvperp,nz,nr,n_ion_species)
+    gEr = allocate_shared_float(nvperp,nz,nr,n_ion_species)
+    gEz = allocate_shared_float(nvperp,nz,nr,n_ion_species)
+    return em_fields_struct(phi, phi0, Er, Ez, gphi, gEr, gEz, force_phi, drive_amplitude, drive_frequency, force_Er_zero)
 end
 
 """
 update_phi updates the electrostatic potential, phi
 """
-function update_phi!(fields, fvec, z, r, composition, z_spectral, r_spectral, scratch_dummy)
+function update_phi!(fields, fvec, vperp, z, r, composition, z_spectral, r_spectral, scratch_dummy, gyroavs::gyro_operators)
     n_ion_species = composition.n_ion_species
     @boundscheck size(fields.phi,1) == z.n || throw(BoundsError(fields.phi))
     @boundscheck size(fields.phi,2) == r.n || throw(BoundsError(fields.phi))
@@ -137,6 +140,20 @@ function update_phi!(fields, fvec, z, r, composition, z_spectral, r_spectral, sc
             fields.Ez[:,:] .= 0.0
         end
     end
+    
+    # get gyroaveraged field arrays for distribution function advance
+    gkions = composition.gyrokinetic_ions
+    if gkions
+        gyroaverage_field!(fields.gphi,fields.phi,gyroavs,vperp,z,r,composition)
+        gyroaverage_field!(fields.gEz,fields.Ez,gyroavs,vperp,z,r,composition)
+        gyroaverage_field!(fields.gEr,fields.Er,gyroavs,vperp,z,r,composition)
+    else # use the drift-kinetic form of the fields in the kinetic equation
+        @loop_s_r_z_vperp is ir iz ivperp begin
+            fields.gphi[ivperp,iz,ir,is] = fields.phi[iz,ir]
+            fields.gEz[ivperp,iz,ir,is] = fields.Ez[iz,ir]
+            fields.gEr[ivperp,iz,ir,is] = fields.Er[iz,ir]
+        end
+    end
 
 end
 
diff --git a/moment_kinetics/src/gyroaverages.jl b/moment_kinetics/src/gyroaverages.jl
new file mode 100644
index 000000000..f76edc742
--- /dev/null
+++ b/moment_kinetics/src/gyroaverages.jl
@@ -0,0 +1,353 @@
+"""
+module for supporting gyroaverages at 
+fixed guiding centre R and fixed position r
+"""
+module gyroaverages
+
+export gyro_operators
+export init_gyro_operators
+export gyroaverage_field!
+export gyroaverage_pdf!
+
+using ..type_definitions: mk_float, mk_int
+using ..array_allocation: allocate_float, allocate_shared_float
+using ..array_allocation: allocate_int, allocate_shared_int
+using ..looping
+using ..communication: MPISharedArray, comm_block, _block_synchronize
+
+struct gyro_operators
+    # matrix for applying a gyroaverage to a function F(r,vpa,vperp) at fixed r, with R = r - rhovec and rhovec = b x v / Omega
+    # the matrix can also be used ring-averaging a function F(R,vpa,vperp) at fixed R, since F(R,vpa,vperp) is independent of gyrophase
+    # other matrices are required for gyroaveraging functions that depend on gyrophase.
+    gyromatrix::MPISharedArray{mk_float,6}
+    gyroloopsizes::MPISharedArray{mk_int,4}
+    izpgyroindex::MPISharedArray{mk_int,5}
+    irpgyroindex::MPISharedArray{mk_int,5}
+end
+
+"""
+Function to initialise the gyroaverage matrix.
+Gyroaverages are assumed to have 
+no contribution from outside of the domain for
+non-periodic boundary conditions. Periodic boundary
+conditions are supported by taking data from the 
+the appropriate part of the domain, determined by
+the path of the particle into a periodic copy of the domain. 
+This behaviour should be modified if we aim to support 
+other nonzero boundary conditions for the z and r domains.
+"""
+
+function init_gyro_operators(vperp,z,r,gyrophase,geometry,composition;print_info=false)
+    gkions = composition.gyrokinetic_ions
+    if !gkions
+        gyromatrix =  allocate_shared_float(1,1,1,1,1,1)
+        gyroloopsizes = allocate_shared_int(1,1,1,1)
+        izpgyroindex = allocate_shared_int(1,1,1,1,1)
+        irpgyroindex = allocate_shared_int(1,1,1,1,1)
+    else
+       if print_info
+           println("Begin: init_gyro_operators")
+       end
+       gyromatrix = allocate_shared_float(z.n,r.n,vperp.n,z.n,r.n,composition.n_ion_species)
+       # an array to store the value for the number of points in the z', r' sum for each gyroaveraged field index 
+       gyroloopsizes = allocate_shared_int(vperp.n,z.n,r.n,composition.n_ion_species)
+       
+       # init the matrix!
+       # the first two indices are to be summed over
+       # the other indices are the "field" positions of the resulting gyroaveraged quantity
+       begin_serial_region()
+       @serial_region begin
+           zlist = allocate_float(gyrophase.n)
+           rlist = allocate_float(gyrophase.n)
+           zelementlist = allocate_int(gyrophase.n)
+           relementlist = allocate_int(gyrophase.n)
+           
+           @loop_s_r_z_vperp is ir iz ivperp begin
+               #println("ivperp, iz, ir: ",ivperp," ",iz," ",ir)
+               r_val = r.grid[ir]
+               z_val = z.grid[iz]
+               vperp_val = vperp.grid[ivperp]
+               rhostar = geometry.rhostar
+               # Bmag at the centre of the gyroaveraged path
+               Bmag = geometry.Bmag[iz,ir]
+               # bzeta at the centre of the gyroaveraged path
+               bzeta = geometry.bzeta[iz,ir]
+               # rho at the centre of the gyroaveraged path 
+               # (modify to include different mass or reference temperatures)
+               rho_val = vperp_val*rhostar/Bmag
+               # convert these angles to a list of z'(gphase) r'(gphase)
+               zrcoordinatelist!(gyrophase,zlist,rlist,rho_val,r_val,z_val,bzeta)
+               # determine which elements contain these z', r'
+               elementlist!(zelementlist,zlist,z)
+               elementlist!(relementlist,rlist,r)
+               #println(z_val,zlist)
+               #println(r_val,rlist)
+               #println(zelementlist)
+               #println(relementlist)
+               # initialise matrix to zero
+               @. gyromatrix[:,:,ivperp,iz,ir,is] = 0.0
+               for igyro in 1:gyrophase.n
+                   # integration weight from gyroaverage (1/2pi)* int d gyrophase
+                   gyrowgt = gyrophase.wgts[igyro]/(2.0*pi)
+                   # get information about contributing element
+                   
+                   izel = zelementlist[igyro]
+                   if izel < 1
+                       # z' point is outside of the grid, skip this point
+                       # if simply ignore contributions from outside of the domain
+                       #continue
+                       # if set to zero any <field> where the path exits the domain
+                       @. gyromatrix[:,:,ivperp,iz,ir,is] = 0.0
+                       if z.bc == "periodic"
+                           print("ERROR: -1 in zelementlist")
+                       end
+                       break
+                   end
+                   izmin, izmax = z.igrid_full[1,izel], z.igrid_full[z.ngrid,izel]
+                   znodes = z.grid[izmin:izmax]
+                   
+                   irel = relementlist[igyro]
+                   if irel < 1
+                       # r' point is outside of the grid, skip this point
+                       # if simply ignore contributions from outside of the domain
+                       #continue
+                       # if set to zero any <field> where the path exits the domain
+                       @. gyromatrix[:,:,ivperp,iz,ir,is] = 0.0
+                       if r.bc == "periodic"
+                           print("ERROR: -1 in relementlist")
+                       end
+                       break
+                   end
+                   irmin, irmax = r.igrid_full[1,irel], r.igrid_full[r.ngrid,irel]
+                   rnodes = r.grid[irmin:irmax]
+                   
+                   #println("igyro ",igyro)
+                   #println("izel ",izel," znodes ",znodes)
+                   #println("irel ",irel," rnodes ",rnodes)
+                   # sum over all contributing Lagrange polynomials from each
+                   # collocation point in the element
+                   icounter = 0
+                   for irgrid in 1:r.ngrid
+                       irp = r.igrid_full[irgrid,irel]
+                       rpoly = lagrange_poly(irgrid,rnodes,rlist[igyro])
+                       for izgrid in 1:z.ngrid
+                           izp = z.igrid_full[izgrid,izel]
+                           zpoly = lagrange_poly(izgrid,znodes,zlist[igyro])
+                           # add the contribution from this z',r'
+                           gyromatrix[izp,irp,ivperp,iz,ir,is] += gyrowgt*rpoly*zpoly
+                           icounter +=1
+                       end
+                   end
+                   #println("counter: ",icounter)
+               end
+              # count the number of nonzero (izp,irp) elements in gyromatrix for this ivperp, iz, ir, is
+              zero = 1.0e-14 
+              nsum = 0 
+              for irp in 1:r.n 
+                  for izp in 1:z.n
+                      if abs(gyromatrix[izp,irp,ivperp,iz,ir,is]) > zero
+                          nsum += 1
+                      end
+                  end                  
+              end
+              # assign this value to the gyroloopsizes array
+              gyroloopsizes[ivperp,iz,ir,is] = nsum
+           end
+        end
+           
+        # Broadcast the values in gyroloopsizes across the shared-memory block
+        _block_synchronize()
+        # initialise the arrays containing the indexing information
+        # use the fact that the first index cannot be larger than the size of z.n*r.n
+        # and accept that we are storing undefined values in exchange for storing the useful
+        # data in shared-memory.
+        izpgyroindex = allocate_shared_int(z.n*r.n,vperp.n,z.n,r.n,composition.n_ion_species)
+        irpgyroindex = allocate_shared_int(z.n*r.n,vperp.n,z.n,r.n,composition.n_ion_species)
+
+        # compute the indices on the root process  
+        @serial_region begin
+            zero = 1.0e-14  
+            @loop_s_r_z_vperp is ir iz ivperp begin
+                # fill these arrays with the index locations using the same
+                # conditions as used to create the gyroloopsizes array
+                # note that values of the array only up to nsum = gyroloopsizes[ivperp,iz,ir,is] will be filled
+                # any access to unassigned values of izpgyroindex or irpgyroindex will result in undefined behaviour
+                isum = 0 
+                for irp in 1:r.n 
+                    for izp in 1:z.n
+                        if abs(gyromatrix[izp,irp,ivperp,iz,ir,is]) > zero
+                            isum += 1
+                            izpgyroindex[isum,ivperp,iz,ir,is] = izp
+                            irpgyroindex[isum,ivperp,iz,ir,is] = irp
+                        end
+                    end                  
+                end
+            end
+        end
+        _block_synchronize()
+        if print_info
+            println("Finished: init_gyro_operators")
+        end
+    end
+    
+    gyro = gyro_operators(gyromatrix,gyroloopsizes,izpgyroindex,irpgyroindex)
+    return gyro
+end
+
+function zrcoordinatelist!(gyrophase,zlist,rlist,rho_val,r_val,z_val,bzeta)
+    ngyro = gyrophase.n
+    for i in 1:ngyro
+        gphase = gyrophase.grid[i]
+        zlist[i] = z_val + rho_val*cos(gphase)*bzeta
+        rlist[i] = r_val + rho_val*sin(gphase)
+    end
+end
+
+"""
+for a given list of coordinate values, determine in which elements they are found
+-1 indicates that the required element would be outside of the existing grid
+
+-- assume here that the coordinates are fully local in memory
+
+"""
+function elementlist!(elist,coordlist,coord)
+    zero = 1.0e-14
+    ngyro = size(coordlist,1)
+    nelement = coord.nelement_global
+    xebs = coord.element_boundaries
+    bc = coord.bc
+    L = coord.L
+    for i in 1:ngyro
+        if bc=="periodic"
+            x = coordlist[i]
+            # if x is outside the domain, shift x to the appropriate periodic copy within the domain
+            x0 = xebs[1] # the lower endpoint
+            y = x-x0
+            r = rem(y,L) + 0.5*L*(1.0 - sign(y)) # get the remainder of x - x0 w.r.t. the domain length L, noting that r should be r > 0
+            x = r + x0 # shift so that x is bounded below by x0
+            coordlist[i] = x # update coordlist for later use
+        end
+        # determine which element contains the position x
+        x = coordlist[i]
+        elist[i] = -1
+        for j in 1:nelement
+            # check internal locations
+            if (x - xebs[j])*(xebs[j+1] - x) > zero
+                elist[i] = j
+                break
+            # check element boundary 
+            # (lower or upper, force a choice of element for boundary values)
+            elseif (abs(x-xebs[j]) < 100*zero) || (abs(x-xebs[j+1]) < 100*zero && j == nelement)
+                elist[i] = j
+                break
+            end            
+        end
+    end
+    return 
+end
+
+"""
+Copy of function in fokker_planck_calculus.jl
+Lagrange polynomial
+args: 
+j - index of l_j from list of nodes
+x_nodes - array of x node values
+x - point where interpolated value is returned
+"""
+function lagrange_poly(j,x_nodes,x)
+    # get number of nodes
+    n = size(x_nodes,1)
+    # location where l(x0) = 1
+    x0 = x_nodes[j]
+    # evaluate polynomial
+    poly = 1.0
+    for i in 1:j-1
+            poly *= (x - x_nodes[i])/(x0 - x_nodes[i])
+    end
+    for i in j+1:n
+            poly *= (x - x_nodes[i])/(x0 - x_nodes[i])
+    end
+    return poly
+end
+
+"""
+function for gyroaveraging a field of shape (z,r)
+and filling the result into an array of shape (vperp,z,r,s)
+"""
+function gyroaverage_field!(gfield_out,field_in,gyro,vperp,z,r,composition)
+    @boundscheck z.n == size(field_in, 1) || throw(BoundsError(field_in))
+    @boundscheck r.n == size(field_in, 2) || throw(BoundsError(field_in))
+    @boundscheck vperp.n == size(gfield_out, 1) || throw(BoundsError(gfield))
+    @boundscheck z.n == size(gfield_out, 2) || throw(BoundsError(gfield))
+    @boundscheck r.n == size(gfield_out, 3) || throw(BoundsError(gfield))
+    @boundscheck composition.n_ion_species == size(gfield_out, 4) || throw(BoundsError(gfield))
+    
+    nr = r.n
+    nz = z.n
+    gyromatrix = gyro.gyromatrix
+    gyroloopsizes = gyro.gyroloopsizes
+    izpgyroindex = gyro.izpgyroindex
+    irpgyroindex = gyro.irpgyroindex
+    
+    begin_s_r_z_vperp_region()
+    @loop_s_r_z_vperp is ir iz ivperp begin
+        nsum = gyroloopsizes[ivperp,iz,ir,is]
+        @views izplist = izpgyroindex[1:nsum,ivperp,iz,ir,is]
+        @views irplist = irpgyroindex[1:nsum,ivperp,iz,ir,is]
+        
+        gfield_out[ivperp,iz,ir] = 0.0
+        # sum over all the contributions in the gyroaverage
+        # use list of indices here instead of simply summing over
+        # irp in 1:nr and izp in 1:nz to try to make use of the sparsity of the gyromatrix
+        for isum in 1:nsum
+            izp, irp = izplist[isum], irplist[isum]
+            gfield_out[ivperp,iz,ir,is] += gyromatrix[izp,irp,ivperp,iz,ir,is]*field_in[izp,irp]
+        end
+    end
+
+end
+
+"""
+function for gyroaveraging a charge particle pdf of shape (vpa,vperp,z,r,s)
+and filling the result into an of the same shape
+"""
+function gyroaverage_pdf!(gpdf_out,pdf_in,gyro,vpa,vperp,z,r,composition)
+    @boundscheck vpa.n == size(pdf_in, 1) || throw(BoundsError(pdf_in))
+    @boundscheck vperp.n == size(pdf_in, 2) || throw(BoundsError(pdf_in))
+    @boundscheck z.n == size(pdf_in, 3) || throw(BoundsError(pdf_in))
+    @boundscheck r.n == size(pdf_in, 4) || throw(BoundsError(pdf_in))
+    @boundscheck composition.n_ion_species == size(pdf_in, 5) || throw(BoundsError(pdf_in))
+    @boundscheck vpa.n == size(gpdf_out, 1) || throw(BoundsError(gpdf_out))
+    @boundscheck vperp.n == size(gpdf_out, 2) || throw(BoundsError(gpdf_out))
+    @boundscheck z.n == size(gpdf_out, 3) || throw(BoundsError(gpdf_out))
+    @boundscheck r.n == size(gpdf_out, 4) || throw(BoundsError(gpdf_out))
+    @boundscheck composition.n_ion_species == size(gpdf_out, 5) || throw(BoundsError(gpdf_out))
+    
+    nr = r.n
+    nz = z.n
+    gyromatrix = gyro.gyromatrix
+    gyroloopsizes = gyro.gyroloopsizes
+    izpgyroindex = gyro.izpgyroindex
+    irpgyroindex = gyro.irpgyroindex
+    
+    begin_s_r_z_vperp_vpa_region()
+    @loop_s_r_z_vperp is ir iz ivperp begin
+        nsum = gyroloopsizes[ivperp,iz,ir,is]
+        @views izplist = izpgyroindex[1:nsum,ivperp,iz,ir,is]
+        @views irplist = irpgyroindex[1:nsum,ivperp,iz,ir,is]
+        @loop_vpa ivpa begin 
+           gpdf_out[ivpa,ivperp,iz,ir,is] = 0.0
+           # sum over all the contributions in the gyroaverage
+           # use list of indices here instead of simply summing over
+           # irp in 1:nr and izp in 1:nz to try to make use of the sparsity of the gyromatrix
+           for isum in 1:nsum
+               izp, irp = izplist[isum], irplist[isum]
+               gpdf_out[ivpa,ivperp,iz,ir,is] += gyromatrix[izp,irp,ivperp,iz,ir,is]*pdf_in[ivpa,ivperp,izp,irp,is]
+           end
+        end
+    end
+
+end
+
+
+end
diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl
index d0500c986..a0bf9151f 100644
--- a/moment_kinetics/src/initial_conditions.jl
+++ b/moment_kinetics/src/initial_conditions.jl
@@ -41,7 +41,7 @@ using MPI
 """
 struct pdf_substruct{n_distribution}
     norm::MPISharedArray{mk_float,n_distribution}
-    buffer::MPISharedArray{mk_float,n_distribution} # for collision operator terms when pdfs must be interpolated onto different velocity space grids
+    buffer::MPISharedArray{mk_float,n_distribution} # for collision operator terms when pdfs must be interpolated onto different velocity space grids, and for gyroaveraging
 end
 
 # struct of structs neatly contains i+n info?
@@ -127,10 +127,12 @@ Allocate arrays for pdfs
 function create_pdf(composition, r, z, vperp, vpa, vzeta, vr, vz)
     # allocate pdf arrays
     pdf_charged_norm = allocate_shared_float(vpa.n, vperp.n, z.n, r.n, composition.n_ion_species)
+    # buffer array is for ion-neutral collisions, not for storing charged pdf
     pdf_charged_buffer = allocate_shared_float(vpa.n, vperp.n, z.n, r.n, composition.n_neutral_species) # n.b. n_species is n_neutral_species here
     pdf_neutral_norm = allocate_shared_float(vz.n, vr.n, vzeta.n, z.n, r.n, composition.n_neutral_species)
+    # buffer array is for neutral-ion collisions, not for storing neutral pdf
     pdf_neutral_buffer = allocate_shared_float(vz.n, vr.n, vzeta.n, z.n, r.n, composition.n_ion_species)
-
+    
     return pdf_struct(pdf_substruct(pdf_charged_norm, pdf_charged_buffer),
                       pdf_substruct(pdf_neutral_norm, pdf_neutral_buffer))
 
diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl
index 8e8e67bbe..6c3c848d8 100644
--- a/moment_kinetics/src/input_structs.jl
+++ b/moment_kinetics/src/input_structs.jl
@@ -278,6 +278,10 @@ mutable struct species_composition
     # The ion flux reaching the wall that is recycled as neutrals is reduced by
     # `recycling_fraction` to account for ions absorbed by the wall.
     recycling_fraction::mk_float
+    # gyrokinetic_ions is a flag determining if the ion species is gyrokinetic
+    # gyrokinetic_ions = true -> use gyroaveraged fields at fixed guiding centre and moments of the pdf computed at fixed r
+    # gyrokinetic_ions = false -> use drift kinetic approximation
+    gyrokinetic_ions::Bool
     # scratch buffer whose size is n_species
     scratch::Vector{mk_float}
 end
diff --git a/moment_kinetics/src/moment_kinetics.jl b/moment_kinetics/src/moment_kinetics.jl
index 1296bcf34..35f939cab 100644
--- a/moment_kinetics/src/moment_kinetics.jl
+++ b/moment_kinetics/src/moment_kinetics.jl
@@ -33,6 +33,7 @@ include("reference_parameters.jl")
 include("coordinates.jl")
 include("file_io.jl")
 include("geo.jl")
+include("gyroaverages.jl")
 include("velocity_moments.jl")
 include("velocity_grid_transforms.jl")
 include("em_fields.jl")
@@ -399,8 +400,8 @@ function setup_moment_kinetics(input_dict::AbstractDict;
     # the main time advance loop -- including normalisation of f by density if requested
 
     moments, fields, spectral_objects, advect_objects,
-    scratch, advance, fp_arrays, scratch_dummy, manufactured_source_list =
-        setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral,
+    scratch, advance, fp_arrays, gyroavs, scratch_dummy, manufactured_source_list =
+        setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, gyrophase, vz_spectral,
             vr_spectral, vzeta_spectral, vpa_spectral, vperp_spectral, z_spectral,
             r_spectral, composition, drive_input, moments, t_input, collisions, species,
             geometry, boundary_distributions, external_source_settings, num_diss_params,
@@ -430,7 +431,7 @@ function setup_moment_kinetics(input_dict::AbstractDict;
 
     return pdf, scratch, code_time, t_input, vz, vr, vzeta, vpa, vperp, gyrophase, z, r,
            moments, fields, spectral_objects, advect_objects,
-           composition, collisions, geometry, boundary_distributions,
+           composition, collisions, geometry, gyroavs, boundary_distributions,
            external_source_settings, num_diss_params, advance, fp_arrays, scratch_dummy,
            manufactured_source_list, ascii_io, io_moments, io_dfns
 end
diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl
index 570740787..f77368236 100644
--- a/moment_kinetics/src/moment_kinetics_input.jl
+++ b/moment_kinetics/src/moment_kinetics_input.jl
@@ -107,7 +107,10 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     if !(0.0 <= composition.recycling_fraction <= 1.0)
         error("recycling_fraction must be between 0 and 1. Got $recycling_fraction.")
     end
-
+    # gyrokinetic_ions = True -> use gyroaveraged fields at fixed guiding centre and moments of the pdf computed at fixed r
+    # gyrokinetic_ions = False -> use drift kinetic approximation
+    composition.gyrokinetic_ions = get(scan_input, "gyrokinetic_ions", false)
+    
     # Reference parameters that define the conversion between physical quantities and
     # normalised values used in the code.
     reference_params = setup_reference_parameters(scan_input)
@@ -936,9 +939,10 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # The ion flux reaching the wall that is recycled as neutrals is reduced by
     # `recycling_fraction` to account for ions absorbed by the wall.
     recycling_fraction = 1.0
+    gyrokinetic_ions = false
     composition = species_composition(n_species, n_ion_species, n_neutral_species,
         electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant,
-        mn_over_mi, me_over_mi, recycling_fraction, allocate_float(n_species))
+        mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species))
     
     species_charged = Array{species_parameters_mutable,1}(undef,n_ion_species)
     species_neutral = Array{species_parameters_mutable,1}(undef,n_neutral_species)
diff --git a/moment_kinetics/src/moment_kinetics_structs.jl b/moment_kinetics/src/moment_kinetics_structs.jl
index f2c16d4b0..25d23be1f 100644
--- a/moment_kinetics/src/moment_kinetics_structs.jl
+++ b/moment_kinetics/src/moment_kinetics_structs.jl
@@ -35,6 +35,12 @@ struct em_fields_struct
     Er::MPISharedArray{mk_float,2}
     # Ez is the parallel electric field
     Ez::MPISharedArray{mk_float,2}
+    # gphi is the gyroaveraged electrostatic potential
+    gphi::MPISharedArray{mk_float,4}
+    # gEr is the gyroaveraged radial electric field
+    gEr::MPISharedArray{mk_float,4}
+    # gEz is the gyroaveraged parallel electric field
+    gEz::MPISharedArray{mk_float,4}
     # if including an external forcing for phi, it is of the form
     # phi_external = phi0*drive_amplitude*sinpi(t*drive_frequency)
     force_phi::Bool
diff --git a/moment_kinetics/src/r_advection.jl b/moment_kinetics/src/r_advection.jl
index d205b04a1..e3f23aef7 100644
--- a/moment_kinetics/src/r_advection.jl
+++ b/moment_kinetics/src/r_advection.jl
@@ -22,7 +22,7 @@ function r_advection!(f_out, fvec_in, moments, fields, advect, r, z, vperp, vpa,
         # get the updated speed along the r direction using the current f
         @views update_speed_r!(advect[is], fvec_in.upar[:,:,is],
                                moments.charged.vth[:,:,is], fields, moments.evolve_upar,
-                               moments.evolve_ppar, vpa, vperp, z, r, geometry)
+                               moments.evolve_ppar, vpa, vperp, z, r, geometry, is)
         # advance r-advection equation
         @loop_z_vpa iz ivpa begin
         end
@@ -84,7 +84,7 @@ end
 calculate the advection speed in the r-direction at each grid point
 """
 function update_speed_r!(advect, upar, vth, fields, evolve_upar, evolve_ppar, vpa, vperp,
-                         z, r, geometry)
+                         z, r, geometry, is)
     @boundscheck z.n == size(advect.speed,4) || throw(BoundsError(advect))
     @boundscheck vperp.n == size(advect.speed,3) || throw(BoundsError(advect))
     @boundscheck vpa.n == size(advect.speed,2) || throw(BoundsError(advect))
@@ -105,7 +105,7 @@ function update_speed_r!(advect, upar, vth, fields, evolve_upar, evolve_ppar, vp
             @loop_z_vperp_vpa iz ivperp ivpa begin
                 # ExB drift
                 @. geofac = bzeta[iz,:]*jacobian[iz,:]/Bmag[iz,:]
-                @views @. advect.speed[:,ivpa,ivperp,iz] = ExBfac*geofac*fields.Ez[iz,:]
+                @views @. advect.speed[:,ivpa,ivperp,iz] = ExBfac*geofac*fields.gEz[ivperp,iz,:,is]
                 # magnetic curvature drift
                 @. @views advect.speed[:,ivpa,ivperp,iz] += rhostar*(vpa.grid[ivpa]^2)*cvdriftr[iz,:]
                 # magnetic grad B drift
diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl
index acfcb6a43..a3a46e3ea 100644
--- a/moment_kinetics/src/time_advance.jl
+++ b/moment_kinetics/src/time_advance.jl
@@ -55,6 +55,7 @@ using ..force_balance: force_balance!, neutral_force_balance!
 using ..energy_equation: energy_equation!, neutral_energy_equation!
 using ..em_fields: setup_em_fields, update_phi!
 using ..fokker_planck: init_fokker_planck_collisions_weak_form, explicit_fokker_planck_collisions_weak_form!
+using ..gyroaverages: init_gyro_operators, gyroaverage_pdf!
 using ..manufactured_solns: manufactured_sources
 using ..advection: advection_info
 using ..utils: to_minutes
@@ -180,7 +181,7 @@ this includes creating and populating structs
 for Chebyshev transforms, velocity space moments,
 EM fields, and advection terms
 """
-function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral,
+function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, gyrophase, vz_spectral,
                              vr_spectral, vzeta_spectral, vpa_spectral, vperp_spectral,
                              z_spectral, r_spectral, composition, drive_input, moments,
                              t_input, collisions, species, geometry,
@@ -216,12 +217,16 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral,
     else
         fp_arrays = nothing
     end
+    # create gyroaverage matrix arrays
+    gyroavs = init_gyro_operators(vperp,z,r,gyrophase,geometry,composition)
+
     # create the "fields" structure that contains arrays
     # for the electrostatic potential phi and eventually the electromagnetic fields
-    fields = setup_em_fields(z.n, r.n, drive_input.force_phi, drive_input.amplitude, drive_input.frequency, drive_input.force_Er_zero_at_wall)
+    fields = setup_em_fields(vperp.n, z.n, r.n, composition.n_ion_species, drive_input.force_phi,
+      drive_input.amplitude, drive_input.frequency, drive_input.force_Er_zero_at_wall)
     # initialize the electrostatic potential
     begin_serial_region()
-    update_phi!(fields, scratch[1], z, r, composition, z_spectral, r_spectral, scratch_dummy)
+    update_phi!(fields, scratch[1], vperp, z, r, composition, z_spectral, r_spectral, scratch_dummy, gyroavs)
     @serial_region begin
         # save the initial phi(z) for possible use later (e.g., if forcing phi)
         fields.phi0 .= fields.phi
@@ -250,7 +255,7 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral,
         @loop_s is begin
             @views update_speed_r!(r_advect[is], moments.charged.upar[:,:,is],
                                    moments.charged.vth[:,:,is], fields, moments.evolve_upar,
-                                   moments.evolve_ppar, vpa, vperp, z, r, geometry)
+                                   moments.evolve_ppar, vpa, vperp, z, r, geometry, is)
         end
         # enforce prescribed boundary condition in r on the distribution function f
     end
@@ -267,7 +272,7 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral,
             @views update_speed_z!(z_advect[is], moments.charged.upar[:,:,is],
                                    moments.charged.vth[:,:,is], moments.evolve_upar,
                                    moments.evolve_ppar, fields, vpa, vperp, z, r, 0.0,
-                                   geometry)
+                                   geometry, is)
         end
     end
     begin_serial_region()
@@ -388,10 +393,8 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral,
         # update moments in case they were affected by applying boundary conditions or
         # constraints to the pdf
         reset_moments_status!(moments)
-        update_moments!(moments, pdf.charged.norm, vpa, vperp, z, r, composition)
-        # update the Chodura diagnostic -- note that the pdf should be the unnormalised one
-        # so this will break for the split moments cases
-        update_chodura!(moments,pdf.charged.norm,vpa,vperp,z,r,r_spectral,composition,geometry,scratch_dummy,z_advect)
+        update_moments!(moments, pdf.charged.norm, gyroavs, vpa, vperp, z, r, composition,
+           r_spectral,geometry,scratch_dummy,z_advect)
         # enforce boundary conditions in r and z on the neutral particle distribution function
         if n_neutral_species > 0
             # Note, so far vr and vzeta do not need advect objects, so pass `nothing` for
@@ -431,8 +434,8 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral,
         end
     end
 
-    update_phi!(fields, scratch[1], z, r, composition, z_spectral, r_spectral,
-                scratch_dummy)
+    update_phi!(fields, scratch[1], vperp, z, r, composition, z_spectral, r_spectral,
+                scratch_dummy, gyroavs)
     calculate_moment_derivatives!(moments, scratch[1], scratch_dummy, z, z_spectral, num_diss_params)
     calculate_moment_derivatives_neutral!(moments, scratch[1], scratch_dummy, z,
                                           z_spectral, num_diss_params)
@@ -441,7 +444,7 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral,
     _block_synchronize()
 
     return moments, fields, spectral_objects, advect_objects,
-    scratch, advance, fp_arrays, scratch_dummy, manufactured_source_list
+    scratch, advance, fp_arrays, gyroavs, scratch_dummy, manufactured_source_list
 end
 
 """
@@ -824,7 +827,7 @@ time integrator can be used without severe CFL condition
 """
 function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase, z, r,
            moments, fields, spectral_objects, advect_objects,
-           composition, collisions, geometry, boundary_distributions, 
+           composition, collisions, geometry, gyroavs, boundary_distributions, 
            external_source_settings, num_diss_params, advance, fp_arrays, scratch_dummy,
            manufactured_source_list, ascii_io, io_moments, io_dfns)
 
@@ -858,6 +861,13 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro
     iwrite_dfns = 2
     finish_now = false
     for i ∈ 1:t_input.nstep
+        
+        if i == t_input.nstep
+            # Ensure all output is written at the final step
+            finish_now = true
+        end
+        diagnostic_checks = mod(i,t_input.nwrite_moments) == 0 || mod(i,t_input.nwrite_dfns) == 0 || finish_now
+        
         if t_input.split_operators
             # MRH NOT SUPPORTED
             time_advance_split_operators!(pdf, scratch, t, t_input, vpa, z,
@@ -867,22 +877,14 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro
         else
             time_advance_no_splitting!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase, z, r,
                 moments, fields, spectral_objects, advect_objects,
-                composition, collisions, geometry, boundary_distributions,
+                composition, collisions, geometry, gyroavs, boundary_distributions,
                 external_source_settings, num_diss_params, advance, fp_arrays,  scratch_dummy,
-                manufactured_source_list, i)
+                manufactured_source_list, diagnostic_checks, i)
         end
         # update the time
         t += t_input.dt
 
-        if i == t_input.nstep
-            # Ensure all output is written at the final step
-            finish_now = true
-        end
-
-        if mod(i,t_input.nwrite_moments) == 0 || mod(i,t_input.nwrite_dfns) == 0 || finish_now
-            # update the diagnostic chodura condition
-            update_chodura!(moments,scratch[end].pdf,vpa,vperp,z,r,spectral_objects.r_spectral,composition,geometry,scratch_dummy,advect_objects.z_advect)
-
+        if diagnostic_checks
             # Always synchronise here, regardless of if we changed region or not
             begin_serial_region(no_synchronize=true)
             _block_synchronize()
@@ -1206,15 +1208,15 @@ end
 """
 function time_advance_no_splitting!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase, z, r,
            moments, fields, spectral_objects, advect_objects,
-           composition, collisions, geometry, boundary_distributions,
+           composition, collisions, geometry, gyroavs, boundary_distributions,
            external_source_settings, num_diss_params, advance, fp_arrays, scratch_dummy,
-           manufactured_source_list, istep)
+           manufactured_source_list, diagnostic_checks, istep)
 
     if t_input.n_rk_stages > 1
         ssp_rk!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase, z, r,
             moments, fields, spectral_objects, advect_objects, composition, collisions,
-            geometry, boundary_distributions, external_source_settings, num_diss_params,
-            advance, fp_arrays, scratch_dummy, manufactured_source_list, istep)
+            geometry, gyroavs, boundary_distributions, external_source_settings, num_diss_params,
+            advance, fp_arrays, scratch_dummy, manufactured_source_list, diagnostic_checks, istep)
     else
         euler_time_advance!(scratch, scratch, pdf, fields, moments,
             advect_objects, vz, vr, vzeta, vpa, vperp, gyrophase, z, r, t,
@@ -1237,7 +1239,7 @@ or update them by taking the appropriate velocity moment of the evolved pdf
 """
 function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, vr, vzeta,
                     vpa, vperp, z, r, spectral_objects, advect_objects, rk_coefs, istage, composition,
-                    geometry, num_diss_params, advance, scratch_dummy)
+                    geometry, gyroavs, num_diss_params, advance, scratch_dummy, diagnostic_moments)
     begin_s_r_z_region()
 
     new_scratch = scratch[istage+1]
@@ -1285,30 +1287,9 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v
         end
     end
     # update remaining velocity moments that are calculable from the evolved pdf
-    update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composition)
-    # update the diagnostic chodura condition
-    # update_chodura!(moments,new_scratch.pdf,vpa,vperp,z,r,r_spectral,composition,geometry,scratch_dummy,z_advect)
-    # update the thermal speed
-    begin_s_r_z_region()
-    try #below block causes DomainError if ppar < 0 or density, so exit cleanly if possible
-        update_vth!(moments.charged.vth, new_scratch.ppar, new_scratch.pperp, new_scratch.density, vperp, z, r, composition)
-    catch e
-        if global_size[] > 1
-            println("ERROR: error calculating vth in time_advance.jl")
-            println(e)
-            display(stacktrace(catch_backtrace()))
-            flush(stdout)
-            flush(stderr)
-            MPI.Abort(comm_world, 1)
-        end
-        rethrow(e)
-    end
-    # update the parallel heat flux
-    update_qpar!(moments.charged.qpar, moments.charged.qpar_updated, new_scratch.density,
-                 new_scratch.upar, moments.charged.vth, new_scratch.pdf, vpa, vperp, z, r,
-                 composition, moments.evolve_density, moments.evolve_upar,
-                 moments.evolve_ppar)
-
+    update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composition,
+        r_spectral, geometry, gyroavs, scratch_dummy, z_advect, diagnostic_moments)
+    
     calculate_moment_derivatives!(moments, new_scratch, scratch_dummy, z, z_spectral,
                                   num_diss_params)
 
@@ -1376,7 +1357,7 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v
     end
 
     # update the electrostatic potential phi
-    update_phi!(fields, scratch[istage+1], z, r, composition, z_spectral, r_spectral, scratch_dummy)
+    update_phi!(fields, scratch[istage+1], vperp, z, r, composition, z_spectral, r_spectral, scratch_dummy, gyroavs)
     if !(( moments.evolve_upar || moments.evolve_ppar) &&
               istage == length(scratch)-1)
         # _block_synchronize() here because phi needs to be read on different ranks than
@@ -1441,23 +1422,61 @@ end
 """
 update velocity moments that are calculable from the evolved charged pdf
 """
-function update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composition)
+function update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composition,
+    r_spectral, geometry, gyroavs, scratch_dummy, z_advect, diagnostic_moments)
+    
+    if composition.gyrokinetic_ions
+        ff = scratch_dummy.buffer_vpavperpzrs_1
+        # fill buffer with ring-averaged F (gyroaverage at fixed position)
+        gyroaverage_pdf!(ff,new_scratch.pdf,gyroavs,vpa,vperp,z,r,composition)
+    else
+        ff = new_scratch.pdf
+    end
+    
     if !moments.evolve_density
         update_density!(new_scratch.density, moments.charged.dens_updated,
-                        new_scratch.pdf, vpa, vperp, z, r, composition)
+                        ff, vpa, vperp, z, r, composition)
     end
     if !moments.evolve_upar
         update_upar!(new_scratch.upar, moments.charged.upar_updated, new_scratch.density,
-                     new_scratch.ppar, new_scratch.pdf, vpa, vperp, z, r, composition,
+                     new_scratch.ppar, ff, vpa, vperp, z, r, composition,
                      moments.evolve_density, moments.evolve_ppar)
     end
     if !moments.evolve_ppar
         # update_ppar! calculates (p_parallel/m_s N_e c_s^2) + (n_s/N_e)*(upar_s/c_s)^2 = (1/√π)∫d(vpa/c_s) (vpa/c_s)^2 * (√π f_s c_s / N_e)
         update_ppar!(new_scratch.ppar, moments.charged.ppar_updated, new_scratch.density,
-                     new_scratch.upar, new_scratch.pdf, vpa, vperp, z, r, composition,
+                     new_scratch.upar, ff, vpa, vperp, z, r, composition,
                      moments.evolve_density, moments.evolve_upar)
     end 
-    update_pperp!(new_scratch.pperp, new_scratch.pdf, vpa, vperp, z, r, composition)
+    update_pperp!(new_scratch.pperp, ff, vpa, vperp, z, r, composition)
+    
+    # if diagnostic time step/RK stage
+    # update the diagnostic chodura condition
+    if diagnostic_moments
+        update_chodura!(moments,ff,vpa,vperp,z,r,r_spectral,composition,geometry,scratch_dummy,z_advect)
+    end
+    # update the thermal speed
+    begin_s_r_z_region()
+    try #below block causes DomainError if ppar < 0 or density, so exit cleanly if possible
+        update_vth!(moments.charged.vth, new_scratch.ppar, new_scratch.pperp, new_scratch.density, vperp, z, r, composition)
+    catch e
+        if global_size[] > 1
+            println("ERROR: error calculating vth in time_advance.jl")
+            println(e)
+            display(stacktrace(catch_backtrace()))
+            flush(stdout)
+            flush(stderr)
+            MPI.Abort(comm_world, 1)
+        end
+        rethrow(e)
+    end
+    # update the parallel heat flux
+    update_qpar!(moments.charged.qpar, moments.charged.qpar_updated, new_scratch.density,
+                 new_scratch.upar, moments.charged.vth, ff, vpa, vperp, z, r,
+                 composition, moments.evolve_density, moments.evolve_upar,
+                 moments.evolve_ppar)
+    # add further moments to be computed here
+    
 end
 
 """
@@ -1487,8 +1506,8 @@ end
 """
 function ssp_rk!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase, z, r,
            moments, fields, spectral_objects, advect_objects, composition, collisions,
-           geometry, boundary_distributions, external_source_settings, num_diss_params,
-           advance, fp_arrays, scratch_dummy, manufactured_source_list,  istep)
+           geometry, gyroavs, boundary_distributions, external_source_settings, num_diss_params,
+           advance, fp_arrays, scratch_dummy, manufactured_source_list, diagnostic_checks, istep)
 
     begin_s_r_z_region()
 
@@ -1534,10 +1553,11 @@ function ssp_rk!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase,
             t_input, spectral_objects, composition,
             collisions, geometry, scratch_dummy, manufactured_source_list,
             external_source_settings, num_diss_params, advance, fp_arrays, istage)
+        diagnostic_moments = diagnostic_checks && istage == n_rk_stages
         @views rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, vr,
                           vzeta, vpa, vperp, z, r, spectral_objects, advect_objects,
                           advance.rk_coefs[:,istage], istage, composition, geometry,
-                          num_diss_params, advance, scratch_dummy)
+                          gyroavs, num_diss_params, advance, scratch_dummy, diagnostic_moments)
     end
 
     istage = n_rk_stages+1
diff --git a/moment_kinetics/src/velocity_moments.jl b/moment_kinetics/src/velocity_moments.jl
index 881502886..3da3e251f 100644
--- a/moment_kinetics/src/velocity_moments.jl
+++ b/moment_kinetics/src/velocity_moments.jl
@@ -41,7 +41,7 @@ using ..communication
 using ..derivatives: derivative_z!
 using ..derivatives: derivative_r!
 using ..looping
-
+using ..gyroaverages: gyro_operators, gyroaverage_pdf!
 #global tmpsum1 = 0.0
 #global tmpsum2 = 0.0
 #global dens_hist = zeros(17,1)
@@ -486,8 +486,18 @@ end
 
 """
 calculate the updated density (dens) and parallel pressure (ppar) for all species
-"""
-function update_moments!(moments, ff, vpa, vperp, z, r, composition)
+this function is only used once after initialisation
+the function used to update moments at run time is update_derived_moments! in time_advance.jl
+"""
+function update_moments!(moments, ff_in, gyroavs::gyro_operators, vpa, vperp, z, r, composition,
+        r_spectral, geometry, scratch_dummy, z_advect)
+    if composition.gyrokinetic_ions
+        ff = scratch_dummy.buffer_vpavperpzrs_1 # the buffer array for the charged pdf -> make sure not to reuse this array below
+        # fill buffer with ring-averaged F (gyroaverage at fixed position)
+        gyroaverage_pdf!(ff,ff_in,gyroavs,vpa,vperp,z,r,composition)
+    else
+        ff = ff_in
+    end
     begin_s_r_z_region()
     n_species = size(ff,5)
     @boundscheck n_species == size(moments.charged.dens,3) || throw(BoundsError(moments))
@@ -528,7 +538,12 @@ function update_moments!(moments, ff, vpa, vperp, z, r, composition)
             moments.charged.qpar_updated[is] = true
         end
     end
+
     update_vth!(moments.charged.vth, moments.charged.ppar, moments.charged.pperp, moments.charged.dens, vperp, z, r, composition)
+    # update the Chodura diagnostic -- note that the pdf should be the unnormalised one
+    # so this will break for the split moments cases
+    update_chodura!(moments,ff,vpa,vperp,z,r,r_spectral,composition,geometry,scratch_dummy,z_advect)
+        
     return nothing
 end
 
@@ -850,8 +865,9 @@ in a single species plasma with Z = 1
 function update_chodura!(moments,ff,vpa,vperp,z,r,r_spectral,composition,geometry,scratch_dummy,z_advect)
     @boundscheck composition.n_ion_species == size(ff, 5) || throw(BoundsError(ff))
     begin_s_z_vperp_vpa_region()
-    dffdr = scratch_dummy.buffer_vpavperpzrs_1
-    ff_dummy = scratch_dummy.buffer_vpavperpzrs_2
+    # use buffer_vpavperpzrs_2 here as buffer_vpavperpzrs_1 is in use storing ff
+    dffdr = scratch_dummy.buffer_vpavperpzrs_2 
+    ff_dummy = scratch_dummy.dummy_vpavperp
     if r.n > 1
     # first compute d f / d r using centred reconciliation and place in dummy array #1
     derivative_r!(dffdr, ff[:,:,:,:,:],
@@ -869,7 +885,7 @@ function update_chodura!(moments,ff,vpa,vperp,z,r,r_spectral,composition,geometr
     if z.irank == 0
         @loop_s_r is ir begin
             @views moments.charged.chodura_integral_lower[ir,is] = update_chodura_integral_species!(ff[:,:,1,ir,is],dffdr[:,:,1,ir,is],
-            ff_dummy[:,:,1,ir,is],vpa,vperp,z,r,composition,geometry,z_advect[is].speed[1,:,:,ir],moments.charged.dens[1,ir,is],del_vpa,1,ir)
+            ff_dummy[:,:],vpa,vperp,z,r,composition,geometry,z_advect[is].speed[1,:,:,ir],moments.charged.dens[1,ir,is],del_vpa,1,ir)
         end
     else # we do not save this Chodura integral to the output file
         @loop_s_r is ir begin
@@ -879,7 +895,7 @@ function update_chodura!(moments,ff,vpa,vperp,z,r,r_spectral,composition,geometr
     if z.irank == z.nrank - 1
         @loop_s_r is ir begin
             @views moments.charged.chodura_integral_upper[ir,is] = update_chodura_integral_species!(ff[:,:,end,ir,is],dffdr[:,:,end,ir,is],
-            ff_dummy[:,:,end,ir,is],vpa,vperp,z,r,composition,geometry,z_advect[is].speed[end,:,:,ir],moments.charged.dens[end,ir,is],del_vpa,z.n,ir)
+            ff_dummy[:,:],vpa,vperp,z,r,composition,geometry,z_advect[is].speed[end,:,:,ir],moments.charged.dens[end,ir,is],del_vpa,z.n,ir)
         end
     else # we do not save this Chodura integral to the output file
         @loop_s_r is ir begin
diff --git a/moment_kinetics/src/vpa_advection.jl b/moment_kinetics/src/vpa_advection.jl
index 47da74aff..ed8afd2fa 100644
--- a/moment_kinetics/src/vpa_advection.jl
+++ b/moment_kinetics/src/vpa_advection.jl
@@ -90,7 +90,7 @@ function update_speed_default!(advect, fields, fvec, moments, vpa, vperp, z, r,
                 # mu, the adiabatic invariant
                 mu = 0.5*(vperp.grid[ivperp]^2)/Bmag[iz,ir]
                 # bzed = B_z/B
-                advect[is].speed[ivpa,ivperp,iz,ir] = (0.5*bzed[iz,ir]*fields.Ez[iz,ir] - 
+                advect[is].speed[ivpa,ivperp,iz,ir] = (0.5*bzed[iz,ir]*fields.gEz[ivperp,iz,ir,is] - 
                                                        mu*bzed[iz,ir]*dBdz[iz,ir])
             end
         end
diff --git a/moment_kinetics/src/z_advection.jl b/moment_kinetics/src/z_advection.jl
index 8a4a3f273..b88808a3e 100644
--- a/moment_kinetics/src/z_advection.jl
+++ b/moment_kinetics/src/z_advection.jl
@@ -22,7 +22,7 @@ function z_advection!(f_out, fvec_in, moments, fields, advect, z, vpa, vperp, r,
         # get the updated speed along the z direction using the current f
         @views update_speed_z!(advect[is], fvec_in.upar[:,:,is],
                                moments.charged.vth[:,:,is], moments.evolve_upar,
-                               moments.evolve_ppar, fields, vpa, vperp, z, r, t, geometry)
+                               moments.evolve_ppar, fields, vpa, vperp, z, r, t, geometry, is)
         # update adv_fac
         @loop_r_vperp_vpa ir ivperp ivpa begin
             @views adjust_advection_speed!(advect[is].speed[:,ivpa,ivperp,ir],
@@ -85,7 +85,7 @@ end
 calculate the advection speed in the z-direction at each grid point
 """
 function update_speed_z!(advect, upar, vth, evolve_upar, evolve_ppar, fields, vpa, vperp,
-                         z, r, t, geometry)
+                         z, r, t, geometry, is)
     @boundscheck r.n == size(advect.speed,4) || throw(BoundsError(advect))
     @boundscheck vperp.n == size(advect.speed,3) || throw(BoundsError(advect))
     @boundscheck vpa.n == size(advect.speed,2) || throw(BoundsError(advect))
@@ -106,7 +106,7 @@ function update_speed_z!(advect, upar, vth, evolve_upar, evolve_ppar, fields, vp
                 # vpa bzed
                 @. @views advect.speed[:,ivpa,ivperp,ir] = vpa.grid[ivpa]*bzed[:,ir]
                 # ExB drift
-                @. @views advect.speed[:,ivpa,ivperp,ir] += ExBfac*bzeta[:,ir]*jacobian[:,ir]/Bmag[:,ir]*fields.Er[:,ir]
+                @. @views advect.speed[:,ivpa,ivperp,ir] += ExBfac*bzeta[:,ir]*jacobian[:,ir]/Bmag[:,ir]*fields.gEr[ivperp,:,ir,is]
                 # magnetic curvature drift
                 @. @views advect.speed[:,ivpa,ivperp,ir] += rhostar*(vpa.grid[ivpa]^2)*cvdriftz[:,ir]
                 # magnetic grad B drift
diff --git a/moment_kinetics/test/gyroaverage_tests.jl b/moment_kinetics/test/gyroaverage_tests.jl
new file mode 100644
index 000000000..8b1740d02
--- /dev/null
+++ b/moment_kinetics/test/gyroaverage_tests.jl
@@ -0,0 +1,315 @@
+module GyroAverageTests
+
+include("setup.jl")
+
+export gyroaverage_test
+
+using MPI
+using SpecialFunctions: besselj0
+
+import moment_kinetics
+using moment_kinetics.input_structs
+using moment_kinetics.coordinates: define_coordinate
+using moment_kinetics.geo: init_magnetic_geometry
+using moment_kinetics.communication
+using moment_kinetics.looping
+using moment_kinetics.array_allocation: allocate_float, allocate_shared_float
+using moment_kinetics.gyroaverages: gyroaverage_pdf!
+using moment_kinetics.gyroaverages: gyroaverage_field!, init_gyro_operators
+using moment_kinetics.type_definitions: mk_float, mk_int
+
+print_test_results = false
+
+function runtests()
+    # basic tests using periodic boundary conditions for the test function
+    # functionality to change rhostar and vperp_max = Lvperp and the pitch = bzed is provided
+    # the test could be sped up by moving lists of kz and kr inside the gyroaverage_test() function, 
+    # so that the gyromatrix does not have to be reinitialised
+
+    # Only needed to save FFTW 'wisdom'
+    test_output_directory = get_MPI_tempdir()
+
+    @testset "Gyroaverage tests" verbose=use_verbose begin
+        println("Gyroaverages test")
+        @testset " - test real-space path-integral gyroaverage (periodic functions)" begin
+            ngrid_vperp = 3
+            ngrid = 5; nelement = 4; ngrid_gyrophase = 100
+            z_bc = "periodic"; r_bc = "periodic"
+            kr = 1; kz = 1; Lvperp = 3.0; pitch = 0.5; rhostar = 0.1; phaser = 0.0; phasez = 0.0
+            @testset "kr $kr kz $kz vperpmax $Lvperp rhostar $rhostar phaser $phaser phasez $phasez" begin
+                absolute_error = 2.0e-4
+                gyroaverage_test(absolute_error; rhostar=rhostar, pitch=pitch, ngrid=ngrid, kr=kr, kz=kz, phaser=phaser, phasez=phasez, nelement=nelement, ngrid_vperp=ngrid_vperp, nelement_vperp=1, Lvperp=Lvperp, ngrid_gyrophase=ngrid_gyrophase, discretization="chebyshev_pseudospectral", r_bc=r_bc, z_bc=z_bc, test_output_directory=test_output_directory)
+            end
+            kr = 1; kz = 3; Lvperp = 3.0; pitch = 0.5; rhostar = 0.1; phaser = 0.0; phasez = 0.0
+            @testset "kr $kr kz $kz vperpmax $Lvperp rhostar $rhostar phaser $phaser phasez $phasez" begin
+                absolute_error = 2.0e-2
+                gyroaverage_test(absolute_error; rhostar=rhostar, pitch=pitch, ngrid=ngrid, kr=kr, kz=kz, phaser=phaser, phasez=phasez, nelement=nelement, ngrid_vperp=ngrid_vperp, nelement_vperp=1, Lvperp=Lvperp, ngrid_gyrophase=ngrid_gyrophase, discretization="chebyshev_pseudospectral", r_bc=r_bc, z_bc=z_bc, test_output_directory=test_output_directory)
+            end
+            kr = 3; kz = 1; Lvperp = 3.0; pitch = 0.5; rhostar = 0.1; phaser = 0.0; phasez = 0.0
+            @testset "kr $kr kz $kz vperpmax $Lvperp rhostar $rhostar phaser $phaser phasez $phasez" begin
+                absolute_error = 2.0e-2
+                gyroaverage_test(absolute_error; rhostar=rhostar, pitch=pitch, ngrid=ngrid, kr=kr, kz=kz, phaser=phaser, phasez=phasez, nelement=nelement, ngrid_vperp=ngrid_vperp, nelement_vperp=1, Lvperp=Lvperp, ngrid_gyrophase=ngrid_gyrophase, discretization="chebyshev_pseudospectral", r_bc=r_bc, z_bc=z_bc, test_output_directory=test_output_directory)
+            end
+            ngrid = 5; nelement = 8; ngrid_gyrophase = 100
+            z_bc = "periodic"; r_bc = "periodic"
+            kr = 1; kz = 1; Lvperp = 3.0; pitch = 0.5; rhostar = 0.1; phaser = 0.0; phasez = 0.0
+            @testset "kr $kr kz $kz vperpmax $Lvperp rhostar $rhostar phaser $phaser phasez $phasez" begin
+                absolute_error = 4.0e-6
+                gyroaverage_test(absolute_error; rhostar=rhostar, pitch=pitch, ngrid=ngrid, kr=kr, kz=kz, phaser=phaser, phasez=phasez, nelement=nelement, ngrid_vperp=ngrid_vperp, nelement_vperp=1, Lvperp=Lvperp, ngrid_gyrophase=ngrid_gyrophase, discretization="chebyshev_pseudospectral", r_bc=r_bc, z_bc=z_bc, test_output_directory=test_output_directory)
+            end
+            ngrid = 5; nelement = 8; ngrid_gyrophase = 100
+            z_bc = "periodic"; r_bc = "periodic"
+            kr = 1; kz = 3; Lvperp = 3.0; pitch = 0.5; rhostar = 0.1; phaser = 0.0; phasez = 0.0
+            @testset "kr $kr kz $kz vperpmax $Lvperp rhostar $rhostar phaser $phaser phasez $phasez" begin
+                absolute_error = 3.0e-3
+                gyroaverage_test(absolute_error; rhostar=rhostar, pitch=pitch, ngrid=ngrid, kr=kr, kz=kz, phaser=phaser, phasez=phasez, nelement=nelement, ngrid_vperp=ngrid_vperp, nelement_vperp=1, Lvperp=Lvperp, ngrid_gyrophase=ngrid_gyrophase, discretization="chebyshev_pseudospectral", r_bc=r_bc, z_bc=z_bc, test_output_directory=test_output_directory)
+            end
+            ngrid = 5; nelement = 8; ngrid_gyrophase = 100
+            z_bc = "periodic"; r_bc = "periodic"
+            kr = 3; kz = 1; Lvperp = 3.0; pitch = 0.5; rhostar = 0.1; phaser = 0.0; phasez = 0.0
+            @testset "kr $kr kz $kz vperpmax $Lvperp rhostar $rhostar phaser $phaser phasez $phasez" begin
+                absolute_error = 3.0e-3
+                gyroaverage_test(absolute_error; rhostar=rhostar, pitch=pitch, ngrid=ngrid, kr=kr, kz=kz, phaser=phaser, phasez=phasez, nelement=nelement, ngrid_vperp=ngrid_vperp, nelement_vperp=1, Lvperp=Lvperp, ngrid_gyrophase=ngrid_gyrophase, discretization="chebyshev_pseudospectral", r_bc=r_bc, z_bc=z_bc, test_output_directory=test_output_directory)
+            end
+        end
+    end
+end
+
+function gyroaverage_test(absolute_error; rhostar=0.1, pitch=0.5, ngrid=5, kr=2, kz=2, phaser=0.0, phasez=0.0, nelement=4, ngrid_vperp=3, nelement_vperp=1, Lvperp=3.0, ngrid_gyrophase=100, discretization="chebyshev_pseudospectral", r_bc="periodic", z_bc = "wall", print_test_results=print_test_results, test_output_directory)
+
+        #ngrid = 17
+        #nelement = 4
+        r_ngrid = ngrid #number of points per element 
+        r_nelement_local = nelement # number of elements per rank
+        r_nelement_global = r_nelement_local # total number of elements 
+        r_L = 1.0
+
+        z_ngrid = ngrid #number of points per element 
+        z_nelement_local = nelement # number of elements per rank
+        z_nelement_global = z_nelement_local # total number of elements 
+        z_L = 1.0
+
+        vperp_ngrid = ngrid_vperp #number of points per element 
+        vperp_nelement_local = nelement_vperp # number of elements per rank
+        vperp_nelement_global = vperp_nelement_local # total number of elements 
+        vperp_L = Lvperp
+        vperp_bc = "zero"
+        
+        vpa_ngrid = 1 #number of points per element 
+        vpa_nelement_local = 1 # number of elements per rank
+        vpa_nelement_global = vpa_nelement_local # total number of elements 
+        vpa_L = 1.0
+        vpa_bc = "" # should not be used
+        
+        gyrophase_ngrid = ngrid_gyrophase #number of points per element 
+        gyrophase_nelement_local = 1 # number of elements per rank
+        gyrophase_nelement_global = gyrophase_nelement_local # total number of elements 
+        gyrophase_discretization = "finite_difference"
+        gyrophase_L = 2.0*pi
+        
+        # fd_option and adv_input not actually used so given values unimportant
+        fd_option = "fourth_order_centered"
+        cheb_option = "matrix"
+        adv_input = advection_input("default", 1.0, 0.0, 0.0)
+        nrank = 1
+        irank = 0#1
+        comm = MPI.COMM_NULL
+        element_spacing_option = "uniform"
+        # create the 'input' struct containing input info needed to create a
+        # coordinate
+        
+        # Set up MPI
+        initialize_comms!()
+        setup_distributed_memory_MPI(1,1,1,1)
+        irank_z, nrank_z, comm_sub_z, irank_r, nrank_r, comm_sub_r = setup_distributed_memory_MPI(z_nelement_global,z_nelement_local,r_nelement_global,r_nelement_local)
+        
+        r_input = grid_input("r", r_ngrid, r_nelement_global, r_nelement_local, 
+                nrank_r, irank_r, r_L, discretization, fd_option, cheb_option, r_bc, adv_input,comm_sub_r,element_spacing_option)
+        z_input = grid_input("z", z_ngrid, z_nelement_global, z_nelement_local, 
+                nrank_z, irank_z, z_L, discretization, fd_option, cheb_option, z_bc, adv_input,comm_sub_z,element_spacing_option)
+        vperp_input = grid_input("vperp", vperp_ngrid, vperp_nelement_global, vperp_nelement_local, 
+                nrank, irank, vperp_L, discretization, fd_option, cheb_option, vperp_bc, adv_input,comm,element_spacing_option)
+        vpa_input = grid_input("vpa", vpa_ngrid, vpa_nelement_global, vpa_nelement_local, 
+                nrank, irank, vpa_L, discretization, fd_option, cheb_option, vpa_bc, adv_input,comm,element_spacing_option)
+        gyrophase_input = grid_input("gyrophase", gyrophase_ngrid, gyrophase_nelement_global, gyrophase_nelement_local, 
+                nrank, irank, gyrophase_L, gyrophase_discretization, fd_option, cheb_option, "periodic", adv_input,comm,element_spacing_option)
+        
+        # create the coordinate structs
+        r, r_spectral = define_coordinate(r_input; init_YY=false, run_directory=test_output_directory)
+        z, z_spectral = define_coordinate(z_input; init_YY=false, run_directory=test_output_directory)
+        vperp, vperp_spectral = define_coordinate(vperp_input; init_YY=false, run_directory=test_output_directory)
+        vpa, vpa_spectral = define_coordinate(vpa_input; init_YY=false, run_directory=test_output_directory)
+        gyrophase, gyrophase_spectral = define_coordinate(gyrophase_input; init_YY=false, run_directory=test_output_directory)
+        
+        # create test geometry
+        #rhostar = 0.1 #rhostar of ions for ExB drift
+        option = "constant-helical"
+        #pitch = 1.0
+        DeltaB = 1.0
+        geometry_in = geometry_input(rhostar,option,pitch,DeltaB)
+        geometry = init_magnetic_geometry(geometry_in,z,r)
+        
+        # create test composition
+        composition = create_test_composition()
+            
+        # setup shared-memory MPI ranges
+        looping.setup_loop_ranges!(block_rank[], block_size[];
+                                       s=composition.n_ion_species, sn=1,
+                                       r=r.n, z=z.n, vperp=vperp.n, vpa=vpa.n,
+                                       vzeta=1, vr=1, vz=1)
+                                       
+        # initialise the matrix for the gyroaverages
+        gyro = init_gyro_operators(vperp,z,r,gyrophase,geometry,composition,print_info=print_test_results)
+        # initialise a test field
+        phi = allocate_shared_float(z.n,r.n)
+        gphi = allocate_shared_float(vperp.n,z.n,r.n,composition.n_ion_species)
+        gphi_exact = allocate_float(vperp.n,z.n,r.n,composition.n_ion_species)
+        gphi_err = allocate_float(vperp.n,z.n,r.n)
+        begin_serial_region()
+        @serial_region begin 
+            fill_test_arrays!(phi,gphi_exact,vperp,z,r,geometry,composition,kz,kr,phasez,phaser)
+        end
+        
+        # gyroaverage phi
+        gyroaverage_field!(gphi,phi,gyro,vperp,z,r,composition)
+        
+        # compute errors
+        begin_serial_region()
+        @serial_region begin
+            @. gphi_err = abs(gphi - gphi_exact)
+            if print_test_results
+                println("Test gyroaverage_field!()")
+            end
+            for ivperp in 1:vperp.n
+                if print_test_results
+                    println("ivperp: ",ivperp," max(abs(gphi_err)): ",maximum(gphi_err[ivperp,:,:])," max(abs(gphi)): ",maximum(gphi[ivperp,:,:]))
+                end
+                @test maximum(gphi_err[ivperp,:,:]) < absolute_error
+            end
+            if print_test_results
+                println("")
+            end
+        end
+        
+        # repeat the test for a pdf
+        # initialise a test field
+        n_ion_species = composition.n_ion_species
+        pdf = allocate_shared_float(vpa.n,vperp.n,z.n,r.n,n_ion_species)
+        gpdf = allocate_shared_float(vpa.n,vperp.n,z.n,r.n,n_ion_species)
+        gpdf_exact = allocate_float(vpa.n,vperp.n,z.n,r.n,n_ion_species)
+        gpdf_err = allocate_float(vpa.n,vperp.n,z.n,r.n,n_ion_species)
+        begin_serial_region()
+        @serial_region begin
+            fill_pdf_test_arrays!(pdf,gpdf_exact,vpa,vperp,z,r,composition,geometry,kz,kr,phasez,phaser)
+        end
+        
+        gyroaverage_pdf!(gpdf,pdf,gyro,vpa,vperp,z,r,composition)
+        # compute errors
+        begin_serial_region()
+        @serial_region begin
+            @. gpdf_err = abs(gpdf - gpdf_exact)
+            if print_test_results
+                println("Test gyroaverage_pdf!()")
+            end
+            for is in 1:n_ion_species
+                for ivperp in 1:vperp.n
+                    for ivpa in 1:vpa.n
+                        if print_test_results
+                            println("ivpa: ",ivpa," ivperp: ",ivperp," is: ",is," max(abs(gpdf_err)): ",maximum(gpdf_err[ivpa,ivperp,:,:,is]),
+                                " max(abs(gpdf)): ",maximum(gpdf[ivpa,ivperp,:,:,is]), " max(abs(gpdf_exact)): ",maximum(gpdf_exact[ivpa,ivperp,:,:,is]))
+                        end
+                        @test maximum(gpdf_err[ivpa,ivperp,:,:,is]) < absolute_error
+                    end
+                end
+            end
+            if print_test_results
+                println("")
+            end
+        end
+        
+        finalize_comms!()
+end
+
+function create_test_composition()
+    electron_physics = boltzmann_electron_response
+    n_ion_species = 1
+    n_neutral_species = 0
+    n_species = n_ion_species + n_neutral_species
+    use_test_neutral_wall_pdf = false
+    # electron temperature over reference temperature
+    T_e = 1.0
+    # temperature at the entrance to the wall in terms of the electron temperature
+    T_wall = 1.0
+    # wall potential at z = 0
+    phi_wall = 0.0
+    # constant to test nonzero Er
+    Er_constant = 0.0
+    # ratio of the neutral particle mass to the ion particle mass
+    mn_over_mi = 1.0
+    # ratio of the electron particle mass to the ion particle mass
+    me_over_mi = 1.0/1836.0
+    # The ion flux reaching the wall that is recycled as neutrals is reduced by
+    # `recycling_fraction` to account for ions absorbed by the wall.
+    recycling_fraction = 1.0
+    gyrokinetic_ions = true
+    return composition = species_composition(n_species, n_ion_species, n_neutral_species,
+            electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant,
+            mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species))
+end
+
+function fill_test_arrays!(phi,gphi,vperp,z,r,geometry,composition,kz,kr,phasez,phaser)
+   n_ion_species = composition.n_ion_species
+   for ir in 1:r.n
+      for iz in 1:z.n
+         Bmag = geometry.Bmag[iz,ir] 
+         bzeta = geometry.bzeta[iz,ir] 
+         rhostar = geometry.rhostar 
+         # convert integer "wavenumbers" to actual wavenumbers 
+         kkr = 2.0*pi*kr/r.L
+         kkz = 2.0*pi*kz/z.L
+         kperp = sqrt(kkr^2 + (bzeta*kkz)^2)
+         
+         phi[iz,ir] = sin(r.grid[ir]*kkr + phaser)*sin(z.grid[iz]*kkz + phasez)
+         for is in 1:n_ion_species
+             for ivperp in 1:vperp.n
+                krho = kperp*vperp.grid[ivperp]*rhostar/Bmag
+                # use that phi is a sum of Kronecker deltas in k space to write gphi
+                gphi[ivperp,iz,ir,is] = besselj0(krho)*phi[iz,ir]
+             end
+         end
+      end
+   end  
+   return nothing
+end
+
+function fill_pdf_test_arrays!(pdf,gpdf,vpa,vperp,z,r,composition,geometry,kz,kr,phasez,phaser)
+   for is in 1:composition.n_ion_species
+       for ir in 1:r.n
+          for iz in 1:z.n
+             Bmag = geometry.Bmag[iz,ir] 
+             bzeta = geometry.bzeta[iz,ir] 
+             rhostar = geometry.rhostar 
+             # convert integer "wavenumbers" to actual wavenumbers 
+             kkr = 2.0*pi*kr/r.L
+             kkz = 2.0*pi*kz/z.L
+             kperp = sqrt(kkr^2 + (bzeta*kkz)^2)
+             
+             for ivperp in 1:vperp.n
+                for ivpa in 1:vpa.n
+                    pdf[ivpa,ivperp,iz,ir,is] = sin(r.grid[ir]*kkr + phaser)*sin(z.grid[iz]*kkz + phasez)
+                    krho = kperp*vperp.grid[ivperp]*rhostar/Bmag
+                    # use that pdf is a sum of Kronecker deltas in k space to write gpdf
+                    gpdf[ivpa,ivperp,iz,ir,is] = besselj0(krho)*pdf[ivpa,ivperp,iz,ir,is]
+                end
+             end
+          end
+       end
+   end  
+   return nothing
+end
+
+end #GyroAverageTests
+
+using .GyroAverageTests
+
+GyroAverageTests.runtests()
diff --git a/moment_kinetics/test/runtests.jl b/moment_kinetics/test/runtests.jl
index 26a8b4e67..1b78fca36 100644
--- a/moment_kinetics/test/runtests.jl
+++ b/moment_kinetics/test/runtests.jl
@@ -17,6 +17,7 @@ function runtests()
         include(joinpath(@__DIR__, "recycling_fraction_tests.jl"))
         include(joinpath(@__DIR__, "fokker_planck_tests.jl"))
         include(joinpath(@__DIR__, "fokker_planck_time_evolution_tests.jl"))
+        include(joinpath(@__DIR__, "gyroaverage_tests.jl"))
     end
 end
 
diff --git a/test_scripts/gyroaverage_test.jl b/test_scripts/gyroaverage_test.jl
new file mode 100644
index 000000000..eb1c3bcb6
--- /dev/null
+++ b/test_scripts/gyroaverage_test.jl
@@ -0,0 +1,276 @@
+export gyroaverage_test
+
+using Printf
+using Plots
+using LaTeXStrings
+using MPI
+using Measures
+using SpecialFunctions: besselj0
+
+import moment_kinetics
+using moment_kinetics.input_structs
+using moment_kinetics.coordinates: define_coordinate
+using moment_kinetics.geo: init_magnetic_geometry
+using moment_kinetics.communication
+using moment_kinetics.looping
+using moment_kinetics.array_allocation: allocate_float, allocate_shared_float
+using moment_kinetics.gyroaverages: gyroaverage_pdf!
+using moment_kinetics.gyroaverages: gyroaverage_field!, init_gyro_operators
+using moment_kinetics.type_definitions: mk_float, mk_int
+
+function print_matrix(matrix,name::String,n::mk_int,m::mk_int)
+        println("\n ",name," \n")
+        for i in 1:n
+            for j in 1:m
+                @printf("%.2f ", matrix[i,j])
+            end
+            println("")
+        end
+        println("\n")
+    end
+    
+    function print_vector(vector,name::String,m::mk_int)
+        println("\n ",name," \n")
+        for j in 1:m
+            @printf("%.3f ", vector[j])
+        end
+        println("")
+        println("\n")
+    end
+
+
+function gyroaverage_test(;rhostar=0.1, pitch=0.5, ngrid=5, kr=2, kz=2, phaser=0.0, phasez=0.0, nelement=4, ngrid_vperp=3, nelement_vperp=1, Lvperp=3.0, ngrid_gyrophase=100, discretization="chebyshev_pseudospectral", r_bc="periodic", z_bc = "wall", plot_test_results=false)
+
+        #ngrid = 17
+        #nelement = 4
+        r_ngrid = ngrid #number of points per element 
+        r_nelement_local = nelement # number of elements per rank
+        r_nelement_global = r_nelement_local # total number of elements 
+        r_L = 1.0
+
+        z_ngrid = ngrid #number of points per element 
+        z_nelement_local = nelement # number of elements per rank
+        z_nelement_global = z_nelement_local # total number of elements 
+        z_L = 1.0
+
+        vperp_ngrid = ngrid_vperp #number of points per element 
+        vperp_nelement_local = nelement_vperp # number of elements per rank
+        vperp_nelement_global = vperp_nelement_local # total number of elements 
+        vperp_L = Lvperp
+        vperp_bc = "zero"
+        
+        vpa_ngrid = 1 #number of points per element 
+        vpa_nelement_local = 1 # number of elements per rank
+        vpa_nelement_global = vpa_nelement_local # total number of elements 
+        vpa_L = 1.0
+        vpa_bc = "" # should not be used
+        
+        gyrophase_ngrid = ngrid_gyrophase #number of points per element 
+        gyrophase_nelement_local = 1 # number of elements per rank
+        gyrophase_nelement_global = gyrophase_nelement_local # total number of elements 
+        gyrophase_discretization = "finite_difference"
+        gyrophase_L = 2.0*pi
+        
+        # fd_option and adv_input not actually used so given values unimportant
+        fd_option = "fourth_order_centered"
+        cheb_option = "matrix"
+        adv_input = advection_input("default", 1.0, 0.0, 0.0)
+        nrank = 1
+        irank = 0#1
+        comm = MPI.COMM_NULL
+        element_spacing_option = "uniform"
+        # create the 'input' struct containing input info needed to create a
+        # coordinate
+        
+        r_input = grid_input("r", r_ngrid, r_nelement_global, r_nelement_local, 
+                nrank, irank, r_L, discretization, fd_option, cheb_option, r_bc, adv_input,comm,element_spacing_option)
+        z_input = grid_input("z", z_ngrid, z_nelement_global, z_nelement_local, 
+                nrank, irank, z_L, discretization, fd_option, cheb_option, z_bc, adv_input,comm,element_spacing_option)
+        vperp_input = grid_input("vperp", vperp_ngrid, vperp_nelement_global, vperp_nelement_local, 
+                nrank, irank, vperp_L, discretization, fd_option, cheb_option, vperp_bc, adv_input,comm,element_spacing_option)
+        vpa_input = grid_input("vpa", vpa_ngrid, vpa_nelement_global, vpa_nelement_local, 
+                nrank, irank, vpa_L, discretization, fd_option, cheb_option, vpa_bc, adv_input,comm,element_spacing_option)
+        gyrophase_input = grid_input("gyrophase", gyrophase_ngrid, gyrophase_nelement_global, gyrophase_nelement_local, 
+                nrank, irank, gyrophase_L, gyrophase_discretization, fd_option, cheb_option, "periodic", adv_input,comm,element_spacing_option)
+        
+        # create the coordinate structs
+        r, r_spectral = define_coordinate(r_input,init_YY=false)
+        z, z_spectral = define_coordinate(z_input,init_YY=false)
+        vperp, vperp_spectral = define_coordinate(vperp_input,init_YY=false)
+        vpa, vpa_spectral = define_coordinate(vpa_input,init_YY=false)
+        gyrophase, gyrophase_spectral = define_coordinate(gyrophase_input,init_YY=false)
+        
+        # create test geometry
+        #rhostar = 0.1 #rhostar of ions for ExB drift
+        option = "constant-helical"
+        #pitch = 1.0
+        DeltaB = 1.0
+        geometry_in = geometry_input(rhostar,option,pitch,DeltaB)
+        geometry = init_magnetic_geometry(geometry_in,z,r)
+        
+        # create test composition
+        composition = create_test_composition()
+            
+        # Set up MPI
+        initialize_comms!()
+        setup_distributed_memory_MPI(1,1,1,1)
+        looping.setup_loop_ranges!(block_rank[], block_size[];
+                                       s=composition.n_ion_species, sn=1,
+                                       r=r.n, z=z.n, vperp=vperp.n, vpa=vpa.n,
+                                       vzeta=1, vr=1, vz=1)
+                                       
+        # initialise the matrix for the gyroaverages
+        gyro = init_gyro_operators(vperp,z,r,gyrophase,geometry,composition)
+        # initialise a test field
+        phi = allocate_shared_float(z.n,r.n)
+        gphi = allocate_shared_float(vperp.n,z.n,r.n)
+        gphi_exact = allocate_float(vperp.n,z.n,r.n)
+        gphi_err = allocate_float(vperp.n,z.n,r.n)
+        begin_serial_region()
+        @serial_region begin 
+            fill_test_arrays!(phi,gphi_exact,vperp,z,r,geometry,kz,kr,phasez,phaser)
+        end
+        
+        # gyroaverage phi
+        gyroaverage_field!(gphi,phi,gyro,vperp,z,r)
+        
+        # compute errors
+        begin_serial_region()
+        @serial_region begin
+            @. gphi_err = abs(gphi - gphi_exact)
+            println("Test gyroaverage_field!()")
+            for ivperp in 1:vperp.n
+                println("ivperp: ",ivperp," max(abs(gphi_err)): ",maximum(gphi_err[ivperp,:,:])," max(abs(gphi)): ",maximum(gphi[ivperp,:,:]))
+            end
+            println("")
+            if plot_test_results
+                @views heatmap(r.grid, z.grid, phi[:,:], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic,
+                    windowsize = (360,240), margin = 15pt)
+                outfile = "phi_vs_r_z.pdf"
+                savefig(outfile)
+                println("Saved outfile: "*outfile)
+                for ivperp in 1:vperp.n
+                    @views heatmap(r.grid, z.grid, gphi[ivperp,:,:], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic,
+                        windowsize = (360,240), margin = 15pt)
+                    outfile = "gphi_ivperp_"*string(ivperp)*"_vs_r_z.pdf"
+                    savefig(outfile)
+                    println("Saved outfile: "*outfile)
+                end
+            end
+        end
+        
+        # repeat the test for a pdf
+        # initialise a test field
+        nvpa = 1
+        n_ion_species = composition.n_ion_species
+        pdf = allocate_shared_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
+        gpdf = allocate_shared_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
+        gpdf_exact = allocate_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
+        gpdf_err = allocate_float(nvpa,vperp.n,z.n,r.n,n_ion_species)
+        begin_serial_region()
+        @serial_region begin
+            fill_pdf_test_arrays!(pdf,gpdf_exact,vpa,vperp,z,r,composition,geometry,kz,kr,phasez,phaser)
+        end
+        
+        gyroaverage_pdf!(gpdf,pdf,gyro,vpa,vperp,z,r,composition)
+        # compute errors
+        begin_serial_region()
+        @serial_region begin
+            @. gpdf_err = abs(gpdf - gpdf_exact)
+            println("Test gyroaverage_pdf!()")
+            for is in 1:n_ion_species
+                for ivperp in 1:vperp.n
+                    for ivpa in 1:nvpa
+                        println("ivpa: ",ivpa," ivperp: ",ivperp," is: ",is," max(abs(gphi_err)): ",maximum(gpdf_err[ivpa,ivperp,:,:,is]),
+                         " max(abs(gpdf)): ",maximum(gpdf[ivpa,ivperp,:,:,is]))
+                    end
+                end
+            end
+            println("")
+        end
+        
+        finalize_comms!()
+end
+
+function create_test_composition()
+    electron_physics = boltzmann_electron_response
+    n_ion_species = 1
+    n_neutral_species = 0
+    n_species = n_ion_species + n_neutral_species
+    use_test_neutral_wall_pdf = false
+    # electron temperature over reference temperature
+    T_e = 1.0
+    # temperature at the entrance to the wall in terms of the electron temperature
+    T_wall = 1.0
+    # wall potential at z = 0
+    phi_wall = 0.0
+    # constant to test nonzero Er
+    Er_constant = 0.0
+    # ratio of the neutral particle mass to the ion particle mass
+    mn_over_mi = 1.0
+    # ratio of the electron particle mass to the ion particle mass
+    me_over_mi = 1.0/1836.0
+    # The ion flux reaching the wall that is recycled as neutrals is reduced by
+    # `recycling_fraction` to account for ions absorbed by the wall.
+    recycling_fraction = 1.0
+    gyrokinetic_ions = true
+    return composition = species_composition(n_species, n_ion_species, n_neutral_species,
+            electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant,
+            mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species))
+end
+
+function fill_test_arrays!(phi,gphi,vperp,z,r,geometry,kz,kr,phasez,phaser)
+   for ir in 1:r.n
+      for iz in 1:z.n
+         Bmag = geometry.Bmag[iz,ir] 
+         bzeta = geometry.bzeta[iz,ir] 
+         rhostar = geometry.rhostar 
+         # convert integer "wavenumbers" to actual wavenumbers 
+         kkr = 2.0*pi*kr/r.L
+         kkz = 2.0*pi*kz/z.L
+         kperp = sqrt(kkr^2 + (bzeta*kkz)^2)
+         
+         phi[iz,ir] = sin(r.grid[ir]*kkr + phaser)*sin(z.grid[iz]*kkz + phasez)
+         for ivperp in 1:vperp.n
+            krho = kperp*vperp.grid[ivperp]*rhostar/Bmag
+            # use that phi is a sum of Kronecker deltas in k space to write gphi
+            gphi[ivperp,iz,ir] = besselj0(krho)*phi[iz,ir]
+         end
+      end
+   end  
+   return nothing
+end
+
+function fill_pdf_test_arrays!(pdf,gpdf,vpa,vperp,z,r,composition,geometry,kz,kr,phasez,phaser)
+   for is in 1:composition.n_ion_species
+       for ir in 1:r.n
+          for iz in 1:z.n
+             Bmag = geometry.Bmag[iz,ir] 
+             bzeta = geometry.bzeta[iz,ir] 
+             rhostar = geometry.rhostar 
+             # convert integer "wavenumbers" to actual wavenumbers 
+             kkr = 2.0*pi*kr/r.L
+             kkz = 2.0*pi*kz/z.L
+             kperp = sqrt(kkr^2 + (bzeta*kkz)^2)
+             
+             for ivperp in 1:vperp.n
+                for ivpa in 1:vpa.n
+                    pdf[ivpa,ivperp,iz,ir,is] = sin(r.grid[ir]*kkr + phaser)*sin(z.grid[iz]*kkz + phasez)
+                    krho = kperp*vperp.grid[ivperp]*rhostar/Bmag
+                    # use that pdf is a sum of Kronecker deltas in k space to write gpdf
+                    gpdf[ivpa,ivperp,iz,ir,is] = besselj0(krho)*pdf[ivpa,ivperp,iz,ir,is]
+                end
+             end
+          end
+       end
+   end  
+   return nothing
+end
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    using Pkg
+    Pkg.activate(".")
+    # example function call with arguments for a successful test
+    #gyroaverage_test()
+    gyroaverage_test(rhostar=0.01,pitch=0.5,kr=2,kz=3,phaser=0.25*pi,phasez=0.5*pi,ngrid=9,nelement=8,ngrid_vperp=5,nelement_vperp=3,Lvperp=18.0,ngrid_gyrophase=100,r_bc="periodic",z_bc="periodic")
+end