/*
 *  =================================================================== 
 *  polygon_properties.c : compute engineering properties of polygons.
 *                                                                     
 *  Copyright (C) 1994-96 by Mark Austin and David Mazzoni.
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software on any computer system,
 *  and to redistribute it freely, subject to the following restrictions:
 * 
 *  1. The authors are not responsible for the consequences of use of
 *     this software, even if they arise from defects in the software.
 *  2. The origin of this software must not be misrepresented, either
 *     by explicit claim or by omission.
 *  3. Altered versions must be plainly marked as such, and must not
 *     be misrepresented as being the original software.
 *  4. This notice is to remain intact.
 *                                                                    
 *  Written by: Mark Austin                                  March 1994
 *  ===================================================================
 */

#include <stdio.h>

#include <stddef.h>
#include <stdlib.h>
#include <stdarg.h>

#include "polygon.h"
#include "miscellaneous.h"

/* 
 *  ==================================================================================== 
 *  Compute engineering Property of exterior polygon and interior loops
 *  
 *  Input  :  float (*EngPropertyPtr)( float, float, float, float ) -- pointer to
 *            function containing four arguments of type float.
 *         :  POLYGON *spPolygon  -- Pointer to polygon.
 *  Output :  float -- Engineering poperty.
 *  ==================================================================================== 
 */ 

float
polyProperty( float (*EngPropertyPtr)( float, float, float, float ), POLYGON *spPolygon ) {
POLYGON *spPoly1, *spPoly2;
float dProperty = 0.0;

  /* [a] : Walk along list of polygons and compute total area */

  for(spPoly1 = spPolygon; spPoly1 != NULL; spPoly1 = spPoly1->spNext) {

     /* [a.1] : Compute area of exterior polygon */

     dProperty += polyPropertyOfLoop( EngPropertyPtr, spPoly1->spEdge );

     /* [a.2] : Subtract area of holes in polygon */

     spPoly2 = spPoly1->uPolygon.spChild;
     while(spPoly2 != NULL ) {
           dProperty += polyPropertyOfLoop( EngPropertyPtr, spPoly2->spEdge );
           spPoly2 = spPoly2->spNext;
     }
  }

  return (dProperty);
}

float
polyPropertyOfLoop( float (*EngPropertyPtr)( float, float, float, float ), EDGE *spFirst ) {
EDGE *spTemp;
float fX1, fY1, fX2, fY2;
float fX1save,   fY1save;
float dProperty = 0.0;

   spTemp = spFirst;
   fX1save = fX1 = spFirst->spVertex->fXcoord;
   fY1save = fY1 = spFirst->spVertex->fYcoord;

   while( spTemp->spNext != spFirst ) {
       spTemp = spTemp->spNext;

       fX2 = spTemp->spVertex->fXcoord;
       fY2 = spTemp->spVertex->fYcoord;

       dProperty += (*EngPropertyPtr)( fX1, fY1, fX2, fY2);

       fX1 = fX2; fY1 = fY2;
   }

   dProperty += (*EngPropertyPtr)( fX1, fY1, fX1save, fY1save);

   return (dProperty);
}

/* 
 *  ==================================================== 
 *  Compute Engineering Property of Trapezoid Components 
 *  ==================================================== 
 */ 

/* Area of Trapezoid */

float polyPropertyArea( float fX1, float fY1, float fX2, float fY2 ) {
float dArea;

     dArea = 0.5*(fY1 + fY2)*(fX2 - fX1);
     return (dArea);
}

/* I_xx of Trapezoid */

float polyPropertyIxx( float fX1, float fY1, float fX2, float fY2 ) {
float dI_xx;

     dI_xx = (1.0/12.0)*(fY1 + fY2)*(fY1*fY1 + fY2*fY2)*(fX2 - fX1);
     return (dI_xx);
}

/* I_yy of Trapezoid */

float polyPropertyIyy( float fX1, float fY1, float fX2, float fY2 ) {
float dI_yy;

     dI_yy = (1.0/12.0)*(fX2*fX2*(3*fY2 + fY1) + 2*fX1*fX2*(fY1 + fY2) +
                         fX1*fX1*(3*fY1 + fY2))*(fX2 - fX1);
     return (dI_yy);
}

/* I_xy of Trapezoid */

float polyPropertyIxy( float fX1, float fY1, float fX2, float fY2 ) {
float dI_xy;

     dI_xy = (1.0/12.0)*(fY2*fY2*(3*fX2 + fX1) + 2*fY1*fY2*(fX1 + fX2) +
                         fY1*fY1*(3*fX1 + fX2))*(fX2 - fX1);
     return (dI_xy);
}