Source: public/javascript/modules/utils.js

import { Constants } from "./constants.js";
import {getProjection} from "./olqv_projections.js";

/*
 ** A set of classes and function definitions utilized by the
 ** differents flavours of OLQV viewers.
 **
 ** Author : M. Caillat
 ** Date : 06th December 2018
 **        07th April 2021
 */

 let id = 0; 
 function getUniqueId() {
    return id++ + '';
  }

 class CoordinatesFormatter {
    static enter(what) {
        console.group(this.name + "." + what);
    }

    static exit() {
        console.groupEnd();
    }


    constructor(canvas, pixelValueUnit, projection) {
        CoordinatesFormatter.enter(this.constructor.name);
        this.canvas = canvas;
        this.pixelValueUnit = pixelValueUnit;
        this.dataSteps = null;
        this.projection = projection;
        CoordinatesFormatter.exit();
    }

    /**
     *
     * @param {integer[2]} olc pixel integer coordinates
     * @returns {String} a string built after the pixel integer and astronomical coordinates and pixel value/unit when possible
     */
    format(olc) {
        //CoordinatesFormatter.enter(this.format.name);
        var result;
        const radec = this.getRaDec(Math.round(olc[0]), Math.round(olc[1]));
        var ctx = this.canvas.getContext('2d');
        var pixelAtPosition = ctx.getImageData(olc[0], this.canvas.height - olc[1], 1, 1).data;
        if (pixelAtPosition) {
            var dataStepsIndex = pixelAtPosition.slice(0, 3).join('_');
            result = 'iRA=' + Math.round(olc[0]) + ' iDEC=' + Math.round(olc[1]) + ' RA=' + radec['ra'] + ' DEC=' + radec['dec'];
            if (this.dataSteps && this.dataSteps[dataStepsIndex]) {
                result += ` - Pixel value : ${this.dataSteps[dataStepsIndex].toExponential(4)} ${this.pixelValueUnit}`;
            }
        } else {
            result = "???";
        }
        return result;
    }

    /**
     * Returns RA/DEC values in HMS/DMS format corresponding to given x/y value
     *
     * @param {float} x 
     * @param {float} y 
     * @returns {dict}
     */
    getRaDec(x, y){
        return this.projection.iRaiDecToHMSDMS(x, y);
    }

    setDataSteps(dataSteps) {
        this.dataSteps = dataSteps;
    }

}


/*
 ** Converts a decimal number expected to represent an angle in degree
 ** into a string expressing a right ascension ( H:M:S)
 */
var DecDeg2HMS = function(deg, sep = ':') {
    //if(any(deg< 0 | deg>360)){stop('All deg values should be 0<=d<=360')}
    //if (deg < 0)
    //deg[deg < 0] = deg[deg < 0] + 360
    const HRS = Math.floor(deg / 15);
    const MIN = Math.floor((deg / 15 - HRS) * 60);
    let SEC = (deg / 15 - HRS - MIN / 60) * 3600;
    SEC = Math.floor(SEC * 1000) / 1000.;
    return HRS + sep + MIN + sep + SEC.toFixed(3);
};

/**
 * Converts a H:M:S value into decimal degrees
 * @param {*} value 
 * @returns 
 */
var HMS2DecDeg = function(value){
    const parts = value.split(":");
    if(parts.length !== 3){
        throw("Invalid format for H:M:S value");
    }else{
        const a = parseInt(parts[0]);
        const b = parseFloat(parts[1])*(1/60);
        const c = parseFloat(parts[2]) * (1/3600);
        return (a + b + c)*15;
    }
};

/*
 ** Converts a decimal number expected to represent an angle in degree
 ** into a string expressing a declination ( D:M:S)
 */
var DecDeg2DMS = function(deg, sep = ':') {
    var sign = deg < 0 ? '-' : '+';
    deg = Math.abs(deg);
    var DEG = Math.floor(deg);
    var MIN = Math.floor((deg - DEG) * 60);
    var SEC = (deg - DEG - MIN / 60) * 3600;
    SEC = Math.floor(SEC * 1000.) / 1000.;
    if (SEC < 0.) SEC = 0.;
    if (SEC > 60) SEC = 60.;

    return (sign + DEG + sep + MIN + sep + SEC.toFixed(3));
};

