//
//  $Id$
//-----------------------------------------------------------------------
//
//            (c) Copyright 2009 by GATS, Inc.,
//     11864 Canon Blvd, Suite 101, Newport News VA 23606
//
//  All Rights Reserved. No part of this software or publication may be
//  reproduced, stored in a retrieval system, or transmitted, in any form
//  or by any means, electronic, mechanical, photocopying, recording, or
//  otherwise without the prior written permission of GATS, Inc.
//
//-----------------------------------------------------------------------
//
// Module:	FOVcorrection.cpp
//
// Author:	Lance Deaver
//
// Date:  	Thu Sep 17 17:24:50 EDT 2009
//
//-----------------------------------------------------------------------
//
// Modification History:
//
//     $Log$
//
//-----------------------------------------------------------------------
//
//-----------------------------------------------------------------------
// Include Files:
//-----------------------------------------------------------------------
//
#include "FOVcorrection.h"
#include "ConfigFile.h"
#include "Event.h"
#include "sofiefov.h"

#include "F77Radtran.h"
#include "F77Functions.h"

#include <iostream>
//
//-----------------------------------------------------------------------
// Defines and Macros:
//-----------------------------------------------------------------------
//



//
//-----------------------------------------------------------------------
// Global Variables:
//-----------------------------------------------------------------------
//



//
//-----------------------------------------------------------------------
// Utility Routines:
//-----------------------------------------------------------------------
//

