module RangeClass use MathClass use ArrayClass implicit none real(real64), parameter :: PF_RANGE_INFTY = dble(1.0e+14) type :: Range_ logical :: rect_is_active = .false. real(real64) :: x_range(1:2) real(real64) :: y_range(1:2) real(real64) :: z_range(1:2) real(real64) :: t_range(1:2) logical :: x_substituted(1:2) logical :: y_substituted(1:2) logical :: z_substituted(1:2) logical :: t_substituted(1:2) ! for sphere logical :: sphere_is_active = .false. real(real64) :: sphere_center(1:3) real(real64) :: sphere_radius ! for cylinder logical :: cylinder_is_active = .false. real(real64) :: cylinder_p1(1:3) real(real64) :: cylinder_p2(1:3) real(real64) :: cylinder_radius ! for cone logical :: cone_is_active = .false. real(real64) :: cone_p1(1:3) real(real64) :: cone_p2(1:3) real(real64) :: cone_radius(1:2) ! for multiple type(Range_),allocatable :: and_ranges(:) type(Range_),allocatable :: or_ranges(:) contains procedure,public :: init => initRange procedure,public :: set => setRange procedure,public :: get => getRange procedure,public :: set_x => set_xRange procedure,public :: set_y => set_yRange procedure,public :: set_z => set_zRange procedure,public :: set_t => set_tRange procedure,public :: inside => insideRange procedure,public :: lx => lx_RangeClass procedure,public :: ly => ly_RangeClass procedure,public :: lz => lz_RangeClass !procedure,public :: getSubrange => getSubrangeRange !procedure,public :: getSubrangeIdx => getSubrangeIdxRange end type interface to_range module procedure :: to_range_int32, to_range_real64, to_range_real64_rect,& to_range_real64_sphere,to_range_real64_cylinder,to_range_real64_Cone end interface interface print module procedure :: printRange end interface print interface operator(.in.) module procedure :: in_detect_by_range_int32, in_detect_by_range_real64, & in_int32_int32vector,in_detect_by_range_xyz_real64 end interface interface operator(.and.) module procedure :: and_rect_real64 end interface interface operator(.or.) module procedure :: or_xyz_range_real64 end interface contains ! ############################################################### function in_int32_int32vector(intval, intlist) result(ret) integer(int32), intent(in) :: intval, intlist(:) integer(int32) :: i logical :: ret ret = .false. do i = 1, size(intlist) if (intlist(i) == intval) then ret = .true. return end if end do end function ! ############################################################### ! ############################################################### function in_detect_by_range_int32(intval, in_range) result(ret) integer(int32), intent(in) :: intval type(Range_), intent(in) :: in_range logical :: ret ret = (in_range%x_range(1) <= dble(intval)) .and. (dble(intval) <= in_range%x_range(2)) end function ! ############################################################### ! ############################################################### function in_detect_by_range_real64(real64val, in_range) result(ret) real(real64), intent(in) :: real64val type(Range_), intent(in) :: in_range logical :: ret ret = (in_range%x_range(1) <= real64val) .and. (real64val <= in_range%x_range(2)) end function ! ############################################################### ! ############################################################### recursive function in_detect_by_range_xyz_real64(real64val, in_range) result(ret) real(real64), intent(in) :: real64val(:) real(real64) :: x(1:3),y(1:3),tot_l,l,nvec(1:3),this_R,cone_r_limit,theta type(Range_), intent(in) :: in_range logical :: ret x(:) = 0.0d0 ! <Boolean :: and> if (allocated(in_range%and_ranges))then ret = (real64val .in. in_range%and_ranges(1) ) .and. (real64val .in. in_range%and_ranges(2)) return endif ! <Boolean :: or> if (allocated(in_range%or_ranges))then ret = (real64val .in. in_range%or_ranges(1) ) .or. (real64val .in. in_range%or_ranges(2)) return endif ! <SPHERE> if(in_range%sphere_is_active )then x(1:size(real64Val)) = real64Val(:) ret = (in_range%sphere_radius & >= sqrt(dot_product(in_range%sphere_center-x,in_range%sphere_center-x) ) ) return endif ! <CYLINDER> if(in_range%cylinder_is_active )then x(1:size(real64Val)) = real64Val(:) ! p1 p2 ! *-------------------------* ! ! * x ! tot_l = sqrt(dot_product(& in_range%cylinder_p2-in_range%cylinder_p1,& in_range%cylinder_p2-in_range%cylinder_p1) ) if (tot_l == 0.0d0)then if (dot_product(x-in_range%cylinder_p1,x-in_range%cylinder_p1)==0.0d0 )then ret = .True. else ret = .False. endif else nvec = (in_range%cylinder_p2-in_range%cylinder_p1)/tot_l l = dot_product(x-in_range%cylinder_p1,in_range%cylinder_p2-in_range%cylinder_p1)/tot_L if (l < 0.0d0)then ret = .false. elseif (l > tot_L)then ret = .false. else y = in_range%cylinder_p1(:) + l*nvec(:) this_R = sqrt(dot_product(x-y,x-y)) ret = (in_range%cylinder_radius >= this_R ) endif endif return endif !<CONE> if(in_range%cone_is_active )then x(1:size(real64Val)) = real64Val(:) ! p1 p2 ! *-------------------------* ! ! * x ! tot_l = sqrt(dot_product(& in_range%cone_p2-in_range%cone_p1,& in_range%cone_p2-in_range%cone_p1) ) if (tot_l == 0.0d0)then if (dot_product(x-in_range%cone_p1,x-in_range%cone_p1)==0.0d0 )then ret = .True. else ret = .False. endif else nvec = (in_range%cone_p2-in_range%cone_p1)/tot_l l = dot_product(x-in_range%cone_p1,in_range%cone_p2-in_range%cone_p1)/tot_L theta = l/tot_l cone_r_limit = (1.0d0-theta)*in_range%cone_radius(1) + (theta)*in_range%cone_radius(2) if (l < 0.0d0)then ret = .false. elseif (l > tot_L)then ret = .false. else y = in_range%cone_p1(:) + l*nvec(:) this_R = sqrt(dot_product(x-y,x-y)) ret = (cone_r_limit >= this_R ) endif endif return endif if(size(real64val)==2 )then ret = ((in_range%x_range(1) <= real64val(1)) .and. (real64val(1) <= in_range%x_range(2))) & .and. ((in_range%y_range(1) <= real64val(2)) .and. (real64val(2) <= in_range%y_range(2))) elseif(size(real64val)==3 )then ret = ((in_range%x_range(1) <= real64val(1)) .and. (real64val(1) <= in_range%x_range(2))) & .and. ((in_range%y_range(1) <= real64val(2)) .and. (real64val(2) <= in_range%y_range(2))) & .and. ((in_range%z_range(1) <= real64val(3)) .and. (real64val(3) <= in_range%z_range(2))) else ret = (in_range%x_range(1) <= real64val(1)) .and. (real64val(1) <= in_range%x_range(2)) endif end function ! ############################################################### function to_range_int32(from_to) result(ret_range) type(Range_) :: ret_range integer(int32), intent(in) :: from_to(1:2) call ret_range%init() ret_range%x_range(1:2) = dble(from_to) ret_range%x_substituted(:) = .True. end function function to_range_real64(from_to) result(ret_range) type(Range_) :: ret_range real(real64), intent(in) :: from_to(:) call ret_range%init() ret_range%x_range(1:2) = from_to ret_range%x_substituted(:) = .True. end function function to_range_real64_rect(x_min, y_min, z_min, x_max, y_max, z_max) result(ret_range) type(Range_) :: ret_range real(real64), optional, intent(in) :: x_min, y_min, z_min, x_max, y_max, z_max call ret_range%init() if (present(x_min)) then ret_range%x_range(1) = x_min ret_range%x_substituted(1) = .True. end if if (present(y_min)) then ret_range%y_range(1) = y_min ret_range%y_substituted(1) = .True. end if if (present(z_min)) then ret_range%z_range(1) = z_min ret_range%z_substituted(1) = .True. end if if (present(x_max)) then ret_range%x_range(2) = x_max ret_range%x_substituted(2) = .True. end if if (present(y_max)) then ret_range%y_range(2) = y_max ret_range%y_substituted(2) = .True. end if if (present(z_max)) then ret_range%z_range(2) = z_max ret_range%z_substituted(2) = .True. end if end function ! ######################################################### subroutine initRange(this, MaxRange) class(Range_), intent(inout) :: this real(real64), optional, intent(in) :: MaxRange this%rect_is_active = .true. if (present(MaxRange)) then this%x_range = [-MaxRange, MaxRange] this%y_range = [-MaxRange, MaxRange] this%z_range = [-MaxRange, MaxRange] this%t_range = [-MaxRange, MaxRange] else this%x_range = [-PF_RANGE_INFTY, PF_RANGE_INFTY] this%y_range = [-PF_RANGE_INFTY, PF_RANGE_INFTY] this%z_range = [-PF_RANGE_INFTY, PF_RANGE_INFTY] this%t_range = [-PF_RANGE_INFTY, PF_RANGE_INFTY] end if this%x_substituted = [.false., .false.] this%y_substituted = [.false., .false.] this%z_substituted = [.false., .false.] this%t_substituted = [.false., .false.] end subroutine ! ######################################################### subroutine setRange(this, x_min, x_max, y_min, y_max, z_min, z_max) class(Range_), intent(inout) :: this real(real64), optional, intent(in) :: x_min, x_max real(real64), optional, intent(in) :: y_min, y_max real(real64), optional, intent(in) :: z_min, z_max this%rect_is_active = .true. call this%set_x(x_min=x_min, x_max=x_max) call this%set_y(y_min=y_min, y_max=y_max) call this%set_z(z_min=z_min, z_max=z_max) end subroutine ! ######################################################### subroutine set_xRange(this, x_min, x_max) class(Range_), intent(inout) :: this real(real64), optional, intent(in) :: x_min, x_max this%rect_is_active = .true. if (present(x_min)) then this%x_range(1) = x_min this%x_substituted(1) = .True. end if if (present(x_max)) then this%x_range(2) = x_max this%x_substituted(2) = .True. end if end subroutine ! ######################################################### ! ######################################################### subroutine set_yRange(this, y_min, y_max) class(Range_), intent(inout) :: this real(real64), optional, intent(in) :: y_min, y_max this%rect_is_active = .true. if (present(y_min)) then this%y_range(1) = y_min this%y_substituted(1) = .True. end if if (present(y_max)) then this%y_range(2) = y_max this%y_substituted(2) = .True. end if end subroutine ! ######################################################### ! ######################################################### subroutine set_zRange(this, z_min, z_max) class(Range_), intent(inout) :: this real(real64), optional, intent(in) :: z_min, z_max this%rect_is_active = .true. if (present(z_min)) then this%z_range(1) = z_min this%z_substituted(1) = .True. end if if (present(z_max)) then this%z_range(2) = z_max this%z_substituted(2) = .True. end if end subroutine ! ######################################################### ! ######################################################### subroutine set_tRange(this, t_min, t_max) class(Range_), intent(inout) :: this real(real64), optional, intent(in) :: t_min, t_max if (present(t_min)) then this%t_range(1) = t_min this%t_substituted(1) = .True. end if if (present(t_max)) then this%t_range(2) = t_max this%t_substituted(2) = .True. end if end subroutine ! ######################################################### ! ######################################################### function getRange(this, range_type) result(min_and_max) class(Range_), intent(inout) :: this character(1), intent(in) :: range_type real(real64) :: min_and_max(2) if(this%rect_is_active)then if (range_type == "x" .or. range_type == "X") then min_and_max = this%x_range end if if (range_type == "y" .or. range_type == "Y") then min_and_max = this%y_range end if if (range_type == "z" .or. range_type == "Z") then min_and_max = this%z_range end if if (range_type == "t" .or. range_type == "T") then min_and_max = this%t_range end if endif end function ! ######################################################### pure function insideRange(this, point) result(inside_is_true) class(Range_), intent(in) :: this real(real64), intent(in) :: point(:) integer(int32) :: i logical :: inside_is_true if (this%rect_is_active)then inside_is_true = .true. do i = 1, size(point) if (i == 1) then if (this%x_range(1) >= point(i) .or. this%x_range(2) <= point(i)) then inside_is_true = .false. return end if elseif (i == 2) then if (this%y_range(1) >= point(i) .or. this%y_range(2) <= point(i)) then inside_is_true = .false. return end if elseif (i == 3) then if (this%z_range(1) >= point(i) .or. this%z_range(2) <= point(i)) then inside_is_true = .false. return end if elseif (i == 4) then if (this%t_range(1) >= point(i) .or. this%t_range(2) <= point(i)) then inside_is_true = .false. return end if end if end do endif end function ! ######################################################### subroutine printRange(range) type(Range_),intent(in) :: range if(range%rect_is_active)then print *, "[",range%x_range(1),",",range%x_range(2),"]" print *, "[",range%y_range(1),",",range%y_range(2),"]" print *, "[",range%z_range(1),",",range%z_range(2),"]" endif ! cylinder等も追記 end subroutine ! ######################################################### elemental function and_rect_real64(range1,range2) result(ret) type(Range_),intent(in) :: range1,range2 type(Range_) :: ret allocate(ret%and_ranges(1:2)) ret%and_ranges(1) = range1 ret%and_ranges(2) = range2 ! regacy code >>> if (range1%rect_is_active .and. range2%rect_is_active)then ret%x_range(1) = maxval([range1%x_range(1),range2%x_range(1)]) ret%x_range(2) = minval([range1%x_range(2),range2%x_range(2)]) ret%y_range(1) = maxval([range1%y_range(1),range2%y_range(1)]) ret%y_range(2) = minval([range1%y_range(2),range2%y_range(2)]) ret%z_range(1) = maxval([range1%z_range(1),range2%z_range(1)]) ret%z_range(2) = minval([range1%z_range(2),range2%z_range(2)]) if(ret%x_range(1) > ret%x_range(2))then ret%x_range(:) = 0.0d0 endif if(ret%y_range(1) > ret%y_range(2))then ret%y_range(:) = 0.0d0 endif if(ret%z_range(1) > ret%z_range(2))then ret%z_range(:) = 0.0d0 endif endif end function ! ######################################################### ! ######################################################### elemental function or_xyz_range_real64(range1,range2) result(ret) type(Range_),intent(in) :: range1,range2 type(Range_) :: ret allocate(ret%or_ranges(1:2)) ret%or_ranges(1) = range1 ret%or_ranges(2) = range2 end function ! ######################################################### !function getSubrange(this,dim,n) result(ret) ! class(Range_),intent(in) :: this ! integer(int32),intent(in) :: dim ! number of dimension (e.g. x-y-z => 3) ! integer(int32),intent(in) :: n ! number of division (e.g. n=2 for binomial search) ! type(Range_),allocatable :: ret(:) ! real(real64) :: dr ! integer(int32) :: i,j,k,itr,n,idx ! ! if(dim==3)then ! allocate(ret(8)) ! idx = 0 ! do x_id=1,n ! do y_id=1,n ! do z_id=1,n ! idx = idx + 1 ! dr = (this%x_range(2)-this%x_range(1))/dble(n) ! ret(idx)%x_range(1) = this%x_range(1) + dr*dble(x_id-1) ! ret(idx)%x_range(2) = this%x_range(2) + dr*dble(x_id ) ! dr = (this%y_range(2)-this%y_range(1))/dble(n) ! ret(idx)%y_range(1) = this%y_range(1) + dr*dble(y_id-1) ! ret(idx)%y_range(2) = this%y_range(2) + dr*dble(y_id ) ! dr = (this%z_range(2)-this%z_range(1))/dble(n) ! ret(idx)%z_range(1) = this%z_range(1) + dr*dble(z_id-1) ! ret(idx)%z_range(2) = this%z_range(2) + dr*dble(z_id ) ! enddo ! enddo ! enddo ! return ! else ! ! not implemented yet. ! return ! endif ! ! ! !end function ! ###################################################### function lx_RangeClass(this) result(ret) class(Range_),intent(in) :: this real(real64) :: ret ret = this%x_range(2) - this%x_range(1) end function ! ###################################################### ! ###################################################### function ly_RangeClass(this) result(ret) class(Range_),intent(in) :: this real(real64) :: ret ret = this%y_range(2) - this%y_range(1) end function ! ###################################################### ! ###################################################### function lz_RangeClass(this) result(ret) class(Range_),intent(in) :: this real(real64) :: ret ret = this%z_range(2) - this%z_range(1) end function ! ###################################################### ! ###################################################### function to_range_real64_sphere(center,radius) result(ret) real(real64),intent(in) :: center(1:3),radius type(Range_) :: ret ret%sphere_is_active = .true. ret%sphere_center = center ret%sphere_radius = abs(radius) end function ! ###################################################### ! ###################################################### function to_range_real64_cylinder(p1,p2,radius) result(ret) real(real64),intent(in) :: p1(1:3),p2(1:3),radius type(Range_) :: ret ret%cylinder_is_active = .true. ret%cylinder_p1 = p1 ret%cylinder_p2 = p2 ret%cylinder_radius = abs(radius) end function ! ###################################################### ! ###################################################### function to_range_real64_Cone(p1,p2,radius) result(ret) real(real64),intent(in) :: p1(1:3),p2(1:3),radius(1:2) type(Range_) :: ret ret%cone_is_active = .true. ret%cone_p1 = p1 ret%cone_p2 = p2 ret%cone_radius = abs(radius) end function ! ###################################################### end module RangeClass