SimpleGMap2::SimpleGMap2()
{
	 position = myMap.addAttribute<VEC3, VERTEX>("position");

     Dart d = Algo::Modelisation::createTetrahedron<PFP>(myMap);
     position[d] = VEC3(0,0,0);
     position[myMap.phi1(d)] = VEC3(10,0,15);
     position[myMap.phi_1(d)] = VEC3(10,20,15);
     position[myMap.phi_1(myMap.phi2(d))] = VEC3(0,0,30);

     VEC3 mid = (position[d] + position[myMap.phi1(d)]) / 2.0f;
     myMap.cutEdge(d);
     position[myMap.phi1(d)] = mid;

     Algo::Modelisation::Polyhedron<PFP> poly(myMap, position);

     d = poly.cylinder_topo(5, 1, false, false);

     poly.embedCylinder(10, 10, 5);

     d = myMap.phi1(d);
     Dart dd = myMap.beta2(d);
     myMap.unsewFaces(d);
     myMap.sewFaces(d, dd);

     position[d][1] += 3.0f;
}

 void read_polygon_xyz(boost::ptr_vector<geometry_type> & paths)
 {
     int num_rings = read_integer();
     if (num_rings > 0)
     {
         std::unique_ptr<geometry_type> poly(new geometry_type(geometry_type::types::Polygon));
         for (int i = 0; i < num_rings; ++i)
         {
             int num_points = read_integer();
             if (num_points > 0)
             {
                 CoordinateArray ar(num_points);
                 read_coords_xyz(ar);
                 poly->move_to(ar[0].x, ar[0].y);
                 for (int j = 1; j < num_points; ++j)
                 {
                     poly->line_to(ar[j].x, ar[j].y);
                 }
                 poly->close_path();
             }
         }
         if (poly->size() > 2) // ignore if polygon has less than 3 vertices
             paths.push_back(poly.release());
     }
 }

/*
 * drawpoly
 *
 *	draw some polygons
 */
void
drawpoly()
{
	float	vec[3];
	short	val;

	color(YELLOW);

	/*
	 * Draw a polygon using poly, parray is our array of
	 * points and 4 is the number of points in it.
	 */
	poly(4L, parray);

	color(GREEN);

	/*
	 * Draw a 5 sided figure by using bgnpolygon, v3d, and endpolygon
	 */
	polymode(PYM_LINE);
	bgnpolygon();
		vec[0] = 0.0;
		vec[1] = 0.0;
		vec[2] = 0.0;
		v3f(vec);
		vec[0] = 3.0;
		vec[1] = 0.0;
		vec[2] = 0.0;
		v3f(vec);
		vec[0] = 3.0;
		vec[1] = 4.0;
		vec[2] = 0.0;
		v3f(vec);
		vec[0] = -1.0;
		vec[1] = 5.0;
		vec[2] = 0.0;
		v3f(vec);
		vec[0] = -2.0;
		vec[1] = 2.0;
		vec[2] = 0.0;
		v3f(vec);
	endpolygon();

	color(MAGENTA);

	/*
	 * draw a sector representing a 1/4 circle
	 */
	arc(1.5, -7.0, 3.0, 0, 900);

	move2(1.5, -7.0);
	draw2(1.5, -4.0);

	move2(1.5, -7.0);
	draw2(4.5, -7.0);

	qread(&val);
}

