/* "UNTIL", a graphics editor,
   Copyright (C) 1985, 1990 California Institute of Technology.
   Original authors: Glenn Gribble, port by Steve DeWeerth
   Unix Port Maintainer: John Lazzaro
   Maintainers's address: lazzaro@hobiecat.cs.caltech.edu;
                          CB 425 CU Boulder/Boulder CO 91125. 

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (Version 1, 1989).

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; see the file COPYING.  If not, write to the
Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139,
USA. */


/*******************************************************************************/
/*                                                                             */ 
/*  procedures to mangle matrices                                              */
/*  cleaned up by steve - 9 May 1990                                           */
/*                                                                             */ 
/*******************************************************************************/


#ifdef mips
#define const
#endif

#include <p2c/p2c.h>

#include "mat_stuff.h"

#ifndef SYSGLOBALS_H
#include <p2c/sysglobals.h>
#endif

#ifndef MYLIB_H
#include <p2c/mylib.h>
#endif


/* exported global variables */

const mat_matrix mat_ident_mat = {
  1, 0, 0, 1, 0, 0, 1
};

mat_matrix ctm;
double *rSinI;



/* to get m_across/m_down */
void mat_init(void)
{
  long i;

  rSinI = Malloc(sizeof(mat_r_array));
  for (i = 0; i <= 90; i++)
    rSinI[i] = sin(i * pi / 180.0);
}


/* Find the amount of rotation in the current transformation matrix */
long mat_findRot(void)
{
  long x, y;

  mat_rtransform(1024, 0, &x, &y);
  return (aTan2I(x, y));
}


/* Compute the new length of something (works for same X and Y scaling) */
long mat_newLength(long oldLen)
{
  long rx, ry;

  mat_rtransform(oldLen, 0, &rx, &ry);
  if (ry == 0)
    return labs(rx);
  else if (rx == 0)
    return labs(ry);
  else
    return ((long)floor(sqrt((double)(rx * rx + ry * ry)) + 0.5));
}


/* Return the Sine of D (in degrees) */
double sinD(long d)
{
  long i;
  double m;

  d %= 360;
  i = d % 90;
  d /= 90;

  switch (d) {
  case 0:
    m = rSinI[i];
    break;
  case 1:
    m = rSinI[90 - i];
    break;
  case 2:
    m = -rSinI[i];
    break;
  case 3:
    m = -rSinI[90 - i];
    break;
  }

  return m;
}


double cosD(long d)
{
  return (sinD(d + 90));
}


#define pi_             3.14159265358979

long aTan2I(double s, double c)
{
  long t;

  if (c == 0)
    t = 90;
  else
    t = (long)floor(180 * atan(s / c) / pi_ + 0.5);

  /* Now t is in range -90<t<=90, check for quads 2 and 3 and adjust */
  if (c < 0) {
    if (s < 0)   /* quad 3 */
      t -= 180;
    else  /* quad 2 */
      t += 180;
  }

  t = (90 - t) % 360;   /* The minus is because things are screwed up. */

  return t;
}

#undef pi_


/* Save the matrix state */
void mat_push(mat_state *m)
{
  m->m = ctm;
}


/* Restore the matrix state */
void mat_pop(mat_state *m)
{
  ctm = m->m;
}


const mat_matrix ident_mat = {
  1, 0, 0, 1, 0, 0, 1
};


/* This is not really an identity matrix because we like to have 0,0 at the center of the screen */
void mat_ident(void)
{
  ctm = ident_mat;
}


/* What is the current matrix??? */
void mat_print(void)
{
  printf("|%6ld%6ld|    1\n", ctm.a, ctm.b);
  printf("|%6ld%6ld|   ---\n", ctm.c, ctm.d);
  printf("|%6ld%6ld|%6ld\n", ctm.e, ctm.f, ctm.g);
}


