# include "System.h"
# include "Syslog.h"
# include "HS2Cal.h"

# include <stdio.h> 
# include <stdlib.h>
# include <math.h>

double HS2Cal::s_calTemp = 25;

double HS2Cal::s_gain[] = {
  1, 9.4573, 94.65, 932.1, 9497.5,
  1, 9.4627, 94.348, 926.66, 9684.3,
  1, 10, 100, 1000, 10000
};

double HS2Cal::s_mu[] = {52.13, 16.37, 10};

double HS2Cal::s_tempCoeff[] = {0, 0, 0};

double HS2Cal::s_Rnominal[] = {8000, 8000, 8000};

double HS2Cal::s_beta2bb[] = {6.79, 6.79, 6.79};

double HS2Cal::s_sigmaExp[] = {.146, .159, 0};

double const HS2Cal::C = 0.1;
double const HS2Cal::gamma_y = 0.014;
double const HS2Cal::a_d_400 = 0.01;
double const HS2Cal::gamma_d = 0.011;
double const HS2Cal::a_star[] = {
  400,0.687,
  410,0.828,
  420,0.913,
  430,0.973,
  440,1,
  450,0.944,
  460,0.917,
  470,0.87,
  480,0.798,
  490,0.75,
  500,0.668,
  510,0.618,
  520,0.528,
  530,0.474,
  540,0.416,
  550,0.357,
  560,0.294,
  570,0.276,
  580,0.291,
  590,0.282,
  600,0.236,
  610,0.252,
  620,0.276,
  630,0.317,
  640,0.334,
  650,0.356,
  660,0.441,
  670,0.595,
  680,0.502,
  690,0.329,
  700,0.215,
  0, 0
};

double const HS2Cal::beta_w_ref = 2.18e-4;
double const HS2Cal::b_bw_ref = 1.17e-3;
double const HS2Cal::lambda_ref = 525;
double const HS2Cal::gamma = 4.32;
double const HS2Cal::b_b_tilde = 0.015;
double const HS2Cal::k_1 = 1.0;

HS2Cal::HS2Cal() {
  loadCfg(System::configurationFile("hs2Calibration.dat"));
}

double HS2Cal::typ_absorption(double freq) {
  double pfreq = 0, pfact = 0;
  
  for(double const *i = a_star; 0!=*i && *i<freq; i+=2) {
    pfreq = *i;
    pfact = i[1];
  }
  double star_a;
  if( *i==freq )
    star_a = i[1];
  else {
    star_a = (freq-pfreq)*(i[1]-pfact)/(*i-pfreq);
    star_a += pfact;
  }
  return (0.06*star_a*std::pow(C, 0.65))*(1+0.2*std::exp(-gamma_y*(freq-440)))
    + (a_d_400*std::exp(-gamma_d*(freq-400)));
}

double HS2Cal::calc_uncorr(short channel, short snorm, short gainStatus, 
			   double temp) {
#if 0
  if( channel>3 || channel<1 ) {
    Syslog::write("Cluster:HS2Cal: channel has to be in [1 3].");
    return NAN;
  }
  if( gainStatus>5 || gainStatus<1 ) {
    Syslog::write("Cluster:HS2Cal: gainStatus has to be in [1 5].");
    return NAN;
  }
#endif // DEBUG
  double beta_uncorr = snorm*s_mu[channel-1];

  beta_uncorr /= (1+s_tempCoeff[channel-1]*(temp-s_calTemp)) 
    *s_gain[5*(channel-1)+(gainStatus-1)]*s_Rnominal[channel-1];
  return beta_uncorr;
} // HS2Cal::calc_uncorr

double HS2Cal::calc_bb(short channel, short snorm, short gainStatus, 
		       double temp, short freq) {
  double beta_uncorr = calc_uncorr(channel, snorm, gainStatus, temp), 
    b_b_uncorr;
#if 0
  if( isnan(beta_uncorr) ) {
    Syslog::write("Cluster:HS2Cal: NAN value returned by calc_uncorr");
    return NAN;
  }
#endif // DEBUG
  double tmp = std::pow(lambda_ref/freq, gamma);
  double beta_w = beta_w_ref*tmp;
  double b_bw = b_bw_ref*tmp;
  
  b_b_uncorr = s_beta2bb[channel-1]*(beta_uncorr-beta_w)+b_bw;
  
  double a = typ_absorption(freq);
  double b = (b_b_uncorr - b_bw)/b_b_tilde;
  
  double k_bb = a+0.4*b;
  double k_exp = s_sigmaExp[channel-1];
  double sigma = k_1*std::exp(k_exp*k_bb);

  return sigma*b_b_uncorr-b_bw;
} // HS2Cal::calc_bb

void HS2Cal::patchShort(short &value) {
  if( value>32768 ) 
    value = 65536-value;
} // HS2Cal::patchShort