bool monster::make_fungus()
{
    char polypick = 0;
    std::string tid = type->id;
    if (type->in_species("FUNGUS")) { // No friendly-fungalizing ;-)
        return true;
    }
    if (tid == "mon_ant" || tid == "mon_ant_soldier" || tid == "mon_ant_queen" || tid == "mon_fly" ||
      tid == "mon_bee" || tid == "mon_dermatik") {
        polypick = 1;
    } else if (tid == "mon_zombie" || tid == "mon_zombie_shrieker" || tid == "mon_zombie_electric" ||
      tid == "mon_zombie_spitter" || tid == "mon_zombie_dog" || tid == "mon_zombie_brute" ||
      tid == "mon_zombie_hulk" || tid == "mon_zombie_soldier" || tid == "mon_zombie_tough" ||
      tid == "mon_zombie_scientist" || tid == "mon_zombie_hunter" || tid == "mon_zombie_child"||
      tid == "mon_zombie_bio_op" || tid == "mon_zombie_survivor" || tid == "mon_zombie_fireman" ||
      tid == "mon_zombie_cop" || tid == "mon_zombie_fat" || tid == "mon_zombie_rot" ||
      tid == "mon_zombie_swimmer" || tid == "mon_zombie_grabber" || tid == "mon_zombie_technician" ||
      tid == "mon_zombie_brute_shocker") {
        polypick = 2; // Necro and Master have enough Goo to resist conversion.
        // Firefighter, hazmat, and scarred/beekeeper have the PPG on.
    } else if (tid == "mon_zombie_necro" || tid == "mon_zombie_master" || tid == "mon_zombie_fireman" ||
      tid == "mon_zombie_hazmat" || tid == "mon_beekeeper") {
        return true;
    } else if (tid == "mon_boomer" || tid == "mon_zombie_gasbag" || tid == "mon_zombie_smoker") {
        polypick = 3;
    } else if (tid == "mon_triffid" || tid == "mon_triffid_young" || tid == "mon_triffid_queen") {
        polypick = 4;
    }
    switch (polypick) {
        case 1: // bugs, why do they all turn into fungal ants?
            poly(GetMType("mon_ant_fungus"));
            return true;
        case 2: // zombies, non-boomer
            poly(GetMType("mon_zombie_fungus"));
            return true;
        case 3:
            poly(GetMType("mon_boomer_fungus"));
            return true;
        case 4:
            poly(GetMType("mon_fungaloid"));
            return true;
        default:
            return false;
    }
}

Real DistLine3Circle3<Real>::GetSquared ()
{
    Vector3<Real> diff = mLine->Origin - mCircle->Center;
    Real diffSqrLen = diff.SquaredLength();
    Real MdM = mLine->Direction.SquaredLength();
    Real DdM = diff.Dot(mLine->Direction);
    Real NdM = mCircle->Normal.Dot(mLine->Direction);
    Real DdN = diff.Dot(mCircle->Normal);

    Real a0 = DdM;
    Real a1 = MdM;
    Real b0 = DdM - NdM*DdN;
    Real b1 = MdM - NdM*NdM;
    Real c0 = diffSqrLen - DdN*DdN;
    Real c1 = b0;
    Real c2 = b1;
    Real rsqr = mCircle->Radius*mCircle->Radius;

    Real a0sqr = a0*a0;
    Real a1sqr = a1*a1;
    Real twoA0A1 = ((Real)2)*a0*a1;
    Real b0sqr = b0*b0;
    Real b1Sqr = b1*b1;
    Real twoB0B1 = ((Real)2)*b0*b1;
    Real twoC1 = ((Real)2)*c1;

    // The minimum point B+t*M occurs when t is a root of the quartic
    // equation whose coefficients are defined below.
    Polynomial1<Real> poly(4);
    poly[0] = a0sqr*c0 - b0sqr*rsqr;
    poly[1] = twoA0A1*c0 + a0sqr*twoC1 - twoB0B1*rsqr;
    poly[2] = a1sqr*c0 + twoA0A1*twoC1 + a0sqr*c2 - b1Sqr*rsqr;
    poly[3] = a1sqr*twoC1 + twoA0A1*c2;
    poly[4] = a1sqr*c2;

    PolynomialRoots<Real> polyroots(Math<Real>::ZERO_TOLERANCE);
    polyroots.FindB(poly, 6);
    int count = polyroots.GetCount();
    const Real* roots = polyroots.GetRoots();

    Real minSqrDist = Math<Real>::MAX_REAL;
    for (int i = 0; i < count; ++i)
    {
        // Compute distance from P(t) to circle.
        Vector3<Real> P = mLine->Origin + roots[i]*mLine->Direction;
        DistPoint3Circle3<Real> query(P, *mCircle);
        Real sqrDist = query.GetSquared();
        if (sqrDist < minSqrDist)
        {
            minSqrDist = sqrDist;
            mClosestPoint0 = query.GetClosestPoint0();
            mClosestPoint1 = query.GetClosestPoint1();
        }
    }

    return minSqrDist;
}

