/****************************************************************************/
/* (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 aagraph.h */

#include "aagraph.h"
#include "fltgraph.h"
#include <js.h>
#include <flt.h>
#include <ia.h>
#include <aa.h>
#include <ps.h>
#include <math.h>
#include <ioprotos.h>
#include <stdio.h>
#include <malloc.h>

/*** INTERNAL PROTOTYPES ***/

void aagraph_plot_aa_graph(
    FILE *psfile,
    AAP ff (AAP x),
    Interval xd,
    Interval yd,
    int n
  );

void aagraph_plot_aa_range_graph(
    FILE *psfile,
    AAP ff (AAP x),
    Interval xd,
    Interval yd,
    int n
  );

void aagraph_plot_ia_graph(
    FILE *psfile,
    Interval fv (Interval x),
    Interval xd,
    Interval yd,
    int n
  );
  
char *aagraph_format_parms(
    Interval xd,
    Interval yd,
    int n
  );
  
/*** IMPLEMENTATIONS ***/

void aagraph_plot_aa_graph(
    FILE *psfile,
    AAP ff (AAP x),
    Interval xd,
    Interval yd,
    int n
  )
  {
    Interval xv, yvlo, yvhi;
    int xi;
    AAP xf, yf;
    double xp[4], yp[4];
    AATerm eps[1];
    double gray = 0.75;

    ps_begin_section(psfile, "Plot of AA graph");

    for (xi=0; xi<n; xi++)
      {
        MemP frame = aa_top();

        ROUND_DOWN;
        xv.lo = xd.lo + ((xd.hi - xd.lo)*xi)/n;

        ROUND_UP;
        xv.hi = xd.lo + ((xd.hi - xd.lo)*(xi+1))/n;

        xf = aa_from_interval(xv);
        assert((xf->nterms == 1), "aagraph_plot_aa_graph: xf->nterms != 1");
        eps[0].id = ((AATermP)(xf + 1))->id;

        yf = ff(xf);
        
        eps[0].coef = -One;
        yvlo = aa_range(aa_fix_eps(yf, 1, eps));
        eps[0].coef = One;
        yvhi = aa_range(aa_fix_eps(yf, 1, eps));

        if (ia_is_full(&yvlo) || ia_is_full(&yvhi)) 
          { yvlo = yd; yvhi = yd; }

        ROUND_NEAR;
        xp[0] = xv.lo;  yp[0] = yvlo.lo;
        xp[1] = xv.hi;  yp[1] = yvhi.lo;
        xp[2] = xv.hi;  yp[2] = yvhi.hi;
        xp[3] = xv.lo;  yp[3] = yvlo.hi;
        
        ps_fill_and_draw_polygon(psfile, xp, yp, 4, gray);

        aa_flush(frame);
      }

    fprintf(stderr, "\n");
    ps_end_section(psfile);
  }

void aagraph_plot_aa_range_graph(
    FILE *psfile,
    AAP ff (AAP x),
    Interval xd,
    Interval yd,
    int n
  )
  {
    Interval xv, yv;
    int xi;
    AAP xf, yf;
    double gray = 0.75;

    ps_begin_section(psfile, "Plot of AA range graph");

    for (xi=0; xi<n; xi++)
      {
        MemP frame = aa_top();

        ROUND_DOWN;
        xv.lo = xd.lo + ((xd.hi - xd.lo)*xi)/n;

        ROUND_UP;
        xv.hi = xd.lo + ((xd.hi - xd.lo)*(xi+1))/n;

        xf = aa_from_interval(xv);

        yf = ff(xf);
        yv = aa_range(yf);

        if (ia_is_full(&yv)) yv = yd;

        ROUND_NEAR;
        ps_fill_and_draw_rectangle(psfile, xv.lo, xv.hi, yv.lo, yv.hi, gray);

        aa_flush(frame);
      }

    fprintf(stderr, "\n");
    ps_end_section(psfile);
  }

void aagraph_plot_ia_graph(
    FILE *psfile,
    Interval fv (Interval x),
    Interval xd,
    Interval yd,
    int n
  )
  {
    Interval xv, yv;
    int xi;
    double gray = 0.75;

    ps_begin_section(psfile, "Plot of IA graph");

    for (xi=0; xi<n; xi++)
      {
        ROUND_DOWN;
        xv.lo = xd.lo + ((xd.hi - xd.lo)*xi)/n;

        ROUND_UP;
        xv.hi = xd.lo + ((xd.hi - xd.lo)*(xi+1))/n;

        yv = fv(xv);

        if (ia_is_full(&yv)) yv = yd;

        ROUND_NEAR;
        ps_fill_and_draw_rectangle(psfile, xv.lo, xv.hi, yv.lo, yv.hi, gray);

      }

    fprintf(stderr, "\n");
    ps_end_section(psfile);
  }
  
