#include <math.h>
#include "mytypes.h"
#include "revolve.h"       /* need to get scrnpair from here */
#include "mapstuff.h"
#include "menuexp.h"

#ifdef TEST
#undef DebugOn
#define DebugOn 1
#endif TEST


#define DEBUG
/*
 * this section of code derived from:
 * "The essential algorithms of ray tracing" by Eric Haines
 * presented in Sigraph proceedings on Ray Tracing
 * my major change has been to simplify it for two dimensions
 */

typedef struct {
    float x, y;
} Vector;

static float DotVector(a,b)
    Vector *a, *b;
{
    return( a->x * b->x + a->y * b->y);
}

static void DivVector(in, scale, out)
    Vector *in, *out;
    float scale;
{
    if ( fabs(scale) < SingleTinyVal ) {
        out->x = SingleLargeVal;
        out->y = SingleLargeVal;
    }
    else {
        out->x = in->x / scale;
        out->y = in->y / scale;
    }
}



static Vector Na, Nb, Nc;
static float Du0, Du1, Du2,
             Dv0, Dv1, Dv2;
static Vector Qux, Quy,
              Qvx, Qvy;
static float Dux, Duy,
             Dvx, Dvy;
static bool IsQuadu, IsQuadv;

void CalcMapConsts(vp)
    register ScrnPair *vp;
#define p00 vp[0]
#define p01 vp[1]
#define p11 vp[2]
#define p10 vp[3]
{
    Vector Pa, Pb, Pc, Pd;

    Pa.x = p00.x - p10.x + p11.x - p01.x;
    Pa.y = p00.y - p10.y + p11.y - p01.y;

    Pb.x = p10.x - p00.x;
    Pb.y = p10.y - p00.y;

    Pc.x = p01.x - p00.x;
    Pc.y = p01.y - p00.y;

    Pd.x = p00.x;
    Pd.y = p00.y;

    Na.x = Pa.y; Na.y = -Pa.x;
    Nc.x = Pc.y; Nc.y = -Pc.x;
    Nb.x = Pb.y; Nb.y = -Pb.x;

    Du0 = DotVector(&Nc, &Pd);
    Du1 = DotVector(&Na, &Pd) + DotVector(&Nc, &Pb);
    Du2 = DotVector( &Na, &Pb);

    if( fabs( Du2 ) > SingleTinyVal ) {
        float TwoDu2;

        IsQuadu = true;
        TwoDu2 = 2.0 * Du2;
        DivVector( &Na, TwoDu2, &Qux );
        DivVector( &Nc, -Du2, &Quy );
        Duy = Du0/Du2;
        Dux = -Du1/TwoDu2;
    }
    else {
        IsQuadu = false;
    }

    Dv0 = DotVector( &Nb, &Pd);
    Dv1 = DotVector(&Na, &Pd) + DotVector(&Nb, &Pc);
    Dv2 = DotVector( &Na, &Pc);
    if( fabs( Dv2 ) > SingleTinyVal ) {
        float TwoDv2;

        IsQuadv = true;
        TwoDv2 = 2.0 * Dv2;
        DivVector( &Na, TwoDv2, &Qvx);
        DivVector( &Nb, -Dv2, &Qvy);
/*      DivVector( &Nc, -Dv2, &Qvy); */
        Dvx = - Dv1/TwoDv2;
        Dvy = Dv0/Dv2;
    }
    else {
        IsQuadv = false;
    }
#ifdef DEBUG
    if( DebugOn ) {
        printf("du2 %f, du1 %f, du0 %f\n", Du2, Du1, Du0 );
        printf("dv2 %f, dv1 %f, dv0 %f\n", Dv2, Dv1, Dv0 );
        printf("Na = (%f, %f), Nb = (%f,%f), Nc = (%f,%f)\n",
        Na.x, Na.y, Nb.x, Nb.y, Nc.x, Nc.y );
        printf("IsQuad =(%c, %c)\n", IsQuadu?'t':'f', IsQuadv? 't': 'f' );
    }
#endif DEBUG

}


/*
 * given points px,py in the quadrilateral, map them to points inside a
 * unit square
 */
void  MapXYRatio(px, py, outx, outy, SweepCode)
    float px, py;
    float *outx, *outy;
    short SweepCode;
{
    float resu, resv;
    Vector Ri;


    Ri.x = px; Ri.y = py;

    if( !IsQuadu ) {
        float denom;
        denom = (Du1 - DotVector(&Na, &Ri));
        if( fabs(denom) < SingleTinyVal )
            resu = 2.0;
        else
            resu = (DotVector(&Nc, &Ri) - Du0)/denom;
    } else {
        float Ka, Kb;
        float discrim;

        Ka = Dux + DotVector( &Qux, &Ri);
        Kb = Duy + DotVector( &Quy, &Ri);
        discrim = sqrt(fabs(Ka * Ka - Kb));
        resu = Ka + ((discrim > Ka)? discrim: (-discrim));
#ifdef DEBUG
        if( DebugOn ) {
            printf("dux=%f, duy = %f, ka = %f, kb = %f\n",
                Dux, Duy, Ka, Kb );
        }
#endif DEBUG

    }

    if( !IsQuadv ) {
        float denom;
        denom = (Dv1 - DotVector(&Na, &Ri));
        if( fabs(denom) < SingleTinyVal )
            resv = 2.0;
        else
            resv = (DotVector(&Nb, &Ri) - Dv0)/denom;
    } else {
        float Ka, Kb;
        float discrim;

        Ka = Dvx + DotVector( &Qvx, &Ri);
        Kb = Dvy + DotVector( &Qvy, &Ri);
        discrim = sqrt(fabs( Ka * Ka - Kb));
        resv = Ka + ((discrim > Ka)? discrim: (-discrim));
#ifdef DEBUG
        if( DebugOn ) {
            printf("dvx=%f, dvy = %f, ka = %f, kb = %f\n",
                Dvx, Dvy, Ka, Kb );
        }
#endif DEBUG
    }

#ifdef DEBUG
    if( DebugOn ) {
        printf("(%f,%f) -> (%f, %f)\n", px, py, resu, resv );
    }
#endif DEBUG

    if( resu > 1.0 || resu < 0.0 ) {
        resu = ( SweepCode & MP_XMAX)? 1.0: 0.0;
    }
    if( resv > 1.0 || resv < 0.0 ) {
        resv = ( SweepCode & MP_YMAX)? 1.0: 0.0;
    }

    *outx = resu; *outy = resv;
}



/*
 * here abides testcode
 */
#ifdef TEST
#include <stdio.h>

ReadScrnPair(set, a)
    char *set;
    ScrnPair *a;
{
    int tx, ty;
    printf("enter screen pair %s\n",set);
    scanf("%d %d", &tx, &ty);
    a->x = tx; a->y = ty;
}

ReadLimits(a)
ScrnPair a[];
{
    ReadScrnPair("a", &a[0]);
    ReadScrnPair("b", &a[1]);
    ReadScrnPair("c", &a[2]);
    ReadScrnPair("d", &a[3]);
}

main() {
    float inx, iny;
    float outy, outx;
    ScrnPair pts[4];

    ReadLimits(pts);
    CalcMapConsts(pts);
    while(!feof(stdin)) {
        printf("enter quadrilateral points\n");
        scanf("%f %f", &inx, &iny );
        MapXYRatio( inx, iny, &outx, &outy, 0);
        printf("p(%f,%f)->p(%f,%f)\n", inx, iny, outx, outy);
    }
}
#endif TEST