poly makePoly( vecRef pos, scalar radius, int segments )
{
	vertexList verts;
	for (int i=0; i<segments; i++) {
		scalar a = TWO_PI - i/scalar(segments) * TWO_PI; // ccw winding
		verts.push_back( pos + radius * vec::polar(a) );
	}
	return poly( verts );
}

// -----------------------------------------------------------------------------
double
dNtcResToTemp(double dR, double dCoeff[])
{
  double ti;

  ti = poly(log(dR), 3, dCoeff);
  ti = 1.0 / ti + TABS;
  return ti;
}

int main (void)

{
    
    static int i;
    
    int n;
    
    printf("Please enter the maximum degree of the polynomial:\n");
    
    scanf("%d", &n);
    
    if(n>0){                            //if power > 0 do followings
    
    printf("Please enter the coefficients:\n");

    
    
    for(i=0; i<=n; i++)                 //saving array a
    {
        
        scanf("%d", &a[i]);
        
    }
    
    
    for(i=0; i<=n; i++)                 //saving array b
    {
        
        b[9-i]=a[i];
        
    }


    
    printf("The polynomial is ");
    
    poly(n);
    
    printf("\nIts reciprocal is ");
    
    print_reci(n);
    
    printf("\nThe polynomial is");
    
    compare(n);
        
    }
    
    else                           //if power < = 0 do nothing
    {
      
    }
    
    
}

GeomShape *gpcPolygonToGeomShape( gpc_polygon *gpcPolygon ) {
	GeomShape *geomShape = new GeomShape;

	for( int i=0; i<gpcPolygon->num_contours; i++ ) {
		GeomPoly poly( gpcPolygon->contour[i].num_vertices, (DVec2*)gpcPolygon->contour[i].vertex );
		geomShape->addPoly( &poly, gpcPolygon->hole[i] );
	}

	return geomShape;
}

/*
 * Генерация и печать битовой последовательности длины `len',
 * используя сдвиговый регистр длины `N' с обратной связью `tap'.
 */
void makeseq (int tap, int len)
{
	unsigned char seq [2*P];
	int mem, i;

	mem = 0;
	for (i=0; i<P; ++i)
		poly (&mem, tap, 1);
	for (i=0; i<2*P; ++i)
		seq [i] = poly (&mem, tap, 1);

	printf ("\nx^%d + x^%d + 1 =\n", N, tap);
	seqprint (seq, len);
/*	hexprint (seq, len);*/

	for (i=0; i<2*P; ++i)
		seq [i] ^= 1;
	seqprint (seq, len);
}

void test_polygon_Bbox()
{
    Point_E3d P1(-5,-4,-3);
    Point_E3d P2( 5, 4, 3);

    Point_E3d P[] = { P1, P2, Point_E3d(-2, 2, 2) };
    Polygon_E3d poly(P, P+3);
    Bbox_E3d box = poly.get_Bbox();
    assert( box == Bbox_E3d(P1, P2) );
}

