module params
  implicit none

  !  Useful constants
  integer, parameter :: &
       fdp=kind(0.0d0), variablelength=256, errorlength=256

  real(fdp), parameter :: pi = 3.141592653589793238d0, pi2 = pi/2.0d0, &
       pi4 = pi/4.0d0, eps = 1.0d-16, inf = huge(inf), &
       el=0.5772156649015329d0, sqrt2 = 1.4142135623730950488d0, zero=tiny(zero)


  character(len=*), parameter :: &
       operatorchars='+-/*^><=!', numberchars='1234567890.', &
       exponentchars='ed', lowerchars='abcdefghijklmnopqrstuvwxyz', &
       upperchars='ABCDEFGHIJKLMNOPQRSTUVWXYZ', bracketchars='()[]{}', &
       specialchars=',_', funcvarchars=lowerchars//upperchars, &
       allowedchars=operatorchars//numberchars//funcvarchars//bracketchars//specialchars

  ! Controlling definitions
  integer, parameter :: &
       number_=10, bracket_=11, nop_=12, comma_=13, end_=14

  !  Bracket definitions
  integer, parameter :: &
       bracketstart=20, &
       roundbracket_=20, squarebracket_=21, squigbracket_=22, &
       bracketend=22

  !  Expression definitions
  integer, parameter :: &
       expressionlength=8, &
       expressionstart=100, &
       anint_=100, log10_=101, sqrt_=102, acos_=103, asin_=104, &
       atan_=105, cosh_=106, sinh_=107, tanh_=108, aint_=109, &
       cos_=110, sin_=111, tan_=112, exp_=113, log_=114, abs_=115, &
       delta_=116, step_=117, hat_=118, max_=119, min_=120, &
       besj0_=121, besj1_=122, besjn_=123, besy0_=124, besy1_=125, &
       besyn_=126, logn_=127, erf_=128, erfc_=129, lgamma_=130, &
       gamma_=131, csch_=132, sech_=133, coth_=134, if_=135, &
       gauss_=136, sinc_=137, besi0_=138, besi1_=139, besk0_=140, &
       besk1_=141, ierfc_=142, besin_=143, beskn_=144, cbrt_=145, &
       fresc_=146, fress_=147, expi_=148, sini_=149, cosi_=150, &
       logi_=151, ierf_=152, elle_=153, ellk_=154, ielle_=155, iellf_=156, &
       expressionend=156

  character(len=expressionlength), dimension(expressionstart:expressionend), parameter :: &
       expressions=(/ 'anint(  ', 'log10(  ', 'sqrt(   ', 'acos(   ', &
       'asin(   ', 'atan(   ', 'cosh(   ', 'sinh(   ', 'tanh(   ', 'aint(   ', &
       'cos(    ', 'sin(    ', 'tan(    ', 'exp(    ', 'log(    ', 'abs(    ', &
       'delta(  ', 'step(   ', 'hat(    ', 'max(    ', 'min(    ', 'besj0(  ', &
       'besj1(  ', 'besjn(  ', 'besy0(  ', 'besy1(  ', 'besyn(  ', 'logn(   ', &
       'erf(    ', 'erfc(   ', 'lgamma( ', 'gamma(  ', 'csch(   ', 'sech(   ', &
       'coth(   ', 'if(     ', 'gauss(  ', 'sinc(   ', 'besi0(  ', 'besi1(  ', &
       'besk0(  ', 'besk1(  ', 'ierfc(  ', 'besin(  ', 'beskn(  ', 'cbrt(   ', &
       'fresc(  ', 'fress(  ', 'expi(   ', 'sini(   ', 'cosi(   ', 'logi(   ', &
       'ierf(   ', 'elle(   ', 'ellk(   ', 'ielle(  ', 'iellf(  ' &
       /)
      

  integer, dimension(expressionstart:expressionend), parameter :: &
       expressionslen=len_trim(expressions), &
       commandcomps=(/ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
       1, 1, 1, 1, 1, 2, 2, 1, 1, 2, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, &
       1, 3, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, &
       1, 1, 1, 2, 2 /)
  integer, parameter :: &
       maxcommandcomps=4

  !  Operators definitions
  integer, parameter :: &
       operatorlength=2, &
       operatorstart=500, &
       conditionalstart=500, &
       neq_=500, eqeq_=501, gteq_=502, eqgt_=503, lteq_=504, eqlt_=505, eq_=506, gt_=507, lt_=508,&
       conditionalend=508, &
       opsstart=509, &
       pow_=509, caret_=510, mult_=511, div_=512, &
       opsend=512, &
       pmstart=513, &
       plus_=513, minus_=514, &
       pmend=514, &
       operatorend=514

  character(len=operatorlength), dimension(operatorstart:operatorend), parameter :: &
       operators=(/ '!=', '==', '>=', '=>', '<=', '=<', '= ', &
       '> ', '< ', '**', '^ ', '* ', '/ ', '+ ', '- ' /)

  integer, dimension(operatorstart:operatorend), parameter :: &
       operatorsorder=(/ 10, 10, 10, 10, 10, 10, 10, 10, 10, &
       10000, 10000, 1000, 1000, 100, 100 /), &
       operatorslen=len_trim(operators)

  !  Variable definitions
  integer, parameter :: &
       variablestart=1000, &
       variableend=9999

  ! Separates logical blocks of code
  integer, dimension(10), parameter :: blockseparator = (/ comma_, neq_, eqeq_, gteq_, &
       eqgt_, lteq_, eqlt_, eq_, gt_, lt_ /) 

