Added angle distribution

src/matcha/distribution_m.f90

@@ -8,20 +8,21 @@ module distribution_m
 
   type distribution_t
     private
-    double precision, allocatable, dimension(:) :: vel_, cumulative_distribution_
+    double precision, allocatable, dimension(:) :: vel_,cumulative_distribution_,angle_,cumulative_angle_distribution_
   contains
     procedure :: cumulative_distribution
+    procedure :: cumulative_angle_distribution
     procedure :: velocities
   end type  
 
   interface distribution_t
 
-   pure module function construct(sample_distribution) result(distribution)
+   pure module function construct(sample_distribution,sample_angle_distribution) result(distribution)
       implicit none
-      double precision, intent(in) :: sample_distribution(:,:)
+      double precision, intent(in) :: sample_distribution(:,:),sample_angle_distribution(:,:)
       type(distribution_t) distribution
-    end function
-    
+    end function 
+
   end interface
 
   interface
@@ -30,13 +31,20 @@ pure module function cumulative_distribution(self) result(my_cumulative_distribu
       implicit none
       class(distribution_t), intent(in) :: self
       double precision, allocatable :: my_cumulative_distribution(:)
-    end function
+    end function cumulative_distribution
+
+    pure module function cumulative_angle_distribution(self) result(my_cumulative_angle_distribution)
+      implicit none
+      class(distribution_t), intent(in) :: self
+      double precision, allocatable :: my_cumulative_angle_distribution(:)
+    end function cumulative_angle_distribution
+
 
-    pure module function velocities(self, speeds, directions) result(my_velocities)
+    pure module function velocities(self, speeds, angles, directions) result(my_velocities)
       !! Return the t_cell_collection_t object's velocity vectors
       implicit none
       class(distribution_t), intent(in) :: self
-      double precision, intent(in) :: speeds(:,:), directions(:,:,:)
+      double precision, intent(in) :: speeds(:,:),angles(:,:),directions(:,:,:)
       double precision, allocatable :: my_velocities(:,:,:)
     end function velocities
 

src/matcha/distribution_s.F90

@@ -2,7 +2,7 @@
 ! Terms of use are as specified in LICENSE.tx
 submodule(distribution_m) distribution_s
   use intrinsic_array_m, only : intrinsic_array_t
-  use do_concurrent_m, only : do_concurrent_sampled_speeds, do_concurrent_my_velocities
+  use do_concurrent_m, only : do_concurrent_sample, do_concurrent_my_velocities
 
 #ifdef USE_CAFFEINE
    use caffeine_assert_m, only : assert
@@ -35,44 +35,137 @@ pure function monotonically_increasing(f) result(monotonic)
         "distribution_t: monotonically_increasing(distribution%cumulative_distribution_)", &
         intrinsic_array_t(distribution%cumulative_distribution_))
     end associate
+
+    call assert(all(sample_angle_distribution(:,2)>=0.D0), "distribution_t%construct: sample_distribution>=0.", &
+      intrinsic_array_t(sample_angle_distribution))
+
+    associate(nintervals => size(sample_angle_distribution,1))      
+      distribution%angle_ = [(sample_angle_distribution(i,1), i =1, nintervals)]  ! Assign speeds to each distribution bin         
+      distribution%cumulative_angle_distribution_ = [0.D0, [(sum(sample_angle_distribution(1:i,2)), i=1, nintervals)]]
+
+      call assert(monotonically_increasing(distribution%cumulative_angle_distribution_), &
+        "distribution_t: monotonically_increasing(distribution%cumulative_distribution_)", &
+        intrinsic_array_t(distribution%cumulative_angle_distribution_))
+    end associate
 
+
   end procedure construct
 
   module procedure cumulative_distribution
     call assert(allocated(self%cumulative_distribution_), &
       "distribution_t%cumulative_distribution: allocated(cumulative_distribution_)")
     my_cumulative_distribution = self%cumulative_distribution_
