source: EEW/trunk/mcast2eew/EEW.cc @ 3747

#include "EEW.h"
#include "eewPropertiesST.h"
#include "Duration.h"
#include "../eewserver/udp.h"
#include "../eewserver/MessageBuffer.h"
#include <iomanip>
#include <arpa/inet.h>
#include <netdb.h>

unsigned int EEW::serverip = 0;
bool EEW::udpsetup = false;

EEW::EEW(Channel Z,Channel E,Channel N){

  chanZ = Z;
  chanE = E;
  chanN = N;

  samplerate = chanZ.samprate;

  gcbrZ = NULL;
  gcbrN = NULL;
  gcbrE = NULL;
  
  picker = NULL;
  buthfilter = NULL;
  
  buf = NULL;
  pickbuf = NULL;
  vdbuf = NULL;
  trigger = false;
  firsttimeE = true;
  firsttimeN = true;
  firsttimeZ = true;
  
  TC = 0.0;
  Pd3 = 0.0;
  Q = 0.0;
  mag = 0.0;
  pgv = 0.0; 
  clipped = 0; 
  
  init();
}

void EEW::init(){

  if(udpsetup==false){
    const int  hostlen = 128;
    char hostname[hostlen];

    int udp_port = 30000+(getpid()%1000);
    if(gethostname(hostname,hostlen)<0){
      perror("Can't get hostname");
      exit(-1);
    }

    udp::getInstance()->init_udp(udp_port,(char*)gethostnametoip(hostname));
    eewPropertiesST* prop = eewPropertiesST::getInstance();
    const char* dst = gethostnametoip((char*)prop->getServerName().c_str());
    if(dst){
      serverip = ntohl(inet_addr(dst));
    }
    else{
      cout <<"Error: Can't find IP Address for the Server Host Name"<<endl;   
      exit(-1);
    }
    udpsetup = true;
  }

  gcbrN = new GCBRTSPModule(chanN);
  gcbrN->load("");

  gcbrE = new GCBRTSPModule(chanE);
  gcbrE->load("");

  gcbrZ = new GCBRTSPModule(chanZ);
  gcbrZ->load("");

  picker = new PPickerTSPModule(chanZ);
  picker->setEnvHandle(this);
  picker->load("");
  
  buthfilter = new ButhpinFilterTSPModule(chanZ);
  buthfilter->load("");

  cout <<"SampleRate for "<<chanZ<<" is "<<samplerate<<endl;

  buf = new float[samplerate*10];
  pickbuf = new float[samplerate*10];
  vdbuf = new struct vd[samplerate*10];
}


void EEW::process(Channel chn,timeval stime,int* samps,int size){
  // cout <<"Processing "<<chn<<endl;
  unsigned int outsize = 0;
  eewPropertiesST* prop = eewPropertiesST::getInstance();
  TimeStamp starttime = TimeStamp(stime);

  if((double)(TimeStamp::current_time() - starttime) > 20.0){
    return;
  }

  if(string(chn.channel) == string(chanZ.channel)){
    checkgap(starttime,lasttimeZ,size,firsttimeZ);

    gcbrZ->process(stime,size,samps,outsize,buf);
    maintainGCBRData(zq,zmap,stime,size);
    buthfilter->process(stime,size,buf,outsize,vdbuf);  
    picker->process(stime,size,buf,outsize,pickbuf);
    maintainVDQueue(vdq,stime,size);    
  }
  else if(string(chn.channel) == string(chanN.channel)){
    checkgap(starttime,lasttimeN,size,firsttimeN);
    gcbrN->process(stime,size,samps,outsize,buf);
    maintainGCBRData(nq,nmap,stime,size);
  }
  else if(string(chn.channel) == string(chanE.channel)){
    checkgap(starttime,lasttimeE,size,firsttimeE);
    gcbrE->process(stime,size,samps,outsize,buf);
    maintainGCBRData(eq,emap,stime,size);
  }

  //TimeStamp endtime = TimeStamp::current_time();
  //cout <<"Processing Delay : "<<(double)(endtime - begtime)<<endl;
}

void EEW::checkgap(TimeStamp firstsamptime,TimeStamp& lastsamptime,int nsamps,bool& firsttime){
  float dt = 1.0/samplerate;

  TimeStamp temptime;
  float offset = dt*nsamps;
  struct timeval time_offset = {0,0};
  time_offset.tv_sec = 0;
  time_offset.tv_usec = offset*1000000;
  temptime = firstsamptime + Duration(time_offset); 

  if(firsttime){
    lastsamptime = temptime;
    firsttime = false;
  }
  else{
    double timediff = (double)(firstsamptime - lastsamptime);
    //    cout <<"TimeDiff: "<<timediff<<" dt:"<<dt<<endl;
    if(timediff > dt){
      //gap
      lastsamptime = temptime;
      reset();
      firsttime = true;
    }
    else{
      lastsamptime = temptime;
    }
  }
}

