#include "ppSim.h"

bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config) {

    bool mdok;
    assert (stars);

    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe

    bool starsReal = psMetadataLookupBool(&mdok, recipe, "STARS.REAL"); // Density of fakes
    if (!starsReal) return true;

    // Read catalogue stars using psastro
    psMetadata *astroRecipe = psMetadataLookupPtr(NULL, config->recipes, PSASTRO_RECIPE);
    if (!astroRecipe) {
        psError(PSASTRO_ERR_CONFIG, false, "Can't find recipe %s", PSASTRO_RECIPE);
        return NULL;
    }

    float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point
    float ra0 = psMetadataLookupF32(NULL, config->arguments, "RA"); // Boresight RA (radians)
    float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians)
    float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians)
    float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
    float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
    float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time

    // Size of FPA
    psRegion *bounds = ppSimFPABounds (fpa);
    float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scale;
    psFree(bounds);

    float x0fpa = 0.5*(bounds->x0 + bounds->x1);
    float y0fpa = 0.5*(bounds->y0 + bounds->y1);

    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MIN",  PS_DATA_F32 | PS_META_REPLACE, "", ra0 - radius);
    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MAX",  PS_DATA_F32 | PS_META_REPLACE, "", ra0 + radius);
    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MIN", PS_DATA_F32 | PS_META_REPLACE, "", dec0 - radius);
    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MAX", PS_DATA_F32 | PS_META_REPLACE, "", dec0 + radius);
    psArray *refStars = psastroLoadRefstars(config);
    if (!refStars || refStars->n == 0) {
        psError(PS_ERR_UNKNOWN, false, "Unable to find reference stars.");
        psFree(refStars);
        return NULL;
    }
    psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", refStars->n);

    long oldSize = stars->n;
    stars = psArrayRealloc (stars, refStars->n);

    // Conversion loop
    for (long i = 0; i < refStars->n; i++) {
        ppSimStar *star = ppSimStarAlloc ();

        pmAstromObj *ref = refStars->data[i]; // Reference star
        star->ra  = ref->sky->r; // RA of star
        star->dec = ref->sky->d; // Dec of star
        star->mag = ref->Mag;       // Magnitude of star

        // Convert to x,y position on tangent plane, in pixels
        float div = (sin(star->ra) * sin(dec0) + cos(star->dec) * cos(dec0) * cos(star->ra - ra0)) * scale;
        float xi = cos(star->dec) * sin(star->ra - ra0) / div;
        float eta = (sin(star->dec) * cos(dec0) - cos(star->dec) * sin(dec0) * cos(star->ra - ra0)) / div;

        // Apply rotation, make FPA center of boresite
        star->x = cos(pa) * xi - sin(pa) * eta + x0fpa;
        star->y = sin(pa) * xi + cos(pa) * eta + y0fpa;

        // Convert magnitude to peak flux
        star->flux = powf(10.0, -0.4 * (star->mag - zp)) * expTime;
        star->peak = star->flux / (2.0*M_PI * PS_SQR(seeing));
        stars->data[oldSize + i] = star;

        psTrace("ppSim", 10, "Adding catalogue star: %.1f,%.1f --> %.2f\n", star->x, star->y, star->flux);
    }
    stars->n = oldSize + refStars->n;
    psFree(refStars);

    psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_META_REPLACE, "Right ascension", ra0);
    psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_META_REPLACE, "Declination", dec0);

    return stars;
}