/**
 * Converts a declination D:M:S into decimal degrees
 * @param {*} value 
 * @returns result in degrees
 */
var DMS2DecDeg = function(value, sep=":"){
    const parts = value.split(sep);
    if(parts.length !== 3){
        throw("Invalid format for D:M:S value");
    }else if(isNaN(Number(parts[0])) || isNaN(Number(parts[1])) || isNaN(Number(parts[2]))){
        throw("Invalid D:M:S value");
    }else  {
        //can start by "-"
        const a = parseFloat(parts[0]);  
        const b = parseFloat(parts[1]) * (1/60);
        const c = parseFloat(parts[2]) * (1/3600);
        return (a + b + c);
    }
};

/*
 ** A function which returns the units of a 2D array resulting from
 ** the summation of a range of RA-DEC planes in a FITS cube.
 */
function summedPixelsUnit(unit) {
    switch (unit) {
        case "Jy/beam":
            return "Jy/beam * km/s";


        case "erg/s/cm^2/A/arcsec^2":
            return "erg/s/cm^2/arcsec^2";

        case "K (Ta*)":
            return "K.km/s";

        default:
            return "";
    }
}

/*
 ** A function which returns a factor
 ** for the display of physical quantity
 ** The returned values are based on experience rather
 ** than on a purely rational approach
 */
function unitRescale(unit) {
    switch (unit.toLowerCase()) {
        case "k (ta*)":
            return 1.0;

        case "jy/beam":
            return 1.0;

        case "erg/s/cm^2/a/arcsec^2":
            //     return 1e18;
            return 1.;

        case "erg/s/cm^2/a":
            //      return 1e12;
            return 1.;

        case "m/s":
            return 0.001;

        case "km/s":
            return 1.;

        case "mhz":
            return 1.;

        case "hz":
            return 0.000001;

        default:
            return 1.0;
    }
}

/*
 ** A function which sums the values of a array
 ** between two indices i0 (included ) and i1 ( excluded )
 ** and returns the sum multiplied by a coefficient coeff.
 */
function sumArr(arr, i0, i1, coeff) {
    i0 = Math.max(0, i0);
    i1 = Math.min(arr.length - 1, i1);
    if (i0 > i1)[i0, i1] = [i1, i0];

    try{
        return coeff * arr.slice(i0, i1 + 1).reduceRight(function(a, b) { return a + b; });
    }catch(exception){
        if(arr.slice(i0, i1 + 1).length === 0)
            alert(`Selected interval [${i0}, ${i1+1}] is out of bound`);
        else
            alert(exception);
    }
}

/**
 * Utility.
 * A function which parses a string into an array of floats. The string is expected
 * to be a comma separated list of textual representation of decimal numbers.
 * A range [min, max[ can be provided,  then the values are considered valif if and only if they lie in that range.
 *
 * @param {string} s a comma separated list of textual numbers representations
 * @param {number[2]} range if provided values are checked to lie in the range
 * 
 * @returns {number[]|undefined} an array of float numbers or undefined.
*/
function str2FloatArray(s, range = false) {
    let x = s.split(",");
    let w = undefined;
    let result = undefined;

    let y = x.map(function(z) { return parseFloat(z) });
    if (range) {
        w = y.map(function(t) { return (!isNaN(t) && (range[0] <= t) && (t < range[1])) });
    } else {
        w = y.map(function(t) { return true; });
    }

    if (w.reduce(function(accum, u) { return accum && u }, true)) {
        result = y;
    }

    return result;
}

/*

 ** A function which creates a document fragment out of an HTML string and appends it to the content of an existing element.
 ** The HTML string is assumed to describe a single element ( e.g. one signle div, p, etc. ).
 ** Returns the created element.
 */
function createAndAppendFromHTML(html, element) {
    var template = document.createElement('template');
    template.innerHTML = html.trim();
    $(element).append(template.content.cloneNode(true));
    return element.lastChild;
}

