/****************************************************************************/
/* (C) Copyright 1993 Universidade Estadual de Campinas (UNICAMP)           */
/*                    Campinas, SP, Brazil                                  */
/*                                                                          */
/* This file can be freely distributed, modified, and used for any          */
/*   non-commercial purpose, provided that this copyright and authorship    */
/*   notice be included in any copy or derived version of this file.        */
/*                                                                          */
/* DISCLAIMER: This software is offered ``as is'', without any guarantee    */
/*   as to fitness for any particular purpose.  Neither the copyright       */
/*   holder nor the authors or their employers can be held responsible for  */
/*   any damages that may result from its use.                              */
/****************************************************************************/

/* See fltzeros2.h */

#include "fltzeros2.h"
#include <flt.h>
#include <js.h>
#include <ia.h>
#include <ps.h>
#include <math.h>
#include <ioprotos.h>
#include <stdio.h>

/*** PROTOTYPES FOR INTERNAL ROUTINES ***/

void fltzeros2_in_triangle(
    FILE *psfile,
    double xa, double ya, Float za,
    double xb, double yb, Float zb,
    double xc, double yc, Float zc
  );
  /* Draws intersection of triangle abc with plane z=0. */
  /* Includes edges ab and bc but not ac.               */

/*** IMPLEMENTATIONS ***/

void fltzeros2_plot(
    FILE *psfile,
    Float f (Float x, Float y),
    Interval xd,
    Interval yd,
    int m
  )
  {
    Interval xv, yv;
    int xi, yi;
    Float xm, ym;
    Float f00, f01, f10, f11, fmm;

    ROUND_NEAR;
    ps_begin_section(psfile, "Plot of actual zeros");

    ps_set_pen(psfile, 0.0, 0.30, 0.0, 0.0);

    for (yi=0; yi<m; yi++)
      {
        for (xi=0; xi<m; xi++)
          {
            ROUND_DOWN;
            xv.lo = xd.lo + ((xd.hi - xd.lo)*xi)/m;
            yv.lo = yd.lo + ((yd.hi - yd.lo)*yi)/m;

            ROUND_UP;
            xv.hi = xd.lo + ((xd.hi - xd.lo)*(xi+1))/m;
            yv.hi = yd.lo + ((yd.hi - yd.lo)*(yi+1))/m;

            ROUND_NEAR;
            xm = (xv.lo + xv.hi)/Two;
            ym = (yv.lo + yv.hi)/Two;

            f00 = f(xv.lo, yv.lo);
            f01 = f(xv.lo, yv.hi);
            f10 = f(xv.hi, yv.lo);
            f11 = f(xv.hi, yv.hi);
            fmm = f(xm, ym);

            fltzeros2_in_triangle(psfile,
              xv.lo,  yv.lo,  f00,
              xv.hi,  yv.lo,  f10,
              xm,     ym,     fmm
            );

            fltzeros2_in_triangle(psfile,
              xv.hi,  yv.lo,  f10,
              xv.hi,  yv.hi,  f11,
              xm,     ym,     fmm
            );

            fltzeros2_in_triangle(psfile,
              xv.hi,  yv.hi,  f11,
              xv.lo,  yv.hi,  f01,
              xm,     ym,     fmm
            );

            fltzeros2_in_triangle(psfile,
              xv.lo,  yv.hi,  f01,
              xv.lo,  yv.lo,  f00,
              xm,     ym,     fmm
            );

          }
      }

    ps_end_section(psfile);
  }

void fltzeros2_in_triangle(
    FILE *psfile,
    double xa, double ya, Float za,
    double xb, double yb, Float zb,
    double xc, double yc, Float zc
  )
  { if (za == Zero && zb == Zero && zc == Zero)
      {
        ps_fill_triangle(psfile, xa, ya, xb, yb, xc, yc, 0.0);
      }
    if (za == Zero && zb == Zero)
      {
        ps_draw_segment(psfile, xa, ya, xb, yb);
      }
    if (zb == Zero && zc == Zero)
      {
        ps_draw_segment(psfile, xb, yb, xc, yc);
      }
    if ((za == Zero) + (zb == Zero) + (zc == Zero) >= 2) return;

    if (za >= Zero && zb >= Zero && zc >= Zero) return;
    if (za <= Zero && zb <= Zero && zc <= Zero) return;

    if ((zc == Zero) || ((zb > Zero) != (zc > Zero)))
      {
        double t; Float s;
        t = xa; xa = xc; xc = t;
        t = ya; ya = yc; yc = t;
        s = za; za = zc; zc = s;
      }
    /* Now zc != 0. */
    /* Also za == 0, or zb == 0, or (sign(za) != sign(zb)) */
    if ((zb == Zero) || ((zb > Zero) != (zc > Zero)))
      {
        double t; Float s;
        t = xa; xa = xb; xb = t;
        t = ya; ya = yb; yb = t;
        s = za; za = zb; zb = s;
      }
    /* Now zb != 0, zc != 0. */
    /* Also za = 0, or (sign(za) != sign(zb) and sign(za) != sign(zc)) */
    if (za != Zero)
      { 
        double sab = zb/(zb - za);
	double sba = za/(za - zb);
	double xab = xa*sab + xb*sba;
	double yab = ya*sab + yb*sba;
	double sac = zc/(zc - za);
	double sca = za/(za - zc);
	double xac = xa*sac + xc*sca;
	double yac = ya*sac + yc*sca;
	ps_draw_segment(psfile, xab, yab, xac, yac);
      }
    else 
      {
	double sbc = zc/(zc - zb);
	double scb = zb/(zb - zc);
	double xbc = xb*sbc + xc*scb;
	double ybc = yb*sbc + yc*scb;
        assert((zb > 0) != (zc > 0), "fltzeros2_in_triangle: sign confusion");
	ps_draw_segment(psfile, xa, ya, xbc, ybc);
      }
  }