#include <except.hh>
#include <su.h>
#include <segy.h>
#include <header.h>
#include <parser.h>
#include <math.h>

const char * sdoc[] = { 
  "Usage: my_wlt.x CWPROOT= wavelet= dt=",
  " ",
  "Required parameters:",
  "  CWPROOT  = path to SU root directory",
  "  output   = filename of output wavelet",
  "  type     = waveform type (0=gauss,1=dgauss,2=ricker,3=cinf)",
  "  nt       = number of time grid",
  "  dt       = time increment",
  "  cit      = time index of center",
  "  scal     = scaling factor",
  "  width    = width of cinf support",
  "  fpeak    = peak freq for gauss,dgauss,ricker",
NULL};  

using RVL::RVLException;
using RVL::valparse;
using std::string;

int xargc;
char **xargv;

//--------------------------------------------------------------------
float my_cinf(float dt, int it, int it0, int itF){
//--------------------------------------------------------------------
  if( it<=it0 || it>=itF )
    return 0.0;
  else{
    float h = (itF+it0)*0.5;
    float tmp = exp( -1/( dt*dt*(it-it0)*(itF-it) ) );
    tmp *= exp( 1/( dt*dt*(h-it0)*(itF-h) ) );
    return tmp;
  }
}

//--------------------------------------------------------------------
float my_gauss(float dt, int it, int cit, int fpeak){
//--------------------------------------------------------------------
  return exp( -dt*dt*(it-cit)*(it-cit)*PI*PI*fpeak*fpeak );
}

//--------------------------------------------------------------------
float my_dgauss(float dt, int it, int cit, int fpeak){
//--------------------------------------------------------------------
  return -2*dt*(it-cit)*PI*PI*fpeak*fpeak*my_gauss(dt,it,cit,fpeak);
}

//--------------------------------------------------------------------
float my_ricker(float dt, int it, int cit, int fpeak){
//--------------------------------------------------------------------
  return (1-2*dt*dt*(it-cit)*(it-cit)*PI*PI*fpeak*fpeak)*my_gauss(dt,it,cit,fpeak);
}

//--------------------------------------------------------------------
int main(int argc, char ** argv) {
//--------------------------------------------------------------------
  try {

    xargc=argc; xargv=argv;
    requestdoc(1);

    PARARRAY * par = ps_new();
    if ( ps_createargs(par, argc - 1, argv + 1) ) {
      printf("Error parsing input data. ABORT.\n");
      exit(1);
    }
    //ps_printall(*par,stderr);

    string cwp    = valparse<string>(*par,"CWPROOT");
    string output = valparse<string>(*par,"output");
    int    type   = valparse<int>(*par,"type");
    int    nt     = valparse<int>(*par,"nt");  
    float  dt     = valparse<float>(*par,"dt"); //units [s]
    
    int cit;
    try{ cit = valparse<int>(*par,"cit") -1;}
    catch(RVLException &e){cit = 0;}

    float scal;
    try{ scal = valparse<float>(*par,"scal");}
    catch(RVLException &e){scal = 1.0;}

    int fpeak;
    try{ fpeak = valparse<int>(*par,"fpeak");}
    catch(RVLException &e){
      if(type>=0&&type<=2){
	e << "fpeak is required for type="<< type <<"\n";
	throw;
      }
    }

    float width; //units [s]
    try{ width = valparse<float>(*par,"width");}
    catch(RVLException &e){width=0.5;}
    int it0 = cit - (int)(width*0.5/dt)-1;
    int itF = cit + (int)(width*0.5/dt)+1;

    //creating su file    
    string sunull = cwp + "/bin/sunull";
    stringstream cmd;
    cmd << sunull << " ntr=1 nt="<< nt <<" dt="<< dt
	<< " >" << output; 
    system(cmd.str().c_str());

    //reading in trace
    FILE* fp;
    segy tr;
   
    fp = fopen(output.c_str(),"r"); 
    fvgettr(fp,&tr);
    fclose(fp);

    fp = fopen(output.c_str(),"w");

    //writing out traces
    if(type==0){
      for(int it=0; it<nt; it++) tr.data[it]=my_gauss(dt,it,cit,fpeak)*scal;
    } 
    else
    if(type==1){
      for(int it=0; it<nt; it++) tr.data[it]=my_dgauss(dt,it,cit,fpeak)*scal;
    }
    else
    if(type==2){
      for(int it=0; it<nt; it++) tr.data[it]=my_ricker(dt,it,cit,fpeak)*scal;
    }
    else
    if(type==3){
      for(int it=0; it<nt; it++) tr.data[it]=my_cinf(dt,it,it0,itF)*scal;
    }
    else{
      RVLException e;
      e << "type must be <= 3!\n";
      throw e;
    }

    fvputtr(fp,&tr);
    fclose(fp);
  }
  catch (RVLException & e) {
    e << "Error caught in main for my_waveform.cc\n";
    e.write(cerr);
    exit(1);
  }
}