/**
 * 
 * @param {*} z 
 * @returns velocity in km/s
 */
function redshift2Velocity(z) {
    // Constants.SPEED_OF_LIGHT in m/s
    //return z * (Constants.SPEED_OF_LIGHT / 1000);
    return  Constants.SPEED_OF_LIGHT / 1000 * (z/(1+z))
}

/**
 * 
 * @param {number} v in km/s 
 * @returns redshift
 */
function velocity2Redshift(v) {
    // Constants.SPEED_OF_LIGHT in m/s
    //return v / (Constants.SPEED_OF_LIGHT / 1000);
    return v / ((Constants.SPEED_OF_LIGHT / 1000) - v);
}

/**
 * Returns standard deviation of an array
 * (https://codesource.io/how-to-do-standard-deviation-in-javascript/)
 * @param {array} arr 
 * @returns float
 */
function standardDeviation(arr){
    // removes null values from array before counting number of elements
    let tmp_array = arr.filter(val => val !== null);
    let mean = tmp_array.reduce((acc, curr)=>{
        return acc + curr
    }, 0) / tmp_array.length;    
    
    tmp_array = tmp_array.map((el)=>{
        return (el - mean) ** 2
    })
    
    let total = tmp_array.reduce((acc, curr)=> acc + curr, 0);
    
    return Math.sqrt(total / tmp_array.length)
}


/*
 ** Two functions to log when a function is entered and exited
 */
function ENTER() {
    var caller = ENTER.caller;
    if (caller == null) {
        result = "_TOP_";
    } else {
        result = caller.name + ": entering";
    }
    console.log(result + ": entering");
}

function EXIT() {
    var caller = EXIT.caller;
    if (caller == null) {
        result = "_TOP_";
    } else {
        result = caller.name + ": exiting";
    }
    console.log(result + ": exiting");
}

function inRange(x, xmin, xmax) {
    return ((x - xmin) * (x - xmax) <= 0);
}

/*
 
  ** Frequency <-> Velocity derivations.
  **
  ** From https://www.iram.fr/IRAMFR/ARN/may95/node4.html
  */

/*
  ** No verification is done on the values of restfreq and frequency. Both are expected to have realistic values.
  ** frequency and restfreq are expected to by in the same unit (i.e. both in HZ or both in GHZ)
  ** The result is in m/s.
  
 var f2v = function (frequency, restfreq, crval3) {
  return (Constants.SPEED_OF_LIGHT * (restfreq - frequency) / restfreq) + crval3;
}*/

var f2v = function(frequency, restfreq, vcenter) {
    //SPEED_OF_LIGHT in m/s
    return (Constants.SPEED_OF_LIGHT * (restfreq - frequency) / restfreq) + vcenter;
}


/*
 ** No verification is done on the values of restfreq and velocity. Both are expected to have realistic values.
 */
var v2f = function(velocity, restfreq, vcenter) {
    //SPEED_OF_LIGHT in m/s
    return restfreq * (1 - (velocity - vcenter) / Constants.SPEED_OF_LIGHT)
}


/*
 ** Frequency to wavelength
 */
var f2lambda = function(frequency) {
    return Constants.SPEED_OF_LIGHT / frequency;
}

/*
 ** Revert 1D array.
 */
var revertArray = function(a) {
    var result = [];
    for (let i = 0; i < a.length; i++) {
        result.push(a[a.length - i - 1]);
    }
    return result;
}

/*
 ** Round *original* number to *round* numbers after 0.
 */
var round = function(original, round) {
    var i = 0;
    var r = 1;
    while (i < round) {
        ++i;
        r *= 10;
    }
    return Math.round(original * r) / r;
};

/*
 ** Coordinate tabulators
 **
 **  - crval : value at reference point
 **  - cdelt : increment of abscissa
 **  - crpix : coordinate of reference point
 **  - n : tabulate at point of coordinate n
 */
var linearTabulator = function(crval, cdelt, crpix, n) {
    return crval + (n - crpix) * cdelt;
}

