/*
 * Common definitions
 * Made by Daniel E. Severin (CONICET) in 2015-2017
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

/*************************************
 * GENERAL DEFINITIONS AND FUNCTIONS *
 *************************************/
#define PI 3.1415926535897932384
#define ONEOVERSQRTTWOPI 1.0/sqrt(2.0*PI)
#define DISTTOL 2.0/60.0
#define MAGTOL 1.5
#define SIGMAALPHA 10.15/3600.0 /* sigma_{alpha*} in degrees */
#define SIGMADELTA 22.74/3600.0 /* sigma_{delta} in degrees */
#define SIGMAMAG 0.2759 /* sigma_{mag} in "magnitudes" */
#define WCS_J2000 1 /* J2000(FK5) right ascension and declination */
#define WCS_B1950 2 /* B1950(FK4) right ascension and declination */
extern "C" void wcsconp(int sys1, int sys2, double eq1, double eq2, double ep1, double ep2,
	double *dtheta, double *dphi, double *ptheta, double *pphi);
//SYNOPSIS:   wcsconp(sys1, sys2, eq1, eq2, ep1, ep2, dtheta, dphi, ptheta, pphi)
//  int     sys1;   /* Input coordinate system (J2000, B1950, ECLIPTIC, GALACTIC) */
//  int     sys2;   /* Output coordinate system (J2000, B1950, ECLIPTIC, GALACTIC) */
//  double  eq1;    /* Input equinox (default of sys1 if 0.0) */
//  double  eq2;    /* Output equinox (default of sys2 if 0.0) */
//  double  ep1;    /* Input Besselian epoch in years (for proper motion) */
//  double  ep2;    /* Output Besselian epoch in years (for proper motion) */
//  double  *dtheta; /* Right ascension in degrees: Input in sys1, returned in sys2 */
//  double  *dphi;  /* Declination in degrees: Input in sys1, returned in sys2 */
//  double  *ptheta; /* Right ascension proper motion in RA degrees/year: Input in sys1, returned in sys2 */
//  double  *pphi;  /* Declination proper motion in Dec degrees/year: Input in sys1, returned in sys2 */

/*
 * bye - show error and exit
 */
void bye(char *string)
{
	printf(string);
	exit(1);
}

/*
 * readfield - read a field from "initial" position (first = 1)
 */
void readfield(char *buffer, char *cell, int initial, int bytes)
{
	int length = strlen(buffer);
	if (initial > length) {
		for (int i = 0; i < bytes; i++) cell[i] = 0;
	}
	else {
		for (int i = 0; i < bytes; i++) cell[i] = buffer[i + initial - 1];
	}
	cell[bytes] = 0;
}

/*
 * Compute cos of angle in degrees
 */
double dcos(double angle)
{
	return cos(angle * PI / 180.0);
}

/*
 * Compute sin of angle in degrees
 */
double dsin(double angle)
{
	return sin(angle * PI / 180.0);
}

/*
 * Compute arc tan
 */
double datan2(double y, double x)
{
	if (x < 0.000001 && x > -0.000001) {
		/* caso x = 0 */
		if (y > 0.0) return 90.0;
		else return 270.0;
	}
	double angle = atan(y / x) * 180.0 / PI;
	if (x < 0.0) {
		if (y < 0.0) angle -= 180.0;
		else angle += 180.0;
	}
	if (angle < 0.0) angle += 360.0;
	return angle;
}

/*
 * Convert (RA,DECL) to (X,Y,Z)
 */
void sph2rec(double ra, double decl, double *x, double *y, double *z)
{
	*x = dcos(ra) * dcos(decl);
	*y = dsin(ra) * dcos(decl);
	*z = dsin(decl);
}

/*
 * calcscalar - compute scalar product between two rectangular coord. in the unit sphere
 */
double calcscalar(double x1, double y1, double z1, double x2, double y2, double z2)
{
	return x1*x2 + y1*y2 + z1*z2;
}

/*
 * calcdist - compute angular distance between two points, given scalar product of it (0 to 180 degrees)
 */
double calcdist(double cosdist)
{
	if (cosdist < -1.0) cosdist = -1.0;
	if (cosdist > 1.0) cosdist = 1.0;
	return acos(cosdist) * 180.0 / PI;
}

/*
 * formula for converting Johnson V to estimated visual mag. in CD scale
 * (according to quadratic fit due to paper Severin, Sevilla 2015)
 */
float convert_V(float vmag) {
	float A = -0.01335368, B = 1.076636, C = 0.2249828;
	return A*vmag*vmag + B*vmag + C;
}

/*
 * compute the difference alpha1* - alpha2*
 */
double alpha_difference(double alpha1, double alpha2, double delta) {
	double dif = fabs(alpha1 - alpha2);
	double dif1 = fabs(alpha1 - alpha2 + 360.0);
	if (dif1 < dif) dif = dif1;
	double dif2 = fabs(alpha1 - alpha2 - 360.0);
	if (dif2 < dif) dif = dif2;
	return dif*dcos(delta);
}

/*
 * pdf of normal distribution
 */
double pdfn(double x, double sigma) {
	float xoversigma = x / sigma;
	return ONEOVERSQRTTWOPI * exp(-0.5 * xoversigma * xoversigma) / sigma;
}

/************************
 * CD CATALOG STRUCTURE *
 ************************/

struct CDstar_struct {
	int decl_ref, num_ref; /* CD ident. */
	char suppl_ref; /* CD ident. suppl. (a,b,c) */
	float rah, ramin, raseg, decldeg, declmin; /* right ascension and declination, in parts for B1875 */
	double RA1875, DE1875, RA2000, DE2000; /* right ascension and declination for B1875 and J2000 (in degrees) */
	float mag; /* original estimated visual magnitude */
	bool dpl, color; /* = true if it is double doble, or color */
	int deleted; /* 0 or >0 if it should not be considered for the following reasons: 1) variable, 2) nebulae, 3) differs with PPM position or mag. */
	char reg[29]; /* complete register */
};

/*************************
* PPMX CATALOG STRUCTURE *
*************************/

struct PPMXstar_struct {
	char ppmx_name[16]; /* name of star */
	double RA1875, DE1875, RA2000, DE2000; /* right ascension and declination for B1875 and J2000 (in degrees) */
	float cmag, vmag, mag; /* magnitudes in C, Johnson V and CD scale */
	bool Vmag_assigned; /* false if Vmag is absent and not found in APASS */
};

/*************************
* STRUCTURE OF SETS P_a  *
*************************/

struct Pa_struct {
	int cardinal; /* number of elements of P_a */
	int *b1, *b2; /* arrays of "cardinal" elements with elements of B: 
				   *    emptyset -> b1 = -1, b2 = -1
				   * single star -> b1 in B, b2 = -1
				   * double star -> b1 in B, b2 in B */
	long double *weight; /* array of "cardinal" elements with weights */
	int chosen; /* the one chosen after cross identification */
};