/********************************************************/
/* code for finding the spec for each of the star  */
/* within the catalogue                            */
/* compilation:                                    */
/* gcc -g -O0 -Wall -I/usr/include -I/usr/include/cfitsio -L/lib64 -lcfitsio -lm main.c -o main */
/* usage: */
/* code/main file.par cat.list spec/Ext_odonnell94_R3.1_CK04W.fits */
/*********************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <fitsio.h>
#include <time.h>
#include "utils.c"
// #include "csst_spec_interp.h"


//#include <stdio.h>
//#include <stdlib.h>
//#include <math.h>
#define MAXSTRLEN 200
//#define G_const 6.67384E-8
#define lgG log10(6.67384E-8)
#define lgMsun log10(1.9891E33)
#define lgExp log10(exp(1.))
#define twolgpc2cm 2.*log10(3.086e+18)

struct SPECLIB
{
  float Z;
  char name[MAXSTRLEN];
  int nspec;
  int nwv;
  float *logts;
  float *UniqueLogts;
  int *iLogts;
  int NULogts;
  float *loggs;
  float *UniqueLoggs;
  int *iLoggs;
  int NULoggs;
  float *wvs;
  float **specs;
};

struct STAR
{
  float logte;
  float logg;
  //float logL;
  float Mass;
  float Av;
  float mu0;
  float Z;
  //float mbolmag;
};

void replace(char * o_string, char * s_string, char * r_string);
void SortedUniqueElements(float in_array[], float out_array[], int size, int *NU);
int interp1D (float x, float *v, int *n, int nmax);
float locateitin1Dtab (float x, float *v, int n, int *i1, int *i2,
		       int flagrange, int flaglimits, int flaglog);

float locateitin1Dtab(float, float *, int, int *, int *, int, int, int);

void printerror(int);

int loadSpecLibs(char *, char *);
void loadExts(char *);
void interpSingleStar(struct STAR, float *, float *);

int NZ=0,NWV=0;
struct SPECLIB *speclibs=NULL;
float *Exts=NULL;
//float *outspec=NULL;


void printerror( int status)
{
  /*****************************************************/
  /* Print out cfitsio error messages and exit program */
  /*****************************************************/
  if (status)
    {
      fits_report_error(stderr, status); /* print error report */
      exit( status );    /* terminate the program, returning error status */
    }
  return;
}




