// pntsrc_appx.cc
// Author: Mario J. Bencomo
// last modified: 11/10/16

#include "pntsrc_appx.hh"

#define MINUS_TO_POW(p) ((p)%2)?(-1):(1)

namespace TSOpt {

  //-----------------------------------------------------------------------------------
  int trunc_fact(int m,int s){
  //-----------------------------------------------------------------------------------
    if( m<0 || s<0 ){
      RVL::RVLException e;
      e << "ERROR from trunc_fact: s="<<s<<" or m="<<m<<" is negative!\n";
      throw e;
    }
    if( s>m ){
      RVL::RVLException e;
      e << "ERROR from trunc_fact: s="<<s<<" is greater than m="<<m<<"!\n";
      throw e;
    }
    int prod = 1;
    for(int i=m; i>(m-s); i--) prod*=i;
    return prod;
  }
  
  //-----------------------------------------------------------------------------------
  float pntsrc_weight( int s, int p, int indx, float alpha, float h){
  //-----------------------------------------------------------------------------------
    int Q = p + s -1;
    if(Q>Q_MAX){
      RVL::RVLException e;
      e << "ERROR from pntsrc_weight: Q="<<Q<<" exceeds maxQ="<<Q_MAX<<"!\n";
      throw e;
    }
    
    if( indx<0 || indx>Q ){
      return 0.0;
    }
    else{
      //compensating for shift in Q_mat for starting index  
      int start = 0;
      for(int i=1; i<=Q; i++) start += i*i;
      //compensating for row shift for starting index
      start+=indx*(Q+1);
      
      //computing weight
      float w   = 0.0;
      
      for(int col=0; col<=Q; col++){
	//computing b term
	float b = 0.0;
	if(col>=s){
	  //b =  aux * trunc_fact(col,s) * pow(alpha,col-s);
	  b = trunc_fact(col,s) * pow(alpha,col-s);
	}
	w += Q_MAT[start+col]*b;
      }

      //Modification: MJB 07-11-16
      //  In order to compensate for an extra h^-1 factor introduced via IWave
      //  during adjoint sampling of point source into FD spatial grid.
      //  OLD VERSION: 
      //w *= pow(h,-s-1)* pow(-1,s);
      w *= pow(h,-s) * pow(-1,s);//MINUS_TO_POW(s);

      return w;
    }
  }
  
  //-----------------------------------------------------------------------------------
  float pntsrc_weight( int sten_n, int s, int p, int indx, float alpha, float h){
    //-----------------------------------------------------------------------------------
    if(sten_n==1){
      return 1;
    }
    else{
      return pntsrc_weight(s,p,indx,alpha,h);
    }
  }
  
}