/* See r2.h */
/* Last edited on 2005-01-16 17:15:14 by stolfi */
/* Copyright © 2005 Jorge Stolfi, UNICAMP. See note at end of file. */

#define _GNU_SOURCE
#include <r2.h>
#include <math.h>
#include <stdio.h>
#include <js.h>
#include <rn.h>

#define N 2

void r2_zero (r2_t *r)
  { r->c[0] = 0.0;
    r->c[1] = 0.0;
  }

void r2_all (double x, r2_t *r)
  { r->c[0] = x;
    r->c[1] = x;
  }

void r2_axis (int i, r2_t *r)
  { affirm((i >= 0) && (i < N), "r2_axis: bad index");
    r->c[0] = 0.0;
    r->c[1] = 0.0;

    r->c[i] = 1.0;
  }

void r2_add (r2_t *a, r2_t *b, r2_t *r)
  { r->c[0] = a->c[0] + b->c[0];
    r->c[1] = a->c[1] + b->c[1];
  }

void r2_sub (r2_t *a, r2_t *b, r2_t *r)
  { r->c[0] = a->c[0] - b->c[0];
    r->c[1] = a->c[1] - b->c[1];
  }

void r2_neg (r2_t *a, r2_t *r)
  { r->c[0] = - a->c[0];
    r->c[1] = - a->c[1];
  }

void r2_scale (double s, r2_t *a, r2_t *r)
  { r->c[0] = s * a->c[0];
    r->c[1] = s * a->c[1];
  }

void r2_mix (double s, r2_t *a, double t, r2_t *b, r2_t *r)
  { r->c[0] = s * a->c[0] + t * b->c[0];
    r->c[1] = s * a->c[1] + t * b->c[1];
  }

void r2_mix_in (double s, r2_t *a, r2_t *r)
  { r->c[0] += s * a->c[0];
    r->c[1] += s * a->c[1];
  }

void r2_weigh (r2_t *a, r2_t *w, r2_t *r)
  { r->c[0] = a->c[0] * w->c[0];
    r->c[1] = a->c[1] * w->c[1];
  }

double r2_norm (r2_t *a)
  { return hypot(a->c[0], a->c[1]); }

double r2_norm_sqr (r2_t *a)
  { double a0 = a->c[0];
    double a1 = a->c[1];
    return a0*a0 + a1*a1;
  }

double r2_L_inf_norm (r2_t *a)
  { double d = 0.0;
    double a0 = fabs(a->c[0]);
    double a1 = fabs(a->c[1]);
    if (a0 > d) d = a0;
    if (a1 > d) d = a1;
    return (d);
  }

double r2_dist (r2_t *a, r2_t *b)
  { double d0 = (a->c[0] - b->c[0]);
    double d1 = (a->c[1] - b->c[1]);
    return hypot(d0, d1);
  }

double r2_dist_sqr (r2_t *a, r2_t *b)
  { double d0 = (a->c[0] - b->c[0]);
    double d1 = (a->c[1] - b->c[1]);
    return d0*d0 + d1*d1;
  }

double r2_L_inf_dist (r2_t *a, r2_t *b)
  { double d = 0.0;
    double d0 = fabs(a->c[0] - b->c[0]);
    double d1 = fabs(a->c[1] - b->c[1]);
    if (d0 > d) d = d0;
    if (d1 > d) d = d1;
    return (d);
  }

double r2_dir (r2_t *a, r2_t *r)
  { double d = sqrt(a->c[0]*a->c[0] + a->c[1]*a->c[1]);
    r->c[0] = a->c[0]/d;
    r->c[1] = a->c[1]/d;
    return (d);
  }

double r2_L_inf_dir (r2_t *a, r2_t *r)
  { double d = 0.0;
    double a0 = fabs(a->c[0]);
    double a1 = fabs(a->c[1]);
    if (a0 > d) d = a0;
    if (a1 > d) d = a1;
    r->c[0] = a->c[0]/d;
    r->c[1] = a->c[1]/d;
    return (d);
  }

double r2_dot (r2_t *a, r2_t *b)
  { return a->c[0]*b->c[0] + a->c[1]*b->c[1]; }

double r2_cos (r2_t *a, r2_t *b)
  { double a0 = a->c[0];
    double a1 = a->c[1];
    double b0 = b->c[0];
    double b1 = b->c[1];
    double ab = a0*b0 + a1*b1;
    double aa = a0*a0 + a1*a1;
    double bb = b0*b0 + b1*b1;
    return ab/(sqrt(aa)*sqrt(bb));
  }

double r2_sin (r2_t *a, r2_t *b)
  { return rn_sin(N, &(a->c[0]), &(b->c[0])); }

void r2_cross (r2_t *a, r2_t *r)
  { double a0 = a->c[0];
    double a1 = a->c[1];
    r->c[0] = - a1;
    r->c[1] =   a0;
  }

double r2_det (r2_t *a, r2_t *b)
  { return a->c[0]*b->c[1] - a->c[1]*b->c[0]; }

double r2_decomp (r2_t *a, r2_t *u, r2_t *para, r2_t *perp)
  { double u0 = u->c[0];
    double u1 = u->c[1];
    double sau = a->c[0]*u0 + a->c[1]*u1;
    if (sau == 0.0) 
      { if (para != NULL)
          { para->c[0] = 0.0; 
            para->c[1] = 0.0;
          }
        if (perp != NULL)
          { perp->c[0] = a->c[0]; 
            perp->c[1] = a->c[1];
          }
        return 0.0;
      }
    else
      { double suu = u0*u0 + u1*u1;
        double c = sau / suu;
        double p0 = c*u0; 
        double p1 = c*u1; 

        if (para != NULL)
          { para->c[0] = p0;
            para->c[1] = p1;
          }

        if (perp != NULL)
          { perp->c[0] = a->c[0] - p0; 
            perp->c[1] = a->c[1] - p1;
          }
        return c;
      }
  }

void r2_throw_cube (r2_t *r)
  { r->c[0] = 2.0 * drandom() - 1.0;
    r->c[1] = 2.0 * drandom() - 1.0;
  }

void r2_throw_ball (r2_t *r)
  { double x, y;
    do
      { x = 2.0 * drandom() - 1.0; r->c[0] = x;
        y = 2.0 * drandom() - 1.0; r->c[1] = y;
      }
    while (x*x + y*y >= 1.0);
  }

void r2_throw_normal (r2_t *r)
  { r->c[0] = dgaussrand();
    r->c[1] = dgaussrand();
  }

void r2_print (FILE *f, r2_t *a)
  { r2_gen_print(f, a, NULL, NULL, NULL, NULL); }

void r2_gen_print (FILE *f, r2_t *a, char *fmt, char *lp, char *sep, char *rp)
  { rn_gen_print(f, N, &(a->c[0]), fmt, lp, sep, rp); }

/* COPYRIGHT AND AUTHORSHIP NOTICE

  Copyright © 2005 Jorge Stolfi, Universidade Estadual de Campinas (UNICAMP).
  Created by Jorge Stolfi in 1992--2005.
  
  This source file can be freely distributed, used, and modified,
  provided that this copyright and authorship notice is preserved in
  all copies, and that any modified versions of this file are clearly
  marked as such.
  
  This software has NO WARRANTY of correctness or applicability for
  any purpose. Neither the author nor his employers shall be held
  responsible for any losses or damages that may result from its use.

  END OF NOTICE */