Main Program
!-----------------------------------------------------------------------
! HPF_RAYTRACE.F90
! Parallel Ray Tracing program using HPF language
! Author : Kayasidh Kasyapanun
! Course : CPS615 Fall96 at Syracuse University, Syracuse, NY, USA
!-----------------------------------------------------------------------
PROGRAM ray_tracing
USE hpf_raytrace_type
IMPLICIT NONE
INTEGER :: sub_img_x, sub_img_y ! subimage size to trace the ray
INTEGER :: num_img_x, num_img_y ! no of subimages in each axis
TYPE(ivector) :: img ! total image size to trace the ray
TYPE(vector) :: eye ! view point
TYPE(rgb_vector) :: light_color ! color of the light
TYPE(vector) :: light ! light source location
INTEGER :: num_of_object ! number of objects
TYPE(object),ALLOCATABLE,DIMENSION(:) :: obj ! objects
TYPE(rgb_vector),ALLOCATABLE,DIMENSION(:,:) :: ray_color
TYPE(vector),ALLOCATABLE,DIMENSION(:,:) :: ray_direction,ray_origin
INTEGER,ALLOCATABLE,DIMENSION (:,:) :: scr_x,scr_y
LOGICAL,ALLOCATABLE,DIMENSION (:,:) :: proc_active
CHARACTER,ALLOCATABLE,DIMENSION (:,:) :: raw_r,raw_g,raw_b
INTEGER :: i,j,l,depth,time1,time2,time_rate,stepx,stepy,ierr
REAL(8) :: xx,yy,zz,aa,bb,cc,dd,ee,ff,gg,hh,ii,jj
REAL(8) :: r1,g1,b1,r2,g2,b2,low_dot,low_thr
!--------------------------------------------------------------------------
! HPF directives
!--------------------------------------------------------------------------
!HPF$ DISTRIBUTE ray_color(CYCLIC,CYCLIC)
!HPF$ ALIGN (:,:) WITH ray_color(:,:) :: ray_direction,ray_origin
!HPF$ ALIGN (:,:) WITH ray_color(:,:) :: scr_x,scr_y,proc_active
!HPF$ ALIGN (:,:) WITH ray_color(:,:) :: raw_r,raw_g,raw_b
!--------------------------------------------------------------------------
! Read the ray tracing model parameters from file 'data.image'
!--------------------------------------------------------------------------
OPEN(1,FILE="data.image",status="old")
READ(1,*) sub_img_x,sub_img_y ! size of subimage
READ(1,*) img%x,img%y ! size of the entire image to be created
READ(1,*) img%z ! z-axis distance of the image from 0,0,0
ALLOCATE (ray_color(sub_img_y,sub_img_x),STAT=ierr)
IF (ierr .NE. 0) THEN
PRINT*, "ray_color: can't allocate"
STOP
ENDIF
ALLOCATE (ray_direction(sub_img_y,sub_img_x),STAT=ierr)
IF (ierr .NE. 0) THEN
PRINT*, "ray_direction: can't allocate"
STOP
ENDIF
ALLOCATE (ray_origin(sub_img_y,sub_img_x),STAT=ierr)
IF (ierr .NE. 0) THEN
PRINT*, "ray_origin: can't allocate"
STOP
ENDIF
allocate (scr_x(sub_img_y,sub_img_x),STAT=ierr)
IF (ierr .NE. 0) THEN
PRINT*, "scr_x: can't allocate"
STOP
ENDIF
ALLOCATE (scr_y(sub_img_y,sub_img_x),STAT=ierr)
IF (ierr .NE. 0) THEN
PRINT*, "scr_y: can't allocate"
STOP
ENDIF
ALLOCATE (proc_active(sub_img_y,sub_img_x),STAT=ierr)
IF (ierr .NE. 0) THEN
PRINT*, "proce_active: can't allocate"
STOP
ENDIF
ALLOCATE (raw_r(img%y,img%x),STAT=ierr) ! to store entire image
IF (ierr .NE. 0) THEN ! in 24bit rgb raw format
PRINT*, "raw_r: can't allocate"
STOP
ENDIF
ALLOCATE (raw_g(img%y,img%x),STAT=ierr)
IF (ierr .NE. 0) THEN
PRINT*, "raw_g : can't allocate"
STOP
ENDIF
ALLOCATE (raw_b(img%y,img%x),STAT=ierr)
IF (ierr .NE. 0) THEN
PRINT*, "raw_b: can't allocate"
STOP
ENDIF
READ(1,*) light%x,light%y,light%z
READ(1,*) light_color%r,light_color%g,light_color%b
READ(1,*) eye%x, eye%y, eye%z
READ(1,*) num_of_object
ALLOCATE (obj(num_of_object),STAT=ierr)
IF (ierr .NE. 0) THEN
PRINT*, "obj: can't allocate"
STOP
ENDIF
DO i=1,num_of_object
READ(1,*) aa,bb,cc,dd,ee,ff,gg,hh,ii,jj,r1,g1,b1,r2,g2,b2,xx,l
obj(i)%a=aa; obj(i)%b=bb; obj(i)%c=cc
obj(i)%d=dd; obj(i)%e=ee; obj(i)%f=ff
obj(i)%g=gg; obj(i)%h=hh; obj(i)%i=ii
obj(i)%j=jj; obj(i)%are=r1; obj(i)%agr=g1
obj(i)%abl=b1; obj(i)%dre=r2; obj(i)%dgr=g2
obj(i)%dbl=b2; obj(i)%s=xx; obj(i)%n=l
END DO
READ(1,*) depth,low_thr,low_dot
CLOSE(1)
PRINT *, "Finish reading input file"
!--------------------------------------------------------------------------
! Start the computation
!--------------------------------------------------------------------------
proc_active = .TRUE. ! initially, all processors are acive
ray_color = rgb_vector(0.0,0.0,0.0) ! init color to black
ray_origin%x = eye%x
ray_origin%y = eye%y
ray_origin%z = eye%z
print *, "Start"
CALL SYSTEM_CLOCK(time1)
! create coordinate (x,y) of each screen pixel on subimage plane
scr_x = SPREAD( (/(l,l=0,sub_img_x - 1)/), 2, &
& sub_img_y )
scr_y = SPREAD( (/(l,l=0,sub_img_y - 1)/), 1, &
& sub_img_x )
! loop doing every subimages and put it in a result raw-format array
DO i = 1, CEILING(REAL(img%y)/REAL(sub_img_y))
stepy = sub_img_y * (i-1)
DO j = 1, CEILING(real(img%x)/real(sub_img_x))
stepx = sub_img_x * (j-1)
! compute the (x,y) of each pixel on the subimage plane
ray_direction%x = scr_x - img%x/2 + stepx
ray_direction%y = img%y/2 - stepy - scr_y - 1.0
ray_direction%z = img%z
! compute a vector from eye to each pixel on the subimage plane
ray_direction%x = ray_direction%x - eye%x
ray_direction%y = ray_direction%y - eye%y
ray_direction%z = ray_direction%z - eye%z
CALL normalize(ray_direction)
! start the recursive ray tracing
CALL trace_ray(ray_origin,ray_direction,depth,proc_active,ray_color)
! store the subimage rgb at the right place in the whole image
raw_r( stepx+1:stepx+sub_img_x , stepy+1:stepy+sub_img_y ) = &
& char(int(ray_color%r*255.0))
raw_g( stepx+1:stepx+sub_img_x , stepy+1:stepy+sub_img_y ) = &
& char(int(ray_color%g*255.0))
raw_b( stepx+1:stepx+sub_img_x , stepy+1:stepy+sub_img_y ) = &
& char(int(ray_color%b*255.0))
END DO
END DO
CALL SYSTEM_CLOCK(time2,time_rate)
WRITE(*,*) "=============================================="
WRITE(*,*) "TIMING"
WRITE(*,*) "Computing : ",(real(time2-time1))/real(time_rate)," sec."
WRITE(*,*) "=============================================="
PRINT *, "Writing image to file ... please wait ..."
! write the whole image to the file 'pic.raw'
OPEN (1,FILE='pic.raw',STATUS='unknown')
DO i=1,img%y
DO j=1,img%x
WRITE(1,'(3A)',ADVANCE='no') raw_r(j,i),raw_g(j,i),raw_b(j,i)
END DO
END DO
CLOSE(1)
PRINT *, "Finish. Image saved in pic.raw file"
DEALLOCATE(ray_direction)
DEALLOCATE(ray_origin)
DEALLOCATE(ray_color)
DEALLOCATE(scr_x)
DEALLOCATE(scr_y)
DEALLOCATE(raw_r)
DEALLOCATE(raw_g)
DEALLOCATE(raw_b)
DEALLOCATE(proc_active)
DEALLOCATE(obj)
STOP
CONTAINS
!=========================================================================
! TRACE_RAY: Main subroutine of the recursive ray-tracing algorithm
! IN: ray_ori = starting point of the rays (x,y,z)
! ray_dir = unit vectors of ray direction (x,y,z) (normalized)
! r_depth = recursion depth
! active = processers that are active
! OUT: ray_col = pixels' color
!=========================================================================
RECURSIVE SUBROUTINE trace_ray(ray_ori,ray_dir,r_depth, &
& active,ray_col)
TYPE(vector),INTENT(IN),DIMENSION(:,:) :: ray_ori,ray_dir
TYPE(rgb_vector),INTENT(INOUT),DIMENSION(:,:) :: ray_col
INTEGER,INTENT(IN) :: r_depth ! recursion depth
LOGICAL,INTENT(IN),DIMENSION(:,:) :: active
REAL(8),DIMENSION(SIZE(active,1),SIZE(active,2)) :: t
TYPE(vector),DIMENSION(SIZE(active,1),SIZE(active,2)) :: R_vec, &
& intersect_point
TYPE(rgb_vector),DIMENSION(SIZE(active,1),SIZE(active,2)) :: &
& local_col,pixel_col,amb_col,dif_col,spe_col,ref_col
INTEGER,DIMENSION(SIZE(active,1),SIZE(active,2)) :: nearest_object
LOGICAL,DIMENSION(SIZE(active,1),SIZE(active,2)) :: active_ray
INTEGER :: i
!HPF$ ALIGN (:,:) WITH ray_ori(:,:) :: t,intersect_point,R_vec,local_col, &
!HPF$ & nearest_object,active_ray,amb_col,pixel_col,dif_col, &
!HPF$ & spe_col,ref_col
!HPF$ DISTRIBUTE *(CYCLIC,CYCLIC) :: ray_ori,ray_dir,ray_col,active
nearest_object = 0 ! object is numbered from 1 (0=not intersect)
pixel_col = rgb_vector(0.0,0.0,0.0)
local_col = rgb_vector(0.0,0.0,0.0)
amb_col = rgb_vector(0.0,0.0,0.0)
dif_col = rgb_vector(0.0,0.0,0.0)
spe_col = rgb_vector(0.0,0.0,0.0)
ref_col = rgb_vector(0.0,0.0,0.0)
t = 0.0 ! distance from view point to object
intersect_point = vector(0.0,0.0,0.0)
R_vec = vector(0.0,0.0,0.0)
active_ray = active ! processors that are active
! find where the ray intersect the object
CALL Find_intersect(ray_ori,ray_dir,t,nearest_object,active_ray)
! where the ray intersect nothing that ray is inactive
WHERE ( t == 0 )
active_ray = .FALSE.
END WHERE
! start finding the color of each pixel
WHERE ( t == 0.0 ) ! intersect nothing
pixel_col%r = 0.0; pixel_col%g = 0.0; pixel_col%b = 0.0
ELSEWHERE ! intersect something
! find location of the intersect point
intersect_point%x = ray_ori%x + (t * ray_dir%x)
intersect_point%y = ray_ori%y + (t * ray_dir%y)
intersect_point%z = ray_ori%z + (t * ray_dir%z)
END WHERE
! *** find ambient color component
CALL find_ambient(nearest_object,amb_col,active_ray)
! *** find diffuse color component
CALL find_diffuse(intersect_point,nearest_object,dif_col,active_ray);
! *** find specular color component
CALL find_specular(intersect_point,ray_ori,nearest_object &
& ,spe_col,active_ray)
! combine all color components into one color
WHERE (active_ray .EQV. .TRUE.)
pixel_col%r = (amb_col%r + dif_col%r + spe_col%r)/3.0
pixel_col%g = (amb_col%g + dif_col%g + spe_col%g)/3.0
pixel_col%b = (amb_col%b + dif_col%b + spe_col%b)/3.0
local_col = pixel_col
END WHERE
! continue the ray if more recursion is needed
IF (r_depth > 0) THEN
CALL find_reflect(intersect_point,ray_ori,nearest_object &
& ,R_vec,active_ray)
CALL normalize(R_vec)
CALL trace_ray(intersect_point,R_vec,r_depth-1,active_ray,ref_col)
END IF
! weight local and reflect color, otherwise reflecting colors are too much
IF (r_depth > 0) THEN
WHERE ((active_ray .EQV. .TRUE.).and.(ref_col%r <> 0.0).and. &
& (ref_col%g<>0.0).and.(ref_col%b<>0.0))
pixel_col%r = ((2*local_col%r) + ref_col%r) / 3.0
pixel_col%g = ((2*local_col%g) + ref_col%g) / 3.0
pixel_col%b = ((2*local_col%b) + ref_col%b) / 3.0
END WHERE
END IF
ray_col = pixel_col
END SUBROUTINE trace_ray
!=======================================================================
! Find_intersect : find if the ray intersects any object
! IN : ray_ori, origin of the ray
! : ray_dir, direction of the ray (unit vector of ray direction)
! : active_ray, active ray or inactive ray
! OUT : distance, how far is the intersection point from the origin
! : obj_no, which object is intersected
!=======================================================================
SUBROUTINE Find_intersect(ray_ori,ray_dir,distance, &
& obj_no,active_ray)
type(vector),INTENT(IN),dimension(:,:) :: ray_ori,ray_dir
real(8),INTENT(OUT),dimension(:,:) :: distance
integer,INTENT(OUT),dimension(:,:) :: obj_no
logical,INTENT(IN),DIMENSION(:,:) :: active_ray
real(8),dimension(size(ray_ori,1),size(ray_ori,2)) :: min_t
integer :: i
!HPF$ ALIGN (:,:) WITH ray_ori(:,:) :: min_t
!HPF$ DISTRIBUTE *(CYCLIC,CYCLIC) :: ray_ori,ray_dir,distance, &
!HPF$ & obj_no,active_ray
min_t = 999999.0
! for each object, find if the ray intersect it or not
DO i=1,num_of_object
CALL Find_intersect_Quad(ray_ori,ray_dir,i,distance,active_ray)
WHERE ((distance0.0)) ! closest in front
min_t = distance ! the nearest distance
obj_no = i ! that nearest object
END WHERE
END DO
distance = min_t
WHERE (min_t == 999999.0) ! intersect nothing
distance = 0.0
obj_no = 0
END WHERE
END SUBROUTINE Find_intersect
!=======================================================================
!
!=======================================================================
SUBROUTINE Find_intersect_Quad(ray_ori,ray_dir,obj_no, &
& distance, active_ray)
type(vector),INTENT(IN),dimension(:,:) :: ray_ori,ray_dir
real(8),INTENT(OUT),dimension(:,:) :: distance
integer,INTENT(IN) :: obj_no
logical,INTENT(IN),dimension(:,:) :: active_ray
REAL(8),dimension(size(ray_ori,1),size(ray_ori,2)) :: a,b,c,d, &
& e,f,g,h,i,j,aq,bq,cq,dq,t0,t1
!HPF$ ALIGN (:,:) WITH ray_ori(:,:) :: a,b,c,d,e,f,g,h,i,j,aq,bq,cq,dq, &
!HPF$ & t0,t1
!HPF$ DISTRIBUTE *(CYCLIC,CYCLIC) :: ray_ori,ray_dir,distance,active_ray
t0 = 0.0; t1 = 0.0; aq = 0.0
bq = 0.0; cq = 0.0; dq = 0.0
a = obj(obj_no)%a
b = obj(obj_no)%b
c = obj(obj_no)%c
d = obj(obj_no)%d
e = obj(obj_no)%e
f = obj(obj_no)%f
g = obj(obj_no)%g
h = obj(obj_no)%h
i = obj(obj_no)%i
j = obj(obj_no)%j
where (active_ray .EQV. .TRUE.)
aq = a * (ray_dir%x**2) + (2 * b * ray_dir%x * ray_dir%y) &
& + (2 * c * ray_dir%x * ray_dir%z) + (e * (ray_dir%y**2)) &
& + (2 * f * ray_dir%y * ray_dir%z) + (h * (ray_dir%z**2))
bq = a * (ray_ori%x * ray_dir%x) + (b * ((ray_ori%x * &
& ray_dir%y) + (ray_dir%x * ray_ori%y))) + (c * &
& ((ray_ori%x * ray_dir%z)+(ray_dir%x * ray_ori%z))) &
& + (d*ray_dir%x) + (e * ray_ori%y * ray_dir%y) + &
& (f*((ray_ori%y*ray_dir%z)+(ray_dir%y*ray_ori%z))) &
& + (g * ray_dir%y) + (h * ray_ori%z * ray_dir%z) &
& + (i * ray_dir%z)
bq = 2.0 * bq
cq = a * (ray_ori%x**2) + (2*b*ray_ori%x*ray_ori%y) &
& + (2*c*ray_ori%x*ray_ori%z) + (2*d*ray_ori%x) &
& + (e*(ray_ori%y**2)) + (2*f*ray_ori%y*ray_ori%z) &
& + (2*g*ray_ori%y) + (h*(ray_ori%z**2)) + (2*i* &
& ray_ori%z) + j
end where
where ((aq <> 0.0) .and. (active_ray .EQV. .TRUE.))
dq = (bq**2) - (4*aq*cq)
end where
where ((aq == 0.0) .and. (active_ray .EQV. .TRUE.))
t0 = -cq/bq
end where
where ((aq <> 0.0) .and. (dq < 0.0) .and. (active_ray .EQV. .TRUE.))
distance = 0.0
end where
where ((aq <> 0.0) .and. (dq >= 0.0) .and. (active_ray .EQV. .TRUE.))
aq = 2*aq; dq = sqrt(dq)
t0 = ((-bq-dq)/aq)
end where
where ((aq == 0.0).and.(t0 < low_thr).and.(active_ray .EQV. .TRUE.))
distance = 0.0
end where
where ((aq == 0.0).and.(t0 >= low_thr).and.(active_ray .EQV. .TRUE.))
distance = t0
end where
where ((aq <> 0.0) .and. (dq >= 0.0) .and. (t0 >= low_thr) .and. &
& (active_ray .EQV. .TRUE.))
distance = t0
end where
where ((aq <> 0.0) .and. (dq >= 0.0) .and. (t0 < low_thr) .and. &
& (active_ray .EQV. .TRUE.))
t1 = (-bq+dq)/aq
end where
WHERE ((aq <> 0.0) .AND. (dq >= 0.0) .AND. (t0=low_thr) .AND. (active_ray .EQV. .TRUE.))
distance = t0
END WHERE
WHERE ((aq <> 0.0).AND.(dq >= 0.0).AND.(t00.0))
dif_col%r = light_color%r * LdotN
dif_col%g = light_color%g * LdotN
dif_col%b = light_color%b * LdotN
elsewhere
dif_col%r = 0.0; dif_col%g = 0.0; dif_col%b = 0.0
end where
do i=1,num_of_object
where ((intersect == 0).and.(i==obj_no))
dif_col%r = dif_col%r * obj(i)%dre
dif_col%g = dif_col%g * obj(i)%dgr
dif_col%b = dif_col%b * obj(i)%dbl
end where
end do
END SUBROUTINE find_diffuse
!=======================================================================
SUBROUTINE intersect_object(point,L_vec,intersect,active_ray)
TYPE(vector),intent(in),dimension(:,:) :: point,L_vec
INTEGER,intent(out),dimension(:,:) :: intersect
LOGICAL,intent(in),dimension(:,:) :: active_ray
INTEGER :: i
REAL(8),DIMENSION(size(point,1),size(point,2)) :: intersect_t
!HPF$ ALIGN (:,:) WITH point(:,:) :: intersect_t
!HPF$ DISTRIBUTE *(CYCLIC,CYCLIC) :: point,L_vec,intersect,active_ray
DO i=1,num_of_object
CALL find_intersect_quad(point,L_vec,i,intersect_t,active_ray)
WHERE ((intersect_t <> 0).and.(active_ray .EQV. .TRUE.))
intersect = 1
END WHERE
END DO
END SUBROUTINE intersect_object
!=======================================================================
SUBROUTINE find_specular(point,ray_ori,obj_no,spe_col,active_ray)
INTEGER,INTENT(IN),DIMENSION(:,:) :: obj_no
TYPE(vector),INTENT(IN),DIMENSION(:,:) :: point,ray_ori
TYPE(rgb_vector),INTENT(OUT),DIMENSION(:,:) :: spe_col
LOGICAL,INTENT(IN),DIMENSION(:,:) :: active_ray
TYPE(vector),DIMENSION(size(point,1),size(point,2)) :: L_vec, &
& N_vec,R_vec,V_vec,light_vec
REAL(8),DIMENSION(size(point,1),size(point,2)) :: LdotN,VdotR
INTEGER,DIMENSION(size(point,1),size(point,2)) :: intersect
INTEGER :: i
!HPF$ ALIGN (:,:) with point(:,:) :: LdotN,VdotR,intersect,light_vec, &
!HPF$ & L_vec,N_vec,R_vec,V_vec
!HPF$ DISTRIBUTE *(CYCLIC,CYCLIC) :: obj_no,point,ray_ori, &
!HPF$ & spe_col,active_ray
L_vec = vector(0.0,0.0,0.0)
N_vec = vector(0.0,0.0,0.0)
R_vec = vector(0.0,0.0,0.0)
V_vec = vector(0.0,0.0,0.0)
light_vec = light ! location of the light source
LdotN = 0.0
VdotR = 0.0
intersect = 0
where (active_ray .EQV. .TRUE.)
L_vec%x = light%x - point%x
L_vec%y = light%y - point%y
L_vec%z = light%z - point%z
end where
CALL normalize(L_vec)
where (active_ray .EQV. .TRUE.)
V_vec%x = ray_ori%x - point%x
V_vec%y = ray_ori%y - point%y
V_vec%z = ray_ori%z - point%z
end where
CALL normalize(V_vec)
CALL find_reflect(point,light_vec,obj_no,R_vec,active_ray) !get R_vec
CALL normalize(R_vec)
CALL calc_Normal_vector(N_vec,point,obj_no,active_ray)
! do i=1,num_of_object
! where ((active_ray .EQV. .TRUE.).and.(i==obj_no))
! N_vec%x = 2 * (obj(i)%a * point%x + &
! & obj(i)%b * point%y + &
! & obj(i)%c * point%z + &
! & obj(i)%d)
! N_vec%y = 2 * (obj(i)%b * point%x + &
! & obj(i)%e * point%y + &
! & obj(i)%f * point%z + &
! & obj(i)%g)
! N_vec%z = 2 * (obj(i)%c * point%x + &
! & obj(i)%f * point%y + &
! & obj(i)%h * point%z + &
! & obj(i)%i)
! END WHERE
! END DO
CALL normalize(N_vec)
CALL vec_dot_product(L_vec,N_vec,LdotN)
WHERE ((active_ray .EQV. .TRUE.).AND.(LdotN < low_dot))
spe_col%r = 0.0; spe_col%g = 0.0; spe_col%b = 0.0
END WHERE
CALL intersect_object(point,L_vec,intersect,active_ray)
WHERE ((intersect == 1) .AND. (active_ray .EQV. .TRUE.))
spe_col%r = 0.0; spe_col%g = 0.0; spe_col%b = 0.0
END WHERE
CALL vec_dot_product(V_vec,R_vec,VdotR)
DO i=1,num_of_object
WHERE ((intersect == 0) .AND. (VdotR > 0.0) .AND. &
& (active_ray .EQV. .TRUE.) .AND. (i == obj_no))
spe_col%r = (( (VdotR**obj(i)%n) * &
& light_color%r ) * obj(i)%s)
spe_col%g = (( (VdotR**obj(i)%n) * &
& light_color%g ) * obj(i)%s)
spe_col%b = (( (VdotR**obj(i)%n) * &
& light_color%b ) * obj(i)%s)
END WHERE
WHERE ((intersect==0) .AND. (VdotR<=0.0) .AND. &
& (active_ray .EQV. .TRUE.) .AND. (i==obj_no))
spe_col%r = 0.0; spe_col%g = 0.0; spe_col%b = 0.0
END WHERE
END DO
END SUBROUTINE find_specular
!=======================================================================
SUBROUTINE find_reflect(point,ray_ori,obj_no,R_vec,active_ray)
INTEGER,INTENT(IN),DIMENSION(:,:) :: obj_no
TYPE(vector),INTENT(IN),DIMENSION(:,:) :: point,ray_ori
TYPE(vector),INTENT(OUT),DIMENSION(:,:) :: R_vec
LOGICAL,INTENT(IN),DIMENSION(:,:) :: active_ray
REAL(8),DIMENSION(SIZE(point,1),SIZE(point,2)) :: IdotN
TYPE(vector),DIMENSION(SIZE(point,1),SIZE(point,2)) :: I_vec,N_vec
INTEGER :: i
!HPF$ ALIGN (:,:) WITH point(:,:) :: IdotN,I_vec,N_vec
!HPF$ DISTRIBUTE *(CYCLIC,CYCLIC) :: point,ray_ori,obj_no,R_vec, &
!HPF$ & active_ray
I_vec = vector(0.0,0.0,0.0)
N_vec = vector(0.0,0.0,0.0)
IdotN = 0.0
WHERE (active_ray .EQV. .TRUE.)
I_vec%x = ray_ori%x - point%x
I_vec%y = ray_ori%y - point%y
I_vec%z = ray_ori%z - point%z
end where
CALL normalize(I_vec)
CALL calc_Normal_vector(N_vec,point,obj_no,active_ray)
! DO i=1,num_of_object
! WHERE ((active_ray .EQV. .TRUE.) .AND. (i==obj_no))
! N_vec%x = 2 * (obj(i)%a * point%x + &
! & obj(i)%b * point%y + &
! & obj(i)%c * point%z + &
! & obj(i)%d)
! N_vec%y = 2 * (obj(i)%b * point%x + &
! & obj(i)%e * point%y + &
! & obj(i)%f * point%z + &
! & obj(i)%g)
! N_vec%z = 2 * (obj(i)%c * point%x + &
! & obj(i)%f * point%y + &
! & obj(i)%h * point%z + &
! & obj(i)%i)
! END WHERE
! END DO
CALL normalize(N_vec)
CALL vec_dot_product(I_vec,N_vec,IdotN)
WHERE (active_ray .EQV. .TRUE.)
IdotN = 2 * IdotN
N_vec%x = N_vec%x * IdotN
N_vec%y = N_vec%y * IdotN
N_vec%z = N_vec%z * IdotN
R_vec%x = N_vec%x - I_vec%x
R_vec%y = N_vec%y - I_vec%y
R_vec%z = N_vec%z - I_vec%z
END WHERE
END SUBROUTINE find_reflect
!=======================================================================
SUBROUTINE calc_Normal_vector(N_vec,point,obj_no,active_ray)
TYPE(vector),DIMENSION(:,:),INTENT(IN) :: point
TYPE(vector),DIMENSION(:,:),INTENT(INOUT) :: N_vec
INTEGER,DIMENSION(:,:),INTENT(IN) :: obj_no
LOGICAL,DIMENSION(:,:),INTENT(IN) :: active_ray
INTEGER :: i
!HPF$ DISTRIBUTE *(CYCLIC,CYCLIC) :: point,N_vec,obj_no,active_ray
DO i=1,num_of_object
WHERE ((active_ray .EQV. .TRUE.) .AND. (i == obj_no))
N_vec%x = 2 * (obj(i)%a * point%x + &
& obj(i)%b * point%y + &
& obj(i)%c * point%z + &
& obj(i)%d)
N_vec%y = 2 * (obj(i)%b * point%x + &
& obj(i)%e * point%y + &
& obj(i)%f * point%z + &
& obj(i)%g)
N_vec%z = 2 * (obj(i)%c * point%x + &
& obj(i)%f * point%y + &
& obj(i)%h * point%z + &
& obj(i)%i)
END WHERE
END DO
END SUBROUTINE
!=======================================================================
END PROGRAM ray_tracing