end module params

module extrafunc
  !
  !  Extra functions, they must take real input and return real
  !

  !
  !  Many of these functions were taken from Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp)
  !  with modifications http://momonga.t.u-tokyo.ac.jp/~ooura/index.html
  !
  !  Some of these functions were programmed with ideas from Numerical Recipes
  !
  !  Other functions came from Shanjie Zhang and Jianming Jin and their specfun 
  !  functions http://jin.ece.uiuc.edu/routines/routines.html with modifications
  !

  use params

  implicit none
  
  private
  public :: f_besj0, f_besj1, f_besjn, f_besy0, f_besy1, f_besyn,&
       & f_erf, f_gamma, f_lgamma, f_cbrt, f_besi0, f_besi1, f_besin,&
       & f_besk0, f_besk1, f_beskn, f_ierfc, f_fresc, f_fress, &
       & f_expi, f_sini, f_cosi, f_logi, f_ierf, f_elle, f_ellk, f_ielle, f_iellf

contains

  recursive function f_besj0(x) result(f)
    !
    ! Bessel function 1st kind, order 0
    !

    real(fdp), intent(in) :: x
    real(fdp) :: f

    integer :: k
    real(fdp) :: t, w, v, y, theta
    real(fdp), dimension(0:7), parameter :: a = (/ &
         &-0.0000000000023655394d0, 0.0000000004708898680d0, &
         &-0.0000000678167892231d0, 0.0000067816840038636d0, &
         &-0.0004340277777716935d0, 0.0156249999999992397d0, &
         &-0.2499999999999999638d0, 0.9999999999999999997d0 /)
    real(fdp), dimension(0:64), parameter :: b = (/&
         & 0.0000000000626681117d0, -0.0000000022270614428d0,&
         & 0.0000000662981656302d0, -0.0000016268486502196d0,&
         & 0.0000321978384111685d0, -0.0005005237733315830d0,&
         & 0.0059060313537449816d0, -0.0505265323740109701d0,&
         & 0.2936432097610503985d0, -1.0482565081091638637d0,&
         & 1.9181123286040428113d0, -1.1319199475221700100d0, &
         &-0.1965480952704682000d0, 0.0000000000457457332d0, &
         &-0.0000000015814772025d0, 0.0000000455487446311d0, &
         &-0.0000010735201286233d0, 0.0000202015179970014d0, &
         &-0.0002942392368203808d0, 0.0031801987726150648d0, &
         &-0.0239875209742846362d0, 0.1141447698973777641d0, &
         &-0.2766726722823530233d0, 0.1088620480970941648d0,&
         & 0.5136514645381999197d0, -0.2100594022073706033d0,&
         & 0.0000000000331366618d0, -0.0000000011119090229d0,&
         & 0.0000000308823040363d0, -0.0000006956602653104d0,&
         & 0.0000123499947481762d0, -0.0001662951945396180d0,&
         & 0.0016048663165678412d0, -0.0100785479932760966d0,&
         & 0.0328996815223415274d0, -0.0056168761733860688d0, &
         &-0.2341096400274429386d0, 0.2551729256776404262d0,&
         & 0.2288438186148935667d0, 0.0000000000238007203d0, &
         &-0.0000000007731046439d0, 0.0000000206237001152d0, &
         &-0.0000004412291442285d0, 0.0000073107766249655d0, &
         &-0.0000891749801028666d0, 0.0007341654513841350d0, &
         &-0.0033303085445352071d0, 0.0015425853045205717d0,&
         & 0.0521100583113136379d0, -0.1334447768979217815d0, &
         &-0.1401330292364750968d0, 0.2685616168804818919d0,&
         & 0.0000000000169355950d0, -0.0000000005308092192d0,&
         & 0.0000000135323005576d0, -0.0000002726650587978d0,&
         & 0.0000041513240141760d0, -0.0000443353052220157d0,&
         & 0.0002815740758993879d0, -0.0004393235121629007d0, &
         &-0.0067573531105799347d0, 0.0369141914660130814d0,&
         & 0.0081673361942996237d0, -0.2573381285898881860d0,&
         & 0.0459580257102978932d0 /) 
    real(fdp), dimension(0:69), parameter :: c = (/ &
         &-0.00000000003009451757d0, -0.00000000014958003844d0,&
         & 0.00000000506854544776d0, 0.00000001863564222012d0, &
         &-0.00000060304249068078d0, -0.00000147686259937403d0,&
         & 0.00004714331342682714d0, 0.00006286305481740818d0, &
         &-0.00214137170594124344d0, -0.00089157336676889788d0,&
         & 0.04508258728666024989d0, -0.00490362805828762224d0, &
         &-0.27312196367405374426d0, 0.04193925184293450356d0, &
         &-0.00000000000712453560d0, -0.00000000041170814825d0,&
         & 0.00000000138012624364d0, 0.00000005704447670683d0, &
         &-0.00000019026363528842d0, -0.00000533925032409729d0,&
         & 0.00001736064885538091d0, 0.00030692619152608375d0, &
         &-0.00092598938200644367d0, -0.00917934265960017663d0,&
         & 0.02287952522866389076d0, 0.10545197546252853195d0, &
         &-0.16126443075752985095d0, -0.19392874768742235538d0,&
         & 0.00000000002128344556d0, -0.00000000031053910272d0, &
         &-0.00000000334979293158d0, 0.00000004507232895050d0,&
         & 0.00000036437959146427d0, -0.00000446421436266678d0, &
         &-0.00002523429344576552d0, 0.00027519882931758163d0,&
         & 0.00097185076358599358d0, -0.00898326746345390692d0, &
         &-0.01665959196063987584d0, 0.11456933464891967814d0,&
         & 0.07885001422733148815d0, -0.23664819446234712621d0,&
         & 0.00000000003035295055d0, 0.00000000005486066835d0, &
         &-0.00000000501026824811d0, -0.00000000501246847860d0,&
         & 0.00000058012340163034d0, 0.00000016788922416169d0, &
         &-0.00004373270270147275d0, 0.00001183898532719802d0,&
         & 0.00189863342862291449d0, -0.00113759249561636130d0, &
         &-0.03846797195329871681d0, 0.02389746880951420335d0,&
         & 0.22837862066532347461d0, -0.06765394811166522844d0,&
         & 0.00000000001279875977d0, 0.00000000035925958103d0, &
         &-0.00000000228037105967d0, -0.00000004852770517176d0,&
         & 0.00000028696428000189d0, 0.00000440131125178642d0, &
         &-0.00002366617753349105d0, -0.00024412456252884129d0,&
         & 0.00113028178539430542d0, 0.00708470513919789080d0, &
         &-0.02526914792327618386d0, -0.08006137953480093426d0,&
         & 0.16548380461475971846d0, 0.14688405470042110229d0 /)
    real(fdp), dimension(0:51), parameter :: d = (/&
         & 1.059601355592185731d-14, -2.71150591218550377d-13,&
         & 8.6514809056201638d-12, -4.6264028554286627d-10,&
         & 5.0815403835647104d-8, -1.76722552048141208d-5,&
         & 0.16286750396763997378d0, 2.949651820598278873d-13, &
         &-8.818215611676125741d-12, 3.571119876162253451d-10, &
         &-2.631924120993717060d-8, 4.709502795656698909d-6, &
         &-5.208333333333283282d-3, 7.18344107717531977d-15, &
         &-2.51623725588410308d-13, 8.6017784918920604d-12, &
         &-4.6256876614290359d-10, 5.0815343220437937d-8, &
         &-1.76722551764941970d-5, 0.16286750396763433767d0,&
         & 2.2327570859680094777d-13, -8.464594853517051292d-12,&
         & 3.563766464349055183d-10, -2.631843986737892965d-8,&
         & 4.709502342288659410d-6, -5.2083333332278466225d-3,&
         & 5.15413392842889366d-15, -2.27740238380640162d-13,&
         & 8.4827767197609014d-12, -4.6224753682737618d-10,&
         & 5.0814848128929134d-8, -1.76722547638767480d-5,&
         & 0.16286750396748926663d0, 1.7316195320192170887d-13, &
         &-7.971122772293919646d-12, 3.544039469911895749d-10, &
         &-2.631443902081701081d-8, 4.709498228695400603d-6, &
         &-5.2083333315143653610d-3, 3.84653681453798517d-15, &
         &-2.04464520778789011d-13, 8.3089298605177838d-12, &
         &-4.6155016158412096d-10, 5.0813263696466650d-8, &
         &-1.76722528311426167d-5, 0.16286750396650065930d0,&
         & 1.3797879972460878797d-13, -7.448089381011684812d-12,&
         & 3.512733797106959780d-10, -2.630500895563592722d-8,&
         & 4.709483934775839193d-6, -5.2083333227940760113d-3 /)

    w = abs(x)

    if (w  <  1) then
       t = w * w
       y = ((((((a(0) * t + a(1)) * t + a(2)) * t + a(3)) * t +&
            & a(4)) * t + a(5)) * t + a(6)) * t + a(7)
    else if (w  <  8.5d0) then
       t = w * w * 0.0625d0
       k = int(t)
       t = t - (k + 0.5d0)
       k = k * 13
       y = (((((((((((b(k) * t + b(k + 1)) * t + b(k + 2)) * t +&
            & b(k + 3)) * t + b(k + 4)) * t + b(k + 5)) * t + b(k&
            & + 6)) * t + b(k + 7)) * t + b(k + 8)) * t + b(k +&
            & 9)) * t + b(k + 10)) * t + b(k + 11)) * t + b(k +&
            & 12)
    else if (w  <  12.5d0) then
       k = int(w)
       t = w - (k + 0.5d0)