void EEW::reset(){
  //  cout <<"Gap Found. RESETTING EEW"<<endl;

  zq.clear();
  eq.clear();
  nq.clear();

  zmap.clear();
  emap.clear();
  nmap.clear();

  vdq.clear();

  struct timeval dummy = {0,0};

  gcbrZ->restart(dummy);
  picker->restart(dummy);
  buthfilter->restart(dummy);

  gcbrE->restart(dummy);
  gcbrN->restart(dummy);
}

void EEW::clear(){
  cout <<"EEW STOPPED"<<endl;
  delete gcbrZ;
  delete picker;
}

void EEW::maintainGCBRData(GCBRQueue& dataq,GCBRMap& datamap,timeval stime,int size){
       float dt = (float)1.0/(samplerate*1.0);
       
       for(int j=0;j<size;j++){
         float offset = dt*j;
         struct timeval time_offset = {0,0};
         time_offset.tv_sec = 0;
         time_offset.tv_usec = offset*1000000;
         
         TimeStamp timestamp = TimeStamp(stime) + Duration(time_offset);
         dataq.push_back(make_pair(timestamp,buf[j]));
         if(dataq.size() > MAXQ){
           dataq.pop_front();
         }
       }
}


void EEW::maintainVDQueue(VDQueue&,timeval stime,int size){
       float dt = (float)1.0/(samplerate*1.0);
       //       cout <<"In maintainVDQueue .."<<endl;
       for(int j=0;j<size;j++){
         float offset = dt*j;
         struct timeval time_offset = {0,0};
         time_offset.tv_sec = 0;
         time_offset.tv_usec = offset*1000000;
         TimeStamp timestamp = TimeStamp(stime) + Duration(time_offset);    
         vdq.push_back(make_pair(timestamp,vdbuf[j]));
       }

       if(trigger){
         // cout <<"In maintainVDQueue.. in Trigger...";
         //pop samples while we find trigger
         int remcount = 0;
         while(vdq.size() > 0){
           TimeStamp t = (*(vdq.begin())).first;
           //cout <<"VDQ TimeStamp: "<<t<<" TrigTime: "<<trigtime<<endl;
           if(t.seconds() == trigtime.seconds() && t.milliseconds() == trigtime.milliseconds()){
             break;
           }
           else{
             remcount++;
             vdq.pop_front();
             //cout <<"Removing Sample from VDQ"<<endl;
           }
         }
         //cout <<" VDQ size = "<<vdq.size()<<" and removed "<<remcount<<" elements"<<endl;
         if(vdq.size() >= samplerate*3.0){
           //Calculate TauC/PD
           cout <<"VDQ size = "<<vdq.size()<<" Calculate TauC/PD"<<endl;
           calcTCPD();
         }
       }
       else{
         vdq.clear();
       }
}

void EEW::setTriggerOn(TimeStamp timestamp){
  trigger = true;
  trigtime = timestamp;
}