/* call me whenever there seems to be an overflow error */
Static void smash_matrix(void)
{
  if (P_escapecode != -4)
    _Escape(P_escapecode);

  printf("\201o\200");

  ctm.a = (ctm.a + 1) / 2;
  ctm.b = (ctm.b + 1) / 2;
  ctm.c = (ctm.c + 1) / 2;
  ctm.d = (ctm.d + 1) / 2;
  ctm.e = (ctm.e + 1) / 2;
  ctm.f = (ctm.f + 1) / 2;
  ctm.g = (ctm.g + 1) / 2;
}


#define numprimes       20
typedef short prime_list[numprimes];
const prime_list primes = {
  2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 71, 79
};

Static void reduce_matrix(void)
{
  long i;
  short p;
  
  for (i = 0; i < numprimes; i++) {
    p = primes[i];
    /* assembly type thing to do is */
    /* if element is zero, ignore it */
    /* if element is < current_prime, we are done reducing */
    /* if element is not divisible by p, skip this prime */
    /* if element was divisible, be sure to save intermediate result */
    while (ctm.g % p == 0 && ctm.f % p == 0 && ctm.e % p == 0 && ctm.d % p == 0 &&
	   ctm.c % p == 0 && ctm.b % p == 0 && ctm.a % p == 0) {
      ctm.a /= p;
      ctm.b /= p;
      ctm.c /= p;
      ctm.d /= p;
      ctm.e /= p;
      ctm.f /= p;
      ctm.g /= p;
    }
  }
}

#undef numprimes


void mat_tran(long dx, long dy)
{
  ctm.e += ctm.a * dx + ctm.c * dy;
  ctm.f += ctm.b * dx + ctm.d * dy;
  reduce_matrix();
}


void mat_scale(long x, long y, long num)
{
  mat_matrix omat;
  long count;

  /* First, reduce x and y vs NUM */
  count = 10;

  do {
    omat = ctm;
    TRY(try1);
      ctm.a *= x;
      ctm.b *= x;
      ctm.c *= y;
      ctm.d *= y;
      ctm.e *= num;
      ctm.f *= num;
      ctm.g *= num;
      count = 0;
    RECOVER(try1);
      if (P_escapecode != -4)
	_Escape(P_escapecode);
      ctm = omat;
      printf("\201s\200");
      x = (x + 1) / 2;
      y = (y + 1) / 2;
      num = (num + 1) / 2;
      count--;
    ENDTRY(try1);
  } while (count > 0);
  reduce_matrix();
}


/* pre-multiply the ctm */
void mat_premult(mat_matrix *m)
{
  mat_matrix tmpmat;

  tmpmat.a = m->a * ctm.a + m->b * ctm.c;
  tmpmat.b = m->a * ctm.b + m->b * ctm.d;
  tmpmat.c = m->c * ctm.a + m->d * ctm.c;
  tmpmat.d = m->c * ctm.b + m->d * ctm.d;
  tmpmat.e = m->e * ctm.a + m->f * ctm.c + m->g * ctm.e;
  tmpmat.f = m->e * ctm.b + m->f * ctm.d + m->g * ctm.f;
  tmpmat.g = m->g * ctm.g;

  ctm = tmpmat;
  reduce_matrix();
}


void mat_rot(long c, long s)
{
  /* cosine, sine of rotation */
  long r;
  mat_matrix new_ctm;

  /* this should be done some faster way in the future */
  r = (long)floor(sqrt((double)(s * s + c * c)) + 0.5);

  new_ctm.a = c * ctm.a + s * ctm.c;
  new_ctm.b = c * ctm.b + s * ctm.d;
  new_ctm.c = c * ctm.c - s * ctm.a;
  new_ctm.d = c * ctm.d - s * ctm.b;
  new_ctm.e = r * ctm.e;
  new_ctm.f = r * ctm.f;
  new_ctm.g = r * ctm.g;
  ctm = new_ctm;
  reduce_matrix();
}