/**
 * Generates a random string 
 * (from https://stackoverflow.com/questions/7616461/generate-a-hash-from-string-in-javascript)
 * @param {*} s 
 * @returns 
 */
var hashCode = function(s) {
    return s.split("").reduce(function(a, b) { a = ((a << 5) - a) + b.charCodeAt(0); return a & a }, 0);
}

var degToArcsec = function(deg) {
    return deg * 3600;
}

function degToRad(deg) {
    return deg * (Math.PI / 180);
}

function radToDeg(rad) {
    return rad * ( 180 / Math.PI);
}

/**
 * get the rest value for a frequency according to a redshift or a velocity
 * returns a deep-copy of the array with shifted line values
 * @param {*} frequency 
 * @param {*} value  value used to calculate shift 
 * @param {*} type  redshift or velocity
 */
function unshift(frequency, velocity, velolsr) {
    let result = null;
    /*if (type == Constants.REDSHIFT_SHIFT_TYPE) {
        result = frequency * (1 + value);
    } else if (type == Constants.VELOCITY_SHIFT_TYPE) {*/
        result = frequency * (1 + (velocity / (Constants.SPEED_OF_LIGHT / 10 ** 3)));
    /*} else {
        result = frequency;
    }*/
    return result;
}

/**
 * calculate the shifted value of all the lines in the given array
 * returns a deep-copy of the array with shifted line values
 * @param {*} lines  array of lines
 * @param {*} velocity  value used to calculate shift in km/s
 * @param {*} velolsr  velolsr in km/s
 */
function shift(frequency, velocity, velolsr) {
    let result = null;
    let vel = velocity;
    /*if (type == Constants.REDSHIFT_SHIFT_TYPE) {
        result = frequency / (1 + value);
    } else if (type == Constants.VELOCITY_SHIFT_TYPE) {*/
    /*if(velolsr !== undefined){
        vel = vel - velolsr
    }*/
    result = frequency/* / (1 + (vel / (Constants.SPEED_OF_LIGHT / 10 ** 3)))*/;
    /*} else {
        alert("Unknown data type");
    }*/
    return result;
}

/**
 * Computes a shifted frequency value for a given z
 * @param {*} frequency 
 * @param {*} z 
 * @returns 
 */
function shiftFrequencyByZ(frequency, z){
    return frequency * (1+z);
}

/**
 * Computes a shifted frequency value for a given z
 * @param {*} frequency 
 * @param {*} z 
 * @returns 
 */
function unshiftFrequencyByZ(frequency, z){
    return frequency / (1+z);
}

/**
 * Computes a shifted frequency value for a given velocity in km/s
 * @param {*} frequency 
 * @param {*} z 
 * @returns 
 */
function shiftFrequencyByV(frequency, velocity){
    let result = frequency / (1 + (velocity / (Constants.SPEED_OF_LIGHT / 10 ** 3)));
    return frequency * (1 + (velocity / (Constants.SPEED_OF_LIGHT / 10 ** 3)));
}

/**
 * Remove shift from a frequency value with a given velocity in km/s
 * @param {*} frequency 
 * @param {*} z 
 * @returns 
 */
function unshiftFrequencyByV(frequency, velocity){
    return frequency / (1 + (velocity / (Constants.SPEED_OF_LIGHT / 10 ** 3)));
}


/**
 * Compute the search radius in degrees from coordinates
 * to get matching sources in NED
 * 
 * if radius_value is FOV : Returns the min value between abs(RAmax-RAmin) and abs(Decmax-Decmin)
 * else : return radius_value
 * 
 * radius_value : radius in arcmin or "fov"
 * fits_header : FITS_HEADER object
 */
function getRadiusInDegrees(radius_value, fits_header) {
    let proj = getProjection(fits_header.projectionType);
    if (radius_value === "fov") {
        let minvals = proj.iRaiDecToRaDec(0, 0);
        let maxvals = proj.iRaiDecToRaDec(fits_header.naxis1, fits_header.naxis2);

        // search radius
        let dist_in_ra = Math.abs(maxvals["ra"] - minvals["ra"]);
        let dist_in_dec = Math.abs(maxvals["dec"] - minvals["dec"]);

        return Math.min(dist_in_ra, dist_in_dec)/2;
    } else {
        return radius_value;
    }
}