Polynomial1<Real> IntrEllipse2Ellipse2<Real>::GetQuartic (
    const Ellipse2<Real>& ellipse0, const Ellipse2<Real>& ellipse1)
{
    Real p0[6], p1[6];
    ellipse0.ToCoefficients(p0);
    ellipse1.ToCoefficients(p1);

    // The polynomials are
    //   P0 = a0 + a1*x + a2*y + a3*x^2 + a4*x*y + a5*y^2
    //      = (a0 + a2*y + a5*y^2) + (a1 + a4*y)*x + (a3)*x^2
    //      = u0(y) + u1(y)*x + u2(y)*x^2
    //   P1 = b0 + b1*x + b2*y + b3*x^2 + b4*x*y + b5*y^2
    //      = (b0 + b2*y + b5*y^2) + (b1 + b4*y)*x + (b3)*x^2
    //      = v0(y) + v1(y)*x + v2(y)*x^2
    // The Bezout determinant eliminates the variable x when solving the
    // equations P0(x,y) = 0 and P1(x,y) = 0.  We have
    //   0 = P0 = u0 + u1*x + u2*x^2
    //   0 = P1 = v0 + v1*x + v2*x^2
    //   0 = v2*P0 - u2*P1 = (u0*v2 - u2*v0) + (u1*v2 - u2*v1)*x
    //   0 = v1*P0 - u1*P1 = (u0*v1 - u1*v0) + (u2*v1 - u1*v2)*x^2
    // Solve the equation 0 = v2*P0-u2*P1 for x and substitute in the other
    // equation and simplify to
    //   Q(y) = (u0*v1-v1*u0)*(u1*v2-u2*v1) - (u0*v2-u2*v0)^2 = 0
    //        = c0 + c1*y + c2*y^2 + c3*y^3 + c4*y^4
    // Define dij = ai*bj - aj*bi for various indices i and j.  For example,
    // d01 = a0*b1-b1*a0.  The coefficients of Q(y) are
    //   c0 = d01*d13 - d30^2
    //   c1 = d01*d43 + (d04+d21)*d13 - 2*d30*d32
    //   c2 = (d04+d21)*d43 + (d24+d51)*d13 - 2*d30*d35 - d32^2
    //   c3 = (d24+d51)*d43 + d54*d13 - 2*d32*d35
    //   c4 = d54*d43 - d35^2

    Real d01 = p0[0]*p1[1] - p0[1]*p1[0];
    Real d04 = p0[0]*p1[4] - p0[4]*p1[0];
    Real d13 = p0[1]*p1[3] - p0[3]*p1[1];
    Real d21 = p0[2]*p1[1] - p0[1]*p1[2];
    Real d24 = p0[2]*p1[4] - p0[4]*p1[2];
    Real d30 = p0[3]*p1[0] - p0[0]*p1[3];
    Real d32 = p0[3]*p1[2] - p0[2]*p1[3];
    Real d35 = p0[3]*p1[5] - p0[5]*p1[3];
    Real d43 = p0[4]*p1[3] - p0[3]*p1[4];
    Real d51 = p0[5]*p1[1] - p0[1]*p1[5];
    Real d54 = p0[5]*p1[4] - p0[4]*p1[5];
    Real d04p21 = d04 + d21;
    Real d24p51 = d24 + d51;

    Polynomial1<Real> poly(4);
    poly[0] = d01*d13 - d30*d30;
    poly[1] = d01*d43 + d04p21*d13 - ((Real)2)*d30*d32;
    poly[2] = d04p21*d43 + d24p51*d13 - ((Real)2)*d30*d35 - d32*d32;
    poly[3] = d24p51*d43 + d54*d13 - ((Real)2)*d32*d35;
    poly[4] = d54*d43 - d35*d35;

    return poly;
}

