/*
 *  $Id: shape.c 769 2006-01-02 20:03:23Z denneau $
 *  SHAPE
 *
 *  Original Fortran code by Mikko Kaasalainen (2005)
 *  C-code by Tommy Grav on 7/27/05.
 *  Copyright 2005. All rights reserved.
 *  For use in Pan-STARRS Moving Object Pipeline System
 */

// Computes the magnitude of an object based on a tri-axial ellipsoid model rotating in space

/* Input variables
 *   nrows - number of rows that make up the triaxial ellipsoid shape of the object (most often 8-12)
 *   orbfile - file containing the orbital and physical paramters of the file. Single line format:
 *             q e incl node peri T_peri epoch H G period T_rot amplitude axis_a axis_b axis_c lat long phase_1 phase_2 phase_3
 *   ephfile - file containing the observatioal data from the survey. List of observations, each line:
 *             t_obs ra dec mag_jpl
 *
 * Output variables:
 *   List of detections, each line:
 *   t_obs  
 */

#include <getopt.h>
#include "shape.h"
#include "mopsephem.h"

// Global variables
int debug = 0;


void usage(void) {
    /*! \brief Display usage message */
    char *msg = "\
usage: shape [--debug] inputfile outputfile\n\
\n\
  --debug : enable debug mode\n\
  inputfile : contains orbit and shape def on single line\n\
  outputfile : output filename\n\
";
    fprintf(stderr, msg);
}


static int get_opts(
        int argc,
        char **argv,
        char **infile,
        char **outfile,
        char **station)
    {

    int c;
    int option_index = 0;
    while (1) {
        static struct option long_opts[] = {
            {"debug", no_argument, NULL, 0},
            {0, 0, 0, 0},
        };

        c = getopt_long(argc, argv, "", long_opts, &option_index);
        if (c == -1)
            break;

        switch(c) {
            case 0:
                if (strcmp(long_opts[option_index].name, "debug") == 0) {
                    debug = 1;
                }
                else if (strcmp(long_opts[option_index].name, "obscode") == 0) {
                    *station = optarg;
                }
                break;

            default:
                usage();
                exit(100);
        }
    }

    if (optind < argc) {
        *infile = argv[optind++];
    }
    if (optind < argc) {
        *outfile = argv[optind++];
    }

    return 0;
}


