/*
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 *
 *<LicenseText>
 *
 * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
 * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
 * Copyright (C) 1994-2005, California Institute of Technology.
 *
 * 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; either version 2 of the License, or
 * (at your option) any later version.
 *
 * 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; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 *</LicenseText>
 *
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 */


#include <vector>
#include <list>
#include <math.h>

/* forward declaration */
struct All_variables;

#ifndef _TRACER_DEFS_H_
#define _TRACER_DEFS_H_

typedef int ElementID;
#define UNDEFINED_ELEMENT		((ElementID)-99)

class CartesianCoord;

// Position or vector in spherical coordinates
class SphericalCoord {
public:
	double	_theta, _phi, _rad;
	SphericalCoord(void) : _theta(0), _phi(0), _rad(0) {};
	SphericalCoord(double theta, double phi, double rad) : _theta(theta), _phi(phi), _rad(rad) {};
	
	size_t size(void) const { return 3; };
	double *writeToMem(double *mem) const;
	double *readFromMem(double *mem);
	CartesianCoord toCartesian(void) const;
	
	void constrainThetaPhi(void);
	double constrainAngle(const double angle) const;
};

// Position or vector in Cartesian coordinates
class CartesianCoord {
public:
	double	_x, _y, _z;
	CartesianCoord(void) : _x(0), _y(0), _z(0) {};
	CartesianCoord(double x, double y, double z) : _x(x), _y(y), _z(z) {};
	
	size_t size(void) const { return 3; };
	double *writeToMem(double *mem) const;
	double *readFromMem(double *mem);
	SphericalCoord toSpherical(void) const;
	CartesianCoord crossProduct(const CartesianCoord &b) const;
	double dist(const CartesianCoord &o) const {
		double xd=_x-o._x, yd=_y-o._y, zd=_z-o._z;
		return sqrt(xd*xd+yd*yd+zd*zd);
	};
	
	// Operations that return a new object
	const CartesianCoord operator+(const CartesianCoord &other) const;
	const CartesianCoord operator-(const CartesianCoord &other) const;
	const CartesianCoord operator*(const double &val) const;
	const CartesianCoord operator/(const double &val) const;
	
	// Operations that modify this object
	CartesianCoord & operator+=(const CartesianCoord &other);
	CartesianCoord & operator-=(const CartesianCoord &other);
	CartesianCoord & operator*=(const double &val);
	CartesianCoord & operator/=(const double &val);
	void operator=(const CartesianCoord &other) { _x = other._x; _y = other._y; _z = other._z; };
};

class CoordUV {
public:
	double u, v;
	CoordUV(void) : u(NAN), v(NAN) {};
	CoordUV(double new_u, double new_v) : u(new_u), v(new_v) {};
};

class BoundaryPoint {
public:
	CartesianCoord	cartesian_pt;
	SphericalCoord	spherical_pt;
	double			cos_theta;
	double			sin_theta;
	double			cos_phi;
	double			sin_phi;
	
	BoundaryPoint(void) : cartesian_pt(), spherical_pt(), cos_theta(NAN), sin_theta(NAN), cos_phi(NAN), sin_phi(NAN) {};
	BoundaryPoint(CartesianCoord cc, SphericalCoord sc, double cost, double sint, double cosf, double sinf) :
		cartesian_pt(cc), spherical_pt(sc), cos_theta(cost), sin_theta(sint), cos_phi(cosf), sin_phi(sinf) {};
	
	CartesianCoord cartesian(void) const { return cartesian_pt; };
};

class CapBoundary {
private:
	BoundaryPoint	bounds[4];
	
public:
	BoundaryPoint operator[] (unsigned int i) const { assert(i<4); return bounds[i]; };
	void setBoundary(unsigned int bnum, SphericalCoord sc);
	void setCartTrigBounds(unsigned int bnum, CartesianCoord cc, double cost, double sint, double cosf, double sinf);
};

class TriElemLinearShapeFunc {
private:
	double		sf[3][3];
public:
	TriElemLinearShapeFunc(CoordUV xy1, CoordUV xy2, CoordUV xy3);
	double applyShapeFunc(const unsigned int i, const CoordUV uv) const;
};