+  end procedure
+
+  module procedure cumulative_angle_distribution
+    call assert(allocated(self%cumulative_angle_distribution_), &
+      "distribution_t%cumulative_angle_distribution: allocated(cumulative_angle_distribution_)")
+    my_cumulative_angle_distribution = self%cumulative_angle_distribution_
   end procedure 
 
+
   module procedure velocities
 
-    double precision, allocatable :: sampled_speeds(:,:),  dir(:,:,:)
-    integer step
+    double precision, allocatable :: sampled_speeds(:,:), sampled_angles(:,:), dir(:,:,:)
+    integer cell,step
+    double precision dir_mag_
+
+    double precision pi,twopi
+    double precision a,b,c
+    double precision x1,y1,z1,x2,y2,z2,x3,y3,z3,vec1mag
+    double precision vxp,vyp,vzp,eps,vec3mag
+    double precision random_phi
 
     call assert(allocated(self%cumulative_distribution_), &
       "distribution_t%cumulative_distribution: allocated(cumulative_distribution_)")
     call assert(allocated(self%vel_), "distribution_t%cumulative_distribution: allocated(vel_)")
 
+    call assert(allocated(self%cumulative_angle_distribution_), &
+      "distribution_t%cumulative_angle_distribution: allocated(cumulative_angle_distribution_)")
+    call assert(allocated(self%angle_), "distribution_t%cumulative_angle_distribution: allocated(angle_)")    
+
      ! Sample from the distribution
-     call do_concurrent_sampled_speeds(speeds, self%vel_, self%cumulative_distribution(), sampled_speeds)
+     call do_concurrent_sample(speeds, self%vel_, self%cumulative_distribution(), sampled_speeds)
+     call do_concurrent_sample(angles, self%angle_, self%cumulative_angle_distribution(), sampled_angles)
 
-     associate(nsteps => size(speeds,2))
-
-       ! Create unit vectors
-       dir = directions(:,1:nsteps,:)
-
-       associate(dir_mag => sqrt(dir(:,:,1)**2 +dir(:,:,2)**2 + dir(:,:,3)**2))
-         associate(dir_mag_ => merge(dir_mag, epsilon(dir_mag), dir_mag/=0.))
-           dir(:,:,1) = dir(:,:,1)/dir_mag_
-           dir(:,:,2) = dir(:,:,2)/dir_mag_
-           dir(:,:,3) = dir(:,:,3)/dir_mag_
-         end associate
-       end associate
-
-       call do_concurrent_my_velocities(nsteps, dir, sampled_speeds, my_velocities)
+     associate(ncells => size(speeds,1), nsteps => size(speeds,2))
+
+!       ! Create unit vectors
+!       dir = directions(:,1:nsteps,:)
+!
+!       associate(dir_mag => sqrt(dir(:,:,1)**2 +dir(:,:,2)**2 + dir(:,:,3)**2))
+!         associate(dir_mag_ => merge(dir_mag, epsilon(dir_mag), dir_mag/=0.))
+!           dir(:,:,1) = dir(:,:,1)/dir_mag_
+!           dir(:,:,2) = dir(:,:,2)/dir_mag_
+!           dir(:,:,3) = dir(:,:,3)/dir_mag_
+!         end associate
+!       end associate
 
