///////////////////////////////////////////////////////////////////////////
//    Copyright (C) 2015 Wizardry and Steamworks - License: GNU GPLv3    //
///////////////////////////////////////////////////////////////////////////
vector wasEllipsePoint(float a, float b) {
    float x = llPow(-1, 1 + (integer) llFrand(2)) * llFrand(a);
    float y = llPow(-1, 1 + (integer) llFrand(2)) * llFrand(b);
    if(llPow(x/a,2) + llPow(y/b,2) <= 1)
        return <x, y, 0>;
    return wasEllipsePoint(a, b);
}

///////////////////////////////////////////////////////////////////////////
//    Copyright (C) 2015 Wizardry and Steamworks - License: GNU GPLv3    //
///////////////////////////////////////////////////////////////////////////
// the angle (in radians) o represents a rotation in the trigonometric 
// sense around the Ox axis where O represents the ellipse center
vector wasRotatedEllipsePoint(float a, float b, float o) {
    float x = llPow(-1, 1 + (integer) llFrand(2)) * llFrand(a);
    float y = llPow(-1, 1 + (integer) llFrand(2)) * llFrand(b);
    if(llPow(x/a,2) + llPow(y/b,2) <= 1)
        return <x*llCos(o) + y*llSin(o), x*llSin(o)-y*llCos(o), 0>;
    return wasRotatedEllipsePoint(a, b, o);
}

///////////////////////////////////////////////////////////////////////////
//    Copyright (C) 2015 Wizardry and Steamworks - License: GNU GPLv3    //
///////////////////////////////////////////////////////////////////////////
// given a set of points, determines the largest segment that can be 
// created between any of those points
list wasLargestSegment(list points) {
    list dists = [];
    list xi = [];
    list xj = [];
    integer i = llGetListLength(points) - 1;
    do {
        integer j = llGetListLength(points) - 1;
        do {        
            if(i == j) jump continue; 
            float d = llVecDist(
                llList2Vector(points, j), 
                llList2Vector(points, i)
            );
            dists += d;
            xi += i;
            xj += j;
@continue;
        } while(--j > -1);
    } while(--i > -1);
    
    i = llGetListLength(dists) - 1;
    float d = llList2Float(dists, i);
    --i;
    integer x = 0;
    integer y = 0;
    do {
        float t = llList2Float(dists, i);
        if(d < t) {
            x = llList2Integer(xi, i);
            y = llList2Integer(xj, i);
            d = t;
        }
    } while(--i > -1);
    
    return [llList2Vector(points, x), llList2Vector(points, y)];
}

///////////////////////////////////////////////////////////////////////////
//    Copyright (C) 2015 Wizardry and Steamworks - License: GNU GPLv3    //
///////////////////////////////////////////////////////////////////////////
// given a set of points, determines the smallest segment that can be 
// created between any of those points
list wasSmallestSegment(list points) {
    list dists = [];
    list xi = [];
    list xj = [];
    integer i = llGetListLength(points) - 1;
    do {
        integer j = llGetListLength(points) - 1;
        do {        
            if(i == j) jump continue; 
            float d = llVecDist(
                llList2Vector(points, j), 
                llList2Vector(points, i)
            );
            dists += d;
            xi += i;
            xj += j;
@continue;
        } while(--j > -1);
    } while(--i > -1);
    
    i = llGetListLength(dists) - 1;
    float d = llList2Float(dists, i);
    --i;
    integer x = 0;
    integer y = 0;
    do {
        float t = llList2Float(dists, i);
        if(d > t) {
            x = llList2Integer(xi, i);
            y = llList2Integer(xj, i);
            d = t;
        }
    } while(--i > -1);
    
    return [llList2Vector(points, x), llList2Vector(points, y)];
}

///////////////////////////////////////////////////////////////////////////
//    Copyright (C) 2015 Wizardry and Steamworks - License: GNU GPLv3    //
///////////////////////////////////////////////////////////////////////////
// takes two points defining a segment and returns the angle created 
// between the Ox axis and the segement (the angle will be positive for the
// trigonometric quadrants I and II and negative for III and IV)
float wasOXSegmentAngle(vector A, vector B) {
    float alpha = llAtan2(
        llVecDist(
            B, 
            <B.x, A.y, 0>
        ) / 
        llVecDist(
            A, 
            <B.x, A.y, 0>
        ), 
        1
    );
    if(B.y < A.y) return -alpha;
    return alpha;
}

