//#include "F77Functions.h"

#include "LoadLinesFromDB.h"
#include "GATS_Utilities.hpp"
#include "F77Radtran.h"
#include "Molecule.h"


#include <numeric>
#include <cassert>
#include <iostream>

using namespace GATS_Utilities;
using namespace gatsDBpp;

int LoadLinesFromDB( const std::vector<int>& igas, 
	const double nu1, const double nu2, const float wnmax, 
	const int mcmix , const int mlmix, std::vector<float>& afgl,
	GATS_DB* dbConnection, const std::string& tableName   )
	
{
	static char r[] = "2003";
	int ifact, igs, iso, nl, ngas, nlines;
	std::string sqlstring;
	float selfHW, foreignHW, *estart, *estop;
	std::vector<int> linecount;
	GATS_DB_Results R;
	
	nl = 0;
	ngas = igas.size();
//count up the total number of lines that will be queried from DB
	for(int i = 0; i < ngas ; ++i ) {
    		assert(igas[i] >= 100 );
       	ifact = (igas[i] > 10000) ? 10000 : 1;
     	igs = igas[i] / (ifact*100);
    		iso=igas[i]/ifact - igs * 100; 	
    		                
		sqlstring = std::string("SELECT COUNT(*) FROM ") + tableName ;
		sqlstring += "  WHERE frequency BETWEEN ";
		sqlstring += ConvertToStringPrec( nu1-wnmax );
		sqlstring += " AND ";
		sqlstring += ConvertToStringPrec( nu2+wnmax );
		sqlstring += " AND molecule=";
		sqlstring += ConvertToString( igs );
		if(iso > 0 ) {
			sqlstring += " AND isotope=";
			sqlstring += ConvertToString( iso );
		}
//		sqlstring += " AND intensity > 1.e-32 ";  // need to convert intensity to double	

// call the database to get the line count
		R = dbConnection->Query(sqlstring.c_str() );
		linecount.push_back( R[0][0].asInt() );  
		std::cout << "Found " << linecount.back() << " lines for Gas ID " << igas[i]
		<< std::endl;
	}
// total number of lines for all gases	
	nl = std::accumulate(linecount.begin(),linecount.end(), 0);
//	std::cout << " ready for memory allocation " << std::endl;
// Allocate the AFGL array
//
//      NAFGL=2+13*NL+MCMIX*(1+MLMIX)+NGAS*(13+MQSET+MQSET*MISO*
 //    &(1+2*MCOEF))
//std::cout << "need " <<  nl * 13 + 2 + mcmix * (mlmix + 1) + ngas * 369  << std::endl;
	afgl.resize( nl * 13 + 2 + mcmix * (mlmix + 1) + ngas * 369  );
//std::cout << "done with allocation " << std::endl;
// Set the common block pointers
    ::lp_loca__.lv0 = 1;
// making sl a real*8
    ::lp_loca__.lsl = ::lp_loca__.lv0 + (nl << 1) ;
//
    ::lp_loca__.lqcoef = ::lp_loca__.lsl +  (nl << 1);
    ::lp_loca__.liranacc = ::lp_loca__.lqcoef + ngas * 320;
    ::lp_loca__.lils = ::lp_loca__.liranacc + 1;
    ::lp_loca__.ljgas = ::lp_loca__.lils + ngas;
    ::lp_loca__.lnlines = ::lp_loca__.ljgas + ngas;
    ::lp_loca__.lilines = ::lp_loca__.lnlines + ngas;
    ::lp_loca__.lindex = ::lp_loca__.lilines + ngas;
    ::lp_loca__.liset = ::lp_loca__.lindex + ngas;
    ::lp_loca__.lnmode = ::lp_loca__.liset + ngas;
    ::lp_loca__.lwtt = ::lp_loca__.lnmode + ngas;
    ::lp_loca__.lgamma = ::lp_loca__.lwtt + ngas;
    ::lp_loca__.lalmin = ::lp_loca__.lgamma + ngas;
    ::lp_loca__.lalmax = ::lp_loca__.lalmin + ngas;
    ::lp_loca__.lalptem = ::lp_loca__.lalmax + ngas;
    ::lp_loca__.lrelease = ::lp_loca__.lalptem + ngas;
    ::lp_loca__.lfound = ::lp_loca__.lrelease + 1;
    ::lp_loca__.lifit = ::lp_loca__.lfound + ngas;
    ::lp_loca__.lnq = ::lp_loca__.lifit +  (ngas << 2);
    ::lp_loca__.llset = ::lp_loca__.lnq + (ngas << 5);
    ::lp_loca__.ljol = ::lp_loca__.llset + nl;
//    ::lp_loca__.lsl = ::lp_loca__.ljol + nl;
    ::lp_loca__.laln = ::lp_loca__.ljol + nl;
    ::lp_loca__.lals = ::lp_loca__.laln + nl;
    ::lp_loca__.lel = ::lp_loca__.lals + nl;
    ::lp_loca__.lwt = ::lp_loca__.lel + nl;
    ::lp_loca__.lalgam = ::lp_loca__.lwt + nl;
    ::lp_loca__.lpsc = ::lp_loca__.lalgam + nl;
    ::lp_loca__.lipmix = ::lp_loca__.lpsc + nl;
    ::lp_loca__.lcfmix = ::lp_loca__.lipmix + nl;
    ::lp_loca__.lpadd = ::lp_loca__.lcfmix + 2;

//std::cout << "setting commonblock " << std::endl;
// set afgl to use the 2003 TIPS release. 
	for(int i=0; i< 4; i++)
    		*(i+ (char*)(&afgl[0]+::lp_loca__.lrelease-1) ) = r[i];

// set the use random access file flag to false or do not use
	*( (int*)(&afgl[0]+::lp_loca__.liranacc-1) ) = -1;

// now read in the lines

	nlines = 0;
	for(int i = 0; i < ngas ; ++i ) {
		*(i + (int*)(&afgl[0]+::lp_loca__.lnlines-1) ) = linecount[i];
		*(i + (int*)(&afgl[0]+::lp_loca__.lilines-1) ) = linecount[i];
		*(i + (int*)(&afgl[0]+::lp_loca__.ljgas-1) ) = nlines+1;
		
		if (linecount[i] > 0) {
    	    		ifact = (igas[i] > 10000) ? 10000 : 1;
     			igs = igas[i] / (ifact*100);
    			iso=igas[i]/ifact - igs * 100;
    			Molecule M(igs, iso);
    			*(i+ &afgl[0]+::lp_loca__.lwtt-1) = static_cast<float>(M.Weight(1));  //Wt of primary isotope
  		                
			sqlstring = std::string("SELECT ");
//                                       0           1         2                 3
			sqlstring += "frequency, intensity, air_halfwidth, self_halfwidth, ";
//                                       4              5         6       7
			sqlstring += "lower_energy, t_exponent, p_shift, isotope ";

			sqlstring += " FROM " + tableName ;
			sqlstring += "  WHERE frequency BETWEEN ";
			sqlstring += ConvertToStringPrec( nu1-wnmax );
			sqlstring += " AND ";
			sqlstring += ConvertToStringPrec( nu2+wnmax );
			sqlstring += " AND molecule=";
			sqlstring += ConvertToString( igs );
			if(iso > 0 ) {
				sqlstring += " AND isotope=";
				sqlstring += ConvertToString( iso );
			}
//			sqlstring += " AND intensity > 1.e-32 ";  // need to convert intensity to double	
			sqlstring += " Order by frequency";

std::cout << sqlstring << std::endl;
			R = dbConnection->Query(sqlstring.c_str() );
//std::cout << "called query " << std::endl;
			GATS_DB_Results::const_iterator resIter = R.begin();

//std::cout << "finished query have " << R.size() << "  Lines " << std::endl;

			for ( ; resIter != R.end(); ++resIter) {

				foreignHW = (*resIter)[2].asFloat();
				selfHW = (*resIter)[3].asFloat() ;
// in the vain of LP_RAFGL it resets a <= 0 halfwidth to the other halfwidth
				if(selfHW <= 0.0) { selfHW = foreignHW; }
				if(foreignHW <= 0.0) { foreignHW = selfHW ; }
				
				*(nlines + (int*)(&afgl[0]+::lp_loca__.ljol-1) ) = (*resIter)[7].asInt()  ;		//isotope
				*(nlines+ (double*)(&afgl[0]+::lp_loca__.lv0-1) ) =  (*resIter)[0].asDouble() ;	//frequency
				*(nlines+ (double*)(&afgl[0]+::lp_loca__.lsl-1) )= (*resIter)[1].asDouble() ; 	//intensity
				*(nlines+ &afgl[0]+::lp_loca__.laln-1) = foreignHW ;     //Air Broadened HW
				*(nlines+ &afgl[0]+::lp_loca__.lals-1) = selfHW ; 			//Self Broadened HW
				*(nlines + &afgl[0]+::lp_loca__.lel-1) =  (*resIter)[4].asFloat() ;	//Lower State Energy
				*(nlines + &afgl[0]+::lp_loca__.lalgam-1) =  (*resIter)[5].asFloat() ;  // temperature dependence of air hw
				*(nlines + &afgl[0]+::lp_loca__.lpsc-1) =(*resIter)[6].asDouble() ;	   // pressure shift
				*(nlines + &afgl[0]+::lp_loca__.lwt-1) = M.Weight( (*resIter)[7].asInt() ) ;  // molecular weight
				*(nlines + (int*)(&afgl[0]+::lp_loca__.lipmix-1) ) = 0;   // No Mixing allowed
	 			nlines++;		
			}
			estart =   &afgl[0]+ ::lp_loca__.laln-1 + nlines - linecount[i] ;
			estop =    &afgl[0]+ ::lp_loca__.laln-1 + nlines   ;
			*(i+ &afgl[0]+::lp_loca__.lalmax-1) = *std::max_element ( estart, estop );  //ALMIN
			*(i+ &afgl[0]+::lp_loca__.lalmin-1) = *std::min_element ( estart, estop );  //ALMax

			estart =   &afgl[0]+ ::lp_loca__.lalgam-1 + nlines - linecount[i] ;
			estop  =   &afgl[0]+ ::lp_loca__.lalgam-1 + nlines  ;
			std::vector<float> xx(estart, estop );
			std::sort(xx.begin(), xx.end() );
			*(i+ &afgl[0]+::lp_loca__.lalptem-1)  =  xx[ xx.size()/2 ];    // ALPTEM parameter
		} else {   // Line count is 0
			*(i+ &afgl[0]+::lp_loca__.lalptem-1) = 0. ;
			*(i+ &afgl[0]+::lp_loca__.lalmax-1) = 0. ;
			*(i+ &afgl[0]+::lp_loca__.lalmin-1) =  0. ;
		}
		
	}
// replace non physical (<= 0.) lower state energies with 300
       std::replace_if(&afgl[0]+ ::lp_loca__.lel-1, &afgl[0]+ ::lp_loca__.lel-1 + nl,
              std::bind2nd(std::less_equal<float>(),0.0),300.0);

//std::cout << "all lines loaded " << std::endl;
	return nl;
}