void EEW::calcTCPD(){

  trigger = false;

  eewPropertiesST* prop = eewPropertiesST::getInstance();
  Channel channel = chanZ; //Channel(prop->getNetwork().c_str(),prop->getStation().c_str(),"  ",string(string(prop->getChannelType())+"Z").c_str());
  
  clipped = 0;

  //Get Trigger Time Samples from E and N Channels
  float eout = findSample(eq,trigtime);
  float nout = findSample(nq,trigtime);
  float zout = findSample(zq,trigtime);
  if(eout == 0.0 || nout == 0.0 || zout == 0.0){
    double delay = TimeStamp::current_time() - TimeStamp(trigtime);
    cout <<"ERROR: E/N/Z Channels "<<chanZ<<" DO NOT HAVE TRIGGER SAMPLES. ABORTING TauC/PD Calculation For now, After waiting "<<delay<<" seconds"<<endl;
    cout <<"Amp for N = "<<nout<<" and E = "<<eout<<" and Z = "<<zout<<endl; 
    cout <<"Size for GCBRMaps : N = "<<nmap.size()<<" E = "<<emap.size()<<" Z = "<<zmap.size()<<endl;
    trigger = false;
    return;
  }
  //

    int nsamples = samplerate * 3.0;
    struct vd* out = new struct vd[nsamples];

    for(int i=0;i<nsamples;i++){
      out[i] = (*(vdq.begin())).second;
      vdq.pop_front();
    }
    
//     //print data
//     cout <<"--- Z Displacement ---"<<endl;
//     for(int i=0;i<nsamples;i++){
//       cout <<out[i].d<<",";
//     }
//     cout <<endl;

    //Check clipping
    if(isClipping(out,nsamples)){
      clipped = 1;
      clipbegin = TimeStamp(trigtime);
    }
    else{
      double lastclipdur = (double)(TimeStamp(trigtime) - clipbegin);
      double clipwaitsec = (double)prop->getClipWaitSeconds();
      
      cout <<channel<< " lastclipdur:"<<lastclipdur<<" clipwaitsec:"<<clipwaitsec<<endl;
      
      if(lastclipdur <= clipwaitsec){
        clipped = 1;
      }
    }
    
    TC =  calTauC(out,nsamples);
    Pd3 = calPd3(out,nsamples);
    delete [] out;

    cout <<"Total Latency for "<<channel<<" is "<<(double)(TimeStamp::current_time() - TimeStamp(trigtime))<<" sec"<<endl;
    cout <<"Amp for N = "<<nout<<" and E = "<<eout<<" and Z = "<<zout<<endl;


    //Check if E,N,Z components satisfy the criteria//
    double X = fabs(sqrt(pow(eout,2) + pow(nout,2))/zout);
    if(X > 1.5){
      cout <<"REJECTED "<<channel<<" Reason: Not a p-wave."<<endl;
      cout <<"The trigger ignored for "<<channel<<" PD: "<<Pd3<<", X = "<<X<<", Amp for N = "<<nout<<" and E = "<<eout<<" and Z = "<<zout<<endl;
      return;
    }

    Q = tauc_pd_trigger(TC,Pd3);
    if(Q == 0){
      cout <<"REJECTED "<<channel<<" Reason: Q = 0  Trigger Time: "<<TimeStamp(trigtime)<<" TauC = "<<TC<<" PD = "<<Pd3<<endl;
      return;
    }

    mag = 4.218*log10(TC)+6.166;
    pgv = pow(10,0.920*log10(Pd3)+1.642);

    cout <<channel<<": TauC = "<<TC<<" and Pd3 = "<<Pd3<<" Q: "<<Q<<" Clipped: "<<clipped<<" Mag: "<<mag<<" PGV: "<<pgv<<endl;
    sendToServer();
}

double EEW::calTauC(struct vd* data,unsigned int size){
    double s1=0.0;
    double s2=0.0;
    double tauc[300]; //temp debug

    for(int i=0;i<size;i++){
        s1=s1+pow(data[i].d,2);
        s2=s2+pow(data[i].v,2);
        tauc[i] = 2.0*PI*sqrt(s1/s2);
    }

    double TauC = 2.0*PI*sqrt(s1/s2);
    return TauC;
}


double EEW::calPd3(struct vd* data,unsigned int size){
    double pd3 = 0.0;
    for(int i=0;i<size;i++){
        double temp = fabs(data[i].d);
        if(temp > pd3)
            pd3 = temp;
    }
    return pd3;
}

bool EEW::isClipping(struct vd* data, unsigned int size){
  eewPropertiesST* prop = eewPropertiesST::getInstance(); 
  float valthresh = prop->getVelocityThreshold();
  for(int i=0;i<size;i++){
      if(abs(data[i].v) >= 0.8){
          return true;
      }
  }
  return false;
}

