/*
 * $Smake: cc -O -o %F %f -lm
 * 
 * xe-to-xe2.c
 * 
 * Jonathan Senning <senning@gordon.edu>
 * December 23, 2000
 * 
 * Compile with:
 * 
 *     cc -O -o xe-to-xe2 xe-to-xe2.c -lm
 * 
 * Usage (example):
 * 
 *     xe-to-xe2 < ppm.xe > ppm.xe2
 *
 * Description:
 * 
 * This program reads files in XE format and writes out the same data
 * converted to XE2 format.  XE is a raw data format used by XEphem
 * (http://www.clearskyinstitute.com/xephem/xephem.html) version 3.2.3
 * and earlier.  XE2 is a new format used on XEphem version 3.4.
 * 
 * The conversion done here was based on the information in the
 * xephem-3.2.3/GUI/xephem/ppm.c and xephem-3.4/GUI/xephem/xe2.c files
 * from the two distributions.  The new data file format reflects more
 * closely the format in the Tycho-2 data that became available early 2000.
 * The conversion is fairly straight-forward, but several points are worth
 * mentioning:
 * 
 * 1. Double precision floating point numbers are stored in 1 to 4 bytes,
 *    depending on the dynamic range of the numbers.  In some cases
 *    (e.g. the visual magnitude) the mapping between the floating-point
 *    version of the number and the byte(s) that is stored changed from
 *    XE format to XE2 format.  This means that numbers like the magnitude
 *    will be slightly different between XEphem versions using the XE format
 *    files and those using the XE2 format files.
 *
 * 2. The Right Ascension PPM correction information in the XE2 format
 *    contains a factor of cos(dec) where "dec" is the declination.  In
 *    the xephem-3.4/GUI/xephem/xe2.c file the value of "dec" that appears
 *    in this calculation has been adjusted with the PPM correction for
 *    declination.  In this program, however, the uncorrected value (i.e.
 *    correct for Epoch 2000) is used.  This may or may not be correct; I'm
 *    just not sure at this time.  As of December 2000 it seems to produce
 *    correct results.
 * 
 * No error checking is performed in this program.
 */

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

#define XE1PKTSZ	19
#define XE2PKTSZ	18

typedef unsigned char UC;
typedef unsigned long UL;

typedef UC XE1PKT[XE1PKTSZ];
typedef UC XE2PKT[XE2PKTSZ];

/*--------------------------------------------------------------------------*/

void convert( XE1PKT in, XE2PKT out )
{
    UL t;
    double dec;

    /*
     * Type
     * Encoding is the same, but location in first byte of record has changed.
     */
    out[0] = ( in[0] >> 3 ) & 0x38;	/* Select bits 3,4 and 5 */

    /*
     * Assume PPM or Tycho unless data is SAO or HD
     */
    t = ( (UL) in[3] << 16 ) | ( (UL) in[4] << 8 ) | (UL) in[5];
    if ( t != 0L )
    {
	out[0] = out[0] | ( ( in[3] & 0x80 ) ? 1 /* SAO */ : 2 /* HD */ );
	out[1] = 0;
	out[2] = in[3] & 0x7f;
	out[3] = in[4];
	out[4] = in[5];
    }
    else
    {
	/* PPM or Tycho */
	t = ( (UL) in[0] << 16 ) | ( (UL) in[1] << 8 ) | (UL) in[2];
	
	/*
	 * If the PPM, SAO and HD numbers are all zero then we don't
	 * really know the source, at least with the ppm.xe file that
	 * was available for XEphem 3.2.3.  Because I don't know a better
	 * way to handle this here, I'll give a type code that, according
	 * to the xe2.c file in XEphem 3.4, will assign the name "<Bogus>"
	 * to the star.
	 */
	if ( ( t & 0x3fffff ) == 0 ) out[0] = out[0] | 7 /* Bogus */;
	out[1] = 0;
	out[2] = in[0] & 0x3f;
	out[3] = in[1];
	out[4] = in[2];
    }
    
    /*
     * RA
     * Just need to change the location in the record.
     */
    out[5] = in[6];
    out[6] = in[7];
    out[7] = in[8];
    
    /*
     * Dec
     * Just need to change the location in the record.
     */
    out[8] = in[9];
    out[9] = in[10];
    out[10] = in[11];
    
    /*
     * Mag
     * Encoding has changed so we need to convert.  The "+ 0.5" is used
     * so the corresponding integer value is rounded not truncated.
     */
    out[11] = (UC) ( (double) in[12] * 255.0 / 190.0 + 0.5 );
    
    /*
     * Spectral class
     * Just need to change the location in the record.
     */
    out[12] = in[13];
    out[13] = in[14];

    /*
     * PM Dec
     * Encoding has changed.
     */
    t = ( (UL) in[17] << 8 ) | (UL) in[18];
    t = 2L * ( t - 10000L ) + 32768L;
    out[16] = (UC) ( ( t >> 8 ) & 0xff );
    out[17] = (UC) ( t & 0xff );

    /*
     * Compute Dec, which is needed for PM RA.
     */
    t = ( (UL) in[9] << 16 ) | ( (UL) in[10] << 8 ) | (UL) in[11];
    dec = t * M_PI / ( ( 1L << 24 ) - 1 ) - M_PI / 2.0;

    /*
     * PM RA
     * Encoding has changed.
     */
    t = ( (UL) in[15] << 8 ) | (UL) in[16];
    t = (UL) ( 3.0 * cos( dec ) * ( t - 5000.0 ) + 32768.0 + 0.5 );
    out[14] = (UC) ( ( t >> 8 ) & 0xff );
    out[15] = (UC) ( t & 0xff );
}

/*--------------------------------------------------------------------------*/

int main( void )
{
    XE1PKT xe1pkt;
    XE2PKT xe2pkt;

    write( 1, "XE2.\n", 5 );
    while ( XE1PKTSZ == read( 0, xe1pkt, XE1PKTSZ ) )
    {
	convert( xe1pkt, xe2pkt );
	write( 1, xe2pkt, XE2PKTSZ );
    }

    return 0;
}