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

#include "io-C/lisaxml.h"
#include "io-C/readxml.h"
#include "io-C/xmlbinary.h"

#include "LISAconstants.h"

void KILL(char*);

int main(int argc,char **argv) {
    writeXML *myxml;
    char buffer[256];
    char inputname[50];
    char outputname[50];
    char binaryname[50];
    char inname[50];
    char ErrorMessage[100];
    int cs;
    size_t ns;
    int i;
    double Ttot;
    int Ntot;
    double t, hp, hc;
    double Ap, Ac, cinc, cpsi, spsi, fgw, iphase;
    FILE* AmpsOut;
    GBParameters *gbparams;

  if (argc < 2)
  {
    printf("The correct call of this program is:\n\n");
    printf("   ./galactic_barycenter name\n");
    KILL("Please try again.\n");
  }

  strcpy(inname, argv[1]);

  printf("\nProcessing %s\n", inname);

   sprintf(inputname, "XML/%s.xml", inname);

   sprintf(outputname, "XML/%s_Barycenter.xml", inname);
   sprintf(binaryname, "Binary/%s.bin", inname);

    gbparams = getGBParameters(inputname);

    Ntot = NFFT+2*(int)ceil(Tpad/dt);

    Ttot = dt*(double)(Ntot);
  
    cinc = cos(gbparams->Inclination);

  Ap = (gbparams->Amplitude)*(1.0+pow(cinc,2.0));
  Ac = -2.0*(gbparams->Amplitude)*cinc;

  cpsi = cos(2.0*gbparams->Polarization);
  spsi = sin(2.0*gbparams->Polarization);

  fgw = gbparams->Frequency;
  iphase = gbparams->InitialPhase;

    if ((AmpsOut = fopen(binaryname, "wb")) == NULL)
    {
      sprintf(ErrorMessage, "Cannot open the file %s.\n", binaryname);
      KILL(ErrorMessage);
    }

  /* to use centered times: t = -Tpad+dt/2.0; */

  t = -Tpad+dt/2.0;

  for (i = 0 ; i < Ntot ; i++)
  {

   /* Calculate the polarized signals h_plus and h_cross */

    hp = Ap*cpsi*cos(2.0*pi*fgw*t+iphase) 
          + Ac*spsi*sin(2.0*pi*fgw*t+iphase);
   
    hc = -Ap*spsi*cos(2.0*pi*fgw*t+iphase) 
           + Ac*cpsi*sin(2.0*pi*fgw*t+iphase);

    /* hp and hc are interleaved */
      fwrite(&hp, sizeof(double), 1, AmpsOut);
      fwrite(&hc, sizeof(double), 1, AmpsOut); 

   /* Increment the time */

    t += dt;

   /* User update */

    if (i == (int)((double)Ntot/5.0)) printf("  20 percent done\n");
    if (i == 2*(int)((double)Ntot/5.0)) printf("  40 percent done\n");
    if (i == 3*(int)((double)Ntot/5.0)) printf("  60 percent done\n");
    if (i == 4*(int)((double)Ntot/5.0)) printf("  80 percent done\n");
    if (i == Ntot-1) printf("  100 percent done\n\n");
  }

  fclose(AmpsOut);
    
    myxml = XMLopen(outputname);

    /* Standard beginning of the file */

    XMLcontentstring(myxml,"<?xml version=\"1.0\"?>");
    XMLcontentstring(myxml,"<!DOCTYPE XSIL SYSTEM \"http://www.vallis.org/lisa-xml.dtd\">");

    XMLcontentstring(myxml,"<?xml-stylesheet type=\"text/xsl\" href=\"http://www.vallis.org/lisa-xml.xsl\"?>");

    XMLopentag(myxml,"XSIL","");

        /* Prolog */

        XMLopentag(myxml,"Param","Name=\"Author\"");
            XMLcontentstring(myxml,"Neil Cornish");
        XMLclosetag(myxml,"Param");
    
        XMLopentag(myxml,"Param","Name=\"GenerationDate\" Type=\"ISO-8601\"");
            getISOtime(buffer);
            XMLcontentstring(myxml,buffer);
        XMLclosetag(myxml,"Param");

        XMLopentag(myxml,"Comment","");
            XMLcontentstring(myxml,"This file generated by galactic_barycenter.c (NJC 2006/06/06)");
        XMLclosetag(myxml,"Comment");

        /* LISAData section (not required for simulator input files) */

        /* XMLopentag(myxml,"XSIL","Type=\"LISAData\"");

        XMLopentag(myxml,"XSIL","Name=\"LISA\" Type=\"PseudoLISA\"");

        XMLparamdouble(myxml,"InitialPosition","Radian",kappa);

        XMLparamdouble(myxml,"InitialRotation","Radian",lambda);

        XMLclosetag(myxml,"XSIL");

        XMLclosetag(myxml,"XSIL"); */

        /* SourceData section */

        XMLopentag(myxml,"XSIL","Type=\"SourceData\"");

            XMLopentag(myxml,"XSIL","Name=\"%s\" Type=\"SampledPlaneWave\"",gbparams->Name);

                /* Give parameters as needed */

                XMLparamdouble(myxml,"EclipticLatitude","Radian",gbparams->EclipticLatitude);
                XMLparamdouble(myxml,"EclipticLongitude","Radian",gbparams->EclipticLongitude);

                /* Since polarization is already included in the hp and hc */

                XMLparamdouble(myxml,"Polarization","Radian",0.0);

                XMLopentag(myxml,"XSIL","Name=\"%s\" Type=\"TimeSeries\"","hp,hc");

                    XMLparamdouble(myxml,"TimeOffset","Second",-Tpad);
                    XMLparamdouble(myxml,"Cadence","Second",dt);
                    XMLparamdouble(myxml,"Duration","Second",Ttot);

                    if (testbyteorder() == BIGENDIAN) {
                        XMLarray(myxml,"hp,hc","double","1",Ntot,2,"Remote","Binary,BigEndian",binaryname);
                    } else {
                        XMLarray(myxml,"hp,hc","double","1",Ntot,2,"Remote","Binary,LittleEndian",binaryname);                
                    }

                XMLclosetag(myxml,"XSIL");

                XMLparamstring(myxml,"Interpolator","String","Lagrange");
                XMLparamstring(myxml,"InterpolatorWindow","1","2");

            XMLclosetag(myxml,"XSIL");
            
        XMLclosetag(myxml,"XSIL");

    XMLclosetag(myxml,"XSIL");

    XMLclose(myxml);
    
    freeGBParameters(gbparams);
}

void KILL(char* Message)
{
  printf("\a\n");
  printf(Message);
  printf("Terminating the program.\n\n\n");
  exit(1);

 
  return;
}