int loadSpecLibs(char *file_par, char *specDir)
{
  fitsfile *fp_spec;
  FILE *fp_par;
  char tmpstr[MAXSTRLEN];
  int iZ,icol,imod,iwv, nspec_table, itable, icol_start, ispec, Nbyte=4;
  int status=0,nfound,hdutype,nullval,anynulls;
  long naxes[2];
  
  //load the spec libs
  fp_par=fopen(file_par,"r");
  fscanf(fp_par, "%d", &NZ);
  printf("%d different metallicity tables of spec lib\n",NZ);
  fgets(tmpstr, MAXSTRLEN, fp_par);
  speclibs=malloc(sizeof(struct SPECLIB)*NZ);
  for (iZ=0; iZ<NZ; iZ++)
    {
      fscanf(fp_par, "%f %s", &speclibs[iZ].Z, speclibs[iZ].name);
      printf("loading spec file: %s\n",speclibs[iZ].name);
      fgets(tmpstr, MAXSTRLEN, fp_par);
      char tmpName[MAXSTRLEN] = "";
      strcat(tmpName, specDir);
      strcat(tmpName, speclibs[iZ].name);
      printf("%s\n",tmpName);
      if (fits_open_file(&fp_spec, tmpName, READONLY, &status)) printerror(status);
      if (fits_movabs_hdu(fp_spec, 2, &hdutype, &status)) printerror(status);
      if (fits_read_keys_lng(fp_spec, "NAXIS", 1, 2, naxes, &nfound, &status)) printerror(status);
      printf("first table of each metallicity: naxes[0:1]: %ld, %ld\n", naxes[0], naxes[1]);
      speclibs[iZ].nspec=naxes[1];
      speclibs[iZ].logts=malloc(sizeof(float)*naxes[1]);
      speclibs[iZ].UniqueLogts=malloc(sizeof(float)*naxes[1]);
      speclibs[iZ].iLogts=malloc(sizeof(int)*naxes[1]);
      speclibs[iZ].loggs=malloc(sizeof(float)*naxes[1]);
      speclibs[iZ].UniqueLoggs=malloc(sizeof(float)*naxes[1]);
      speclibs[iZ].iLoggs=malloc(sizeof(int)*naxes[1]);
      icol=1;
      fits_read_col(fp_spec, TFLOAT, icol, 1, 1, naxes[1], &nullval, speclibs[iZ].logts, &anynulls, &status);
      for (imod=0; imod<speclibs[iZ].nspec; imod++) speclibs[iZ].logts[imod]=log10(speclibs[iZ].logts[imod]);
      fits_read_col(fp_spec, TFLOAT, icol+1, 1, 1, naxes[1], &nullval, speclibs[iZ].loggs, &anynulls, &status);

      SortedUniqueElements(speclibs[iZ].logts, speclibs[iZ].UniqueLogts, speclibs[iZ].nspec, &speclibs[iZ].NULogts);
      //printf("%d unique Teff grid points of spec lib\n", speclibs[iZ].NULogts);
      for (imod=0; imod<speclibs[iZ].nspec; imod++)
	{
	  for (speclibs[iZ].iLogts[imod]=0;
	       speclibs[iZ].iLogts[imod]<speclibs[iZ].NULogts;
	       speclibs[iZ].iLogts[imod]++)
	    {
	      if (speclibs[iZ].logts[imod]==speclibs[iZ].UniqueLogts[speclibs[iZ].iLogts[imod]])
		break;
	    }
	}
      SortedUniqueElements(speclibs[iZ].loggs, speclibs[iZ].UniqueLoggs, speclibs[iZ].nspec, &speclibs[iZ].NULoggs);
      //printf("%d unique logg grid points of spec lib\n", speclibs[iZ].NULoggs);
      for (imod=0; imod<speclibs[iZ].nspec; imod++)
	{
	  for (speclibs[iZ].iLoggs[imod]=0;
	       speclibs[iZ].iLoggs[imod]<speclibs[iZ].NULoggs;
	       speclibs[iZ].iLoggs[imod]++)
	    {
	      if (speclibs[iZ].loggs[imod]==speclibs[iZ].UniqueLoggs[speclibs[iZ].iLoggs[imod]])
		break;
	    }
	  //printf("speclibs[iZ].iLoggs[imod]: %d\n",speclibs[iZ].iLoggs[imod]);
	}




      speclibs[iZ].specs=(float **) malloc(sizeof(float *)*speclibs[iZ].nspec);
      ispec=0;
      nspec_table=ceil((speclibs[iZ].nspec+1.)/999.);
      printf("nspec_table=%d\n", nspec_table);
      for (itable=0; itable<nspec_table; itable++)
	{
	  printf("itable=%d\n", itable);
	  if (fits_movabs_hdu(fp_spec, 3+itable, &hdutype, &status)) printerror(status);
	  //if (fits_read_keys_lng(fp_spec, "NAXIS", 1, 2, naxes, &nfound, &status)) printerror(status);
	  if (fits_read_keys_lng(fp_spec, "NAXIS", 1, 2, naxes, &nfound, &status)) printerror(status);
	  printf("naxes[0:1]=%ld, %ld\n", naxes[0]/Nbyte, naxes[1]);
	  icol_start=1;
	  if (itable==0)
	    {
	      speclibs[iZ].nwv=naxes[1];
	      speclibs[iZ].wvs=malloc(sizeof(float)*naxes[1]);
	      //fits_read_col(fp_spec, TFLOAT, icol, 1, 1, speclibs[iZ].nwv, &nullval, speclibs[iZ].wvs, &anynulls, &status);
	      fits_read_col(fp_spec, TFLOAT, 1, 1, 1, speclibs[iZ].nwv, &nullval, speclibs[iZ].wvs, &anynulls, &status);
	      printf("wavelengths: speclibs[iZ].nwv, speclibs[iZ].wvs[0], speclibs[iZ].wvs[speclibs[iZ].nwv-1]: %d %f %f\n",
		     speclibs[iZ].nwv, speclibs[iZ].wvs[0], speclibs[iZ].wvs[speclibs[iZ].nwv-1]);
	      //printf("nspecs: %d\n",speclibs[iZ].nspec);
	      //if (iZ>1)	{      fits_close_file(fp_spec, &status);      continue;}
	      icol_start=2;
	    }

	  //for (icol=icol_start; icol<speclibs[iZ].nspec+2; icol++)
	  //for (icol=icol_start; icol<naxes[0]/4+1; icol++)
	  for (icol=icol_start; icol<naxes[0]/Nbyte+1; icol++)
	    {
	      //if (ispec==speclibs[iZ].nspec) break;
	      speclibs[iZ].specs[ispec] = (float *) malloc(sizeof(float)*speclibs[iZ].nwv);
	      //printf("add: %d",malloc(sizeof(float)*speclibs[iZ].nwv));
	      fits_read_col(fp_spec, TFLOAT, icol, 1, 1, speclibs[iZ].nwv,
			    &nullval, speclibs[iZ].specs[ispec], &anynulls, &status);
	      //printf("%f %f\n",speclibs[iZ].specs[icol-2][0],speclibs[iZ].specs[icol-2][speclibs[iZ].nwv-1]);
	      for (iwv=0; iwv < speclibs[iZ].nwv; iwv++)
		{
		  //speclibs[iZ].specs[icol-2][iwv] = (speclibs[iZ].specs[icol-2][iwv] < 1E-32) ?
		  //                                   -32 : log10(speclibs[iZ].specs[icol-2][iwv]);
		  speclibs[iZ].specs[ispec][iwv] = (speclibs[iZ].specs[ispec][iwv] < 1E-20) ?
		    -20 : log10(speclibs[iZ].specs[ispec][iwv]);
		}
	      ispec++;
	    }
	  printf("ispec=%d, nspec=%d\n",ispec, speclibs[iZ].nspec);
	}





      
      fits_close_file(fp_spec, &status);
      fp_spec=NULL;
     }
  NWV=speclibs[0].nwv; //NWV set here and only here
  //printf("Done loading spec lib files: %d\n",speclibs[0].nwv);
  return NWV;
}




