/* MetaContour, version 0.1 -------------------------------------- */
/* Copyright(C) 2004, Brooks Moses                                 */
/*                                                                 */
/* This version of MetaContour is made available under the Gnu     */
/* Public License; see metacontour_main.cc for details.            */
/*                                                                 */
/* This is a very pre-release version of MetaContour, distributed  */
/* primarily as an example of a use of MetaPlot.  It can be        */
/* compiled with gcc using the line:                               */
/*                                                                 */
/*   g++ metacontour.cc cpoint.cc metacontour_main.cc              */
/*                                                                 */
/*-----------------------------------------------------------------*/

/* Class for 2D points with one data value each */  
  
#include<sstream>
#include"cpoint.h"
using namespace std;


//-----cpoint--------------------------------------------------

cpoint::cpoint() {
  x = 0;
  y = 0;
  z = 0;
}

cpoint::cpoint(double xx, double yy, double zz) {
  x = xx;
  y = yy;
  z = zz;
}

cpoint cpoint::interpolateto(cpoint x2, double z3) {
  double w = (z3-z)/(x2.z-z);
  if (w > 1) w=1;
  if (w < 0) w=0;
  double x3 = x + (x2.x-x)*w;
  double y3 = y + (x2.y-y)*w;
  return cpoint(x3,y3,z3);
}

cpoint cpoint::partwayto(cpoint x2, double w) {
  if (w > 1) w=1;
  if (w < 0) w=0;
  double x3 = x + (x2.x-x)*w;
  double y3 = y + (x2.y-y)*w;
  double z3 = z + (x2.z-z)*w;
  return cpoint(x3,y3,z3);
}

// Derivation:
//   F = ax + by + cxy + d     x,y in [0,1]
//   Fx = a + by, Fy = b + cx
//   Thus: saddle point S is at (-b/c, -a/c)
//   Plugging this in, and solving, we get
//   F(S) = (f00*f11 - f01*f10) / (f00 + f11 - f01 - f10)
//
//   Note: mapping x1 = (0,0), x2 = (1,0), x3 = (1,1), x4 = (0,1)
//
double cpoint::saddlepoint(cpoint x1, cpoint x2, cpoint x3, cpoint x4) {
  return (x1.z * x3.z - x2.z * x4.z) / (x1.z + x3.z - x2.z - x4.z);
}

string cpoint::metapoint() {
  ostringstream pointstring;
  pointstring.setf(ios_base::fixed,ios_base::floatfield);
  pointstring.precision(5);
  pointstring << '(' << x << ',' << y << ')';  
  return pointstring.str();
}


//-----cparray-------------------------------------------------

cparray::cparray(int x, int y) {
  if (x<1) x=1;
  if (y<1) y=1;
  i = x;
  j = y;
  a = new valarray<cpoint>(x*y);
}

cparray::~cparray() {
  delete a;
}

cpoint& cparray::operator()(int x, int y)
{
  return (*a)[x*j+y];
}