subroutine disk_grid ( n, r, c, ng, cg ) !*****************************************************************************80 ! !! DISK_GRID computes grid points inside a disk. ! ! Discussion: ! ! The grid is defined by specifying the radius and center of the circle, ! and the number of subintervals N into which the horizontal radius ! should be divided. Thus, a value of N = 2 will result in 5 points ! along that horizontal line. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 09 November 2011 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of subintervals. ! ! Input, real ( kind = rk ) R, the radius of the circle. ! ! Input, real ( kind = rk ) C(2), the coordinates of the center of the circle. ! ! Input, integer NG, the number of grid points, as determined by ! DISK_GRID_COUNT. ! ! Output, real ( kind = rk ) CG(2,NG), the grid points inside the circle. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ng real ( kind = rk ) c(2) real ( kind = rk ) cg(2,ng) integer i integer j integer n integer p real ( kind = rk ) r real ( kind = rk ) x real ( kind = rk ) y p = 0 do j = 0, n i = 0 x = c(1) y = c(2) + r * real ( 2 * j, kind = rk ) / real ( 2 * n + 1, kind = rk ) p = p + 1 cg(1,p) = x cg(2,p) = y if ( 0 < j ) then p = p + 1 cg(1,p) = x cg(2,p) = 2.0D+00 * c(2) - y end if do i = i + 1 x = c(1) + r * real ( 2 * i, kind = rk ) / real ( 2 * n + 1, kind = rk ) if ( r * r < ( x - c(1) )**2 + ( y - c(2) )**2 ) then exit end if p = p + 1 cg(1,p) = x cg(2,p) = y p = p + 1 cg(1,p) = 2.0D+00 * c(1) - x cg(2,p) = y if ( 0 < j ) then p = p + 1 cg(1,p) = x cg(2,p) = 2.0D+00 * c(2) - y p = p + 1 cg(1,p) = 2.0D+00 * c(1) - x cg(2,p) = 2.0D+00 * c(2) - y end if end do end do return end subroutine disk_grid_count ( n, r, c, ng ) !*****************************************************************************80 ! !! DISK_GRID_COUNT counts the grid points inside a disk. ! ! Discussion: ! ! The grid is defined by specifying the radius and center of the circle, ! and the number of subintervals N into which the horizontal radius ! should be divided. Thus, a value of N = 2 will result in 5 points ! along that horizontal line. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 09 November 2011 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of subintervals. ! ! Input, real ( kind = rk ) R, the radius of the circle. ! ! Input, real ( kind = rk ) C(2), the coordinates of the center of the circle. ! ! Output, integer NG, the number of grid points inside ! the circle. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) c(2) integer i integer j integer n integer ng real ( kind = rk ) r real ( kind = rk ) x real ( kind = rk ) y ng = 0 do j = 0, n i = 0 x = c(1) y = c(2) + r * real ( 2 * j, kind = rk ) / real ( 2 * n + 1, kind = rk ) ng = ng + 1 if ( 0 < j ) then ng = ng + 1 end if do i = i + 1 x = c(1) + r * real ( 2 * i, kind = rk ) / real ( 2 * n + 1, kind = rk ) if ( r * r < ( x - c(1) )**2 + ( y - c(2) )**2 ) then exit end if ng = ng + 1 ng = ng + 1 if ( 0 < j ) then ng = ng + 1 ng = ng + 1 end if end do end do return end subroutine disk_grid_fibonacci ( n, r, c, g ) !*****************************************************************************80 ! !! DISK_GRID_FIBONACCI computes Fibonacci grid points inside a disk. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 20 October 2013 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Richard Swinbank, James Purser, ! Fibonacci grids: A novel approach to global modelling, ! Quarterly Journal of the Royal Meteorological Society, ! Volume 132, Number 619, July 2006 Part B, pages 1769-1793. ! ! Parameters: ! ! Input, integer N, the number of points desired. ! ! Input, real ( kind = rk ) R, the radius of the circle. ! ! Input, real ( kind = rk ) C(2), the coordinates of the center of the circle. ! ! Output, real ( kind = rk ) G(2,N), the grid points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) c(2) real ( kind = rk ) g(2,n) real ( kind = rk ) gr real ( kind = rk ) gt integer i real ( kind = rk ) phi real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 real ( kind = rk ) r real ( kind = rk ) r0 r0 = r / sqrt ( real ( n, kind = rk ) - 0.5D+00 ) phi = ( 1.0D+00 + sqrt ( 5.0D+00 ) ) / 2.0D+00 do i = 1, n gr = r0 * sqrt ( real ( i, kind = rk ) - 0.5D+00 ) gt = 2.0D+00 * pi * real ( i, kind = rk ) / phi g(1,i) = c(1) + gr * cos ( gt ) g(2,i) = c(2) + gr * sin ( gt ) end do return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is a value between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is a value between 1 and 99, representing a ! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 October 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IUNIT, the free unit number. ! implicit none integer i integer ios integer iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end subroutine r82vec_print_part ( n, a, max_print, title ) !*****************************************************************************80 ! !! R82VEC_PRINT_PART prints "part" of an R82VEC. ! ! Discussion: ! ! The user specifies MAX_PRINT, the maximum number of lines to print. ! ! If N, the size of the vector, is no more than MAX_PRINT, then ! the entire vector is printed, one entry per line. ! ! Otherwise, if possible, the first MAX_PRINT-2 entries are printed, ! followed by a line of periods suggesting an omission, ! and the last entry. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 09 November 2011 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries of the vector. ! ! Input, real ( kind = rk ) A(2,N), the vector to be printed. ! ! Input, integer MAX_PRINT, the maximum number of lines ! to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a(2,n) integer i integer max_print character ( len = * ) title if ( max_print <= 0 ) then return end if if ( n <= 0 ) then return end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' if ( n <= max_print ) then do i = 1, n write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(1:2,i) end do else if ( 3 <= max_print ) then do i = 1, max_print - 2 write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(1:2,i) end do write ( *, '(a)' ) ' ........ .............. ..............' i = n write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(1:2,i) else do i = 1, max_print - 1 write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a(1:2,i) end do i = max_print write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6,2x,a)' ) i, ':', a(1:2,i), & '...more entries...' end if return end subroutine r8mat_write ( output_filename, m, n, table ) !*****************************************************************************80 ! !! R8MAT_WRITE writes an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! Input, integer M, the spatial dimension. ! ! Input, integer N, the number of points. ! ! Input, real ( kind = rk ) TABLE(M,N), the data. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer j character ( len = * ) output_filename integer output_status integer output_unit character ( len = 30 ) string real ( kind = rk ) table(m,n) ! ! Open the file. ! call get_unit ( output_unit ) open ( unit = output_unit, file = output_filename, & status = 'replace', iostat = output_status ) if ( output_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_WRITE - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the output file "' // & trim ( output_filename ) // '" on unit ', output_unit output_unit = -1 stop end if ! ! Create a format string. ! ! For less precision in the output file, try: ! ! '(', m, 'g', 14, '.', 6, ')' ! if ( 0 < m .and. 0 < n ) then write ( string, '(a1,i8,a1,i8,a1,i8,a1)' ) '(', m, 'g', 24, '.', 16, ')' ! ! Write the data. ! do j = 1, n write ( output_unit, string ) table(1:m,j) end do end if ! ! Close the file. ! close ( unit = output_unit ) return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end