//-----------------------------------------------------------------------------------------
//  Subroutine deals with the nonlinear response of SOFIE radiometer measurements.
//
//  There are 2 options:
//
//    1) Correct (remove) nonlinearity in the input radiometer signal.
//
//    2) Make the input linear signal be nonlinear.
//
//  Input:
//     band.......band number (integer from 1-16)
//     resp.......the input response (counts), can range from 0 - 32768
//     GA.........balance attenuator setting (gain) for the measurement (can range from 0-1)
//     option.....if 1 then correct nonlinearity in resp, if 2 then make resp nonlinear
//
//  Output:
//     resp_out...option 1: response corrected for nonlinearity
//                option 2: response with nonlinearity applied 
//
//  Source: Mark Hervig, GATS
//-----------------------------------------------------------------------------------------
#include "NonLinearity.h"

#include <cassert>
#include <iostream>

//using namespace std;

double sofie_nonlinearity (const int band, const int option, const double resp, const double GA)

{

// check the input for problems

assert(band >= 1 && band <= 16);
//   if (band < 1 || band > 16) {
 //    cout << "band number out of range in sofie_nonlinearity.cpp:  " << band << endl;
  //   return resp;
//   }

//assert( resp >= 0. && resp <= 32768.0);
//   if (resp < 0 || resp > 32768.0) {
//     cout << "resp out of range in sofie_nonlinearity.cpp:  " << resp << endl;
//     return resp;
//   }
	
assert( GA >= 0. && GA <= 1.);
//   if (GA < 0.0 || GA > 1.0) {
//     cout << "ballance attenuator setting out of range in sofie_nonlinearity.cpp:  " << GA << endl;
//     return resp;
//   }

assert(option == 1 || option ==2);
//   if (option < 1 || option > 2) {
//     cout << "option is out of range in sofie_nonlinearity.cpp:  " << option << endl;
//     return resp;
//   }
	
// Cnl are the nonlinearity constants derived in ground calibration,
// see sofiedata.org for further explanation.  The balance attenuator
// setting (electronic gain, GA) for these calibration results was GAcal = 0.83

   static const double GAcal = 0.83;
   static const double Cnl[16] =   // the nonlinearity constants,  1/counts
		{
		    0.,  // band 1
		    0.,
		    0.,
		    0.,
		    1.79e-6, // band 5
		    1.56e-6,
		    9.58e-6,
		    8.55e-6,
		    0.80e-6, // band 9
		    1.60e-6,
		    1.74e-6,
		    2.56e-6,
		    5.01e-6, // band 13
		    3.15e-6,
		    1.93e-6,
		    2.20e-6 // band 16
		};

// correct nonlinearity or apply nonlinearity

   double resp_out;  // the output response

   if (option == 1) {
     resp_out = resp / (1.0 -  Cnl[band-1] * resp * GAcal / GA);  // correct nonlinearity
   }

   if (option == 2) {
     resp_out = resp / (1.0 +  Cnl[band-1] * resp * GAcal / GA);  // apply nonlinearity
   }

   return resp_out; 
}


//-----------------------------------------------------------------------------------------
//  Subroutine calculates the transmission that you would get from using nonlinear data.
//
//  The input exoatmospheric response (Rexo) should nonlinear as measured.
//
//  Input:
//     band.......band number (integer from 1-16)
//     t..........simulated atmospheric transmission
//     GA.........balance attenuator setting (can range from 0-1) 
//     Rexo.......exoatmospheric response (counts), nonlinear 
//
//  Output:
//     tsim.......simulated transmission for nonlinear data
//
//  Source: Mark Hervig, GATS
//-----------------------------------------------------------------------------------------

//#include <iostream>

//using namespace std;

double sofie_sim_nonlin_tran (const int band, const double t, const double GA, const double Rexo)
{

 //  double sofie_nonlinearity (int band, int option, double resp, double GA);

    int    option;
   double RL, R, RexoL, tsim;

// check the input for problems


assert(band >= 1 && band <= 16);
//   if (band < 1 || band > 16) {
 //    cout << "band number out of range in sofie_nonlinearity.cpp:  " << band << endl;
  //   return resp;
//   }

assert(t >=0 && t <= 1.0);

assert( Rexo >= 0. && Rexo <= 32768.0);
//   if (resp < 0 || resp > 32768.0) {
//     cout << "resp out of range in sofie_nonlinearity.cpp:  " << resp << endl;
//     return resp;
//   }
	
assert( GA >= 0. && GA <= 1.);
//   if (GA < 0.0 || GA > 1.0) {
//     cout << "ballance attenuator setting out of range in sofie_nonlinearity.cpp:  " << GA << endl;
//     return resp;
//   }

//assert(option == 1 || option ==2);
//   if (option < 1 || option > 2) {
//     cout << "option is out of range in sofie_nonlinearity.cpp:  " << option << endl;

 
// The linear exo response

   option = 1;  // 1 means correct nonlinearity  in "sofie_nonlinearity.cpp"

   RexoL  = sofie_nonlinearity(band, option, Rexo, GA);   //  linear system exo response, counts

// Calculate dV/Vs for a nonlinear system       

   RL = RexoL * t;    //  linear system response, counts

   option = 2;  // 2 means apply nonlinearity in "sofie_nonlinearity.cpp"

   R  = sofie_nonlinearity(band, option, RL, GA);   //  nonlinear response, counts

   tsim = R / Rexo;

   return tsim;
}


