#include <XSUtil/Utils/XSutility.h>
#include <XSUtil/FunctionUtils/FunctionUtility.h>
#include <XSUtil/FunctionUtils/xsFortran.h>
#include <fitsio.h>
#include <fstream>

// Translated from ldaped.f 12/08

typedef XSutility::auto_array_ptr<float> FloatAP;
typedef XSutility::auto_array_ptr<int> IntAP;

extern "C" void rdaped_(int& iocont, int& ioline, int& nmtval, int& nmelt, int& ncont,
                   int& npseudo, int& nlines, float* tval, float* cont, float* cenerg,
                   int* conEl, int* conpt, int* consiz, float* pseudo, float* penerg, 
                   int* psupt, int* psusiz, float* lineE, float* lineEm, int* lineEl,
                   int* linpt, int* linsiz, int* status); 

void ldaped(int itype, FloatAP& tval, int& nmtval, int& nmelt, FloatAP& cont,
        FloatAP& cenerg, IntAP& conEl, IntAP& conpt, IntAP& consiz, 
        FloatAP& pseudo, FloatAP& penerg, IntAP& psupt, IntAP& psusiz,
        FloatAP& lineE, FloatAP& lineEm, IntAP& lineEl, IntAP& linpt,
        IntAP& linsiz, int& status)
{
   // Subroutine to read the APED style plasma model files.
   // Arguments :
   //    itype      I        i: type of plasma emission file
   //    tval      ptr       r: pointer to array of tabulated temperatures
   //    nmtval     I        r: number of tabulated temperatures
   //    nmelt      I        r: number of tabulated elements
   //    cont      ptr       r: Tabulated continua
   //    cenerg    ptr       r: Energies for tabulated continua
   //    conl      ptr       r: Atomic numbers for tabulated elements
   //    conpt     ptr       r: Offsets in the cont array for each Z and kT
   //    consiz    ptr       r: Sizes in the cont array for each Z and kT
   //    pseudo    ptr       r: Tabulated pseudo-continua
   //    penerg    ptr       r: Energies for tabulated pseudo-continua
   //    psupt     ptr       r: Offsets in the pseudo array for each Z and kT
   //    psusiz    ptr       r: Sizes in the pseudo array for each Z and kT
   //    lineE     ptr       r: Tabulated line energies
   //    lineEm    ptr       r: Tabulated line emissivities
   //    lineEl    ptr       r: Tabulated line elements
   //    linpt     ptr       r: Offsets in the line arrays for each kT
   //    linsiz    ptr       r: Sizes in the line arrays for each kT

   // Set the names of the input files and open them. Note that for APEC
   // we check for the APECROOT string model parameter to allow the user
   // to substitute different _coco and _line tables.

   const string& datadir = FunctionUtility::modelDataPath();
   static string ovalue;
   string pname("APECROOT");
   string pvalue(FunctionUtility::getModelString(pname));

   string contfil;
   string linefil;
   if (itype == 1)
   {
      contfil = datadir + "RS93_newx.fits";
      linefil = datadir + "RS93_line.fits";
   }
   else if (!pvalue.length() || pvalue == FunctionUtility::NOT_A_KEY())
   {
      contfil = datadir + "apec_v2.0.2_coco.fits";
      linefil = datadir + "apec_v2.0.2_line.fits";
   }
   else
   {
      // if APECROOT has been set first test whether only a version number has
      // changed. Then assume that only a root filename has been given. If neither
      // these work then assume that an entire directory path was included.

      contfil = datadir + "apec_v" + pvalue + "_coco.fits";
      linefil = datadir + "apec_v" + pvalue + "_line.fits";
      std::ifstream test1(contfil.c_str());
      if (!test1)
      {
         contfil = datadir + pvalue + "_coco.fits";
         linefil = datadir + pvalue + "_line.fits";
         std::ifstream test2(contfil.c_str());
         if (!test2)
         {
            contfil = pvalue + "_coco.fits";
            linefil = pvalue + "_line.fits";            
         }
      }
      if (ovalue != pvalue)
      {
         ovalue = pvalue;
         string msg("Reading APEC data from " + pvalue);
         xs_write(const_cast<char*>(msg.c_str()), 10);
      }
   }

   fitsfile *pContFil=0;
   fits_open_file(&pContFil, const_cast<char*>(contfil.c_str()), 0, &status);
   if (status)
   {
      std::ostringstream err;
      err << "Failed to open " << contfil <<"\nError in LDAPED : status = "
          << status;
      xs_write(const_cast<char*>(err.str().c_str()), 10);
      return;
   }
   else
   {
      string msg("Loading contents of " + contfil);
      xs_write(const_cast<char*>(msg.c_str()),15);
   }

   fitsfile *pLineFil=0;
   fits_open_file(&pLineFil, const_cast<char*>(linefil.c_str()), 0, &status);
   if (status)
   {
      std::ostringstream err;
      err << "Failed to open " << linefil <<"\nError in LDAPED : status = "
          << status;
      xs_write(const_cast<char*>(err.str().c_str()), 10);
      int status2=0;
      fits_close_file(pContFil, &status2);
      return;
   } 
   else
   {
      string msg("Loading contents of " + linefil);
      xs_write(const_cast<char*>(msg.c_str()),15);
   }

   // First scan the files to get the number of temperatures and elements
   // and to find the sizes of arrays that will be required. First the
   // continuum data

   int ncont=0;
   int npseudo=0;
   nmtval = 0;
   int hdutyp=0;
   std::ostringstream context;

   fits_movrel_hdu(pContFil, 1, &hdutyp, &status);
   bool isOk = true;
   while (status == 0)
   {
      char extnam[FLEN_VALUE+1];
      if (fits_read_key_str(pContFil, "EXTNAME", extnam, 0, &status))
      {
         context << "Failed to read cont file EXTNAME for temperature # " << nmtval;
         isOk = false;
         break;
      }
      if (string(extnam) == "EMISSIVITY")
      {
         ++nmtval;
         long tmpNmelt=0;
         if (fits_read_key_lng(pContFil, "NAXIS2", &tmpNmelt, 0, &status))
         {
            context << "Failed to read cont file NAXIS2 for temperature # " << nmtval;
            isOk = false;
            break;
         }
         nmelt = static_cast<int>(tmpNmelt);

         // Sum up the continuum elements required

         int icol=0;
         if (fits_get_colnum(pContFil, CASEINSEN, (char *)"N_Cont", &icol, &status))
         {
            context << "Failed to find N_Cont column for temperature # " 
                    << nmtval;
            isOk = false;
            break;
         }

         for (int i=1; i<=nmelt; ++i)
         {
            long tmpBuf=0;
            int anynul=0;
            if (fits_read_col_lng(pContFil, icol, i, 1, 1, 0, &tmpBuf, 
                        &anynul, &status))
            {
               context << "Failed to read N_cont : " << nmtval << " " << i;
               isOk = false;
               break;
            }
            ncont += static_cast<int>(tmpBuf);
         }
         if (!isOk) break;

         // Sum up the pseudo-continuum elements required

         if (fits_get_colnum(pContFil, CASEINSEN, (char *)"N_Pseudo", &icol, &status))
         {
            context << "Failed to find N_Pseudo column for temperature # " << nmtval;
            isOk = false;
            break;
         }

         for (int i=1; i<=nmelt; ++i)
         {
            long tmpBuf=0;
            int anynul=0;
            if (fits_read_col_lng(pContFil, icol, i, 1, 1, 0, &tmpBuf,
                        &anynul, &status))
            {
               context << "Failed to read N_Pseudo : " << nmtval << " " << i;
               isOk = false;
               break;
            }
            npseudo += static_cast<int>(tmpBuf);
         }
         if (!isOk) break;
      } // end this EMISSIVITY extension, move on to next

      fits_movrel_hdu(pContFil, 1, &hdutyp, &status);

   } // end contfil extension loop

   if (isOk)
   {
      status = 0;

      // Now the line data - this is easier since we only need the NAXIS2 for
      // each extension

      int nlines = 0;
      int iCheck = 0;
      fits_movrel_hdu(pLineFil, 1, &hdutyp, &status);
      while (status == 0)
      {
         char extnam[FLEN_VALUE+1];
         if (fits_read_key_str(pLineFil, "EXTNAME", extnam, 0, &status))
         {
            context << "Failed to read line file EXTNAME for temperature # " << iCheck;
            isOk = false;
            break;
         }
         if (string(extnam) == "EMISSIVITY")
         {
            ++iCheck;
            long tmpNlines=0;
            if (fits_read_key_lng(pLineFil, "NAXIS2", &tmpNlines, 0, &status))
            {
               context << "Failed to read line file NAXIS2 for temperature # " << iCheck;
               isOk = false;
               break;
            }
            nlines += static_cast<int>(tmpNlines);            
         }
         fits_movrel_hdu(pLineFil, 1, &hdutyp, &status);
      }
      if (isOk && iCheck != nmtval)
      {
         context << "EMISSIVITY extensions mismatch between line and continuum files.";
         isOk = false;
      }
      if (isOk)
      {
         status = 0;

         // Reset both files to the primary HDU

         if (fits_movabs_hdu(pContFil, 1, &hdutyp, &status))
         {
            context << "Failed to move back to continuum file primary header.";
            isOk = false;
         }
         else if (fits_movabs_hdu(pLineFil, 1, &hdutyp, &status))
         {
            context << "Failed to move back to line file primary header.";
            isOk = false;
         }
         if (isOk)
         {
            // Grab the memory that will be needed

            float *oldF = tval.reset(new float[nmtval]);
            delete [] oldF;
            oldF = cont.reset(new float[ncont]);
            delete [] oldF;
            oldF = cenerg.reset(new float[ncont]);
            delete [] oldF;
            int * oldI = conEl.reset(new int[nmelt]);
            delete [] oldI;
            oldI = conpt.reset(new int[nmtval*nmelt]);
            delete [] oldI;
            oldI = consiz.reset(new int[nmtval*nmelt]);
            delete [] oldI;
            oldF = pseudo.reset(new float[npseudo]);
            delete [] oldF;
            oldF = penerg.reset(new float[npseudo]);
            delete [] oldF;
            oldI = psupt.reset(new int[nmtval*nmelt]);
            delete [] oldI;
            oldI = psusiz.reset(new int[nmtval*nmelt]);
            delete [] oldI;
            oldF = lineE.reset(new float[nlines]);
            delete [] oldF;
            oldF = lineEm.reset(new float[nlines]);
            delete [] oldF;
            oldI = lineEl.reset(new int[nlines]);
            delete [] oldI;
            oldI = linpt.reset(new int[nmtval]);
            delete [] oldI;
            oldI = linsiz.reset(new int[nmtval]);
            delete [] oldI;

            // Read the files and load the arrays
            int iocont = CFITS2Unit(pContFil);
            int ioline = CFITS2Unit(pLineFil);
            rdaped_(iocont, ioline, nmtval, nmelt, ncont, npseudo, nlines,
               tval.get(), cont.get(), cenerg.get(), conEl.get(), conpt.get(),
               consiz.get(), pseudo.get(), penerg.get(), psupt.get(), psusiz.get(),
               lineE.get(), lineEm.get(), lineEl.get(), linpt.get(), linsiz.get(),
               &status);
            if (status)
            {
               context << "Failure in RDAPED";
               isOk = false;
            }
         }
      }
   }

   // Close the files. This is after the break-out point from errors to ensure 
   // that everything gets shut down correctly.

   fits_close_file(pContFil, &status);
   fits_close_file(pLineFil, &status);

   if (status || !isOk)
   {
      xs_write(const_cast<char*>(context.str().c_str()), 10);
      std::ostringstream ossErr;
      ossErr << "Error in LDAPED : status = " << status;
      xs_write(const_cast<char*>(ossErr.str().c_str()), 10);      
   }
}