char *aagraph_format_parms(
    Interval xd,
    Interval yd,
    int n
  )
  {
    char *s = (char *) malloc(100);
    sprintf(s, 
      "x in [%f _ %f]  y in [%f _ %f]  %d intervals",
      xd.lo, xd.hi, yd.lo, yd.hi, n
    );
    return(s);
  }

void aagraph_plots(
    char *fileprefix,
    int epsformat,
    char *title,
    Float f(Float x),
    Interval fv(Interval x),
    AAP ff(AAP x),
    Interval xd,
    Interval yd,
    int n,
    int m
  )
  {
    char *parmstr = aagraph_format_parms(xd, yd, n);
    FILE *psfile;
    Float wx = xd.hi - xd.lo;
    Float wy = yd.hi - yd.lo;
    Float wm = FMAX(wx, wy);
    double scale = 6.0 * 72.0 / wm;
    double xc = 4.25 * 72.0;
    double hmin = xc - scale * wx / 2.0;
    double hmax = xc + scale * wx / 2.0;
    double yc = 6.50 *72.0;
    double vmin = yc - scale * wy / 2.0;
    double vmax = yc + scale * wy / 2.0;

    /*** Ordinary interval arithmetic plot ***/
    
    if (epsformat)
      {
        psfile = open_write(txtcat(fileprefix, "-ia.eps"));
        ps_begin_figure(
	  psfile,
	  xd.lo, xd.hi, yd.lo, yd.hi,
	  hmin,  hmax,  vmin,  vmax,
	  n, 1
	);
      }
    else
      {
        psfile = open_write(txtcat(fileprefix, ".ps"));
        ps_begin_document(psfile);
        ps_begin_page(
	  psfile, 1,
	  xd.lo, xd.hi, yd.lo, yd.hi,
	  hmin,  hmax,  vmin,  vmax,
	  n, 1
	);
        ps_add_caption(psfile, title);
        ps_add_caption(psfile, "");
        ps_add_caption(psfile, parmstr);
        ps_add_caption(psfile, "");
        ps_add_caption(psfile, "Ordinary interval arithmetic");
      }

    aagraph_plot_ia_graph(psfile, fv, xd, yd, n);
    fltgraph_draw_axes(psfile, xd, yd);
    fltgraph_plot(psfile, f, xd, yd, m);
    ps_draw_frame(psfile);

    if (epsformat)
      { ps_end_figure(psfile); fclose(psfile); }
    else
      { ps_end_page(psfile); }
      
    /*** Affine arithmetic plot ***/

    if (epsformat)
      {
        psfile = open_write(txtcat(fileprefix, "-aa.eps"));
	ps_begin_figure(
	  psfile,
	  xd.lo, xd.hi, yd.lo, yd.hi,
	  hmin,  hmax,  vmin,  vmax,
	  n, 1
	);
      }
    else
      {
        ps_begin_page(
	  psfile, 2,
	  xd.lo, xd.hi, yd.lo, yd.hi,
	  hmin,  hmax,  vmin,  vmax,
	  n, 1
	);
	ps_add_caption(psfile, title);
	ps_add_caption(psfile, "");
	ps_add_caption(psfile, parmstr);
	ps_add_caption(psfile, "");
	ps_add_caption(psfile, "Affine arithmetic");
      }

    aagraph_plot_aa_graph(psfile, ff, xd, yd, n);
    fltgraph_draw_axes(psfile, xd, yd);
    fltgraph_plot(psfile, f, xd, yd, m);
    ps_draw_frame(psfile);
    
    if (epsformat) 
      { ps_end_figure(psfile); fclose(psfile); }
    else
      { ps_end_page(psfile); }

    /*** Affine arithmetic range plot ***/

    if (epsformat)
      {
        psfile = open_write(txtcat(fileprefix, "-ar.eps"));
	ps_begin_figure(
	  psfile,
	  xd.lo, xd.hi, yd.lo, yd.hi,
	  hmin,  hmax,  vmin,  vmax,
	  n, 1
	);
      }
    else
      {
        ps_begin_page(
	  psfile, 3,
	  xd.lo, xd.hi, yd.lo, yd.hi,
	  hmin,  hmax,  vmin,  vmax,
	  n, 1
	);
	ps_add_caption(psfile, title);
	ps_add_caption(psfile, "");
	ps_add_caption(psfile, parmstr);
	ps_add_caption(psfile, "");
	ps_add_caption(psfile, "Affine arithmetic (converted to intervals)");
      }

    aagraph_plot_aa_range_graph(psfile, ff, xd, yd, n);
    fltgraph_draw_axes(psfile, xd, yd);
    fltgraph_plot(psfile, f, xd, yd, m);
    ps_draw_frame(psfile);
    
    if (epsformat) 
      { ps_end_figure(psfile); fclose(psfile); }
    else
      { ps_end_page(psfile); ps_end_document(psfile, 3); fclose(psfile); }
  }