void mat_transform(long x, long y, long *rx, long *ry)
{
  TRY(try2);
    *rx = (x * ctm.a + y * ctm.c + ctm.e) / ctm.g;
    *ry = (x * ctm.b + y * ctm.d + ctm.f) / ctm.g;
  RECOVER(try2);
    if (P_escapecode != -4)
      _Escape(P_escapecode);
    _Escape(-4);
  ENDTRY(try2);
}


/* same as mat_transform, except leave out the e and f */
void mat_rtransform(long x, long y, long *rx, long *ry)
{
  TRY(try3);
    *rx = (x * ctm.a + y * ctm.c) / ctm.g;
    *ry = (x * ctm.b + y * ctm.d) / ctm.g;
  RECOVER(try3);
    if (P_escapecode != -4)
      _Escape(P_escapecode);
    _Escape(-4);
  ENDTRY(try3);
}


Static void mat_untransformR(long x, long y, long *rx, long *ry)
{
  double xt, yt, det, ra, rb, rc, rd, re, rf, rg, rx0, ry0;

  ra = ctm.a;
  rb = ctm.b;
  rc = ctm.c;
  rd = ctm.d;
  re = ctm.e;
  rf = ctm.f;
  rg = ctm.g;

  TRY(try4);
    det = rb * rc - ra * rd;
    if (det == 0) {
      *rx = 0;
      *ry = 0;
      printf("\201mat_untransformR: transformation matrix is singular\200\n");
    } 
    else {
      det = 1.0 / det;
      xt = rg * x - re;
      yt = rg * y - rf;
      rx0 = (rc * yt - rd * xt) * det;
      ry0 = (rb * xt - ra * yt) * det;

      if (rx0 >= LONG_MAX)
	*rx = LONG_MAX;
      else if (rx0 <= -LONG_MAX)
	*rx = -LONG_MAX;
      else
	*rx = (long)floor(rx0 + 0.5);

      if (ry0 >= LONG_MAX)
	*ry = LONG_MAX;
      else if (ry0 <= -LONG_MAX)
	*ry = -LONG_MAX;
      else
	*ry = (long)floor(ry0 + 0.5);
    }
  RECOVER(try4);
    if (P_escapecode != -4)
      _Escape(P_escapecode);
    printf("\201mat_untransform: overflow (with real #s)!\200\007\n");
    _Escape(-4);
  ENDTRY(try4);
}


void mat_untransform(long x, long y, long *rx, long *ry)
{
  long det, xt, yt;

  TRY(try5);
    det = ctm.b * ctm.c - ctm.a * ctm.d;
    if (det == 0) {
      *rx = 0;
      *ry = 0;
      printf("\201mat_untransform: transformation matrix is singular\200\n");
    } 
    else {
      xt = ctm.g * x - ctm.e;
      yt = ctm.g * y - ctm.f;
      *rx = (ctm.c * yt - ctm.d * xt) / det;
      *ry = (ctm.b * xt - ctm.a * yt) / det;
    }
  RECOVER(try5);
    if (P_escapecode != -4)
      _Escape(P_escapecode);
    printf("\201untransformR called.\200");
    mat_untransformR(x, y, rx, ry);
  ENDTRY(try5);
}


void mat_runtransform(long x, long y, long *rx, long *ry)
{
  long det, xt, yt;

  TRY(try6);
    det = ctm.b * ctm.c - ctm.a * ctm.d;
    if (det == 0) {
      *rx = 0;
      *ry = 0;
      printf("\201mat_untransform: transformation matrix is singular\200\n");
    }
    else {
      xt = ctm.g * x;
      yt = ctm.g * y;
      *rx = (ctm.c * yt - ctm.d * xt) / det;
      *ry = (ctm.b * xt - ctm.a * yt) / det;
    }
  RECOVER(try6);
    if (P_escapecode != -4)
      _Escape(P_escapecode);
    _Escape(-4);
  ENDTRY(try6);
}