float unwrapphase(float ph1, float ph2, float ph3, int indx )
{
    // unwrap phase
    // 6 sections:  (biggest to smallest)
    //     MAX  MED  MIN
    // 0 = ph1, ph3, ph2
    // 1 = ph3, ph1, ph2
    // 2 = ph3, ph2, ph1
    // 3 = ph2, ph3, ph1
    // 4 = ph2, ph1, ph3
    // 5 = ph1, ph2, ph3
    int N;
    float minph = min(ph1,min(ph2,ph3));
    float maxph = max(ph1,max(ph2,ph3));
    float medph;
    if (maxph==ph1 && minph==ph2)
    {
        N=0;
        medph=ph3;
    } 
    else if (maxph==ph3 && minph==ph2)
    {
        N=1;
        medph=ph1;
    } 
    else if (maxph==ph3 && minph==ph1)
    {
        N=2;
        medph=ph2;
    } 
    else  if (maxph==ph2 && minph==ph1)
    {
        N=3;
        medph=ph3;
    } 
    else if (maxph==ph2 && minph==ph3)
    {
        N=4;
        medph=ph1;
    } 
    else // max==ph1, min==ph3	
    {
        N=5;
        medph=ph2;
    }
    
    float r=(medph-minph)/(maxph-minph+0.0001f);

    // correct for sinwavey-ness
    if (r < 1) 
        r = -0.5 - 0.5 / PI * ( 6.0 * atan( scalefac * ( 1+r ) / ( r-1 ) ) );

    if ( (N & 1) == 1 ) // every uneven 
        r = -r;

    r = (2 * ((N+1) / 2) + r);

    gr.pixels[indx] = color( (ph1+ph2+ph3) / 3.0 * 255 );

    return r/6.0;
}

int v4_to_u32 ( v4 vec )
{
    int r=int(vec.x*255.f+0.5f);
    int g=int(vec.y*255.f+0.5f);
    int b=int(vec.z*255.f+0.5f);
    int a=int(vec.w*255.f+0.5f);
    if (r<0) r=0; 
    else if (r>255) r=255;
    if (g<0) g=0; 
    else if (g>255) g=255;
    if (b<0) b=0; 
    else if (b>255) b=255;
    if (a<0) a=0; 
    else if (a>255) a=255;
    return (r<<0) | (g<<8) | (b<<16) | (a<<24);
}

class v4 {
    float x,y,z,w;

    v4 (float _x, float _y, float _z, float _w)
    {
        x = _x; 
        y = _y; 
        z = _z; 
        w = _w;
    }   
}

float frac(float f)
{
    return f-floor(f);
}

//#define IDX(x,y) _isconf[y][x],_isconf[y][x]
int[] IDX (int x, int y)
{
    return new int[]{ _isconf[y][x],_isconf[y][x] };
}

void unwrap( float r, int x, int y )
{
    float myr = _wrapphase[y][x] - frac(r);
    //if (myr>0.5 && myr<0.75) return;
    //if (myr<-0.5 && myr>-0.75) return;
    if (myr >  0.5f)
        myr -= 1.f;
        
    if (myr < -0.5f)
        myr += 1.f;
        
    _wrapphase[y][x] = myr + r;
    _unwraporder[y][x] = idxx++;

    pix.add( new int[]{x,y} );
    _var[y][x]=-1;
}