void loadExts(char *file_ext)
{
  int status=0,hdutype,nullval,anynulls;
  fitsfile *fp_ext;
  //float *Exts=NULL;
  printf("file_ext name is %s\n", file_ext);
  Exts=malloc(sizeof(float)*NWV);
  if (fits_open_file(&fp_ext, file_ext, READONLY, &status)) printerror(status);
  if (fits_movabs_hdu(fp_ext, 2, &hdutype, &status)) printerror(status);
  fits_read_col(fp_ext, TFLOAT, 2, 1, 1, NWV, &nullval, Exts, &anynulls, &status);
  printf("Exts[500], Exts[509]: %f %f\n", Exts[500], Exts[509]);
  fits_close_file(fp_ext, &status);
  fp_ext=NULL;
  return;
}




void interpSingleStar(struct STAR star, float *outspec, float *outspecWave)
{
  int iZ,iZsub,iZup,imod;
  int it1, it2, ig1, ig2;
  float ftsub, ftup, fgsub, fgup;
  int iwv;
  float specsub[NWV], specup[NWV];
  float wht, totalwht, Zwht;
  //printf("%f, %f, %f, %f, %f, %f \n",star.Av,star.logg,star.logte, star.Mass, star.mu0, star.Z);
  //find the two Zs comprise the star
  iZsub=-99; iZup=-99;
  if (star.Z<speclibs[0].Z)
    {
      iZsub=0; iZup=0;
    }
  else if (star.Z>speclibs[NZ-1].Z)
    {
      iZsub=NZ-1;iZup=NZ-1;
    }
  else
    {
      for (iZ=1;iZ<NZ;iZ++)
	{
	  if (star.Z<speclibs[iZ].Z)
	    {
	      iZsub=iZ-1; iZup=iZ; break;
	    }
	}
    }

  //interp within each Z speclib
  for (iwv=0; iwv < NWV; iwv++) specsub[iwv]=0.;
  ftsub = locateitin1Dtab (star.logte, speclibs[iZsub].UniqueLogts,
			   speclibs[iZsub].NULogts, &it1, &it2, 0, 0, 0);
  fgsub = locateitin1Dtab (star.logg,  speclibs[iZsub].UniqueLoggs,
			   speclibs[iZsub].NULoggs, &ig1, &ig2, 0, 0, 0);
  //printf("sub: it1=%d, it2=%d\n", it1, it2);
  //printf("sub: ig1=%d, ig2=%d\n", ig1, ig2);
  totalwht=0.;
  for (imod=0; imod < speclibs[iZsub].nspec; imod++)
    {
      if ( (speclibs[iZsub].iLogts[imod] == it1 || speclibs[iZsub].iLogts[imod] == it2) &&
	   (speclibs[iZsub].iLoggs[imod] == ig1 || speclibs[iZsub].iLoggs[imod] == ig2) )
	{
	  wht = (speclibs[iZsub].iLogts[imod] == it1 ? ftsub : 1.-ftsub) *
	        (speclibs[iZsub].iLoggs[imod] == ig1 ? fgsub : 1.-fgsub);
	  totalwht += wht;
	  for (iwv=0; iwv < speclibs[iZsub].nwv; iwv++)
	    specsub[iwv] += wht * speclibs[iZsub].specs[imod][iwv];
	}
    }
  //printf("totalwht:%f\n",totalwht);
  if (totalwht <= 0.) {printf("zero total weight, quit!"); abort();}
  for (iwv=0; iwv < NWV; iwv++) specsub[iwv] /= totalwht;

  //printf("up-1: it1=%d, it2=%d\n", it1, it2);
  //printf("up-1: ig1=%d, ig2=%d\n", ig1, ig2);
  for (iwv=0; iwv < NWV; iwv++) specup[iwv]=0.;
  ftup = locateitin1Dtab (star.logte, speclibs[iZup].UniqueLogts,
			  speclibs[iZup].NULogts, &it1, &it2, 0, 0, 0);
  //printf("up0: it1=%d, it2=%d\n", it1, it2);
  fgup = locateitin1Dtab (star.logg,  speclibs[iZup].UniqueLoggs,
			  speclibs[iZup].NULoggs, &ig1, &ig2, 0, 0, 0);
  //printf("up1: it1=%d, it2=%d\n", it1, it2);
  //printf("up1: ig1=%d, ig2=%d\n", ig1, ig2);
  totalwht=0.;
  for (imod=0; imod < speclibs[iZup].nspec; imod++)
    {
      if ( (speclibs[iZup].iLogts[imod] == it1 || speclibs[iZup].iLogts[imod] == it2) &&
	   (speclibs[iZup].iLoggs[imod] == ig1 || speclibs[iZup].iLoggs[imod] == ig2) )
	{
	  wht = (speclibs[iZup].iLogts[imod] == it1 ? ftup : 1.-ftup) *
	        (speclibs[iZup].iLoggs[imod] == ig1 ? fgup : 1.-fgup);
	  totalwht += wht;
	  for (iwv=0; iwv < speclibs[iZup].nwv; iwv++)
	    specup[iwv] += wht * speclibs[iZup].specs[imod][iwv];
	}
    }
  if (totalwht <= 0.) abort();
  for (iwv=0; iwv < speclibs[iZup].nwv; iwv++) specup[iwv] /= totalwht;

  Zwht = iZsub == iZup ? 0. : (star.Z - speclibs[iZsub].Z)/(speclibs[iZup].Z - speclibs[iZsub].Z);
  for (iwv=0; iwv < speclibs[iZsub].nwv; iwv++) //assuming same wv grids
    outspec[iwv] = specsub[iwv] * Zwht + specup[iwv] * (1.-Zwht);

  //correct for distance and extinction
  for (iwv=0; iwv < speclibs[iZsub].nwv; iwv++)
    {
      outspec[iwv] = outspec[iwv] + lgG + log10(star.Mass) + lgMsun
	             - star.logg - 0.4*(star.mu0 + 5) - twolgpc2cm;
      outspec[iwv] = outspec[iwv] - lgExp * Exts[iwv] * star.Av;
    }

  for (iwv=0; iwv < speclibs[iZsub].nwv; iwv++)
    {
      outspecWave[iwv] = speclibs[iZsub].wvs[iwv];
    }
  return;
}