//-----------------------------------------------------------------------------------------
//  Subroutine calculates the SOFIE dV/Vs ([weak - strong signal] / [strong band signal])
//  based on input strong and weak band transmissions.  dV/Vs is computed for gain-corrected
//  zero-normalized dV signals,  that is, the dV is balanced to zero exoatmospherically and 
//  has no dV gain.  
//
//  The input is referenced as "band 1 and band 2", which is not literal, but rather
//  the odd then even # bands in a channel. This code computes the difference signals 
//  as weak band - strong band,  which is not always odd - even band or even - odd band.  
//
//  The routine includes the effects of nonlinearity, and returns the dV/Vs response for a 
//  nonlinear system.   The input exoatmospheric response values (Rexo1 & Rexo2)
//  should nonlinear as measured.
//
//  Input:
//     chan.......channel number (integer from 1-8)
//     t1.........atmospheric transmission in odd band
//     t2.........atmospheric transmission in even band
//     G1.........balance attenuator setting (gain) for odd band (can range from 0-1)
//     G2.........balance attenuator setting (gain) for even band (can range from 0-1)
//     Rexo1......exoatmospheric response in odd band (counts) 
//     Rexo2......exoatmospheric response in even band (counts) 
//
//  Output:
//     dV.........dV/Vs
//
//  Source: Mark Hervig, GATS
//-----------------------------------------------------------------------------------------

//#include <iostream>

//using namespace std;

double sofie_sim_nonlin_dv (const int chan, const double t1, const double t2, 
const double G1,const  double G2, const double Rexo1, const double Rexo2)
{

//   double sofie_nonlinearity (int band, int option, double resp, double GA);

   int    ibw, ibs, option;
   int    iweak[8], istrong[8];

   double tw,ts,Gw,Gs,Rexow,Rexos,tsims, tsimw;
   double RwL, RsL, Rw, Rs, dV, RexowL, RexosL;

// check the input for problems

assert(chan >=1 && chan <=8);
	
assert(G1 >=0. && G1 <= 1.0);
assert(G2 >=0. && G2 <= 1.0);
/*
   if (chan < 1 || chan >8) {
     cout << "channel number out of range in sofie_sim_dv.cpp:  " << chan << endl;
     return 0.;
   }
	
   if (G1 < 0.0 || G1 > 1.0) {
     cout << "odd band ballance attenuator setting out of range in sofie_sim_dv.cpp:  " << G1 << endl;
     return 0.;
   }

   if (G2 < 0.0 || G2 > 1.0) {
     cout << "even band ballance attenuator setting out of range in sofie_sim_dv.cpp:  " << G2 << endl;
     return 0.;
   }
*/
// tag the weak & strong band numbers (1-16) for each channel   

   iweak[0] =  2; //channel 1
   iweak[1] =  4;
   iweak[2] =  5; //channel 3
   iweak[3] =  8;
   iweak[4] = 10; //channel 5
   iweak[5] = 12;
   iweak[6] = 14; //channel 7
   iweak[7] = 15;

   istrong[0] =  1; //channel 1
   istrong[1] =  3;
   istrong[2] =  6; //channel 3
   istrong[3] =  7;
   istrong[4] =  9; //channel 5
   istrong[5] = 11;
   istrong[6] = 13; //channel 7
   istrong[7] = 16;

   ibw = iweak[chan-1];    // weak band # (1 - 16)
   ibs = istrong[chan-1];  // strong band # (1 - 16)

// assign the appropriate band data to weak and strong notation

   tw = t2;
   ts = t1;
   Rexow = Rexo2;
   Rexos = Rexo1;
   Gw = G2;
   Gs = G1;

   if (ibw % 2) {
     tw = t1;
     ts = t2;
     Rexow = Rexo1;
     Rexos = Rexo2;
     Gw = G1;
     Gs = G2;
   }

// The linear exo response

   option = 1;  // 1 means correct nonlinearity  in "sofie_nonlinearity.cpp"

   RexowL = sofie_nonlinearity(ibw, option, Rexow, Gw);   //  linear system exo response, counts
   RexosL = sofie_nonlinearity(ibs, option, Rexos, Gs);

// Calculate the Vw and Vs for a nonlinear system   

   RwL = RexowL * tw;    // weaak band linear system response, counts
   RsL = RexosL * ts;    // strong band linear system response, counts

   option = 2;  // 2 means apply nonlinearity in "sofie_nonlinearity.cpp"

std::cout << RwL << "    " << RsL <<std::endl;

   Rw = sofie_nonlinearity(ibw, option, RwL, Gw);   //  nonlinear weak band response, counts
   Rs = sofie_nonlinearity(ibs, option, RsL, Gs);   //  nonlinear strong band response, counts

// make the signals such that the initial balance condition was dV = 0

   tsimw = Rw / Rexow ; //  simulated weak band transmission, nonlinear sigs
   tsims = Rs / Rexos ; //  simulated strong   "

// calculate dV/Vs

   dV = (tsimw - tsims) / tsims;  //  dV/Vs 

   return dV;
}