int EEW::tauc_pd_trigger(double tau_c,double pd){

  // Tauc-Pd trigger criterion after Boese et al. (in prep., Feb. 2008)://
  // Earthquake Early Warning in California: One Year Real-Time Testing of 

  double rmin           =       1.0;            //km
  double rmax           =       100.0;          //km
  double min_tauc       =       0.2;
  double max_tauc       =       8.00;    
  double Pd_thres       =       0.0005; // cm
  double q = 0.0;
  double pi = 3.141592;

 // Coefficients taken from Cua and Heaton (in prep.), for rock sites
double a                =    0.89;
double b                =   -4.3*pow(10.0,-4.0);
double c1               =    1.11;
double c2               =    1.11;
double d                =   -1.44;
double e                =   -2.60;
double f                =    1.10;              // rms_to_max_PGV      
double R_min            =   sqrt(rmin*rmin+9); 
double R_max            =   sqrt(rmax*rmax+9);

// Uncertainties in the parameters/relations:
 double sigma_IM        =   0.280; // Cua and Heaton (in prep.)
 double sigma_M     =   0.385; // Wu et al., 2007: GJI
 double sigma_PGV       =   0.326; // Wu and Kanamori, 2008: Sensors
    
 // Determine trigger quality Q:
 if((tau_c<min_tauc)||(tau_c>max_tauc)){
        return 0.0;
 }
 else{
   // Estimate magnitude from tauc (Wu et al., 2007: GJI)
   double M_est         =   4.218*log10(tau_c)+6.166;  
    
   // Mean values of the Tauc-Pd trigger criterion
   double CM1         =   c1*exp(c2*(M_est-5))*(atan(M_est-5)+pi/2);
   double Pd_min1     =   pow(10,(log10(pow(10,(a*M_est+b*(R_max+CM1)+d*log10(R_max+CM1)+e))*f) -1.642)/0.92);
   double Pd_max1     =   pow(10,(log10(pow(10,(a*M_est+b*(R_min+CM1)+d*log10(R_min+CM1)+e))*f) -1.642)/0.92);

   // Mean + std dev. values of the Tauc-Pd trigger criterion   
   double CM2         = c1*exp(c2*(M_est-sigma_M-5))*(atan(M_est-sigma_M-5)+pi/2);
   double CM3         = c1*exp(c2*(M_est+sigma_M-5))*(atan(M_est+sigma_M-5)+pi/2);
   double Pd_min2     = pow(10,(log10(pow(10,(a*(M_est-sigma_M)+b*(R_max+CM2)+d*log10(R_max+CM2) +e-sigma_IM))*f)-1.642-sigma_PGV)/0.92);
   double Pd_max2     = pow(10,(log10(pow(10,(a*(M_est+sigma_M)+b*(R_min+CM3)+d*log10(R_min+CM3) +e+sigma_IM))*f)-1.642+sigma_PGV)/0.92);

   // Determine quality Q of the trigger:
    if ((pd>=Pd_thres)&&(pd>=Pd_min1)&&(pd<=Pd_max1)){
        q       =       1.0;
    }
    else if((pd>=Pd_thres)&&((pd>=Pd_min2)&&(pd<=Pd_min1)) ||((pd>=Pd_max1)&&(pd<=Pd_max2))){
        q       =       0.5;
    }
    else{
        q       =       0.0;
    }
  }
  
  return q;
}


void EEW::generateRandomTrigger(TimeStamp stime){
  static TimeStamp lasttrig = stime;

  cout <<"********  Trigger In "<<10 - (double)(stime - lasttrig)<<endl;
  if((double)(stime - lasttrig) >= 10.0){
    lasttrig = stime;
    setTriggerOn(stime);
  }
}

float EEW::findSample(GCBRQueue& q,TimeStamp timestamp){
  for(GCBRQueue::iterator it=q.begin();it!=q.end();++it){
    TimeStamp t = (*it).first;
    if(t.seconds() == timestamp.seconds() && t.milliseconds() == timestamp.milliseconds()){
      return (*it).second;
    }
  } 
  return 0.0;
}

void EEW::sendToServer(){
  eewPropertiesST* prop = eewPropertiesST::getInstance();
    char str[256];    
    long sec = trigtime.seconds();
    long usec = trigtime.micro_seconds();

    sprintf(str,"%s %s %s -- %ld %ld %f %f %f %f %f %d",chanZ.network,
                                                        chanZ.station,
                                                        chanZ.channel,
                                                        sec,usec,TC,Pd3,mag,pgv,Q,clipped);  
    MessageBuffer *mbuf = new MessageBuffer();
    mbuf->clear();
    mbuf->append(str,strlen(str));
    
    try{
      udp* udp_hnd = udp::getInstance();
      udp_hnd->send(mbuf,serverip,prop->getServerPort());
      cout <<"UDP Message Sent"<<endl;
    }
    catch(Error e){
      cout<<"Error while sending message"<<endl;
    }
    delete mbuf;   
}

const char* EEW::gethostnametoip(const char* hostname){
  if(!hostname)
    return NULL;

  struct hostent* ent = gethostbyname(hostname);
  if(!ent)
    return NULL;

  if(ent->h_addrtype == AF_INET){
    struct in_addr adr;
    bcopy(ent->h_addr_list[0],(char *) &adr, sizeof(adr));
    return inet_ntoa(adr);
  }
  else{
    return NULL;
  }
}