!> CSV 応答テーブルの読み込み・格子検証を担う。 submodule(bem_matching_plane_response) bem_matching_plane_response_io use bem_string_utils, only: is_decimal_real_token implicit none integer, parameter :: response_column_count = 11 integer, parameter :: response_line_length = 4096 integer, parameter :: initial_row_capacity = 64 character(len=*), parameter :: matching_plane_z_prefix = '# matching_plane_z_m=' type :: matching_plane_axis_type real(dp), allocatable :: values(:) end type matching_plane_axis_type contains module procedure read_matching_plane_response_csv real(dp), allocatable :: rows(:, :) real(dp) :: parsed_row(response_column_count), matching_plane_z_m character(len=response_line_length) :: line character(len=:), allocatable :: left_adjusted integer :: unit_id, ios, line_number, row_count, capacity logical :: header_found, matching_plane_z_found table = matching_plane_response_table_type() call accept(status, message) capacity = initial_row_capacity allocate (rows(response_column_count, capacity)) row_count = 0 line_number = 0 header_found = .false. matching_plane_z_found = .false. matching_plane_z_m = 0.0_dp open (newunit=unit_id, file=trim(path), status='old', action='read', iostat=ios) if (ios /= 0) then call reject( & matching_plane_response_io_error, 'could not open matching-plane response CSV: '//trim(path), status, message & ) return end if do read (unit_id, '(a)', iostat=ios) line if (ios < 0) exit line_number = line_number + 1 if (ios > 0) then close (unit_id) call reject_line( & matching_plane_response_io_error, line_number, 'failed to read matching-plane response CSV.', status, message & ) return end if if (len_trim(line) == 0) cycle left_adjusted = adjustl(trim(line)) if (left_adjusted(1:1) == '#') then if (index(left_adjusted, '# matching_plane_z_m') == 1) then if (header_found) then close (unit_id) call reject_line( & matching_plane_response_invalid_metadata, line_number, & 'matching_plane_z_m metadata must precede the CSV header.', status, message & ) return end if if (matching_plane_z_found) then close (unit_id) call reject_line( & matching_plane_response_invalid_metadata, line_number, & 'matching_plane_z_m metadata must appear exactly once.', status, message & ) return end if call parse_matching_plane_z(left_adjusted, matching_plane_z_m, status, message) if (status /= matching_plane_response_ok) then close (unit_id) call prefix_line_number(line_number, message) return end if matching_plane_z_found = .true. end if cycle end if if (.not. header_found) then if (trim(line) /= matching_plane_response_csv_header) then close (unit_id) call reject_line( & matching_plane_response_invalid_header, line_number, & 'matching-plane response CSV header does not match the exact v1 contract.', status, message & ) return end if if (.not. matching_plane_z_found) then close (unit_id) call reject_line( & matching_plane_response_invalid_metadata, line_number, & 'matching_plane_z_m metadata is required before the CSV header.', status, message & ) return end if header_found = .true. cycle end if call parse_numeric_csv_row(line, parsed_row, status, message) if (status /= matching_plane_response_ok) then close (unit_id) call prefix_line_number(line_number, message) return end if if (row_count >= capacity) call grow_row_buffer(rows, capacity) row_count = row_count + 1 rows(:, row_count) = parsed_row end do close (unit_id) if (.not. header_found) then call reject( & matching_plane_response_invalid_header, 'matching-plane response CSV contains no v1 header.', status, message & ) return end if if (row_count <= 0) then call reject( & matching_plane_response_invalid_grid, 'matching-plane response CSV contains no numeric rows.', status, message & ) return end if call build_matching_plane_response_table( & path, matching_plane_z_m, rows(:, :row_count), table, status, message & ) end procedure read_matching_plane_response_csv subroutine build_matching_plane_response_table(path, matching_plane_z_m, rows, table, status, message) character(len=*), intent(in) :: path real(dp), intent(in) :: matching_plane_z_m real(dp), intent(in) :: rows(:, :) type(matching_plane_response_table_type), intent(out) :: table integer(i32), intent(out) :: status character(len=*), intent(out) :: message type(matching_plane_axis_type) :: axis_data(matching_plane_response_input_count) logical, allocatable :: filled(:) integer(i32) :: coordinate_index(matching_plane_response_input_count) integer(i64) :: expected_points, linear_index, stride integer :: axis, row, max_axis_size, point_count table = matching_plane_response_table_type() call accept(status, message) if (size(rows, 1) /= response_column_count .or. size(rows, 2) <= 0) then call reject( & matching_plane_response_invalid_grid, 'matching-plane response rows have an invalid shape.', status, message & ) return end if if (.not. ieee_is_finite(matching_plane_z_m)) then call reject( & matching_plane_response_invalid_metadata, 'matching_plane_z_m metadata must be finite.', status, message & ) return end if if (any(.not. ieee_is_finite(rows))) then call reject( & matching_plane_response_invalid_row, 'matching-plane response values must all be finite.', status, message & ) return end if if (any(rows(2:5, :) < 0.0_dp) .or. any(rows(7:8, :) < 0.0_dp)) then call reject( & matching_plane_response_invalid_row, & 'matching-plane response fluxes and normal energies must be non-negative.', status, message & ) return end if expected_points = 1_i64 do axis = 1, matching_plane_response_input_count call unique_sorted_values(rows(axis, :), axis_data(axis)%values) table%axis_sizes(axis) = int(size(axis_data(axis)%values), i32) if (axis >= matching_plane_input_photoelectron_outward_flux .and. & table%axis_sizes(axis) > 1_i32 .and. & .not. any(axis_data(axis)%values == 0.0_dp)) then call reject( & matching_plane_response_invalid_grid, & 'every matching-plane feedback axis must include zero.', status, message & ) return end if if (expected_points > huge(expected_points)/int(table%axis_sizes(axis), i64)) then call reject( & matching_plane_response_invalid_grid, 'matching-plane Cartesian grid size overflow.', status, message & ) return end if expected_points = expected_points*int(table%axis_sizes(axis), i64) end do if (expected_points > int(huge(0), i64)) then call reject( & matching_plane_response_invalid_grid, 'matching-plane Cartesian grid is too large.', status, message & ) return end if if (expected_points /= int(size(rows, 2), i64)) then call reject( & matching_plane_response_invalid_grid, & 'matching-plane response rows do not form one complete Cartesian product.', status, message & ) return end if max_axis_size = maxval(table%axis_sizes) point_count = int(expected_points) allocate (table%axes(max_axis_size, matching_plane_response_input_count)) table%axes = 0.0_dp do axis = 1, matching_plane_response_input_count table%axes(:table%axis_sizes(axis), axis) = axis_data(axis)%values end do allocate (table%response_values(matching_plane_response_output_count, point_count)) allocate (filled(point_count), source=.false.) table%response_values = 0.0_dp do row = 1, size(rows, 2) do axis = 1, matching_plane_response_input_count coordinate_index(axis) = exact_axis_index(axis_data(axis)%values, rows(axis, row)) if (coordinate_index(axis) <= 0_i32) then call reject( & matching_plane_response_invalid_grid, & 'matching-plane row could not be mapped to its canonical axis.', status, message & ) return end if end do linear_index = 1_i64 stride = 1_i64 do axis = 1, matching_plane_response_input_count linear_index = linear_index + int(coordinate_index(axis) - 1_i32, i64)*stride stride = stride*int(table%axis_sizes(axis), i64) end do if (filled(int(linear_index))) then call reject( & matching_plane_response_invalid_grid, & 'matching-plane response grid contains a duplicate input coordinate.', status, message & ) return end if filled(int(linear_index)) = .true. table%response_values(:, int(linear_index)) = rows(6:11, row) end do if (size(rows, 2) /= point_count .or. .not. all(filled)) then call reject( & matching_plane_response_invalid_grid, & 'matching-plane response rows must form one complete Cartesian product.', status, message & ) return end if table%source_path = trim(path) table%matching_plane_z_m = matching_plane_z_m table%loaded = .true. end subroutine build_matching_plane_response_table subroutine parse_matching_plane_z(line, matching_plane_z_m, status, message) character(len=*), intent(in) :: line real(dp), intent(out) :: matching_plane_z_m integer(i32), intent(out) :: status character(len=*), intent(out) :: message character(len=:), allocatable :: value_text integer :: ios matching_plane_z_m = 0.0_dp call accept(status, message) if (index(line, matching_plane_z_prefix) /= 1) then call reject( & matching_plane_response_invalid_metadata, & 'matching_plane_z_m metadata must use "# matching_plane_z_m=<finite>".', status, message & ) return end if value_text = trim(adjustl(line(len(matching_plane_z_prefix) + 1:))) if (len(value_text) == 0 .or. scan(value_text, ' ,'//achar(9)) > 0 .or. & .not. is_decimal_real_token(value_text)) then call reject( & matching_plane_response_invalid_metadata, 'matching_plane_z_m metadata has an invalid value.', status, message & ) return end if read (value_text, *, iostat=ios) matching_plane_z_m if (ios /= 0 .or. .not. ieee_is_finite(matching_plane_z_m)) then matching_plane_z_m = 0.0_dp call reject( & matching_plane_response_invalid_metadata, 'matching_plane_z_m metadata must be finite.', status, message & ) end if end subroutine parse_matching_plane_z subroutine parse_numeric_csv_row(line, values, status, message) character(len=*), intent(in) :: line real(dp), intent(out) :: values(response_column_count) integer(i32), intent(out) :: status character(len=*), intent(out) :: message character(len=:), allocatable :: record, token integer :: column, comma, first, last, ios values = 0.0_dp call accept(status, message) record = trim(line) first = 1 do column = 1, response_column_count if (column < response_column_count) then comma = index(record(first:), ',') if (comma == 0) then call reject( & matching_plane_response_invalid_row, & 'matching-plane response row must contain exactly 11 comma-separated values.', status, message & ) return end if last = first + comma - 2 else if (index(record(first:), ',') /= 0) then call reject( & matching_plane_response_invalid_row, & 'matching-plane response row must contain exactly 11 comma-separated values.', status, message & ) return end if last = len(record) end if if (last < first) then call reject( & matching_plane_response_invalid_row, 'matching-plane response CSV values must not be empty.', status, message & ) return end if token = trim(adjustl(record(first:last))) if (len(token) == 0 .or. scan(token, ' '//achar(9)) > 0 .or. & .not. is_decimal_real_token(token)) then call reject( & matching_plane_response_invalid_row, 'matching-plane response CSV contains an invalid numeric token.', & status, message & ) return end if read (token, *, iostat=ios) values(column) if (ios /= 0) then call reject( & matching_plane_response_invalid_row, 'matching-plane response CSV contains an invalid numeric token.', & status, message & ) return end if first = last + 2 end do if (any(.not. ieee_is_finite(values))) then call reject( & matching_plane_response_invalid_row, 'matching-plane response values must all be finite.', status, message & ) end if end subroutine parse_numeric_csv_row subroutine grow_row_buffer(rows, capacity) real(dp), allocatable, intent(inout) :: rows(:, :) integer, intent(inout) :: capacity real(dp), allocatable :: grown(:, :) integer :: new_capacity if (capacity > huge(capacity)/2) error stop 'matching-plane response row capacity overflow.' new_capacity = 2*capacity allocate (grown(response_column_count, new_capacity)) grown(:, :capacity) = rows call move_alloc(grown, rows) capacity = new_capacity end subroutine grow_row_buffer subroutine unique_sorted_values(values, unique) real(dp), intent(in) :: values(:) real(dp), allocatable, intent(out) :: unique(:) real(dp), allocatable :: sorted(:), compact(:) integer :: i, count allocate (sorted, source=values) call merge_sort_real(sorted) allocate (compact(size(sorted))) count = 1 compact(1) = sorted(1) do i = 2, size(sorted) if (sorted(i) /= compact(count)) then count = count + 1 compact(count) = sorted(i) end if end do allocate (unique(count), source=compact(:count)) end subroutine unique_sorted_values subroutine merge_sort_real(values) real(dp), intent(inout) :: values(:) real(dp), allocatable :: work(:) integer :: width, left, middle, right, i, j, k, n n = size(values) if (n <= 1) return allocate (work(n)) width = 1 do while (width < n) left = 1 do while (left <= n) middle = min(left + width, n + 1) right = min(left + 2*width - 1, n) i = left j = middle do k = left, right if (i >= middle) then work(k) = values(j) j = j + 1 else if (j > right) then work(k) = values(i) i = i + 1 else if (values(i) <= values(j)) then work(k) = values(i) i = i + 1 else work(k) = values(j) j = j + 1 end if end do left = left + 2*width end do values = work if (width > n/2) exit width = 2*width end do end subroutine merge_sort_real pure integer(i32) function exact_axis_index(axis_values, value) result(index_value) real(dp), intent(in) :: axis_values(:), value integer :: low, high, middle index_value = 0_i32 low = 1 high = size(axis_values) do while (low <= high) middle = (low + high)/2 if (value == axis_values(middle)) then index_value = int(middle, i32) return else if (value < axis_values(middle)) then high = middle - 1 else low = middle + 1 end if end do end function exact_axis_index subroutine reject_line(code, line_number, text, status, message) integer(i32), intent(in) :: code integer, intent(in) :: line_number character(len=*), intent(in) :: text integer(i32), intent(out) :: status character(len=*), intent(out) :: message character(len=512) :: contextual write (contextual, '(a,i0,a,a)') 'line ', line_number, ': ', trim(text) call reject(code, trim(contextual), status, message) end subroutine reject_line subroutine prefix_line_number(line_number, message) integer, intent(in) :: line_number character(len=*), intent(inout) :: message character(len=512) :: contextual write (contextual, '(a,i0,a,a)') 'line ', line_number, ': ', trim(message) message = trim(contextual) end subroutine prefix_line_number module procedure accept status = matching_plane_response_ok message = '' end procedure accept module procedure reject status = code message = text end procedure reject end submodule bem_matching_plane_response_io