* SUBROUTINE SLICE
* Created by M.L. Nagurka and W.C. Hayes, University of Pennsylvania
* from Journal of Biomechanics, Volume 13 pp 59-64, 1980
* adapted by Amy J Stevenson, University of Calgary
*
* This subroutine calculates the area and area inertial properties
* of any cross section containing holes, given x,y co-ordinate  
* points of outer and inner perimeters, either by manual data input
* or data file
*
* **** VARIABLES ******
*
* AREA - CROSS SECTIONAL AREA
* I - X,Y CO-ORDINATE INDEX
* IX - MOMENT OF INERTIA ABOUT X AXIS
* IXBAR - MOMENT OF INERTIA ABOUT X-AXIS TRANSLATED TO CENTROID
* IXBPH - MOMENT OF INERTIA ABOUT TRANSLATED, ROTATED PRICIPLE AXIS
* IXY - PRODUCT OF INERTIA
* IY - MOMENT OF INERTIA ABOUT Y AXIS
* IYBAR - MOMENT OF INERTIA ABOUT Y AXIS TRANSLATED TO CENTROID
* IYBPH - MOMENT OF INERTIA ABOUT TRANSLATED, ROTATED PRINCIPLE  
*         AXIS
* NPER - NUMBER OF PERIMETERS
* NPOINT(IPER) - NUMBER OF POINTS IN SECTION (IPER)
* PHI - ANGLE BETWEEN TRANSLATED AXIS AND PRINCIPLE AXIS
* X(I) - X COORDINATE OF PREVIOUS VERTEX POINT
* X(I+1) - X COORDINATE OF CURRENT VERTEX POINT
* XBAR - X COORDINATE OF CENTROID
* XBARA - X COORDINATE OF CENTROID TIMES AREA
* XVECT(IPER,I) - ARRAY OF X(I) FOR SECTION IPER
* Y(I) - Y COORDINATE OF PREVIOUS VERTEX POINT
* Y(I+1) - Y COORDINATE OF CURRENT VERTEX POINT
* YBAR - Y COORDINATE OF CENTROID
* YBARA - Y COORDINATE OF CENTROID TIMES AREA
* YVECT(IPER,I) - ARRAY OF Y(I) FOR SECTION IPER
*
*
* ----------------------------------------------------------------
*
* SUBROUTINE SLICE REDUCES SECTION PERIMETER INFO
* TO AREA AND AREA INERTIAL PROPERTIES OF SECTION.
*
      SUBROUTINE SLICE(XVECT, YVECT, NPER, maxpnt, NPOINT, 
     *  AREA, XBAR, YBAR, PHI, IXBPH, IYBPH)
      implicit double precision (a-h,o-z)
c
      dimension xvect(maxpnt, 1), yvect(maxpnt, 1)
c
      DIMENSION X(500), Y(500)
      DIMENSION NPOINT(5)
      double precision IX, IY, IXY, IXBAR, IYBAR, IXBYB, IXBPH, iybph
*
* INITIALIZE VARIABLES
          AREA = 0.0
          XBARA=0.0
          YBARA=0.0
          IX=0.0
          IY=0.0
          IXY=0.0
*
          DO 50 IPER= 1, NPER
*
* FILL X,Y VECTORS FOR EACH SECTION PERIMETER
* ORIGINAL LINE N1=NPOINT(IPER)+1
	  N1=NPOINT(IPER)
          DO 20 I=1, N1
          X(I)=XVECT(I,iper)
          Y(I)=YVECT(I,iper)
20        CONTINUE
* NEED TO ADD ANOTHER VALUE BUT NOT INCREMENT NPOINT
* TO SATISFY (I+1) LOOP CONDITION
	  X(N1+1)=X(1)
	  Y(N1+1)=Y(1)
*
* CALCULATIONS
          DO 40 I=1, NPOINT(IPER)
          AREA=AREA-((Y(I+1)-Y(I))*(X(I+1)+X(I)))/2
          XBARA=XBARA-(((Y(I+1)-Y(I))/8)*((X(I+1)+X(I))**2
     1    + ((X(I+1)-X(I))**2)/3))
          YBARA=YBARA+(((X(I+1)-X(I))/8)*((Y(I+1)+Y(I))**2
     1    + ((Y(I+1)-Y(I))**2)/3))
          IX=IX+((X(I+1)-X(I))*(Y(I+1)+Y(I))/24)*
     1    ((Y(I+1)+Y(I))**2 + (Y(I+1)-Y(I))**2)
          IY=IY-((Y(I+1)-Y(I))*(X(I+1)+X(I))/24)*
     1    ((X(I+1)+X(I))**2+(X(I+1)-X(I))**2)
          IF ((X(I+1)-X(I)) .EQ. 0) GOTO 40
          IXY=IXY+(((Y(I+1)-Y(I))**2)*(X(I+1)+X(I))*
     1    ((X(I+1))**2+(X(I))**2)/8+
     2    (Y(I+1)-Y(I))*(X(I+1)*Y(I)-X(I)*Y(I+1))*
     3    ((X(I+1))**2+X(I+1)*X(I)+(X(I))**2)/3+
     4    ((X(I+1)*Y(I)-X(I)*Y(I+1))**2)*
     5    (X(I+1)+X(I))/4)/(X(I+1)-X(I))
40        CONTINUE
50        CONTINUE
*
* EXTRA CALCULATIONS
          XBAR=XBARA/AREA
          YBAR=YBARA/AREA
***
          IXBAR=IX-AREA*YBAR**2
          IYBAR=IY-AREA*XBAR**2
          IXBYB=IXY-AREA*XBAR*YBAR
          BIG=ABS(IXBAR)
          IF (IYBAR .GT. IXBAR) BIG=ABS(IYBAR)
	IF ((ABS(IXBYB)) .LE. 1d-12) THEN
	  PHI=0.0
	ELSE
	  IF ((ABS(IXBAR-IYBAR)/BIG) .LE. 1d-5) THEN
		PHI=-0.78539816
	  ELSE
		PHI=(ATAN(-2d0*IXBYB/(IXBAR-IYBAR)))/2d0
	  ENDIF
	ENDIF
          IXBPH=IXBAR*(COS(PHI))**2+IYBAR*(SIN(PHI))**2-
     1    IXBYB*SIN(2*PHI)
          IYBPH=IYBAR*(COS(PHI))**2+IXBAR*(SIN(PHI))**2+
     1    IXBYB*SIN(2*PHI)
c        print *,ix,iy,ixy
c        print *,ixbar,iybar,ixbyb
c        print *,ixbph, iybph, phi
*
* CALCULATIONS ABOUT ARBITRARY AXIS
* NOT NEEDED FOR OUR PURPOSES, SECTION ELIMINATED
          RETURN
          END