// grid.cc -- classes for a 2D uniform tensor-product grid
// $Id$
//
// grid_arrays::grid_arrays
//
// grid::grid
// grid::print_gridfn
//

#include <stdio.h>
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <float.h>

#include "jt/stdc.h"
#include "jt/util++.hh"		// jtutil::how_many_in_range(),
				// round<fp>::, fuzzy<fp>::
#include "jt/array.hh"
#include "jt/linear_map.hh"

#include "fp.hh"
#include "coords.hh"
#include "grid.hh"

//*****************************************************************************
//*****************************************************************************
//*****************************************************************************

//
// This function constructs a  grid_arrays  object.
//
grid_arrays::grid_arrays(const grid_array_pars& grid_array_pars_in,
			 int N_gridfns_in)
	: gridfn_data_(0, N_gridfns_in-1,
		       grid_array_pars_in.min_irho
			- grid_array_pars_in.min_rho_N_ghost_zones,
		       grid_array_pars_in.max_irho
			+ grid_array_pars_in.max_rho_N_ghost_zones,
		       grid_array_pars_in.min_isigma
			- grid_array_pars_in.min_sigma_N_ghost_zones,
		       grid_array_pars_in.max_isigma
			+ grid_array_pars_in.max_sigma_N_ghost_zones),

	  N_gridfns_(N_gridfns_in),

	  min_irho_(grid_array_pars_in.min_irho),
	  max_irho_(grid_array_pars_in.max_irho),
	  min_isigma_(grid_array_pars_in.min_isigma),
	  max_isigma_(grid_array_pars_in.max_isigma),

	  full_grid__min_irho_(grid_array_pars_in.min_irho
				- grid_array_pars_in.min_rho_N_ghost_zones),
	  full_grid__max_irho_(grid_array_pars_in.max_irho
				+ grid_array_pars_in.max_rho_N_ghost_zones),
	  full_grid__min_isigma_(grid_array_pars_in.min_isigma
				  - grid_array_pars_in.min_sigma_N_ghost_zones),
	  full_grid__max_isigma_(grid_array_pars_in.max_isigma
				  + grid_array_pars_in.max_sigma_N_ghost_zones)
								     // no comma
{ }

//*****************************************************************************
//*****************************************************************************
//*****************************************************************************

//
// This function constructs a  grid  object.
//
grid::grid(const grid_array_pars& grid_array_pars_in,
	   const grid_pars& grid_pars_in,
	   int N_gridfns_in)

	: grid_arrays(grid_array_pars_in,
		      N_gridfns_in),

	  rho_map_(grid_array_pars_in.min_irho
		    - grid_array_pars_in.min_rho_N_ghost_zones,
		   grid_array_pars_in.max_irho
		    + grid_array_pars_in.max_rho_N_ghost_zones,
		   jtutil::radians_of_degrees(
		      grid_pars_in.min_drho
		       - grid_array_pars_in.min_rho_N_ghost_zones
			 * grid_pars_in.delta_drho
					     ),
		   jtutil::radians_of_degrees(grid_pars_in.delta_drho),
		   jtutil::radians_of_degrees(
		      grid_pars_in.max_drho
		       + grid_array_pars_in.max_rho_N_ghost_zones
			 * grid_pars_in.delta_drho
					     )),

	  sigma_map_(grid_array_pars_in.min_isigma
		      - grid_array_pars_in.min_sigma_N_ghost_zones,
		     grid_array_pars_in.max_isigma
		      + grid_array_pars_in.max_sigma_N_ghost_zones,
		     jtutil::radians_of_degrees(
			grid_pars_in.min_dsigma
			 - grid_array_pars_in.min_sigma_N_ghost_zones
			   * grid_pars_in.delta_dsigma
					       ),
		     jtutil::radians_of_degrees(grid_pars_in.delta_dsigma),
		     jtutil::radians_of_degrees(
			grid_pars_in.max_dsigma
			 + grid_array_pars_in.max_sigma_N_ghost_zones
			   * grid_pars_in.delta_dsigma
					       )),

	  min_rho_(jtutil::radians_of_degrees(grid_pars_in.min_drho)),
	  max_rho_(jtutil::radians_of_degrees(grid_pars_in.max_drho)),
	  min_sigma_(jtutil::radians_of_degrees(grid_pars_in.min_dsigma)),
	  max_sigma_(jtutil::radians_of_degrees(grid_pars_in.max_dsigma))
								     // no comma
{ }

//******************************************************************************

//
// This function prints a gridfn in ASCII format to an already-open
// stdio stream.  The output format is suitable for a gnuplot 'splot'
// command.
//
void grid::print_gridfn(int gfn, FILE *output_fp = stdout)
{
fprintf(output_fp, "# gfn=%d\n", gfn);

fprintf(output_fp, "# drho\t\tdsigma\t\tgridfn\n");
	for (int irho = min_irho() ; irho <= max_irho() ; ++irho)
	{
	const fp drho = drho_of_irho(irho);
		for (int isigma = min_isigma() ;
		     isigma <= max_isigma() ;
		     ++isigma)
		{
		const fp dsigma = dsigma_of_isigma(isigma);
		fprintf(output_fp,
			"%g\t%g\t%.15g\n",
			drho, dsigma, gridfn(gfn, irho,isigma));
		}
	printf("\n");
	}
}