void HS2Cal::loadCfg(char const *filename) {
  Syslog::write("hydroscat:hs2cal parsing %s", filename);
  FILE *cfg = fopen(filename, "r");
  char buff[256];
  buff[255] = '\0';

  if( NULL!=cfg ) {
    fgets(buff, 256, cfg); // [General] - Trash
    fgets(buff, 256, cfg); // DeviceType=
// #ifdef DEBUG
//     Syslog::write("cluster:HS2Cal %s", buff);
// #endif 
    fgets(buff, 256, cfg); // Serial=
    fgets(buff, 256, cfg); // Label=
    fgets(buff, 256, cfg); // Config=
    fgets(buff, 256, cfg); // MaxDepth=
    fgets(buff, 256, cfg); // CalTime=
// #ifdef DEBUG
//     Syslog::write("cluster:HS2Cal %s", buff);
// #endif
    fgets(buff, 256, cfg); // DepthCal=
    fgets(buff, 256, cfg); // DepthOff=
    fgets(buff, 256, cfg); // CalTemp=
    s_calTemp = strtod(buff+8, NULL);
// #ifdef DEBUG
//     Syslog::write("cluster::HS2Cal CalTemp=%g", s_calTemp);
// #endif // DEBUG
    fgets(buff, 256, cfg); // 
    for( int chan=0; chan<3; ++chan ) {
      fgets(buff, 256, cfg); // [Channel 1:3]
      fgets(buff, 256, cfg); // Name=
// #ifdef DEBUG
//       Syslog::write("cluster::HS2Cal Channel.%s", buff);
// #endif 
      for(int gain=0; gain<5; ++gain) {
	fgets(buff, 256, cfg); // GainX=
	s_gain[5*chan+gain] = strtod(buff+6, NULL);
// #ifdef DEBUG
// 	Syslog::write("cluster::HS2Cal    Gain[%d]=%g", 
// 		      gain+1, s_gain[5*chan+gain]);
// #endif
      }
      fgets(buff, 256, cfg); // Mu=
#if 1
      if(chan<2) 
#endif // patch to have mu set to 50 for fl676
      s_mu[chan] = strtod(buff+3, NULL);
#if 1
      else
	s_mu[chan] = 50;
#endif // patch to have mu set to 50 for fl676
// #ifdef DEBUG
//       Syslog::write("cluster::HS2Cal    Mu=%g", s_mu[chan]);
// #endif
      for(int i=0; i<3; ++i )
	fgets(buff, 256, cfg); // SigmaX=
      fgets(buff, 256, cfg); // SigmaExp=
      s_sigmaExp[chan] = strtod(buff+9, NULL);
// #ifdef DEBUG
//       Syslog::write("cluster::HS2Cal    SigmaExp=%g", s_sigmaExp[chan]);
// #endif
      fgets(buff, 256, cfg); // RNominal=
      s_Rnominal[chan] = strtod(buff+9, NULL);
// #ifdef DEBUG
//       Syslog::write("cluster::HS2Cal    RNominal=%g", s_Rnominal[chan]);
// #endif
      fgets(buff, 256, cfg); // Beta2Bb=
      s_beta2bb[chan] = strtod(buff+8, NULL);
// #ifdef DEBUG
//       Syslog::write("cluster::HS2Cal    Beta2Bb= %g", s_beta2bb[chan]);
// #endif
      fgets(buff, 256, cfg); // MuRho=
      fgets(buff, 256, cfg); // TempCoeff=
      s_tempCoeff[chan] = strtod(buff+10, NULL);
// #ifdef DEBUG
//       Syslog::write("cluster::HS2Cal    TempCoeff=%g", s_tempCoeff[chan]);
// #endif
      for(i=0; i<5; ++i) 
	fgets(buff, 256, cfg); // OffsetX=
      for(i=0; i<5; ++i) 
	fgets(buff, 256, cfg); // NomOffsetX=
      fgets(buff, 256, cfg); //
    }
    fclose(cfg);
  }
}

HS2Cal::Output HS2Cal::calibrate(HydroscatIF::Data &data) const {
  Output res;
  
  res.ready = data.deviceReady;

  if( res.ready ) {
    patchShort(data.snorm1);
    patchShort(data.snorm2);
    patchShort(data.snorm3);
#ifdef DEBUG
    Syslog::write("Cluster:HS2Cal : input is");
    Syslog::write("\tbb470 = [%d %d]", data.snorm1, data.gainstatus1);
    Syslog::write("\tbb676 = [%d %d]", data.snorm2, data.gainstatus2);
    Syslog::write("\tfl676 = [%d %d]", data.snorm3, data.gainstatus3);
    Syslog::write("\ttemperature = %g", data.temp);
#endif // DEBUG
    
    res.bb470 = calc_bb(1, data.snorm1, data.gainstatus1, data.temp, 470);
    res.bb676 = calc_bb(2, data.snorm2, data.gainstatus2, data.temp, 676);
    res.fl676_uncorr = calc_uncorr(3, data.snorm3, data.gainstatus3, data.temp);

#ifdef DEBUG
    Syslog::write("Cluster:HS2Cal : calibrated output is");
    Syslog::write("\tbb470 = %g", res.bb470);
    Syslog::write("\tbb676 = %g", res.bb676);
    Syslog::write("\tfl676_uncorr = %g", res.fl676_uncorr);
#endif // DEBUG

  }
  return res;
} // HS2Cal::operator()