#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "galex_footprint.h"


void radec2xyz(double ra, double dec, double *x, double *y, double *z)
{
    double a, d, ca, cd, sa, sd;
    a = ra*d2r;
    d = dec*d2r;
    ca = cos(a); sa = sin(a);
    cd = cos(d); sd = sin(d);
    *x = cd * ca;
    *y = cd * sa;
    *z = sd;
}

void normalize(double *x, double *y, double *z)
{
    double len = sqrt( (*x)*(*x) + (*y)*(*y) + (*z)*(*z) );
    *x /= len;
    *y /= len;
    *z /= len;
}

point* point_new(double ra, double dec)
{
    point* p = (point*) malloc(sizeof(point));
    radec2xyz(ra,dec, &p->x,&p->y,&p->z);
    return p;
}

point* point_set(point *p, double ra, double dec)
{
  radec2xyz(ra,dec, &p->x,&p->y,&p->z);
    return p;
}

void fprint_point(FILE* stream, point* p)
{
    fprintf(stream, "%lf %lf %lf", p->x, p->y, p->z);
}


void halfspace_set(halfspace* h, double x,double y,double z,bool norm,double c)
{
    h->x = x;
    h->y = y;
    h->z = z;
    h->c = c;
    if (norm) 
        normalize(&h->x,&h->y,&h->z);
}

halfspace* halfspace_new(double x, double y, double z, bool norm, double c)
{
    halfspace* h = (halfspace*) malloc(sizeof(halfspace));
    halfspace_set(h, x,y,z, norm, c);
    return h;
}

void fprint_halfspace(FILE *stream, halfspace *h)
{
    fprintf(stream, "%lf %lf %lf %lf", h->x, h->y, h->z, h->c);
}

bool halfspace_covers(halfspace* h, point* p, double tol)
{
    double dot = h->x*p->x + h->y*p->y + h->z*p->z;
    return dot > h->c - tol;
}

void convex_set(convex *c, int num)
{
    c->halfspacelist = (halfspace*) malloc(num*sizeof(halfspace));
    c->nh = num;
}    

convex* convex_new(int num)
{
    convex *c = (convex*) malloc(sizeof(convex));
    convex_set(c, num);
    return c;
}

bool convex_covers(convex* c, point* p, double tol)
{
    int i;
    for (i=0; i<c->nh; i++)
        if (! halfspace_covers(c->halfspacelist+i,p,tol)) return FALSE;
    return TRUE;
}

void convex_free(convex *c)
{
    if (c->halfspacelist != NULL) free(c->halfspacelist);
    if (c != NULL) { free(c); c = NULL; }
}

region* region_new(int num)
{
    region *r = (region*) malloc(sizeof(region));
    r->convexlist = (convex*) malloc(num*sizeof(convex));
    r->nc = num;
    return r;
}

void region_free(region *r)
{
    int i;
    if (r->convexlist != NULL)
    {
        for (i=0; i<r->nc; i++) 
	    if (r->convexlist[i].halfspacelist != NULL)
	        free(r->convexlist[i].halfspacelist);
	free(r->convexlist);
	r->convexlist = NULL;
    }
    if (r != NULL) { free(r); r = NULL; }
}

int region_firstcovers(region *r, point *p, double tol)
{
    int i;
    for (i=0; i<r->nc; i++)
        if (convex_covers(r->convexlist+i, p, tol))
	    return i;
    return -1;
}