import { FITS_HEADER } from "./fitsheader.js";
import {DecDeg2DMS, DecDeg2HMS, degToRad} from "./utils.js";
let ProjEnum = new Enum(['AZP', 'SZP', 'TAN', 'STG', 'SIN', 'ARC', 'ZPN', 'ZEA', 'AIR'
, 'CYP', 'CEA', 'CAR', 'MER', 'COP', 'COE', 'COD', 'COO', 'SFL', 'PAR'
, 'MOL', 'AIT', 'BON', 'PCO', 'TSC', 'CSC', 'QSC', 'HPX', 'XPH', 'GLS']);
//p_azimutal pour RA--ARC : cl.fits, p_radio RA--GLS : m33, p_ortho RA--SIN (mrc.fits), p_gnomonic RA--TAN (nenufar.fits)
class ProjUtils {
/**
* Applies arcsin function and returns an angle in radian
* @param {number} x (float)
* @returns {float}
*/
static arcsin(x) {
var result;
if (x < -1.e0) {
result = -Math.PI * 0.5e0;
} else if (x > 1.e0) {
result = pi * 0.5e0
} else {
result = Math.asin(x);
}
return result;
}
/**
* Applies arccos function and returns an angle in radian
* @param {number} x (float)
* @returns {float}
*/
static arccos(x) {
var result;
if (x < - 1.e0) {
result = Math.PI
} else if (x > 1.e0) {
result = 0.0e0;
} else {
result = Math.acos(x);
}
return result;
}
}
/**
* Base class for oject providing projection functions
*/
class Projection {
static precision = 1e-10;
static print = false;
/**
* Factory function creating a projection object
* @param {number} a0 crval1 in radians (float)
* @param {number} d0 crval2 in radians (float)
* @param {number} angle rotation angle (float)
* @param {ProjEnum} type type of projection
* @returns
*/
static create(a0, d0, angle, type) {
let result = null;
if(!ProjEnum.isDefined(type)){
throw('Undefined projection type');
}
switch (type) {
case ProjEnum.TAN:
result = new GnomonicProjection(a0, d0, angle);
break;
case ProjEnum.SIN:
result = new OrthoGraphicProjection(a0, d0, angle);
break;
case ProjEnum.ARC:
result = new AzimutalProjection(a0, d0, angle);
break;
case ProjEnum.GLS:
result = new RadioProjection(a0, d0, angle);
break;
default:
throw "Projection of type " + type.toString() + " is not yet implemented.";
}
return result;
}
/**
* @constructor
* @param {number} a0 reference coordinate on axis 1 in radians (float)
* @param {number} d0 reference coordinate on axis 2 in radians (float)
* @param {number} angle rotation in radians (float)
*/
constructor(a0, d0, angle) {
this._a0 = a0;
this._d0 = d0;
this._angle = angle;
this._sind0 = Math.sin(d0);
this._cosd0 = Math.cos(d0);
this._sina0 = Math.sin(a0);
this._cosa0 = Math.cos(a0);
}
/**
* Returns RA/Dec in degrees from a coordinate in pixels
* @param {number} ra pixel coordinates on ra axis (int)
* @param {number} dec pixel coordinates on dec axis (int)
* @returns {Object} Object containing ra, dec coordinates
*/
iRaiDecToRaDec(iRa, iDec) {
const xi = (iRa - FITS_HEADER.crpix1) * FITS_HEADER.cdelt1 * Math.PI / 180;
const yi = (iDec - FITS_HEADER.crpix2) * FITS_HEADER.cdelt2 * Math.PI / 180;
// pixels to coordinates in radians
const radec = this.relToAbs([xi], [yi]);
//const radec = this.relToAbs([0.005274772851999852], [-0.014117774397999605]);
/*console.log("### iRaiDecToRaDec");
console.log("iRA : "+iRa + ", iDec :" + iDec);
console.log("radians");
console.log("ra : " + radec['ra'] + " dec : " + radec['dec']);*/
// ra, dec in degrees
let ra = radec['ra'] / Math.PI * 180;
let dec = radec['dec'] / Math.PI * 180;
/*console.log("degrees");
console.log("ra : " + ra + " dec : " + dec);*/
return { 'ra': ra, 'dec': dec };
}
/**
* Returns a coordinate in pixels from RA/DEC in degrees
* @param {number} ra RA in degrees (float)
* @param {number} dec DEC in degreess (float)
* @returns {Object} Object containing iRa, iDec coordinates
*/
raDecToiRaiDec(ra, dec) {
// absToRel expects radians
const iRaiDec = this.absToRel([degToRad(ra)], [degToRad(dec)]);
const x = iRaiDec['x'];
const y = iRaiDec['y'];
const i_xp=Math.round(FITS_HEADER.crpix1+x/(degToRad(FITS_HEADER.cdelt1))-1);
const i_yp=Math.round(FITS_HEADER.crpix2+y/(degToRad(FITS_HEADER.cdelt2))-1);
// pixels to coordinates in radians
return { 'iRa': i_xp, 'iDec': i_yp };
}
/**
* returns RA/Dec in HMS/DMS
* @param {number} ra pixel coordinates on ra axis (int)
* @param {number} dec pixel coordinates on dec axis (int)
* @returns {Object} Object containing ra, dec coordinates
*/
iRaiDecToHMSDMS(iRa, iDec) {
let raDec = this.iRaiDecToRaDec(iRa, iDec);
/*console.log("### iRaiDecToHMSDMS");
console.log("RA : " + DecDeg2HMS(raDec['ra']));
console.log("DEC : " + DecDeg2DMS(raDec['dec']));*/
return { 'ra': DecDeg2HMS(raDec['ra']), 'dec': DecDeg2DMS(raDec['dec']) };
}
/**
* RA is between 0 and 2PI
* @param {number} ra RA in radians (float)
* @returns {number} (float)
*/
moduloRa(ra) {
if (ra > 2 * Math.PI) {
return ra - 2 * Math.PI;
} else if (ra < 0) {
return ra + 2 * Math.PI;
} else {
return ra;
}
}
/**
* DEC is between -PI/2 and PI/2
* @param {number} dec DEC in radians (float)
* @returns {number} (float)
*/
moduloDec(dec) {
if (dec > Math.PI / 2) {
return dec - Math.PI;
} else if (dec < -Math.PI / 2) {
return dec + Math.PI;
} else {
return dec;
}
}
/**
* Computes the absolute (RA, DEC) from the relative linear offset given in the header
*
* @abstract
* @param {array} x array of offsets in X
* @param {array} y array of offsets in Y
* @returns {Object} an object containing ra/dec coordinates
*/
relToAbs(x, y) {
//Projection.enter(this.relToAbs.name);
//Projection.exit();
throw 'Use this method on a derived class !'
}
/**
* Computes the relative offset (xp, yp) in a linear plane, from an absolute (RA, DEC).
* Use xp and yp to retrieve the i_yp, i_yp pixel values
* @abstract
* @param {array} a right ascensions (float)
* @param {array} d declinations (float)
* @returns {Object} array of offsets x and y
*/
absToRel(a, d) {
//Projection.enter(this.absToRel.name);
//Projection.exit();
throw 'Use this method on a derived class !'
}
}
/**
* Class representing a Radio projection
* @extends Projection
*/
class RadioProjection extends Projection {
/**
* @constructor
* @param {number} a0 reference coordinate on axis 1 in radians (float)
* @param {number} d0 reference coordinate on axis 2 in radians (float)
* @param {*} angle
*/
constructor(a0, d0, angle) {
super(a0, d0, angle);
if (angle != 0.0) {
//alert("RADIO does not support a projection angle, angle is set to 0.");
//angle = 0.0;
throw("RadioProjection does not support a projection angle, value must be 0");
}
this._npole = Math.PI * 0.5 - this._d0;
this._spole = -this.d0 - Math.PI * 0.5;
}
/**
* Computes the absolute (RA, DEC) from the relative linear offset given in the header
*
* @abstract
* @param {array} x array of offsets in X
* @param {array} y array of offsets in Y
* @returns {Object} an object containing ra/dec coordinates
*/
relToAbs(x, y) {
const n = x.length;
let a = new Array(n);
let d = new Array(n);
for (let i = 0; i < n; i++) {
if (y[i] <= this._spole) {
a[i] = this._a0;
d[i] = -Math.PI * 0.5;
} else if (y[i] >= this._npole) {
a[i] = this._a0;
d[i] = Math.PI * 0.5;
} else {
const p = this._d0 + y[i];
d[i] = p;
const r = x[i] / Math.cos(p);
if (r <= -Math.PI) {
a[i] = this._a0 - Math.PI;
} else if (r >= Math.PI) {
a[i] = this._a0 + Math.PI;
} else {
a[i] = this._a0 + r;
}
}
}
for(let i = 0; i < n; i++){
a[i] = this.moduloRa(a[i]);
d[i] = this.moduloDec(d[i]);
}
return { 'ra': a, 'dec': d };
}
/**
* Computes the relative offset (xp, yp) in a linear plane, from an absolute (RA, DEC).
* Use xp and yp to retrieve the i_yp, i_yp pixel values
* @abstract
* @param {array} a right ascensions (float)
* @param {array} d declinations (float)
* @returns {Object} array of offsets x and y
*/
absToRel(a, d) {
const n = a.length;
let x = new Array(n);
let y = new Array(n);
for (let i = 0; i < n; i++) {
const p = d[i];
let da = a[i] - this._a0;
while (da < -Math.PI) {
da = da + 2. * Math.PI;
}
while (da > Math.PI) {
da = da - 2. * Math.PI;
}
x[i] = da * Math.cos(p);
y[i] = p - this._d0;
}
return { 'x': x, 'y': y };
}
}
/**
* Class representing a Gnomonic projection
* @extends Projection
*/
class GnomonicProjection extends Projection {
/**
* @constructor
* @param {number} a0 reference coordinate on axis 1 in radians (float)
* @param {number} d0 reference coordinate on axis 2 in radians (float)
* @param {*} angle
*/
constructor(a0, d0, angle) {
super(a0, d0, angle);
if (this._d0 < -1.e-30) {
this._spole = 1. / Math.tan(this._d0);
this._npole = 1.e38;
}
else if (this._d0 > 1.e-30) {
this._npole = 1. / Math.tan(this._d0);
this._npole = 1.e38;
}
else {
this._npole = 1.e38;
this._spole = 1.e38;
}
}
/**
* Computes the absolute (RA, DEC) from the relative linear offset given in the header
*
* @abstract
* @param {array} x array of offsets in X
* @param {array} y array of offsets in Y
* @returns {Object} an object containing ra/dec coordinates
*/
relToAbs(x, y) {
let n = x.length;
let a = new Array(n);
let d = new Array(n);
if (this._sind0 > Projection.precision) {
for (let i = 0; i < n; i++) {
const r = Math.atan(Math.sqrt(x[i] ** 2 + y[i] ** 2));
const sinr = Math.sin(r);
const cosr = Math.cos(r);
if (r > Projection.precision) {
const p = Math.atan2(x[i], y[i]) - this._angle;
const sinp = Math.sin(p);
const cosp = Math.cos(p);
const sindec = (this._sind0 * cosr + this._cosd0 * sinr * cosp);
if (y[i] < this._npole) {
d[i] = Math.asin(sindec);
a[i] = this._a0 + Math.asin(sinr * sinp / Math.cos(d[i]));
} else {
d[i] = Math.asin(sindec);
a[i] = this._a0 - Math.asin(sinr * sinp / Math.cos(d[i])) + Math.PI;
}
} else {
a[i] = this._a0;
d[i] = this._d0;
}
}
} else if (this._sind0 < -Projection.precision) {
for (let i = 0; i < n; i++) {
const r = Math.atan(Math.sqrt(x[i] ** 2 + y[i] ** 2));
const sinr = Math.sin(r);
const cosr = Math.cos(r);
if (r > Projection.precision) {
const p = Math.atan2( x[i],y[i] ) - this._angle;
const sinp = Math.sin(p);
const cosp = Math.cos(p);
const sindec = (this._sind0 * cosr + this._cosd0 * sinr * cosp);
if (y[i] > this._spole) {
d[i] = Math.asin(sindec);
a[i] = this._a0 + Math.asin(sinr * sinp / Math.cos(d[i]));
} else {
d[i] = Math.asin(sindec);
a[i] = this._a0 - Math.asin(sinr * sinp / Math.cos(d[i])) + Math.PI;
}
} else {
a[i] = this._a0;
d[i] = this._d0;
}
}
} else {
for (let i = 0; i < n; i++) {
r = Math.atan(Math.sqrt(x[i] ** 2 + y[i] ** 2));
const sinr = Math.sin(r);
const cosr = Math.cos(r);
if (r > Projection.precision) {
const p = Math.atan2(x[i], y[i]) - this._angle;
const sinp = Math.sin(p);
const cosp = Math.cos(p);
const sindec = (this._sind0 * cosr + this._cosd0 * sinr * cosp);
d[i] = Math.asin(sindec);
a[i] = this._a0 + Math.asin(sinr * sinp / mth.cos(d[i]));
} else {
a[i] = this._a0;
d[i] = this._d0;
}
}
}
//Projection.exit();
return { 'ra': a, 'dec': d };
}
/**
* Computes the relative offset (xp, yp) in a linear plane, from an absolute (RA, DEC).
* Use xp and yp to retrieve the i_yp, i_yp pixel values
* @abstract
* @param {array} a right ascensions (float)
* @param {array} d declinations (float)
* @returns {Object} array of offsets x and y
*/
absToRel(a, d) {
const n = a.length;
let x = new Array(n)
let y = new Array(n);
//console.log("### absToRel");
//console.log("a0 : " + this._a0);
//console.log("sind0 : " + this._sind0);
//console.log("cosd0 : " + this._cosd0);
//console.log("angle : " + this._angle);
if (Math.abs(this._sind0) >= Projection.precision) {
for (let i = 0; i < n; i++) {
const sindec = Math.sin(d[i]);
const cosdec = Math.cos(d[i]);
const sina = Math.sin(a[i] - this._a0);
const cosa = Math.cos(a[i] - this._a0);
const product = sindec * this._cosd0 - cosdec * this._sind0 * cosa;
const r = Math.acos((sindec - this._cosd0 * product) / this._sind0);
//console.log("product : " + product);
//console.log("r : " + r);
if (r > Projection.precision) {
const p = Math.atan2(cosdec * sina, product) + this._angle;
const expr = Math.tan(r);
x[i] = expr * Math.sin(p);
y[i] = expr * Math.cos(p);
} else {
x[i] = 0.0;
y[i] = 0.0;
}
}
} else {
for (let i = 0; i < n; i++) {
const sina = Math.sin(a[i] - this._a0);
const cosa = Math.cos(a[i] - this._a0);
const cosdec = Math.cos(d[i]);
const sindec = Math.sin(d[i]);
const r = Math.acos(cosdec * cosa);
if (r >= Projection.precision) {
const p = Math.atan2(cosdec * sina, sindec) + this._angle;
const expr = Math.tan(r);
x[i] = expr * Math.sin(p);
y[i] = expr * Math.cos(p);
} else {
x[i] = 0.0;
y[i] = 0.0;
}
}
}
//Projection.exit();
return { 'x': x, 'y': y };
}
}
/**
* Class representing an OrthoGraphicProjection projection
* @extends Projection
*/
class OrthoGraphicProjection extends Projection {
constructor(a0, d0, angle) {
super(a0, d0, angle);
if (this._d0 < 1.e-30) {
this._spole = -this._cosd0;
this._npole = 1.e38;
}
else if (this._d0 > 1.e-30) {
this._npole = this._cosd0;
this._spole = 1.e38;
}
else {
this._npole = 1.;
this._spole = -1.;
}
}
/**
* Computes the absolute (RA, DEC) from the relative linear offset given in the header
*
* @abstract
* @param {array} x array of offsets in X
* @param {array} y array of offsets in Y
* @returns {Object} an object containing ra/dec coordinates
*/
relToAbs(x, y) {
//Projection.enter(this.relToAbs.name);
let n = x.length;
let a = new Array(n);
let d = new Array(n);
if (this._sind0 > Projection.precision) {
for (let i = 0; i < n; i++) {
let r = x[i] ** 2 + y[i] ** 2;
if (r < 1) {
const sinr = Math.sqrt(r);
r = Math.asin(sinr);
if (sinr > Projection.precision) {
const p = Math.atan2(x[i], y[i]) - this._angle;
const cosr = Math.cos(r);
const sinp = Math.sin(p);
const cosp = Math.cos(p);
const sindec = (this._sind0 * cosr + this._cosd0 * sinr * cosp);
d[i] = Math.asin(sindec);
a[i] = this._a0 + Math.atan2(sinr * sinp * this._sind0, sindec * this._cosd0 - sinr * cosp);
} else {
a[i] = this._a0;
d[i] = this._d0;
}
} else {
a[i] = 0.;
d[i] = 0.;
}
}
} else if (this._sind0 < -Projection.precision) {
for (let i = 0; i < n; i++) {
let r = x[i] ** 2 + y[i] ** 2;
if (r < 1) {
const sinr = Math.sqrt(r);
r = Math.asin(sinr);
if (sinr > Projection.precision) {
const p = Math.atan2(x[i], y[i]) - this._angle;
const cosr = Math.cos(r);
const sinp = Math.sin(p);
const cosp = Math.cos(p);
const sindec = (this._sind0 * cosr + this._cosd0 * sinr * cosp);
d[i] = Math.asin(sindec);
a[i] = this._a0 + Math.atan2(sinr * sinp * this._sind0, sindec * this._cosd0 - sinr * cosp) + Math.PI;
} else {
a[i] = this._a0;
d[i] = this._d0;
}
} else {
a[i] = 0.;
d[i] = 0.;
}
}
} else {
for (let i = 0; i < n; i++) {
let r = x[i] ** 2 + y[i] ** 2;
if (r < 1.) {
const sinr = Math.sqrt(r);
r = Math.asin(sinr);
if (sinr > Projection.precision) {
const p = Math.atan2(x[i], y[i]) - angle;
const cosr = Math.cos(r);
const sinp = Math.sin(p);
const cosp = Math.cos(p);
const sindec = (sind0 * cosr + cosd0 * sinr * cosp);
d[i] = Math.asin(sindec);
a[i] = a0 + Math.atan2(sinr * sinp, cosr + sind0 * sinr * cosp);
} else {
a[i] = a0;
d[i] = d0;
}
} else {
a[i] = 0.;
d[i] = 0.;
}
}
}
for(let i = 0; i < n; i++){
a[i] = this.moduloRa(a[i]);
d[i] = this.moduloDec(d[i]);
}
return { 'ra': a, 'dec': d };
};
/**
* Computes the relative offset (xp, yp) in a linear plane, from an absolute (RA, DEC).
* Use xp and yp to retrieve the i_yp, i_yp pixel values
* @abstract
* @param {array} a right ascensions (float)
* @param {array} d declinations (float)
* @returns {Object} array of offsets x and y
*/
absToRel(a, d) {
//Projection.enter(this.absToRel.name);
let n = a.length;
let x = new Array(n);
let y = new Array(n);
if (Math.abs(this._sind0) >= Projection.precision) {
for (let i = 0; i < n; i++) {
const sindec = Math.sin(d[i]);
const cosdec = Math.cos(d[i]);
const sina = Math.sin(a[i] - this._a0);
const cosa = Math.cos(a[i] - this._a0);
const diff = sindec * this._cosd0 - cosdec * this._sind0 * cosa;
const r = Math.acos((sindec - this._cosd0 * diff) / this._sind0);
if (r > Projection.precision) {
const p = Math.atan2(cosdec * sina, diff) + this._angle;
const expr = Math.sin(r);
x[i] = expr * Math.sin(p);
y[i] = expr * Math.cos(p);
} else {
x[i] = 0;
y[i] = 0;
}
}
} else {
for (let i = 0; i < n; i++) {
const sina = Math.sin(a[i] - this._a0);
const cosa = Math.cos(a[i] - this._a0);
const cosdec = Math.cos(d[i]);
const sindec = Math.sin(d[i]);
const r = Math.acos(cosdec * cosa);
if (r > Projection.precision) {
const p = Math.atan2(cosdec * sina, sindec) + this._angle;
const expr = Math.sin(r);
x[i] = expr * Math.sin(p);
y[i] = expr * Math.cos(p);
} else {
x[i] = 0;
y[i] = 0;
}
}
}
//Projection.exit();
return { 'x': x, 'y': y };
}
}
/**
* Class representing an AzimutalProjection projection
* @extends Projection
*/
class AzimutalProjection extends Projection {
constructor(a0, d0, angle) {
super(a0, d0, angle);
if (d0 < 0) {
this._spole = -d0 - Math.PI * 0.5;
this._npole = Math.PI * 0.5 + d0;
} else {
this._spole = d0 - Math.PI * 0.5;
this._npole = Math.PI * 0.5 - d0;
}
}
/**
* Computes the absolute (RA, DEC) from the relative linear offset given in the header
*
* @abstract
* @param {array} x array of offsets in X
* @param {array} y array of offsets in Y
* @returns {Object} an object containing ra/dec coordinates
*/
relToAbs(x, y) {
const n = x.length;
let a = new Array(n);
let d = new Array(n);
if (this._sind0 > Projection.precision) {
for (let i = 0; i < n; i++) {
const r = Math.sqrt(x[i] ** 2 + y[i] ** 2);
//if(Projection.print)
//console.log(`x:${x}, y:${y}, r:${r}`);
if (r <= Projection.precision) {
a[i] = this._a0;
d[i] = this._d0;
} else if (r <= Math.PI - Projection.precision) {
let p = Math.atan2(x[i], y[i]) - this._angle;
//if(Projection.print)
// console.log(`angle: ${this._angle}, p: ${p}`);
let sinr = Math.sin(r);
let cosr = Math.cos(r);
let sinp = Math.sin(p);
let cosp = Math.cos(p);
let sindec = (this._sind0 * cosr + this._cosd0 * sinr * cosp);
//if(Projection.print)
// console.log(`sindec:${sindec}, cosr:${cosr}, sinr:${sinr}, sind0:${this._sind0}, cosd0:${this._cosd0}, cosp:${cosp}`);
d[i] = Math.asin(sindec);
a[i] = this._a0 + Math.atan2(sinr * sinp * this._sind0, sindec * this._cosd0 - sinr * cosp);
//if(Projection.print)
// console.log(`a[i]:${a[i]}, d[i]:${d[i]}`);
} else if (r < Math.PI + Projection.precision) {
a[i] = this._a0;
d[i] = -d0;
} else {
a[i] = 0;
d[i] = 0;
}
}
} else if (this._sind0 < -Projection.precision) {
for (let i = 0; i < n; i++) {
const r = Math.sqrt(x[i] ** 2 + y[i] ** 2);
let sinr = Math.sin(r);
let cosr = Math.cos(r);
if (r <= Projection.precision) {
a[i] = this._a0;
d[i] = this._d0;
} else if (r <= Math.PI - Projection.precision) {
let p = Math.atan2(x[i], y[i]) - this._angle;
let sinp = Math.sin(p);
let cosp = Math.cos(p);
let sindec = (this._sind0 * cosr + this._cosd0 * sinr * cosp);
d[i] = Math.asin(sindec);
a[i] = this._a0 + Math.atan2(sinr * sinp * this._sind0, sindec * this._cosd0 - sinr * cosp) + Math.PI;
} else if (r < Math.PI + Projection.precision) {
a[i] = this._a0;
d[i] = -this._d0;
} else {
a[i] = 0;
d[i] = 0;
}
}
} else {
for (let i = 0; i < n; i++) {
const r = Math.sqrt(x[i] ** 2 + y[i] ** 2);
let sinr = Math.sin(r);
let cosr = Math.cos(r);
if (r <= Projection.precision) {
a[i] = this._a0;
d[i] = this._d0;
} else if (r <= Math.PI - Projection.precision) {
let p = Math.atan2(x[i], y[i]) - this._angle;
let sinp = Math.sin(p);
let cosp = Math.cos(p);
let sindec = (this._sind0 * cosr + this._cosd0 * sinr * cosp);
d[i] = Math.asin(sindec);
a[i] = this._a0 + Math.atan2(sinr * sinp, cosr + this._sind0 * sinr * cosp);
} else if (r < Math.PI + Projection.precision) {
a[i] = this._a0;
d[i] = -this._d0;
} else {
a[i] = 0;
d[i] = 0;
}
}
}
for(let i = 0; i < n; i++){
a[i] = this.moduloRa(a[i]);
d[i] = this.moduloDec(d[i]);
}
//Projection.exit();
return { 'ra': a, 'dec': d };
}
/**
* Computes the relative offset (xp, yp) in a linear plane, from an absolute (RA, DEC).
* Use xp and yp to retrieve the i_yp, i_yp pixel values
* @abstract
* @param {array} a right ascensions (float)
* @param {array} d declinations (float)
* @returns {Object} array of offsets x and y
*/
absToRel(a, d) {
const n = a.length;
let x = new Array(n);
let y = new Array(n);
if (Math.abs(this._sind0) >= Projection.precision) {
for (let i = 0; i < n; i++) {
let sindec = Math.sin(d[i]);
let cosdec = Math.cos(d[i]);
let sina = Math.sin(a[i] - this._a0);
let cosa = Math.cos(a[i] - this._a0);
let expr = sindec * this._cosd0 - cosdec * this._sind0 * cosa;
let r = Math.acos((sindec - this._cosd0 * expr) / this._sind0);
if (r > Projection.precision) {
const p = Math.atan2(cosdec * sina, expr) + this._angle;
x[i] = r * Math.sin(p);
y[i] = r * Math.cos(p);
} else {
x[i] = 0;
y[i] = 0;
}
}
} else {
for (let i = 0; i < n; i++) {
let sina = Math.sin(a[i] - this._a0);
let cosa = Math.cos(a[i] - this._a0);
let cosdec = Math.cos(d[i]);
let sindec = Math.sin(d[i]);
let r = acos(cosdec * cosa);
if (r >= Projection.precision) {
const p = Math.atan2(cosdec * sina, sindec) + this._angle;
x[i] = r * Math.sin(p);
y[i] = r * Math.cos(p);
} else {
x[i] = 0;
y[i] = 0;
}
}
}
//Projection.exit();
return { 'x': x, 'y': y };
}
}
/**
* Returns a Projection object according to the type passed in parameter
* @param {string} type
* @returns {Projection}
*/
function getProjection(type) {
if(!ProjEnum.isDefined(type)){
throw("Projection type is not implemented : " + type);
}else{
const proj_type = ProjEnum.get(type);
let pc1_1 = FITS_HEADER.pc1_1 != undefined ? FITS_HEADER.pc1_1 : 1.;
let pc2_1 = FITS_HEADER.pc2_1 != undefined ? FITS_HEADER.pc2_1 : 0.;
const ROTA2 = Math.atan(pc2_1 / pc1_1);
let angle;
if(FITS_HEADER.crota2 !== undefined){
angle=FITS_HEADER.crota2 * Math.PI/180.;
}else{
angle = ROTA2;
}
//angle = 0;
const CRVAL1 = FITS_HEADER.crval1 * Math.PI / 180;
const CRVAL2 = FITS_HEADER.crval2 * Math.PI / 180;
let projection = Projection.create(CRVAL1, CRVAL2, angle, proj_type);
return projection;
}
}
export {
ProjUtils, Projection, GnomonicProjection, OrthoGraphicProjection,
AzimutalProjection, RadioProjection, getProjection, ProjEnum
}