diff --git a/README.md b/README.md
index 84d791e..017222d 100644
--- a/README.md
+++ b/README.md
@@ -47,3 +47,74 @@ plasma_type = 1
 ```
 
 For a better understanding of the varibles take a look at the [C. Fausser et al](https://www.sciencedirect.com/science/article/pii/S0920379612000853) paper.
+
+## C API for use with Fortran
+
+A C API is provided for linking of the plasma source routine to Fortran. This is particularly useful for compilation with MCNP. To compile a static library and a test program a build script is provided in the ```parametric-plasma-source/parametric_plasma_source/fortran_api``` folder. This can be run in the following manner:
+
+```[bash]
+cd parametric-plasma-source/parametric_plasma_source/fortran_api
+./build_lib.sh intel
+```
+for use with intel ifort and icpc compilers or
+```[bash]
+cd parametric-plasma-source/parametric_plasma_source/fortran_api
+./build_lib.sh gnu
+```
+for use with the gnu gfortran and g++ compilers.
+
+### Use with MCNP
+
+In order to use the library with MCNPv6.2 the ```plasma_source_module.F90``` and ```mcnp_pp.F90``` should be placed in the ```src``` folder. The ```source.F90``` provided with MCNP should then be modified to:
+
+```[fortran]
+subroutine source
+
+  ! .. Use Statements ..
+  use mcnp_interfaces_mod, only : expirx
+  use mcnp_debug
+  use pp_source_mk2_mod
+
+  implicit none
+
+  call parametric_plasma_2
+
+  return
+end subroutine source
+```
+The MCNP Makefile should also be updated to point to the library during linking. This can be done by adding the first line and modifying the second line in the Makefile:
+```
+PPLIB = -lplasmasource -L$(PLASMA_SOURCE)
+COMPILE_LINE=$(LD) $(OUT_EXE)$(EXEC) $(F_OBJS) $(C_OBJS) $(ALL_LDFLAGS) $(PLOTLIBS) $(LIB_DMMP) $(EXTRALIBS) \
+             $(PPLIB)
+```
+When compiling MCNP the PLASMA_SOURCE variable then needs to be set to the folder containing ```libplasmasource.a```. An example of the compliation line would be:
+
+```
+make build CONFIG='intel openmpi omp plot' PLASMA_SOURCE=/plasma/source/dir/
+```
+
+The source parameters are passed using the rdum and idum cards in the mcnp input file
+
+```[fortran]
+ion_density_pedistal           = rdum(1)
+ion_density_seperatrix         = rdum(2)
+ion_density_origin             = rdum(3)
+ion_temperature_pedistal       = rdum(4)
+ion_temperature_seperatrix     = rdum(5)
+ion_temperature_origin         = rdum(6)
+ion_density_peaking_factor     = rdum(7)
+ion_temperature_peaking_factor = rdum(8)
+ion_temperature_beta           = rdum(9)
+minor_radius                   = rdum(10)
+major_radius                   = rdum(11)
+pedistal_radius                = rdum(12)
+elongation                     = rdum(13)
+triangularity                  = rdum(14)
+min_toroidal_angle             = rdum(15)
+max_toroidal_angle             = rdum(16)
+        
+number_of_bins                 = idum(2)
+plasma_id                      = idum(3)
+```
+Note that idum(1) is intentionally left unused. This can be used for source selection if multiple user defined sources are to be compiled in the same executable. 
diff --git a/parametric_plasma_source/Plasma_source.h b/parametric_plasma_source/Plasma_source.h
new file mode 100644
index 0000000..4d412c6
--- /dev/null
+++ b/parametric_plasma_source/Plasma_source.h
@@ -0,0 +1,45 @@
+#ifdef __cplusplus // Are we compiling this with a C++ compiler ?
+extern "C" {
+    namespace plasma_source{
+        class PlasmaSource;
+    }
+    using namespace plasma_source;
+    typedef PlasmaSource PLASMASOURCE;
+#else
+    // From the C side, we use an opaque pointer.
+    typedef struct PLASMASOURCE PLASMASOURCE;
+#endif
+
+//
+// PLASMASOURCE* create_emptyplasmasource();
+
+// Large Constructor
+PLASMASOURCE* create_plasmasource(const double ion_density_ped, const double ion_density_sep,
+                                  const double ion_density_origin, const double ion_temp_ped,
+                                  const double ion_temp_sep, const double ion_temp_origin, 
+                                  const double pedistal_rad, const double ion_density_peak,
+                                  const double ion_temp_peak, const double ion_temp_beta,
+                                  const double minor_radius, const double major_radius,
+                                  const double elongation, const double triangularity,
+                                  const double shafranov, char* plasma_type_c,
+                                  const int plasma_id, const int number_of_bins,
+                                  const double min_toroidal_angle,
+                                  const double max_toroidal_angle );
+                                  
+// Destructor
+void delete_plasmasource(PLASMASOURCE* source);
+
+// Sample Function
+void Sample_Plasma_Source(PLASMASOURCE* source, 
+                          double (&random_numbers_c)[8],
+                          double &x,
+                          double &y,
+                          double &z,
+                          double &u,
+                          double &v,
+                          double &w,
+                          double &E);
+                          
+#ifdef __cplusplus
+}
+#endif 
\ No newline at end of file
diff --git a/parametric_plasma_source/fortran_api/build_lib.sh b/parametric_plasma_source/fortran_api/build_lib.sh
new file mode 100755
index 0000000..f5fe35f
--- /dev/null
+++ b/parametric_plasma_source/fortran_api/build_lib.sh
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+if [ $1 == "gnu" ] ; then
+  FC="gfortran"
+  CXX="g++"
+elif [ $1 == "intel" ] ; then
+  FC="ifort"
+  CXX="icpc"
+else 
+  echo "you must provide a compiler name - either gnu or intel"
+  exit
+fi
+echo Building fortran api plasma source library with $FC and $CXX
+
+rm -rf *.o libplasmasource.a *.mod testprog; $CXX -c ../plasma_source.cpp ../plasma_source_api.cpp -std=c++11; ar cr libplasmasource.a *.o; $FC -o testprog plasma_source_module.F90 testprog.F90 -lplasmasource -L./ -lstdc++
+
+
diff --git a/parametric_plasma_source/fortran_api/mcnp_pp.F90 b/parametric_plasma_source/fortran_api/mcnp_pp.F90
new file mode 100644
index 0000000..e5aa32c
--- /dev/null
+++ b/parametric_plasma_source/fortran_api/mcnp_pp.F90
@@ -0,0 +1,153 @@
+module pp_source_mk2_mod
+
+! Modules for use with the parametric plasma source mk2
+use cf_plasma_source
+use iso_c_binding
+
+! MCNP modules
+use mcnp_debug
+use pblcom, only   : pbl
+use mcnp_params, only   : dknd, i4knd
+use mcnp_random, only   : rang
+
+implicit none
+
+contains
+
+    subroutine parametric_plasma_2()
+    
+        type(plasma_source_f), save :: source
+        
+        real(c_double) :: ion_density_pedistal 
+        real(c_double) :: ion_density_seperatrix
+        real(c_double) :: ion_density_origin 
+        real(c_double) :: ion_temperature_pedistal
+        real(c_double) :: ion_temperature_seperatrix
+        real(c_double) :: ion_temperature_origin
+        real(c_double) :: ion_density_peaking_factor
+        real(c_double) :: ion_temperature_peaking_factor
+        real(c_double) :: ion_temperature_beta
+        real(c_double) :: minor_radius
+        real(c_double) :: major_radius
+        real(c_double) :: pedistal_radius
+        real(c_double) :: elongation
+        real(c_double) :: triangularity
+        real(c_double) :: shafranov_shift
+        character (len=100) :: plasma_name
+        integer(c_int) :: number_of_bins
+        integer(c_int) :: plasma_id
+        real(c_double) :: min_toroidal_angle
+        real(c_double) :: max_toroidal_angle
+        
+        logical, save :: bInit=.FALSE. !Start uninitialised
+        real(c_double), allocatable :: random_numbers(:)
+        integer(i4knd), parameter :: num_of_rands = 8
+        
+        real(c_double) :: xxx
+        real(c_double) :: yyy
+        real(c_double) :: zzz
+        real(c_double) :: uuu
+        real(c_double) :: vvv
+        real(c_double) :: www
+        real(c_double) :: EEE
+        
+        
+        
+        ! Set up parameters from the rdum and idum cards
+        plasma_name = "parametric_plasma_source"
+        ion_density_pedistal           = rdum(1)
+        ion_density_seperatrix         = rdum(2)
+        ion_density_origin             = rdum(3)
+        ion_temperature_pedistal       = rdum(4)
+        ion_temperature_seperatrix     = rdum(5)
+        ion_temperature_origin         = rdum(6)
+        ion_density_peaking_factor     = rdum(7)
+        ion_temperature_peaking_factor = rdum(8)
+        ion_temperature_beta           = rdum(9)
+        minor_radius                   = rdum(10)
+        major_radius                   = rdum(11)
+        pedistal_radius                = rdum(12)
+        elongation                     = rdum(13)
+        triangularity                  = rdum(14)
+        min_toroidal_angle             = rdum(15)
+        max_toroidal_angle             = rdum(16)
+        
+        number_of_bins                 = idum(2)
+        plasma_id                      = idum(3) ! 1 is default; //0 = L mode anything else H/A mode
+        
+        ! If this is the first time through init the source
+        if(.not. bInit) then
+            bInit = .True.
+            
+            source = plasma_source_f(ion_density_pedistal, ion_density_seperatrix,         &
+                                     ion_density_origin, ion_temperature_pedistal,         &
+                                     ion_temperature_seperatrix, ion_temperature_origin,   &
+                                     pedistal_radius, ion_density_peaking_factor,          &
+                                     ion_temperature_peaking_factor, ion_temperature_beta, &
+                                     minor_radius, major_radius,                           &
+                                     elongation, triangularity,                            &
+                                     shafranov_shift, plasma_name,                         &
+                                     plasma_id, number_of_bins,                            &
+                                     min_toroidal_angle,                                   &
+                                     max_toroidal_angle )
+            
+        end if
+        
+        ! Start here if this isn't the first time this routine has been called and sample source
+        
+        !source is a cell 
+        pbl%i%jsu=0
+        !source is neutrons    
+        pbl%i%ipt=1
+        !time independent   
+        pbl%r%tme=0
+        !initialise to weight 1.0
+        pbl%r%wgt=1.0_dknd
+        
+        ! Get random numbers
+        random_numbers = get_random_numbers(num_of_rands)
+        
+        ! init variables
+        xxx=0.0_c_double
+        yyy=0.0_c_double
+        zzz=0.0_c_double
+        uuu=0.0_c_double
+        vvv=0.0_c_double
+        www=0.0_c_double
+        EEE=0.0_c_double
+        
+        ! Call C++ sampling routine
+        call source%sample(random_numbers,xxx,yyy,zzz,uuu,vvv,www,EEE)
+        
+        ! Convert into cm
+        pbl%r%x = xxx * 100 
+        pbl%r%y = yyy * 100
+        pbl%r%z = zzz * 100
+        pbl%r%erg = EEE
+        
+        ! Setup the rest of the particle information by calling findlv
+        call findlv()
+    
+    end subroutine parametric_plasma_2
+    
+    !----------------------------------------------------------------------------------
+    ! Function to fill an array with a given quantity of random numbers
+    !----------------------------------------------------------------------------------
+    function get_random_numbers(quantity) result(random_nums)
+        implicit none
+        
+        integer(i4knd), intent(in) :: quantity
+        real(dknd), allocatable :: random_nums(:)
+        
+        integer(i4knd) :: i
+        
+        if(allocated(random_nums))deallocate(random_nums)
+        allocate(random_nums(quantity))
+        do i = 1, quantity
+            random_nums(i) = rang()
+        end do
+    
+    end function get_random_numbers
+
+
+end module pp_source_mk2_mod
\ No newline at end of file
diff --git a/parametric_plasma_source/fortran_api/plasma_source_module.F90 b/parametric_plasma_source/fortran_api/plasma_source_module.F90
new file mode 100644
index 0000000..28f9ce2
--- /dev/null
+++ b/parametric_plasma_source/fortran_api/plasma_source_module.F90
@@ -0,0 +1,190 @@
+module cf_plasma_source
+
+    use iso_c_binding
+
+    implicit none
+    
+    !private
+    public :: plasma_source_f
+
+    ! Interfaces to C routines    
+  interface
+
+    function create_plasmasource_c(ion_density_ped, ion_density_sep, &
+                                   ion_density_origin, ion_temp_ped, &
+                                   ion_temp_sep, ion_temp_origin,    &
+                                   pedistal_rad, ion_density_peak,   &
+                                   ion_temp_peak, ion_temp_beta,     &
+                                   minor_radius, major_radius,       &
+                                   elongation, triangularity,        &
+                                   shafranov, plasma_type,           &
+                                   plasma_id, number_of_bins,        &
+                                   min_toroidal_angle,               &
+                                   max_toroidal_angle ) result(plasmasource_c) bind(C, name="create_plasmasource")
+                                 
+                                 
+        use iso_c_binding
+        implicit none
+
+        type(c_ptr) :: plasmasource_c
+        
+        real(c_double), intent(in), value :: ion_density_ped
+        real(c_double), intent(in), value :: ion_density_sep
+        real(c_double), intent(in), value :: ion_density_origin
+        real(c_double), intent(in), value :: ion_temp_ped
+        real(c_double), intent(in), value :: ion_temp_sep
+        real(c_double), intent(in), value :: ion_temp_origin
+        real(c_double), intent(in), value :: pedistal_rad
+        real(c_double), intent(in), value :: ion_density_peak
+        real(c_double), intent(in), value :: ion_temp_peak
+        real(c_double), intent(in), value :: ion_temp_beta
+        real(c_double), intent(in), value :: minor_radius
+        real(c_double), intent(in), value :: major_radius
+        real(c_double), intent(in), value :: elongation
+        real(c_double), intent(in), value :: triangularity
+        real(c_double), intent(in), value :: shafranov
+        character(len=1, kind=c_char), intent(in) :: plasma_type(*)
+        integer(c_int), intent(in), value :: plasma_id
+        integer(c_int), intent(in), value :: number_of_bins
+        real(c_double), intent(in), value :: min_toroidal_angle
+        real(c_double), intent(in), value :: max_toroidal_angle
+        
+    end function create_plasmasource_c
+
+    subroutine delete_plasmasource_c(source) bind(C, name="delete_plasmasource")
+    
+        use iso_c_binding
+        implicit none
+    
+        type(c_ptr), value :: source
+    
+    end subroutine delete_plasmasource_c
+    
+    subroutine Sample_Plasma_Source_c(source, random_numbers, x, y, z, u, v, w, E) bind(C, name="Sample_Plasma_Source")
+
+        use iso_c_binding
+        implicit none
+        
+        type(c_ptr), intent(in), value :: source
+        real(c_double), dimension(8), intent(in) :: random_numbers
+        real(c_double), intent(out) :: x
+        real(c_double), intent(out) :: y
+        real(c_double), intent(out) :: z
+        real(c_double), intent(out) :: u
+        real(c_double), intent(out) :: v
+        real(c_double), intent(out) :: w
+        real(c_double), intent(out) :: E
+        
+    end subroutine Sample_Plasma_Source_c
+
+  end interface
+  
+    type plasma_source_f
+        type(c_ptr) :: ptr
+  
+        contains
+        
+#ifdef __GNUC__
+            procedure :: delete => delete_plasmasource_polymorph_f ! Destructor for gfortran
+#else
+            final :: delete_plasmasource_f ! Destructor
+#endif
+        
+            procedure :: sample => samplesource_f
+    end type
+    
+    interface plasma_source_f
+        procedure create_plasmasource_f
+    end interface
+   
+    
+    contains
+    
+        function create_plasmasource_f(ion_density_ped, ion_density_sep, &
+                                     ion_density_origin, ion_temp_ped, &
+                                     ion_temp_sep, ion_temp_origin,    &
+                                     pedistal_rad, ion_density_peak,   &
+                                     ion_temp_peak, ion_temp_beta,     &
+                                     minor_radius, major_radius,       &
+                                     elongation, triangularity,        &
+                                     shafranov, f_plasma_type,         &
+                                     plasma_id, number_of_bins,        &
+                                     min_toroidal_angle,               &
+                                     max_toroidal_angle ) result(plasmasource)
+        
+            type(plasma_source_f) :: plasmasource
+        
+            real(c_double), intent(in) :: ion_density_ped
+            real(c_double), intent(in) :: ion_density_sep
+            real(c_double), intent(in) :: ion_density_origin
+            real(c_double), intent(in) :: ion_temp_ped
+            real(c_double), intent(in) :: ion_temp_sep
+            real(c_double), intent(in) :: ion_temp_origin
+            real(c_double), intent(in) :: pedistal_rad
+            real(c_double), intent(in) :: ion_density_peak
+            real(c_double), intent(in) :: ion_temp_peak
+            real(c_double), intent(in) :: ion_temp_beta
+            real(c_double), intent(in) :: minor_radius
+            real(c_double), intent(in) :: major_radius
+            real(c_double), intent(in) :: elongation
+            real(c_double), intent(in) :: triangularity
+            real(c_double), intent(in) :: shafranov
+            character(len=*), intent(in) :: f_plasma_type
+            integer(c_int), intent(in) :: plasma_id
+            integer(c_int), intent(in) :: number_of_bins
+            real(c_double), intent(in) :: min_toroidal_angle
+            real(c_double), intent(in) :: max_toroidal_angle
+    
+            character(len=1, kind=C_CHAR) :: c_str(len_trim(f_plasma_type) + 1)
+            integer :: n, i
+            
+            do i = 1, len_trim(f_plasma_type)
+              c_str(i)=f_plasma_type(i:i)
+            end do
+            c_str(len_trim(f_plasma_type)+1) = C_NULL_CHAR
+            
+            plasmasource%ptr = create_plasmasource_c(ion_density_ped, ion_density_sep, &
+                                                     ion_density_origin, ion_temp_ped, &
+                                                     ion_temp_sep, ion_temp_origin,    &
+                                                     pedistal_rad, ion_density_peak,   &
+                                                     ion_temp_peak, ion_temp_beta,     &
+                                                     minor_radius, major_radius,       &
+                                                     elongation, triangularity,        &
+                                                     shafranov, c_str,                 &
+                                                     plasma_id, number_of_bins,        &
+                                                     min_toroidal_angle,               &
+                                                     max_toroidal_angle )              
+            
+    
+        end function create_plasmasource_f
+        
+        subroutine delete_plasmasource_f(this)
+            implicit none
+            type(plasma_source_f) :: this
+            call delete_plasmasource_c(this%ptr)
+        end subroutine delete_plasmasource_f
+
+        ! Bounds procedure needs to take a polymorphic (class) argument
+        subroutine delete_plasmasource_polymorph_f(this)
+            implicit none
+            class(plasma_source_f) :: this
+            call delete_plasmasource_c(this%ptr)
+        end subroutine delete_plasmasource_polymorph_f
+        
+        subroutine samplesource_f(this, random_numbers, x, y, z, u, v, w, E)
+            
+            class(plasma_source_f) :: this
+            real(c_double), dimension(8), intent(in) :: random_numbers
+            real(c_double), intent(out) :: x
+            real(c_double), intent(out) :: y
+            real(c_double), intent(out) :: z
+            real(c_double), intent(out) :: u
+            real(c_double), intent(out) :: v
+            real(c_double), intent(out) :: w
+            real(c_double), intent(out) :: E
+        
+            call Sample_Plasma_Source_c(this%ptr, random_numbers, x, y, z, u, v, w, E)
+        
+        end subroutine samplesource_f
+
+end module cf_plasma_source
diff --git a/parametric_plasma_source/fortran_api/testprog.F90 b/parametric_plasma_source/fortran_api/testprog.F90
new file mode 100644
index 0000000..4334a86
--- /dev/null
+++ b/parametric_plasma_source/fortran_api/testprog.F90
@@ -0,0 +1,74 @@
+program test
+
+use cf_plasma_source
+use iso_c_binding
+
+implicit none
+
+    type(plasma_source_f), save :: source
+    
+    real(c_double), parameter :: ion_density_pedistal = 1.09e+20
+    real(c_double), parameter :: ion_density_seperatrix = 3e+19
+    real(c_double), parameter :: ion_density_origin = 1.09e+20
+    real(c_double), parameter :: ion_temperature_pedistal = 6.09
+    real(c_double), parameter :: ion_temperature_seperatrix = 0.1
+    real(c_double), parameter :: ion_temperature_origin = 45.9
+    real(c_double), parameter :: ion_density_peaking_factor = 1
+    real(c_double), parameter :: ion_temperature_peaking_factor = 8.06
+    real(c_double), parameter :: ion_temperature_beta = 6.0
+    real(c_double), parameter :: minor_radius = 1.56
+    real(c_double), parameter :: major_radius = 2.5
+    real(c_double), parameter :: pedistal_radius = 0.8 * minor_radius
+    real(c_double), parameter :: elongation = 2.0
+    real(c_double), parameter :: triangularity = 0.55
+    real(c_double), parameter :: shafranov_shift = 0.0
+    character (len=*), parameter :: plasma_name="parametric_plasma_source"
+    integer(c_int), parameter :: number_of_bins  = 100
+    integer(c_int), parameter :: plasma_id = 1
+    real(c_double), parameter :: min_toroidal_angle = 0
+    real(c_double), parameter :: max_toroidal_angle = 360
+
+    real(c_double) :: x
+    real(c_double) :: y
+    real(c_double) :: z
+    real(c_double) :: u 
+    real(c_double) :: v
+    real(c_double) :: w
+    real(c_double) :: E
+    
+    real(c_double), dimension(8) :: rand_nums
+    
+    integer :: i
+    integer, parameter :: num_points = 1000
+    character (len=*), parameter :: filename="starting_points.o"
+    
+    write(*,*)"Initialising source"
+    ! Create source type
+    source = plasma_source_f(ion_density_pedistal, ion_density_seperatrix, &
+                           ion_density_origin, ion_temperature_pedistal, &
+                           ion_temperature_seperatrix, ion_temperature_origin,    &
+                           pedistal_radius, ion_density_peaking_factor,   &
+                           ion_temperature_peaking_factor, ion_temperature_beta,     &
+                           minor_radius, major_radius,       &
+                           elongation, triangularity,        &
+                           shafranov_shift, plasma_name,         &
+                           plasma_id, number_of_bins,        &
+                           min_toroidal_angle,               &
+                           max_toroidal_angle )
+
+
+    ! Open a file to output results to
+    open(unit=100, file=filename)
+    
+    write(*,*)"Sampling ", num_points," points and writing them to ", filename
+    ! Sample the source and return the starting variables
+    do i=1,1000
+      ! write(*,*)"calling random number"
+      call random_number(rand_nums)
+      ! write(*,*)rand_nums
+      call source%sample(rand_nums,x,y,z,u,v,w,E)
+      write(100,*)x,y,z,u,v,w,E
+    end do
+    close(100)
+
+end program test
diff --git a/parametric_plasma_source/plasma_source_api.cpp b/parametric_plasma_source/plasma_source_api.cpp
new file mode 100644
index 0000000..98fc9e3
--- /dev/null
+++ b/parametric_plasma_source/plasma_source_api.cpp
@@ -0,0 +1,64 @@
+#include <algorithm>
+#include <array>
+#include <iterator>
+#include <string>
+#include <iostream>
+#include <vector>
+
+#include "Plasma_source.h"
+#include "plasma_source.hpp"
+
+using namespace std;
+
+PLASMASOURCE* create_plasmasource(const double ion_density_ped, const double ion_density_sep,
+                                  const double ion_density_origin, const double ion_temp_ped,
+                                  const double ion_temp_sep, const double ion_temp_origin, 
+                                  const double pedistal_rad, const double ion_density_peak,
+                                  const double ion_temp_peak, const double ion_temp_beta,
+                                  const double minor_radius, const double major_radius,
+                                  const double elongation, const double triangularity,
+                                  const double shafranov, char* plasma_type_c,
+                                  const int plasma_id, const int number_of_bins,
+                                  const double min_toroidal_angle,
+                                  const double max_toroidal_angle){                               
+     
+    std::string plasma_type;
+    plasma_type = plasma_type_c;                                   
+
+    return new plasma_source::PlasmaSource(ion_density_ped, ion_density_sep,
+                            ion_density_origin, ion_temp_ped,
+                            ion_temp_sep, ion_temp_origin, 
+                            pedistal_rad, ion_density_peak,
+                            ion_temp_peak, ion_temp_beta,
+                            minor_radius, major_radius,
+                            elongation, triangularity,
+                            shafranov, plasma_type,
+                            plasma_id, number_of_bins,
+                            min_toroidal_angle,
+                            max_toroidal_angle);
+           
+}
+
+void delete_plasmasource(PLASMASOURCE* source){
+    
+    // delete source;
+    
+}
+
+void Sample_Plasma_Source(PLASMASOURCE* source, 
+                          double (&random_numbers_c)[8],
+                          double &x,
+                          double &y,
+                          double &z,
+                          double &u,
+                          double &v,
+                          double &w,
+                          double &E){
+  
+    std::array<double, 8> random_numbers;
+    std::copy_n(std::begin(random_numbers_c),8,std::begin(random_numbers));
+    
+    
+    source->SampleSource(random_numbers,x,y,z,u,v,w,E);
+    
+}