int
MAIN__(int argc, char **argv)
{

	int i, j;					// Counters 
	int nmax;
	int ch_vertices;			// To check that number of vertices generated in ell2conv() is correct 
	int ch_facets;				// To check that number of facets generated in trifac() is correct  
	double rmagnitude;
	double cmagnitude;
	double h;					// The absolute magnitude of the object
	double ell[3];				// The three semi_major axis of the ellipsoid
	double phase[4];			// The three Kaasalainen phase function parameter and the phase angle (p[3]) 
	double sinlat;
	double len_ps, len_pe, ps_dot_pe;
	double absmag;
	double beta, lambda;
	double dec, ra, ra1, ra2;
	double hcorr;
	double scattering[2];
	double angpos_earth[2], angpos_sun[2];
	double ast_pos_earth[3], ast_pos_sun[3];
	double rotstate[4];
	double pos_earth[3], pos_sun[3], epoch, mjd;
	double ini_rotstate[4];
	double t_start, t_end, dt;
	double a, g, p, amp;
	double diam;
    char objectName[100];
	ORBIT_BLK orb;
	XYZ_BLK xyz;
	EPHEM_BLK obs;
	char label[256];

	char *station = "568";
    char *infilename = NULL;
	FILE *in_fp = NULL;
    char *outfilename = NULL;
    FILE *out_fp = NULL;
	int nrows = 10;

	VERTEX *ifp;
	double *vertex_pos;

    /* Setup. */
    get_opts(argc, argv, &infilename, &outfilename, &station);
    if (NULL == infilename || NULL == outfilename) {
        usage();        ///< not all inputs specified
        exit(100);
    }


	/* Initializations. */
#define MAX_VERTICES 100000
	ifp = (VERTEX *) malloc(MAX_VERTICES * sizeof(VERTEX));
	vertex_pos = (double *) malloc(MAX_VERTICES * 3 * sizeof(double));
	mopsephem_init();

	in_fp = fopen(infilename, "r");
	if (in_fp == NULL) {
		fprintf(stderr, "Error opening input file %s\n", infilename);
		exit(-1);
	}

	out_fp = fopen(outfilename, "w");
	if (out_fp == NULL) {
		fprintf(stderr, "Error opening output file %s\n", outfilename);
		exit(-1);
	}

    /* Read in the orbital and shape parameters */
    int go = 1;
    int num = 0;
    while (go) {
        if ((num = fscanf(in_fp, "%lf", &orb.q_au)) <= 0) {
            go = 0;
            break;
        }

        fscanf(in_fp, "%lf", &orb.e);
        fscanf(in_fp, "%lf", &orb.i_deg);
        fscanf(in_fp, "%lf", &orb.node_deg);
        fscanf(in_fp, "%lf", &orb.arg_peri_deg);
        fscanf(in_fp, "%lf", &orb.time_peri_mjd);
        fscanf(in_fp, "%lf", &orb.h_v);             // Absolute magnitude in V
        fscanf(in_fp, "%lf", &orb.t_epoch_mjd);
        fscanf(in_fp, "%10s", &objectName);         // objectname

        fscanf(in_fp, "%lf", &g);                   // g; Phase angle parameter G
        fscanf(in_fp, "%lf", &p);                   // p; Albedo
        fscanf(in_fp, "%lf", &rotstate[2]);         // period_d; Rotational period in days
        fscanf(in_fp, "%lf", &amp);                 // amp_mag; Maximum amplitude 
        fscanf(in_fp, "%lf", &ell[0]);              // a_km; Triaxial diameter
        fscanf(in_fp, "%lf", &ell[1]);              // b_km;
        fscanf(in_fp, "%lf", &ell[2]);              // c_km;
        fscanf(in_fp, "%lf", &rotstate[0]);         // beta_deg; Rotational pole orientation
        fscanf(in_fp, "%lf", &rotstate[1]);         // lambda_deg
        fscanf(in_fp, "%lf", &phase[0]);            // pha_a
        fscanf(in_fp, "%lf", &phase[1]);            // pha_d
        fscanf(in_fp, "%lf", &phase[2]);            // pha_k

        fscanf(in_fp, "%lf", &epoch);               // epoch for calculation

        if (debug) {
            printf("# orbit %6.2lf %6.2lf %6.2lf\n", orb.q_au, orb.e, orb.i_deg);
            printf("#       %6.2lf %6.2lf %12.2lf\n", orb.node_deg, orb.arg_peri_deg,
                   orb.time_peri_mjd);
            printf("# H,G,p %6.2lf %6.2lf %6.2lf\n", orb.h_v, g, p);
            printf("# shape %3d %6.2lf %6.2lf %6.2lf %10.4lf\n", nrows, ell[0],
                   ell[1], ell[2], amp);
            printf("# phy   %10.5lf %10.5lf %10.5lf\n", rotstate[0], rotstate[1],
                   rotstate[2]);
            printf("# phase %6.2lf %6.2lf %6.2lf\n", phase[0], phase[1], phase[2]);
        }


        int nvertices = 4 * nrows * nrows + 2;
        int nfacets = 8 * nrows * nrows;

        int facet_vertices[nfacets][3];
        double longitude[nvertices];
        double latitude[nvertices];

        // Setting up input variable
        // convert input pole direction from ecliptic to from pole
        lambda = rotstate[0] * DEG2RAD;
        beta = rotstate[1] * DEG2RAD;
        rotstate[1] =
            PI / 2. - asin(sin(beta) * cos(EPSILON) +
                           cos(beta) * sin(EPSILON) * sin(lambda));
        //printf(" GG %lf %lf %lf %lf GG ",lambda,beta,EPSILON,dec) ;
        ra1 = cos(lambda) * cos(beta);
        ra2 = -sin(beta) * sin(EPSILON) + cos(beta) * cos(EPSILON) * sin(lambda);
        rotstate[0] = atan2(ra2, ra1);
        if (debug) {
            printf("# equatorial rotation state: %6.2lf  %6.2lf\n",
               rotstate[0] * 180. / PI, rotstate[1] * 180. / PI);
        }

        // Initializing global variables 
        memset(ifp, 0, 1000 * sizeof(ifp));
#if 0
        for (i = 0; i < 1000; i++) {
            for (j = 0; j < 3; j++) {
                ifp[i].v[0] = 0.;
                ifp[i].v[1] = 0.;
                ifp[i].v[2] = 0.;
            }
        }
#endif

        // Vertices 
        ch_vertices = ell2conv(nrows, longitude, latitude);
        //printf("Vertices:\n") ;
        //printf("%3d %3d \n",ch_vertices,nvertices) ;
        for (i = 0; i < nvertices; i++) {
            sinlat = sin(latitude[i]);
            vertex_pos[i * 3 + 0] = (ell[0] / 2.) * sinlat * cos(longitude[i]);
            vertex_pos[i * 3 + 1] = (ell[1] / 2.) * sinlat * sin(longitude[i]);
            vertex_pos[i * 3 + 2] = (ell[2] / 2.) * cos(latitude[i]);
            //printf("#%3d  %8.2lf %8.2lf - %8.2lf %8.2lf %8.2lf \n",
            //   i,latitude[i],longitude[i],vertice_pos[i][0],vertice_pos[i][1],vertice_pos[i][2]) ;
        }

        // Facets
        ch_facets = trifac(ifp, nrows);
        //printf("Facets:\n") ;
        //printf("%3d %3d \n",ch_facets,nfacets) ;
        for (i = 0; i < nfacets; i++) {
            //printf("#%3d - ",i) ;
            for (j = 0; j < 3; j++) {
                facet_vertices[i][j] = ifp[i].v[0];
                facet_vertices[i][j] = ifp[i].v[1];
                facet_vertices[i][j] = ifp[i].v[2];
                //printf("[%8.2lf , %8.2lf] %3d ",latitude[facet_vertices[i][j]],longitude[facet_vertices[i][j]],ifp[i][j]) ; 
            }
            //printf("\n") ;
        }
        //Find the inital offset to get the absolute magnitude
        ini_rotstate[0] = rotstate[0];
        ini_rotstate[1] = rotstate[1];
        ini_rotstate[2] = rotstate[2];
        ini_rotstate[3] = 0.;

        diam = pow(ell[0] * ell[1] * ell[2], 1. / 3.);

        for (i = 0; i < nvertices; i++) {
            vertex_pos[i * 3 + 0] *= diam / ell[0];
            vertex_pos[i * 3 + 1] *= diam / ell[1];
            vertex_pos[i * 3 + 2] *= diam / ell[2];
        }

        pos_earth[0] = 1.;
        pos_earth[1] = 0.;
        pos_earth[2] = 0.;
        pos_sun[0] = 1.;
        pos_sun[1] = 0.;
        pos_sun[2] = 0.;
        ecl2ast(ini_rotstate, pos_earth, pos_sun, epoch, ast_pos_earth,
                ast_pos_sun);

        absmag = conv2lc(ifp, vertex_pos, nfacets, scattering, ast_pos_earth, ast_pos_sun);
        if (debug) {
            printf("# %lf \n", absmag);
        }
        hcorr = orb.h_v - absmag;
        if (debug) {
            printf("#\n");
        }

        for (i = 0; i < nvertices; i++) {
            vertex_pos[i * 3 + 0] *= ell[0] / diam;
            vertex_pos[i * 3 + 1] *= ell[1] / diam;
            vertex_pos[i * 3 + 2] *= ell[2] / diam;
        }

        // Get the position of the earth and the target at time mjd
        mopsephem_compute_xyz(&orb, &xyz, epoch, station);
        mopsephem_compute_ephemeris(&orb, &obs, epoch, station);

        pos_earth[0] = xyz.ex - xyz.ox;
        pos_earth[1] = xyz.ey - xyz.oy;
        pos_earth[2] = xyz.ez - xyz.oz;
        pos_sun[0] = -xyz.ox;
        pos_sun[1] = -xyz.oy;
        pos_sun[2] = -xyz.oz;

        ps_dot_pe =
            pos_earth[0] * pos_sun[0] + pos_earth[1] * pos_sun[1] +
            pos_earth[2] * pos_sun[2];
        len_ps =
            sqrt(pos_sun[0] * pos_sun[0] + pos_sun[1] * pos_sun[1] +
                 pos_sun[2] * pos_sun[2]);
        len_pe =
            sqrt(pos_earth[0] * pos_earth[0] + pos_earth[1] * pos_earth[1] +
                 pos_earth[2] * pos_earth[2]);
        phase[3] = acos(ps_dot_pe / (len_ps * len_pe)) * 180. / PI;

        ecl2ast(rotstate, pos_earth, pos_sun, epoch, ast_pos_earth,
                ast_pos_sun);
        rmagnitude =
            conv2lc(ifp, vertex_pos, nfacets, scattering, ast_pos_earth,
                    ast_pos_sun) + hcorr;

                // Note that the following function also writes output.
        cmagnitude =
            magnitude_corrections(rmagnitude, pos_earth, pos_sun, phase);
        //    printf("\n") ;

        // Real output here.  Woohoo.
        fprintf(out_fp, "%s ", objectName);
        fprintf(out_fp, "%.10f ", epoch);
        fprintf(out_fp, "%.10f %10.5f ", obs.ra_deg, obs.dec_deg);
        fprintf(out_fp, "%.10f ", obs.mag) ;
        fprintf(out_fp, "%.10f ", cmagnitude);
        fprintf(out_fp, "\n");

    }

	close(in_fp);
	close(out_fp);
}
