program create_ss_grid

	! reads in the file with polygons "gmt_file.txt" and creates a grid with a given resolution set by the user
	implicit none

	double precision :: latitude_spacing, longitude_spacing, longitude, latitude, minimum_latitude, minimum_longitude
	double precision :: maximum_latitude, maximum_longitude
	double precision :: local_minimum_longitude, local_maximum_longitude, local_minimum_latitude, local_maximum_latitude
	double precision :: local_longitude, local_latitude

	character(len=255) :: input_parameter, dummy
	character(len=255), parameter :: input_file = "gmt_file.txt", output_file = "domains.txt"

	character(len=1) :: bracket

	integer :: istat, number_polygons, polygon_counter, point_counter, maximum_points, number_latitude, number_longitude
	integer :: start_latitude_index, start_longitude_index, end_latitude_index, end_longitude_index
	integer :: local_latitude_counter, local_longitude_counter, latitude_counter, longitude_counter
	integer, parameter :: input_file_unit=10, max_polygons = 10000, output_file_unit=20

	integer, dimension(max_polygons) :: polygon_point_size ! if you have more than 10000 polygons, you are ambitious!

	double precision, allocatable, dimension(:,:) :: latitude_array, longitude_array
	integer, allocatable, dimension(:) :: shear_stess_id_array
	integer, allocatable, dimension(:,:) :: grid

	logical :: inside

	character(len=50), parameter :: output_format = "(F8.3,1X,F7.3,1X,I5)"

	call getarg(1,input_parameter)
	read(input_parameter,*) latitude_spacing
	call getarg(2,input_parameter)
	read(input_parameter,*) longitude_spacing

	! first, find the number of polygons, and the maximum amount of points

	number_polygons = 0
	polygon_point_size = 0



	open(unit=input_file_unit, file=input_file, access="sequential", form="formatted", status="old")

	initial_read: do

		read(input_file_unit,'(A1,A)', iostat=istat) bracket, dummy
		if(istat /=0) THEN !hopefully end of file
			exit initial_read
		endif

		if(bracket == ">") then
			number_polygons = number_polygons + 1

			if(number_polygons > max_polygons) THEN
				write(6,*) "you have to increase the max_polygons parameter in create_ss_grid.f90, and recompile"
				stop
			endif

		else
			polygon_point_size(number_polygons) = polygon_point_size(number_polygons) + 1
		endif

	end do initial_read

	rewind(unit=input_file_unit)

	maximum_points = maxval(polygon_point_size)

	allocate(latitude_array(number_polygons,maximum_points), longitude_array(number_polygons,maximum_points),&
		shear_stess_id_array(number_polygons), stat=istat)
	if(istat /=0) THEN
		write(6,*) "problem allocating arrays"
		stop
	endif


	minimum_latitude = 90.
	minimum_longitude = 180.
	maximum_latitude = -90.
	maximum_longitude = -180.

	do polygon_counter = 1, number_polygons, 1

		read(input_file_unit,*) bracket, shear_stess_id_array(polygon_counter)


		do point_counter = 1, polygon_point_size(polygon_counter), 1

			read(input_file_unit,*) longitude_array(polygon_counter,point_counter), &
						      latitude_array(polygon_counter,point_counter)

			if(latitude_array(polygon_counter,point_counter) < minimum_latitude) THEN
				minimum_latitude = latitude_array(polygon_counter,point_counter)
			endif

			if(latitude_array(polygon_counter,point_counter) > maximum_latitude) THEN
				maximum_latitude = latitude_array(polygon_counter,point_counter)
			endif

			if(longitude_array(polygon_counter,point_counter) < minimum_longitude) THEN
				minimum_longitude = longitude_array(polygon_counter,point_counter)
			endif

			if(longitude_array(polygon_counter,point_counter) > maximum_longitude) THEN
				maximum_longitude = longitude_array(polygon_counter,point_counter)
			endif

		end do
	end do

	close(input_file_unit)


	! expand the region a bit

	minimum_longitude = dble(floor(minimum_longitude)) - 2. * longitude_spacing
	maximum_longitude = dble(ceiling(maximum_longitude)) + 2. * longitude_spacing

	minimum_latitude = dble(floor(minimum_latitude)) - 2. * latitude_spacing
	maximum_latitude = dble(ceiling(maximum_latitude)) + 2. * latitude_spacing




	number_latitude = nint((maximum_latitude - minimum_latitude) / latitude_spacing) + 1
	number_longitude = nint((maximum_longitude - minimum_longitude) / longitude_spacing) + 1


	allocate(grid(number_longitude,number_latitude), stat=istat)
	if(istat /=0) THEN
		write(6,*) "problem allocating arrays"
		stop
	endif

	! note that it is entirely possible that polgons will be skipped if your lat/long spacing is not dense enough

	! also, if your input shapefile has any gaps between polygons, they could also be skipped. 
	! If they polygons overlap, the values could be overwritten

	grid = 0

	do polygon_counter = 1, number_polygons, 1


		local_minimum_longitude = minval(longitude_array(polygon_counter,1:polygon_point_size(polygon_counter)))
		local_maximum_longitude = maxval(longitude_array(polygon_counter,1:polygon_point_size(polygon_counter)))

		local_minimum_latitude = minval(latitude_array(polygon_counter,1:polygon_point_size(polygon_counter)))
		local_maximum_latitude = maxval(latitude_array(polygon_counter,1:polygon_point_size(polygon_counter)))


		local_minimum_longitude = dble(floor(local_minimum_longitude/longitude_spacing)) * longitude_spacing
		local_maximum_longitude = dble(ceiling(local_maximum_longitude/longitude_spacing)) * longitude_spacing

		local_minimum_latitude = dble(floor(local_minimum_latitude/latitude_spacing)) * latitude_spacing
		local_maximum_latitude = dble(ceiling(local_maximum_latitude/latitude_spacing)) * latitude_spacing

		start_latitude_index = nint((local_minimum_latitude - minimum_latitude)/latitude_spacing) + 1
		start_longitude_index = nint((local_minimum_longitude - minimum_longitude)/longitude_spacing) + 1
		end_latitude_index = nint((local_maximum_latitude - minimum_latitude)/latitude_spacing) + 1
		end_longitude_index = nint((local_maximum_longitude - minimum_longitude)/longitude_spacing) + 1
		
		do  local_longitude_counter = start_longitude_index, end_longitude_index
			do local_latitude_counter = start_latitude_index, end_latitude_index

				
				local_longitude = minimum_longitude + dble(local_longitude_counter-1) * longitude_spacing
				local_latitude = minimum_latitude + dble(local_latitude_counter-1) * latitude_spacing

				inside = point_in_polygon(longitude_array(polygon_counter,1:polygon_point_size(polygon_counter)), &
				  latitude_array(polygon_counter,1:polygon_point_size(polygon_counter)), local_longitude, &
				  local_latitude, polygon_point_size(polygon_counter))

				if(inside) THEN
					grid(local_longitude_counter,local_latitude_counter) = shear_stess_id_array(polygon_counter)
				endif

			end do
		end do


	end do


	! output

	open(unit=output_file_unit, file=output_file, access="sequential", form="formatted", status="replace")
	do longitude_counter = 1, number_longitude, 1
		do latitude_counter = 1, number_latitude, 1
			local_longitude = minimum_longitude + dble(longitude_counter-1) * longitude_spacing
			local_latitude = minimum_latitude + dble(latitude_counter-1) * latitude_spacing

			write(output_file_unit,output_format) local_longitude, local_latitude, grid(longitude_counter,latitude_counter) 

		end do
	end do

	close(output_file_unit)


	deallocate(latitude_array, longitude_array,shear_stess_id_array, stat=istat)
	if(istat /=0) THEN
		write(6,*) "problem deallocating arrays"
		stop
	endif