void AverageAreaGridDensityProvider::initialize(const real proscenium[4], real sizeFactor)
{
	float prosceniumWidth = (proscenium[1] - proscenium[0]);
	float prosceniumHeight = (proscenium[3] - proscenium[2]);

	real cellArea = 0.0;
	unsigned numFaces = 0;
	for (source.begin(); source.isValid(); source.next()) {
		Polygon3r& poly(source.getGridSpacePolygon());
		Vec3r min, max;
		poly.getBBox(min, max);
		cellArea += (max[0] - min[0]) * (max[1] - min[1]);
		++numFaces;
	}
	if (G.debug & G_DEBUG_FREESTYLE) {
		cout << "Total area: " << cellArea << ".  Number of faces: " << numFaces << "." << endl;
	}
	cellArea /= numFaces;
	cellArea *= sizeFactor;
	if (G.debug & G_DEBUG_FREESTYLE) {
		cout << "Building grid with average area " << cellArea << endl;
	}

	_cellSize = sqrt(cellArea);
	unsigned maxCells = 931; // * 1.1 = 1024
	if (std::max(prosceniumWidth, prosceniumHeight) / _cellSize > maxCells) {
		if (G.debug & G_DEBUG_FREESTYLE) {
			cout << "Scene-dependent cell size (" << _cellSize << " square) is too small." << endl;
		}
		_cellSize = std::max(prosceniumWidth, prosceniumHeight) / maxCells;
	}
	// Now we know how many cells make each side of our grid
	_cellsX = ceil(prosceniumWidth / _cellSize);
	_cellsY = ceil(prosceniumHeight / _cellSize);
	if (G.debug & G_DEBUG_FREESTYLE) {
		cout << _cellsX << "x" << _cellsY << " cells of size " << _cellSize << " square." << endl;
	}

	// Make sure the grid exceeds the proscenium by a small amount
	float safetyZone = 0.1f;
	if (_cellsX * _cellSize < prosceniumWidth * (1.0 + safetyZone)) {
		_cellsX = ceil(prosceniumWidth * (1.0 + safetyZone) / _cellSize);
	}
	if (_cellsY * _cellSize < prosceniumHeight * (1.0 + safetyZone)) {
		_cellsY = ceil(prosceniumHeight * (1.0 + safetyZone) / _cellSize);
	}
	if (G.debug & G_DEBUG_FREESTYLE) {
		cout << _cellsX << "x" << _cellsY << " cells of size " << _cellSize << " square." << endl;
	}

	// Find grid origin
	_cellOrigin[0] = ((proscenium[0] + proscenium[1]) / 2.0) - (_cellsX / 2.0) * _cellSize;
	_cellOrigin[1] = ((proscenium[2] + proscenium[3]) / 2.0) - (_cellsY / 2.0) * _cellSize;
}

void
PolynomialFitTest::sample()
{
  PolynomialFit poly( *_x, *_y, 2 );
  poly.generate();

  CPPUNIT_ASSERT( std::abs(poly.sample( -2. ) - 4) < _tol );
  CPPUNIT_ASSERT( std::abs(poly.sample( -1. ) - 1) < _tol );
  CPPUNIT_ASSERT( std::abs(poly.sample(  0. ) - 0) < _tol );
  CPPUNIT_ASSERT( std::abs(poly.sample(  1. ) - 1) < _tol );
  CPPUNIT_ASSERT( std::abs(poly.sample(  2. ) - 4) < _tol );
}

int main(int argc, char const *argv[]) {

  int i;
  double a[MAXN];
  for (i = 0; i < MAXN; i++) a[i]=(double)i;
  start=clock();
  poly(MAXN-1,a,1.1);
  stop=clock();
  duration = ((double)(stop-start))/CLK_TCK;
  printf("%f\n",(double)duration );
  return 0;
}

int drawpolys(Polys *p)
{
    int i;
    int l;

    for (i=0; i<p->polyc;i++){
        l = *(p->polyl+i);
        poly(*(p->x+i), *(p->y+i), l);
    }

    return 0;
}

void El_Window::hach
    (
          std::vector<Pt2dr> & aVPoly,
          Pt2dr           dir,
          REAL            esp,
          Line_St         lst
    )
{
    ElFifo<Pt2dr> poly(aVPoly,true);
    HWind   HW(*this,lst);
    HPoly(HW,poly,dir,esp);
}

