!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
!written by Xinyu Li(202335098@mail.sdu.edu.cn) and Qi Zhao(zhaoqi_326326@163.com)
SUBROUTINE anomalous_conformal
!This procedure is used to compute anomalous body conformal
    USE RES_MODEL_PARAMETER
    USE CONSTANTPARAMETERS
	use omp_lib
    IMPLICIT NONE
    CHARACTER(200) FinenameOfAnomalous
	CHARACTER(200) print_vert0,print_vert1,print_vert2,Complex_Inclusion_x,Complex_Inclusion_y,Complex_Inclusion_z,temp
    CHARACTER(300) :: line	!used to read one line of the ASCII STL file
    INTEGER(KIND=4)::i,j,k,t,ii,jj,kk,n,i1,i2,i3,j1,j2,j3,k1,k2,k3,l1,l2,l3,i0
    INTEGER(KIND=4):: ios_stl,i_face,iv_face,nv_tmp,idx_v,ip_stl
    INTEGER(KIND=4), DIMENSION(:,:), ALLOCATABLE :: tmp_face
    REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: tmp_vert,tmp_norm
    REAL(KIND=8):: vx,vy,vz
    REAL(KIND=8):: ex1x,ex1y,ex1z,ex2x,ex2y,ex2z,crossx,crossy,crossz

	INTEGER(KIND=4)::iz,jz,kz,ix,jx,kx,iy,jy,ky,iii,jjj,kkk
	INTEGER(KIND=4)::index_z,index_x,index_y,dex_z,dex_x,dex_y
	INTEGER(KIND=4)::num_z,num_x,num_y,dex_dd_z,dex_dd_x,dex_dd_y,dex_max_z,dex_max_x,dex_max_y,dex_d_z,dex_d_x,dex_d_y
	REAL(KIND=8)::dist_d_z,dist_d_x,dist_d_y,dist_z,dist_x,dist_y,threshold
	INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: XX_min,XX_max,YY_min,YY_max,ZZ_min,ZZ_max
	
	REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: coor_z_max,coor_z_min,coor_x_max,coor_x_min,coor_y_max,coor_y_min
    TYPE Triangular_Coordinates
        INTEGER(KIND=4)::point_number
		REAL(KIND=8):: Coord_X,Coord_Y,Coord_Z                   
    END TYPE Triangular_Coordinates
    TYPE face
        INTEGER(KIND=4)::node_face
        TYPE(Triangular_Coordinates)             			::node_point1
        TYPE(Triangular_Coordinates)             			::node_point2
        TYPE(Triangular_Coordinates)             			::node_point3
    END TYPE face
    TYPE(face), DIMENSION(:), ALLOCATABLE 					:: Triangular_face_element
    TYPE(Triangular_Coordinates), DIMENSION(:), ALLOCATABLE :: Vert
    
    print*,'Start conformal processing of the anomalous body'
    mmx = 0
    mmy = 0
    mmz = 0
	X_max=1
	X_min=NXB
	Y_max=1
	Y_min=NYB
	Z_max=1
	Z_min=NZB
	threshold=1.0e-8
    
	!===============================Read triangle face element information====================================
	IF(Logic_AnomalousDat)THEN
		!>Original .dat format:
		!!  Line 1 : "Number of Nodes and Elements:"
		!!  Line 2 : n_point (number of nodes)
		!!  Line 3 : n_face  (number of triangular elements)
		!!  Line 4 : "Nodes Coordinates:"
		!!  next n_point lines : label, Coord_X, Coord_Y, Coord_Z
		!!  then 2 title lines, then n_face lines : label, node1, node2, node3
		FinenameOfAnomalous='Complex_anomalous.dat'
		OPEN(20240506,FILE=FinenameOfAnomalous,STATUS='OLD')
		Read(20240506, *) temp
        Read(20240506, *) n_point !get total number of  Node
        Read(20240506, *) n_face !get total number of  Element
        Read(20240506, *) temp

        ALLOCATE(Triangular_face_element(n_face))
        ALLOCATE(Vert(n_point))
        DO j=1,n_point
            READ(20240506,*)Vert(j)%point_number,Vert(j)%Coord_X,Vert(j)%Coord_Y,Vert(j)%Coord_Z
        ENDDO
		Read(20240506,*) temp
        Read(20240506,*) temp
		DO i=1,n_face
            READ(20240506,*)Triangular_face_element(i)%node_face,Triangular_face_element(i)%node_point1%point_number,Triangular_face_element(i)%node_point2%point_number,Triangular_face_element(i)%node_point3%point_number
        ENDDO
        CLOSE(20240506)
	ELSEIF(Logic_AnomalousStl)THEN
		!>ASCII STL format:
		!!  solid <name>
		!!  facet normal nx ny nz
		!!    outer loop
		!!      vertex x y z
		!!      vertex x y z
		!!      vertex x y z
		!!    endloop
		!!  endfacet
		!!  ...
		!!  endsolid <name>
		!!In an STL file the vertices are written once per facet, so the duplicated
		!!vertices are merged into unique nodes before filling Vert/Triangular_face_element.
		!!The vertex order of each facet is also checked against the facet normal so that
		!!the normal direction convention is the same as the .dat format.
		FinenameOfAnomalous='Complex_anomalous.stl'
		OPEN(20240506,FILE=FinenameOfAnomalous,STATUS='OLD')
		!>First pass: count the number of facets.
		n_face=0
		DO
			READ(20240506,'(A)',IOSTAT=ios_stl) line
			IF(ios_stl/=0) EXIT
			!>convert the line into lower case for keyword matching
			DO ii=1,LEN_TRIM(line)
				IF(line(ii:ii)>='A'.AND.line(ii:ii)<='Z') line(ii:ii)=ACHAR(IACHAR(line(ii:ii))+32)
			ENDDO
			IF(INDEX(line,'facet')>0 .AND. INDEX(line,'endfacet')==0) n_face=n_face+1
		ENDDO
		IF(n_face==0)THEN
			WRITE(*,*)'Error: no facet is found in Complex_anomalous.stl!'
			STOP
		ENDIF
		REWIND(20240506)
		ALLOCATE(tmp_face(3,n_face),tmp_norm(3,n_face),tmp_vert(3,3*n_face))
		nv_tmp=0
		i_face=0
		iv_face=0
		DO
			READ(20240506,'(A)',IOSTAT=ios_stl) line
			IF(ios_stl/=0) EXIT
			DO ii=1,LEN_TRIM(line)
				IF(line(ii:ii)>='A'.AND.line(ii:ii)<='Z') line(ii:ii)=ACHAR(IACHAR(line(ii:ii))+32)
			ENDDO
			IF(INDEX(line,'facet')>0 .AND. INDEX(line,'endfacet')==0)THEN
				!>a new facet begins
				i_face=i_face+1
				iv_face=0
				IF(INDEX(line,'normal')>0)THEN
					READ(line(INDEX(line,'normal')+6:),*,IOSTAT=ios_stl) tmp_norm(1,i_face),tmp_norm(2,i_face),tmp_norm(3,i_face)
					IF(ios_stl/=0)THEN
						WRITE(*,*)'Error: failed to read the facet normal line in Complex_anomalous.stl!'
						STOP
					ENDIF
				ENDIF
			ELSEIF(INDEX(line,'endfacet')>0)THEN
				!>a facet is finished, check that it has exactly 3 vertices
				IF(iv_face/=3)THEN
					WRITE(*,*)'Error: a facet with',iv_face,'vertices (instead of 3) is found in Complex_anomalous.stl!'
					STOP
				ENDIF
			ELSEIF(INDEX(line,'vertex')>0)THEN
				iv_face=iv_face+1
				IF(iv_face>3)THEN
					WRITE(*,*)'Error: a facet with more than 3 vertices is found in Complex_anomalous.stl!'
					STOP
				ENDIF
				READ(line(INDEX(line,'vertex')+6:),*,IOSTAT=ios_stl) vx,vy,vz
				IF(ios_stl/=0)THEN
					WRITE(*,*)'Error: failed to read a vertex line in Complex_anomalous.stl!'
					STOP
				ENDIF
				!>merge the duplicated vertices
				idx_v=0
				DO ip_stl=1,nv_tmp
					IF(ABS(tmp_vert(1,ip_stl)-vx)<threshold.AND.ABS(tmp_vert(2,ip_stl)-vy)<threshold.AND.ABS(tmp_vert(3,ip_stl)-vz)<threshold)THEN
						idx_v=ip_stl
						EXIT
					ENDIF
				ENDDO
				IF(idx_v==0)THEN
					nv_tmp=nv_tmp+1
					tmp_vert(1,nv_tmp)=vx
					tmp_vert(2,nv_tmp)=vy
					tmp_vert(3,nv_tmp)=vz
					idx_v=nv_tmp
				ENDIF
				tmp_face(iv_face,i_face)=idx_v
			ENDIF
		ENDDO
		!>Correct the vertex order of each facet: compare the cross-product normal
		!!with the facet normal stored in the STL file, swap node2/node3 if they
		!!point in opposite directions, so that the normal direction convention
		!!is the same as in the .dat format.
		DO i=1,n_face
			ex1x=tmp_vert(1,tmp_face(2,i))-tmp_vert(1,tmp_face(1,i))
			ex1y=tmp_vert(2,tmp_face(2,i))-tmp_vert(2,tmp_face(1,i))
			ex1z=tmp_vert(3,tmp_face(2,i))-tmp_vert(3,tmp_face(1,i))
			ex2x=tmp_vert(1,tmp_face(3,i))-tmp_vert(1,tmp_face(1,i))
			ex2y=tmp_vert(2,tmp_face(3,i))-tmp_vert(2,tmp_face(1,i))
			ex2z=tmp_vert(3,tmp_face(3,i))-tmp_vert(3,tmp_face(1,i))
			crossx=ex1y*ex2z-ex1z*ex2y
			crossy=ex1z*ex2x-ex1x*ex2z
			crossz=ex1x*ex2y-ex1y*ex2x
			IF(crossx*tmp_norm(1,i)+crossy*tmp_norm(2,i)+crossz*tmp_norm(3,i)<0.0D0)THEN
				idx_v=tmp_face(2,i)
				tmp_face(2,i)=tmp_face(3,i)
				tmp_face(3,i)=idx_v
			ENDIF
		ENDDO
		n_point=nv_tmp
        ALLOCATE(Triangular_face_element(n_face))
        ALLOCATE(Vert(n_point))
		DO j=1,n_point
			Vert(j)%point_number=j
			Vert(j)%Coord_X=tmp_vert(1,j)
			Vert(j)%Coord_Y=tmp_vert(2,j)
			Vert(j)%Coord_Z=tmp_vert(3,j)
		ENDDO
		DO i=1,n_face
			Triangular_face_element(i)%node_face=i
			Triangular_face_element(i)%node_point1%point_number=tmp_face(1,i)
			Triangular_face_element(i)%node_point2%point_number=tmp_face(2,i)
			Triangular_face_element(i)%node_point3%point_number=tmp_face(3,i)
		ENDDO
        CLOSE(20240506)
		DEALLOCATE(tmp_face,tmp_norm,tmp_vert)
		WRITE(*,*)'Complex_anomalous.stl read: n_point=',n_point,' n_face=',n_face
	ELSE
		WRITE(*,*)'Error: neither Complex_anomalous.dat nor Complex_anomalous.stl exists, anomalous_conformal can not run!'
		RETURN
	ENDIF
	!=========================================================================================================	
	!=====================================The first range reduction===========================================
	!Find out the maximum and minimum values of the abnormal volume triangular mesh in the three ranges.
		max_coord_x=maxval(Vert(:)%Coord_X)
		min_coord_x=minval(Vert(:)%Coord_X)
		max_coord_y=maxval(Vert(:)%Coord_Y)
		min_coord_y=minval(Vert(:)%Coord_Y)
		max_coord_z=maxval(Vert(:)%Coord_Z)
		min_coord_z=minval(Vert(:)%Coord_Z)
		!The anomalous volume is delimited in the hexahedron.
		DO ii=1,NX
            IF(coordinates_x(ii)>min_coord_x)THEN
			    X_min=ii-1
			    EXIT
			ENDIF
        ENDDO
		DO ii=1,NX
            IF(coordinates_x(ii)>max_coord_x)THEN
			    X_max=ii
			    EXIT
			ENDIF
        ENDDO
        DO jj=1,NY
            IF(coordinates_y(jj)>min_coord_y)THEN
			    Y_min=jj-1
			    EXIT
			ENDIF
        ENDDO
		DO jj=1,NY
            IF(coordinates_y(jj)>max_coord_y)THEN
			    Y_max=jj
			    EXIT
			ENDIF
        ENDDO
		DO kk=1,NZ
            IF(coordinates_z(kk)>min_coord_z)THEN
			    Z_min=kk-1
			    EXIT
			ENDIF
        ENDDO
		DO kk=1,NZ
            IF(coordinates_z(kk)>max_coord_z)THEN
			    Z_max=kk
			    EXIT
			ENDIF
        ENDDO
    !=========================================================================================================
		ALLOCATE(orig_z(3,NXB*NYB),orig_y(3,NXB*NZB),orig_x(3,NXB*NZB))
		ALLOCATE(vert0(3,n_face),vert1(3,n_face),vert2(3,n_face),edge1(3,n_face),edge2(3,n_face))
        ALLOCATE(det_z(NXB*NYB),det_x(NYB*NZB),det_y(NXB*NZB))
		ALLOCATE(u_z(NXB*NYB),u_x(NYB*NZB),u_y(NXB*NZB))
		ALLOCATE(v_z(NXB*NYB),v_x(NYB*NZB),v_y(NXB*NZB))
		ALLOCATE(t_z(NXB*NYB),t_x(NYB*NZB),t_y(NXB*NZB))
		ALLOCATE(coor_z(NXB,NYB,n_face),coor_y(NXB,NZB,n_face),coor_x(NYB,NZB,n_face))
		ALLOCATE(pvec_z(3,NXB*NYB),pvec_y(3,NXB*NZB),pvec_x(3,NYB*NZB))
        ALLOCATE(tvec_z(3,NXB*NYB),tvec_y(3,NXB*NZB),tvec_x(3,NYB*NZB))
		ALLOCATE(ZZ_min(X_max),ZZ_max(X_max),XX_min(Y_max),XX_max(Y_max),YY_min(Z_max),YY_max(Z_max))
		ALLOCATE(coor_z_min(X_max),coor_z_max(X_max),coor_x_min(Y_max),coor_x_max(Y_max),coor_y_min(Z_max),coor_y_max(Z_max))
		ALLOCATE(mmz_per(NXB*NYB),mmy_per(NXB*NZB),mmx_per(NYB*NZB))
		ALLOCATE(crosspoint_ZZ(X_max,NXB*NYB),crosspoint_YY(X_max,NXB*NZB),crosspoint_XX(Y_max,NYB*NZB))
        ALLOCATE( Face_Triangle_NormVect(3,n_face))
        orig_z=0.0D0
		orig_y=0.0D0
		orig_x=0.0D0
		vert0=0.0D0
		vert1=0.0D0
		vert2=0.0D0
		edge1=0.0D0
		edge2=0.0D0
		det_z=0.0D0
		det_x=0.0D0
		det_y=0.0D0
		u_z=0.0D0
		u_x=0.0D0
		u_y=0.0D0
		v_z=0.0D0
		v_x=0.0D0
		v_y=0.0D0
		t_z=0.0D0
		t_x=0.0D0
		t_y=0.0D0
		pvec_z=0.0D0
		pvec_x=0.0D0
		pvec_y=0.0D0
		tvec_z=0.0D0
		tvec_x=0.0D0
		tvec_y=0.0D0
        mmz=0
		mmx=0
		mmy=0
		mmz_per=0
		mmx_per=0
		mmy_per=0
	dir_z   = [0.D0,0.D0,1.D0]
	dir_y   = [0.D0,1.D0,0.D0]
	dir_x   = [1.D0,0.D0,0.D0]
	!=====================================Möller-Trumbore algorithm===========================================	
	vert0(1,:) = Vert(Triangular_face_element(:)%node_point1%point_number)%Coord_X
    vert0(2,:) = Vert(Triangular_face_element(:)%node_point1%point_number)%Coord_Y
    vert0(3,:) = Vert(Triangular_face_element(:)%node_point1%point_number)%Coord_Z

    vert1(1,:) = Vert(Triangular_face_element(:)%node_point2%point_number)%Coord_X
    vert1(2,:) = Vert(Triangular_face_element(:)%node_point2%point_number)%Coord_Y
    vert1(3,:) = Vert(Triangular_face_element(:)%node_point2%point_number)%Coord_Z

    vert2(1,:) = Vert(Triangular_face_element(:)%node_point3%point_number)%Coord_X
    vert2(2,:) = Vert(Triangular_face_element(:)%node_point3%point_number)%Coord_Y
    vert2(3,:) = Vert(Triangular_face_element(:)%node_point3%point_number)%Coord_Z

    edge1 = vert1 - vert0
    edge2 = vert2 - vert0
    DO i0=1,n_face
        Face_Triangle_NormVect(1,i0)=edge1(2,i0) * edge2(3,i0)-edge1(3,i0) * edge2(2,i0)
        Face_Triangle_NormVect(2,i0)=edge1(3,i0) * edge2(1,i0)-edge1(1,i0) * edge2(3,i0)
        Face_Triangle_NormVect(3,i0)=edge1(1,i0) * edge2(2,i0)-edge1(2,i0) * edge2(1,i0)
    ENDDO 
	print*,'Ray tracing begins'
	!call OMP_SET_NUM_THREADS(16)
	!$OMP PARALLEL DO PRIVATE(i0,j1,i1,Rz,kk,kkk,verts_Dotmultp)			
    !get the intersaction of ray and z-face and save it into "crosspoint_z"	
    Do j1=Y_min,Y_max
        Do i1=X_min,X_max
            Rz=(j1-1)*NXB+i1
            DO i0=1,n_face
                orig_z(1,Rz)=coordinates_x(i1)
            	orig_z(2,Rz)=coordinates_y(j1)
        		orig_z(3,Rz)=coordinates_z(1)
    			tvec_z(1:3,Rz) = orig_z(1:3,Rz) - vert0(1:3,i0)
				pvec_z(1,Rz) = dir_z(2)*edge2(3,i0) - dir_z(3)*edge2(2,i0)
                pvec_z(2,Rz) = dir_z(3)*edge2(1,i0) - dir_z(1)*edge2(3,i0)
            	pvec_z(3,Rz) = dir_z(1)*edge2(2,i0) - dir_z(2)*edge2(1,i0)
        		det_z(Rz)=edge1(1,i0)*pvec_z(1,Rz)+edge1(2,i0)*pvec_z(2,Rz)+edge1(3,i0)*pvec_z(3,Rz)
                IF (abs(det_z(Rz)) < eps105) THEN
            		CYCLE
        		ENDIF
    			u_z(Rz) = (tvec_z(1,Rz)*pvec_z(1,Rz)+tvec_z(2,Rz)*pvec_z(2,Rz)+tvec_z(3,Rz)*pvec_z(3,Rz))/det_z(Rz)
			    IF (u_z(Rz) < 0.0 .or. u_z(Rz) > 1.0) THEN
                    CYCLE
                ENDIF 
            	pvec_z(1,Rz) = tvec_z(2,Rz)*edge1(3,i0) - tvec_z(3,Rz)*edge1(2,i0)
        		pvec_z(2,Rz) = tvec_z(3,Rz)*edge1(1,i0) - tvec_z(1,Rz)*edge1(3,i0)
    			pvec_z(3,Rz) = tvec_z(1,Rz)*edge1(2,i0) - tvec_z(2,Rz)*edge1(1,i0)
				v_z(Rz) = (dir_z(1)*pvec_z(1,Rz)+dir_z(2)*pvec_z(2,Rz)+dir_z(3)*pvec_z(3,Rz))/det_z(Rz)
                IF (v_z(Rz) < 0.0 .or. u_z(Rz) + v_z(Rz) > 1.0) THEN
            		CYCLE
        		ENDIF
                t_z(Rz)=(edge2(1,i0)*pvec_z(1,Rz)+edge2(2,i0)*pvec_z(2,Rz)+edge2(3,i0)*pvec_z(3,Rz))/det_z(Rz)
            	coor_z(i1,j1,i0)%Global_Coord%Coord_X = orig_z(1,Rz) + t_z(Rz) * dir_z(1)
        		coor_z(i1,j1,i0)%Global_Coord%Coord_Y = orig_z(2,Rz) + t_z(Rz) * dir_z(2)
    			coor_z(i1,j1,i0)%Global_Coord%Coord_Z = orig_z(3,Rz) + t_z(Rz) * dir_z(3)
				mmz_per(Rz) = mmz_per(Rz) + 1    !The number of z-direction intersections of each facet element is stored
            	crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_X = coor_z(i1,j1,i0)%Global_Coord%Coord_X
        		crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Y = coor_z(i1,j1,i0)%Global_Coord%Coord_Y
    			crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z = coor_z(i1,j1,i0)%Global_Coord%Coord_Z
				IF(mmz_per(Rz)>1)THEN
                    IF(ABS(crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z-crosspoint_ZZ((mmz_per(Rz)-1),Rz)%Global_Coord%Coord_Z)<eps105)THEN
                        mmz_per(Rz) = mmz_per(Rz) - 1
                        CYCLE
                    ENDIF
                ENDIF
                !*********************************Determine the intersection_Z attribute*******************************
            	verts_Dotmultp = dir_z(1) * Face_Triangle_NormVect(1,i0) + dir_z(2) * Face_Triangle_NormVect(2,i0) + dir_z(3) * Face_Triangle_NormVect(3,i0)
        		if(verts_Dotmultp > 0.D0) THEN
    			    crosspoint_ZZ(mmz_per(Rz),Rz)%Log_In=.TRUE.
				ELSE
                	crosspoint_ZZ(mmz_per(Rz),Rz)%Log_In=.FALSE.
                ENDIF
            ENDDO
            IF(mmz_per(Rz)>1)THEN
                DO kk=2,mmz_per(Rz)
                    DO kkk=1,kk-1
                        IF(crosspoint_ZZ(kk,Rz)%Global_Coord%Coord_Z<crosspoint_ZZ(kkk,Rz)%Global_Coord%Coord_Z)THEN
                            CALL SWAP(crosspoint_ZZ(kk,Rz), crosspoint_ZZ(kkk,Rz))
                        ENDIF
                    ENDDO
                ENDDO
            ENDIF
        ENDDO
    ENDDO
    !$OMP END PARALLEL DO
    !get the intersaction of ray and y-face and save it into "crosspoint_y"
    !=======================================================================================================
    !$OMP PARALLEL DO PRIVATE(i0,k2,i2,Ry,jj,jjj,verts_Dotmultp)
    Do i2=X_min,X_max
        Do k2=Z_min,Z_max
            Ry=(i2-1)*NZB+k2
            DO i0=1,n_face 
			    orig_y(1,Ry)=coordinates_x(i2)
				orig_y(2,Ry)=coordinates_y(1)
				orig_y(3,Ry)=coordinates_z(k2)
                tvec_y(1:3,Ry) = orig_y(1:3,Ry) - vert0(1:3,i0)
				pvec_y(1,Ry) = dir_y(2)*edge2(3,i0) - dir_y(3)*edge2(2,i0)
                pvec_y(2,Ry) = dir_y(3)*edge2(1,i0) - dir_y(1)*edge2(3,i0)
                pvec_y(3,Ry) = dir_y(1)*edge2(2,i0) - dir_y(2)*edge2(1,i0)
				det_y(Ry)=edge1(1,i0)*pvec_y(1,Ry)+edge1(2,i0)*pvec_y(2,Ry)+edge1(3,i0)*pvec_y(3,Ry)
				IF (abs(det_y(Ry)) < eps105) THEN
					CYCLE
                END IF
				u_y(Ry) = (tvec_y(1,Ry)*pvec_y(1,Ry)+tvec_y(2,Ry)*pvec_y(2,Ry)+tvec_y(3,Ry)*pvec_y(3,Ry))/det_y(Ry)
                IF (u_y(Ry) < 0.0 .or. u_y(Ry) > 1.0) THEN
					CYCLE
                END IF
				pvec_y(1,Ry) = tvec_y(2,Ry)*edge1(3,i0) - tvec_y(3,Ry)*edge1(2,i0)
                pvec_y(2,Ry) = tvec_y(3,Ry)*edge1(1,i0) - tvec_y(1,Ry)*edge1(3,i0)
                pvec_y(3,Ry) = tvec_y(1,Ry)*edge1(2,i0) - tvec_y(2,Ry)*edge1(1,i0)
				v_y(Ry) = (dir_y(1)*pvec_y(1,Ry)+dir_y(2)*pvec_y(2,Ry)+dir_y(3)*pvec_y(3,Ry))/det_y(Ry)
                IF (v_y(Ry) < 0.0 .or. u_y(Ry) + v_y(Ry) > 1.0) THEN
					CYCLE
                END IF
                t_y(Ry)=(edge2(1,i0)*pvec_y(1,Ry)+edge2(2,i0)*pvec_y(2,Ry)+edge2(3,i0)*pvec_y(3,Ry))/det_y(Ry)
                coor_y(i2,k2,i0)%Global_Coord%Coord_X = orig_y(1,Ry) + t_y(Ry) * dir_y(1)
                coor_y(i2,k2,i0)%Global_Coord%Coord_Y = orig_y(2,Ry) + t_y(Ry) * dir_y(2)
                coor_y(i2,k2,i0)%Global_Coord%Coord_Z = orig_y(3,Ry) + t_y(Ry) * dir_y(3)
                mmy_per(Ry) = mmy_per(Ry) + 1    !The number of y-direction intersections of each facet element is stored
                crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_X = coor_y(i2,k2,i0)%Global_Coord%Coord_X
                crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y = coor_y(i2,k2,i0)%Global_Coord%Coord_Y
                crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Z = coor_y(i2,k2,i0)%Global_Coord%Coord_Z
				IF(mmy_per(Ry)>1)THEN
                    IF(ABS(crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y-crosspoint_YY((mmy_per(Ry)-1),Ry)%Global_Coord%Coord_Y)<eps105)THEN
                        mmy_per(Ry) = mmy_per(Ry) - 1
                        CYCLE
                    ENDIF
                ENDIF
                !*********************************Determine the intersection_Y attribute*******************************
                verts_Dotmultp = dir_y(1) * Face_Triangle_NormVect(1,i0) + dir_y(2) * Face_Triangle_NormVect(2,i0) + dir_y(3) * Face_Triangle_NormVect(3,i0)
                if(verts_Dotmultp > 0.D0)THEN
            		crosspoint_YY(mmy_per(Ry),Ry)%Log_In=.TRUE.
                ELSE
    				crosspoint_YY(mmy_per(Ry),Ry)%Log_In=.FALSE.
                ENDIF
            ENDDO
            IF(mmy_per(Ry)>1)THEN
                DO jj=2,mmy_per(Ry)
                    DO jjj=1,jj-1
                        IF(crosspoint_YY(jj,Ry)%Global_Coord%Coord_Y<crosspoint_YY(jjj,Ry)%Global_Coord%Coord_Y)THEN
                            CALL SWAP(crosspoint_YY(jj,Ry), crosspoint_YY(jjj,Ry))
                        ENDIF
                    ENDDO
                ENDDO
            ENDIF
		ENDDO
    ENDDO
    !$OMP END PARALLEL DO
    !get the intersaction of ray and y-face and save it into "crosspoint_x"
	!=======================================================================================================
    !$OMP PARALLEL DO PRIVATE(i0,k3,j3,ii,iii,Rx,verts_Dotmultp)
    Do k3=Z_min,Z_max
        Do j3=Y_min,Y_max
            Rx=(k3-1)*NYB+j3
            DO i0=1,n_face
                orig_x(1,Rx)=coordinates_x(1)
				orig_x(2,Rx)=coordinates_y(j3)
				orig_x(3,Rx)=coordinates_z(k3)
				tvec_x(1:3,Rx) = orig_x(1:3,Rx) - vert0(1:3,i0)
				pvec_x(1,Rx) = dir_x(2)*edge2(3,i0) - dir_x(3)*edge2(2,i0)
                pvec_x(2,Rx) = dir_x(3)*edge2(1,i0) - dir_x(1)*edge2(3,i0)
                pvec_x(3,Rx) = dir_x(1)*edge2(2,i0) - dir_x(2)*edge2(1,i0)
				det_x(Rx)=edge1(1,i0)*pvec_x(1,Rx)+edge1(2,i0)*pvec_x(2,Rx)+edge1(3,i0)*pvec_x(3,Rx)
				IF (abs(det_x(Rx)) < eps105) THEN
					CYCLE
                END IF
				u_x(Rx) = (tvec_x(1,Rx)*pvec_x(1,Rx)+tvec_x(2,Rx)*pvec_x(2,Rx)+tvec_x(3,Rx)*pvec_x(3,Rx))/det_x(Rx)
                IF (u_x(Rx) < 0.0 .or. u_x(Rx) > 1.0) THEN
					CYCLE
                END IF
				pvec_x(1,Rx) = tvec_x(2,Rx)*edge1(3,i0) - tvec_x(3,Rx)*edge1(2,i0)
                pvec_x(2,Rx) = tvec_x(3,Rx)*edge1(1,i0) - tvec_x(1,Rx)*edge1(3,i0)
                pvec_x(3,Rx) = tvec_x(1,Rx)*edge1(2,i0) - tvec_x(2,Rx)*edge1(1,i0)
				v_x(Rx) = (dir_x(1)*pvec_x(1,Rx)+dir_x(2)*pvec_x(2,Rx)+dir_x(3)*pvec_x(3,Rx))/det_x(Rx)
                IF (v_x(Rx) < 0.0 .or. u_x(Rx) + v_x(Rx) > 1.0) THEN
					CYCLE
                END IF
                t_x(Rx)=(edge2(1,i0)*pvec_x(1,Rx)+edge2(2,i0)*pvec_x(2,Rx)+edge2(3,i0)*pvec_x(3,Rx))/det_x(Rx)
                coor_x(j3,k3,i0)%Global_Coord%Coord_X = orig_x(1,Rx) + t_x(Rx) * dir_x(1)
                coor_x(j3,k3,i0)%Global_Coord%Coord_Y = orig_x(2,Rx) + t_x(Rx) * dir_x(2)
                coor_x(j3,k3,i0)%Global_Coord%Coord_Z = orig_x(3,Rx) + t_x(Rx) * dir_x(3)
                mmx_per(Rx) = mmx_per(Rx) + 1    !The number of x-direction intersections of each facet element is stored
				crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X = coor_x(j3,k3,i0)%Global_Coord%Coord_X
                crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Y = coor_x(j3,k3,i0)%Global_Coord%Coord_Y
                crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Z = coor_x(j3,k3,i0)%Global_Coord%Coord_Z
				IF(mmx_per(Rx)>1)THEN
                    IF(ABS(crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X-crosspoint_XX((mmx_per(Rx)-1),Rx)%Global_Coord%Coord_X)<eps105)THEN
                        mmx_per(Rx) = mmx_per(Rx) - 1
                        CYCLE
                    ENDIF
                ENDIF
                !*********************************Determine the intersection_X attribute*******************************
                verts_Dotmultp = dir_x(1) * Face_Triangle_NormVect(1,i0) + dir_x(2) * Face_Triangle_NormVect(2,i0) + dir_x(3) * Face_Triangle_NormVect(3,i0)
                if(verts_Dotmultp > 0.D0)THEN
                    crosspoint_XX(mmx_per(Rx),Rx)%Log_In=.TRUE.
                ELSE
                    crosspoint_XX(mmx_per(Rx),Rx)%Log_In=.FALSE.
                ENDIF
            ENDDO
            IF(mmx_per(Rx)>1)THEN
                DO ii=2,mmx_per(Rx)
                    DO iii=1,mmx_per(Rx)-1
                        IF(crosspoint_XX(ii,Rx)%Global_Coord%Coord_X<crosspoint_XX(iii,Rx)%Global_Coord%Coord_X)THEN
                            CALL SWAP(crosspoint_XX(ii,Rx), crosspoint_XX(iii,Rx))
                        ENDIF
                    ENDDO
                ENDDO
            ENDIF
        ENDDO
    ENDDO
    !$OMP END PARALLEL DO
    !=======================================================================================================
	!=========================================================================================================
	!$OMP PARALLEL
    !$OMP DO PRIVATE(ii)			
	Do j1=Y_min,Y_max
		coor_x_min(j1)=minval(coor_x(j1,:,:)%Global_Coord%Coord_X)
	    coor_x_max(j1)=maxval(coor_x(j1,:,:)%Global_Coord%Coord_X)
		XX_min(j1) = -1
        XX_max(j1) = -1
		DO ii=1,NX
			IF(coordinates_x(ii)>=coor_x_min(j1).and.XX_min(j1) == -1)THEN
				XX_min(j1)=ii-1
			ENDIF
			IF(coordinates_x(ii)>=coor_x_max(j1).and.XX_max(j1) == -1)THEN
				XX_max(j1)=ii
				EXIT
			ENDIF
		ENDDO
	ENDDO	
	!$OMP END  DO 
	!'Calculate the projection of the anomalous volume on the y-plane'
	!$OMP DO PRIVATE(kk)	
	Do i2=X_min,X_max
		coor_z_min(i2)=minval(coor_z(i2,:,:)%Global_Coord%Coord_Z)
	    coor_z_max(i2)=maxval(coor_z(i2,:,:)%Global_Coord%Coord_Z)
		ZZ_min(i2) = -1
        ZZ_max(i2) = -1
		DO kk=1,NZ
			IF(coordinates_z(kk)>=coor_z_min(i2).and.ZZ_min(i2) == -1)THEN
				ZZ_min(i2)=kk-1
			ENDIF
			IF(coordinates_z(kk)>=coor_z_max(i2).and.ZZ_max(i2) == -1)THEN
				ZZ_max(i2)=kk
				EXIT
			ENDIF
		ENDDO
	ENDDO
	!$OMP END  DO 
	!'Calculate the projection of the anomalous volume on the x-plane'
	!$OMP DO PRIVATE(jj)
	DO k3=Z_min,Z_max
	    coor_y_min(k3)=minval(coor_y(:,k3,:)%Global_Coord%Coord_Y)
	    coor_y_max(k3)=maxval(coor_y(:,k3,:)%Global_Coord%Coord_Y)
		YY_min(k3) = -1
        YY_max(k3) = -1
		DO jj=1,NY
			IF(coordinates_y(jj)>=coor_y_min(k3).and.YY_min(k3) == -1)THEN
				YY_min(k3)=jj-1
			ENDIF
			IF(coordinates_y(jj)>=coor_y_max(k3).and.YY_max(k3) == -1)THEN
				YY_max(k3)=jj
				EXIT
			ENDIF
		ENDDO
	ENDDO
	!$OMP END DO 
	!$OMP END PARALLEL
	print*,'The projection calculation of the model in three directions is completed.'			
	!=========================================================================================================
	!=====================================Calculate the electric conductivity===========================================	
    !Calculate the electric conductivity of z-dection
    DO j=Y_min,Y_max
		DO i=XX_min(j),XX_max(j)
			Rz=(j-1)*NXB+i
			IF (mod(mmz_per(Rz),2)==0 .and. mmz_per(Rz)/=0) THEN  !The case of an even number of intersection points
                idx_start=0
                idx_end=0
                Logic_1=.false.
                Logic_2=.false.
                DO KK=1,mmz_per(Rz)
                    IF(crosspoint_ZZ(KK,Rz)%Log_In)THEN  !The intersection point is the entry point
                        idx_start=KK
                    ELSEIF(.NOT.crosspoint_ZZ(KK,Rz)%Log_In)THEN  !The intersection point is the exit point
                        idx_end=KK
                    ENDIF
                    IF((idx_start>0).AND.(idx_end>0))THEN   !There are both entry and exit points on the ray simultaneously
                        DO k=1,NZ
                            Logic_1=((coordinates_z(k)<=crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z).AND.&
                                    &(coordinates_z(k+1)>=crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z))
							Logic_2=((coordinates_z(k)<=crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z).AND.&
                                    &(coordinates_z(k+1)>=crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z))
							IF(Logic_1) KIdx_1=k
                            IF(Logic_2) KIdx_2=k
                            IF(Logic_1.and.Logic_2)	EXIT
                        ENDDO
                        IF(idx_start<idx_end) THEN  !The ray first penetrate the stratum and enter the anomalous body
							IF(KIdx_1==KIdx_2) THEN 
								LenRatio_CCSIGZ(i,j,KIdx_1) = (crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z-crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z)/Cdelz(KIdx_1)
                                CCSIGZ(i,j,KIdx_1) = tao_abnormal * LenRatio_CCSIGZ(i,j,KIdx_1)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGZ(i,j,KIdx_1))
                            ELSEIF(KIdx_2 > KIdx_1)THEN
                                LenRatio_CCSIGZ(i,j,KIdx_1) = (crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z-coordinates_z(KIdx_1))/Cdelz(KIdx_1)
                        		LenRatio_CCSIGZ(i,j,KIdx_2) = (crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z-coordinates_z(KIdx_2))/Cdelz(KIdx_2)
                                CCSIGZ(i,j,KIdx_1) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGZ(i,j,KIdx_1) + tao_abnormal * (1-LenRatio_CCSIGZ(i,j,KIdx_1))
                                CCSIGZ(i,j,(KIdx_1+1):(KIdx_2-1)) = tao_abnormal
                                CCSIGZ(i,j,KIdx_2) = tao_abnormal * LenRatio_CCSIGZ(i,j,KIdx_2)+TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGZ(i,j,KIdx_2))
                            ENDIF
                        ELSE  !The ray first emerges from the anomaly and penetrates the stratum
                            IF(KIdx_1==KIdx_2) THEN
                                LenRatio_CCSIGZ(i,j,KIdx_1) = (crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z-crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z)/Cdelz(KIdx_1)
                                CCSIGZ(i,j,KIdx_1) = tao_abnormal * LenRatio_CCSIGZ(i,j,KIdx_1)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGZ(i,j,KIdx_1))
                            ELSEIF(KIdx_2 < KIdx_1)THEN
                                LenRatio_CCSIGZ(i,j,KIdx_1) = (coordinates_z(KIdx_1+1)-crosspoint_ZZ(idx_start,Rz)%Global_Coord%Coord_Z)/Cdelz(KIdx_1)
                        		LenRatio_CCSIGZ(i,j,KIdx_2) = (crosspoint_ZZ(idx_end,Rz)%Global_Coord%Coord_Z-coordinates_z(KIdx_2))/Cdelz(KIdx_2)
                            ENDIF
                            CCSIGZ(i,j,KIdx_2) = LenRatio_CCSIGZ(i,j,KIdx_2)*tao_abnormal+(1-LenRatio_CCSIGZ(i,j,KIdx_2))*TAR_CONDUCTIVITY(2)
                            CCSIGZ(i,j,(KIdx_2+1):(KIdx_1-1)) = TAR_CONDUCTIVITY(2)
                            CCSIGZ(i,j,KIdx_1) = TAR_CONDUCTIVITY(2)*LenRatio_CCSIGZ(i,j,KIdx_1)+(1-LenRatio_CCSIGZ(i,j,KIdx_1))*tao_abnormal
                        ENDIF
                    ENDIF
                ENDDO
            ELSEIF(mod(mmz_per(Rz),2)==1)THEN
                print*,'z-dection!!!ERROR!!!ERROR!!!ERROR!!!'
				print*,i,j
            ENDIF
        ENDDO
    ENDDO
	print*,'The equivalent conductivity calculation in the z direction is completed!'
    !Calculate the electric conductivity of x-dection
	DO k=Z_min,Z_max
		DO j=YY_min(k),YY_max(k)
			Rx=(k-1)*NYB+j
			IF (mod(mmx_per(Rx),2)==0.and.mmx_per(Rx)/=0) THEN  !The case of an even number of intersection points
                idx_start=0
                idx_end=0
                Logic_1=.false.
                Logic_2=.false.
                DO II=1,mmx_per(Rx)
                    IF(crosspoint_XX(II,Rx)%Log_In)THEN  !The intersection point is the entry point
                        idx_start=II
                    ELSEIF(.NOT.crosspoint_XX(II,Rx)%Log_In)THEN  !The intersection point is the exit point
                        idx_end=II
                    ENDIF
                    IF((idx_start>0).AND.(idx_end>0))THEN   !There are both entry and exit points on the ray simultaneously
                        DO i=1,NX
                            Logic_1=((coordinates_x(i)<=crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X).AND.&
                                    &(coordinates_x(i+1)>=crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X))
							Logic_2=((coordinates_x(i)<=crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X).AND.&
                                    &(coordinates_x(i+1)>=crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X))
							IF(Logic_1) KIdx_1=i
                            IF(Logic_2) KIdx_2=i
                            IF(Logic_1.and.Logic_2)	EXIT
                        ENDDO
                        IF(idx_start<idx_end) THEN  !The ray first penetrate the stratum and enter the anomalous body
							IF(KIdx_1==KIdx_2) THEN
								LenRatio_CCSIGX(KIdx_1,j,k) = (crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X)/Cdelx(KIdx_1)
                                CCSIGX(KIdx_1,j,k) = tao_abnormal * LenRatio_CCSIGX(KIdx_1,j,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGX(KIdx_1,j,k))
                            ELSEIF(KIdx_2 > KIdx_1)THEN
                                LenRatio_CCSIGX(KIdx_1,j,k) = (crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_1))/Cdelx(KIdx_1)
                        		LenRatio_CCSIGX(KIdx_2,j,k) = (crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_2))/Cdelx(KIdx_2)
                                CCSIGX(KIdx_1,j,k) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGX(KIdx_1,j,k) + tao_abnormal * (1-LenRatio_CCSIGX(KIdx_1,j,k))
                                CCSIGX((KIdx_1+1):(KIdx_2-1),j,k) = tao_abnormal
                                CCSIGX(KIdx_2,j,k) = tao_abnormal * LenRatio_CCSIGX(KIdx_2,j,k)+TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGX(KIdx_2,j,k))
                            ENDIF
                        ELSE  !The ray first emerges from the anomaly and penetrates the stratum
                            IF(KIdx_1==KIdx_2) THEN
                                LenRatio_CCSIGX(KIdx_1,j,k) = (crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X-crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X)/Cdelx(KIdx_1)
                                CCSIGX(KIdx_1,j,k) = tao_abnormal * LenRatio_CCSIGX(KIdx_1,j,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGX(KIdx_1,j,k))
                            ELSEIF(KIdx_2 < KIdx_1)THEN
                                LenRatio_CCSIGX(KIdx_1,j,k) = (coordinates_x(KIdx_1+1)-crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X)/Cdelx(KIdx_1)
                        		LenRatio_CCSIGX(KIdx_2,j,k) = (crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_2))/Cdelx(KIdx_2)
                            ENDIF
                            CCSIGX(KIdx_2,j,k) = tao_abnormal*LenRatio_CCSIGX(KIdx_2,j,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGX(KIdx_2,j,k))
                            CCSIGX((KIdx_2+1):(KIdx_1-1),j,k) = TAR_CONDUCTIVITY(2)
                            CCSIGX(KIdx_1,j,k) = TAR_CONDUCTIVITY(2)*LenRatio_CCSIGX(KIdx_1,j,k)+tao_abnormal*(1-LenRatio_CCSIGX(KIdx_1,j,k))
                        ENDIF
                    ENDIF
                ENDDO
            ELSEIF(mod(mmx_per(Rx),2)==1)THEN
                print*,'x-dection!!!ERROR!!!ERROR!!!ERROR!!!'
				print*,j,k,coordinates_y(j),coordinates_z(k)
            ENDIF
        ENDDO
    ENDDO
	print*,'The equivalent conductivity calculation in the x direction is completed!'
    !===========================================================================================================================================================
    !Calculate the electric conductivity of y-dection
    DO i=X_min,X_max
		DO k=ZZ_min(i),ZZ_max(i)
			Ry=(i-1)*NZB+k
			IF (mod(mmy_per(Ry),2)==0.and. mmy_per(Ry)/=0) THEN  !The case of an even number of intersection points
                idx_start=0
                idx_end=0
                Logic_1=.false.
                Logic_2=.false.
                DO JJ=1,mmy_per(Ry)
                    IF(crosspoint_YY(JJ,Ry)%Log_In)THEN  !The intersection point is the entry point
                        idx_start=JJ
                    ELSEIF(.NOT.crosspoint_YY(JJ,Ry)%Log_In)THEN  !The intersection point is the exit point
                        idx_end=JJ
                    ENDIF
                    IF((idx_start>0).AND.(idx_end>0))THEN   !There are both entry and exit points on the ray simultaneously
                        DO j=1,NY
                            Logic_1=((coordinates_y(j)<=crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y).AND.&
                                    &(coordinates_y(j+1)>=crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y))
							Logic_2=((coordinates_y(j)<=crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y).AND.&
                                    &(coordinates_y(j+1)>=crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y))
							IF(Logic_1) KIdx_1=j
                            IF(Logic_2) KIdx_2=j
                            IF(Logic_1.and.Logic_2)	EXIT
                        ENDDO
                        IF(idx_start<idx_end) THEN  !The ray first penetrate the stratum and enter the anomalous body
							IF(KIdx_1==KIdx_2) THEN 
								LenRatio_CCSIGY(i,KIdx_1,k) = (crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y)/Cdely(KIdx_1)
                                CCSIGY(i,KIdx_1,k) = tao_abnormal * LenRatio_CCSIGY(i,KIdx_1,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGY(i,KIdx_1,k))
                            ELSEIF(KIdx_2 > KIdx_1)THEN
                                LenRatio_CCSIGY(i,KIdx_1,k) = (crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_1))/Cdely(KIdx_1)
                        		LenRatio_CCSIGY(i,KIdx_2,k) = (crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_2))/Cdely(KIdx_2)
                                CCSIGY(i,KIdx_1,k) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGY(i,KIdx_1,k) + tao_abnormal * (1-LenRatio_CCSIGY(i,KIdx_1,k))
                                CCSIGY(i,(KIdx_1+1):(KIdx_2-1),k) = tao_abnormal
                                CCSIGY(i,KIdx_2,k) = tao_abnormal * LenRatio_CCSIGY(i,KIdx_2,k)+TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGY(i,KIdx_2,k))
                            ENDIF
                        ELSE  !The ray first emerges from the anomaly and penetrates the stratum
                            IF(KIdx_1==KIdx_2) THEN
                                LenRatio_CCSIGY(i,KIdx_1,k) = (crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y-crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y)/Cdely(KIdx_1)
                                CCSIGY(i,KIdx_1,k) = tao_abnormal * LenRatio_CCSIGY(i,KIdx_1,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGY(i,KIdx_1,k))
                            ELSEIF(KIdx_2 < KIdx_1)THEN 
                                LenRatio_CCSIGY(i,KIdx_1,k) = (coordinates_y(KIdx_1+1)-crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y)/Cdely(KIdx_1)
                        		LenRatio_CCSIGY(i,KIdx_2,k) = (crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_2))/Cdely(KIdx_2)
                            ENDIF
                            CCSIGY(i,KIdx_2,k) = tao_abnormal*LenRatio_CCSIGY(i,KIdx_2,k)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGY(i,KIdx_2,k))
                            CCSIGY(i,(KIdx_2+1):(KIdx_1-1),k) = TAR_CONDUCTIVITY(2)
                            CCSIGY(i,KIdx_1,k) = TAR_CONDUCTIVITY(2)*LenRatio_CCSIGY(i,KIdx_1,k)+tao_abnormal*(1-LenRatio_CCSIGY(i,KIdx_1,k))
                        ENDIF
                    ENDIF
                ENDDO
            ELSEIF(mod(mmy_per(Ry),2)==1)THEN
                print*,'y-dection!!!ERROR!!!ERROR!!!ERROR!!!'
				print*,i,k,coordinates_x(i),coordinates_z(k)
            ENDIF
        ENDDO
    ENDDO       
	print*,'The equivalent conductivity calculation in the y direction is completed!'
    
	DEALLOCATE(orig_z,orig_y,orig_x)
	DEALLOCATE(vert0,vert1,vert2,edge1,edge2)
    DEALLOCATE(det_z,det_x,det_y)
	DEALLOCATE(u_z,u_x,u_y)
	DEALLOCATE(v_z,v_x,v_y)
	DEALLOCATE(t_z,t_x,t_y)
	DEALLOCATE(coor_z,coor_y,coor_x)
	DEALLOCATE(pvec_z,pvec_y,pvec_x)
    DEALLOCATE(tvec_z,tvec_y,tvec_x)
	DEALLOCATE(ZZ_min,ZZ_max,XX_min,XX_max,YY_min,YY_max)
	DEALLOCATE(coor_z_min,coor_z_max,coor_x_min,coor_x_max,coor_y_min,coor_y_max)
	DEALLOCATE(mmz_per,mmy_per,mmx_per)
	DEALLOCATE(crosspoint_ZZ,crosspoint_YY,crosspoint_XX)
		
    
	RETURN
END SUBROUTINE anomalous_conformal