contains


! Written by Evan Gowan, last updated August 24, 2015


logical function point_in_polygon(x_boundary, y_boundary, x, y, number_points)

! this function determines whether a give point (x,y) is within a polygon defined by x_boundary and y_boundary
! number_points is the number of points in the polygon

	
	implicit none

	integer, intent(in) :: number_points
	double precision, intent(in) :: x, y
	double precision, dimension(number_points), intent(in) :: x_boundary, y_boundary

	integer :: current_point, next_point, last_point, crossover_counter
	logical :: found_first, found_last, inside

	found_first = .false.
	found_last = .false.
	inside = .false.

	current_point = 1
	search_boundary: do

		next_point = current_point + 1
	
		if (next_point == number_points) THEN

			found_last = .true.
		endif

! even-odd rule algorithm to determine if the point is inside or outside

		if (min(y_boundary(current_point), y_boundary(next_point)) < y .and.&
		    max(y_boundary(current_point), y_boundary(next_point)) >= y) THEN

			if (x_boundary(current_point) + (y - y_boundary(current_point)) /&
			    (y_boundary(next_point) - y_boundary(current_point)) * &
			    (x_boundary(next_point) - x_boundary(current_point)) < x) THEN

				inside = .not.(inside)

			endif

		endif

		current_point = current_point + 1

		if (found_last) THEN
			exit search_boundary
		endif
		
	
	end do search_boundary

	point_in_polygon = inside

	return
end function point_in_polygon
	


end program create_ss_grid
