//-----------------------------------------------------------------------------------------
//  Function:   read_fov_file (string fovfile)
//
//  Reads the FOV data file.
//
//  Input:
//    file.....the complete path + name of the FOV data file
//
//  Notes:
//     The FOV file format is 5 header lines, then a matrix (# of bands+1, # elevation angles),
//     where each row is: elevation angle, then FOV response for bands 1-nband
//
//  Source: Mark Hervig, GATS
//-----------------------------------------------------------------------------------------

#include "sofiefov.h"
#include "ConvFunctions.h"

#include <fstream>
#include <iostream>
#include <cassert>
#include <algorithm>


void sofiefov::read_fov_file (const std::string& fovfile, int reverse) {
	fov_response.resize(16);
/* 
THis function is expecting to read a file with the following convention.
+ elevations values mean towards the ground    - elevation angles are toward space
*/

   char   header[512];   // for header lines in the file         

   std::ifstream inputFile;   //the input file

   inputFile.open(fovfile.c_str(), std::ios_base::in);  // open the file

   inputFile.getline(header,511);  // read a header line
   inputFile.getline(header,511);
   inputFile >> nel;               // read the number of elevation angles
   inputFile.getline(header,511);
   inputFile >> nband;             // read the number of bands
   inputFile.getline(header,511);
   inputFile.getline(header,511);

   std::cout << nel << " elevation angles and " << nband << " bands" << std::endl;

//   std::vector < double > temp_ar;
   double temp_sc;

   for (int i = 0; i < nel; i++) {  // loop over elevation

 //    temp_ar.clear();

//       if (j == 0) {  // the elevation
         inputFile >> temp_sc; 
         elevation.push_back(temp_sc);
   //      std::cout << temp_sc;
 //      }
     for (int j = 0; j < nband; j++) {  // loop over bands
 
//       if (j > 0) {  // the FOV response data
         inputFile >> temp_sc;
         fov_response[j].push_back(temp_sc);
    //     std::cout << " " << temp_sc;
 //      }

     }

 //    fov_response.push_back(temp_ar);
   //  std::cout << std::endl;
   }

   inputFile.close();
 
/* file sofie_fov_elevation_V1.3.txt has + angles toward space - toward ground
   we reverse that here so positive is toward ground and - is toward space */


if(reverse != 0) { 
   std::reverse(elevation.begin(),elevation.end() );
   std::transform(elevation.begin(), elevation.end(), 
		elevation.begin(),std::negate<double>() );
   for (int j = 0; j < nband; j++) {  // loop over bands
     std::reverse(fov_response[j].begin(), fov_response[j].end()) ;
   }
std::cout << " reversed FOV arrays " << std::endl;
}


 
/* we calculate the equal area weighted center of the band 3 FOV repsonse function and create a new
    elevation grid Eprime which is zero at the weighted center value.
*/
   	Eprime= elevation;
	std::transform(elevation.begin(), elevation.end(), Eprime.begin(), 
		std::bind2nd( std::plus<double>(), -midpoint(elevation,fov_response[2])  ));
		
std::cout << "new midpoint " << midpoint(Eprime, fov_response[2]) << std::endl;
//exit(23);
//   return;
}

std::vector<double>  sofiefov::get_fov_elev(const int band) 
{
	assert(band > 0 && band <= 16);
	return elevation;
}

std::vector<double>  sofiefov::get_fov_response(const int band) 
{
	assert(band > 0 && band <= 16);
	return fov_response.at(band-1);
}


 double sofiefov::midpoint( std::vector<double>& elev,  std::vector<double>& f)
{
	std::vector<double> S;
	S.push_back(0.0);

	for(unsigned j = 1; j < f.size(); ++j) 
		S.push_back( S.back() + 0.5*(f[j-1]+f[j])*(elev[j]-elev[j-1]) );

	double Half = 0.5* S.back(), Of;
	int npts = (int)S.size(), c__1 = 1;

	::linerd_(&S[0], &elev[0], &Half, &Of, &npts, &c__1);	
	return Of;
}



std::vector<double> sofiefov::get_fov_elev2(const int band, const double offset)
{
/* offset is under the sign convention that + elevation angles are towards space  and - elevation angles are
   towards the ground so this function will account for this
*/
	assert(band > 0 && band <= 16);

/* now we want each band's center weighted FOV to be at a particular angular location (offset)  */
	std::vector<double> E = Eprime;

//std::cout << "Band " << band << "  " << midpoint(Eprime,fov_response[band-1]) << std::endl;

// Notice the -offset  this accounts for the sign convention of offset

	std::transform(Eprime.begin(), Eprime.end(), E.begin(), 
		std::bind2nd( std::plus<double>(), -offset - midpoint(Eprime,fov_response[band-1])  ));


//std::cout << "band " << band << " new midpoint " << midpoint(E, fov_response[band-1]) << std::endl;
	return E;
}




//-----------------------------------------------------------------------------------------
//  Function:   get_fov_response(int band, double el_in) 
//
//  Returns the FOV relative response for the given band # and elevation angle.
//  The routine does 2-point interpolation of the FOV response vs. angle.
//  A response of 0 is returned if the input angle is outside the data range.
//
//  Input:
//    el_in....elevation angle, arcmin
//    band.....band number (1-nband)
//
//  Output:
//    relative FOV response (0-1)
//
//  Source: Mark Hervig, GATS
//-----------------------------------------------------------------------------------------
/*
double fov::get_fov_response(int band, double el_in) {

  double resp=0., term;
  	
// check the input for problems

  if (band < 1 || band > nband) {
    cout << "band number out of range in fov::get_fov_response:  " << band << endl;
    cout << "there are " << nband << " bands available" << endl;
    return resp;
  }

// interpolate FOV in angle

  for (int i = 0; i < nel; i++) {   // loop over nel angles

   if ( ((elevation[i] <= el_in) && (elevation[i+1] >= el_in)) || ((elevation[i] >= el_in) && (elevation[i+1] <= el_in)) ) {
	  
    term  = (el_in - elevation[i]) / (elevation[i+1] - elevation[i]);    
    resp  = term * (fov_response[i+1][band-1] - fov_response[i][band-1]) + fov_response[i][band-1];    
    return resp;
   }

  }    
  return resp;  // return a response of 0 if outside the data range
}
*/