!=====================================================================
!
!               S p e c f e m 3 D  V e r s i o n  3 . 0
!               ---------------------------------------
!
!     Main historical authors: Dimitri Komatitsch and Jeroen Tromp
!                        Princeton University, USA
!                and CNRS / University of Marseille, France
!                 (there are currently many more authors!)
! (c) Princeton University and CNRS / University of Marseille, July 2012
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along
! with this program; if not, write to the Free Software Foundation, Inc.,
! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
!
!=====================================================================

  double precision function comp_source_time_function(t,hdur)

  use constants

  implicit none

  double precision t,hdur

  double precision, external :: netlib_specfun_erf

  ! quasi Heaviside, small Gaussian moment-rate tensor with hdur
  comp_source_time_function = 0.5d0*(1.0d0 + netlib_specfun_erf(t/hdur))

  end function comp_source_time_function


!
!-------------------------------------------------------------------------------------------------
!
!Carene

  double precision function comp_source_time_function_mono(t,hdur,ncycles)

  implicit none

  include "constants.h"

  double precision t,hdur,f0,ncycles

  double precision, external :: netlib_specfun_erf
!Carene
  double precision:: risingt,stopt

  risingt = 1.25d0/hdur
! cycles
  stopt = ncycles/hdur
  f0 = hdur

  !Mono chromatic
  if (t<=0.d0) then 
    comp_source_time_function_mono = 0.d0
  else
    !erf function stops at -4 theoretically but already "enough" at 0 at x=-3
    comp_source_time_function_mono = 0.5d0*(1.d0+netlib_specfun_erf((t-3.d0*risingt)/risingt))*dsin(2*pi*f0*t)
    comp_source_time_function_mono = 0.5d0*(1.d0 + netlib_specfun_erf((stopt-t)/risingt)) * comp_source_time_function_mono
  endif

  end function comp_source_time_function_mono

!
!-------------------------------------------------------------------------------------------------
!
!
!-------------------------------------------------------------------------------------------------
!
!Carene

  double precision function comp_source_time_function_mono2(t,hdur)

  implicit none

  include "constants.h"

  double precision t,hdur,f0

  double precision, external :: netlib_specfun_erf
!Carene
  double precision:: risingt,stopt,df,tmp
  double precision, parameter :: fmin=0.1d0,fmax=5.d0
  integer :: if
  integer, parameter :: nf=50

  risingt = 1.25d0/hdur
! 10 cycles
  stopt = 10.d0/hdur
  f0 = hdur

  !Mono chromatic
  if (t<=0.d0) then 
    comp_source_time_function_mono2 = 0.d0
  else
    !erf function stops at -4 theoretically but already "enough" at 0 at x=-3
    comp_source_time_function_mono2 = 0.5d0*(1.d0+netlib_specfun_erf((t-3.d0*risingt)/risingt))!*dsin(2*pi*f0*t)
    comp_source_time_function_mono2 = 0.5d0*(1.d0 + netlib_specfun_erf((stopt-t)/risingt)) * comp_source_time_function_mono2
  endif
!  df=(fmax/fmin)**(1/10.d0)
   df=(fmax-fmin)/dble(nf)
  do if = 0,nf
   !f0=fmin*(df**if)
   f0=fmin+if*df
   tmp = tmp+dsin(2*pi*f0*t)
  enddo
  comp_source_time_function_mono2 = comp_source_time_function_mono2*tmp

  end function comp_source_time_function_mono2

!
!-------------------------------------------------------------------------------------------------
!

  double precision function comp_source_time_function_gauss(t,hdur)

  use constants

  implicit none

  double precision :: t,hdur
  double precision :: hdur_decay
  double precision,parameter :: SOURCE_DECAY_STRONG = 2.0d0/SOURCE_DECAY_MIMIC_TRIANGLE

  ! note: hdur given is hdur_gaussian = hdur/SOURCE_DECAY_MIMIC_TRIANGLE
  !           and SOURCE_DECAY_MIMIC_TRIANGLE ~ 1.68
  hdur_decay = hdur

  ! this here uses a stronger gaussian decay rate (empirical value) to avoid non-zero onset times;
  ! however, it should mimik a triangle source time function...
  !hdur_decay = hdur  / SOURCE_DECAY_STRONG

  ! note: a nonzero time to start the simulation with would lead to more high-frequency noise
  !          due to the (spatial) discretization of the point source on the mesh

  ! gaussian
  comp_source_time_function_gauss = exp(-(t/hdur_decay)**2)/(sqrt(PI)*hdur_decay)

  end function comp_source_time_function_gauss

!
!-------------------------------------------------------------------------------------------------
!

  double precision function comp_source_time_function_rickr(t,f0)

  use constants

  implicit none

  double precision t,f0

  ! ricker
  comp_source_time_function_rickr = (1.d0 - 2.d0*PI*PI*f0*f0*t*t ) &
                                    * exp( -PI*PI*f0*f0*t*t )

  !!! another source time function they have called 'ricker' in some old papers,
  !!! e.g., 'Finite-Frequency Kernels Based on Adjoint Methods' by Liu & Tromp, BSSA (2006)
  !!! in order to benchmark those simulations, the following formula is needed.
  ! comp_source_time_function_rickr = -2.d0*PI*PI*f0*f0*f0*t * exp(-PI*PI*f0*f0*t*t)

  end function comp_source_time_function_rickr

!-------------------------------------------------------------------------------------------------
!
!Carene

  double precision function comp_source_time_function_custom(t,hdur,ISOURCE)
  
  use specfem_par, only : DT,nsteps_stf_custom,stf_custom,CUSTOM_REAL,SIZE_REAL

  implicit none

  double precision,intent(in) :: t,hdur
  integer, intent(in) :: ISOURCE

  integer :: it

  it = ceiling(t/DT)
  if (it<=nsteps_stf_custom) then 
    if (CUSTOM_REAL==SIZE_REAL) then 
      comp_source_time_function_custom=dble(stf_custom(ISOURCE,it))
    else
      comp_source_time_function_custom=stf_custom(ISOURCE,it)
    endif
  else
    comp_source_time_function_custom=stf_custom(ISOURCE,nsteps_stf_custom)
  endif

  end function comp_source_time_function_custom 