///////////////////////////////////////////////////////////////////////////
//    Copyright (C) 2015 Wizardry and Steamworks - License: GNU GPLv3    //
///////////////////////////////////////////////////////////////////////////
// takes a set of points describing a polygon and returns the radiuses of
// the minimal ellipse that encloses all the points of the polygon
list wasMinimalPolygonEllipse(
    list points, 
    vector C, 
    float a, 
    float b, 
    float angle, 
    float increment
) {
    integer i = llGetListLength(polygonPoints)-1;
    do {
        vector p = llList2Vector(polygonPoints, i);
        if(llPow(
            (
                (p.x-C.x)*llCos(a)+(p.y-C.y)*llSin(a)
            )/a, 
            2
        ) + 
            llPow(
            (
                (p.x-C.x)*llSin(a)-(p.y-C.y)*llCos(a)
            )/b, 
            2) <= 1) jump continue;
        b += increment;
        return wasMinimalPolygonEllipse(points, C, a, b, angle, increment);
@continue;
    } while(--i>-1);
    return [a, b];
}

///////////////////////////////////////////////////////////////////////////
//    Copyright (C) 2013 Wizardry and Steamworks - License: GNU GPLv3    //
//    www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html    //
///////////////////////////////////////////////////////////////////////////
integer wasPointInPolygon(vector p, list polygon) {
    integer inside = FALSE;
    integer i = 0;
    integer nvert = llGetListLength(polygon);
    integer j = nvert-1;
    do {
        vector pi = llList2Vector(polygon, i);
        vector pj = llList2Vector(polygon, j);
        if ((pi.y > p.y) != (pj.y > p.y))
            if(p.x < (pj.x - pi.x) * (p.y - pi.y) / (pj.y - pi.y) + pi.x)
                inside = !inside;
        j = i++;
    } while(i<nvert);
    return inside;
}

///////////////////////////////////////////////////////////////////////////
//    Copyright (C) 2015 Wizardry and Steamworks - License: GNU GPLv3    //
///////////////////////////////////////////////////////////////////////////
// takes as input the points describing a polygon and returns a point 
// within the polygon
//
// note that this function allows to dynamically change the polygon points
// and will adjust to the newly specified polygon if the list has changed
// otherwise, consider calculating all the variables within the marked 
// blocks only once and storing them since the performance will improve!
vector wasPolygonPoint(list polygon) {
    ////////////////// CAN BE CALCULATED ONLY ONCE ////////////////////////
    ///////////////// AND STORED IN GLOBAL VARIABLES //////////////////////
    // Get the largest segment.
    list largeSegment = wasLargestSegment(polygonPoints);
    vector A = llList2Vector(largeSegment, 0);
    A.z = 0;
    vector B = llList2Vector(largeSegment, 1);
    B.z = 0;
    // Get the center of the largest segement.
    vector C = <(A.x + B.x)/2, (A.y + B.y)/2, 0>;
    // Get the size of the largest segment.
    float largeSegmentSize = llVecDist(A, B);
    // Get the angle between 0x and the segment
    float a = wasOXSegmentAngle(A, B);
    // Get the smallest segment.
    list smallSegement = wasSmallestSegment(polygonPoints);
    // Get the size of the smallest segment.
    float smallSegmentSize = llVecDist(
        llList2Vector(smallSegement, 0), 
        llList2Vector(smallSegement, 1)
    );
    // Now find the minimal ellipse around the polygon.
    list m = wasMinimalPolygonEllipse(
        polygonPoints,
        C,
        largeSegmentSize,
        smallSegmentSize,
        a,
        0.1
    );
    ////////////////// CAN BE CALCULATED ONLY ONCE ////////////////////////
    ///////////////// AND STORED IN GLOBAL VARIABLES //////////////////////

    vector d = ZERO_VECTOR;
    do {
        d = C + wasRotatedEllipsePoint(llList2Float(m, 0), llList2Float(m, 1), a);
    } while(wasPointInPolygon(d, polygon) == FALSE);
    return d;
}

///////////////////////////////////////////////////////////////////////////
// these points represent the coordinates of the points that define the  //
// contour (permimeter) of the polgyon                                   //
///////////////////////////////////////////////////////////////////////////
list polygonPoints = [
    <163.349121, 171.559189, 3400.894287>,
    <165.56424, 170.422821, 3400.894287>,
    <164.364594, 169.194382, 3400.894287>,
    <165.171875, 167.453903, 3400.894287>,
    <162.345123, 167.61731, 3400.894287>,
    <163.349121, 171.659189, 3400.894287>
];

default
{
    state_entry()
    {
        integer i = 1000;
        do {
            llOwnerSay((string)wasPolygonPoint(polygonPoints));
        } while(--i > -1);
    }
}