program main c*********************************************************************72 c cc ORBITAL uses DISLIN to display a contour plot of Z(X,Y) data. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 20 May 2011 c c Author: c c John Burkardt c c Reference: c c Helmut Michels, c The Data Plotting Software DISLIN - version 10.4, c Shaker Media GmbH, January 2010, c ISBN13: 978-3-86858-517-9. c implicit none integer m parameter ( m = 3 ) integer n parameter ( n = 101 ) integer nn parameter ( nn = 101 * 101 ) integer i integer j integer k integer level integer level_num real level_value real x(n) real xmax real xmin real xyz(m,nn) real y(n) real ymax real ymin real z(n,n) real zmax real zmin write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ORBITAL:' write ( *, '(a)' ) ' FORTRAN77 version' write ( *, '(a)' ) & ' Use DISLIN to make a contour plot of Z(X,Y) data.' c c Read the data. c call r4mat_data_read ( 'orbital.txt', m, nn, xyz ) c c Split the data. c The contouring routine expects that data is along fixed X and Y coordinates, c and so the X and Y data is to be given as vectors, not arrays. c k = 1 do i = 1, n x(i) = xyz(1,k) k = k + 1 end do call r4vec_max ( n, x, xmax ) call r4vec_min ( n, x, xmin ) k = 1 do i = 1, n y(i) = xyz(2,k) k = k + n end do call r4vec_max ( n, y, ymax ) call r4vec_min ( n, y, ymin ) c c Z is a table. c The first dimension should contain values for constant X. c k = 0 do j = 1, n do i = 1, n k = k + 1 z(i,j) = xyz(3,k) end do end do call r4mat_max ( n, n, z, zmax ) call r4mat_min ( n, n, z, zmin ) c c Specify the format of the output file. c call metafl ( 'png' ) c c Indicate that new data overwrites old data. c call filmod ( 'delete' ) c c Specify the name of the output graphics file. c call setfil ( 'orbital.png' ) c c Choose the page size and orientation. c 'USA' is 2160 plot units wide and 2790 plot units high. c 'P' requests LANDSCAPE orientation. c call setpag ( 'usap' ) c c For PNG output, use reverse the default black background to white. c call scrmod ( 'reverse' ) c c Open DISLIN. c call disini ( ) c c Plot a border around the page. c call pagera ( ) c c Use the SIMPLX font. c call simplx ( ) c c Set the axis origin in plot units to the right, and plot units DOWN. c call axspos ( 230, 2500 ) c c Define the X and Y sizes of the axis system in plot units. c call axslen ( 1700, 1700 ) c c Label the X and Y axes. c call name ( 'X axis', 'X' ) call name ( 'Y axis', 'Y' ) c c Relate the physical coordinates to the axes. c call graf ( xmin, xmax, xmin, 1.0, ymin, ymax, ymin, 1.0 ) c c BEGIN LEVEL 2 COMMANDS. c c Define the title. c call titlin ( 'Orbital contour plot', 1 ) call title ( ) c c Set the color to blue. c call color ( 'blue' ) c c Draw the contour plot. c level_num = 10 do level = 1, level_num level_value = ( real ( level_num + 1 - level ) * zmin & + real ( level ) * zmax ) & / real ( level_num + 1 ) call contur ( x, n, y, n, z, level_value ) end do c c End this graph. c call endgrf ( ) c c RETURN FROM LEVEL 2 TO LEVEL 1. c c Close DISLIN. c call disfin ( ) c c Terminate. c write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ORBITAl:' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine ch_cap ( ch ) c*********************************************************************72 c cc CH_CAP capitalizes a single character. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 03 January 2007 c c Author: c c John Burkardt c c Parameters: c c Input/output, character CH, the character to capitalize. c implicit none character ch integer itemp itemp = ichar ( ch ) if ( 97 .le. itemp .and. itemp .le. 122 ) then ch = char ( itemp - 32 ) end if return end function ch_eqi ( c1, c2 ) c*********************************************************************72 c cc CH_EQI is a case insensitive comparison of two characters for equality. c c Example: c c CH_EQI ( 'A', 'a' ) is TRUE. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 03 January 2007 c c Author: c c John Burkardt c c Parameters: c c Input, character C1, C2, the characters to compare. c c Output, logical CH_EQI, the result of the comparison. c implicit none character c1 character c1_cap character c2 character c2_cap logical ch_eqi c1_cap = c1 c2_cap = c2 call ch_cap ( c1_cap ) call ch_cap ( c2_cap ) if ( c1_cap .eq. c2_cap ) then ch_eqi = .true. else ch_eqi = .false. end if return end subroutine ch_to_digit ( c, digit ) c*********************************************************************72 c cc CH_TO_DIGIT returns the integer value of a base 10 digit. c c Example: c c C DIGIT c --- ----- c '0' 0 c '1' 1 c ... ... c '9' 9 c ' ' 0 c 'X' -1 c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 04 August 1999 c c Author: c c John Burkardt c c Parameters: c c Input, character C, the decimal digit, '0' through '9' or blank c are legal. c c Output, integer DIGIT, the corresponding integer value. If C was c 'illegal', then DIGIT is -1. c implicit none character c integer digit if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then digit = ichar ( c ) - 48 else if ( c .eq. ' ' ) then digit = 0 else digit = -1 end if return end subroutine file_column_count ( input_filename, column_num ) c*********************************************************************72 c cc FILE_COLUMN_COUNT counts the number of columns in the first line of a file. c c Discussion: c c The file is assumed to be a simple text file. c c Most lines of the file is presumed to consist of COLUMN_NUM words, c separated by spaces. There may also be some blank lines, and some c comment lines, c which have a "#" in column 1. c c The routine tries to find the first non-comment non-blank line and c counts the number of words in that line. c c If all lines are blanks or comments, it goes back and tries to analyze c a comment line. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 28 April 2008 c c Author: c c John Burkardt c c Parameters: c c Input, character * ( * ) INPUT_FILENAME, the name of the file. c c Output, integer COLUMN_NUM, the number of columns in the file. c implicit none integer column_num logical got_one character * ( * ) input_filename integer input_unit character * ( 255 ) line integer s_len_trim c c Open the file. c call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, & status = 'old', form = 'formatted', access = 'sequential' ) c c Read one line, but skip blank lines and comment lines. c got_one = .false. 10 continue read ( input_unit, '(a)', err = 20 ) line if ( s_len_trim ( line ) .eq. 0 ) then go to 10 end if if ( line(1:1) .eq. '#' ) then go to 10 end if got_one = .true. go to 20 go to 10 20 continue if ( .not. got_one ) then rewind ( input_unit ) 30 continue read ( input_unit, '(a)', err = 40 ) line if ( s_len_trim ( line ) .eq. 0 ) then go to 30 end if got_one = .true. go to 40 go to 30 40 continue end if close ( unit = input_unit ) if ( .not. got_one ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Warning.' write ( *, '(a)' ) ' The file does not contain any data.' column_num = -1 return end if call s_word_count ( line, column_num ) return end subroutine file_row_count ( input_filename, row_num ) c*********************************************************************72 c cc FILE_ROW_COUNT counts the number of row records in a file. c c Discussion: c c It does not count lines that are blank, or that begin with a c comment symbol '#'. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 28 April 2008 c c Author: c c John Burkardt c c Parameters: c c Input, character * ( * ) INPUT_FILENAME, the name of the input file. c c Output, integer ROW_NUM, the number of rows found. c implicit none integer bad_num integer comment_num integer ierror character * ( * ) input_filename integer input_status integer input_unit character * ( 255 ) line integer record_num integer row_num integer s_len_trim call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, & status = 'old' ) comment_num = 0 row_num = 0 record_num = 0 bad_num = 0 10 continue read ( input_unit, '(a)', err = 20, end = 20 ) line record_num = record_num + 1 if ( line(1:1) .eq. '#' ) then comment_num = comment_num + 1 go to 10 end if if ( s_len_trim ( line ) .eq. 0 ) then comment_num = comment_num + 1 go to 10 end if row_num = row_num + 1 go to 10 20 continue close ( unit = input_unit ) return end subroutine get_unit ( iunit ) c*********************************************************************72 c cc GET_UNIT returns a free FORTRAN unit number. c c Discussion: c c A "free" FORTRAN unit number is a value between 1 and 99 which c is not currently associated with an I/O device. A free FORTRAN unit c number is needed in order to open a file with the OPEN command. c c If IUNIT = 0, then no free FORTRAN unit could be found, although c all 99 units were checked (except for units 5, 6 and 9, which c are commonly reserved for console I/O). c c Otherwise, IUNIT is a value between 1 and 99, representing a c free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 c are special, and will never return those values. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 02 September 2013 c c Author: c c John Burkardt c c Parameters: c c Output, integer IUNIT, the free unit number. c implicit none integer i integer iunit logical value iunit = 0 do i = 1, 99 if ( i .ne. 5 .and. i .ne. 6 .and. i .ne. 9 ) then inquire ( unit = i, opened = value, err = 10 ) if ( .not. value ) then iunit = i return end if end if 10 continue end do return end subroutine r4mat_data_read ( input_filename, m, n, table ) c*********************************************************************72 c cc R4MAT_DATA_READ reads data from an R4MAT file. c c Discussion: c c An R4MAT is an array of R4's. c c The file may contain more than N points, but this routine will c return after reading N of them. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 28 April 2008 c c Author: c c John Burkardt c c Parameters: c c Input, character * ( * ) INPUT_FILENAME, the name of the input file. c c Input, integer M, the spatial dimension. c c Input, integer N, the number of points. c c Output, real TABLE(M,N), the data. c implicit none integer m integer n integer i integer ierror character * ( * ) input_filename integer input_status integer input_unit integer j character * ( 255 ) line integer s_len_trim real table(m,n) real x(m) ierror = 0 call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, & status = 'old' ) j = 0 10 continue if ( j .lt. n ) then read ( input_unit, '(a)' ) line if ( line(1:1) .eq. '#' .or. s_len_trim ( line ) .eq. 0 ) then go to 10 end if call s_to_r4vec ( line, m, x, ierror ) if ( ierror .ne. 0 ) then go to 10 end if j = j + 1 do i = 1, m table(i,j) = x(i) end do go to 10 end if close ( unit = input_unit ) return end subroutine r4mat_header_read ( input_filename, m, n ) c*********************************************************************72 c cc R4MAT_HEADER_READ reads the header from an R4MAT file. c c Discussion: c c An R4MAT is an array of R4's. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 28 April 2008 c c Author: c c John Burkardt c c Parameters: c c Input, character * ( * ) INPUT_FILENAME, the name of the input file. c c Output, integer M, spatial dimension. c c Output, integer N, the number of points. c implicit none character * ( * ) input_filename integer m integer n call file_column_count ( input_filename, m ) if ( m .le. 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R4MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was an I/O problem while trying' write ( *, '(a)' ) ' to count the number of data columns in' write ( *, '(a,a,a)' ) & ' the file "', trim ( input_filename ), '".' stop end if call file_row_count ( input_filename, n ) if ( n .le. 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R4MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was an I/O problem while trying' write ( *, '(a)' ) ' to count the number of data rows in' write ( *, '(a,a,a)' ) & ' the file "', trim ( input_filename ), '".' stop end if return end subroutine r4mat_max ( m, n, a, amax ) c*********************************************************************72 c cc R4MAT_MAX returns the maximum entry of an R4MAT. c c Discussion: c c An R4MAT is an array of R4's. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 19 May 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer M, the number of rows in A. c c Input, integer N, the number of columns in A. c c Input, real A(M,N), the matrix. c c Output, real AMAX, the maximum entry of A. c implicit none integer m integer n real a(m,n) real amax integer i integer j amax = a(1,1) do j = 1, n do i = 1, m amax = max ( amax, a(i,j) ) end do end do return end subroutine r4mat_min ( m, n, a, amin ) c*********************************************************************72 c cc R4MAT_MIN returns the minimum entry of an R4MAT. c c Discussion: c c An R4MAT is an array of R4's. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 19 May 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer M, the number of rows in A. c c Input, integer N, the number of columns in A. c c Input, real A(M,N), the matrix. c c Output, real AMIN, the minimum entry of A. c implicit none integer m integer n real a(m,n) real amin integer i integer j amin = a(1,1) do j = 1, n do i = 1, m amin = min ( amin, a(i,j) ) end do end do return end subroutine r4vec_max ( n, a, amax ) c*********************************************************************72 c cc R4VEC_MAX returns the maximum value in an R4VEC. c c Discussion: c c An R4VEC is a vector of R4's. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 31 May 2009 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of entries in the array. c c Input, real A(N), the array. c c Output, real AMAX, the value of the largest entry. c implicit none integer n real a(n) real amax integer i amax = a(1) do i = 2, n amax = max ( amax, a(i) ) end do return end subroutine r4vec_min ( n, a, amin ) c*********************************************************************72 c cc R4VEC_MIN returns the minimum value in an R4VEC. c c Discussion: c c An R4VEC is a vector of R4's. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 31 May 2009 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of entries in the array. c c Input, real A(N), the array. c c Output, real AMIN, the value of the smallest entry. c implicit none integer n real a(n) real amin integer i amin = a(1) do i = 2, n amin = min ( amin, a(i) ) end do return end function s_len_trim ( s ) c*********************************************************************72 c cc S_LEN_TRIM returns the length of a string to the last nonblank. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 05 March 2004 c c Author: c c John Burkardt c c Parameters: c c Input, character*(*) S, a string. c c Output, integer S_LEN_TRIM, the length of the string to the last nonblank. c implicit none integer i character*(*) s integer s_len_trim do i = len ( s ), 1, -1 if ( s(i:i) .ne. ' ' ) then s_len_trim = i return end if end do s_len_trim = 0 return end subroutine s_to_r4 ( s, dval, ierror, length ) c*********************************************************************72 c cc S_TO_R4 reads an R4 from a string. c c Discussion: c c The routine will read as many characters as possible until it reaches c the end of the string, or encounters a character which cannot be c part of the number. c c Legal input is: c c 1 blanks, c 2 '+' or '-' sign, c 2.5 blanks c 3 integer part, c 4 decimal point, c 5 fraction part, c 6 'E' or 'e' or 'D' or 'd', exponent marker, c 7 exponent sign, c 8 exponent integer part, c 9 exponent decimal point, c 10 exponent fraction part, c 11 blanks, c 12 final comma or semicolon, c c with most quantities optional. c c Example: c c S DVAL c c '1' 1.0 c ' 1 ' 1.0 c '1A' 1.0 c '12,34,56' 12.0 c ' 34 7' 34.0 c '-1E2ABCD' -100.0 c '-1X2ABCD' -1.0 c ' 2E-1' 0.2 c '23.45' 23.45 c '-4.2E+2' -420.0 c '17d2' 1700.0 c '-14e-2' -0.14 c 'e2' 100.0 c '-12.73e-9.23' -12.73 * 10.0**(-9.23) c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 02 January 2009 c c Author: c c John Burkardt c c Parameters: c c Input, character * ( * ) S, the string containing the c data to be read. Reading will begin at position 1 and c terminate at the end of the string, or when no more c characters can be read to form a legal real. Blanks, c commas, or other nonnumeric data will, in particular, c cause the conversion to halt. c c Output, real DVAL, the value read from the string. c c Output, integer IERROR, error flag. c 0, no errors occurred. c 1, 2, 6 or 7, the input number was garbled. The c value of IERROR is the last type of input successfully c read. For instance, 1 means initial blanks, 2 means c a plus or minus sign, and so on. c c Output, integer LENGTH, the number of characters read c to form the number, including any terminating c characters such as a trailing comma or blanks. c implicit none logical ch_eqi character c real dval integer ierror integer ihave integer isgn integer iterm integer jbot integer jsgn integer jtop integer length integer nchar integer ndig real rbot real rexp real rtop character * ( * ) s integer s_len_trim nchar = s_len_trim ( s ) ierror = 0 dval = 0.0E+00 length = -1 isgn = 1 rtop = 0 rbot = 1 jsgn = 1 jtop = 0 jbot = 1 ihave = 1 iterm = 0 10 continue length = length + 1 if ( nchar .lt. length+1 ) then go to 20 end if c = s(length+1:length+1) c c Blank character. c if ( c .eq. ' ' ) then if ( ihave .eq. 2 ) then else if ( ihave .eq. 6 .or. ihave .eq. 7 ) then iterm = 1 else if ( 1 .lt. ihave ) then ihave = 11 end if c c Comma. c else if ( c .eq. ',' .or. c .eq. ';' ) then if ( ihave .ne. 1 ) then iterm = 1 ihave = 12 length = length + 1 end if c c Minus sign. c else if ( c .eq. '-' ) then if ( ihave .eq. 1 ) then ihave = 2 isgn = -1 else if ( ihave .eq. 6 ) then ihave = 7 jsgn = -1 else iterm = 1 end if c c Plus sign. c else if ( c .eq. '+' ) then if ( ihave .eq. 1 ) then ihave = 2 else if ( ihave .eq. 6 ) then ihave = 7 else iterm = 1 end if c c Decimal point. c else if ( c .eq. '.' ) then if ( ihave .lt. 4 ) then ihave = 4 else if ( 6 .le. ihave .and. ihave .le. 8 ) then ihave = 9 else iterm = 1 end if c c Scientific notation exponent marker. c else if ( ch_eqi ( c, 'E' ) .or. ch_eqi ( c, 'D' ) ) then if ( ihave .lt. 6 ) then ihave = 6 else iterm = 1 end if c c Digit. c else if ( ihave .lt. 11 .and. lle ( '0', c ) & .and. lle ( c, '9' ) ) then if ( ihave .le. 2 ) then ihave = 3 else if ( ihave .eq. 4 ) then ihave = 5 else if ( ihave .eq. 6 .or. ihave .eq. 7 ) then ihave = 8 else if ( ihave .eq. 9 ) then ihave = 10 end if call ch_to_digit ( c, ndig ) if ( ihave .eq. 3 ) then rtop = 10.0E+00 * rtop + real ( ndig ) else if ( ihave .eq. 5 ) then rtop = 10.0E+00 * rtop + real ( ndig ) rbot = 10.0E+00 * rbot else if ( ihave .eq. 8 ) then jtop = 10 * jtop + ndig else if ( ihave .eq. 10 ) then jtop = 10 * jtop + ndig jbot = 10 * jbot end if c c Anything else is regarded as a terminator. c else iterm = 1 end if c c If we haven't seen a terminator, and we haven't examined the c entire string, go get the next character. c if ( iterm .eq. 1 ) then go to 20 end if go to 10 20 continue c c If we haven't seen a terminator, and we have examined the c entire string, then we're done, and LENGTH is equal to NCHAR. c if ( iterm .ne. 1 .and. length+1 .eq. nchar ) then length = nchar end if c c Number seems to have terminated. Have we got a legal number? c Not if we terminated in states 1, 2, 6 or 7. c if ( ihave .eq. 1 .or. ihave .eq. 2 .or. & ihave .eq. 6 .or. ihave .eq. 7 ) then ierror = ihave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'S_TO_R4 - Serious error!' write ( *, '(a)' ) ' Illegal or nonnumeric input:' write ( *, '(a,a)' ) ' ', s return end if c c Number seems OK. Form it. c if ( jtop .eq. 0 ) then rexp = 1.0E+00 else if ( jbot .eq. 1 ) then rexp = 10.0E+00 ** ( jsgn * jtop ) else rexp = 10.0E+00 ** ( real ( jsgn * jtop ) / real ( jbot ) ) end if end if dval = real ( isgn ) * rexp * rtop / rbot return end subroutine s_to_r4vec ( s, n, rvec, ierror ) c*********************************************************************72 c cc S_TO_R4VEC reads an R4VEC from a string. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 02 January 2009 c c Author: c c John Burkardt c c Parameters: c c Input, character * ( * ) S, the string to be read. c c Input, integer N, the number of values expected. c c Output, real RVEC(N), the values read from the string. c c Output, integer IERROR, error flag. c 0, no errors occurred. c -K, could not read data for entries -K through N. c implicit none integer n integer i integer ierror integer ilo integer lchar real rvec(n) character * ( * ) s i = 0 ierror = 0 ilo = 1 10 continue if ( i .lt. n ) then i = i + 1 call s_to_r4 ( s(ilo:), rvec(i), ierror, lchar ) if ( ierror .ne. 0 ) then ierror = -i go to 20 end if ilo = ilo + lchar go to 10 end if 20 continue return end subroutine s_word_count ( s, nword ) c*********************************************************************72 c cc S_WORD_COUNT counts the number of "words" in a string. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 28 April 2008 c c Author: c c John Burkardt c c Parameters: c c Input, character * ( * ) S, the string to be examined. c c Output, integer NWORD, the number of "words" in the string. c Words are presumed to be separated by one or more blanks. c implicit none logical blank integer i integer lens integer nword character * ( * ) s nword = 0 lens = len ( s ) if ( lens .le. 0 ) then return end if blank = .true. do i = 1, lens if ( s(i:i) .eq. ' ' ) then blank = .true. else if ( blank ) then nword = nword + 1 blank = .false. end if end do return end