+!       call do_concurrent_my_velocities(nsteps, dir, sampled_speeds, my_velocities)
+
+       allocate(my_velocities(ncells,nsteps,3))
+
+       pi = acos(-1.d0)
+       twopi = 2.d0*pi
+       eps = 1.d-12
+
+!       do cell = 1,ncells
+       do concurrent(cell = 1:ncells)
+         do step = 1,nsteps
+            random_phi = 2.d0*pi*directions(cell,step,1)
+            vxp = sampled_speeds(cell,step)*sin(sampled_angles(cell,step))*cos(random_phi)
+            vyp = sampled_speeds(cell,step)*sin(sampled_angles(cell,step))*sin(random_phi)
+            vzp = sampled_speeds(cell,step)*cos(sampled_angles(cell,step))
+            if (step .eq. 1) then
+               my_velocities(cell,step,1) = vxp
+               my_velocities(cell,step,2) = vyp
+               my_velocities(cell,step,3) = vzp
+            else
+               x3 = my_velocities(cell,step-1,1)
+               y3 = my_velocities(cell,step-1,2)
+               z3 = my_velocities(cell,step-1,3)
+               vec3mag = sqrt(x3**2 + y3**2 + z3**2)+eps
+               x3 = x3/vec3mag
+               y3 = y3/vec3mag
+               z3 = z3/vec3mag
+
+               if (abs(z3) .gt. eps) then
+                  a = 0.d0
+                  b = 1.d0
+                  c = 0.d0
+               else
+                  a = 0.d0
+                  b = 0.d0
+                  c = 1.d0
+               end if
+
+               !curl(x3,y3,z3,a,b,c,x1,y1,z1)               
+               x1 = y3*c - z3*b
+               y1 = z3*a - x3*c
+               z1 = x3*b - y3*a
+               vec1mag = sqrt(x1**2 + y1**2 + z1**2) + eps
+               x1 = x1/vec1mag
+               y1 = y1/vec1mag
+               z1 = z1/vec1mag
+
+               ! curl(x3,y3,z3,x1,y1,z1,x2,y2,z2)
+               x2 = y3*z1 - z3*y1
+               y2 = z3*x1 - x3*z1
+               z2 = x3*y1 - y3*x1               
+
+
+               my_velocities(cell,step,1) = vxp*x1 + vyp*x2 + vzp*x3
+               my_velocities(cell,step,2) = vxp*y1 + vyp*y2 + vzp*y3
+               my_velocities(cell,step,3) = vxp*z1 + vyp*z2 + vzp*z3
+
+            end if
+
+         end do
+       end do
+
      end associate
 
   end procedure
 
 end submodule distribution_s
+

src/matcha/do_concurrent_m.f90

@@ -3,15 +3,15 @@ module do_concurrent_m
   use t_cell_collection_m, only : t_cell_collection_t, t_cell_collection_bind_C_t
   implicit none
   private
-  public :: do_concurrent_sampled_speeds, do_concurrent_my_velocities, do_concurrent_k, do_concurrent_speeds
+  public :: do_concurrent_sample, do_concurrent_my_velocities, do_concurrent_k, do_concurrent_speeds
   public :: do_concurrent_output_distribution
 
   interface
 
-    pure module subroutine do_concurrent_sampled_speeds(speeds, vel, cumulative_distribution, sampled_speeds) bind(C)
+    pure module subroutine do_concurrent_sample(rvalues, val, cumulative_distribution, sampled_values) bind(C)
       implicit none
-      real(c_double), intent(in) :: speeds(:,:), vel(:), cumulative_distribution(:)
-      real(c_double), intent(out), allocatable :: sampled_speeds(:,:)
+      real(c_double), intent(in) :: rvalues(:,:), val(:), cumulative_distribution(:)
+      real(c_double), intent(out), allocatable :: sampled_values(:,:)
     end subroutine
 
     pure module subroutine do_concurrent_my_velocities(nsteps, dir, sampled_speeds, my_velocities) bind(C)

src/matcha/do_concurrent_s.f90

@@ -4,15 +4,15 @@
 
 contains
 
-  module procedure do_concurrent_sampled_speeds
+  module procedure do_concurrent_sample
 
     integer cell, step
 
-    associate(ncells => size(speeds,1), nsteps => size(speeds,2))
-      allocate(sampled_speeds(ncells,nsteps))
+    associate(ncells => size(rvalues,1), nsteps => size(rvalues,2))
+      allocate(sampled_values(ncells,nsteps))
       do concurrent(cell = 1:ncells, step = 1:nsteps)
