!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
!Code distribution @ https://git.em3d.cn/

SUBROUTINE GET_COORDINATES
	!> @brief This subroutine calculates the coordinates of each grid node in the 3D domain,
	!>including the Yee grid nodes and the source-centered coordinate system.
	USE CONSTANTPARAMETERS
	IMPLICIT NONE
	
	INTEGER(KIND=4) :: ii,jj,kk
	!================================================
	Coordiz(NZS) = -Cdelz(NZS)/2
	!>In the Coordiz coordinate system, the coordinates are on the grid edge
	!! Coordinates are defined at cell centers relative to source location
    do kk=NZS-1,1,-1
        Coordiz(kk)=Coordiz(kk+1)-(Cdelz(kk+1)+Cdelz(kk))/2
    end do
    do kk=NZS+1,NZ,1
        Coordiz(kk)=Coordiz(kk-1)+(Cdelz(kk-1)+Cdelz(kk))/2
    end do
    !> Initialize X and Y coordinates based on source grid length being odd or even
    IF(Logi_Sourcelenth) THEN		
		!> For odd source grid length, origin is at the grid center
        Coordix(NXS)=0
        do ii=NXS-1,1,-1
            Coordix(ii)=Coordix(ii+1)-(Cdelx(ii)+Cdelx(ii+1))/2
        end do
        do ii=NXS+1,NX,1
            Coordix(ii)=Coordix(ii-1)+(Cdelx(ii-1)+Cdelx(ii))/2
        end do

        Coordiy(NYS)=0
        do jj=NYS-1,1,-1
            Coordiy(jj)=Coordiy(jj+1)-(Cdely(jj)+Cdely(jj+1))/2
        end do
        do jj=NYS+1,NY,1
            Coordiy(jj)=Coordiy(jj-1)+(Cdely(jj)+Cdely(jj-1))/2
        end do
    ELSE
		!> For even source grid length, origin spans between two grid centers
        Coordix(NXS)=-GridSize/2.0
        Coordix(NXS+1)=GridSize/2.0
        do ii=NXS-1,1,-1
            Coordix(ii)=Coordix(ii+1)-(Cdelx(ii)+Cdelx(ii+1))/2
        end do
        do ii=NXS+2,nx,1
            Coordix(ii)=Coordix(ii-1)+(Cdelx(ii-1)+Cdelx(ii))/2
        end do

        Coordiy(NYS)=-GridSize/2.0
        Coordiy(NYS+1)=GridSize/2.0
        do jj=NYS-1,1,-1
            Coordiy(jj)=Coordiy(jj+1)-(Cdely(jj)+Cdely(jj+1))/2
        end do
        do jj=NYS+2,ny,1
            Coordiy(jj)=Coordiy(jj-1)+(Cdely(jj)+Cdely(jj-1))/2
        end do		
    ENDIF
	
	!===============================Calculate the Yee node coordinates========================================
	!> Allocate and calculate Yee grid node coordinates
	ALLOCATE(coordinates_x(NXB),coordinates_y(NYB),coordinates_z(NZB))

	DO ii=1,NX
        coordinates_x(ii)=Coordix(ii)-Cdelx(ii)/2.0
    ENDDO
    coordinates_x(NXB)=Coordix(NX)+Cdelx(NX)/2.0
		
    DO jj=1,NY
        coordinates_y(jj)=Coordiy(jj)-Cdely(jj)/2.0
    ENDDO
    coordinates_y(NYB)=Coordiy(NY)+Cdely(NY)/2.0
        
	DO kk=1,NZ
        coordinates_z(kk)=Coordiz(kk)-Cdelz(kk)/2
    ENDDO
	coordinates_z(NZB)=coordinates_z(NZ)+Cdelz(NZ)
	
		!---------------------Create a global coordinate system about HZ-----------------------!
	! The HZ grid planes are located at the Yee nodes, whose source-centered
	! coordinates are stored in coordinates_x/y/z (origin at the loop source /
	! ground surface). The receiver coordinates read from input.dat are also
	! source-centered, so Coord_HZ_* must use the SAME origin, otherwise the
	! receiver-to-grid search in Get_Receiver_Gridlabel yields index 0 (or an
	! uninitialized value) and the observer interpolation in Iteration.f90 reads
	! out-of-bounds EX/EY/CDELX/CDELY entries -> NaN in the dBzdt output files.
	ALLOCATE(Coord_HZ_X(Nx),Coord_HZ_Y(NY),Coord_HZ_Z(NZB))
	Coord_HZ_X(1:NX)  = coordinates_x(1:NX)    !Record the HZ coordinate information in the x direction.
	Coord_HZ_Y(1:NY)  = coordinates_y(1:NY)    !Record the HZ coordinate information in the y direction.
	Coord_HZ_Z(1:NZB) = coordinates_z(1:NZB)   !Record the HZ coordinate information in the z direction.
END SUBROUTINE