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

#include "pntsrc_appx.hh"
#include "MPS_includes.hh"

const char * sdoc[] = { 
  "Usage: pntsrc_appx.x CWPROOT= src= wavelet= appx_ord= dervx= dervz=",
  "       nt= ot= dt= sx= sz=",
  "       nx= nz= ox= oz= dx= dz=",
  " ",
  "Required parameters:",
  "  CWPROOT        = path to SU root directory",
  "  src            = filename for output source traces",
  "  wavelet        = filename for input source wavelet (single trace)",
  "  appx_ord       = desired approximation order of point source",
  "  dervx          = number of x-derivatives",
  "  dervz          = number of z-derivatives",
  "  nt             = number of time samples in wavelet",
  "  ot             = time of first sample in wavelet (ms)",
  "  dt             = time step size of wavelet (ms)",
  "  sx,sz          = source coordinates (m)",
  "  nx             = number of grid points in model (x-axis)",
  "  nz             = number of grid points in model (z-axis)",
  "  ox             = x coord of first grid point in model (m)",
  "  oz             = z coord of first grid point in model (m)",
  "  dx             = x coord grid size (m)",
  "  dz             = z coord grid size (m)",
NULL};  

using TSOpt::pntsrc_weight;
using RVL::valparse;

int xargc;
char **xargv;

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

    /*
    //Testing pntsrc_weight//
    float h = 2.0;
    float alpha =0.2;
  
    for( int s=0; s<=1; s++){
      for( int p=1; p<=4; p++){
	cerr << "(p,s)=("<<p<<","<<s<<")\n";
	int Q = p+s-1;
	for (int indx=0; indx<=Q; indx++){
	  float w = pntsrc_weight(s,p,indx,alpha,h);
	  cerr << "d["<<indx<<"]="<<w<<"\n";
	}
      }
    }
    */

    cerr << "\n////////////////////////////////////////////////////////////////////\n"
	 << "RUNNING: PNTSRC_APPX.X\n"
	 << "////////////////////////////////////////////////////////////////////\n";

    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);

    //CWPROOT
    std::string cwp = valparse<std::string>(*par,"CWPROOT");
    
    //wavelet information
    std::string wavelet = valparse<std::string>(*par,"wavelet");
    //int   nt = valparse<int>  (*par,"nt");  
    //float ot = valparse<float>(*par,"ot");
    float dt = valparse<float>(*par,"dt"); dt*=1e-3; //convert to sec.
    float sx = valparse<float>(*par,"sx");
    float sz = valparse<float>(*par,"sz");

    //model information
    float ox = valparse<float>(*par,"ox");
    float oz = valparse<float>(*par,"oz");
    float dx = valparse<float>(*par,"dx");
    float dz = valparse<float>(*par,"dz");
    int   nx = valparse<int>  (*par,"nx");
    int   nz = valparse<int>  (*par,"nz");

    //taking precautions for negative z-coord
    if( oz<0 ) oz *= -1;
    if( sz<0 ) sz *= -1;

    //parameters
    std::string src = valparse<std::string>(*par,"src");
    int appx_ord    = valparse<int>(*par,"appx_ord");
    int dervx       = valparse<int>(*par,"dervx");
    int dervz       = valparse<int>(*par,"dervz");

    /*
    // extract time step from source wavelet, convert unit from mus to s
    FILE * fp = fopen(src.c_str(),"r");
    segy tr;
    fgettr(fp,&tr);
    // dt in ms
    float dt = tr.dt*1.e-3;
    */

    //set up paths to SU commands
    std::string suplane = cwp + "/bin/suplane";
    std::string sushw   = cwp + "/bin/sushw";
    std::string suconv  = cwp + "/bin/suconv";
    std::string sugain  = cwp + "/bin/sugain";

    //finding mdl grid point index closest to src point
    int jx = (int)((sx-ox)/dx);
    int jz = (int)((sz-oz)/dz);
    float ax = (sx-ox-jx*dx)/dx;
    float az = (sz-oz-jz*dz)/dz;

    cerr << "(jx,jz) = ("<<jx<<","<<jz<<")\n";
    cerr << "(ax,az) = ("<<ax<<","<<az<<")\n";
    if( jx<0 || jx>nx ){
      RVLException e;
      e << "ERROR from pntsrc_appx.x: source x-coord is out of bounds!";
      e << " sx="<<sx<<", jx="<<jx<<"\n";
      throw e;
    }
    if( jz<0 || jz>nz ){
      RVLException e;
      e << "ERROR from pntsrc_appx.x: source z-coord is out of bounds!";
      e << " sz="<<sz<<", jz="<<jz<<"\n";
      throw e;
    }
    if( az<0 || az<0 ){
      RVLException e;
      e << "ERROR from pntsrc_appx.x: relative source position is <0!";
      e << " az="<<az<<", ax="<<ax<<"\n";
      throw e;
    }

    //grid point info of approx src point
    int pt_Nx = appx_ord + dervx;
    int pt_Nz = appx_ord + dervz;

    if(pt_Nx<=0||pt_Nz<=0){
      RVLException e;
      e << "ERROR from pntsrc_appx.x: stencil size of pnt appx in x/z-axis is <=0!";
      e << "(pt_Nx,pt_Nz) = (" << pt_Nx <<","<< pt_Nz <<")\n";
      throw e;
    }
      
    int ellx = jx - (int)((pt_Nx-1)/2);
    int ellz = jz - (int)((pt_Nz-1)/2);
    
    float pt_ox = ox + ellx*dx;
    float pt_oz = oz + ellz*dz;
      
    float pt_x  = 0.0;
    float pt_z  = 0.0;
    float w     = 0.0;
    int   count = 0;

    for (int ix=0; ix<pt_Nx; ix++){
      for (int iz=0; iz<pt_Nz; iz++){

	//computing weights
	w  = pntsrc_weight(dervx,appx_ord,ix,ax,dx);
	w *= pntsrc_weight(dervz,appx_ord,iz,az,dz);
       
	pt_x = pt_ox + ix*dx;
	pt_z = pt_oz + iz*dz;
	
	//cmd for writing source .su output
	std::stringstream cmd;
	//cmd << suplane <<" ntr=1 npl=1 nt="<< nt <<"dt="<< dt <<"ln1=1 ct1=1 cx1=1 | ";
       
	cmd << "<" << wavelet << " ";
	cmd << sugain  <<" scale="<< w;
	cmd <<" | "<< sushw <<" key=gelev a=-"<< pt_z;
	cmd <<" | "<< sushw <<" key=gx a="<< pt_x;
	cmd <<" | "<< sushw <<" key=tracr a="<< count;
	cmd <<" | "<< sushw <<" key=tracl a="<< count;
      	cmd <<" >>"<< src <<";";

	//cmd << suconv  <<" sufile="<< wavelet <<" >>"<< src <<";";
	
	cerr << "pt(ix="<<ix<<",iz="<<iz<<") = ("<<pt_x<<","<<pt_z<<")\n";
	cerr << "w=" << w << "\n";

	cerr << cmd.str() << "\n";
	system(cmd.str().c_str());
	
	count++;
      }
    }
  }
  catch (RVLException & e) {
    e.write(cerr);
    exit(1);
  }
}