class Tracer {
private:
	// Tracer position in spherical coordinates
	SphericalCoord	_sc;
	// Tracer position in Cartesian coordinates
	CartesianCoord	_cc;
	// Previous Cartesian position
	CartesianCoord	_cc0;
	// Previous Cartesian velocity
	CartesianCoord	_Vc;
	
	// Tracer flavor (meaning should be application dependent)
	double _flavor;
	
	// ID of element containing this tracer
	ElementID _ielement;
	
public:
	Tracer(void) : _sc(), _cc(), _cc0(), _Vc(), _flavor(0), _ielement(UNDEFINED_ELEMENT) {};
	Tracer(SphericalCoord new_sc, CartesianCoord new_cc) :
		_sc(new_sc), _cc(new_cc), _cc0(), _Vc(), _flavor(0), _ielement(UNDEFINED_ELEMENT) {};
	
	CartesianCoord getCartesianPos(void) const { return _cc; };
	SphericalCoord getSphericalPos(void) const { return _sc; };
	CartesianCoord getOrigCartesianPos(void) const { return _cc0; };
	CartesianCoord getCartesianVel(void) const { return _Vc; };
	
	void setCoords(const SphericalCoord new_sc, const CartesianCoord new_cc) {
		_sc = new_sc;
		_cc = new_cc;
	}
	void setOrigVals(const CartesianCoord new_cc0, const CartesianCoord new_vc) {
		_cc0 = new_cc0;
		_Vc = new_vc;
	}
	
	double theta(void) { return _sc._theta; };
	double phi(void) { return _sc._phi; };
	double rad(void) { return _sc._rad; };
	
	double x(void) { return _cc._x; };
	double y(void) { return _cc._y; };
	double z(void) { return _cc._z; };
	
	ElementID ielement(void) const { return _ielement; };
	void set_ielement(const ElementID ielement) { _ielement = ielement; };
	
	double flavor(void) const { return _flavor; };
	void set_flavor(const double flavor) { _flavor = flavor; };
	
	size_t size(void);
	double *writeToMem(double *mem) const;
	double *readFromMem(double *mem);
};

typedef std::list<Tracer> TracerList;

#endif

struct TRACE{

    FILE *fpt;

    char tracer_file[200];

    int itracer_warnings;
    int ianalytical_tracer_test;
    int ic_method;
    int itperel;
    int itracer_interpolation_scheme;

    double box_cushion;

    /* tracer arrays */
	
	// Sets of tracers organized by cap
    TracerList *tracers;
	
	// Sets of tracers that have escaped a cap, organized by cap
    TracerList *escaped_tracers;

    /* tracer flavors */
    int nflavors;
    int **ntracer_flavor[13];

	int number_of_tracers;
	
    int ic_method_for_flavors;
    double *z_interface;

    char ggrd_file[255];		/* for grd input */
    int ggrd_layers;

    /* statistical parameters */
    int istat_ichoice[13][5];
    int istat_isend;
    int istat_iempty;
    int istat1;
    int istat_elements_checked;
    int ilast_tracer_count;

	// Whether we are doing tracers in a full sphere shell or just a region
	bool full_tracers;
	
    /* timing information */
    double advection_time;
    double find_tracers_time;
    double lost_souls_time;

    /* Mesh information */
	CapBoundary	boundaries[13];

	// Map of node numbers to gnomonic coordinates
	std::map<int, CoordUV> *gnomonic;
    double gnomonic_reference_phi;
	
    /*********************/
    /* for global model  */
    /*********************/

    /* regular mesh parameters */
    int numtheta[13];
    int numphi[13];
    unsigned int numregel[13];
    unsigned int numregnodes[13];
    double deltheta[13];
    double delphi[13];
    double thetamax[13];
    double thetamin[13];
    double phimax[13];
    double phimin[13];
    int *regnodetoel[13];
    int *regtoel[13][5];

    /* gnomonic shape functions */
    TriElemLinearShapeFunc **shape_coefs[13][2];

    /**********************/
    /* for regional model */
    /**********************/

    double *x_space;
    double *y_space;
    double *z_space;

    /*********************/
    /* function pointers */
    /*********************/

    ElementID (* iget_element)(struct All_variables*, int, int,
                         CartesianCoord, SphericalCoord);

    CartesianCoord (* get_velocity)(struct All_variables*, int, int,
                          SphericalCoord);
};