/**
 * Returns value of L'_Line in K.km/s.pc2, then Mgas is a function of L'_Line
 * in case bunit is in Jy/beam
 * @param {*} sline    surface selected in summed pixel spectrum in Jy.km/s
 * @param {*} nu_obs   line observed frequency in GHz
 * @param {*} dl       distance in Mpc (from NED)
 * @param {*} z        redshift
 * @returns 
 */
function getLpLineJ(sline, nu_obs, dl, z){
    return 3.25e7 * sline * (nu_obs**-2) * (dl**2) * ((1+z)**-3);
}

/**
 * Returns value of L'_Line in K.km/s.pc2, then Mgas is a function of L'_Line
 * in case bunit is in K
 * 
 * @param {*} sline    surface selected in summed pixel spectrum in Jy.km/s
 * @param {*} dl       distance in Mpc (from NED)
 * @param {*} z        redshift
 * @param {*} cdelt1   
 * @param {*} cdelt2   
 * @param {*} bmin        
 * @param {*} bmaj        
 * @param {*} box      coordinates of selected box ([xmin, xmax, ymin, ymax])
 * @returns 
 */
function getLpLineK(sline, dl, z, cdelt1, cdelt2, bmin, bmaj, box){
    const npix = (box[3]-box[2]+1)*(box[1]-box[0]+1);
    const omeg1 = npix * Math.abs(cdelt1 * 3600 * cdelt2 * 3600);
    const omeg2 = bmin*3600*bmaj*3600;
    let omega = (omeg1**2 + omeg2 ** 2 ) ** 0.5;
    return 23.5 * omega * dl**2 * sline * (1+z)**-3/1e7
}

/**
*Converts K in cm-1 (for energy values)
*  @param {*} energy_in_K
*
*/
var KToCm = function(energy_in_K){
    return energy_in_K * 0.695;
  }
  
  /**
  *Converts cm-1 in K (for energy values)
*  @param {*} energy_in_cm
  */
  var cmToK = function(energy_in_cm){
    return energy_in_cm / 0.695;
  }

/**
 * Displays current value of variable (console.log sometimes show updated value)
 * @param {*} value 
 */
function printVal(value){
    console.log(JSON.parse(JSON.stringify(value)));
}

/**
 * Computes jansky per kelvin. 
 * @param {number} nu  frequency in Hz
 * @param {number} bmin  bmin in degrees
 * @param {number} bmaj  bmaj in degrees
 * @returns 
 */
function jyperk(nu, bmin, bmaj){
    const lambda = 2.99792458e8 / nu;
    const theta1 = bmin * (Math.PI / 180) / 2 / Math.sqrt(Math.log(2.0));
    const theta2 = bmaj * (Math.PI / 180) / 2 / Math.sqrt(Math.log(2.0));
    return  2.0 * 1.38e3 * Math.PI * theta1 * theta2 / lambda**2;
}

export {
    hashCode,
    degToRad,
    radToDeg,
    degToArcsec,
    linearTabulator,
    round,
    revertArray,
    /*f2lambda,*/ v2f,
    f2v,
    inRange,
    velocity2Redshift,
    redshift2Velocity,
    /*ENTER, EXIT, createAndAppendFromHTML, str2FloatArray,*/ sumArr,
    createAndAppendFromHTML,
    str2FloatArray,
    unitRescale,
    summedPixelsUnit,
    DecDeg2DMS,
    DecDeg2HMS,
    HMS2DecDeg,
    DMS2DecDeg,
    CoordinatesFormatter,
    getRadiusInDegrees,
    shift, 
    unshift,
    shiftFrequencyByZ,
    unshiftFrequencyByZ,
    shiftFrequencyByV,
    unshiftFrequencyByV,
    getLpLineJ,
    getLpLineK,
    KToCm,
    cmToK,
    printVal,
    standardDeviation,
    jyperk,  
    getUniqueId
}