program poly_area * Program to compute the area of an arbitrarily shaped * Polygon but breaking the polygon into triangles. * MOD TAH: Version 2: * Added a check of the sign of the darea calculation. * If the sign changes as we compute each triangle then * if means the figure has not been input correctly. * PARAMETERS integer*4 max_nodes ! Maximum number of nodes in the polygon ! that the user can input. parameter ( max_nodes = 1000 ) * MAIN PROGRAM VARIABLES integer*4 num_nodes ! Number of nodes input by the user. List ! of node coordinates is read until end-of- ! file (EOF or ^D) is reached. integer*4 sign_darea ! Sign of the area of the triangle. If ! this changes then figure is not ordered ! correctly. real*8 nodes_xy(2,max_nodes) ! XY coordinates of nodes of the ! of polygon real*8 triangle_vec(2,2) ! Two vectors that make up the ! current triangle. Lead index is over ! XY and the second index is over side ! 1 and 2. real*8 area, darea ! Total area and increment on area for ! each triangle. These can be positive or ! negative depending whether nodes are ! entered clockwise or anti-clockwise. ! Sign is fixed in output routine. character*8 units ! Units of the coordinates. Used only for ! output. * Miscellaneous variables that are needed integer*4 i ! Counter for looping over the nodes ***** START PROGRAM. * Tell user what this program does write(*,110) 110 format(/,' POLY_AREA: Program to compute the area of a', . ' plane polygon') * Now get the polygon node coordinates call read_nodes(num_nodes, nodes_xy, units, max_nodes) * Make sure we have enough nodes to from a figure. if( num_nodes.le.2 ) then write(*,130) num_nodes 130 format('**ERROR** Only ',i2,' nodes entered. This is', . ' not enough to form figure') * Some programmers considering putting "stops" inside * routines a bad practice, and that in this case, I * should execute the remainder of the program in an * "else" structure. I feel that if the program has reached * a dead-end, it should stop to avoid the problem that with * later modifications, it may continue to run. stop 'Insufficient number of nodes' end if * We have enough nodes, now start to loop over the triangles * in the figure computing the incremental area and summing to * get total. The order here is that the first node will be * apex of all triangles, and in each loop we form the triangle * of apex-node_i and apex-node_i-1. (Could have started at * 2 and gone to num_nodes-1). * Initialize the total area to 0 before starting area = 0.0d0 do i = 3, num_nodes * Form the two vectors that make up the triangle call form_triangle(i, nodes_xy, triangle_vec, num_nodes) * Now compute the increment on the area call triangle_darea(triangle_vec, darea) * MOD TAH Version 2: If this is first triangle set the sign, on * subsequent triangles check to see that it remains the * same call check_dir(i, sign_darea, darea) * Increment into total area area = area + darea end do * We are now complete, write out the results call output_area(num_nodes, area, nodes_xy, units) **** We are finished end CTITLE READ_NODES subroutine read_nodes(num_nodes, nodes_xy, units, max_nodes) * Routine to read in the coordinate of the nodes. Values * are read until End-of-file or too many nodes are entered. * PASSED VARIABLES * Input values: integer*4 max_nodes ! Maximum number of nodes allowed (defined ! in main program) * Returned values integer*4 num_nodes ! Number of nodes input by the user. List ! of node coordinates is read until end-of- ! file (EOF or ^D) is reached. real*8 nodes_xy(2,max_nodes) ! XY coordinates of nodes of the ! of polygon character*(*) units ! Units of the coordinates. Used only for ! output. * LOCAL VARIABLES * Miscellaneous variables that are needed integer*4 ierr ! IOSTAT error returned from reads and ! writes. real*8 input_xy(2) ! User input values of X and Y. If ! values are OK, they are saved in nodes_xy ***** Start: Tell user what they need to enter. Since we will end * list with EOF, the list must be last things entered. write(*,110) 110 format(' What are the units of the coordinates? ') read(*,'(a)',iostat=ierr) units * Report any error on read. if( ierr.ne.0 ) then write(*,120) ierr, units 120 format('**WARNING** IOSTAT error ',i5,' occurred', . ' reading units. Recorded units are ',a,/, . ' Area will still be calculated') end if * MOD TAH Version 2. Change message here to say that warning * is printed if direction changed. Old code has been * commented out. C write(*,200) C200 format(' Input the coordinates of the nodes of the polygon', C . ' in X, Y order, one pair per line, free format',/, C . ' NOTE: No check is made on shape of figure and nodes', C . ' should be entered moving around the figure.',/, C . ' Use ^D (control-D or EOF) to end input',/, C . ' X Y ') write(*,200) 200 format(' Input the coordinates of the nodes of the polygon', . ' in X, Y order, one pair per line, free format',/, . ' NOTE: Warning printed if direction around polygon', . ' changes during area calculation',/, . ' Use ^D (control-D or EOF) to end input',/, . ' X Y ') **** Initialize the number of nodes and start reading until an * error occurrs (EOF will return ierr=-1 and this is the * expected end) num_nodes = 0 ierr = 0 ! Set the error variable so that we enter the ! do while loop. do while ( ierr.eq.0 ) * Read the input coordinates and see if EOF is reached. read(*,*,iostat=ierr) input_xy * If no error occurred in read, save the coordinates of the * node, provided we have not reached maximum allowed. if ( ierr.eq.0 ) then ! No read error. **** Check to see if the number of nodes is too large. * Increment the number of nodes and make sure not too large * (i.e., bigger than max_nodes). if ( num_nodes+1.gt.max_nodes ) then **** Too many nodes are about to be entered. We could * stop at this point or compute using the nodes we * already have. We will do the later so that user * does not need to enter and can complete the nodes * with another program run. write(*,220) max_nodes 220 format('**ERROR** Attempt to input ',i5,' nodes', . ' which exceeds maximum allowed.',/, . ' Computing area with nodes entered. Run', . ' program again with additional nodes and', . ' sum results') ierr = -2 ! Set error to show an error occurred. ! This forces use out of reading loop. else ! Enough space to save num_nodes = num_nodes + 1 nodes_xy(1,num_nodes) = input_xy(1) nodes_xy(2,num_nodes) = input_xy(2) end if * Some error occurred on read. If -1 (EOF) or -2 (to many * entries (already reported) then OK, and * we just continue. If some other error, then report * problem and use just the current nodes. else if ( ierr.ne. -1 .and. ierr.ne.-2 ) then write(*,250) ierr, num_nodes+1 250 format('**ERROR** IOSTAT error ',i5,' occurred reading', . ' input node # ',i5,' coordinates.',/, . ' Truncating list at this point and computing', . ' area') end if end do **** Finally, tell user how nodes were entered write(*,320) num_nodes 320 format(i5,' Node coordinates have been entered.') **** That's all return end CTITLE FORM_TRIANGLES subroutine form_triangle(n, nodes_xy, triangle_vec, . num_nodes) * Routine form the two vectors that make up the triangle * from the node coordinates. * This routine contains a possible stop if n is greater * than the number of nodes (num_nodes) * PASSED VARIABLES * Input values: integer*4 num_nodes ! Maximum number of nodes allowed (defined ! in main program) integer*4 n ! Node number of second side of triangle ! triangle will be formed with 1->(n-1) and ! 1->n real*8 nodes_xy(2,num_nodes) ! XY coordinates of nodes of the ! of polygon * Returned values real*8 triangle_vec(2,2) ! Two vectors that make up the ! current triangle. Lead index is over ! XY and the second index is over side ! 1 and 2. * LOCAL VARIABLES integer*4 j ! Variable for looping over X and Y components **** First check that we are not exceeding the total number * of nodes. This should not happen given way routine is * called but in some future use it might happen if( n.gt.num_nodes) then ! Something is wrong: requested ! node out of range write(*,120) n, num_nodes 120 format(/,'**DISASTER** Triangle requested with node number', . i5,' but only ',i5,' nodes were input',/, . ' Program terminating in form_triangles') * See discussion in main program about internal stops. stop 'Node number too large in form_triangles' end if **** OK, Node number is OK, form the vectors from the apex to * node n and n-1 * Loop over the XY components do j = 1, 2 triangle_vec(j,1) = nodes_xy(j,n-1)-nodes_xy(j,1) triangle_vec(j,2) = nodes_xy(j,n )-nodes_xy(j,1) end do **** Thats all return end CTITLE TRIANGLE_DAREA subroutine triangle_darea(triangle_vec, darea) * Routine to compute the area of the triangle formed * by the two vectors in triangle_vec. A cross product * calculation is used. Since the vectors are 2-d, the * cross product will have only a Z-component and so this * is the only component we compute. * PASSED VARIABLES * Input values: real*8 triangle_vec(2,2) ! Two vectors that make up the ! current triangle. Lead index is over ! XY and the second index is over side ! 1 and 2. * Returned values real*8 darea ! Area of triangle. **** Form the cross product Z-component and divide by two since * cross product a x b = |a||b| sin(theta) in direction normal * to the plane of vectors a and b. theta is the angle between * vectors. * Note: Sign will depend on if we rotate clockwise or * anti-clockwise between vector. We treat this at the end * in the output routine where we use the absolute value. darea = (triangle_vec(1,1)*triangle_vec(2,2) - . triangle_vec(2,1)*triangle_vec(1,2))/2.d0 return end CTITLE OUTPUT_AREA subroutine output_area(num_nodes, area, nodes_xy, units) * Routine to output the final area. Here we take the absolute * value in case the figure is entered anti-clockwise. * PASSED VARIABLES * Input values: integer*4 num_nodes ! Number of nodes input by the user. List ! of node coordinates is read until end-of- ! file (EOF or ^D) is reached. real*8 area ! Total area. This can be positive or ! negative depending whether nodes are ! entered clockwise or anti-clockwise. ! Sign is fixed in output routine. real*8 nodes_xy(2,num_nodes) ! XY coordinates of nodes of the ! of polygon character*(*) units ! Units of the coordinates. Used only for ! output. * LOCAL VARIABLES integer*4 ierr ! IOSTAT error returned from reads and ! writes. integer*4 i, j ! Loop variables for list nodes. **** List the nodes we have used. write(*,120,iostat=ierr) num_nodes, units 120 format(/,' POLY_AREA: The area enclosed by the ',i5, . ' nodes in the following polygon:',/, . ' # X Y ',a) do i = 1, num_nodes write(*,140) i, (nodes_xy(j,i),j=1,2) * Should be careful with format since numbers could be * too large or small for the area. 140 format(i5,1x,2(F12.3,1x)) end do write(*,160,iostat=ierr) abs(area), units 160 format(' Total area is ',f20.6,' square ',a) if( ierr.ne.0 ) then write(*,210) ierr 210 format('**WARNING** IOSTAT error ',i5,' occurred during', . ' output of results') end if **** Thats all return end CTITLE CHECK_DIR subroutine check_dir(n, sign_darea, darea) * Routine to check and warn user if the sign of the * triangle area changes. This means that the direction * around the figure has changed. * PASSED VARIABLES * Input variables integer*4 n ! Node number of second vector in triangle integer*4 sign_darea ! Sign of the area of the first triangle real*8 darea ! Area increment of the triangle formed with ! nodes n and n-1. **** See if first triangle. If it is then set the sign. * For other triangles make sure the same sign. if( n.eq.3 ) then * First triangle. Set the sign of the area. NOTE: * small problem with this code if darea=0 (routine * will set to +ve sign. if( darea.lt.0.d0 ) then sign_darea = -1 else sign_darea = +1 endif else * See if sign has changed. By multiplying the values * we get a positive number of sign has not changed. if( sign_darea*darea.lt.0 ) then write(*,210) n 210 format('**WARNING** Direction around figure ', . ' changes at node number ',i4) end if end if **** Thats all return end