int FOVcorrection(Event& L0, Event& L1, Event& Tmp, Event& SD, ConfigFile& cf)
{

  static sofiefov gFOVdata;
  EventVar<double> lockdownEl, angleSLDC, Za, ApparentAng, Z, P, T, RadCurv, TP_Lat;
  EventVarVect<double> LDcurves, MeasuredTrans;
// bandpass wavenumbers (approximately)
//                           1      2      3        4     5      6      7       8
    static float wave[] = { 33800.,30100., 11600., 9600., 4060., 3810., 3590., 3395.,
//           9        10         11         12        13        14      15      16
          3260.,     3140.,     2950.,    2875.,     2310.,    2150.,  2000.,  1885. };
    static int mcel;
    static int ncel;
    static int mpout;

    static int c_one = 1;
    static int debug = 0;           // 0 is a false (this is the debug flag)
    static int irf =  1;            // 1 is true always do refraction
    static int iplat = 0;       // observer in space
    static int ic1 = -1;
    static int ic2 = 1;
    static  bool firstCall = true;
    

    static std::vector<float> za, zt, pt, tt, gr, pout;
    static std::vector<double> zdub,apparentAng, trueAng,lockang;
//    static std::string  AltEventVarName,
//        RefracEventVarName, sectionName ;

//  float lat83;
    float  re, tearth,xlat;
    int mr, nr;
    std::vector<double> refractionAng;
    std::vector< std::vector<double> > DetectorData;
    
    std::vector<float> zr,pr,tr,tgr;

    if(firstCall) {
// read in the FOV file
// std::string fovFileName = cf.GetStr("AltitudeRegistration", "FOVFile");

      gFOVdata.read_fov_file(cf.GetStr("FOVcorrection", "FOVFile"), 1 );  // the 1 is a reverse flag for 1.3


      float z = 180.;
       while(z > 2.8) {
     za.push_back(z);
         z -= 1.0;
       }
/*
       float z = 150.;
       while(z > 9.) {
	 za.push_back(z);
         z -= 1.0;
       }
 */       
        ncel=mcel=(int)za.size() ;
            int it2 = mcel * mcel;
            int it1 = it2 + mcel;
            mpout = mcel * 11 + 3 + it1 * 12 + it2 * 3;
        pout.resize(mpout);
        zt.resize(mcel);
        pt.resize(mcel);
        tt.resize(mcel);
	apparentAng.resize(mcel);
	trueAng.resize(mcel);
	lockang.resize(mcel);

        gr.resize(2*mcel);
// Call  KP_inita  so we can PROCEED

        ::kp_inita__(&debug, &mcel, (int*)&pout[0]  )  ;
        zdub = std::vector<double>(za.begin(), za.end() );
	
        firstCall = false;
     }

  std::cout << "\nIn FOVcorrection \n";

  L1.getEventVar("L1_SolarModel", LDcurves);
  L1.getEventVar("L1_SLDCAngleFromTopEdge", angleSLDC);

  L1.getEventVar("gridImpact", Za);
  std::vector<double> zimpact(&Za[0], &Za[0]+Za.size() );

  L1.getEventVar("L1_EffectiveElLockdown",lockdownEl);
  std::vector<double> fovLock(&lockdownEl[0], &lockdownEl[0]+lockdownEl.size() );
//  double lockdownExo = lockdownEl[0];

  L1.getEventVar("gridAngleImpact", ApparentAng);
  std::vector<double> appang(&ApparentAng[0], &ApparentAng[0]+ApparentAng.size() );

  L1.getEventVar("GRIDDED_TRANSMISSIONS", MeasuredTrans);
  
  for(int i=0; i<16;++i) {
      EventVar<double> Ev = MeasuredTrans[i];
      DetectorData.push_back( std::vector<double>(&Ev[0], &Ev[0]+Ev.size())  );
  }
//  filter out missing data

  std::vector<double>::iterator res = std::find_if(zimpact.begin(), zimpact.end(), 
	std::bind2nd(std::greater<double>(),0.) );
  unsigned ik = std::distance(zimpact.begin(), res);
  if(ik != 0)  {
        zimpact.erase( zimpact.begin(), zimpact.begin()+ik);
	fovLock.erase(fovLock.begin(), fovLock.begin()+ik);
	appang.erase(appang.begin(), appang.begin()+ik) ;
        for(int i=0; i<16;++i) {
           DetectorData[i].erase( DetectorData[i].begin(), DetectorData[i].begin()+ik );
        } 
  } 

  res = std::find_if(zimpact.begin(), zimpact.end(), 
	std::bind2nd(std::less<double>(),0.) );
  ik = std::distance(zimpact.begin(), res);
  if(ik != zimpact.size() )  {
        zimpact.erase( zimpact.begin()+ik, zimpact.end() );
	fovLock.erase(fovLock.begin()+ik, fovLock.end());
	appang.erase(appang.begin()+ik, appang.end() ) ;
        for(int i=0; i<16;++i) {
           DetectorData[i].erase( DetectorData[i].begin()+ik, DetectorData[i].end() );
        } 
  } 

std::reverse(zimpact.begin(), zimpact.end() );
std::reverse(fovLock.begin(), fovLock.end());
std::reverse(appang.begin(), appang.end() ) ;
for(int i=0; i<16;++i) {
    std::reverse( DetectorData[i].begin(), DetectorData[i].end() );
}
 int numObs=(int)zimpact.size();


  ::linerd_(&zimpact[0],&appang[0],&zdub[0], &apparentAng[0], &numObs, &ncel  );
  ::linerd_(&zimpact[0],&fovLock[0],&zdub[0], &lockang[0], &numObs, &ncel  );
  
  
  L1.getEventVar("L1_CurvatureRadius", RadCurv);
  L1.getEventVar("L1_TPLat", TP_Lat);
  xlat = TP_Lat[0];
  re = RadCurv[0];

// Ok so we have the solar limb darkening curves and measured transmissions. . we need a zpt profile
// to calculate the atmospheric refraction
  L1.getEventVar("Z_Merged", Z);
  L1.getEventVar("P_Merged", P);
  L1.getEventVar("T_Merged", T);
// Now we will Loop over all 16 bands and convolve the transmission over the LDC with refraction using the
// FOV function;

// make floats from doubles
    zr = std::vector<float>(&Z[0], &Z[0]+Z.size() );
    pr = std::vector<float>(&P[0], &P[0]+P.size() );
    tr = std::vector<float>(&T[0], &T[0]+T.size() );
    nr = mr=(int)zr.size();
// no temperature Gradients
    tgr = std::vector<float>(2*mr, 0.);


for(int iband =0; iband < 16; ++iband) {

     EventVar<double> respSLDC = LDcurves[iband];
     int numSLDCs = (int)respSLDC.size();
   

    ::kp_mpath__( &mr,  &mcel,  &mpout,  &irf,  &c_one, &ncel, &iband, &re, &ncel, &nr,
        &wave[iband], &iplat, &za[0], &zr[0], &pr[0], &tr[0], &tgr[0],
        &zt[0], &pt[0], &tt[0], &gr[0], &tearth, &pout[0] );

      refractionAng.clear();
      for(int ir=1; ir<= ncel; ++ir ) {
         refractionAng.push_back( ::refraction_( &ir, &ic1, &ic2, &pout[0] ) );
	 trueAng[ir-1] = apparentAng[ir-1]+refractionAng[ir-1];
      }

      std::vector<double> angleFOV = gFOVdata.get_fov_elev(iband+1);
//for(int i = 0; i< angleFOV.size(); ++i)
// std::cout << angleFOV[i] << std::endl;
//exit(23);
      std::vector<double> respFOV = gFOVdata.get_fov_response(iband+1);
      std::transform( angleFOV.begin(), angleFOV.end(), angleFOV.begin(), 
          std::bind2nd(std::multiplies<double>(), M_PI/10800.0) );  //convert to radians
      int numFOVs = (int)angleFOV.size();
 double lockdownExo = fovLock[0];      
// determine the exo lockdown value
      double simExo;      
//std::cout << "insidde fovcorrections " << std::endl;
/*
for(int i = 0 ; i< numSLDCs; i++) {

std::cout << angleSLDC[i]*10800./M_PI << "  " <<  respSLDC[i] << std::endl;

}
std::cout << " now fov " << std::endl;
for(int i = 0 ; i< numFOVs; i++) {

std::cout << angleFOV[i] << "  " <<  respFOV[i] << std::endl;

}
*/
  ::calc_sim_exo__(&respFOV[0], &angleFOV[0], &numFOVs, &respSLDC[0], 
           &angleSLDC[0], &numSLDCs, &lockdownExo, &simExo);
	   
       std::vector<double> tau(mcel);
std::cout <<" lock " <<  lockdownExo*180.*60./M_PI << "   "  << simExo << std::endl;
//exit(23);

::linerd_(&zimpact[0],&(DetectorData[iband][0]),&zdub[0], &tau[0], &numObs, &ncel  ); 
     //   = DetectorData[iband]; // = MeasuredTrans[iband];
//std::cout << lockdownExo << std::endl << std::endl;
/*
if(iband==0) {
	for(int i=0; i< numObs; ++i) 
	std::cout << DetectorData[iband][i] << "  " <<  zimpact[i] << std::endl;       

std::cout << "&^^^^^" << std::endl; 
for(int i =0; i< ncel; i++) 
   std::cout << tau[i] << "  " << zdub[i] << std::endl;
exit(23);
}
*/
      std::vector<double> respConvolved(mcel);
     
	   
      ::conv_fov_atm_sun__(
                 &ncel, &c_one, &ncel, 
                 &apparentAng[0],  &trueAng[0], &lockang[0], &tau[0],
                 &respFOV[0], &angleFOV[0], &numFOVs, 
		         &respSLDC[0], &angleSLDC[0], &numSLDCs, 
		         &simExo, &respConvolved[0]); 
/*
if(iband==0) {
for(int i =0; i< ncel; i++)
   std::cout << tau[i] << "  " << respConvolved[i] << "  " << zdub[i] << std::endl;

exit(23);
}
*/

std::vector<double> correction(numObs);
::linerd_(&zdub[0], &respConvolved[0], &zimpact[0], &correction[0], &ncel, &numObs);
/*
if(iband==0) {
for(int i =0; i< ncel; i++)
   std::cout <<  respConvolved[i] << "  " << zdub[i] << std::endl;
std::cout << "&" << std::endl;
for(int i =0; i< numObs; i++)
   std::cout <<  correction[i] << "  " << zimpact[i] << std::endl;


exit(23);
}
*/

for(int i=0; i< numObs; ++i) 
 correction[i] = 2.0*DetectorData[iband][i] - correction[i];


 
 
EventVar<double> Ev = MeasuredTrans[iband];

//std::cout << "# Before Band " << iband << std::endl;
//for(int i=Za.size()-1; i>=0; i--) 
//  std::cout << Ev[i] << "  " << Za[i] << std::endl;
//std::cout << "&" << std::endl;

int j=0;
for(int i=Za.size()-1; i>=0; i--) {
if(j < numObs && zimpact[j] == Za[i]  ) {
//	std::cout << Za[i] << std::endl;
	Ev[i] = correction[j];
	j++;
}
}
//std::cout << "#after Band "<< std::endl;
//for(int i=Za.size()-1; i>=0; i--)
//  std::cout << Ev[i] << "  " << Za[i] << std::endl;
//std::cout << "&" << std::endl;



MeasuredTrans[iband] = Ev;


//exit(23);
//std::cout << "# band " << iband+1 << std::endl;
//for(int i=0; i< ncel; ++i) 
//	std::cout << respConvolved[i] << "  " <<  zt[i] << std::endl;
//std::cout << "&" << std::endl;
//exit(23);


/*
    for(int i=0; i< Za.size(); ++i) {
//  if(Za[i] > 0.) 
  std::cout << Ev[i] << "  " << Za[i] <<   std::endl;
  }
std::cout << "&" << std::endl;
    for(int i=0; i< numObs; ++i) {
//  if(Za[i] > 0.)
  std::cout << DetectorData[iband][i] << "  " << zimpact[i] <<   std::endl;
  }
std::cout << "&" << std::endl;

//exit(23);
}
exit(23);
*/
}

L1.replaceEventVar(MeasuredTrans);

 std::cout << " did fov corrections " << std::endl;
//exit(23);
  return 0;
}