void GDSin::Gds2Ted::convert(GDSin::GdsStructure* src, laydata::tdtcell* dst)
{
   //@FIXME!!! ged rid of this silly loop! GDS database is layer oriented now!
   for(int2b laynum = 1 ; laynum < GDS_MAX_LAYER ; laynum++)
   {
      if (src->allLay(laynum))
      {// layers
         GDSin::GdsData *wd = src->fDataAt(laynum);
         while( wd )
         {
            word tdtlaynum;
            if (_theLayMap.getTdtLay(tdtlaynum, laynum, wd->singleType()) )
            {
               // convert only if layer/data type pair is defined
               laydata::tdtlayer* dwl = static_cast<laydata::tdtlayer*>(dst->securelayer(tdtlaynum));
               switch( wd->gdsDataType() )
               {
                  case      gds_BOX: box (static_cast<GDSin::GdsBox*>(wd)     , dwl, laynum);  break;
                  case gds_BOUNDARY: poly(static_cast<GDSin::GdsPolygon*>(wd) , dwl, laynum);  break;
                  case     gds_PATH: wire(static_cast<GDSin::GDSpath*>(wd)    , dwl, laynum);  break;
                  case     gds_TEXT: text(static_cast<GDSin::GdsText*>(wd)    , dwl);  break;
                  default: assert(false); /*Error - unexpected type*/
               }
            }
            //else
            //{
            //   // The message below could be noughty - so put it in place ONLY if the data type
            //   // is different from default
            //   std::ostringstream ost;
            //   ost << "Layer: " << laynum << "; data type: " << wd->singleType() <<
            //         "; found in GDS database, but not in the conversion map.";
            //   tell_log(console::MT_INFO,ost.str());
            //}
            wd = wd->last();
         }
      }
   }
   if (src->allLay(0))
   {// references
      GDSin::GdsData *wd = src->fDataAt(0);
      while( wd )
      {
         switch( wd->gdsDataType() )
         {
            case gds_SREF: ref (static_cast<GDSin::GdsRef*>(wd)   , dst);  break;
            case gds_AREF: aref(static_cast<GDSin::GdsARef*>(wd)  , dst);  break;
                  default: assert(false); /*Error - unexpected type*/
         }
         wd = wd->last();
      }
   }
   dst->resort();
}

// Constructors
View::View() 
  : iOrigin(0.0, 0.0)
{
  iColor[0] = 0.0;
  iColor[1] = 1.0;  
  iColor[2] = 0.0;    
  Polygon2 poly(gPoints, gPoints+4);
  setCollisionPolygon(poly);
 
  // cout << hex << "0x" << (int)this << " view created" << endl;  // DEBUG
  setTag(gNextTag++); // DEBUG
}

int main(int argc, char** argv)
{

  bem2d::freqtype k={10,0};
        int n=3*k.re;




        std::vector<bem2d::Point> triangle;
        triangle.push_back(bem2d::Point(0,0));
        triangle.push_back(bem2d::Point(1,0));
        triangle.push_back(bem2d::Point(.5,.5*sqrt(3)));
        bem2d::Polygon poly(triangle,n);
        bem2d::pGeometry pgeom=poly.GetGeometry();



        bem2d::PolBasis::AddBasis(0,pgeom);


        bem2d::Point direction=bem2d::normalize(bem2d::Point(1,1));
        bem2d::PlaneWave pw(direction,k);
        bem2d::CombinedPlaneWave cpw(direction,k);

        bem2d::SoundSoftScattering<bem2d::PlaneWave,bem2d::CombinedPlaneWave> soundsoft(pgeom,k,pw,cpw);
        soundsoft.SetQuadOption(5,5,0.15);
        soundsoft.set_polygons(pgeom);
        soundsoft.set_plotInterior();



        std::cout << "Discretizing with k=" << bem2d::complex(k.re,k.im) << " and n=" << pgeom->size() << std::endl;
        clock_t start, finish;
        double time;

        start=clock();
        soundsoft.Discretize();
        finish=clock();
        time=(double(finish)-double(start))/CLOCKS_PER_SEC;
        std::cout << "Computing time (minutes): " << time/60 << std::endl;



        int xpts=300;
        int ypts=300;
        bem2d::pOutputHandler pout(new bem2d::GplotOutput(xpts,ypts,-4,8,-4,8,"wedge"));
        soundsoft.SetOutput(pout);
        soundsoft.WriteAll();


}