-        associate(k => findloc(speeds(cell,step) >= cumulative_distribution, value=.false., dim=1)-1)
-          sampled_speeds(cell,step) = vel(k)
+        associate(k => findloc(rvalues(cell,step) >= cumulative_distribution, value=.false., dim=1)-1)
+          sampled_values(cell,step) = val(k)
         end associate
       end do
     end associate

src/matcha/input_m.f90

@@ -8,7 +8,7 @@ module input_m
 
   type input_t
     private
-    integer :: num_cells_ = 6000, num_positions_ = 6000, num_dimensions_ = 3, num_intervals_ = 4
+    integer :: num_cells_ = 6000, num_positions_ = 6000, num_dimensions_ = 3, num_intervals_ = 4, num_angle_intervals_ = 4
     double precision :: time_step_ = 0.1D0
     !double precision, allocatable :: sample_distribution_(:,:)
     !allocate(sample_distribution_(num_intervals_,2))
@@ -17,8 +17,10 @@ module input_m
     procedure :: num_positions
     procedure :: num_dimensions
     procedure :: num_intervals
+    procedure :: num_angle_intervals
     procedure :: time_step
     procedure :: sample_distribution
+    procedure :: sample_angle_distribution
   end type
 
   interface
@@ -46,6 +48,12 @@ pure module function num_intervals(self) result(n)
       class(input_t), intent(in) :: self
       integer n
     end function 
+
+    pure module function num_angle_intervals(self) result(n)
+      implicit none
+      class(input_t), intent(in) :: self
+      integer n
+    end function 
 
     pure module function time_step(self) result(dt)
       implicit none
@@ -57,7 +65,14 @@ pure module function sample_distribution(self) result(empirical_distribution)
       implicit none
       class(input_t), intent(in) :: self
       double precision, allocatable :: empirical_distribution(:,:)
-    end function sample_distribution    
+    end function sample_distribution
+
+    pure module function sample_angle_distribution(self) result(empirical_angle_distribution)
+      implicit none
+      class(input_t), intent(in) :: self
+      double precision, allocatable :: empirical_angle_distribution(:,:)
+    end function sample_angle_distribution    
+
 
   end interface
 

src/matcha/input_s.f90

@@ -21,6 +21,10 @@
       n = self%num_intervals_
   end procedure
 
+  module procedure num_angle_intervals
+      n = self%num_angle_intervals_
+  end procedure   
+
   module procedure time_step
       dt = self%time_step_
   end procedure
@@ -64,5 +68,47 @@
      end do
 
   end procedure    
+
+  module procedure sample_angle_distribution
+     integer i,nintervals
+     double precision angle_lower,angle_upper
+     double precision range,pi,dangle,sumy
+     double precision angle_lower_bin
+     double precision angle_upper_bin
+     double precision angle_middle_bin
+     double precision pi
+     double precision, allocatable :: angles(:),probability(:)
+     nintervals = self%num_angle_intervals_
+
+     pi = acos(-1.d0)
+     angle_lower = 0.d0
+     angle_upper = pi
+     range = angle_upper - angle_lower
+     pi = acos(-1.d0)
+     dangle = range/dble(nintervals)
+     allocate(angles(nintervals),probability(nintervals))
+!    Create normal distribution     
+     sumy = 0.d0
+     do i = 1,nintervals
+        angle_lower_bin = angle_lower + dble(i-1)*dangle
+        angle_upper_bin = angle_lower + dble(i)*dangle
+        angle_middle_bin = 0.5d0*(angle_lower_bin + angle_upper_bin)
+        angles(i) = angle_middle_bin
+        probability(i) = 1.d0/dble(nintervals) ! Use uniform distribution
+        sumy = sumy + probability(i)
+     end do
+
+     do i = 1,nintervals
+        probability(i) = probability(i)/sumy
+     end do
+
+     allocate(empirical_angle_distribution(nintervals,2))
+
+     do i = 1,nintervals
+         empirical_angle_distribution(i,1) = angles(i)
+         empirical_angle_distribution(i,2) = probability(i)
+     end do
+
+  end procedure    
 
 end submodule input_s